MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
estimators.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 "estimators.hpp"
13 
14 namespace mfem
15 {
16 
18 {
19  flux_space->Update(false);
20  // In parallel, 'flux' can be a GridFunction, as long as 'flux_space' is a
21  // ParFiniteElementSpace and 'solution' is a ParGridFunction.
23 
24  if (!anisotropic) { aniso_flags.SetSize(0); }
26  anisotropic ? &aniso_flags : NULL,
28  with_coeff);
29 
31 }
32 
33 
34 #ifdef MFEM_USE_MPI
35 
37 {
38  flux_space->Update(false);
39  smooth_flux_space->Update(false);
40 
41  // TODO: move these parameters in the class, and add Set* methods.
42  const double solver_tol = 1e-12;
43  const int solver_max_it = 200;
46  local_norm_p, solver_tol, solver_max_it);
47 
49 }
50 
51 #endif // MFEM_USE_MPI
52 
54  GridFunction& sol_,
55  FiniteElementSpace& flux_fespace_,
56  const Array<int> &attributes_)
57  : attributes(attributes_)
58  , flux_integrator(&di_)
59  , solution(&sol_)
60  , flux_space(&flux_fespace_)
61  , own_flux_fespace(false)
62 #ifdef MFEM_USE_MPI
63  , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
64 #endif // MFEM_USE_MPI
65 {
67 }
68 
70  GridFunction& sol_,
71  FiniteElementSpace* flux_fespace_,
72  const Array<int> &attributes_)
73  : attributes(attributes_)
74  , flux_integrator(&di_)
75  , solution(&sol_)
76  , flux_space(flux_fespace_)
77  , own_flux_fespace(true)
78 #ifdef MFEM_USE_MPI
79  , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
80 #endif // MFEM_USE_MPI
81 {
83 }
84 
86 {
87  if (own_flux_fespace)
88  {
89  delete flux_space;
90  }
91 }
92 
94 {
95  compute_element_coefficient = [](Mesh* mesh, const int e)
96  {
97  return 1.0;
98  };
99 
100  compute_face_coefficient = [](Mesh* mesh, const int f,
101  const bool shared_face)
102  {
103  auto FT = [&]()
104  {
105 #ifdef MFEM_USE_MPI
106  if (shared_face)
107  {
108  return dynamic_cast<ParMesh*>(mesh)->GetSharedFaceTransformations(f);
109  }
110 #endif // MFEM_USE_MPI
111  return mesh->GetFaceElementTransformations(f);
112  }();
113  const auto order = FT->GetFE()->GetOrder();
114 
115  // Poor man's face diameter.
116  double diameter = 0.0;
117 
118  Vector p1(mesh->SpaceDimension());
119  Vector p2(mesh->SpaceDimension());
120  // NOTE: We have no direct access to vertices for shared faces,
121  // so we fall back to compute the positions from the element.
122  // This can also be modified to compute the diameter for non-linear
123  // geometries by sampling along geometry-specific lines.
124  auto vtx_intrule = Geometries.GetVertices(FT->GetGeometryType());
125  const auto nip = vtx_intrule->GetNPoints();
126  for (int i = 0; i < nip; i++)
127  {
128  // Evaluate flux vector at integration point
129  auto fip1 = vtx_intrule->IntPoint(i);
130  FT->Transform(fip1, p1);
131 
132  for (int j = 0; j < nip; j++)
133  {
134  auto fip2 = vtx_intrule->IntPoint(j);
135  FT->Transform(fip2, p2);
136 
137  diameter = std::max<double>(diameter, p2.DistanceTo(p1));
138  }
139  }
140  return diameter/(2.0*order);
141  };
142 }
143 
144 void KellyErrorEstimator::ComputeEstimates()
145 {
146  // Remarks:
147  // For some context you may have to consult the documentation of
148  // the FaceInfo class [1]. Also, the FaceElementTransformations
149  // documentation [2] may be helpful to grasp what is going on. Note
150  // that the FaceElementTransformations also works in the non-
151  // conforming case to transfer the gauss points from the slave to
152  // the master element.
153  // [1]
154  // https://github.com/mfem/mfem/blob/02d0bfe9c18ce049c3c93a6a4208080fcfc96991/mesh/mesh.hpp#L94
155  // [2]
156  // https://github.com/mfem/mfem/blob/02d0bfe9c18ce049c3c93a6a4208080fcfc96991/fem/eltrans.hpp#L435
157 
158  flux_space->Update(false);
159 
160  auto xfes = solution->FESpace();
161  MFEM_ASSERT(xfes->GetVDim() == 1,
162  "Estimation for vector-valued problems not implemented yet.");
163  auto mesh = xfes->GetMesh();
164 
165  this->error_estimates.SetSize(xfes->GetNE());
166  this->error_estimates = 0.0;
167 
168  // 1. Compute fluxes in discontinuous space
169  GridFunction *flux =
170 #ifdef MFEM_USE_MPI
171  isParallel ? new ParGridFunction(dynamic_cast<ParFiniteElementSpace*>
172  (flux_space)) :
173 #endif // MFEM_USE_MPI
174  new GridFunction(flux_space);
175 
176  *flux = 0.0;
177 
178  // We pre-sort the array to speed up the search in the following loops.
179  if (attributes.Size())
180  {
181  attributes.Sort();
182  }
183 
184  Array<int> xdofs, fdofs;
185  Vector el_x, el_f;
186  for (int e = 0; e < xfes->GetNE(); e++)
187  {
188  auto attr = xfes->GetAttribute(e);
189  if (attributes.Size() && attributes.FindSorted(attr) == -1)
190  {
191  continue;
192  }
193 
194  xfes->GetElementVDofs(e, xdofs);
195  solution->GetSubVector(xdofs, el_x);
196 
197  ElementTransformation* Transf = xfes->GetElementTransformation(e);
198  flux_integrator->ComputeElementFlux(*xfes->GetFE(e), *Transf, el_x,
199  *flux_space->GetFE(e), el_f, true);
200 
201  flux_space->GetElementVDofs(e, fdofs);
202  flux->AddElementVector(fdofs, el_f);
203  }
204 
205  // 2. Add error contribution from local interior faces
206  for (int f = 0; f < mesh->GetNumFaces(); f++)
207  {
208  auto FT = mesh->GetFaceElementTransformations(f);
209 
210  auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * xfes->GetFaceOrder(f));
211  const auto nip = int_rule.GetNPoints();
212 
213  if (mesh->FaceIsInterior(f))
214  {
215  int Inf1, Inf2, NCFace;
216  mesh->GetFaceInfos(f, &Inf1, &Inf2, &NCFace);
217 
218  // Convention
219  // * Conforming face: Face side with smaller element id handles
220  // the integration
221  // * Non-conforming face: The slave handles the integration.
222  // See FaceInfo documentation for details.
223  bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
224  bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
225  if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
226  {
227  if (attributes.Size() &&
228  (attributes.FindSorted(FT->Elem1->Attribute) == -1
229  || attributes.FindSorted(FT->Elem2->Attribute) == -1))
230  {
231  continue;
232  }
233 
234  IntegrationRule eir;
235  Vector jumps(nip);
236 
237  // Integral over local half face on the side of e₁
238  // i.e. the numerical integration of ∫ flux ⋅ n dS₁
239  for (int i = 0; i < nip; i++)
240  {
241  // Evaluate flux at IP
242  auto &fip = int_rule.IntPoint(i);
243  IntegrationPoint ip;
244  FT->Loc1.Transform(fip, ip);
245 
246  Vector val(flux_space->GetVDim());
247  flux->GetVectorValue(FT->Elem1No, ip, val);
248 
249  // And build scalar product with normal
250  Vector normal(mesh->SpaceDimension());
251  FT->Face->SetIntPoint(&fip);
252  if (mesh->Dimension() == mesh->SpaceDimension())
253  {
254  CalcOrtho(FT->Face->Jacobian(), normal);
255  }
256  else
257  {
258  Vector ref_normal(mesh->Dimension());
259  FT->Loc1.Transf.SetIntPoint(&fip);
260  CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
261  auto &e1 = FT->GetElement1Transformation();
262  e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
263  normal /= e1.Weight();
264  }
265  jumps(i) = val * normal * fip.weight * FT->Face->Weight();
266  }
267 
268  // Subtract integral over half face of e₂
269  // i.e. the numerical integration of ∫ flux ⋅ n dS₂
270  for (int i = 0; i < nip; i++)
271  {
272  // Evaluate flux vector at IP
273  auto &fip = int_rule.IntPoint(i);
274  IntegrationPoint ip;
275  FT->Loc2.Transform(fip, ip);
276 
277  Vector val(flux_space->GetVDim());
278  flux->GetVectorValue(FT->Elem2No, ip, val);
279 
280  // And build scalar product with normal
281  Vector normal(mesh->SpaceDimension());
282  FT->Face->SetIntPoint(&fip);
283  if (mesh->Dimension() == mesh->SpaceDimension())
284  {
285  CalcOrtho(FT->Face->Jacobian(), normal);
286  }
287  else
288  {
289  Vector ref_normal(mesh->Dimension());
290  FT->Loc1.Transf.SetIntPoint(&fip);
291  CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
292  auto &e1 = FT->GetElement1Transformation();
293  e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
294  normal /= e1.Weight();
295  }
296 
297  jumps(i) -= val * normal * fip.weight * FT->Face->Weight();
298  }
299 
300  // Finalize "local" L₂ contribution
301  for (int i = 0; i < nip; i++)
302  {
303  jumps(i) *= jumps(i);
304  }
305  auto h_k_face = compute_face_coefficient(mesh, f, false);
306  double jump_integral = h_k_face*jumps.Sum();
307 
308  // A local face is shared between two local elements, so we
309  // can get away with integrating the jump only once and add
310  // it to both elements. To minimize communication, the jump
311  // of shared faces is computed locally by each process.
312  error_estimates(FT->Elem1No) += jump_integral;
313  error_estimates(FT->Elem2No) += jump_integral;
314  }
315  }
316  }
317 
318  current_sequence = solution->FESpace()->GetMesh()->GetSequence();
319 
320 #ifdef MFEM_USE_MPI
321  if (!isParallel)
322 #endif // MFEM_USE_MPI
323  {
324  // Finalize element errors
325  for (int e = 0; e < xfes->GetNE(); e++)
326  {
327  auto factor = compute_element_coefficient(mesh, e);
328  // The sqrt belongs to the norm and hₑ to the indicator.
329  error_estimates(e) = sqrt(factor * error_estimates(e));
330  }
331 
332  total_error = error_estimates.Norml2();
333  delete flux;
334  return;
335  }
336 
337 #ifdef MFEM_USE_MPI
338 
339  // 3. Add error contribution from shared interior faces
340  // Synchronize face data.
341 
342  ParGridFunction *pflux = dynamic_cast<ParGridFunction*>(flux);
343  MFEM_VERIFY(pflux, "flux is not a ParGridFunction pointer");
344 
345  ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
346  MFEM_VERIFY(pmesh, "mesh is not a ParMesh pointer");
347 
348  pflux->ExchangeFaceNbrData();
349 
350  for (int sf = 0; sf < pmesh->GetNSharedFaces(); sf++)
351  {
352  auto FT = pmesh->GetSharedFaceTransformations(sf, true);
353  if (attributes.Size() &&
354  (attributes.FindSorted(FT->Elem1->Attribute) == -1
355  || attributes.FindSorted(FT->Elem2->Attribute) == -1))
356  {
357  continue;
358  }
359 
360  auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * xfes->GetFaceOrder(0));
361  const auto nip = int_rule.GetNPoints();
362 
363  IntegrationRule eir;
364  Vector jumps(nip);
365 
366  // Integral over local half face on the side of e₁
367  // i.e. the numerical integration of ∫ flux ⋅ n dS₁
368  for (int i = 0; i < nip; i++)
369  {
370  // Evaluate flux vector at integration point
371  auto &fip = int_rule.IntPoint(i);
372  IntegrationPoint ip;
373  FT->Loc1.Transform(fip, ip);
374 
375  Vector val(flux_space->GetVDim());
376  flux->GetVectorValue(FT->Elem1No, ip, val);
377 
378  Vector normal(mesh->SpaceDimension());
379  FT->Face->SetIntPoint(&fip);
380  if (mesh->Dimension() == mesh->SpaceDimension())
381  {
382  CalcOrtho(FT->Face->Jacobian(), normal);
383  }
384  else
385  {
386  Vector ref_normal(mesh->Dimension());
387  FT->Loc1.Transf.SetIntPoint(&fip);
388  CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
389  auto &e1 = FT->GetElement1Transformation();
390  e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
391  normal /= e1.Weight();
392  }
393 
394  jumps(i) = val * normal * fip.weight * FT->Face->Weight();
395  }
396 
397  // Subtract integral over non-local half face of e₂
398  // i.e. the numerical integration of ∫ flux ⋅ n dS₂
399  for (int i = 0; i < nip; i++)
400  {
401  // Evaluate flux vector at integration point
402  auto &fip = int_rule.IntPoint(i);
403  IntegrationPoint ip;
404  FT->Loc2.Transform(fip, ip);
405 
406  Vector val(flux_space->GetVDim());
407  flux->GetVectorValue(FT->Elem2No, ip, val);
408 
409  // Evaluate gauss point
410  Vector normal(mesh->SpaceDimension());
411  FT->Face->SetIntPoint(&fip);
412  if (mesh->Dimension() == mesh->SpaceDimension())
413  {
414  CalcOrtho(FT->Face->Jacobian(), normal);
415  }
416  else
417  {
418  Vector ref_normal(mesh->Dimension());
419  CalcOrtho(FT->Loc1.Transf.Jacobian(), ref_normal);
420  auto &e1 = FT->GetElement1Transformation();
421  e1.AdjugateJacobian().MultTranspose(ref_normal, normal);
422  normal /= e1.Weight();
423  }
424 
425  jumps(i) -= val * normal * fip.weight * FT->Face->Weight();
426  }
427 
428  // Finalize "local" L₂ contribution
429  for (int i = 0; i < nip; i++)
430  {
431  jumps(i) *= jumps(i);
432  }
433  auto h_k_face = compute_face_coefficient(mesh, sf, true);
434  double jump_integral = h_k_face*jumps.Sum();
435 
436  error_estimates(FT->Elem1No) += jump_integral;
437  // We skip "error_estimates(FT->Elem2No) += jump_integral"
438  // because the error is stored on the remote process and
439  // recomputed there.
440  }
441  delete flux;
442 
443  // Finalize element errors
444  for (int e = 0; e < xfes->GetNE(); e++)
445  {
446  auto factor = compute_element_coefficient(mesh, e);
447  // The sqrt belongs to the norm and hₑ to the indicator.
448  error_estimates(e) = sqrt(factor * error_estimates(e));
449  }
450 
451  // Finish by computing the global error.
452  auto pfes = dynamic_cast<ParFiniteElementSpace*>(xfes);
453  MFEM_VERIFY(pfes, "xfes is not a ParFiniteElementSpace pointer");
454 
455  double process_local_error = pow(error_estimates.Norml2(),2.0);
456  MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
457  MPI_SUM, pfes->GetComm());
458  total_error = sqrt(total_error);
459 #endif // MFEM_USE_MPI
460 }
461 
463 {
464  MFEM_VERIFY(coef != NULL || vcoef != NULL,
465  "LpErrorEstimator has no coefficient! Call SetCoef first.");
466 
468  if (coef)
469  {
471  }
472  else
473  {
475  }
476 #ifdef MFEM_USE_MPI
478  auto pfes = dynamic_cast<ParFiniteElementSpace*>(sol->FESpace());
479  if (pfes)
480  {
481  auto process_local_error = total_error;
482  MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
483  MPI_SUM, pfes->GetComm());
484  }
485 #endif // MFEM_USE_MPI
488 }
489 
490 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void ComputeEstimates()
Compute the element error estimates.
Definition: estimators.cpp:36
virtual void ComputeElementLpErrors(const double p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:3356
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:465
ParFiniteElementSpace * flux_space
Ownership based on the flag own_flux_fes. Its Update() method is called automatically by this class w...
Definition: estimators.hpp:228
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
VectorCoefficient * vcoef
Definition: estimators.hpp:346
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3284
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
Abstract parallel finite element space.
Definition: pfespace.hpp:28
const FiniteElement * GetFE() const
Get the current element used to compute the transformations.
Definition: eltrans.hpp:389
GridFunction * sol
Definition: estimators.hpp:347
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2288
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:960
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:265
void ComputeEstimates()
Compute the element error estimates.
Definition: estimators.cpp:17
ParGridFunction * solution
Not owned.
Definition: estimators.hpp:226
Geometry Geometries
Definition: fe.cpp:49
long GetSequence() const
Definition: mesh.hpp:1639
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
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
void ResetCoefficientFunctions()
Change the coefficients back to default as described above.
Definition: estimators.cpp:93
GridFunction * solution
Not owned.
Definition: estimators.hpp:99
int local_norm_p
Local L_p norm to use, default is 1.
Definition: estimators.hpp:221
void Sort()
Sorts the array in ascending order. This requires operator&lt; to be defined for T.
Definition: array.hpp:248
int FindSorted(const T &el) const
Do bisection search for &#39;el&#39; in a sorted array; return -1 if not found.
Definition: array.hpp:820
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:557
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:652
KellyErrorEstimator(BilinearFormIntegrator &di_, GridFunction &sol_, FiniteElementSpace &flux_fes_, const Array< int > &attributes_=Array< int >())
Construct a new KellyErrorEstimator object for a scalar field.
Definition: estimators.cpp:53
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:626
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
int SpaceDimension() const
Definition: mesh.hpp:1000
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, const ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, double solver_tol, int solver_max_it)
Definition: pgridfunc.cpp:1137
ParFiniteElementSpace * smooth_flux_space
Ownership based on the flag own_flux_fes. Its Update() method is called automatically by this class w...
Definition: estimators.hpp:231
double ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains, bool with_coeff)
Definition: gridfunc.cpp:4093
FiniteElementSpace * flux_space
Ownership based on own_flux_fes. Its Update() method is called automatically by this class when neede...
Definition: estimators.hpp:101
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
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
BilinearFormIntegrator * integ
Not owned.
Definition: estimators.hpp:225
void ComputeEstimates()
Compute the element error estimates.
Definition: estimators.cpp:462
Vector data type.
Definition: vector.hpp:60
virtual void ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans, Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef=true)
Virtual method required for Zienkiewicz-Zhu type error estimators.
Definition: bilininteg.hpp:217
Class for parallel grid function.
Definition: pgridfunc.hpp:32
BilinearFormIntegrator * integ
Not owned.
Definition: estimators.hpp:98
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:900