MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
linearform.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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  // Linear forms are stored on the device
23  UseDevice(true);
24 
25  fes = f;
26  extern_lfs = 1;
27 
28  // Copy the pointers to the integrators
29  dlfi = lf->dlfi;
30 
31  dlfi_delta = lf->dlfi_delta;
32 
33  blfi = lf->blfi;
34 
35  flfi = lf->flfi;
37 }
38 
40 {
41  DeltaLFIntegrator *maybe_delta =
42  dynamic_cast<DeltaLFIntegrator *>(lfi);
43  if (!maybe_delta || !maybe_delta->IsDelta())
44  {
45  dlfi.Append(lfi);
46  }
47  else
48  {
49  dlfi_delta.Append(maybe_delta);
50  }
51 }
52 
54 {
55  blfi.Append (lfi);
56  blfi_marker.Append(NULL); // NULL -> all attributes are active
57 }
58 
60  Array<int> &bdr_attr_marker)
61 {
62  blfi.Append (lfi);
63  blfi_marker.Append(&bdr_attr_marker);
64 }
65 
67 {
68  flfi.Append(lfi);
69  flfi_marker.Append(NULL); // NULL -> all attributes are active
70 }
71 
73  Array<int> &bdr_attr_marker)
74 {
75  flfi.Append(lfi);
76  flfi_marker.Append(&bdr_attr_marker);
77 }
78 
80 {
81  Array<int> vdofs;
82  ElementTransformation *eltrans;
83  Vector elemvect;
84 
85  int i;
86 
87  Vector::operator=(0.0);
88 
89  // The above operation is executed on device because of UseDevice().
90  // The first use of AddElementVector() below will move it back to host
91  // because both 'vdofs' and 'elemvect' are on host.
92 
93  if (dlfi.Size())
94  {
95  for (i = 0; i < fes -> GetNE(); i++)
96  {
97  fes -> GetElementVDofs (i, vdofs);
98  eltrans = fes -> GetElementTransformation (i);
99  for (int k=0; k < dlfi.Size(); k++)
100  {
101  dlfi[k]->AssembleRHSElementVect(*fes->GetFE(i), *eltrans, elemvect);
102  AddElementVector (vdofs, elemvect);
103  }
104  }
105  }
106  AssembleDelta();
107 
108  if (blfi.Size())
109  {
110  Mesh *mesh = fes->GetMesh();
111 
112  // Which boundary attributes need to be processed?
113  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
114  mesh->bdr_attributes.Max() : 0);
115  bdr_attr_marker = 0;
116  for (int k = 0; k < blfi.Size(); k++)
117  {
118  if (blfi_marker[k] == NULL)
119  {
120  bdr_attr_marker = 1;
121  break;
122  }
123  Array<int> &bdr_marker = *blfi_marker[k];
124  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
125  "invalid boundary marker for boundary integrator #"
126  << k << ", counting from zero");
127  for (int i = 0; i < bdr_attr_marker.Size(); i++)
128  {
129  bdr_attr_marker[i] |= bdr_marker[i];
130  }
131  }
132 
133  for (i = 0; i < fes -> GetNBE(); i++)
134  {
135  const int bdr_attr = mesh->GetBdrAttribute(i);
136  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
137  fes -> GetBdrElementVDofs (i, vdofs);
138  eltrans = fes -> GetBdrElementTransformation (i);
139  for (int k=0; k < blfi.Size(); k++)
140  {
141  if (blfi_marker[k] &&
142  (*blfi_marker[k])[bdr_attr-1] == 0) { continue; }
143 
144  blfi[k]->AssembleRHSElementVect(*fes->GetBE(i), *eltrans, elemvect);
145 
146  AddElementVector (vdofs, elemvect);
147  }
148  }
149  }
150  if (flfi.Size())
151  {
153  Mesh *mesh = fes->GetMesh();
154 
155  // Which boundary attributes need to be processed?
156  Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
157  mesh->bdr_attributes.Max() : 0);
158  bdr_attr_marker = 0;
159  for (int k = 0; k < flfi.Size(); k++)
160  {
161  if (flfi_marker[k] == NULL)
162  {
163  bdr_attr_marker = 1;
164  break;
165  }
166  Array<int> &bdr_marker = *flfi_marker[k];
167  MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
168  "invalid boundary marker for boundary face integrator #"
169  << k << ", counting from zero");
170  for (int i = 0; i < bdr_attr_marker.Size(); i++)
171  {
172  bdr_attr_marker[i] |= bdr_marker[i];
173  }
174  }
175 
176  for (i = 0; i < mesh->GetNBE(); i++)
177  {
178  const int bdr_attr = mesh->GetBdrAttribute(i);
179  if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
180 
181  tr = mesh->GetBdrFaceTransformations(i);
182  if (tr != NULL)
183  {
184  fes -> GetElementVDofs (tr -> Elem1No, vdofs);
185  for (int k = 0; k < flfi.Size(); k++)
186  {
187  if (flfi_marker[k] &&
188  (*flfi_marker[k])[bdr_attr-1] == 0) { continue; }
189 
190  flfi[k] -> AssembleRHSElementVect (*fes->GetFE(tr -> Elem1No),
191  *tr, elemvect);
192  AddElementVector (vdofs, elemvect);
193  }
194  }
195  }
196  }
197 }
198 
200 {
201  fes = f;
202  NewMemoryAndSize(Memory<double>(v.GetMemory(), v_offset, f->GetVSize()),
203  f->GetVSize(), false);
205 }
206 
208 {
209  MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
210  fes = f;
211  v.UseDevice(true);
212  this->Vector::MakeRef(v, v_offset, fes->GetVSize());
213 }
214 
216 {
217  if (dlfi_delta.Size() == 0) { return; }
218 
219  if (!HaveDeltaLocations())
220  {
221  int sdim = fes->GetMesh()->SpaceDimension();
222  Vector center;
223  DenseMatrix centers(sdim, dlfi_delta.Size());
224  for (int i = 0; i < centers.Width(); i++)
225  {
226  centers.GetColumnReference(i, center);
227  dlfi_delta[i]->GetDeltaCenter(center);
228  MFEM_VERIFY(center.Size() == sdim,
229  "Point dim " << center.Size() <<
230  " does not match space dim " << sdim);
231  }
233  }
234 
235  Array<int> vdofs;
236  Vector elemvect;
237  for (int i = 0; i < dlfi_delta.Size(); i++)
238  {
239  int elem_id = dlfi_delta_elem_id[i];
240  // The delta center may be outside of this sub-domain, or
241  // (Par)Mesh::FindPoints() failed to find this point:
242  if (elem_id < 0) { continue; }
243 
244  const IntegrationPoint &ip = dlfi_delta_ip[i];
246  Trans.SetIntPoint(&ip);
247 
248  fes->GetElementVDofs(elem_id, vdofs);
249  dlfi_delta[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id), Trans,
250  elemvect);
251  AddElementVector(vdofs, elemvect);
252  }
253 }
254 
256 {
257  Vector::operator=(value);
258  return *this;
259 }
260 
262 {
263  MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
265  return *this;
266 }
267 
269 {
270  if (!extern_lfs)
271  {
272  int k;
273  for (k=0; k < dlfi_delta.Size(); k++) { delete dlfi_delta[k]; }
274  for (k=0; k < dlfi.Size(); k++) { delete dlfi[k]; }
275  for (k=0; k < blfi.Size(); k++) { delete blfi[k]; }
276  for (k=0; k < flfi.Size(); k++) { delete flfi[k]; }
277  }
278 }
279 
280 }
Array< LinearFormIntegrator * > blfi
Set of Boundary Integrators to be applied.
Definition: linearform.hpp:40
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:400
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1053
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:740
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:173
virtual void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
Make the LinearForm reference external data on a new FiniteElementSpace.
Definition: linearform.cpp:207
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
int extern_lfs
Indicates the LinearFormIntegrators stored in dlfi, dlfi_delta, blfi, and flfi are owned by another L...
Definition: linearform.hpp:31
void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:89
LinearForm & operator=(const LinearForm &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition: linearform.hpp:98
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
Array< Array< int > * > blfi_marker
Entries are not owned.
Definition: linearform.hpp:41
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
void ResetDeltaLocations()
Force (re)computation of delta locations.
Definition: linearform.hpp:57
Array< DeltaLFIntegrator * > dlfi_delta
Separate array for integrators with delta function coefficients.
Definition: linearform.hpp:37
Array< IntegrationPoint > dlfi_delta_ip
The reference coordinates where the centers of the delta functions lie.
Definition: linearform.hpp:51
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:183
Abstract base class LinearFormIntegrator.
Definition: lininteg.hpp:22
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:120
bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:92
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:314
Array< LinearFormIntegrator * > dlfi
Set of Domain Integrators to be applied.
Definition: linearform.hpp:34
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
FiniteElementSpace * fes
FE space on which the LinearForm lives. Not owned.
Definition: linearform.hpp:27
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1003
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:53
bool IsDelta() const
Returns true if the derived class instance uses a delta coefficient.
Definition: lininteg.hpp:70
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:587
int SpaceDimension() const
Definition: mesh.hpp:789
LinearForm()
Create an empty LinearForm without an associated FiniteElementSpace.
Definition: linearform.hpp:82
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:460
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:270
Abstract class for integrators that support delta coefficients.
Definition: lininteg.hpp:47
void AssembleDelta()
Assembles delta functions of the linear form.
Definition: linearform.cpp:215
void Update()
Update the object according to the associated FE space fes.
Definition: linearform.hpp:166
Class for integration point with weight.
Definition: intrules.hpp:25
~LinearForm()
Destroys linear form.
Definition: linearform.cpp:268
Array< int > dlfi_delta_elem_id
The element ids where the centers of the delta functions lie.
Definition: linearform.hpp:48
void NewMemoryAndSize(const Memory< double > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition: vector.hpp:508
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:1798
bool HaveDeltaLocations()
If true, the delta locations are not (re)computed during assembly.
Definition: linearform.hpp:54
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:10561
Vector data type.
Definition: vector.hpp:51
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:517
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
Array< LinearFormIntegrator * > flfi
Set of Boundary Face Integrators to be applied.
Definition: linearform.hpp:44
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:2047
Array< Array< int > * > flfi_marker
Entries are not owned.
Definition: linearform.hpp:45
double f(const Vector &p)