1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "bilininteg.hpp"
13 #include "pfespace.hpp"
14 #include <algorithm>
15
16 namespace mfem
17 {
18
DGDiffusionBR2Integrator(FiniteElementSpace * fes,double e)19 DGDiffusionBR2Integrator::DGDiffusionBR2Integrator(FiniteElementSpace *fes,
20 double e) : eta(e)
21 {
22 // Precompute local mass matrix inverses needed for the lifting operators
23 // First compute offsets and total size needed (e.g. for mixed meshes or
24 // p-refinement)
25 int nel = fes->GetNE();
26 Minv_offsets.SetSize(nel+1);
27 ipiv_offsets.SetSize(nel+1);
28 ipiv_offsets[0] = 0;
29 Minv_offsets[0] = 0;
30 for (int i=0; i<nel; ++i)
31 {
32 int dof = fes->GetFE(i)->GetDof();
33 ipiv_offsets[i+1] = ipiv_offsets[i] + dof;
34 Minv_offsets[i+1] = Minv_offsets[i] + dof*dof;
35 }
36
37 #ifdef MFEM_USE_MPI
38 // When running in parallel, we also need to compute the local mass matrices
39 // of face neighbor elements
40 ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace *>(fes);
41 if (pfes != NULL)
42 {
43 ParMesh *pmesh = pfes->GetParMesh();
44 pfes->ExchangeFaceNbrData();
45 int nel_nbr = pmesh->GetNFaceNeighborElements();
46 Minv_offsets.SetSize(nel+nel_nbr+1);
47 ipiv_offsets.SetSize(nel+nel_nbr+1);
48 for (int i=0; i<nel_nbr; ++i)
49 {
50 int dof = pfes->GetFaceNbrFE(i)->GetDof();
51 ipiv_offsets[nel+i+1] = ipiv_offsets[nel+i] + dof;
52 Minv_offsets[nel+i+1] = Minv_offsets[nel+i] + dof*dof;
53 }
54 nel += nel_nbr;
55 }
56 #endif
57 // The final "offset" is the total size of all the blocks
58 Minv.SetSize(Minv_offsets[nel]);
59 ipiv.SetSize(ipiv_offsets[nel]);
60
61 // Assemble the local mass matrices and compute LU factorization
62 MassIntegrator mi;
63 for (int i=0; i<nel; ++i)
64 {
65 const FiniteElement *fe = NULL;
66 ElementTransformation *tr = NULL;
67 if (i < fes->GetNE())
68 {
69 fe = fes->GetFE(i);
70 tr = fes->GetElementTransformation(i);
71 }
72 else
73 {
74 #ifdef MFEM_USE_MPI
75 int inbr = i - fes->GetNE();
76 fe = pfes->GetFaceNbrFE(inbr);
77 tr = pfes->GetParMesh()->GetFaceNbrElementTransformation(inbr);
78 #endif
79 }
80 int dof = fe->GetDof();
81 double *Minv_el = &Minv[Minv_offsets[i]];
82 int *ipiv_el = &ipiv[ipiv_offsets[i]];
83 DenseMatrix Me(Minv_el, dof, dof);
84 mi.AssembleElementMatrix(*fe, *tr, Me);
85 LUFactors lu(Minv_el, ipiv_el);
86 lu.Factor(dof);
87 }
88 }
89
AssembleFaceMatrix(const FiniteElement & el1,const FiniteElement & el2,FaceElementTransformations & Trans,DenseMatrix & elmat)90 void DGDiffusionBR2Integrator::AssembleFaceMatrix(
91 const FiniteElement &el1, const FiniteElement &el2,
92 FaceElementTransformations &Trans, DenseMatrix &elmat)
93 {
94 int ndof1 = el1.GetDof();
95 shape1.SetSize(ndof1);
96
97 R11.SetSize(ndof1, ndof1);
98 R11 = 0.0;
99 LUFactors M1inv(&Minv[Minv_offsets[Trans.Elem1No]],
100 &ipiv[ipiv_offsets[Trans.Elem1No]]);
101 LUFactors M2inv;
102
103 double factor = Geometries.NumBdr(Trans.Elem1->GetGeometryType());
104
105 int ndof2;
106 if (Trans.Elem2No >= 0)
107 {
108 ndof2 = el2.GetDof();
109 shape2.SetSize(ndof2);
110 R12.SetSize(ndof1, ndof2);
111 R21.SetSize(ndof2, ndof1);
112 R22.SetSize(ndof2, ndof2);
113 M2inv.data = &Minv[Minv_offsets[Trans.Elem2No]];
114 M2inv.ipiv = &ipiv[ipiv_offsets[Trans.Elem2No]];
115
116 R12 = 0.0;
117 R21 = 0.0;
118 R22 = 0.0;
119
120 Geometry::Type geom2 = Trans.Elem2->GetGeometryType();
121 factor = std::max(factor, double(Geometries.NumBdr(geom2)));
122 }
123 else
124 {
125 ndof2 = 0;
126 }
127
128 int ndofs = ndof1 + ndof2;
129
130 Re.SetSize(ndofs, ndofs);
131 MinvRe.SetSize(ndofs, ndofs);
132
133 elmat.SetSize(ndofs);
134 elmat = 0.0;
135
136 const IntegrationRule *ir = IntRule;
137 if (ir == NULL)
138 {
139 int order;
140 if (ndof2)
141 {
142 order = 2*std::max(el1.GetOrder(), el2.GetOrder());
143 }
144 else
145 {
146 order = 2*el1.GetOrder();
147 }
148 ir = &IntRules.Get(Trans.FaceGeom, order);
149 }
150
151 for (int p = 0; p < ir->GetNPoints(); p++)
152 {
153 const IntegrationPoint &ip = ir->IntPoint(p);
154 IntegrationPoint eip1, eip2;
155
156 Trans.Loc1.Transform(ip, eip1);
157 el1.CalcShape(eip1, shape1);
158 if (ndof2)
159 {
160 Trans.Loc2.Transform(ip, eip2);
161 el2.CalcShape(eip2, shape2);
162 }
163
164 double w = factor*sqrt(eta)*ip.weight*Trans.Face->Weight();
165 if (ndof2)
166 {
167 w /= 2;
168 }
169
170 for (int i = 0; i < ndof1; i++)
171 {
172 const double wsi = w*shape1(i);
173 for (int j = 0; j < ndof1; j++)
174 {
175 R11(i, j) += wsi*shape1(j);
176 }
177 }
178
179 if (ndof2)
180 {
181 for (int i = 0; i < ndof2; i++)
182 {
183 const double wsi = w*shape2(i);
184 for (int j = 0; j < ndof1; j++)
185 {
186 R21(i, j) += wsi*shape1(j);
187 R12(j, i) -= wsi*shape1(j);
188 }
189 for (int j = 0; j < ndof2; j++)
190 {
191 R22(i, j) -= wsi*shape2(j);
192 }
193 }
194 }
195 }
196
197 MinvR11 = R11;
198 M1inv.Solve(ndof1, ndof1, MinvR11.Data());
199 for (int i = 0; i < ndof1; i++)
200 {
201 for (int j = 0; j < ndof1; j++)
202 {
203 Re(i, j) = R11(i, j);
204 MinvRe(i, j) = MinvR11(i, j);
205 }
206 }
207
208 if (ndof2)
209 {
210 MinvR12 = R12;
211 MinvR21 = R21;
212 MinvR22 = R22;
213 M1inv.Solve(ndof1, ndof2, MinvR12.Data());
214 M2inv.Solve(ndof2, ndof1, MinvR21.Data());
215 M2inv.Solve(ndof2, ndof2, MinvR22.Data());
216
217 for (int i = 0; i < ndof2; i++)
218 {
219 for (int j = 0; j < ndof1; j++)
220 {
221 Re(ndof1 + i, j) = R21(i, j);
222 MinvRe(ndof1 + i, j) = MinvR21(i, j);
223
224 Re(j, ndof1 + i) = R12(j, i);
225 MinvRe(j, ndof1 + i) = MinvR12(j, i);
226 }
227 for (int j = 0; j < ndof2; j++)
228 {
229 Re(ndof1 + i, ndof1 + j) = R22(i, j);
230 MinvRe(ndof1 + i, ndof1 + j) = MinvR22(i, j);
231 }
232 }
233 }
234
235 // Compute the matrix associated with (r_e([u]), r_e([u])).
236 // The matrix for r_e([u]) is `MinvRe`, and so we need to form the product
237 // `(MinvRe)^T M MinvRe`. Using `Minv^T M = Minv M = I`, we obtain
238 // `Re^T MinvRe`.
239 MultAtB(Re, MinvRe, elmat);
240 }
241
242 }
243