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