error_estimator.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_ERROR_ESTIMATOR_NAMESPACE_HEADER
27 #define OOMPH_ERROR_ESTIMATOR_NAMESPACE_HEADER
28 
29 #include "mesh.h"
30 #include "quadtree.h"
31 #include "nodes.h"
32 #include "algebraic_elements.h"
33 
34 namespace oomph
35 {
36  //========================================================================
37  /// Base class for spatial error estimators
38  //========================================================================
40  {
41  public:
42  /// Default empty constructor
44 
45  /// Broken copy constructor
46  ErrorEstimator(const ErrorEstimator&) = delete;
47 
48  /// Broken assignment operator
49  void operator=(const ErrorEstimator&) = delete;
50 
51  /// Empty virtual destructor
52  virtual ~ErrorEstimator() {}
53 
54  /// Compute the elemental error-measures for a given mesh
55  /// and store them in a vector.
56  void get_element_errors(Mesh*& mesh_pt, Vector<double>& elemental_error)
57  {
58  // Create dummy doc info object and switch off output
59  DocInfo doc_info;
60  doc_info.disable_doc();
61  // Forward call to version with doc.
62  get_element_errors(mesh_pt, elemental_error, doc_info);
63  }
64 
65 
66  /// Compute the elemental error measures for a given mesh
67  /// and store them in a vector. Doc errors etc.
68  virtual void get_element_errors(Mesh*& mesh_pt,
69  Vector<double>& elemental_error,
70  DocInfo& doc_info) = 0;
71  };
72 
73 
74  //========================================================================
75  /// Base class for finite elements that can compute the
76  /// quantities that are required for the Z2 error estimator.
77  //========================================================================
79  {
80  public:
81  /// Default empty constructor
83 
84  /// Broken copy constructor
86 
87  /// Broken assignment operator
88  void operator=(const ElementWithZ2ErrorEstimator&) = delete;
89 
90  /// Number of 'flux' terms for Z2 error estimation
91  virtual unsigned num_Z2_flux_terms() = 0;
92 
93  /// A stuitable error estimator for
94  /// a multi-physics elements may require one Z2 error estimate for each
95  /// field (e.g. velocity and temperature in a fluid convection problem).
96  /// It is assumed that these error estimates will each use
97  /// selected flux terms. The number of compound fluxes returns the number
98  /// of such combinations of the flux terms. Default value is one and all
99  /// flux terms are combined with equal weight.
100  virtual unsigned ncompound_fluxes()
101  {
102  return 1;
103  }
104 
105  /// Z2 'flux' terms for Z2 error estimation
106  virtual void get_Z2_flux(const Vector<double>& s, Vector<double>& flux) = 0;
107 
108  /// Plot the error when compared against a given exact flux.
109  /// Also calculates the norm of the error and that of the exact flux.
111  std::ostream& outfile,
113  double& error,
114  double& norm)
115  {
116  std::string error_message =
117  "compute_exact_Z2_error undefined for this element \n";
118  outfile << error_message;
119 
120  throw OomphLibError(
121  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
122  }
123 
124  /// Return the compound flux index of each flux component
125  /// The default (do nothing behaviour) will mean that all indices
126  /// remain at the default value zero.
127  virtual void get_Z2_compound_flux_indices(Vector<unsigned>& flux_index) {}
128 
129  /// Order of recovery shape functions
130  virtual unsigned nrecovery_order() = 0;
131 
132  /// Return the geometric jacobian (should be overloaded in
133  /// cylindrical and spherical geometries).
134  /// Default value one is suitable for Cartesian coordinates
135  virtual double geometric_jacobian(const Vector<double>& x)
136  {
137  return 1.0;
138  }
139  };
140 
141 
142  //========================================================================
143  /// Z2-error-estimator: Elements that can
144  /// be used with Z2 error estimation should be derived from
145  /// the base class ElementWithZ2ErrorEstimator and implement its
146  /// pure virtual member functions to provide the following functionality:
147  /// - pointer-based access to the vertex nodes in the element
148  /// (this is required to facilitate setup of element patches).
149  /// - the computation of a suitable "flux vector" which represents
150  /// a quantity whose FE representation is discontinuous across
151  /// element boundaries but would become continuous under infinite
152  /// mesh refinement.
153  ///
154  /// As an example consider the finite element solution of the Laplace problem,
155  /// \f$ \partial^2 u/\partial x_i^2 = 0 \f$. If we approximate the
156  /// unknown \f$ u \f$ on a finite element mesh with \f$ N \f$ nodes as
157  /// \f[ u^{[FE]}(x_k) = \sum_{j=1}^{N} U_j \ \psi_j(x_k), \f]
158  /// where the \f$ \psi_j(x_k) \f$ are the (global) \f$ C^0 \f$ basis
159  /// functions, the finite-element representation of the flux,
160  /// \f$ f_i = \partial u/\partial x_i \f$,
161  /// \f[ f_i^{[FE]} = \sum_{j=1}^{N} U_j \ \frac{\partial \psi_j}{\partial x_i} \f]
162  /// is discontinuous between elements but the magnitude of the jump
163  /// decreases under mesh refinement. We denote the number
164  /// of flux terms by \f$N_{flux}\f$, so for a 2D (3D) Laplace problem,
165  /// \f$N_{flux}=2 \ (3).\f$
166  ///
167  /// The idea behind Z2 error estimation is to compute an
168  /// approximation for the flux that is more accurate than the value
169  /// \f$ f_i^{[FE]} \f$ obtained directly from the finite element
170  /// solution. We refer to the improved approximation for the flux
171  /// as the "recovered flux" and denote it by \f$ f_i^{[rec]} \f$. Since
172  /// \f$ f_i^{[rec]} \f$ is more accurate than \f$ f_i^{[FE]} \f$, the
173  /// difference between the two provides a measure of the error.
174  /// In Z2 error estimation, the "recovered flux" is determined by
175  /// projecting the discontinuous, FE-based flux \f$ f_i^{[FE]} \f$
176  /// onto a set of continuous basis functions, known as the "recovery shape
177  /// functions". In principle, one could use the
178  /// finite element shape functions \f$ \psi_j(x_k) \f$ as the
179  /// recovery shape functions but then the determination of the
180  /// recovered flux would involve the solution of a linear system
181  /// that is as big as the original problem. To avoid this, the projection
182  /// is performed over small patches of elements within which
183  /// low-order polynomials (defined in terms of the global, Eulerian
184  /// coordinates) are used as the recovery shape functions.
185  ///
186  /// Z2 error estimation is known to be "optimal" for many self-adjoint
187  /// problems. See, e.g., Zienkiewicz & Zhu's original papers
188  /// "The superconvergent patch recovery and a posteriori error estimates."
189  /// International Journal for Numerical Methods in Engineering \b 33 (1992)
190  /// Part I: 1331-1364; Part II: 1365-1382.
191  /// In non-self adjoint problems, the error estimator only
192  /// analyses the "self-adjoint part" of the differential operator.
193  /// However, in many cases, this still produces good error indicators
194  /// since the error estimator detects under-resolved, sharp gradients
195  /// in the solution.
196  ///
197  /// Z2 error estimation works as follows:
198  /// -# We combine the elements in the mesh into a number of "patches",
199  /// which consist of all elements that share a common vertex node.
200  /// Most elements will therefore be members of multiple patches.
201  /// -# Within each patch \f$p\f$, we expand the "recovered flux" as
202  ///
203  /// \f[ f^{[rec](p)}_i(x_k) = \sum_{j=1}^{N_{rec}} F^{(p)}_{ij} \ \psi^{[rec]}_j(x_k) \mbox{ \ \ \ for $i=1,...,N_{flux}$,} \f]
204  /// where the functions \f$ \psi^{[rec]}_j(x_k)\f$
205  /// are the recovery shape functions, which are functions of the global,
206  /// Eulerian coordinates. Typically, these are chosen to be low-order
207  /// polynomials. For instance, in 2D, a bilinear representation of
208  /// \f$ f^{(p)}_i(x_0,x_1) \f$ involves the \f$N_{rec}=3\f$ recovery shape
209  /// functions \f$ \psi^{[rec]}_0(x_0,x_1)=1, \ \psi^{[rec]}_1(x_0,x_1)=x_0 \f$
210  /// and \f$ \psi^{[rec]}_2(x_0,x_1)=x_1\f$.
211  ///
212  /// We determine the coefficients \f$ F^{(p)}_{ij} \f$ by enforcing
213  /// \f$ f^{(p)}_i(x_k) = f^{[FE]}_i(x_k)\f$ in its weak form:
214  ///
215  /// \f[ \int_{\mbox{Patch $p$}} \left( f^{[FE]}_i(x_k) - \sum_{j=1}^{N_{rec}} F^{(p)}_{ij} \ \psi^{[rec]}_j(x_k) \right) \psi^{[rec]}_l(x_k)\ dv = 0 \mbox{ \ \ \ \ for $l=1,...,N_{rec}$ and $i=1,...,N_{flux}$}. \f]
216  /// Once the \f$ F^{(p)}_{ij} \f$ are determined in a given patch,
217  /// we determine the values of the recovered flux at
218  /// all nodes that are part of the patch. We denote the
219  /// value of the recovered flux at node \f$ k \f$ by
220  /// \f$ {\cal F}^{(p)}_{ik}\f$.
221  ///
222  /// We repeat this procedure for every patch. For nodes that are part of
223  /// multiple patches, the procedure
224  /// will provide multiple, slightly different nodal values for the
225  /// recovered flux. We average these values via
226  /// \f[ {\cal F}_{ij} = \frac{1}{N_p(j)} \sum_{\mbox{Node $j \in $ patch $p$}} {\cal F}^{(p)}_{ij}, \f]
227  /// where \f$N_p(j)\f$ denotes the number of patches that node \f$ j\f$ is
228  /// a member of. This allows us to obtain a globally-continuous,
229  /// finite-element based representation of the recovered flux as
230  /// \f[ f_i^{[rec]} = \sum_{j=1}^{N} {\cal F}_{ij}\ \psi_j, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1) \f]
231  /// where the \f$ \psi_j \f$ are the (global) finite element
232  /// shape functions.
233  /// -# Since the recovered flux in (1), is based on nodal values, we can
234  /// evaluate it locally within each of the \f$ N_e\f$ elements in the mesh
235  /// to obtain a normalised elemental error estimate via
236  ///
237  /// \f[ E_{e} = \sqrt{ \frac{ \int_{\mbox{$e$}} \left( f_i^{[rec]} - f_i^{[FE]} \right)^2 dv} {\sum_{e'=1}^{N_e} \int_{\mbox{$e'$}} \left( f_i^{[rec]} \right)^2 dv} } \mbox{\ \ \ for $e=1,...,N_e$.} \f]
238  /// In this (default) form, mesh refinement, based on this error estimate
239  /// will lead to an equidistribution of the error across all elements.
240  /// Usually, this is the desired behaviour. However, there are
241  /// cases in which the solution evolves towards a state in which
242  /// the flux tends to zero while the solution itself becomes so simple
243  /// that it can be represented exactly on any finite element mesh
244  /// (e.g. in spin-up problems in which the fluid motion approaches
245  /// a rigid body rotation). In that case the mesh adaptation tries
246  /// to equidistribute any small roundoff errors in the solution,
247  /// causing strong, spatially uniform mesh refinement throughout
248  /// the domain, even though the solution is already fully converged.
249  /// For such cases, it is possible to replace the denominator in the
250  /// above expression by a prescribed reference flux, which may be
251  /// specified via the access function
252  /// \code
253  /// reference_flux_norm()
254  /// \endcode
255  ///
256  ///
257  //========================================================================
258  class Z2ErrorEstimator : public virtual ErrorEstimator
259  {
260  public:
261  /// Function pointer to combined error estimator function
262  typedef double (*CombinedErrorEstimateFctPt)(const Vector<double>& errors);
263 
264  /// Constructor: Set order of recovery shape functions
268  Reference_flux_norm(0.0),
270  {
271  }
272 
273 
274  /// Constructor: Leave order of recovery shape functions open
275  /// for now -- they will be read out from first element of the mesh
276  /// when the error estimator is applied
278  : Recovery_order(0),
280  Reference_flux_norm(0.0),
282  {
283  }
284 
285  /// Broken copy constructor
287 
288  /// Broken assignment operator
289  void operator=(const Z2ErrorEstimator&) = delete;
290 
291  /// Empty virtual destructor
292  virtual ~Z2ErrorEstimator() {}
293 
294  /// Compute the elemental error measures for a given mesh
295  /// and store them in a vector.
296  void get_element_errors(Mesh*& mesh_pt, Vector<double>& elemental_error)
297  {
298  // Create dummy doc info object and switch off output
299  DocInfo doc_info;
300  doc_info.disable_doc();
301  // Forward call to version with doc.
302  get_element_errors(mesh_pt, elemental_error, doc_info);
303  }
304 
305  /// Compute the elemental error measures for a given mesh
306  /// and store them in a vector.
307  /// If doc_info.enable_doc(), doc FE and recovered fluxes in
308  /// - flux_fe*.dat
309  /// - flux_rec*.dat
310  void get_element_errors(Mesh*& mesh_pt,
311  Vector<double>& elemental_error,
312  DocInfo& doc_info);
313 
314  /// Access function for order of recovery polynomials
315  unsigned& recovery_order()
316  {
317  return Recovery_order;
318  }
319 
320  /// Access function for order of recovery polynomials (const version)
321  unsigned recovery_order() const
322  {
323  return Recovery_order;
324  }
325 
326  /// Access function: Pointer to combined error estimate function
328  {
329  return Combined_error_fct_pt;
330  }
331 
332  /// Access function: Pointer to combined error estimate function.
333  /// Const version
335  {
336  return Combined_error_fct_pt;
337  }
338 
339  /// Setup patches: For each vertex node pointed to by nod_pt,
340  /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
341  /// contains the pointers to the elements that the node is part of.
342  void setup_patches(Mesh*& mesh_pt,
344  adjacent_elements_pt,
345  Vector<Node*>& vertex_node_pt);
346 
347  /// Access function for prescribed reference flux norm
349  {
350  return Reference_flux_norm;
351  }
352 
353  /// Access function for prescribed reference flux norm (const. version)
354  double reference_flux_norm() const
355  {
356  return Reference_flux_norm;
357  }
358 
359  /// Return a combined error estimate from all compound errors
360  double get_combined_error_estimate(const Vector<double>& compound_error);
361 
362  private:
363  /// Given the vector of elements that make up a patch,
364  /// the number of recovery and flux terms, and the spatial
365  /// dimension of the problem, compute
366  /// the matrix of recovered flux coefficients and return
367  /// a pointer to it.
369  const Vector<ElementWithZ2ErrorEstimator*>& patch_el_pt,
370  const unsigned& num_recovery_terms,
371  const unsigned& num_flux_terms,
372  const unsigned& dim,
373  DenseMatrix<double>*& recovered_flux_coefficient_pt);
374 
375 
376  /// Return number of coefficients for expansion of recovered fluxes
377  /// for given spatial dimension of elements.
378  /// (We use complete polynomials of the specified given order.)
379  unsigned nrecovery_terms(const unsigned& dim);
380 
381 
382  /// Recovery shape functions as functions of the global, Eulerian
383  /// coordinate x of dimension dim.
384  /// The recovery shape functions are complete polynomials of
385  /// the order specified by Recovery_order.
386  void shape_rec(const Vector<double>& x,
387  const unsigned& dim,
388  Vector<double>& psi_r);
389 
390  /// Integation scheme associated with the recovery shape functions
391  /// must be of sufficiently high order to integrate the mass matrix
392  /// associated with the recovery shape functions
393  Integral* integral_rec(const unsigned& dim, const bool& is_q_mesh);
394 
395  /// Order of recovery polynomials
396  unsigned Recovery_order;
397 
398  /// Bool to indicate if recovery order is to be read out from
399  /// first element in mesh or set globally
401 
402  /// Doc flux and recovered flux
403  void doc_flux(Mesh* mesh_pt,
404  const unsigned& num_flux_terms,
406  const Vector<double>& elemental_error,
407  DocInfo& doc_info);
408 
409  /// Prescribed reference flux norm
411 
412  /// Function pointer to combined error estimator function
414  };
415 
416 
417  /// /////////////////////////////////////////////////////////////////////
418  /// /////////////////////////////////////////////////////////////////////
419  /// /////////////////////////////////////////////////////////////////////
420 
421 
422  //========================================================================
423  /// Dummy error estimator, allows manual specification of refinement
424  /// pattern by forcing refinement in regions defined by elements in
425  /// a reference mesh.
426  //========================================================================
427  class DummyErrorEstimator : public virtual ErrorEstimator
428  {
429  public:
430  /// Constructor. Provide mesh and number of the elements that define
431  /// the regions within which elements are to be refined subsequently.
432  /// Also specify the node number of a central node
433  /// within elements -- it's used to determine if an element is
434  /// in the region where refinement is supposed to take place.
435  /// Optional boolean flag (defaulting to false) indicates that
436  /// refinement decision is based on Lagrangian coordinates -- only
437  /// applicable to solid meshes.
439  const Vector<unsigned>& elements_to_refine,
440  const unsigned& central_node_number,
441  const bool& use_lagrangian_coordinates = false)
442  : Use_lagrangian_coordinates(use_lagrangian_coordinates),
443  Central_node_number(central_node_number)
444  {
445 #ifdef PARANOID
446 #ifdef OOMPH_HAS_MPI
447  if (mesh_pt->is_mesh_distributed())
448  {
449  throw OomphLibError(
450  "Can't use this error estimator on distributed meshes!",
451  OOMPH_CURRENT_FUNCTION,
452  OOMPH_EXCEPTION_LOCATION);
453  }
454 #endif
455 #endif
456 
457 #ifdef PARANOID
458  if (mesh_pt->nelement() == 0)
459  {
460  throw OomphLibError(
461  "Can't build error estimator if there are no elements in mesh\n",
462  OOMPH_CURRENT_FUNCTION,
463  OOMPH_EXCEPTION_LOCATION);
464  }
465 #endif
466 
467  unsigned dim = mesh_pt->finite_element_pt(0)->node_pt(0)->ndim();
468  if (use_lagrangian_coordinates)
469  {
470  SolidNode* solid_nod_pt =
471  dynamic_cast<SolidNode*>(mesh_pt->finite_element_pt(0)->node_pt(0));
472  if (solid_nod_pt != 0)
473  {
474  dim = solid_nod_pt->nlagrangian();
475  }
476  }
477  unsigned nregion = elements_to_refine.size();
478  Region_low_bound.resize(nregion);
479  Region_upp_bound.resize(nregion);
480  for (unsigned e = 0; e < nregion; e++)
481  {
482  Region_low_bound[e].resize(dim, 1.0e20);
483  Region_upp_bound[e].resize(dim, -1.0e20);
484  FiniteElement* el_pt =
485  mesh_pt->finite_element_pt(elements_to_refine[e]);
486  unsigned nnod = el_pt->nnode();
487  for (unsigned j = 0; j < nnod; j++)
488  {
489  Node* nod_pt = el_pt->node_pt(j);
490  for (unsigned i = 0; i < dim; i++)
491  {
492  double x = nod_pt->x(i);
493  if (use_lagrangian_coordinates)
494  {
495  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
496  if (solid_nod_pt != 0)
497  {
498  x = solid_nod_pt->xi(i);
499  }
500  }
501  if (x < Region_low_bound[e][i])
502  {
503  Region_low_bound[e][i] = x;
504  }
505  if (x > Region_upp_bound[e][i])
506  {
507  Region_upp_bound[e][i] = x;
508  }
509  }
510  }
511  }
512  }
513 
514 
515  /// Constructor. Provide vectors to "lower left" and "upper right"
516  /// vertices of refinement region
517  /// Also specify the node number of a central node
518  /// within elements -- it's used to determine if an element is
519  /// in the region where refinement is supposed to take place.
520  /// Optional boolean flag (defaulting to false) indicates that
521  /// refinement decision is based on Lagrangian coordinates -- only
522  /// applicable to solid meshes.
524  const Vector<double>& lower_left,
525  const Vector<double>& upper_right,
526  const unsigned& central_node_number,
527  const bool& use_lagrangian_coordinates = false)
528  : Use_lagrangian_coordinates(use_lagrangian_coordinates),
529  Central_node_number(central_node_number)
530  {
531 #ifdef PARANOID
532  if (mesh_pt->nelement() == 0)
533  {
534  throw OomphLibError(
535  "Can't build error estimator if there are no elements in mesh\n",
536  OOMPH_CURRENT_FUNCTION,
537  OOMPH_EXCEPTION_LOCATION);
538  }
539 #endif
540 
541  unsigned nregion = 1;
542  Region_low_bound.resize(nregion);
543  Region_upp_bound.resize(nregion);
544  Region_low_bound[0] = lower_left;
545  Region_upp_bound[0] = upper_right;
546  }
547 
548  /// Broken copy constructor
550 
551  /// Broken assignment operator
552  void operator=(const DummyErrorEstimator&) = delete;
553 
554  /// Empty virtual destructor
555  virtual ~DummyErrorEstimator() {}
556 
557 
558  /// Compute the elemental error measures for a given mesh
559  /// and store them in a vector. Doc errors etc.
560  virtual void get_element_errors(Mesh*& mesh_pt,
561  Vector<double>& elemental_error,
562  DocInfo& doc_info)
563  {
564 #ifdef PARANOID
565  if (doc_info.is_doc_enabled())
566  {
567  std::ostringstream warning_stream;
568  warning_stream
569  << "No output defined in DummyErrorEstimator::get_element_errors()\n"
570  << "Ignoring doc_info flag.\n";
571  OomphLibWarning(warning_stream.str(),
572  "DummyErrorEstimator::get_element_errors()",
573  OOMPH_EXCEPTION_LOCATION);
574  }
575 #endif
576  unsigned nregion = Region_low_bound.size();
577  unsigned nelem = mesh_pt->nelement();
578  for (unsigned e = 0; e < nelem; e++)
579  {
580  elemental_error[e] = 0.0;
581 
582  // Check if element is in the regions to be refined
583  // (based on coords of its central node)
584  Node* nod_pt =
586  for (unsigned r = 0; r < nregion; r++)
587  {
588  bool is_inside = true;
589  unsigned dim = Region_low_bound[r].size();
590  for (unsigned i = 0; i < dim; i++)
591  {
592  double x = nod_pt->x(i);
594  {
595  SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
596  if (solid_nod_pt != 0)
597  {
598  x = solid_nod_pt->xi(i);
599  }
600  }
601  if (x < Region_low_bound[r][i])
602  {
603  is_inside = false;
604  break;
605  }
606  if (x > Region_upp_bound[r][i])
607  {
608  is_inside = false;
609  break;
610  }
611  }
612  if (is_inside)
613  {
614  elemental_error[e] = 1.0;
615  break;
616  }
617  }
618  }
619  }
620 
621  private:
622  /// Use Lagrangian coordinates to decide which element is to be
623  /// refined
625 
626  /// Number of local node that is used to identify if an element
627  /// is located in the refinement region
629 
630  /// Upper bounds for the coordinates of the refinement regions
632 
633  /// Lower bounds for the coordinates of the refinement regions
635  };
636 
637 
638 } // namespace oomph
639 
640 #endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
Information for documentation of results: Directory and file number to enable output in the form RESL...
bool is_doc_enabled() const
Are we documenting?
void disable_doc()
Disable documentation.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
bool Use_lagrangian_coordinates
Use Lagrangian coordinates to decide which element is to be refined.
Vector< Vector< double > > Region_low_bound
Lower bounds for the coordinates of the refinement regions.
DummyErrorEstimator(Mesh *mesh_pt, const Vector< unsigned > &elements_to_refine, const unsigned &central_node_number, const bool &use_lagrangian_coordinates=false)
Constructor. Provide mesh and number of the elements that define the regions within which elements ar...
DummyErrorEstimator(const DummyErrorEstimator &)=delete
Broken copy constructor.
virtual void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)
Compute the elemental error measures for a given mesh and store them in a vector. Doc errors etc.
Vector< Vector< double > > Region_upp_bound
Upper bounds for the coordinates of the refinement regions.
unsigned Central_node_number
Number of local node that is used to identify if an element is located in the refinement region.
DummyErrorEstimator(Mesh *mesh_pt, const Vector< double > &lower_left, const Vector< double > &upper_right, const unsigned &central_node_number, const bool &use_lagrangian_coordinates=false)
Constructor. Provide vectors to "lower left" and "upper right" vertices of refinement region Also spe...
virtual ~DummyErrorEstimator()
Empty virtual destructor.
void operator=(const DummyErrorEstimator &)=delete
Broken assignment operator.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
ElementWithZ2ErrorEstimator(const ElementWithZ2ErrorEstimator &)=delete
Broken copy constructor.
void operator=(const ElementWithZ2ErrorEstimator &)=delete
Broken assignment operator.
virtual unsigned ncompound_fluxes()
A stuitable error estimator for a multi-physics elements may require one Z2 error estimate for each f...
virtual void get_Z2_compound_flux_indices(Vector< unsigned > &flux_index)
Return the compound flux index of each flux component The default (do nothing behaviour) will mean th...
ElementWithZ2ErrorEstimator()
Default empty constructor.
virtual void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)=0
Z2 'flux' terms for Z2 error estimation.
virtual void compute_exact_Z2_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
Plot the error when compared against a given exact flux. Also calculates the norm of the error and th...
virtual double geometric_jacobian(const Vector< double > &x)
Return the geometric jacobian (should be overloaded in cylindrical and spherical geometries)....
virtual unsigned num_Z2_flux_terms()=0
Number of 'flux' terms for Z2 error estimation.
virtual unsigned nrecovery_order()=0
Order of recovery shape functions.
Base class for spatial error estimators.
void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error)
Compute the elemental error-measures for a given mesh and store them in a vector.
ErrorEstimator()
Default empty constructor.
virtual ~ErrorEstimator()
Empty virtual destructor.
ErrorEstimator(const ErrorEstimator &)=delete
Broken copy constructor.
virtual void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error, DocInfo &doc_info)=0
Compute the elemental error measures for a given mesh and store them in a vector. Doc errors etc.
void operator=(const ErrorEstimator &)=delete
Broken assignment operator.
A general Finite Element class.
Definition: elements.h:1317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2214
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1763
Generic class for numerical integration schemes:
Definition: integral.h:49
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
A general mesh class.
Definition: mesh.h:67
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1883
unsigned nlagrangian() const
Return number of lagrangian coordinates.
Definition: nodes.h:1870
Z2-error-estimator: Elements that can be used with Z2 error estimation should be derived from the bas...
void get_recovered_flux_in_patch(const Vector< ElementWithZ2ErrorEstimator * > &patch_el_pt, const unsigned &num_recovery_terms, const unsigned &num_flux_terms, const unsigned &dim, DenseMatrix< double > *&recovered_flux_coefficient_pt)
Given the vector of elements that make up a patch, the number of recovery and flux terms,...
CombinedErrorEstimateFctPt & combined_error_fct_pt()
Access function: Pointer to combined error estimate function.
unsigned nrecovery_terms(const unsigned &dim)
Return number of coefficients for expansion of recovered fluxes for given spatial dimension of elemen...
Z2ErrorEstimator(const unsigned &recovery_order)
Constructor: Set order of recovery shape functions.
CombinedErrorEstimateFctPt Combined_error_fct_pt
Function pointer to combined error estimator function.
Z2ErrorEstimator()
Constructor: Leave order of recovery shape functions open for now – they will be read out from first ...
virtual ~Z2ErrorEstimator()
Empty virtual destructor.
double reference_flux_norm() const
Access function for prescribed reference flux norm (const. version)
double get_combined_error_estimate(const Vector< double > &compound_error)
Return a combined error estimate from all compound errors.
void setup_patches(Mesh *&mesh_pt, std::map< Node *, Vector< ElementWithZ2ErrorEstimator * > * > &adjacent_elements_pt, Vector< Node * > &vertex_node_pt)
Setup patches: For each vertex node pointed to by nod_pt, adjacent_elements_pt[nod_pt] contains the p...
void operator=(const Z2ErrorEstimator &)=delete
Broken assignment operator.
double(* CombinedErrorEstimateFctPt)(const Vector< double > &errors)
Function pointer to combined error estimator function.
double & reference_flux_norm()
Access function for prescribed reference flux norm.
unsigned Recovery_order
Order of recovery polynomials.
double Reference_flux_norm
Prescribed reference flux norm.
void get_element_errors(Mesh *&mesh_pt, Vector< double > &elemental_error)
Compute the elemental error measures for a given mesh and store them in a vector.
CombinedErrorEstimateFctPt combined_error_fct_pt() const
Access function: Pointer to combined error estimate function. Const version.
Z2ErrorEstimator(const Z2ErrorEstimator &)=delete
Broken copy constructor.
void shape_rec(const Vector< double > &x, const unsigned &dim, Vector< double > &psi_r)
Recovery shape functions as functions of the global, Eulerian coordinate x of dimension dim....
Integral * integral_rec(const unsigned &dim, const bool &is_q_mesh)
Integation scheme associated with the recovery shape functions must be of sufficiently high order to ...
unsigned recovery_order() const
Access function for order of recovery polynomials (const version)
void doc_flux(Mesh *mesh_pt, const unsigned &num_flux_terms, MapMatrixMixed< Node *, int, double > &rec_flux_map, const Vector< double > &elemental_error, DocInfo &doc_info)
Doc flux and recovered flux.
unsigned & recovery_order()
Access function for order of recovery polynomials.
bool Recovery_order_from_first_element
Bool to indicate if recovery order is to be read out from first element in mesh or set globally.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...