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 #ifndef MFEM_FESPACE
13 #define MFEM_FESPACE
14
15 #include "../config/config.hpp"
16 #include "../linalg/sparsemat.hpp"
17 #include "../mesh/mesh.hpp"
18 #include "fe_coll.hpp"
19 #include "restriction.hpp"
20 #include <iostream>
21 #include <unordered_map>
22
23 namespace mfem
24 {
25
26 /** @brief The ordering method used when the number of unknowns per mesh node
27 (vector dimension) is bigger than 1. */
28 class Ordering
29 {
30 public:
31 /// %Ordering methods:
32 enum Type
33 {
34 byNODES, /**< loop first over the nodes (inner loop) then over the vector
35 dimension (outer loop); symbolically it can be represented
36 as: XXX...,YYY...,ZZZ... */
37 byVDIM /**< loop first over the vector dimension (inner loop) then over
38 the nodes (outer loop); symbolically it can be represented
39 as: XYZ,XYZ,XYZ,... */
40 };
41
42 template <Type Ord>
43 static inline int Map(int ndofs, int vdim, int dof, int vd);
44
45 template <Type Ord>
46 static void DofsToVDofs(int ndofs, int vdim, Array<int> &dofs);
47 };
48
49 template <> inline int
Map(int ndofs,int vdim,int dof,int vd)50 Ordering::Map<Ordering::byNODES>(int ndofs, int vdim, int dof, int vd)
51 {
52 MFEM_ASSERT(dof < ndofs && -1-dof < ndofs && 0 <= vd && vd < vdim, "");
53 return (dof >= 0) ? dof+ndofs*vd : dof-ndofs*vd;
54 }
55
56 template <> inline int
Map(int ndofs,int vdim,int dof,int vd)57 Ordering::Map<Ordering::byVDIM>(int ndofs, int vdim, int dof, int vd)
58 {
59 MFEM_ASSERT(dof < ndofs && -1-dof < ndofs && 0 <= vd && vd < vdim, "");
60 return (dof >= 0) ? vd+vdim*dof : -1-(vd+vdim*(-1-dof));
61 }
62
63
64 /// Constants describing the possible orderings of the DOFs in one element.
65 enum class ElementDofOrdering
66 {
67 /// Native ordering as defined by the FiniteElement.
68 /** This ordering can be used by tensor-product elements when the
69 interpolation from the DOFs to quadrature points does not use the
70 tensor-product structure. */
71 NATIVE,
72 /// Lexicographic ordering for tensor-product FiniteElements.
73 /** This ordering can be used only with tensor-product elements. */
74 LEXICOGRAPHIC
75 };
76
77 // Forward declarations
78 class NURBSExtension;
79 class BilinearFormIntegrator;
80 class QuadratureSpace;
81 class QuadratureInterpolator;
82 class FaceQuadratureInterpolator;
83
84
85 /** @brief Class FiniteElementSpace - responsible for providing FEM view of the
86 mesh, mainly managing the set of degrees of freedom. */
87 class FiniteElementSpace
88 {
89 friend class InterpolationGridTransfer;
90 friend class PRefinementTransferOperator;
91 friend void Mesh::Swap(Mesh &, bool);
92 friend class LORBase;
93
94 protected:
95 /// The mesh that FE space lives on (not owned).
96 Mesh *mesh;
97
98 /// Associated FE collection (not owned).
99 const FiniteElementCollection *fec;
100
101 /// %Vector dimension (number of unknowns per degree of freedom).
102 int vdim;
103
104 /** Type of ordering of the vector dofs when #vdim > 1.
105 - Ordering::byNODES - first nodes, then vector dimension,
106 - Ordering::byVDIM - first vector dimension, then nodes */
107 Ordering::Type ordering;
108
109 /// Number of degrees of freedom. Number of unknowns is #ndofs * #vdim.
110 int ndofs;
111
112 /** Polynomial order for each element. If empty, all elements are assumed
113 to be of the default order (fec->GetOrder()). */
114 Array<char> elem_order;
115
116 int nvdofs, nedofs, nfdofs, nbdofs;
117 int uni_fdof; ///< # of single face DOFs if all faces uniform; -1 otherwise
118 int *bdofs; ///< internal DOFs of elements if mixed/var-order; NULL otherwise
119
120 /** Variable order spaces only: DOF assignments for edges and faces, see
121 docs in MakeDofTable. For constant order spaces the tables are empty. */
122 Table var_edge_dofs;
123 Table var_face_dofs; ///< NOTE: also used for spaces with mixed faces
124
125 /** Additional data for the var_*_dofs tables: individual variant orders
126 (these are basically alternate J arrays for var_edge/face_dofs). */
127 Array<char> var_edge_orders, var_face_orders;
128
129 // precalculated DOFs for each element, boundary element, and face
130 mutable Table *elem_dof; // owned (except in NURBS FE space)
131 mutable Table *bdr_elem_dof; // owned (except in NURBS FE space)
132 mutable Table *face_dof; // owned; in var-order space contains variant 0 DOFs
133
134 Array<int> dof_elem_array, dof_ldof_array;
135
136 NURBSExtension *NURBSext;
137 int own_ext;
138 mutable Array<int> face_to_be; // NURBS FE space only
139
140 /** Matrix representing the prolongation from the global conforming dofs to
141 a set of intermediate partially conforming dofs, e.g. the dofs associated
142 with a "cut" space on a non-conforming mesh. */
143 mutable SparseMatrix *cP; // owned
144 /// Conforming restriction matrix such that cR.cP=I.
145 mutable SparseMatrix *cR; // owned
146 /// A version of the conforming restriction matrix for variable-order spaces.
147 mutable SparseMatrix *cR_hp; // owned
148 mutable bool cP_is_set;
149
150 /// Transformation to apply to GridFunctions after space Update().
151 OperatorHandle Th;
152
153 /// The element restriction operators, see GetElementRestriction().
154 mutable OperatorHandle L2E_nat, L2E_lex;
155 /// The face restriction operators, see GetFaceRestriction().
156 using key_face = std::tuple<bool, ElementDofOrdering, FaceType, L2FaceValues>;
157 struct key_hash
158 {
operator ()mfem::FiniteElementSpace::key_hash159 std::size_t operator()(const key_face& k) const
160 {
161 return std::get<0>(k)
162 + 2 * (int)std::get<1>(k)
163 + 4 * (int)std::get<2>(k)
164 + 8 * (int)std::get<3>(k);
165 }
166 };
167 using map_L2F = std::unordered_map<const key_face,FaceRestriction*,key_hash>;
168 mutable map_L2F L2F;
169
170 mutable Array<QuadratureInterpolator*> E2Q_array;
171 mutable Array<FaceQuadratureInterpolator*> E2IFQ_array;
172 mutable Array<FaceQuadratureInterpolator*> E2BFQ_array;
173
174 /** Update counter, incremented every time the space is constructed/updated.
175 Used by GridFunctions to check if they are up to date with the space. */
176 long sequence;
177
178 /** Mesh sequence number last seen when constructing the space. The space
179 needs updating if Mesh::GetSequence() is larger than this. */
180 long mesh_sequence;
181
182 /// True if at least one element order changed (variable-order space only).
183 bool orders_changed;
184
185 bool relaxed_hp; // see SetRelaxedHpConformity()
186
187 void UpdateNURBS();
188
189 void Construct();
190 void Destroy();
191
192 void BuildElementToDofTable() const;
193 void BuildBdrElementToDofTable() const;
194 void BuildFaceToDofTable() const;
195
196 /** @brief Generates partial face_dof table for a NURBS space.
197
198 The table is only defined for exterior faces that coincide with a
199 boundary. */
200 void BuildNURBSFaceToDofTable() const;
201
202 /// Bit-mask representing a set of orders needed by an edge/face.
203 typedef std::uint64_t VarOrderBits;
204 static constexpr int MaxVarOrder = 8*sizeof(VarOrderBits) - 1;
205
206 /// Return the minimum order (least significant bit set) in the bit mask.
207 static int MinOrder(VarOrderBits bits);
208
209 /// Return element order: internal version of GetElementOrder without checks.
210 int GetElementOrderImpl(int i) const;
211
212 /** In a variable order space, calculate a bitmask of polynomial orders that
213 need to be represented on each edge and face. */
214 void CalcEdgeFaceVarOrders(Array<VarOrderBits> &edge_orders,
215 Array<VarOrderBits> &face_orders) const;
216
217 /** Build the table var_edge_dofs (or var_face_dofs) in a variable order
218 space; return total edge/face DOFs. */
219 int MakeDofTable(int ent_dim, const Array<int> &entity_orders,
220 Table &entity_dofs, Array<char> *var_ent_order);
221
222 /// Search row of a DOF table for a DOF set of size 'ndof', return first DOF.
223 int FindDofs(const Table &var_dof_table, int row, int ndof) const;
224
225 /** In a variable order space, return edge DOFs associated with a polynomial
226 order that has 'ndof' degrees of freedom. */
FindEdgeDof(int edge,int ndof) const227 int FindEdgeDof(int edge, int ndof) const
228 { return FindDofs(var_edge_dofs, edge, ndof); }
229
230 /// Similar to FindEdgeDof, but used for mixed meshes too.
FindFaceDof(int face,int ndof) const231 int FindFaceDof(int face, int ndof) const
232 { return FindDofs(var_face_dofs, face, ndof); }
233
FirstFaceDof(int face,int variant=0) const234 int FirstFaceDof(int face, int variant = 0) const
235 { return uni_fdof >= 0 ? face*uni_fdof : var_face_dofs.GetRow(face)[variant];}
236
237 /// Return number of possible DOF variants for edge/face (var. order spaces).
238 int GetNVariants(int entity, int index) const;
239
240 /// Helper to encode a sign flip into a DOF index (for Hcurl/Hdiv shapes).
EncodeDof(int entity_base,int idx)241 static inline int EncodeDof(int entity_base, int idx)
242 { return (idx >= 0) ? (entity_base + idx) : (-1-(entity_base + (-1-idx))); }
243
244 /// Helpers to remove encoded sign from a DOF
DecodeDof(int dof)245 static inline int DecodeDof(int dof)
246 { return (dof >= 0) ? dof : (-1 - dof); }
247
DecodeDof(int dof,double & sign)248 static inline int DecodeDof(int dof, double& sign)
249 { return (dof >= 0) ? (sign = 1, dof) : (sign = -1, (-1 - dof)); }
250
251 /// Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
252 int GetEntityDofs(int entity, int index, Array<int> &dofs,
253 Geometry::Type master_geom = Geometry::INVALID,
254 int variant = 0) const;
255
256 // Get degenerate face DOFs: see explanation in method implementation.
257 int GetDegenerateFaceDofs(int index, Array<int> &dofs,
258 Geometry::Type master_geom, int variant) const;
259
260 int GetNumBorderDofs(Geometry::Type geom, int order) const;
261
262 /// Calculate the cP and cR matrices for a nonconforming mesh.
263 void BuildConformingInterpolation() const;
264
265 static void AddDependencies(SparseMatrix& deps, Array<int>& master_dofs,
266 Array<int>& slave_dofs, DenseMatrix& I,
267 int skipfirst = 0);
268
269 static bool DofFinalizable(int dof, const Array<bool>& finalized,
270 const SparseMatrix& deps);
271
272 void AddEdgeFaceDependencies(SparseMatrix &deps, Array<int>& master_dofs,
273 const FiniteElement *master_fe,
274 Array<int> &slave_dofs, int slave_face,
275 const DenseMatrix *pm) const;
276
277 /// Replicate 'mat' in the vector dimension, according to vdim ordering mode.
278 void MakeVDimMatrix(SparseMatrix &mat) const;
279
280 /// GridFunction interpolation operator applicable after mesh refinement.
281 class RefinementOperator : public Operator
282 {
283 const FiniteElementSpace* fespace;
284 DenseTensor localP[Geometry::NumGeom];
285 Table* old_elem_dof; // Owned.
286
287 public:
288 /** Construct the operator based on the elem_dof table of the original
289 (coarse) space. The class takes ownership of the table. */
290 RefinementOperator(const FiniteElementSpace* fespace,
291 Table *old_elem_dof/*takes ownership*/, int old_ndofs);
292 RefinementOperator(const FiniteElementSpace *fespace,
293 const FiniteElementSpace *coarse_fes);
294 virtual void Mult(const Vector &x, Vector &y) const;
295 virtual void MultTranspose(const Vector &x, Vector &y) const;
296 virtual ~RefinementOperator();
297 };
298
299 /// Derefinement operator, used by the friend class InterpolationGridTransfer.
300 class DerefinementOperator : public Operator
301 {
302 const FiniteElementSpace *fine_fes; // Not owned.
303 DenseTensor localR[Geometry::NumGeom];
304 Table *coarse_elem_dof; // Owned.
305 Table coarse_to_fine;
306 Array<int> coarse_to_ref_type;
307 Array<Geometry::Type> ref_type_to_geom;
308 Array<int> ref_type_to_fine_elem_offset;
309
310 public:
311 DerefinementOperator(const FiniteElementSpace *f_fes,
312 const FiniteElementSpace *c_fes,
313 BilinearFormIntegrator *mass_integ);
314 virtual void Mult(const Vector &x, Vector &y) const;
315 virtual ~DerefinementOperator();
316 };
317
318 /** This method makes the same assumptions as the method:
319 void GetLocalRefinementMatrices(
320 const FiniteElementSpace &coarse_fes, Geometry::Type geom,
321 DenseTensor &localP) const
322 which is defined below. It also assumes that the coarse fes and this have
323 the same vector dimension, vdim. */
324 SparseMatrix *RefinementMatrix_main(const int coarse_ndofs,
325 const Table &coarse_elem_dof,
326 const DenseTensor localP[]) const;
327
328 void GetLocalRefinementMatrices(Geometry::Type geom,
329 DenseTensor &localP) const;
330 void GetLocalDerefinementMatrices(Geometry::Type geom,
331 DenseTensor &localR) const;
332
333 /** Calculate explicit GridFunction interpolation matrix (after mesh
334 refinement). NOTE: consider using the RefinementOperator class instead
335 of the fully assembled matrix, which can take a lot of memory. */
336 SparseMatrix* RefinementMatrix(int old_ndofs, const Table* old_elem_dof);
337
338 /// Calculate GridFunction restriction matrix after mesh derefinement.
339 SparseMatrix* DerefinementMatrix(int old_ndofs, const Table* old_elem_dof);
340
341 /** @brief Return in @a localP the local refinement matrices that map
342 between fespaces after mesh refinement. */
343 /** This method assumes that this->mesh is a refinement of coarse_fes->mesh
344 and that the CoarseFineTransformations of this->mesh are set accordingly.
345 Another assumption is that the FEs of this use the same MapType as the FEs
346 of coarse_fes. Finally, it assumes that the spaces this and coarse_fes are
347 NOT variable-order spaces. */
348 void GetLocalRefinementMatrices(const FiniteElementSpace &coarse_fes,
349 Geometry::Type geom,
350 DenseTensor &localP) const;
351
352 /// Help function for constructors + Load().
353 void Constructor(Mesh *mesh, NURBSExtension *ext,
354 const FiniteElementCollection *fec,
355 int vdim = 1, int ordering = Ordering::byNODES);
356
357 /// Updates the internal mesh pointer. @warning @a new_mesh must be
358 /// <b>topologically identical</b> to the existing mesh. Used if the address
359 /// of the Mesh object has changed, e.g. in @a Mesh::Swap.
360 virtual void UpdateMeshPointer(Mesh *new_mesh);
361
362 /// Resize the elem_order array on mesh change.
363 void UpdateElementOrders();
364
365 /// @brief Copies the prolongation and restriction matrices from @a fes.
366 ///
367 /// Used for low order preconditioning on non-conforming meshes. If the DOFs
368 /// require a permutation, it will be supplied by non-NULL @a perm. NULL @a
369 /// perm indicates that no permutation is required.
370 virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
371 const Array<int> *perm);
372
373 public:
374 /** @brief Default constructor: the object is invalid until initialized using
375 the method Load(). */
376 FiniteElementSpace();
377
378 /** @brief Copy constructor: deep copy all data from @a orig except the Mesh,
379 the FiniteElementCollection, ans some derived data. */
380 /** If the @a mesh or @a fec pointers are NULL (default), then the new
381 FiniteElementSpace will reuse the respective pointers from @a orig. If
382 any of these pointers is not NULL, the given pointer will be used instead
383 of the one used by @a orig.
384
385 @note The objects pointed to by the @a mesh and @a fec parameters must be
386 either the same objects as the ones used by @a orig, or copies of them.
387 Otherwise, the behavior is undefined.
388
389 @note Derived data objects, such as the conforming prolongation and
390 restriction matrices, and the update operator, will not be copied, even
391 if they are created in the @a orig object. */
392 FiniteElementSpace(const FiniteElementSpace &orig, Mesh *mesh = NULL,
393 const FiniteElementCollection *fec = NULL);
394
FiniteElementSpace(Mesh * mesh,const FiniteElementCollection * fec,int vdim=1,int ordering=Ordering::byNODES)395 FiniteElementSpace(Mesh *mesh,
396 const FiniteElementCollection *fec,
397 int vdim = 1, int ordering = Ordering::byNODES)
398 { Constructor(mesh, NULL, fec, vdim, ordering); }
399
400 /// Construct a NURBS FE space based on the given NURBSExtension, @a ext.
401 /** @note If the pointer @a ext is NULL, this constructor is equivalent to
402 the standard constructor with the same arguments minus the
403 NURBSExtension, @a ext. */
FiniteElementSpace(Mesh * mesh,NURBSExtension * ext,const FiniteElementCollection * fec,int vdim=1,int ordering=Ordering::byNODES)404 FiniteElementSpace(Mesh *mesh, NURBSExtension *ext,
405 const FiniteElementCollection *fec,
406 int vdim = 1, int ordering = Ordering::byNODES)
407 { Constructor(mesh, ext, fec, vdim, ordering); }
408
409 /// Returns the mesh
GetMesh() const410 inline Mesh *GetMesh() const { return mesh; }
411
GetNURBSext() const412 const NURBSExtension *GetNURBSext() const { return NURBSext; }
GetNURBSext()413 NURBSExtension *GetNURBSext() { return NURBSext; }
414 NURBSExtension *StealNURBSext();
415
Conforming() const416 bool Conforming() const { return mesh->Conforming() && cP == NULL; }
Nonconforming() const417 bool Nonconforming() const { return mesh->Nonconforming() || cP != NULL; }
418
419 /// Sets the order of the i'th finite element.
420 /** By default, all elements are assumed to be of fec->GetOrder(). Once
421 SetElementOrder is called, the space becomes a variable order space. */
422 void SetElementOrder(int i, int p);
423
424 /// Returns the order of the i'th finite element.
425 int GetElementOrder(int i) const;
426
427 /// Return the maximum polynomial order.
GetMaxElementOrder() const428 int GetMaxElementOrder() const
429 { return IsVariableOrder() ? elem_order.Max() : fec->GetOrder(); }
430
431 /// Returns true if the space contains elements of varying polynomial orders.
IsVariableOrder() const432 bool IsVariableOrder() const { return elem_order.Size(); }
433
434 /// The returned SparseMatrix is owned by the FiniteElementSpace.
435 const SparseMatrix *GetConformingProlongation() const;
436
437 /// The returned SparseMatrix is owned by the FiniteElementSpace.
438 const SparseMatrix *GetConformingRestriction() const;
439
440 /** Return a version of the conforming restriction matrix for variable-order
441 spaces with complex hp interfaces, where some true DOFs are not owned by
442 any elements and need to be interpolated from higher order edge/face
443 variants (see also @a SetRelaxedHpConformity()). */
444 /// The returned SparseMatrix is owned by the FiniteElementSpace.
445 const SparseMatrix *GetHpConformingRestriction() const;
446
447 /// The returned Operator is owned by the FiniteElementSpace.
GetProlongationMatrix() const448 virtual const Operator *GetProlongationMatrix() const
449 { return GetConformingProlongation(); }
450
451 /// Return an operator that performs the transpose of GetRestrictionOperator
452 /** The returned operator is owned by the FiniteElementSpace. In serial this
453 is the same as GetProlongationMatrix() */
GetRestrictionTransposeOperator() const454 virtual const Operator *GetRestrictionTransposeOperator() const
455 { return GetConformingProlongation(); }
456
457 /// An abstract operator that performs the same action as GetRestrictionMatrix
458 /** In some cases this is an optimized matrix-free implementation. The
459 returned operator is owned by the FiniteElementSpace. */
GetRestrictionOperator() const460 virtual const Operator *GetRestrictionOperator() const
461 { return GetConformingRestriction(); }
462
463 /// The returned SparseMatrix is owned by the FiniteElementSpace.
GetRestrictionMatrix() const464 virtual const SparseMatrix *GetRestrictionMatrix() const
465 { return GetConformingRestriction(); }
466
467 /// The returned SparseMatrix is owned by the FiniteElementSpace.
GetHpRestrictionMatrix() const468 virtual const SparseMatrix *GetHpRestrictionMatrix() const
469 { return GetHpConformingRestriction(); }
470
471 /// Return an Operator that converts L-vectors to E-vectors.
472 /** An L-vector is a vector of size GetVSize() which is the same size as a
473 GridFunction. An E-vector represents the element-wise discontinuous
474 version of the FE space.
475
476 The layout of the E-vector is: ND x VDIM x NE, where ND is the number of
477 degrees of freedom, VDIM is the vector dimension of the FE space, and NE
478 is the number of the mesh elements.
479
480 The parameter @a e_ordering describes how the local DOFs in each element
481 should be ordered, see ElementDofOrdering.
482
483 For discontinuous spaces, the element restriction corresponds to a
484 permutation of the degrees of freedom, implemented by the
485 L2ElementRestriction class.
486
487 The returned Operator is owned by the FiniteElementSpace. */
488 const Operator *GetElementRestriction(ElementDofOrdering e_ordering) const;
489
490 /// Return an Operator that converts L-vectors to E-vectors on each face.
491 virtual const FaceRestriction *GetFaceRestriction(
492 ElementDofOrdering e_ordering, FaceType,
493 L2FaceValues mul = L2FaceValues::DoubleValued) const;
494
495 /** @brief Return a QuadratureInterpolator that interpolates E-vectors to
496 quadrature point values and/or derivatives (Q-vectors). */
497 /** An E-vector represents the element-wise discontinuous version of the FE
498 space and can be obtained, for example, from a GridFunction using the
499 Operator returned by GetElementRestriction().
500
501 All elements will use the same IntegrationRule, @a ir as the target
502 quadrature points. */
503 const QuadratureInterpolator *GetQuadratureInterpolator(
504 const IntegrationRule &ir) const;
505
506 /** @brief Return a QuadratureInterpolator that interpolates E-vectors to
507 quadrature point values and/or derivatives (Q-vectors). */
508 /** An E-vector represents the element-wise discontinuous version of the FE
509 space and can be obtained, for example, from a GridFunction using the
510 Operator returned by GetElementRestriction().
511
512 The target quadrature points in the elements are described by the given
513 QuadratureSpace, @a qs. */
514 const QuadratureInterpolator *GetQuadratureInterpolator(
515 const QuadratureSpace &qs) const;
516
517 /** @brief Return a FaceQuadratureInterpolator that interpolates E-vectors to
518 quadrature point values and/or derivatives (Q-vectors). */
519 const FaceQuadratureInterpolator *GetFaceQuadratureInterpolator(
520 const IntegrationRule &ir, FaceType type) const;
521
522 /// Returns the polynomial degree of the i'th finite element.
523 /** NOTE: it is recommended to use GetElementOrder in new code. */
GetOrder(int i) const524 int GetOrder(int i) const { return GetElementOrder(i); }
525
526 /** Return the order of an edge. In a variable order space, return the order
527 of a specific variant, or -1 if there are no more variants. */
528 int GetEdgeOrder(int edge, int variant = 0) const;
529
530 /// Returns the polynomial degree of the i'th face finite element
531 int GetFaceOrder(int face, int variant = 0) const;
532
533 /// Returns vector dimension.
GetVDim() const534 inline int GetVDim() const { return vdim; }
535
536 /// Returns number of degrees of freedom.
GetNDofs() const537 inline int GetNDofs() const { return ndofs; }
538
539 /// Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
GetVSize() const540 inline int GetVSize() const { return vdim * ndofs; }
541
542 /// Return the number of vector true (conforming) dofs.
GetTrueVSize() const543 virtual int GetTrueVSize() const { return GetConformingVSize(); }
544
545 /// Returns the number of conforming ("true") degrees of freedom
546 /// (if the space is on a nonconforming mesh with hanging nodes).
547 int GetNConformingDofs() const;
548
GetConformingVSize() const549 int GetConformingVSize() const { return vdim * GetNConformingDofs(); }
550
551 /// Return the ordering method.
GetOrdering() const552 inline Ordering::Type GetOrdering() const { return ordering; }
553
FEColl() const554 const FiniteElementCollection *FEColl() const { return fec; }
555
556 /// Number of all scalar vertex dofs
GetNVDofs() const557 int GetNVDofs() const { return nvdofs; }
558 /// Number of all scalar edge-interior dofs
GetNEDofs() const559 int GetNEDofs() const { return nedofs; }
560 /// Number of all scalar face-interior dofs
GetNFDofs() const561 int GetNFDofs() const { return nfdofs; }
562
563 /// Returns number of vertices in the mesh.
GetNV() const564 inline int GetNV() const { return mesh->GetNV(); }
565
566 /// Returns number of elements in the mesh.
GetNE() const567 inline int GetNE() const { return mesh->GetNE(); }
568
569 /// Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
570 /** The co-dimension 1 entities are those that have dimension 1 less than the
571 mesh dimension, e.g. for a 2D mesh, the faces are the 1D entities, i.e.
572 the edges. */
GetNF() const573 inline int GetNF() const { return mesh->GetNumFaces(); }
574
575 /// Returns number of boundary elements in the mesh.
GetNBE() const576 inline int GetNBE() const { return mesh->GetNBE(); }
577
578 /// Returns the number of faces according to the requested type.
579 /** If type==Boundary returns only the "true" number of boundary faces
580 contrary to GetNBE() that returns "fake" boundary faces associated to
581 visualization for GLVis.
582 Similarly, if type==Interior, the "fake" boundary faces associated to
583 visualization are counted as interior faces. */
GetNFbyType(FaceType type) const584 inline int GetNFbyType(FaceType type) const
585 { return mesh->GetNFbyType(type); }
586
587 /// Returns the type of element i.
GetElementType(int i) const588 inline int GetElementType(int i) const
589 { return mesh->GetElementType(i); }
590
591 /// Returns the vertices of element i.
GetElementVertices(int i,Array<int> & vertices) const592 inline void GetElementVertices(int i, Array<int> &vertices) const
593 { mesh->GetElementVertices(i, vertices); }
594
595 /// Returns the type of boundary element i.
GetBdrElementType(int i) const596 inline int GetBdrElementType(int i) const
597 { return mesh->GetBdrElementType(i); }
598
599 /// Returns ElementTransformation for the @a i-th element.
GetElementTransformation(int i) const600 ElementTransformation *GetElementTransformation(int i) const
601 { return mesh->GetElementTransformation(i); }
602
603 /** @brief Returns the transformation defining the @a i-th element in the
604 user-defined variable @a ElTr. */
GetElementTransformation(int i,IsoparametricTransformation * ElTr)605 void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
606 { mesh->GetElementTransformation(i, ElTr); }
607
608 /// Returns ElementTransformation for the @a i-th boundary element.
GetBdrElementTransformation(int i) const609 ElementTransformation *GetBdrElementTransformation(int i) const
610 { return mesh->GetBdrElementTransformation(i); }
611
GetAttribute(int i) const612 int GetAttribute(int i) const { return mesh->GetAttribute(i); }
613
GetBdrAttribute(int i) const614 int GetBdrAttribute(int i) const { return mesh->GetBdrAttribute(i); }
615
616 /// Returns indices of degrees of freedom of element 'elem'.
617 virtual void GetElementDofs(int elem, Array<int> &dofs) const;
618
619 /// Returns indices of degrees of freedom for boundary element 'bel'.
620 virtual void GetBdrElementDofs(int bel, Array<int> &dofs) const;
621
622 /** @brief Returns the indices of the degrees of freedom for the specified
623 face, including the DOFs for the edges and the vertices of the face. */
624 /** In variable order spaces, multiple variants of DOFs can be returned.
625 See @a GetEdgeDofs for more details.
626 @return Order of the selected variant, or -1 if there are no more
627 variants.*/
628 virtual int GetFaceDofs(int face, Array<int> &dofs, int variant = 0) const;
629
630 /** @brief Returns the indices of the degrees of freedom for the specified
631 edge, including the DOFs for the vertices of the edge. */
632 /** In variable order spaces, multiple sets of DOFs may exist on an edge,
633 corresponding to the different polynomial orders of incident elements.
634 The 'variant' parameter is the zero-based index of the desired DOF set.
635 The variants are ordered from lowest polynomial degree to the highest.
636 @return Order of the selected variant, or -1 if there are no more
637 variants. */
638 int GetEdgeDofs(int edge, Array<int> &dofs, int variant = 0) const;
639
640 void GetVertexDofs(int i, Array<int> &dofs) const;
641
642 void GetElementInteriorDofs(int i, Array<int> &dofs) const;
643
644 void GetFaceInteriorDofs(int i, Array<int> &dofs) const;
645
646 int GetNumElementInteriorDofs(int i) const;
647
648 void GetEdgeInteriorDofs(int i, Array<int> &dofs) const;
649
650 /** @brief Returns the indices of all of the VDofs for the specified
651 dimension 'vd'. */
652 /** The 'ndofs' parameter defines the number of Dofs in the
653 FiniteElementSpace. If 'ndofs' is -1 (the default value), then the
654 number of Dofs is determined by the FiniteElementSpace. */
655 void GetVDofs(int vd, Array<int> &dofs, int ndofs = -1) const;
656
657 void DofsToVDofs(Array<int> &dofs, int ndofs = -1) const;
658
659 void DofsToVDofs(int vd, Array<int> &dofs, int ndofs = -1) const;
660
661 int DofToVDof(int dof, int vd, int ndofs = -1) const;
662
VDofToDof(int vdof) const663 int VDofToDof(int vdof) const
664 { return (ordering == Ordering::byNODES) ? (vdof%ndofs) : (vdof/vdim); }
665
666 static void AdjustVDofs(Array<int> &vdofs);
667
668 /// Returns indexes of degrees of freedom in array dofs for i'th element.
669 void GetElementVDofs(int i, Array<int> &vdofs) const;
670
671 /// Returns indexes of degrees of freedom for i'th boundary element.
672 void GetBdrElementVDofs(int i, Array<int> &vdofs) const;
673
674 /// Returns indexes of degrees of freedom for i'th face element (2D and 3D).
675 void GetFaceVDofs(int i, Array<int> &vdofs) const;
676
677 /// Returns indexes of degrees of freedom for i'th edge.
678 void GetEdgeVDofs(int i, Array<int> &vdofs) const;
679
680 void GetVertexVDofs(int i, Array<int> &vdofs) const;
681
682 void GetElementInteriorVDofs(int i, Array<int> &vdofs) const;
683
684 void GetEdgeInteriorVDofs(int i, Array<int> &vdofs) const;
685
686 /// (@deprecated) Use the Update() method if the space or mesh changed.
687 MFEM_DEPRECATED void RebuildElementToDofTable();
688
689 /** @brief Reorder the scalar DOFs based on the element ordering.
690
691 The new ordering is constructed as follows: 1) loop over all elements as
692 ordered in the Mesh; 2) for each element, assign new indices to all of
693 its current DOFs that are still unassigned; the new indices we assign are
694 simply the sequence `0,1,2,...`; if there are any signed DOFs their sign
695 is preserved. */
696 void ReorderElementToDofTable();
697
698 /** @brief Return a reference to the internal Table that stores the lists of
699 scalar dofs, for each mesh element, as returned by GetElementDofs(). */
GetElementToDofTable() const700 const Table &GetElementToDofTable() const { return *elem_dof; }
701
702 /** @brief Return a reference to the internal Table that stores the lists of
703 scalar dofs, for each boundary mesh element, as returned by
704 GetBdrElementDofs(). */
GetBdrElementToDofTable() const705 const Table &GetBdrElementToDofTable() const
706 { if (!bdr_elem_dof) { BuildBdrElementToDofTable(); } return *bdr_elem_dof; }
707
708 /** @brief Return a reference to the internal Table that stores the lists of
709 scalar dofs, for each face in the mesh, as returned by GetFaceDofs(). In
710 this context, "face" refers to a (dim-1)-dimensional mesh entity. */
711 /** @note In the case of a NURBS space, the rows corresponding to interior
712 faces will be empty. */
GetFaceToDofTable() const713 const Table &GetFaceToDofTable() const
714 { if (!face_dof) { BuildFaceToDofTable(); } return *face_dof; }
715
716 /** @brief Initialize internal data that enables the use of the methods
717 GetElementForDof() and GetLocalDofForDof(). */
718 void BuildDofToArrays();
719
720 /// Return the index of the first element that contains dof @a i.
721 /** This method can be called only after setup is performed using the method
722 BuildDofToArrays(). */
GetElementForDof(int i) const723 int GetElementForDof(int i) const { return dof_elem_array[i]; }
724 /// Return the local dof index in the first element that contains dof @a i.
725 /** This method can be called only after setup is performed using the method
726 BuildDofToArrays(). */
GetLocalDofForDof(int i) const727 int GetLocalDofForDof(int i) const { return dof_ldof_array[i]; }
728
729 /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
730 associated with i'th element in the mesh object. */
731 virtual const FiniteElement *GetFE(int i) const;
732
733 /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
734 associated with i'th boundary face in the mesh object. */
735 const FiniteElement *GetBE(int i) const;
736
737 /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
738 associated with i'th face in the mesh object. Faces in this case refer
739 to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are
740 points.*/
741 const FiniteElement *GetFaceElement(int i) const;
742
743 /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
744 associated with i'th edge in the mesh object. */
745 const FiniteElement *GetEdgeElement(int i, int variant = 0) const;
746
747 /// Return the trace element from element 'i' to the given 'geom_type'
748 const FiniteElement *GetTraceElement(int i, Geometry::Type geom_type) const;
749
750 /** @brief Mark degrees of freedom associated with boundary elements with
751 the specified boundary attributes (marked in 'bdr_attr_is_ess').
752 For spaces with 'vdim' > 1, the 'component' parameter can be used
753 to restricts the marked vDOFs to the specified component. */
754 virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
755 Array<int> &ess_vdofs,
756 int component = -1) const;
757
758 /** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
759 boundary attributes marked in the array bdr_attr_is_ess.
760 For spaces with 'vdim' > 1, the 'component' parameter can be used
761 to restricts the marked tDOFs to the specified component. */
762 virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
763 Array<int> &ess_tdof_list,
764 int component = -1);
765
766 /** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
767 with 'vdim' > 1, the 'component' parameter can be used to restricts the
768 marked tDOFs to the specified component. Equivalent to
769 FiniteElementSpace::GetEssentialTrueDofs with all boundary attributes
770 marked as essential. */
771 void GetBoundaryTrueDofs(Array<int> &boundary_dofs, int component = -1);
772
773 /// Convert a Boolean marker array to a list containing all marked indices.
774 static void MarkerToList(const Array<int> &marker, Array<int> &list);
775
776 /** @brief Convert an array of indices (list) to a Boolean marker array where all
777 indices in the list are marked with the given value and the rest are set
778 to zero. */
779 static void ListToMarker(const Array<int> &list, int marker_size,
780 Array<int> &marker, int mark_val = -1);
781
782 /** @brief For a partially conforming FE space, convert a marker array (nonzero
783 entries are true) on the partially conforming dofs to a marker array on
784 the conforming dofs. A conforming dofs is marked iff at least one of its
785 dependent dofs is marked. */
786 void ConvertToConformingVDofs(const Array<int> &dofs, Array<int> &cdofs);
787
788 /** @brief For a partially conforming FE space, convert a marker array (nonzero
789 entries are true) on the conforming dofs to a marker array on the
790 (partially conforming) dofs. A dof is marked iff it depends on a marked
791 conforming dofs, where dependency is defined by the ConformingRestriction
792 matrix; in other words, a dof is marked iff it corresponds to a marked
793 conforming dof. */
794 void ConvertFromConformingVDofs(const Array<int> &cdofs, Array<int> &dofs);
795
796 /** @brief Generate the global restriction matrix from a discontinuous
797 FE space to the continuous FE space of the same polynomial degree. */
798 SparseMatrix *D2C_GlobalRestrictionMatrix(FiniteElementSpace *cfes);
799
800 /** @brief Generate the global restriction matrix from a discontinuous
801 FE space to the piecewise constant FE space. */
802 SparseMatrix *D2Const_GlobalRestrictionMatrix(FiniteElementSpace *cfes);
803
804 /** @brief Construct the restriction matrix from the FE space given by
805 (*this) to the lower degree FE space given by (*lfes) which
806 is defined on the same mesh. */
807 SparseMatrix *H2L_GlobalRestrictionMatrix(FiniteElementSpace *lfes);
808
809 /** @brief Construct and return an Operator that can be used to transfer
810 GridFunction data from @a coarse_fes, defined on a coarse mesh, to @a
811 this FE space, defined on a refined mesh. */
812 /** It is assumed that the mesh of this FE space is a refinement of the mesh
813 of @a coarse_fes and the CoarseFineTransformations returned by the method
814 Mesh::GetRefinementTransforms() of the refined mesh are set accordingly.
815 The Operator::Type of @a T can be set to request an Operator of the set
816 type. Currently, only Operator::MFEM_SPARSEMAT and Operator::ANY_TYPE
817 (matrix-free) are supported. When Operator::ANY_TYPE is requested, the
818 choice of the particular Operator sub-class is left to the method. This
819 method also works in parallel because the transfer operator is local to
820 the MPI task when the input is a synchronized ParGridFunction. */
821 void GetTransferOperator(const FiniteElementSpace &coarse_fes,
822 OperatorHandle &T) const;
823
824 /** @brief Construct and return an Operator that can be used to transfer
825 true-dof data from @a coarse_fes, defined on a coarse mesh, to @a this FE
826 space, defined on a refined mesh.
827
828 This method calls GetTransferOperator() and multiplies the result by the
829 prolongation operator of @a coarse_fes on the right, and by the
830 restriction operator of this FE space on the left.
831
832 The Operator::Type of @a T can be set to request an Operator of the set
833 type. In serial, the supported types are: Operator::MFEM_SPARSEMAT and
834 Operator::ANY_TYPE (matrix-free). In parallel, the supported types are:
835 Operator::Hypre_ParCSR and Operator::ANY_TYPE. Any other type is treated
836 as Operator::ANY_TYPE: the operator representation choice is made by this
837 method. */
838 virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes,
839 OperatorHandle &T) const;
840
841 /** @brief Reflect changes in the mesh: update number of DOFs, etc. Also, calculate
842 GridFunction transformation operator (unless want_transform is false).
843 Safe to call multiple times, does nothing if space already up to date. */
844 virtual void Update(bool want_transform = true);
845
846 /// Get the GridFunction update operator.
GetUpdateOperator()847 const Operator* GetUpdateOperator() { Update(); return Th.Ptr(); }
848
849 /// Return the update operator in the given OperatorHandle, @a T.
GetUpdateOperator(OperatorHandle & T)850 void GetUpdateOperator(OperatorHandle &T) { T = Th; }
851
852 /** @brief Set the ownership of the update operator: if set to false, the
853 Operator returned by GetUpdateOperator() must be deleted outside the
854 FiniteElementSpace. */
855 /** The update operator ownership is automatically reset to true when a new
856 update operator is created by the Update() method. */
SetUpdateOperatorOwner(bool own)857 void SetUpdateOperatorOwner(bool own) { Th.SetOperatorOwner(own); }
858
859 /// Specify the Operator::Type to be used by the update operators.
860 /** The default type is Operator::ANY_TYPE which leaves the choice to this
861 class. The other currently supported option is Operator::MFEM_SPARSEMAT
862 which is only guaranteed to be honored for a refinement update operator.
863 Any other type will be treated as Operator::ANY_TYPE.
864 @note This operation destroys the current update operator (if owned). */
SetUpdateOperatorType(Operator::Type tid)865 void SetUpdateOperatorType(Operator::Type tid) { Th.SetType(tid); }
866
867 /// Free the GridFunction update operator (if any), to save memory.
UpdatesFinished()868 virtual void UpdatesFinished() { Th.Clear(); }
869
870 /** Return update counter, similar to Mesh::GetSequence(). Used by
871 GridFunction to check if it is up to date with the space. */
GetSequence() const872 long GetSequence() const { return sequence; }
873
874 /// Return whether or not the space is discontinuous (L2)
IsDGSpace() const875 bool IsDGSpace() const
876 {
877 return dynamic_cast<const L2_FECollection*>(fec) != NULL;
878 }
879
880 /** In variable order spaces on nonconforming (NC) meshes, this function
881 controls whether strict conformity is enforced in cases where coarse
882 edges/faces have higher polynomial order than their fine NC neighbors.
883 In the default (strict) case, the coarse side polynomial order is
884 reduced to that of the lowest order fine edge/face, so all fine
885 neighbors can interpolate the coarse side exactly. If relaxed == true,
886 some discontinuities in the solution in such cases are allowed and the
887 coarse side is not restricted. For an example, see
888 https://github.com/mfem/mfem/pull/1423#issuecomment-621340392 */
SetRelaxedHpConformity(bool relaxed=true)889 void SetRelaxedHpConformity(bool relaxed = true)
890 {
891 relaxed_hp = relaxed;
892 orders_changed = true; // force update
893 Update(false);
894 }
895
896 /// Save finite element space to output stream @a out.
897 void Save(std::ostream &out) const;
898
899 /** @brief Read a FiniteElementSpace from a stream. The returned
900 FiniteElementCollection is owned by the caller. */
901 FiniteElementCollection *Load(Mesh *m, std::istream &input);
902
903 virtual ~FiniteElementSpace();
904 };
905
906
907 /// Class representing the storage layout of a QuadratureFunction.
908 /** Multiple QuadratureFunction%s can share the same QuadratureSpace. */
909 class QuadratureSpace
910 {
911 protected:
912 friend class QuadratureFunction; // Uses the element_offsets.
913
914 Mesh *mesh;
915 int order;
916 int size;
917
918 const IntegrationRule *int_rule[Geometry::NumGeom];
919 int *element_offsets; // scalar offsets; size = number of elements + 1
920
921 // protected functions
922
923 // Assuming mesh and order are set, construct the members: int_rule,
924 // element_offsets, and size.
925 void Construct();
926
927 public:
928 /// Create a QuadratureSpace based on the global rules from #IntRules.
QuadratureSpace(Mesh * mesh_,int order_)929 QuadratureSpace(Mesh *mesh_, int order_)
930 : mesh(mesh_), order(order_) { Construct(); }
931
932 /// Read a QuadratureSpace from the stream @a in.
933 QuadratureSpace(Mesh *mesh_, std::istream &in);
934
~QuadratureSpace()935 virtual ~QuadratureSpace() { delete [] element_offsets; }
936
937 /// Return the total number of quadrature points.
GetSize() const938 int GetSize() const { return size; }
939
940 /// Return the order of the quadrature rule(s) used by all elements.
GetOrder() const941 int GetOrder() const { return order; }
942
943 /// Returns the mesh
GetMesh() const944 inline Mesh *GetMesh() const { return mesh; }
945
946 /// Returns number of elements in the mesh.
GetNE() const947 inline int GetNE() const { return mesh->GetNE(); }
948
949 /// Get the IntegrationRule associated with mesh element @a idx.
GetElementIntRule(int idx) const950 const IntegrationRule &GetElementIntRule(int idx) const
951 { return *int_rule[mesh->GetElementBaseGeometry(idx)]; }
952
953 /// Write the QuadratureSpace to the stream @a out.
954 void Save(std::ostream &out) const;
955 };
956
UsesTensorBasis(const FiniteElementSpace & fes)957 inline bool UsesTensorBasis(const FiniteElementSpace& fes)
958 {
959 // TODO: mixed meshes: return true if there is at least one tensor-product
960 // Geometry in the global mesh and the FE collection returns a
961 // TensorBasisElement for that Geometry?
962
963 // Potential issue: empty local mesh --> no element 0.
964 return dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
965 }
966
967 }
968
969 #endif
970