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