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