MFEM  v4.5.2
Finite element discretization library
linearform.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 // Implementation of class LinearForm
13 
14 #include "fem.hpp"
15 
16 namespace mfem
17 {
18 
20  : Vector(f->GetVSize())
21 {
22  ext = nullptr;
23  extern_lfs = 1;
24  fast_assembly = false;
25  fes = f;
26 
27  // Linear forms are stored on the device
28  UseDevice(true);
29 
30  // Copy the pointers to the integrators
32 
34 
36 
39 }
40 
42 {
43  DeltaLFIntegrator *maybe_delta =
44  dynamic_cast<DeltaLFIntegrator *>(lfi);
45  if (!maybe_delta || !maybe_delta->IsDelta())
46  {
47  domain_integs.Append(lfi);
48  }
49  else
50  {
51  domain_delta_integs.Append(maybe_delta);
52  }
53  domain_integs_marker.Append(NULL);
54 }
55 
57  Array<int> &elem_marker)
58 {
59  DeltaLFIntegrator *maybe_delta =
60  dynamic_cast<DeltaLFIntegrator *>(lfi);
61  if (!maybe_delta || !maybe_delta->IsDelta())
62  {
63  domain_integs.Append(lfi);
64  }
65  else
66  {
67  domain_delta_integs.Append(maybe_delta);
68  }
69  domain_integs_marker.Append(&elem_marker);
70 }
71 
73 {
74  boundary_integs.Append (lfi);
75  boundary_integs_marker.Append(NULL); // NULL -> all attributes are active
76 }
77 
79  Array<int> &bdr_attr_marker)
80 {
81  boundary_integs.Append (lfi);
82  boundary_integs_marker.Append(&bdr_attr_marker);
83 }
84 
86 {
87  boundary_face_integs.Append(lfi);
88  // NULL -> all attributes are active
89  boundary_face_integs_marker.Append(NULL);
90 }
91 
93  Array<int> &bdr_attr_marker)
94 {
95  boundary_face_integs.Append(lfi);
96  boundary_face_integs_marker.Append(&bdr_attr_marker);
97 }
98 
100 {
101  interior_face_integs.Append(lfi);
102 }
103 
105 {
106  // return false for NURBS meshes, so we don’t convert it to non-NURBS
107  // through Assemble, AssembleDevice, GetGeometricFactors and EnsureNodes
108  const Mesh &mesh = *fes->GetMesh();
109  if (mesh.NURBSext != nullptr) { return false; }
110 
111  // scan integrators to verify that all can use device assembly
112  auto IntegratorsSupportDevice = [](const Array<LinearFormIntegrator*> &integ)
113  {
114  for (int k = 0; k < integ.Size(); k++)
115  {
116  if (!integ[k]->SupportsDevice()) { return false; }
117  }
118  return true;
119  };
120 
121  if (!IntegratorsSupportDevice(domain_integs)) { return false; }
122  if (!IntegratorsSupportDevice(boundary_integs)) { return false; }
123  if (boundary_face_integs.Size() > 0 || interior_face_integs.Size() > 0 ||
124  domain_delta_integs.Size() > 0) { return false; }
125 
126  if (boundary_integs.Size() > 0)
127  {
128  // Make sure there are no boundary faces that are not boundary elements
130  {
131  return false;
132  }
133  // Make sure every boundary element corresponds to a boundary face
134  for (int be = 0; be < fes->GetNBE(); ++be)
135  {
136  const int f = mesh.GetBdrElementEdgeIndex(be);
137  const auto face_info = mesh.GetFaceInformation(f);
138  if (!face_info.IsBoundary())
139  {
140  return false;
141  }
142  }
143  }
144 
145  // no support for elements with varying polynomial orders
146  if (fes->IsVariableOrder()) { return false; }
147 
148  // no support for 1D and embedded meshes
149  const int mesh_dim = mesh.Dimension();
150  if (mesh_dim == 1 || mesh_dim != mesh.SpaceDimension()) { return false; }
151 
152  // tensor-product finite element space only
153  if (!UsesTensorBasis(*fes)) { return false; }
154 
155  return true;
156 }
157 
159 {
160  fast_assembly = use_fa;
161 
162  if (fast_assembly && SupportsDevice() && !ext)
163  {
164  ext = new LinearFormExtension(this);
165  }
166 }
167 
169 {
170  Array<int> vdofs;
171  ElementTransformation *eltrans;
172  DofTransformation *doftrans;
173  Vector elemvect;
174 
175  Vector::operator=(0.0);
176 
177  // The above operation is executed on device because of UseDevice().
178  // The first use of AddElementVector() below will move it back to host
179  // because both 'vdofs' and 'elemvect' are on host.
180 
181  if (fast_assembly && ext) { return ext->Assemble(); }
182 
183  if (domain_integs.Size())
184  {
185  for (int k = 0; k < domain_integs.Size(); k++)
186  {
187  if (domain_integs_marker[k] != NULL)
188  {
189  MFEM_VERIFY(domain_integs_marker[k]->Size() ==
190  (fes->GetMesh()->attributes.Size() ?
191  fes->GetMesh()->attributes.Max() : 0),
192  "invalid element marker for domain linear form "
193  "integrator #" << k << ", counting from zero");
194  }
195  }
196 
197  for (int i = 0; i < fes -> GetNE(); i++)
198  {
199  int elem_attr = fes->GetMesh()->GetAttribute(i);
200  for (int k = 0; k < domain_integs.Size(); k++)
201  {
202  const Array<int> * const markers = domain_integs_marker[k];
203  if ( markers == NULL || (*markers)[elem_attr-1] == 1 )
204  {
205  doftrans = fes -> GetElementVDofs (i, vdofs);
206  eltrans = fes -> GetElementTransformation (i);
207  domain_integs[k]->AssembleRHSElementVect(*fes->GetFE(i),
208  *eltrans, elemvect);
209  if (doftrans)
210  {
211  doftrans->TransformDual(elemvect);
212  }
213  AddElementVector (vdofs, elemvect);
214  }
215  }
216  }
217  }
218  AssembleDelta();
219 
220  if (boundary_integs.Size())
221  {
222  Mesh *mesh = fes->GetMesh();
223 
224  // Which boundary attributes need to be processed?
225  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
226  mesh->bdr_attributes.Max() : 0);
227  bdr_attr_marker = 0;
228  for (int k = 0; k < boundary_integs.Size(); k++)
229  {
230  if (boundary_integs_marker[k] == NULL)
231  {
232  bdr_attr_marker = 1;
233  break;
234  }
235  Array<int> &bdr_marker = *boundary_integs_marker[k];
236  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
237  "invalid boundary marker for boundary integrator #"
238  << k << ", counting from zero");
239  for (int i = 0; i < bdr_attr_marker.Size(); i++)
240  {
241  bdr_attr_marker[i] |= bdr_marker[i];
242  }
243  }
244 
245  for (int i = 0; i < fes -> GetNBE(); i++)
246  {
247  const int bdr_attr = mesh->GetBdrAttribute(i);
248  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
249  doftrans = fes -> GetBdrElementVDofs (i, vdofs);
250  eltrans = fes -> GetBdrElementTransformation (i);
251  for (int k=0; k < boundary_integs.Size(); k++)
252  {
253  if (boundary_integs_marker[k] &&
254  (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
255 
256  boundary_integs[k]->AssembleRHSElementVect(*fes->GetBE(i),
257  *eltrans, elemvect);
258 
259  if (doftrans)
260  {
261  doftrans->TransformDual(elemvect);
262  }
263  AddElementVector (vdofs, elemvect);
264  }
265  }
266  }
267  if (boundary_face_integs.Size())
268  {
270  Mesh *mesh = fes->GetMesh();
271 
272  // Which boundary attributes need to be processed?
273  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
274  mesh->bdr_attributes.Max() : 0);
275  bdr_attr_marker = 0;
276  for (int k = 0; k < boundary_face_integs.Size(); k++)
277  {
278  if (boundary_face_integs_marker[k] == NULL)
279  {
280  bdr_attr_marker = 1;
281  break;
282  }
283  Array<int> &bdr_marker = *boundary_face_integs_marker[k];
284  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
285  "invalid boundary marker for boundary face integrator #"
286  << k << ", counting from zero");
287  for (int i = 0; i < bdr_attr_marker.Size(); i++)
288  {
289  bdr_attr_marker[i] |= bdr_marker[i];
290  }
291  }
292 
293  for (int i = 0; i < mesh->GetNBE(); i++)
294  {
295  const int bdr_attr = mesh->GetBdrAttribute(i);
296  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
297 
298  tr = mesh->GetBdrFaceTransformations(i);
299  if (tr != NULL)
300  {
301  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
302  for (int k = 0; k < boundary_face_integs.Size(); k++)
303  {
305  (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
306  { continue; }
307 
309  AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
310  *tr, elemvect);
311  AddElementVector (vdofs, elemvect);
312  }
313  }
314  }
315  }
316 
317  if (interior_face_integs.Size())
318  {
319  Mesh *mesh = fes->GetMesh();
320 
321  for (int k = 0; k < interior_face_integs.Size(); k++)
322  {
323  for (int i = 0; i < mesh->GetNumFaces(); i++)
324  {
325  FaceElementTransformations *tr = NULL;
326  tr = mesh->GetInteriorFaceTransformations (i);
327  if (tr != NULL)
328  {
329  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
330  Array<int> vdofs2;
331  fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
332  vdofs.Append(vdofs2);
334  AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
335  *fes->GetFE(tr->Elem2No),
336  *tr, elemvect);
337  AddElementVector (vdofs, elemvect);
338  }
339  }
340  }
341  }
342 }
343 
345 {
347  if (ext) { ext->Update(); }
348 }
349 
351 {
352  MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
353  fes = f;
354  v.UseDevice(true);
355  this->Vector::MakeRef(v, v_offset, fes->GetVSize());
357  if (ext) { ext->Update(); }
358 }
359 
361 {
362  Update(f, v, v_offset);
363 }
364 
366 {
367  if (domain_delta_integs.Size() == 0) { return; }
368 
369  if (!HaveDeltaLocations())
370  {
371  int sdim = fes->GetMesh()->SpaceDimension();
372  Vector center;
373  DenseMatrix centers(sdim, domain_delta_integs.Size());
374  for (int i = 0; i < centers.Width(); i++)
375  {
376  centers.GetColumnReference(i, center);
377  domain_delta_integs[i]->GetDeltaCenter(center);
378  MFEM_VERIFY(center.Size() == sdim,
379  "Point dim " << center.Size() <<
380  " does not match space dim " << sdim);
381  }
384  }
385 
386  Array<int> vdofs;
387  Vector elemvect;
388  for (int i = 0; i < domain_delta_integs.Size(); i++)
389  {
390  int elem_id = domain_delta_integs_elem_id[i];
391  // The delta center may be outside of this sub-domain, or
392  // (Par)Mesh::FindPoints() failed to find this point:
393  if (elem_id < 0) { continue; }
394 
397  Trans.SetIntPoint(&ip);
398 
399  fes->GetElementVDofs(elem_id, vdofs);
400  domain_delta_integs[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id),
401  Trans, elemvect);
402  AddElementVector(vdofs, elemvect);
403  }
404 }
405 
407 {
408  Vector::operator=(value);
409  return *this;
410 }
411 
413 {
414  MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
416  return *this;
417 }
418 
420 {
421  if (!extern_lfs)
422  {
423  int k;
424  for (k=0; k < domain_delta_integs.Size(); k++)
425  { delete domain_delta_integs[k]; }
426  for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
427  for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
428  for (k=0; k < boundary_face_integs.Size(); k++)
429  { delete boundary_face_integs[k]; }
430  for (k=0; k < interior_face_integs.Size(); k++)
431  { delete interior_face_integs[k]; }
432  }
433 
434  delete ext;
435 }
436 
437 }
Class extending the LinearForm class to support assembly on devices.
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:463
LinearFormExtension * ext
Extension for supporting different assembly levels.
Definition: linearform.hpp:33
int GetBdrElementEdgeIndex(int i) const
Definition: mesh.cpp:6235
int Dimension() const
Definition: mesh.hpp:1047
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 MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
Make the LinearForm reference external data on a new FiniteElementSpace.
Definition: linearform.cpp:360
Array< Array< int > * > boundary_face_integs_marker
Entries not owned.
Definition: linearform.hpp:63
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
int extern_lfs
Indicates the LinearFormIntegrators stored in domain_integs, domain_delta_integs, boundary_integs...
Definition: linearform.hpp:42
LinearForm & operator=(const LinearForm &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition: linearform.hpp:121
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
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
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Array< IntegrationPoint > domain_delta_integs_ip
The reference coordinates where the centers of the delta functions lie.
Definition: linearform.hpp:72
void ResetDeltaLocations()
Force (re)computation of delta locations.
Definition: linearform.hpp:79
void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Interior Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:99
Array< int > domain_delta_integs_elem_id
The element ids where the centers of the delta functions lie.
Definition: linearform.hpp:69
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:957
Array< LinearFormIntegrator * > boundary_face_integs
Set of Boundary Face Integrators to be applied.
Definition: linearform.hpp:62
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition: mesh.hpp:1327
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:939
Abstract base class LinearFormIntegrator.
Definition: lininteg.hpp:24
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition: fespace.hpp:631
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:124
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1544
Array< LinearFormIntegrator * > interior_face_integs
Set of Internal Face Integrators to be applied.
Definition: linearform.hpp:66
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
Array< LinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
Definition: linearform.hpp:45
void Update()
Update the linear form extension.
FiniteElementSpace * fes
FE space on which the LinearForm lives. Not owned.
Definition: linearform.hpp:30
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:85
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1095
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:72
virtual bool SupportsDevice()
Return true if assembly on device is supported, false otherwise.
Definition: linearform.cpp:104
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:640
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
LinearForm()
Create an empty LinearForm without an associated FiniteElementSpace.
Definition: linearform.hpp:104
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
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 GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
Abstract class for integrators that support delta coefficients.
Definition: lininteg.hpp:62
void AssembleDelta()
Assembles delta functions of the linear form.
Definition: linearform.cpp:365
Array< DeltaLFIntegrator * > domain_delta_integs
Separate array for integrators with delta function coefficients.
Definition: linearform.hpp:54
int SpaceDimension() const
Definition: mesh.hpp:1048
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.cpp:344
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
Class for integration point with weight.
Definition: intrules.hpp:25
~LinearForm()
Destroys linear form.
Definition: linearform.cpp:419
bool fast_assembly
Should we use the device-compatible fast assembly algorithm (false by default)
Definition: linearform.hpp:37
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:623
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
bool HaveDeltaLocations()
If true, the delta locations are not (re)computed during assembly.
Definition: linearform.hpp:75
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
FaceInformation GetFaceInformation(int f) const
Definition: mesh.cpp:1131
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:120
Array< LinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
Definition: linearform.hpp:57
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
Definition: mesh.cpp:11962
Vector data type.
Definition: vector.hpp:60
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
bool IsDelta() const
Returns true if the derived class instance uses a delta coefficient.
Definition: lininteg.hpp:85
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
Array< Array< int > * > domain_integs_marker
Definition: linearform.hpp:51
virtual void TransformDual(double *v) const =0
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3102
void UseFastAssembly(bool use_fa)
Which assembly algorithm to use: the new device-compatible fast assembly (true), or the legacy CPU-on...
Definition: linearform.cpp:158
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:273
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647