MFEM  v4.4.0 Finite element discretization library
restriction.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
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_RESTRICTION
13 #define MFEM_RESTRICTION
14
15 #include "../linalg/operator.hpp"
16 #include "../mesh/mesh.hpp"
17
18 namespace mfem
19 {
20
21 class FiniteElementSpace;
22 enum class ElementDofOrdering;
23
24 /// Operator that converts FiniteElementSpace L-vectors to E-vectors.
25 /** Objects of this type are typically created and owned by FiniteElementSpace
26  objects, see FiniteElementSpace::GetElementRestriction(). */
28 {
29 private:
30  /** This number defines the maximum number of elements any dof can belong to
31  for the FillSparseMatrix method. */
32  static const int MaxNbNbr = 16;
33
34 protected:
36  const int ne;
37  const int vdim;
38  const bool byvdim;
39  const int ndofs;
40  const int dof;
41  const int nedofs;
45
46 public:
48  void Mult(const Vector &x, Vector &y) const;
49  void MultTranspose(const Vector &x, Vector &y) const;
50
51  /// Compute Mult without applying signs based on DOF orientations.
52  void MultUnsigned(const Vector &x, Vector &y) const;
53  /// Compute MultTranspose without applying signs based on DOF orientations.
54  void MultTransposeUnsigned(const Vector &x, Vector &y) const;
55
56  /// Compute MultTranspose by setting (rather than adding) element
57  /// contributions; this is a left inverse of the Mult() operation
58  void MultLeftInverse(const Vector &x, Vector &y) const;
59
60  /// @brief Fills the E-vector y with boolean values 0.0 and 1.0 such that each
61  /// each entry of the L-vector is uniquely represented in y.
62  /** This means, the sum of the E-vector y is equal to the sum of the
63  corresponding L-vector filled with ones. The boolean mask is required to
64  emulate SetSubVector and its transpose on GPUs. This method is running on
65  the host, since the processed array requires a large shared memory. */
67
68  /// Fill a Sparse Matrix with Element Matrices.
69  void FillSparseMatrix(const Vector &mat_ea, SparseMatrix &mat) const;
70
71  /** Fill the I array of SparseMatrix corresponding to the sparsity pattern
72  given by this ElementRestriction. */
73  int FillI(SparseMatrix &mat) const;
74  /** Fill the J and Data arrays of SparseMatrix corresponding to the sparsity
75  pattern given by this ElementRestriction, and the values of ea_data. */
76  void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const;
77 };
78
79 /// Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
80 /** Objects of this type are typically created and owned by FiniteElementSpace
81  objects, see FiniteElementSpace::GetElementRestriction(). L-vectors
82  corresponding to grid functions in L2 finite element spaces differ from
83  E-vectors only in the ordering of the degrees of freedom. */
85 {
86  const int ne;
87  const int vdim;
88  const bool byvdim;
89  const int ndof;
90  const int ndofs;
91 public:
93  void Mult(const Vector &x, Vector &y) const;
94  void MultTranspose(const Vector &x, Vector &y) const;
95  /** Fill the I array of SparseMatrix corresponding to the sparsity pattern
96  given by this ElementRestriction. */
97  void FillI(SparseMatrix &mat) const;
98  /** Fill the J and Data arrays of SparseMatrix corresponding to the sparsity
99  pattern given by this L2FaceRestriction, and the values of ea_data. */
100  void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const;
101 };
102
103 /** An enum type to specify if only e1 value is requested (SingleValued) or both
104  e1 and e2 (DoubleValued). */
106
107 /** @brief Base class for operators that extracts Face degrees of freedom.
108
109  In order to compute quantities on the faces of a mesh, it is often useful to
110  extract the degrees of freedom on the faces of the elements. This class
111  provides an interface for such operations.
112
113  If the FiniteElementSpace is ordered by Ordering::byVDIM, then the expected
114  format for the L-vector is (vdim x ndofs), otherwise if Ordering::byNODES
115  the expected format is (ndofs x vdim), where ndofs is the total number of
116  degrees of freedom.
117  Since FiniteElementSpace can either be continuous or discontinuous, the
118  degrees of freedom on a face can either be single valued or double valued,
119  this is what we refer to as the multiplicity and is represented by the
120  L2FaceValues enum type.
121  The format of the output face E-vector of degrees of freedom is
122  (face_dofs x vdim x multiplicity x nfaces), where face_dofs is the number of
123  degrees of freedom on each face, and nfaces the number of faces of the
124  requested FaceType (see FiniteElementSpace::GetNFbyType).
125
126  @note Objects of this type are typically created and owned by
127  FiniteElementSpace objects, see FiniteElementSpace::GetFaceRestriction(). */
128 class FaceRestriction : public Operator
129 {
130 public:
132
133  FaceRestriction(int h, int w): Operator(h, w) { }
134
135  virtual ~FaceRestriction() { }
136
137  /** @brief Extract the face degrees of freedom from @a x into @a y.
138
139  @param[in] x The L-vector of degrees of freedom.
140  @param[out] y The degrees of freedom on the face, corresponding to a face
141  E-vector.
142  */
143  void Mult(const Vector &x, Vector &y) const override = 0;
144
145  /** @brief Add the face degrees of freedom @a x to the element degrees of
146  freedom @a y.
147
148  @param[in] x The face degrees of freedom on the face.
149  @param[in,out] y The L-vector of degrees of freedom to which we add the
150  face degrees of freedom.
151  */
152  virtual void AddMultTranspose(const Vector &x, Vector &y) const = 0;
153
154  /** @brief Set the face degrees of freedom in the element degrees of freedom
155  @a y to the values given in @a x.
156
157  @param[in] x The face degrees of freedom on the face.
158  @param[in,out] y The L-vector of degrees of freedom to which we add the
159  face degrees of freedom.
160  */
161  void MultTranspose(const Vector &x, Vector &y) const override
162  {
163  y = 0.0;
165  }
166 };
167
168 /// Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
169 /** Objects of this type are typically created and owned by FiniteElementSpace
170  objects, see FiniteElementSpace::GetFaceRestriction(). */
172 {
173 protected:
175  const int nf; // Number of faces of the requested type
176  const int vdim;
177  const bool byvdim;
178  const int face_dofs; // Number of dofs on each face
179  const int elem_dofs; // Number of dofs in each element
180  const int nfdofs; // Total number of face E-vector dofs
181  const int ndofs; // Total number of dofs
182  Array<int> scatter_indices; // Scattering indices for element 1 on each face
183  Array<int> gather_offsets; // offsets for the gathering indices of each dof
184  Array<int> gather_indices; // gathering indices for each dof
185
186  /** @brief Construct an H1FaceRestriction.
187
188  @param[in] fes The FiniteElementSpace on which this operates
189  @param[in] ordering Request a specific element ordering
190  @param[in] type Request internal or boundary faces dofs
191  @param[in] build Request the NCL2FaceRestriction to compute the
192  scatter/gather indices. False should only be used
193  when inheriting from H1FaceRestriction.
194  */
196  const ElementDofOrdering ordering,
197  const FaceType type,
198  bool build);
199 public:
200  /** @brief Construct an H1FaceRestriction.
201
202  @param[in] fes The FiniteElementSpace on which this operates
203  @param[in] ordering Request a specific element ordering
204  @param[in] type Request internal or boundary faces dofs */
206  const ElementDofOrdering ordering,
207  const FaceType type);
208
209  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
210  face E-Vector.
211
212  @param[in] x The L-vector degrees of freedom.
213  @param[out] y The face E-Vector degrees of freedom with the given format:
214  face_dofs x vdim x nf
215  where nf is the number of interior or boundary faces
216  requested by @a type in the constructor.
217  The face_dofs are ordered according to the given
218  ElementDofOrdering. */
219  void Mult(const Vector &x, Vector &y) const override;
220
221  /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
222  L-Vector.
223
224  @param[in] x The face E-Vector degrees of freedom with the given format:
225  face_dofs x vdim x nf
226  where nf is the number of interior or boundary faces
227  requested by @a type in the constructor.
228  The face_dofs should be ordered according to the given
229  ElementDofOrdering
230  @param[in,out] y The L-vector degrees of freedom. */
231  void AddMultTranspose(const Vector &x, Vector &y) const override;
232
233 private:
234  /** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
235  for the gathering: E-vector to L-vector.
236
237  @param[in] ordering Request a specific element ordering.
238  @param[in] type Request internal or boundary faces dofs.
239  */
240  void ComputeScatterIndicesAndOffsets(const ElementDofOrdering ordering,
241  const FaceType type);
242
243  /** @brief Compute the gather indices: E-vector to L-vector.
244
245  Note: Requires the gather offsets to be computed.
246
247  @param[in] ordering Request a specific element ordering.
248  @param[in] type Request internal or boundary faces dofs.
249  */
250  void ComputeGatherIndices(const ElementDofOrdering ordering,
251  const FaceType type);
252
253 protected:
254  mutable Array<int> face_map; // Used in the computation of GetFaceDofs
255
256  /** @brief Verify that H1FaceRestriction is build from an H1 FESpace.
257
258  @param[in] ordering The FESpace element ordering.
259  */
260  void CheckFESpace(const ElementDofOrdering ordering);
261
262  /** @brief Set the scattering indices of elem1, and increment the offsets for
263  the face described by the @a face.
264
265  @param[in] face The face information of the current face.
266  @param[in] face_index The interior/boundary face index.
267  @param[in] ordering Request a specific element ordering.
268  */
270  const int face_index,
271  const ElementDofOrdering ordering);
272
273  /** @brief Set the gathering indices of elem1 for the interior face described
274  by the @a face.
275
276  @param[in] face The face information of the current face.
277  @param[in] face_index The interior/boundary face index.
278  @param[in] ordering Request a specific element ordering.
279  */
281  const int face_index,
282  const ElementDofOrdering ordering);
283 };
284
285 /// Operator that extracts Face degrees of freedom for L2 spaces.
286 /** Objects of this type are typically created and owned by FiniteElementSpace
287  objects, see FiniteElementSpace::GetFaceRestriction(). */
289 {
290 protected:
292  const int nf; // Number of faces of the requested type
293  const int ne; // Number of elements
294  const int vdim; // vdim
295  const bool byvdim;
296  const int face_dofs; // Number of dofs on each face
297  const int elem_dofs; // Number of dofs in each element
298  const int nfdofs; // Total number of dofs on the faces
299  const int ndofs; // Total number of dofs
300  const FaceType type;
302  Array<int> scatter_indices1; // Scattering indices for element 1 on each face
303  Array<int> scatter_indices2; // Scattering indices for element 2 on each face
304  Array<int> gather_offsets; // offsets for the gathering indices of each dof
305  Array<int> gather_indices; // gathering indices for each dof
306
307  /** @brief Constructs an L2FaceRestriction.
308
309  @param[in] fes The FiniteElementSpace on which this operates
310  @param[in] ordering Request a specific ordering
311  @param[in] type Request internal or boundary faces dofs
312  @param[in] m Request the face dofs for elem1, or both elem1 and
313  elem2
314  @param[in] build Request the NCL2FaceRestriction to compute the
315  scatter/gather indices. False should only be used
316  when inheriting from L2FaceRestriction.
317  */
319  const ElementDofOrdering ordering,
320  const FaceType type,
321  const L2FaceValues m,
322  bool build);
323
324 public:
325  /** @brief Constructs an L2FaceRestriction.
326
327  @param[in] fes The FiniteElementSpace on which this operates
328  @param[in] ordering Request a specific ordering
329  @param[in] type Request internal or boundary faces dofs
330  @param[in] m Request the face dofs for elem1, or both elem1 and
331  elem2 */
333  const ElementDofOrdering ordering,
334  const FaceType type,
336
337  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
338  face E-Vector.
339
340  @param[in] x The L-vector degrees of freedom.
341  @param[out] y The face E-Vector degrees of freedom with the given format:
342  if L2FacesValues::DoubleValued (face_dofs x vdim x 2 x nf)
343  if L2FacesValues::SingleValued (face_dofs x vdim x nf)
344  where nf is the number of interior or boundary faces
345  requested by @a type in the constructor.
346  The face_dofs are ordered according to the given
347  ElementDofOrdering. */
348  void Mult(const Vector &x, Vector &y) const override;
349
350  /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
351  L-Vector.
352
353  @param[in] x The face E-Vector degrees of freedom with the given format:
354  if L2FacesValues::DoubleValued (face_dofs x vdim x 2 x nf)
355  if L2FacesValues::SingleValued (face_dofs x vdim x nf)
356  where nf is the number of interior or boundary faces
357  requested by @a type in the constructor.
358  The face_dofs should be ordered according to the given
359  ElementDofOrdering
360  @param[in,out] y The L-vector degrees of freedom. */
361  void AddMultTranspose(const Vector &x, Vector &y) const override;
362
363  /** @brief Fill the I array of SparseMatrix corresponding to the sparsity
364  pattern given by this L2FaceRestriction.
365
366  @param[in,out] mat The sparse matrix for which we want to initialize the
367  row offsets.
368  @param[in] keep_nbr_block When set to true the SparseMatrix will
369  include the rows (in addition to the columns)
370  corresponding to face-neighbor dofs. The
371  default behavior is to disregard those rows. */
372  virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block = false) const;
373
374  /** @brief Fill the J and Data arrays of the SparseMatrix corresponding to
375  the sparsity pattern given by this L2FaceRestriction, and the values of
376  fea_data.
377
378  @param[in] fea_data The dense matrices representing the local operators
379  on each face. The format is:
380  face_dofs x face_dofs x 2 x nf
381  On each face the first local matrix corresponds to
382  the contribution of elem1 on elem2, and the second to
383  the contribution of elem2 on elem1.
384  @param[in,out] mat The sparse matrix that is getting filled.
385  @param[in] keep_nbr_block When set to true the SparseMatrix will
386  include the rows (in addition to the columns)
387  corresponding to face-neighbor dofs. The
388  default behavior is to disregard those rows. */
389  virtual void FillJAndData(const Vector &fea_data,
390  SparseMatrix &mat,
391  const bool keep_nbr_block = false) const;
392
393  /** @brief This methods adds the DG face matrices to the element matrices.
394
395  @param[in] fea_data The dense matrices representing the local operators
396  on each face. The format is:
397  face_dofs x face_dofs x 2 x nf
398  On each face the first and second local matrices
399  correspond to the contributions of elem1 and elem2 on
400  themselves respectively.
401  @param[in,out] ea_data The dense matrices representing the element local
402  contributions for each element to which will be
404  The format is: dofs x dofs x ne, where dofs is the
405  number of dofs per element and ne the number of
406  elements. */
407  virtual void AddFaceMatricesToElementMatrices(const Vector &fea_data,
408  Vector &ea_data) const;
409
410 private:
411  /** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
412  for the gathering: E-vector to L-vector.
413
414  @param[in] ordering Request a specific element ordering.
415  @param[in] type Request internal or boundary faces dofs.
416  */
417  void ComputeScatterIndicesAndOffsets(const ElementDofOrdering ordering,
418  const FaceType type);
419
420  /** @brief Compute the gather indices: E-vector to L-vector.
421
422  Note: Requires the gather offsets to be computed.
423
424  @param[in] ordering Request a specific element ordering.
425  @param[in] type Request internal or boundary faces dofs.
426  */
427  void ComputeGatherIndices(const ElementDofOrdering ordering,
428  const FaceType type);
429
430 protected:
431  mutable Array<int> face_map; // Used in the computation of GetFaceDofs
432
433  /** @brief Verify that L2FaceRestriction is build from an L2 FESpace.
434
435  @param[in] ordering The FESpace element ordering.
436  */
437  void CheckFESpace(const ElementDofOrdering ordering);
438
439  /** @brief Set the scattering indices of elem1, and increment the offsets for
440  the face described by the @a face. The ordering of the face dofs of elem1
441  is lexicographic relative to elem1.
442
443  @param[in] face The face information of the current face.
444  @param[in] face_index The interior/boundary face index.
445  */
447  const int face_index);
448
449  /** @brief Permute and set the scattering indices of elem2, and increment the
450  offsets for the face described by the @a face. The permutation orders the
451  dofs of elem2 lexicographically as the ones of elem1.
452
453  @param[in] face The face information of the current face.
454  @param[in] face_index The interior/boundary face index.
455  */
457  const int face_index);
458
459  /** @brief Permute and set the scattering indices of elem2 for the shared
460  face described by the @a face. The permutation orders the dofs of elem2 as
461  the ones of elem1.
462
463  @param[in] face The face information of the current face.
464  @param[in] face_index The interior/boundary face index.
465  */
467  const Mesh::FaceInformation &face,
468  const int face_index);
469
470  /** @brief Set the scattering indices of elem2 for the boundary face
471  described by the @a face.
472
473  @param[in] face The face information of the current face.
474  @param[in] face_index The interior/boundary face index.
475  */
477  const int face_index);
478
479  /** @brief Set the gathering indices of elem1 for the interior face described
480  by the @a face.
481
482  Note: This function modifies the offsets.
483
484  @param[in] face The face information of the current face.
485  @param[in] face_index The interior/boundary face index.
486  */
488  const int face_index);
489
490  /** @brief Permute and set the gathering indices of elem2 for the interior
491  face described by the @a face. The permutation orders the dofs of elem2 as
492  the ones of elem1.
493
494  Note: This function modifies the offsets.
495
496  @param[in] face The face information of the current face.
497  @param[in] face_index The interior/boundary face index.
498  */
500  const int face_index);
501
502 public:
503  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
504  face E-Vector. Should only be used with conforming faces and when:
505  m == L2FacesValues::SingleValued
506
507  @param[in] x The L-vector degrees of freedom.
508  @param[out] y The face E-Vector degrees of freedom with the given format:
509  face_dofs x vdim x nf
510  where nf is the number of interior or boundary faces
511  requested by @a type in the constructor.
512  The face_dofs are ordered according to the given
513  ElementDofOrdering. */
514  void SingleValuedConformingMult(const Vector& x, Vector& y) const;
515
516  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
517  face E-Vector. Should only be used with conforming faces and when:
518  m == L2FacesValues::DoubleValued
519
520  @param[in] x The L-vector degrees of freedom.
521  @param[out] y The face E-Vector degrees of freedom with the given format:
522  face_dofs x vdim x 2 x nf
523  where nf is the number of interior or boundary faces
524  requested by @a type in the constructor.
525  The face_dofs are ordered according to the given
526  ElementDofOrdering. */
527  virtual void DoubleValuedConformingMult(const Vector& x, Vector& y) const;
528
529  /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
530  L-Vector. Should only be used with conforming faces and when:
531  m == L2FacesValues::SingleValued
532
533  @param[in] x The face E-Vector degrees of freedom with the given format:
534  face_dofs x vdim x nf
535  where nf is the number of interior or boundary faces
536  requested by @a type in the constructor.
537  The face_dofs should be ordered according to the given
538  ElementDofOrdering
539  @param[in,out] y The L-vector degrees of freedom. */
540  void SingleValuedConformingAddMultTranspose(const Vector& x, Vector& y) const;
541
542  /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
543  L-Vector. Should only be used with conforming faces and when:
544  m == L2FacesValues::DoubleValued
545
546  @param[in] x The face E-Vector degrees of freedom with the given format:
547  face_dofs x vdim x 2 x nf
548  where nf is the number of interior or boundary faces
549  requested by @a type in the constructor.
550  The face_dofs should be ordered according to the given
551  ElementDofOrdering
552  @param[in,out] y The L-vector degrees of freedom. */
553  void DoubleValuedConformingAddMultTranspose(const Vector& x, Vector& y) const;
554 };
555
556 /** This struct stores which side is the master nonconforming side and the
557  index of the interpolator, see InterpolationManager class below. */
559 {
560  uint32_t is_non_conforming : 1;
561  uint32_t master_side : 1;
562  uint32_t index : 30;
563
564  // default constructor, create a conforming face with index 0.
565  InterpConfig() = default;
566
567  // Non-conforming face
568  InterpConfig(int master_side, int nc_index)
569  : is_non_conforming(1), master_side(master_side), index(nc_index)
570  { }
571
572  InterpConfig(const InterpConfig&) = default;
573
574  InterpConfig &operator=(const InterpConfig &rhs) = default;
575 };
576
577 /** @brief This class manages the storage and computation of the interpolations
578  from master (coarse) face to slave (fine) face.
579 */
581 {
582 protected:
585  Array<InterpConfig> interp_config; // interpolator index for each face
586  Vector interpolators; // face_dofs x face_dofs x num_interpolators
587  int nc_cpt; // Counter for interpolators, and used as index.
588
589  /** The interpolators are associated to a key of containing the address of
590  PointMatrix and a local face identifier. */
591  using Key = std::pair<const DenseMatrix*,int>;
592  /// The temporary map used to store the different interpolators.
593  using Map = std::map<Key, std::pair<int,const DenseMatrix*>>;
594  Map interp_map; // The temporary map that stores the interpolators.
595
596 public:
597  InterpolationManager() = delete;
598
599  /** @brief main constructor.
600
601  @param[in] fes The FiniteElementSpace on which this operates
602  @param[in] ordering Request a specific element ordering.
603  @param[in] type Request internal or boundary faces dofs
604  */
607  FaceType type);
608
609  /** @brief Register the face with @a face and index @a face_index as a
610  conforming face for the interpolation of the degrees of freedom.
611
612  @param[in] face The face information of the current face.
613  @param[in] face_index The interior/boundary face index.
614  */
616  int face_index);
617
618  /** @brief Register the face with @a face and index @a face_index as a
619  conforming face for the interpolation of the degrees of freedom.
620
621  @param[in] face The face information of the current face.
622  @param[in] face_index The interior/boundary face index.
623  */
625  int face_index);
626
627  /** @brief Transform the interpolation matrix map into a contiguous memory
628  structure. */
630
631  /// @brief Return the total number of interpolators.
633  {
634  return nc_cpt;
635  }
636
637  /** @brief Return an mfem::Vector containing the interpolators in the
638  following format: face_dofs x face_dofs x num_interpolators. */
639  const Vector& GetInterpolators() const
640  {
641  return interpolators;
642  }
643
644  /** @brief Return an array containing the interpolation configuration for
645  each face registered with RegisterFaceConformingInterpolation and
646  RegisterFaceCoarseToFineInterpolation. */
648  {
649  return interp_config;
650  }
651
652 private:
653  /** @brief Returns the interpolation operator from a master (coarse) face to
654  a slave (fine) face.
655
656  @param[in] face The face information of the current face.
657  @param[in] ptMat The PointMatrix describing the position and orientation
658  of the fine face in the coarse face. This PointMatrix is
659  usually obtained from the mesh through the method
660  GetNCFacesPtMat.
661  @param[in] ordering Request a specific element ordering.
662  @return The dense matrix corresponding to the interpolation of the face
663  degrees of freedom of the master (coarse) face to the slave
664  (fine) face. */
665  const DenseMatrix* GetCoarseToFineInterpolation(
666  const Mesh::FaceInformation &face,
667  const DenseMatrix* ptMat);
668 };
669
670 /** @brief Operator that extracts face degrees of freedom for L2 nonconforming
671  spaces.
672
673  In order to support face restrictions on nonconforming meshes, this
674  operator interpolates master (coarse) face degrees of freedom onto the
675  slave (fine) face. This allows face integrators to treat nonconforming
676  faces just as regular conforming faces. */
678 {
679 protected:
681  mutable Vector x_interp;
682
683  /** @brief Constructs an NCL2FaceRestriction, this is a specialization of a
684  L2FaceRestriction for nonconforming meshes.
685
686  @param[in] fes The FiniteElementSpace on which this operates
687  @param[in] ordering Request a specific ordering
688  @param[in] type Request internal or boundary faces dofs
689  @param[in] m Request the face dofs for elem1, or both elem1 and
690  elem2
691  @param[in] build Request the NCL2FaceRestriction to compute the
692  scatter/gather indices. False should only be used
693  when inheriting from NCL2FaceRestriction.
694  */
696  const ElementDofOrdering ordering,
697  const FaceType type,
698  const L2FaceValues m,
699  bool build);
700 public:
701  /** @brief Constructs an NCL2FaceRestriction, this is a specialization of a
702  L2FaceRestriction for nonconforming meshes.
703
704  @param[in] fes The FiniteElementSpace on which this operates
705  @param[in] ordering Request a specific ordering
706  @param[in] type Request internal or boundary faces dofs
707  @param[in] m Request the face dofs for elem1, or both elem1 and
708  elem2
709  */
711  const ElementDofOrdering ordering,
712  const FaceType type,
714
715  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
716  face E-Vector.
717
718  @param[in] x The L-vector degrees of freedom.
719  @param[out] y The face E-Vector degrees of freedom with the given format:
720  if L2FacesValues::DoubleValued (face_dofs x vdim x 2 x nf),
721  if L2FacesValues::SingleValued (face_dofs x vdim x nf),
722  where nf is the number of interior or boundary faces
723  requested by @a type in the constructor.
724  The face_dofs are ordered according to the given
725  ElementDofOrdering. */
726  void Mult(const Vector &x, Vector &y) const override;
727
728  /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
729  L-Vector.
730
731  @param[in] x The face E-Vector degrees of freedom with the given format:
732  if L2FacesValues::DoubleValued (face_dofs x vdim x 2 x nf),
733  if L2FacesValues::SingleValued (face_dofs x vdim x nf),
734  where nf is the number of interior or boundary faces
735  requested by @a type in the constructor.
736  The face_dofs should be ordered according to the given
737  ElementDofOrdering
738  @param[in,out] y The L-vector degrees of freedom. */
739  void AddMultTranspose(const Vector &x, Vector &y) const override;
740
741  /** @brief Fill the I array of SparseMatrix corresponding to the sparsity
742  pattern given by this NCL2FaceRestriction.
743
744  @param[in,out] mat The sparse matrix for which we want to initialize the
745  row offsets.
746  @param[in] keep_nbr_block When set to true the SparseMatrix will
747  include the rows (in addition to the columns)
748  corresponding to face-neighbor dofs. The
749  default behavior is to disregard those rows.
750
751  @warning This method is not implemented yet. */
752  void FillI(SparseMatrix &mat,
753  const bool keep_nbr_block = false) const override;
754
755  /** @brief Fill the J and Data arrays of the SparseMatrix corresponding to
756  the sparsity pattern given by this NCL2FaceRestriction, and the values of
757  ea_data.
758
759  @param[in] fea_data The dense matrices representing the local operators
760  on each face. The format is:
761  face_dofs x face_dofs x 2 x nf.
762  On each face the first local matrix corresponds to
763  the contribution of elem1 on elem2, and the second to
764  the contribution of elem2 on elem1.
765  @param[in,out] mat The sparse matrix that is getting filled.
766  @param[in] keep_nbr_block When set to true the SparseMatrix will
767  include the rows (in addition to the columns)
768  corresponding to face-neighbor dofs. The
769  default behavior is to disregard those rows.
770
771  @warning This method is not implemented yet. */
772  void FillJAndData(const Vector &fea_data,
773  SparseMatrix &mat,
774  const bool keep_nbr_block = false) const override;
775
776  /** @brief This methods adds the DG face matrices to the element matrices.
777
778  @param[in] fea_data The dense matrices representing the local operators
779  on each face. The format is:
780  face_dofs x face_dofs x 2 x nf.
781  On each face the first and second local matrices
782  correspond to the contributions of elem1 and elem2 on
783  themselves respectively.
784  @param[in,out] ea_data The dense matrices representing the element local
785  contributions for each element to which will be
787  The format is: dofs x dofs x ne, where dofs is the
788  number of dofs per element and ne the number of
789  elements.
790
791  @warning This method is not implemented yet. */
793  Vector &ea_data) const override;
794
795 private:
796  /** @brief Compute the scatter indices: L-vector to E-vector, the offsets
797  for the gathering: E-vector to L-vector, and the interpolators from
798  coarse to fine face for master non-comforming faces.
799
800  @param[in] ordering Request a specific element ordering.
801  @param[in] type Request internal or boundary faces dofs.
802  */
803  void ComputeScatterIndicesAndOffsets(const ElementDofOrdering ordering,
804  const FaceType type);
805
806  /** @brief Compute the gather indices: E-vector to L-vector.
807
808  Note: Requires the gather offsets to be computed.
809
810  @param[in] ordering Request a specific element ordering.
811  @param[in] type Request internal or boundary faces dofs.
812  */
813  void ComputeGatherIndices(const ElementDofOrdering ordering,
814  const FaceType type);
815
816 public:
817  /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
818  face E-Vector. Should only be used with nonconforming faces and when:
819  L2FaceValues m == L2FaceValues::DoubleValued
820
821  @param[in] x The L-vector degrees of freedom.
822  @param[out] y The face E-Vector degrees of freedom with the given format:
823  (face_dofs x vdim x 2 x nf),
824  where nf is the number of interior or boundary faces
825  requested by @a type in the constructor.
826  The face_dofs are ordered according to the given
827  ElementDofOrdering. */
828  virtual void DoubleValuedNonconformingMult(const Vector& x, Vector& y) const;
829
830  /** @brief Apply a change of basis from fine element basis to coarse element
831  basis for the coarse face dofs. Should only be used when:
832  L2FaceValues m == L2FaceValues::SingleValued
833
834  @param[in] x The dofs vector that needs coarse dofs to be express in term
835  of the coarse basis, the result is stored in x_interp.
836  */
838
839  /** @brief Apply a change of basis from fine element basis to coarse element
840  basis for the coarse face dofs. Should only be used when:
841  L2FaceValues m == L2FaceValues::DoubleValued
842
843  @param[in] x The dofs vector that needs coarse dofs to be express in term
844  of the coarse basis, the result is stored in x_interp.
845  */
847 };
848
849 /** @brief Return the face map that extracts the degrees of freedom for the
850  requested local face of a quad or hex, returned in Lexicographic order.
851
852  @param[in] dim The dimension of the space
853  @param[in] face_id The local face identifier
854  @param[in] dof1d The 1D number of degrees of freedom for each dimension
855  @param[out] face_map The map that maps each face dof to an element dof
856 */
857 void GetFaceDofs(const int dim, const int face_id,
858  const int dof1d, Array<int> &face_map);
859
860 /** @brief Convert a dof face index from Native ordering to lexicographic
861  ordering for quads and hexes.
862
863  @param[in] dim The dimension of the element, 2 for quad, 3 for hex
864  @param[in] face_id The local face identifier
865  @param[in] size1d The 1D number of degrees of freedom for each dimension
866  @param[in] index The native index on the face
867  @return The lexicographic index on the face
868 */
869 int ToLexOrdering(const int dim, const int face_id, const int size1d,
870  const int index);
871
872 /** @brief Compute the dof face index of elem2 corresponding to the given dof
873  face index.
874
875  @param[in] dim The dimension of the element, 2 for quad, 3 for hex
876  @param[in] face_id1 The local face identifier of elem1
877  @param[in] face_id2 The local face identifier of elem2
878  @param[in] orientation The orientation of elem2 relative to elem1 on the
879  face
880  @param[in] size1d The 1D number of degrees of freedom for each dimension
881  @param[in] index The dof index on elem1
882  @return The dof index on elem2 facing the dof on elem1
883 */
884 int PermuteFaceL2(const int dim, const int face_id1,
885  const int face_id2, const int orientation,
886  const int size1d, const int index);
887
888 }
889
890 #endif // MFEM_RESTRICTION
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
H1FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering ordering, const FaceType type, bool build)
Construct an H1FaceRestriction.
void MultTranspose(const Vector &x, Vector &y) const override
Set the face degrees of freedom in the element degrees of freedom y to the values given in x...
void MultTransposeUnsigned(const Vector &x, Vector &y) const
Compute MultTranspose without applying signs based on DOF orientations.
void Mult(const Vector &x, Vector &y) const override=0
Extract the face degrees of freedom from x into y.
void PermuteAndSetSharedFaceDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the scattering indices of elem2 for the shared face described by the face...
Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
InterpConfig & operator=(const InterpConfig &rhs)=default
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
Array< int > gather_indices
void AddMultTranspose(const Vector &x, Vector &y) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
void SetBoundaryDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Set the scattering indices of elem2 for the boundary face described by the face.
std::pair< const DenseMatrix *, int > Key
void AddFaceMatricesToElementMatrices(const Vector &fea_data, Vector &ea_data) const override
This methods adds the DG face matrices to the element matrices.
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
virtual void DoubleValuedNonconformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with no...
void MultLeftInverse(const Vector &x, Vector &y) const
void SetFaceDofsGatherIndices1(const Mesh::FaceInformation &face, const int face_index)
Set the gathering indices of elem1 for the interior face described by the face.
FaceRestriction(int h, int w)
void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const override
Fill the I array of SparseMatrix corresponding to the sparsity pattern given by this NCL2FaceRestrict...
const L2FaceValues m
Array< int > gather_offsets
Array< int > gather_indices
void DoubleValuedConformingAddMultTranspose(const Vector &x, Vector &y) const
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector. Should only be used with con...
void PermuteAndSetFaceDofsGatherIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the gathering indices of elem2 for the interior face described by the face...
L2ElementRestriction(const FiniteElementSpace &)
Array< InterpConfig > interp_config
InterpConfig(int master_side, int nc_index)
void FillJAndData(const Vector &fea_data, SparseMatrix &mat, const bool keep_nbr_block=false) const override
Fill the J and Data arrays of the SparseMatrix corresponding to the sparsity pattern given by this NC...
FaceType
Definition: mesh.hpp:45
void RegisterFaceConformingInterpolation(const Mesh::FaceInformation &face, int face_index)
Register the face with face and index face_index as a conforming face for the interpolation of the de...
const FiniteElementSpace & fes
void PermuteAndSetFaceDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the scattering indices of elem2, and increment the offsets for the face described by ...
Data type sparse matrix.
Definition: sparsemat.hpp:46
int FillI(SparseMatrix &mat) const
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void SingleValuedConformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with co...
void FillI(SparseMatrix &mat) const
Operator that extracts Face degrees of freedom for L2 spaces.
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
Fills the E-vector y with boolean values 0.0 and 1.0 such that each each entry of the L-vector is uni...
void CheckFESpace(const ElementDofOrdering ordering)
Verify that H1FaceRestriction is build from an H1 FESpace.
const FiniteElementSpace & fes
void LinearizeInterpolatorMapIntoVector()
Transform the interpolation matrix map into a contiguous memory structure.
void AddMultTranspose(const Vector &x, Vector &y) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
void MultUnsigned(const Vector &x, Vector &y) const
Compute Mult without applying signs based on DOF orientations.
ElementRestriction(const FiniteElementSpace &, ElementDofOrdering)
Definition: restriction.cpp:27
void GetFaceDofs(const int dim, const int face_id, const int dof1d, Array< int > &face_map)
Return the face map that extracts the degrees of freedom for the requested local face of a quad or he...
void AddMultTranspose(const Vector &x, Vector &y) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
This class manages the storage and computation of the interpolations from master (coarse) face to sla...
This structure is used as a human readable output format that decipheres the information contained in...
Definition: mesh.hpp:1333
const ElementDofOrdering ordering
int ToLexOrdering(const int dim, const int face_id, const int size1d, const int index)
Convert a dof face index from Native ordering to lexicographic ordering for quads and hexes...
const FiniteElementSpace & fes
Definition: restriction.hpp:35
void SetFaceDofsGatherIndices(const Mesh::FaceInformation &face, const int face_index, const ElementDofOrdering ordering)
Set the gathering indices of elem1 for the interior face described by the face.
virtual void AddMultTranspose(const Vector &x, Vector &y) const =0
Add the face degrees of freedom x to the element degrees of freedom y.
Array< int > gather_offsets
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
int GetNumInterpolators() const
Return the total number of interpolators.
NCL2FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering ordering, const FaceType type, const L2FaceValues m, bool build)
Constructs an NCL2FaceRestriction, this is a specialization of a L2FaceRestriction for nonconforming ...
void SingleValuedConformingAddMultTranspose(const Vector &x, Vector &y) const
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector. Should only be used with con...
void CheckFESpace(const ElementDofOrdering ordering)
Verify that L2FaceRestriction is build from an L2 FESpace.
Array< int > scatter_indices1
L2FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering ordering, const FaceType type, const L2FaceValues m, bool build)
Constructs an L2FaceRestriction.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:66
int dim
Definition: ex24.cpp:53
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
const FiniteElementSpace & fes
void SetFaceDofsScatterIndices(const Mesh::FaceInformation &face, const int face_index, const ElementDofOrdering ordering)
Set the scattering indices of elem1, and increment the offsets for the face described by the face...
void SingleValuedNonconformingTransposeInterpolation(const Vector &x) const
Apply a change of basis from fine element basis to coarse element basis for the coarse face dofs...
Vector data type.
Definition: vector.hpp:60
void DoubleValuedNonconformingTransposeInterpolation(const Vector &x) const
Apply a change of basis from fine element basis to coarse element basis for the coarse face dofs...
virtual void DoubleValuedConformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with co...
virtual void FillJAndData(const Vector &fea_data, SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the J and Data arrays of the SparseMatrix corresponding to the sparsity pattern given by this L2...
int PermuteFaceL2(const int dim, const int face_id1, const int face_id2, const int orientation, const int size1d, const int index)
Compute the dof face index of elem2 corresponding to the given dof face index.
virtual void AddFaceMatricesToElementMatrices(const Vector &fea_data, Vector &ea_data) const
This methods adds the DG face matrices to the element matrices.
void SetFaceDofsScatterIndices1(const Mesh::FaceInformation &face, const int face_index)
Set the scattering indices of elem1, and increment the offsets for the face described by the face...
void RegisterFaceCoarseToFineInterpolation(const Mesh::FaceInformation &face, int face_index)
Register the face with face and index face_index as a conforming face for the interpolation of the de...
Base class for operators that extracts Face degrees of freedom.
Abstract operator.
Definition: operator.hpp:24
uint32_t is_non_conforming
const Array< InterpConfig > & GetFaceInterpConfig() const
Return an array containing the interpolation configuration for each face registered with RegisterFace...
Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:84
Array< int > scatter_indices
InterpolationManager interpolations
InterpConfig()=default
Array< int > scatter_indices2
virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the I array of SparseMatrix corresponding to the sparsity pattern given by this L2FaceRestrictio...
void FillSparseMatrix(const Vector &mat_ea, SparseMatrix &mat) const
Fill a Sparse Matrix with Element Matrices.
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
const Vector & GetInterpolators() const
Return an mfem::Vector containing the interpolators in the following format: face_dofs x face_dofs x ...
std::map< Key, std::pair< int, const DenseMatrix * >> Map
The temporary map used to store the different interpolators.