MFEM  v4.5.2
Finite element discretization library
linearform_ext.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "linearform.hpp"
13 #include "../general/forall.hpp"
14 
15 namespace mfem
16 {
17 
19 
21 {
22  const FiniteElementSpace &fes = *lf->FESpace();
23  MFEM_VERIFY(lf->SupportsDevice(), "Not supported.");
24  MFEM_VERIFY(lf->Size() == fes.GetVSize(), "LinearForm size does not "
25  "match the number of vector dofs!");
26 
27  const Array<Array<int>*> &domain_integs_marker = *lf->GetDLFI_Marker();
28  const int mesh_attributes_max = fes.GetMesh()->attributes.Size() ?
29  fes.GetMesh()->attributes.Max() : 0;
30  const Array<LinearFormIntegrator*> &domain_integs = *lf->GetDLFI();
31 
32  for (int k = 0; k < domain_integs.Size(); ++k)
33  {
34  // Get the markers for this integrator
35  const Array<int> *domain_integs_marker_k = domain_integs_marker[k];
36 
37  // check if there are markers for this integrator
38  const bool has_markers_k = domain_integs_marker_k != nullptr;
39 
40  if (has_markers_k)
41  {
42  // Element attribute marker should be of length mesh->attributes
43  MFEM_VERIFY(mesh_attributes_max == domain_integs_marker_k->Size(),
44  "invalid element marker for domain linear form "
45  "integrator #" << k << ", counting from zero");
46  }
47 
48  // if there are no markers, just use the whole linear form (1)
49  if (!has_markers_k) { markers.HostReadWrite(); markers = 1; }
50  else
51  {
52  // scan the attributes to set the markers to 0 or 1
53  const int NE = fes.GetNE();
54  const auto attr = attributes.Read();
55  const auto dimk = domain_integs_marker_k->Read();
56  auto markers_w = markers.Write();
57  MFEM_FORALL(e, NE, markers_w[e] = dimk[attr[e]-1] == 1;);
58  }
59 
60  // Assemble the linear form
61  b = 0.0;
62  domain_integs[k]->AssembleDevice(fes, markers, b);
63  if (k == 0) { elem_restrict_lex->MultTranspose(b, *lf); }
64  else { elem_restrict_lex->AddMultTranspose(b, *lf); }
65  }
66 
67  const Array<Array<int>*> &boundary_integs_marker = lf->boundary_integs_marker;
68  const int bdr_attributes_max = fes.GetMesh()->bdr_attributes.Size() ?
69  fes.GetMesh()->bdr_attributes.Max() : 0;
70  const Array<LinearFormIntegrator*> &boundary_integs = lf->boundary_integs;
71 
72  for (int k = 0; k < boundary_integs.Size(); ++k)
73  {
74  // Get the markers for this integrator
75  const Array<int> *boundary_integs_marker_k = boundary_integs_marker[k];
76 
77  // check if there are markers for this integrator
78  const bool has_markers_k = boundary_integs_marker_k != nullptr;
79 
80  if (has_markers_k)
81  {
82  // Element attribute marker should be of length mesh->attributes
83  MFEM_VERIFY(bdr_attributes_max == boundary_integs_marker_k->Size(),
84  "invalid boundary marker for boundary linear form "
85  "integrator #" << k << ", counting from zero");
86  }
87 
88  // if there are no markers, just use the whole linear form (1)
89  if (!has_markers_k) { bdr_markers.HostReadWrite(); bdr_markers = 1; }
90  else
91  {
92  // scan the attributes to set the markers to 0 or 1
93  const int NBE = bdr_attributes.Size();
94  const auto attr = bdr_attributes.Read();
95  const auto attr_markers = boundary_integs_marker_k->Read();
96  auto markers_w = bdr_markers.Write();
97  MFEM_FORALL(e, NBE, markers_w[e] = attr_markers[attr[e]-1] == 1;);
98  }
99 
100  // Assemble the linear form
101  bdr_b = 0.0;
102  boundary_integs[k]->AssembleDevice(fes, bdr_markers, bdr_b);
103  bdr_restrict_lex->AddMultTranspose(bdr_b, *lf);
104  }
105 }
106 
108 {
109  const FiniteElementSpace &fes = *lf->FESpace();
110  const Mesh &mesh = *fes.GetMesh();
112 
113  MFEM_VERIFY(lf->Size() == fes.GetVSize(), "");
114 
115  if (lf->domain_integs.Size() > 0)
116  {
117  const int NE = fes.GetNE();
118  markers.SetSize(NE);
119  //markers.UseDevice(true);
120 
121  // Gather the attributes on the host from all the elements
122  attributes.SetSize(NE);
123  for (int i = 0; i < NE; ++i) { attributes[i] = mesh.GetAttribute(i); }
124 
125  elem_restrict_lex = fes.GetElementRestriction(ordering);
126  MFEM_VERIFY(elem_restrict_lex, "Element restriction not available");
127  b.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
128  b.UseDevice(true);
129  }
130 
131  if (lf->boundary_integs.Size() > 0)
132  {
133  const int nf_bdr = fes.GetNFbyType(FaceType::Boundary);
134  bdr_markers.SetSize(nf_bdr);
135  // bdr_markers.UseDevice(true);
136 
137  // The face restriction will give us "face E-vectors" on the boundary that
138  // are numbered in the order of the faces of mesh. This numbering will be
139  // different than the numbering of the boundary elements. We compute
140  // mappings so that the array `bdr_attributes[i]` gives the boundary
141  // attribute of the `i`th boundary face in the mesh face order.
142  std::unordered_map<int,int> f_to_be;
143  for (int i = 0; i < mesh.GetNBE(); ++i)
144  {
145  const int f = mesh.GetBdrElementEdgeIndex(i);
146  f_to_be[f] = i;
147  }
148  MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");
149  bdr_attributes.SetSize(nf_bdr);
150  int f_ind = 0;
151  for (int f = 0; f < mesh.GetNumFaces(); ++f)
152  {
153  if (f_to_be.find(f) != f_to_be.end())
154  {
155  const int be = f_to_be[f];
156  bdr_attributes[f_ind] = mesh.GetBdrAttribute(be);
157  ++f_ind;
158  }
159  }
160 
161  bdr_restrict_lex =
162  dynamic_cast<const FaceRestriction*>(
165  MFEM_VERIFY(bdr_restrict_lex, "Face restriction not available");
166  bdr_b.SetSize(bdr_restrict_lex->Height(), Device::GetMemoryType());
167  bdr_b.UseDevice(true);
168  }
169 }
170 
171 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
LinearFormExtension(LinearForm *lf)
Create a LinearForm extension of lf.
virtual 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 ...
Definition: operator.hpp:93
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:6235
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5389
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override=0
Add the E-vector degrees of freedom x to the L-vector degrees of freedom y.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:939
FiniteElementSpace * FESpace()
Read+write access to the associated FiniteElementSpace.
Definition: linearform.hpp:129
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1259
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition: fespace.hpp:631
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1544
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override=0
Add the face degrees of freedom x to the element degrees of freedom y.
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition: device.hpp:277
Array< LinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
Definition: linearform.hpp:45
void Update()
Update the linear form extension.
virtual bool SupportsDevice()
Return true if assembly on device is supported, false otherwise.
Definition: linearform.cpp:104
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1550
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
Lexicographic ordering for tensor-product FiniteElements.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
T * Write(bool on_dev=true)
Shortcut for mfem::Write(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:315
Array< LinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
Definition: linearform.hpp:57
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
Definition: linearform.hpp:59
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType, L2FaceValues mul=L2FaceValues::DoubleValued) const
Return an Operator that converts L-vectors to E-vectors on each face.
Definition: fespace.cpp:1292
Array< Array< int > * > * GetDLFI_Marker()
Access all domain markers added with AddDomainIntegrator(). If no marker was specified when the integ...
Definition: linearform.hpp:173
Base class for operators that extracts Face degrees of freedom.
Array< LinearFormIntegrator * > * GetDLFI()
Access all integrators added with AddDomainIntegrator() which are not DeltaLFIntegrators or they are ...
Definition: linearform.hpp:168
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:273