MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
change_basis.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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 "change_basis.hpp"
16
17namespace mfem
18{
19
20/// @brief Compute the inverse of the matrix A and store the result in Ainv.
21///
22/// The input A is an array of size n*n, interpreted as a matrix with column
23/// major ordering.
25{
26 Array<real_t> A2 = A;
27 const int n2 = A.Size();
28 const int n = sqrt(n2);
29 Array<int> ipiv(n);
30 LUFactors lu(A2.GetData(), ipiv.GetData());
31 lu.Factor(n);
32 Ainv.SetSize(n2);
33 lu.GetInverseMatrix(n, Ainv.GetData());
34}
35
36void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array<real_t> &B)
37{
40 Vector u(n);
41 B.SetSize(n*n);
42 B = 0.0;
43
44 for (int i = 0; i < n; ++i)
45 {
46 const real_t h = gll_pts[i+1] - gll_pts[i];
47 // Loop over subcell quadrature points
48 for (int iq = 0; iq < ir.Size(); ++iq)
49 {
50 const IntegrationPoint &ip = ir[iq];
51 const real_t x = gll_pts[i] + h*ip.x;
52 const real_t w = h*ip.weight;
53 basis.Eval(x, u);
54 for (int j = 0; j < n; ++j)
55 {
56 B[i + j*n] += w*u[j];
57 }
58 }
59 }
60}
61
63{
64 const int n = sqrt(B.Size());
65 Bt.SetSize(n*n);
66 for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
67}
68
70 : Operator(fes.GetTrueVSize()),
71 ne(fes.GetNE())
72{
73 auto *fec1 = dynamic_cast<const L2_FECollection*>(fes.FEColl());
74 MFEM_VERIFY(fec1, "Must be L2 finite element space");
75
76 const int btype = fec1->GetBasisType();
77
78 // If the basis types are the same, don't need to perform change of basis.
79 no_op = (btype == BasisType::IntegratedGLL);
80 if (no_op) { return; }
81
82 // Convert from the given basis to the "integrated GLL basis".
83 // The degrees of freedom are integrals over subcells.
84 const FiniteElement *fe = fes.GetFE(0);
85 auto *tbe = dynamic_cast<const TensorBasisElement*>(fe);
86 MFEM_VERIFY(tbe != nullptr, "Must be a tensor element.");
87 const Poly_1D::Basis &basis = tbe->GetBasis1D();
88
89 const int p = fes.GetMaxElementOrder();
90 const int pp1 = p + 1;
91
92 Array<real_t> B_inv;
93 SubcellIntegrals(pp1, basis, B_inv);
94
95 ComputeInverse(B_inv, B_1d);
96 Transpose(B_1d, Bt_1d);
97
98 // Set up the DofToQuad object, used in TensorValues
99 dof2quad.FE = fe;
100 dof2quad.mode = DofToQuad::TENSOR;
101 dof2quad.ndof = pp1;
102 dof2quad.nqpt = pp1;
103}
104
105void ChangeOfBasis_L2::Mult(const Vector &x, Vector &y) const
106{
107 if (no_op) { y = x; return; }
108 using namespace internal::quadrature_interpolator;
109 dof2quad.B.MakeRef(B_1d);
110 TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
111}
112
114{
115 if (no_op) { y = x; return; }
116 using namespace internal::quadrature_interpolator;
117 dof2quad.B.MakeRef(Bt_1d);
118 TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
119}
120
122 : Operator(fes.GetTrueVSize()),
123 fes(fes),
124 dim(fes.GetMesh()->Dimension()),
125 ne(fes.GetNE()),
126 p(fes.GetMaxElementOrder())
127{
129 elem_restr = dynamic_cast<const ElementRestriction*>(op);
130 MFEM_VERIFY(elem_restr != NULL, "Missing element restriction.");
131
132 const auto *rt_fec = dynamic_cast<const RT_FECollection*>(fes.FEColl());
133 MFEM_VERIFY(rt_fec, "Must be RT finite element space.");
134
135 const int cb_type = rt_fec->GetClosedBasisType();
136 const int ob_type = rt_fec->GetOpenBasisType();
137
138 no_op = (cb_type == BasisType::GaussLobatto &&
139 ob_type == BasisType::IntegratedGLL);
140 if (no_op) { return; }
141
142 const int pp1 = p + 1;
143
144 Poly_1D::Basis &cbasis = poly1d.GetBasis(p, cb_type);
145 Poly_1D::Basis &obasis = poly1d.GetBasis(p-1, ob_type);
146
148
149 Bci_1d.SetSize(pp1*pp1);
150 Vector b(pp1);
151 for (int i = 0; i < pp1; ++i)
152 {
153 cbasis.Eval(cpts2[i], b);
154 for (int j = 0; j < pp1; ++j)
155 {
156 Bci_1d[i + j*pp1] = b[j];
157 }
158 }
159 SubcellIntegrals(p, obasis, Boi_1d);
160
161 ComputeInverse(Boi_1d, Bo_1d);
162 Transpose(Bo_1d, Bot_1d);
163 ComputeInverse(Bci_1d, Bc_1d);
164 Transpose(Bc_1d, Bct_1d);
165}
166
167const real_t *ChangeOfBasis_RT::GetOpenMap(Mode mode) const
168{
169 switch (mode)
170 {
171 case NORMAL: return Bo_1d.Read();
172 case TRANSPOSE: return Bot_1d.Read();
173 case INVERSE: return Boi_1d.Read();
174 }
175 return nullptr;
176}
177
178const real_t *ChangeOfBasis_RT::GetClosedMap(Mode mode) const
179{
180 switch (mode)
181 {
182 case NORMAL: return Bc_1d.Read();
183 case TRANSPOSE: return Bct_1d.Read();
184 case INVERSE: return Bci_1d.Read();
185 }
186 return nullptr;
187}
188
189void ChangeOfBasis_RT::MultRT_2D(const Vector &x, Vector &y, Mode mode) const
190{
191 const int DIM = dim;
192 const int NE = ne;
193 const int D1D = p + 1;
194 const int ND = (p+1)*p;
195 const real_t *BC = GetClosedMap(mode);
196 const real_t *BO = GetOpenMap(mode);
197 const auto X = Reshape(x.Read(), DIM*ND, ne);
198 auto Y = Reshape(y.Write(), DIM*ND, ne);
199
200 MFEM_FORALL(e, NE,
201 {
202 for (int c = 0; c < DIM; ++c)
203 {
204 const int nx = (c == 0) ? D1D : D1D-1;
205 const int ny = (c == 1) ? D1D : D1D-1;
206 const real_t *Bx = (c == 0) ? BC : BO;
207 const real_t *By = (c == 1) ? BC : BO;
208
209 for (int i = 0; i < ND; ++i)
210 {
211 Y(i + c*ND, e) = 0.0;
212 }
213 for (int iy = 0; iy < ny; ++ iy)
214 {
215 real_t xx[DofQuadLimits::MAX_D1D];
216 for (int ix = 0; ix < nx; ++ix) { xx[ix] = 0.0; }
217 for (int jx = 0; jx < nx; ++jx)
218 {
219 const real_t val = X(jx + iy*nx + c*nx*ny, e);
220 for (int ix = 0; ix < nx; ++ix)
221 {
222 xx[ix] += val*Bx[ix + jx*nx];
223 }
224 }
225 for (int jy = 0; jy < ny; ++jy)
226 {
227 const real_t b = By[jy + iy*ny];
228 for (int ix = 0; ix < nx; ++ix)
229 {
230 Y(ix + jy*nx + c*nx*ny, e) += xx[ix]*b;
231 }
232 }
233 }
234 }
235 });
236}
237
238void ChangeOfBasis_RT::MultRT_3D(const Vector &x, Vector &y, Mode mode) const
239{
240 const int DIM = dim;
241 const int NE = ne;
242 const int D1D = p + 1;
243 const int ND = (p+1)*p*p;
244 const real_t *BC = GetClosedMap(mode);
245 const real_t *BO = GetOpenMap(mode);
246 const auto X = Reshape(x.Read(), DIM*ND, ne);
247 auto Y = Reshape(y.Write(), DIM*ND, ne);
248
249 MFEM_FORALL(e, NE,
250 {
251 for (int c = 0; c < DIM; ++c)
252 {
253 const int nx = (c == 0) ? D1D : D1D-1;
254 const int ny = (c == 1) ? D1D : D1D-1;
255 const int nz = (c == 2) ? D1D : D1D-1;
256 const real_t *Bx = (c == 0) ? BC : BO;
257 const real_t *By = (c == 1) ? BC : BO;
258 const real_t *Bz = (c == 2) ? BC : BO;
259
260 for (int i = 0; i < ND; ++i)
261 {
262 Y(i + c*ND, e) = 0.0;
263 }
264 for (int iz = 0; iz < nz; ++ iz)
265 {
266 real_t xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D];
267 for (int iy = 0; iy < ny; ++iy)
268 {
269 for (int ix = 0; ix < nx; ++ix)
270 {
271 xy[iy][ix] = 0.0;
272 }
273 }
274 for (int iy = 0; iy < ny; ++iy)
275 {
276 real_t xx[DofQuadLimits::MAX_D1D];
277 for (int ix = 0; ix < nx; ++ix) { xx[ix] = 0.0; }
278 for (int ix = 0; ix < nx; ++ix)
279 {
280 const real_t val = X(ix + iy*nx + iz*nx*ny + c*ND, e);
281 for (int jx = 0; jx < nx; ++jx)
282 {
283 xx[jx] += val*Bx[jx + ix*nx];
284 }
285 }
286 for (int jy = 0; jy < ny; ++jy)
287 {
288 const real_t b = By[jy + iy*ny];
289 for (int jx = 0; jx < nx; ++jx)
290 {
291 xy[jy][jx] += xx[jx] * b;
292 }
293 }
294 }
295 for (int jz = 0; jz < nz; ++jz)
296 {
297 const real_t b = Bz[jz + iz*nz];
298 for (int jy = 0; jy < ny; ++jy)
299 {
300 for (int jx = 0; jx < nx; ++jx)
301 {
302 Y(jx + jy*nx + jz*nx*ny + c*ND, e) += xy[jy][jx] * b;
303 }
304 }
305 }
306 }
307 }
308 });
309}
310
311void ChangeOfBasis_RT::Mult(const Vector &x, Vector &y, Mode mode) const
312{
313 if (no_op) { y = x; return; }
314
315 const Operator *P = fes.GetProlongationMatrix();
316
318 {
319 x_l.MakeRef(const_cast<Vector&>(x), 0, fes.GetVSize());
320 y_l.MakeRef(y, 0, fes.GetVSize());
321 }
322 else
323 {
324 x_l.SetSize(fes.GetVSize());
325 y_l.SetSize(fes.GetVSize());
326 P->Mult(x, x_l);
327 }
328
329 x_e.SetSize(elem_restr->Height());
330 y_e.SetSize(elem_restr->Height());
331
332 elem_restr->Mult(x_l, x_e);
333
334 if (dim == 2) { MultRT_2D(x_e, y_e, mode); }
335 else { MultRT_3D(x_e, y_e, mode); }
336
337 elem_restr->MultLeftInverse(y_e, y_l);
338
339 const Operator *R = fes.GetRestrictionOperator();
340 if (R) { R->Mult(y_l, y); }
341 else { MFEM_VERIFY(P == NULL, "Invalid state."); }
342}
343
344void ChangeOfBasis_RT::Mult(const Vector &x, Vector &y) const
345{
346 Mult(x, y, NORMAL);
347}
348
350{
351 Mult(x, y, TRANSPOSE);
352}
353
355{
356 Mult(x, y, INVERSE);
357}
358
359} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void MakeRef(T *data_, int size_, bool own_data=false)
Make this Array a reference to a pointer.
Definition array.hpp:882
T * GetData()
Returns the data.
Definition array.hpp:118
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:317
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:39
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
ChangeOfBasis_L2(FiniteElementSpace &fes)
ChangeOfBasis_RT(FiniteElementSpace &fes)
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void MultInverse(const Vector &x, Vector &y) const
void MultRT_2D(const Vector &x, Vector &y, Mode mode) const
void MultRT_3D(const Vector &x, Vector &y, Mode mode) const
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition fe_base.hpp:170
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:161
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:189
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:174
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition fe_base.hpp:141
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
Operator that converts FiniteElementSpace L-vectors to E-vectors.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
void MultLeftInverse(const Vector &x, Vector &y) const
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition fespace.hpp:597
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
Definition fespace.hpp:616
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1336
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
int GetMaxElementOrder() const
Return the maximum polynomial order.
Definition fespace.hpp:577
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
Abstract class for all finite elements.
Definition fe_base.hpp:239
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
virtual void GetInverseMatrix(int m, real_t *X) const
Assuming L.U = P.A factored data of size (m x m), compute X <- A^{-1}.
virtual bool Factor(int m, real_t TOL=0.0)
Compute the LU factorization of the current matrix.
Abstract operator.
Definition operator.hpp:25
Operator(int s=0)
Construct a square Operator with given size s (default 0).
Definition operator.hpp:59
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition fe_base.hpp:991
void Eval(const real_t x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition fe_base.cpp:1761
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition fe_base.cpp:2304
const real_t * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition fe_base.cpp:2270
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Vector data type.
Definition vector.hpp:80
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:474
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:482
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
int dim
Definition ex24.cpp:53
Mesh * GetMesh(int type)
Definition ex29.cpp:218
real_t b
Definition lissajous.cpp:42
constexpr int DIM
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array< real_t > &B)
void ComputeInverse(const Array< real_t > &A, Array< real_t > &Ainv)
Compute the inverse of the matrix A and store the result in Ainv.
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:414
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
bool IsIdentityProlongation(const Operator *P)
Definition operator.hpp:720
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
@ LEXICOGRAPHIC
Lexicographic ordering for tensor-product FiniteElements.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
real_t p(const Vector &x, real_t t)