MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
linearform.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
104{
105 interior_face_integs.Append(lfi);
106}
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 == NULL || (*markers)[elem_attr-1] == 1 )
208 {
209 doftrans = fes -> GetElementVDofs (i, vdofs);
210 eltrans = fes -> GetElementTransformation (i);
211 domain_integs[k]->AssembleRHSElementVect(*fes->GetFE(i),
212 *eltrans, elemvect);
213 if (doftrans)
214 {
215 doftrans->TransformDual(elemvect);
216 }
217 AddElementVector (vdofs, elemvect);
218 }
219 }
220 }
221 }
223
224 if (boundary_integs.Size())
225 {
226 Mesh *mesh = fes->GetMesh();
227
228 // Which boundary attributes need to be processed?
229 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
230 mesh->bdr_attributes.Max() : 0);
231 bdr_attr_marker = 0;
232 for (int k = 0; k < boundary_integs.Size(); k++)
233 {
234 if (boundary_integs_marker[k] == NULL)
235 {
236 bdr_attr_marker = 1;
237 break;
238 }
239 Array<int> &bdr_marker = *boundary_integs_marker[k];
240 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
241 "invalid boundary marker for boundary integrator #"
242 << k << ", counting from zero");
243 for (int i = 0; i < bdr_attr_marker.Size(); i++)
244 {
245 bdr_attr_marker[i] |= bdr_marker[i];
246 }
247 }
248
249 for (int i = 0; i < fes -> GetNBE(); i++)
250 {
251 const int bdr_attr = mesh->GetBdrAttribute(i);
252 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
253 doftrans = fes -> GetBdrElementVDofs (i, vdofs);
254 eltrans = fes -> GetBdrElementTransformation (i);
255 for (int k=0; k < boundary_integs.Size(); k++)
256 {
257 if (boundary_integs_marker[k] &&
258 (*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
259
260 boundary_integs[k]->AssembleRHSElementVect(*fes->GetBE(i),
261 *eltrans, elemvect);
262
263 if (doftrans)
264 {
265 doftrans->TransformDual(elemvect);
266 }
267 AddElementVector (vdofs, elemvect);
268 }
269 }
270 }
271 if (boundary_face_integs.Size())
272 {
274 Mesh *mesh = fes->GetMesh();
275
276 // Which boundary attributes need to be processed?
277 Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
278 mesh->bdr_attributes.Max() : 0);
279 bdr_attr_marker = 0;
280 for (int k = 0; k < boundary_face_integs.Size(); k++)
281 {
282 if (boundary_face_integs_marker[k] == NULL)
283 {
284 bdr_attr_marker = 1;
285 break;
286 }
287 Array<int> &bdr_marker = *boundary_face_integs_marker[k];
288 MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
289 "invalid boundary marker for boundary face integrator #"
290 << k << ", counting from zero");
291 for (int i = 0; i < bdr_attr_marker.Size(); i++)
292 {
293 bdr_attr_marker[i] |= bdr_marker[i];
294 }
295 }
296
297 for (int i = 0; i < mesh->GetNBE(); i++)
298 {
299 const int bdr_attr = mesh->GetBdrAttribute(i);
300 if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
301
302 tr = mesh->GetBdrFaceTransformations(i);
303 if (tr != NULL)
304 {
305 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
306 for (int k = 0; k < boundary_face_integs.Size(); k++)
307 {
309 (*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
310 { continue; }
311
313 AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
314 *tr, elemvect);
315 AddElementVector (vdofs, elemvect);
316 }
317 }
318 }
319 }
320
321 if (interior_face_integs.Size())
322 {
323 Mesh *mesh = fes->GetMesh();
324
325 for (int k = 0; k < interior_face_integs.Size(); k++)
326 {
327 for (int i = 0; i < mesh->GetNumFaces(); i++)
328 {
330 tr = mesh->GetInteriorFaceTransformations (i);
331 if (tr != NULL)
332 {
333 fes -> GetElementVDofs (tr -> Elem1No, vdofs);
334 Array<int> vdofs2;
335 fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
336 vdofs.Append(vdofs2);
338 AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
339 *fes->GetFE(tr->Elem2No),
340 *tr, elemvect);
341 AddElementVector (vdofs, elemvect);
342 }
343 }
344 }
345 }
346}
347
349{
351 if (ext) { ext->Update(); }
352}
353
355{
356 MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
357 fes = f;
358 v.UseDevice(true);
359 this->Vector::MakeRef(v, v_offset, fes->GetVSize());
361 if (ext) { ext->Update(); }
362}
363
365{
366 Update(f, v, v_offset);
367}
368
370{
371 if (domain_delta_integs.Size() == 0) { return; }
372
373 if (!HaveDeltaLocations())
374 {
375 int sdim = fes->GetMesh()->SpaceDimension();
376 Vector center;
377 DenseMatrix centers(sdim, domain_delta_integs.Size());
378 for (int i = 0; i < centers.Width(); i++)
379 {
380 centers.GetColumnReference(i, center);
381 domain_delta_integs[i]->GetDeltaCenter(center);
382 MFEM_VERIFY(center.Size() == sdim,
383 "Point dim " << center.Size() <<
384 " does not match space dim " << sdim);
385 }
388 }
389
390 Array<int> vdofs;
391 Vector elemvect;
392 for (int i = 0; i < domain_delta_integs.Size(); i++)
393 {
394 int elem_id = domain_delta_integs_elem_id[i];
395 // The delta center may be outside of this sub-domain, or
396 // (Par)Mesh::FindPoints() failed to find this point:
397 if (elem_id < 0) { continue; }
398
401 Trans.SetIntPoint(&ip);
402
403 fes->GetElementVDofs(elem_id, vdofs);
404 domain_delta_integs[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id),
405 Trans, elemvect);
406 AddElementVector(vdofs, elemvect);
407 }
408}
409
411{
412 Vector::operator=(value);
413 return *this;
414}
415
417{
418 MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
420 return *this;
421}
422
424{
425 if (!extern_lfs)
426 {
427 int k;
428 for (k=0; k < domain_delta_integs.Size(); k++)
429 { delete domain_delta_integs[k]; }
430 for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
431 for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
432 for (k=0; k < boundary_face_integs.Size(); k++)
433 { delete boundary_face_integs[k]; }
434 for (k=0; k < interior_face_integs.Size(); k++)
435 { delete interior_face_integs[k]; }
436 }
437
438 delete ext;
439}
440
441}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition: array.hpp:769
Abstract class for integrators that support delta coefficients.
Definition: lininteg.hpp:63
bool IsDelta() const
Returns true if the derived class instance uses a delta coefficient.
Definition: lininteg.hpp:85
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
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:93
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition: eltrans.hpp:484
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:581
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition: fespace.cpp:3204
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:773
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:749
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:287
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:3168
int GetNFbyType(FaceType type) const
Returns the number of faces according to the requested type.
Definition: fespace.hpp:757
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:559
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:713
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:25
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:25
void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Interior Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:103
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:45
Array< Array< int > * > boundary_face_integs_marker
Entries not owned.
Definition: linearform.hpp:63
void AssembleDelta()
Assembles delta functions of the linear form.
Definition: linearform.cpp:369
virtual void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
Make the LinearForm reference external data on a new FiniteElementSpace.
Definition: linearform.cpp:364
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:76
LinearForm()
Create an empty LinearForm without an associated FiniteElementSpace.
Definition: linearform.hpp:104
Array< LinearFormIntegrator * > boundary_face_integs
Set of Boundary Face Integrators to be applied.
Definition: linearform.hpp:62
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:162
Array< int > domain_delta_integs_elem_id
The element ids where the centers of the delta functions lie.
Definition: linearform.hpp:69
virtual bool SupportsDevice() const
Return true if assembly on device is supported, false otherwise.
Definition: linearform.cpp:108
LinearForm & operator=(const LinearForm &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition: linearform.hpp:121
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
Definition: linearform.hpp:59
bool HaveDeltaLocations()
If true, the delta locations are not (re)computed during assembly.
Definition: linearform.hpp:75
FiniteElementSpace * fes
FE space on which the LinearForm lives. Not owned.
Definition: linearform.hpp:30
Array< Array< int > * > domain_integs_marker
Definition: linearform.hpp:51
Array< LinearFormIntegrator * > interior_face_integs
Set of Internal Face Integrators to be applied.
Definition: linearform.hpp:66
LinearFormExtension * ext
Extension for supporting different assembly levels.
Definition: linearform.hpp:33
void ResetDeltaLocations()
Force (re)computation of delta locations.
Definition: linearform.hpp:79
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.cpp:348
Array< DeltaLFIntegrator * > domain_delta_integs
Separate array for integrators with delta function coefficients.
Definition: linearform.hpp:54
Array< LinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
Definition: linearform.hpp:45
~LinearForm()
Destroys linear form.
Definition: linearform.cpp:423
Array< IntegrationPoint > domain_delta_integs_ip
The reference coordinates where the centers of the delta functions lie.
Definition: linearform.hpp:72
Array< LinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
Definition: linearform.hpp:57
int extern_lfs
Indicates the LinearFormIntegrators stored in domain_integs, domain_delta_integs, boundary_integs,...
Definition: linearform.hpp:42
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:89
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:172
bool fast_assembly
Should we use the device-compatible fast assembly algorithm (false by default)
Definition: linearform.hpp:37
Mesh data type.
Definition: mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:282
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:290
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:6250
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1333
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1339
int GetBdrElementFaceIndex(int be_idx) const
Return the local face (codimension-1) index for the given boundary element index.
Definition: mesh.hpp:1518
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
FaceInformation GetFaceInformation(int f) const
Definition: mesh.cpp:1169
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1163
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition: mesh.cpp:1099
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
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:13081
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition: mesh.cpp:1079
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:280
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:80
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:670
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:139
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:136
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:538
Vector & operator=(const real_t *v)
Copy Size() entries from v.
Definition: vector.cpp:118
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:602
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:1345
float real_t
Definition: config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition: lor_mms.hpp:30