MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
linearform.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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
16namespace 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 and the corresponding marker arrays
33
35
38
41
43}
44
46{
47 DeltaLFIntegrator *maybe_delta =
48 dynamic_cast<DeltaLFIntegrator *>(lfi);
49 if (!maybe_delta || !maybe_delta->IsDelta())
50 {
51 domain_integs.Append(lfi);
52 }
53 else
54 {
55 domain_delta_integs.Append(maybe_delta);
56 }
57 domain_integs_marker.Append(NULL);
58}
59
61 Array<int> &elem_marker)
62{
63 DeltaLFIntegrator *maybe_delta =
64 dynamic_cast<DeltaLFIntegrator *>(lfi);
65 if (!maybe_delta || !maybe_delta->IsDelta())
66 {
67 domain_integs.Append(lfi);
68 }
69 else
70 {
71 domain_delta_integs.Append(maybe_delta);
72 }
73 domain_integs_marker.Append(&elem_marker);
74}
75
77{
78 boundary_integs.Append (lfi);
79 boundary_integs_marker.Append(NULL); // NULL -> all attributes are active
80}
81
83 Array<int> &bdr_attr_marker)
84{
85 boundary_integs.Append (lfi);
86 boundary_integs_marker.Append(&bdr_attr_marker);
87}
88
90{
91 boundary_face_integs.Append(lfi);
92 // NULL -> all attributes are active
93 boundary_face_integs_marker.Append(NULL);
94}
95
97 Array<int> &bdr_attr_marker)
98{
99 boundary_face_integs.Append(lfi);
100 boundary_face_integs_marker.Append(&bdr_attr_marker);
101}
102
107
109{
110 // return false for NURBS meshes, so we don’t convert it to non-NURBS
111 // through Assemble, AssembleDevice, GetGeometricFactors and EnsureNodes
112 const Mesh &mesh = *fes->GetMesh();
113 if (mesh.NURBSext != nullptr) { return false; }
114
115 // scan integrators to verify that all can use device assembly
116 auto IntegratorsSupportDevice = [](const Array<LinearFormIntegrator*> &integ)
117 {
118 for (int k = 0; k < integ.Size(); k++)
119 {
120 if (!integ[k]->SupportsDevice()) { return false; }
121 }
122 return true;
123 };
124
125 if (!IntegratorsSupportDevice(domain_integs)) { return false; }
126 if (!IntegratorsSupportDevice(boundary_integs)) { return false; }
127 if (boundary_face_integs.Size() > 0 || interior_face_integs.Size() > 0 ||
128 domain_delta_integs.Size() > 0) { return false; }
129
130 if (boundary_integs.Size() > 0)
131 {
132 // Make sure there are no boundary faces that are not boundary elements
134 {
135 return false;
136 }
137 // Make sure every boundary element corresponds to a boundary face
138 for (int be = 0; be < fes->GetNBE(); ++be)
139 {
140 const int f = mesh.GetBdrElementFaceIndex(be);
141 const auto face_info = mesh.GetFaceInformation(f);
142 if (!face_info.IsBoundary())
143 {
144 return false;
145 }
146 }
147 }
148
149 // no support for elements with varying polynomial orders
150 if (fes->IsVariableOrder()) { return false; }
151
152 // no support for 1D and embedded meshes
153 const int mesh_dim = mesh.Dimension();
154 if (mesh_dim == 1 || mesh_dim != mesh.SpaceDimension()) { return false; }
155
156 // tensor-product finite element space only
157 if (!UsesTensorBasis(*fes)) { return false; }
158
159 return true;
160}
161
163{
164 fast_assembly = use_fa;
165
166 if (fast_assembly && SupportsDevice() && !ext)
167 {
168 ext = new LinearFormExtension(this);
169 }
170}
171
173{
174 Array<int> vdofs;
175 ElementTransformation *eltrans;
176 DofTransformation *doftrans;
177 Vector elemvect;
178
180
181 // The above operation is executed on device because of UseDevice().
182 // The first use of AddElementVector() below will move it back to host
183 // because both 'vdofs' and 'elemvect' are on host.
184
185 if (fast_assembly && ext) { return ext->Assemble(); }
186
187 if (domain_integs.Size())
188 {
189 for (int k = 0; k < domain_integs.Size(); k++)
190 {
191 if (domain_integs_marker[k] != NULL)
192 {
193 MFEM_VERIFY(domain_integs_marker[k]->Size() ==
194 (fes->GetMesh()->attributes.Size() ?
195 fes->GetMesh()->attributes.Max() : 0),
196 "invalid element marker for domain linear form "
197 "integrator #" << k << ", counting from zero");
198 }
199 }
200
201 for (int i = 0; i < fes -> GetNE(); i++)
202 {
203 int elem_attr = fes->GetMesh()->GetAttribute(i);
204 for (int k = 0; k < domain_integs.Size(); k++)
205 {
206 const Array<int> * const markers = domain_integs_marker[k];
207 if (markers) { markers->HostRead(); }
208 if ( markers == NULL || (*markers)[elem_attr-1] == 1 )
209 {
210 doftrans = fes -> GetElementVDofs (i, vdofs);
211 eltrans = fes -> GetElementTransformation (i);
212 domain_integs[k]->AssembleRHSElementVect(*fes->GetFE(i),
213 *eltrans, elemvect);
214 if (doftrans)
215 {
216 doftrans->TransformDual(elemvect);
217 }
218 AddElementVector (vdofs, elemvect);
219 }
220 }
221 }
222 }
224
225 if (boundary_integs.Size())
226 {
227 Mesh *mesh = fes->GetMesh();
228
229 // Which boundary attributes need to be processed?
230 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
231 mesh->bdr_attributes.Max() : 0);
232 bdr_attr_marker = 0;
233 for (int k = 0; k < boundary_integs.Size(); k++)
234 {
235 if (boundary_integs_marker[k] == NULL)
236 {
237 bdr_attr_marker = 1;
238 break;
239 }
240 Array<int> &bdr_marker = *boundary_integs_marker[k];
241 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
242 "invalid boundary marker for boundary integrator #"
243 << k << ", counting from zero");
244 for (int i = 0; i < bdr_attr_marker.Size(); i++)
245 {
246 bdr_attr_marker[i] |= bdr_marker[i];
247 }
248 }
249
250 for (int i = 0; i < fes -> GetNBE(); i++)
251 {
252 const int bdr_attr = mesh->GetBdrAttribute(i);
253 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
254 doftrans = fes -> GetBdrElementVDofs (i, vdofs);
255 eltrans = fes -> GetBdrElementTransformation (i);
256 for (int k=0; k < boundary_integs.Size(); k++)
257 {
258 if (boundary_integs_marker[k] &&
259 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
260
261 boundary_integs[k]->AssembleRHSElementVect(*fes->GetBE(i),
262 *eltrans, elemvect);
263
264 if (doftrans)
265 {
266 doftrans->TransformDual(elemvect);
267 }
268 AddElementVector (vdofs, elemvect);
269 }
270 }
271 }
272 if (boundary_face_integs.Size())
273 {
275 Mesh *mesh = fes->GetMesh();
276
277 // Which boundary attributes need to be processed?
278 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
279 mesh->bdr_attributes.Max() : 0);
280 bdr_attr_marker = 0;
281 for (int k = 0; k < boundary_face_integs.Size(); k++)
282 {
283 if (boundary_face_integs_marker[k] == NULL)
284 {
285 bdr_attr_marker = 1;
286 break;
287 }
288 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
289 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
290 "invalid boundary marker for boundary face integrator #"
291 << k << ", counting from zero");
292 for (int i = 0; i < bdr_attr_marker.Size(); i++)
293 {
294 bdr_attr_marker[i] |= bdr_marker[i];
295 }
296 }
297
298 for (int i = 0; i < mesh->GetNBE(); i++)
299 {
300 const int bdr_attr = mesh->GetBdrAttribute(i);
301 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
302
303 tr = mesh->GetBdrFaceTransformations(i);
304 if (tr != NULL)
305 {
306 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
307 for (int k = 0; k < boundary_face_integs.Size(); k++)
308 {
310 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
311 { continue; }
312
314 AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
315 *tr, elemvect);
316 AddElementVector (vdofs, elemvect);
317 }
318 }
319 }
320 }
321
322 if (interior_face_integs.Size())
323 {
324 Mesh *mesh = fes->GetMesh();
325
326 for (int k = 0; k < interior_face_integs.Size(); k++)
327 {
328 for (int i = 0; i < mesh->GetNumFaces(); i++)
329 {
331 tr = mesh->GetInteriorFaceTransformations (i);
332 if (tr != NULL)
333 {
334 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
335 Array<int> vdofs2;
336 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
337 vdofs.Append(vdofs2);
339 AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
340 *fes->GetFE(tr->Elem2No),
341 *tr, elemvect);
342 AddElementVector (vdofs, elemvect);
343 }
344 }
345 }
346 }
347}
348
350{
352 if (ext) { ext->Update(); }
353}
354
356{
357 MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
358 fes = f;
359 v.UseDevice(true);
360 this->Vector::MakeRef(v, v_offset, fes->GetVSize());
362 if (ext) { ext->Update(); }
363}
364
366{
367 Update(f, v, v_offset);
368}
369
371{
372 if (domain_delta_integs.Size() == 0) { return; }
373
374 if (!HaveDeltaLocations())
375 {
376 int sdim = fes->GetMesh()->SpaceDimension();
377 Vector center;
378 DenseMatrix centers(sdim, domain_delta_integs.Size());
379 for (int i = 0; i < centers.Width(); i++)
380 {
381 centers.GetColumnReference(i, center);
382 domain_delta_integs[i]->GetDeltaCenter(center);
383 MFEM_VERIFY(center.Size() == sdim,
384 "Point dim " << center.Size() <<
385 " does not match space dim " << sdim);
386 }
389 }
390
391 Array<int> vdofs;
392 Vector elemvect;
393 for (int i = 0; i < domain_delta_integs.Size(); i++)
394 {
395 int elem_id = domain_delta_integs_elem_id[i];
396 // The delta center may be outside of this sub-domain, or
397 // (Par)Mesh::FindPoints() failed to find this point:
398 if (elem_id < 0) { continue; }
399
402 Trans.SetIntPoint(&ip);
403
404 fes->GetElementVDofs(elem_id, vdofs);
405 domain_delta_integs[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id),
406 Trans, elemvect);
407 AddElementVector(vdofs, elemvect);
408 }
409}
410
412{
413 Vector::operator=(value);
414 return *this;
415}
416
418{
419 MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
421 return *this;
422}
423
425{
426 if (!extern_lfs)
427 {
428 int k;
429 for (k=0; k < domain_delta_integs.Size(); k++)
430 { delete domain_delta_integs[k]; }
431 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
432 for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
433 for (k=0; k < boundary_face_integs.Size(); k++)
434 { delete boundary_face_integs[k]; }
435 for (k=0; k < interior_face_integs.Size(); k++)
436 { delete interior_face_integs[k]; }
437 }
438
439 delete ext;
440}
441
442}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition array.hpp:341
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
Abstract class for integrators that support delta coefficients.
Definition lininteg.hpp:60
bool IsDelta() const
Returns true if the derived class instance uses a delta coefficient.
Definition lininteg.hpp:82
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void GetColumnReference(int c, Vector &col)
Definition densemat.hpp:319
void TransformDual(real_t *v) const
Definition doftrans.cpp:81
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:750
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3881
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:924
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition fespace.hpp:900
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:332
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3835
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition fespace.hpp:908
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:848
Class for integration point with weight.
Definition intrules.hpp:35
Class extending the LinearForm class to support assembly on devices.
void Update()
Update the linear form extension.
Abstract base class LinearFormIntegrator.
Definition lininteg.hpp:26
Vector with associated FE space and LinearFormIntegrators.
void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Interior Face Integrator. Assumes ownership of lfi.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Array< Array< int > * > boundary_face_integs_marker
Entries not owned.
void AssembleDelta()
Assembles delta functions of the linear form.
virtual void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
Make the LinearForm reference external data on a new FiniteElementSpace.
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
LinearForm()
Create an empty LinearForm without an associated FiniteElementSpace.
Array< LinearFormIntegrator * > boundary_face_integs
Set of Boundary Face Integrators to be applied.
void UseFastAssembly(bool use_fa)
Which assembly algorithm to use: the new device-compatible fast assembly (true), or the legacy CPU-on...
Array< int > domain_delta_integs_elem_id
The element ids where the centers of the delta functions lie.
virtual bool SupportsDevice() const
Return true if assembly on device is supported, false otherwise.
LinearForm & operator=(const LinearForm &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
bool HaveDeltaLocations()
If true, the delta locations are not (re)computed during assembly.
FiniteElementSpace * fes
FE space on which the LinearForm lives. Not owned.
Array< Array< int > * > domain_integs_marker
Array< LinearFormIntegrator * > interior_face_integs
Set of Internal Face Integrators to be applied.
LinearFormExtension * ext
Extension for supporting different assembly levels.
void ResetDeltaLocations()
Force (re)computation of delta locations.
void Update()
Update the object according to the associated FE space fes.
Array< DeltaLFIntegrator * > domain_delta_integs
Separate array for integrators with delta function coefficients.
Array< LinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
~LinearForm()
Destroys linear form.
Array< IntegrationPoint > domain_delta_integs_ip
The reference coordinates where the centers of the delta functions lie.
Array< LinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
int extern_lfs
Indicates the LinearFormIntegrators stored in domain_integs, domain_delta_integs, boundary_integs,...
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
bool fast_assembly
Should we use the device-compatible fast assembly algorithm (false by default)
Mesh data type.
Definition mesh.hpp:64
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:290
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:298
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6531
int GetAttribute(int i) const
Return the attribute of element i.
Definition mesh.hpp:1389
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1395
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
Definition mesh.hpp:1588
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
FaceInformation GetFaceInformation(int f) const
Definition mesh.cpp:1193
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1123
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1285
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:13406
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition mesh.cpp:1103
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:288
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Vector data type.
Definition vector.hpp:82
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:745
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition vector.hpp:147
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:144
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
Vector & operator=(const real_t *v)
Copy Size() entries from v.
Definition vector.cpp:196
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:622
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1555
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30