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-2022 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
34namespace oomph
35{
36 //========================================================================
37 /// Base class for spatial error estimators
38 //========================================================================
40 {
41 public:
42 /// Default empty constructor
44
45 /// Broken copy constructor
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
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.
128
129 /// Number of vertex nodes in the element
130 virtual unsigned nvertex_node() const = 0;
131
132 /// Pointer to the j-th vertex node in the element. Needed for
133 /// efficient patch assmbly
134 virtual Node* vertex_node_pt(const unsigned& j) const = 0;
135
136 /// Order of recovery shape functions
137 virtual unsigned nrecovery_order() = 0;
138
139 /// Return the geometric jacobian (should be overloaded in
140 /// cylindrical and spherical geometries).
141 /// Default value one is suitable for Cartesian coordinates
142 virtual double geometric_jacobian(const Vector<double>& x)
143 {
144 return 1.0;
145 }
146 };
147
148
149 //========================================================================
150 /// Z2-error-estimator: Elements that can
151 /// be used with Z2 error estimation should be derived from
152 /// the base class ElementWithZ2ErrorEstimator and implement its
153 /// pure virtual member functions to provide the following functionality:
154 /// - pointer-based access to the vertex nodes in the element
155 /// (this is required to facilitate setup of element patches).
156 /// - the computation of a suitable "flux vector" which represents
157 /// a quantity whose FE representation is discontinuous across
158 /// element boundaries but would become continuous under infinite
159 /// mesh refinement.
160 ///
161 /// As an example consider the finite element solution of the Laplace problem,
162 /// \f$ \partial^2 u/\partial x_i^2 = 0 \f$. If we approximate the
163 /// unknown \f$ u \f$ on a finite element mesh with \f$ N \f$ nodes as
164 /// \f[ u^{[FE]}(x_k) = \sum_{j=1}^{N} U_j \ \psi_j(x_k), \f]
165 /// where the \f$ \psi_j(x_k) \f$ are the (global) \f$ C^0 \f$ basis
166 /// functions, the finite-element representation of the flux,
167 /// \f$ f_i = \partial u/\partial x_i \f$,
168 /// \f[ f_i^{[FE]} = \sum_{j=1}^{N} U_j \ \frac{\partial \psi_j}{\partial x_i} \f]
169 /// is discontinuous between elements but the magnitude of the jump
170 /// decreases under mesh refinement. We denote the number
171 /// of flux terms by \f$N_{flux}\f$, so for a 2D (3D) Laplace problem,
172 /// \f$N_{flux}=2 \ (3).\f$
173 ///
174 /// The idea behind Z2 error estimation is to compute an
175 /// approximation for the flux that is more accurate than the value
176 /// \f$ f_i^{[FE]} \f$ obtained directly from the finite element
177 /// solution. We refer to the improved approximation for the flux
178 /// as the "recovered flux" and denote it by \f$ f_i^{[rec]} \f$. Since
179 /// \f$ f_i^{[rec]} \f$ is more accurate than \f$ f_i^{[FE]} \f$, the
180 /// difference between the two provides a measure of the error.
181 /// In Z2 error estimation, the "recovered flux" is determined by
182 /// projecting the discontinuous, FE-based flux \f$ f_i^{[FE]} \f$
183 /// onto a set of continuous basis functions, known as the "recovery shape
184 /// functions". In principle, one could use the
185 /// finite element shape functions \f$ \psi_j(x_k) \f$ as the
186 /// recovery shape functions but then the determination of the
187 /// recovered flux would involve the solution of a linear system
188 /// that is as big as the original problem. To avoid this, the projection
189 /// is performed over small patches of elements within which
190 /// low-order polynomials (defined in terms of the global, Eulerian
191 /// coordinates) are used as the recovery shape functions.
192 ///
193 /// Z2 error estimation is known to be "optimal" for many self-adjoint
194 /// problems. See, e.g., Zienkiewicz & Zhu's original papers
195 /// "The superconvergent patch recovery and a posteriori error estimates."
196 /// International Journal for Numerical Methods in Engineering \b 33 (1992)
197 /// Part I: 1331-1364; Part II: 1365-1382.
198 /// In non-self adjoint problems, the error estimator only
199 /// analyses the "self-adjoint part" of the differential operator.
200 /// However, in many cases, this still produces good error indicators
201 /// since the error estimator detects under-resolved, sharp gradients
202 /// in the solution.
203 ///
204 /// Z2 error estimation works as follows:
205 /// -# We combine the elements in the mesh into a number of "patches",
206 /// which consist of all elements that share a common vertex node.
207 /// Most elements will therefore be members of multiple patches.
208 /// -# Within each patch \f$p\f$, we expand the "recovered flux" as
209 ///
210 /// \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]
211 /// where the functions \f$ \psi^{[rec]}_j(x_k)\f$
212 /// are the recovery shape functions, which are functions of the global,
213 /// Eulerian coordinates. Typically, these are chosen to be low-order
214 /// polynomials. For instance, in 2D, a bilinear representation of
215 /// \f$ f^{(p)}_i(x_0,x_1) \f$ involves the \f$N_{rec}=3\f$ recovery shape
216 /// functions \f$ \psi^{[rec]}_0(x_0,x_1)=1, \ \psi^{[rec]}_1(x_0,x_1)=x_0 \f$
217 /// and \f$ \psi^{[rec]}_2(x_0,x_1)=x_1\f$.
218 ///
219 /// We determine the coefficients \f$ F^{(p)}_{ij} \f$ by enforcing
220 /// \f$ f^{(p)}_i(x_k) = f^{[FE]}_i(x_k)\f$ in its weak form:
221 ///
222 /// \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]
223 /// Once the \f$ F^{(p)}_{ij} \f$ are determined in a given patch,
224 /// we determine the values of the recovered flux at
225 /// all nodes that are part of the patch. We denote the
226 /// value of the recovered flux at node \f$ k \f$ by
227 /// \f$ {\cal F}^{(p)}_{ik}\f$.
228 ///
229 /// We repeat this procedure for every patch. For nodes that are part of
230 /// multiple patches, the procedure
231 /// will provide multiple, slightly different nodal values for the
232 /// recovered flux. We average these values via
233 /// \f[ {\cal F}_{ij} = \frac{1}{N_p(j)} \sum_{\mbox{Node $j \in $ patch $p$}} {\cal F}^{(p)}_{ij}, \f]
234 /// where \f$N_p(j)\f$ denotes the number of patches that node \f$ j\f$ is
235 /// a member of. This allows us to obtain a globally-continuous,
236 /// finite-element based representation of the recovered flux as
237 /// \f[ f_i^{[rec]} = \sum_{j=1}^{N} {\cal F}_{ij}\ \psi_j, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1) \f]
238 /// where the \f$ \psi_j \f$ are the (global) finite element
239 /// shape functions.
240 /// -# Since the recovered flux in (1), is based on nodal values, we can
241 /// evaluate it locally within each of the \f$ N_e\f$ elements in the mesh
242 /// to obtain a normalised elemental error estimate via
243 ///
244 /// \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]
245 /// In this (default) form, mesh refinement, based on this error estimate
246 /// will lead to an equidistribution of the error across all elements.
247 /// Usually, this is the desired behaviour. However, there are
248 /// cases in which the solution evolves towards a state in which
249 /// the flux tends to zero while the solution itself becomes so simple
250 /// that it can be represented exactly on any finite element mesh
251 /// (e.g. in spin-up problems in which the fluid motion approaches
252 /// a rigid body rotation). In that case the mesh adaptation tries
253 /// to equidistribute any small roundoff errors in the solution,
254 /// causing strong, spatially uniform mesh refinement throughout
255 /// the domain, even though the solution is already fully converged.
256 /// For such cases, it is possible to replace the denominator in the
257 /// above expression by a prescribed reference flux, which may be
258 /// specified via the access function
259 /// \code
260 /// reference_flux_norm()
261 /// \endcode
262 ///
263 ///
264 //========================================================================
265 class Z2ErrorEstimator : public virtual ErrorEstimator
266 {
267 public:
268 /// Function pointer to combined error estimator function
269 typedef double (*CombinedErrorEstimateFctPt)(const Vector<double>& errors);
270
271 /// Constructor: Set order of recovery shape functions
277 {
278 }
279
280
281 /// Constructor: Leave order of recovery shape functions open
282 /// for now -- they will be read out from first element of the mesh
283 /// when the error estimator is applied
285 : Recovery_order(0),
289 {
290 }
291
292 /// Broken copy constructor
294
295 /// Broken assignment operator
296 void operator=(const Z2ErrorEstimator&) = delete;
297
298 /// Empty virtual destructor
299 virtual ~Z2ErrorEstimator() {}
300
301 /// Compute the elemental error measures for a given mesh
302 /// and store them in a vector.
303 void get_element_errors(Mesh*& mesh_pt, Vector<double>& elemental_error)
304 {
305 // Create dummy doc info object and switch off output
306 DocInfo doc_info;
307 doc_info.disable_doc();
308 // Forward call to version with doc.
309 get_element_errors(mesh_pt, elemental_error, doc_info);
310 }
311
312 /// Compute the elemental error measures for a given mesh
313 /// and store them in a vector.
314 /// If doc_info.enable_doc(), doc FE and recovered fluxes in
315 /// - flux_fe*.dat
316 /// - flux_rec*.dat
317 void get_element_errors(Mesh*& mesh_pt,
318 Vector<double>& elemental_error,
319 DocInfo& doc_info);
320
321 /// Access function for order of recovery polynomials
322 unsigned& recovery_order()
323 {
324 return Recovery_order;
325 }
326
327 /// Access function for order of recovery polynomials (const version)
328 unsigned recovery_order() const
329 {
330 return Recovery_order;
331 }
332
333 /// Access function: Pointer to combined error estimate function
335 {
337 }
338
339 /// Access function: Pointer to combined error estimate function.
340 /// Const version
342 {
344 }
345
346 /// Setup patches: For each vertex node pointed to by nod_pt,
347 /// adjacent_elements_pt[nod_pt] contains the pointer to the vector that
348 /// contains the pointers to the elements that the node is part of.
349 void setup_patches(Mesh*& mesh_pt,
351 adjacent_elements_pt,
352 Vector<Node*>& vertex_node_pt);
353
354 /// Access function for prescribed reference flux norm
356 {
357 return Reference_flux_norm;
358 }
359
360 /// Access function for prescribed reference flux norm (const. version)
361 double reference_flux_norm() const
362 {
363 return Reference_flux_norm;
364 }
365
366 /// Return a combined error estimate from all compound errors
367 double get_combined_error_estimate(const Vector<double>& compound_error);
368
369 private:
370 /// Given the vector of elements that make up a patch,
371 /// the number of recovery and flux terms, and the spatial
372 /// dimension of the problem, compute
373 /// the matrix of recovered flux coefficients and return
374 /// a pointer to it.
376 const Vector<ElementWithZ2ErrorEstimator*>& patch_el_pt,
377 const unsigned& num_recovery_terms,
378 const unsigned& num_flux_terms,
379 const unsigned& dim,
380 DenseMatrix<double>*& recovered_flux_coefficient_pt);
381
382
383 /// Return number of coefficients for expansion of recovered fluxes
384 /// for given spatial dimension of elements.
385 /// (We use complete polynomials of the specified given order.)
386 unsigned nrecovery_terms(const unsigned& dim);
387
388
389 /// Recovery shape functions as functions of the global, Eulerian
390 /// coordinate x of dimension dim.
391 /// The recovery shape functions are complete polynomials of
392 /// the order specified by Recovery_order.
393 void shape_rec(const Vector<double>& x,
394 const unsigned& dim,
395 Vector<double>& psi_r);
396
397 /// Integation scheme associated with the recovery shape functions
398 /// must be of sufficiently high order to integrate the mass matrix
399 /// associated with the recovery shape functions
400 Integral* integral_rec(const unsigned& dim, const bool& is_q_mesh);
401
402 /// Order of recovery polynomials
404
405 /// Bool to indicate if recovery order is to be read out from
406 /// first element in mesh or set globally
408
409 /// Doc flux and recovered flux
410 void doc_flux(Mesh* mesh_pt,
411 const unsigned& num_flux_terms,
413 const Vector<double>& elemental_error,
414 DocInfo& doc_info);
415
416 /// Prescribed reference flux norm
418
419 /// Function pointer to combined error estimator function
421 };
422
423
424 /// /////////////////////////////////////////////////////////////////////
425 /// /////////////////////////////////////////////////////////////////////
426 /// /////////////////////////////////////////////////////////////////////
427
428
429 //========================================================================
430 /// Dummy error estimator, allows manual specification of refinement
431 /// pattern by forcing refinement in regions defined by elements in
432 /// a reference mesh.
433 //========================================================================
434 class DummyErrorEstimator : public virtual ErrorEstimator
435 {
436 public:
437 /// Constructor. Provide mesh and number of the elements that define
438 /// the regions within which elements are to be refined subsequently.
439 /// Also specify the node number of a central node
440 /// within elements -- it's used to determine if an element is
441 /// in the region where refinement is supposed to take place.
442 /// Optional boolean flag (defaulting to false) indicates that
443 /// refinement decision is based on Lagrangian coordinates -- only
444 /// applicable to solid meshes.
446 const Vector<unsigned>& elements_to_refine,
447 const unsigned& central_node_number,
448 const bool& use_lagrangian_coordinates = false)
449 : Use_lagrangian_coordinates(use_lagrangian_coordinates),
450 Central_node_number(central_node_number)
451 {
452#ifdef PARANOID
453#ifdef OOMPH_HAS_MPI
454 if (mesh_pt->is_mesh_distributed())
455 {
456 throw OomphLibError(
457 "Can't use this error estimator on distributed meshes!",
458 OOMPH_CURRENT_FUNCTION,
459 OOMPH_EXCEPTION_LOCATION);
460 }
461#endif
462#endif
463
464#ifdef PARANOID
465 if (mesh_pt->nelement() == 0)
466 {
467 throw OomphLibError(
468 "Can't build error estimator if there are no elements in mesh\n",
469 OOMPH_CURRENT_FUNCTION,
470 OOMPH_EXCEPTION_LOCATION);
471 }
472#endif
473
474 unsigned dim = mesh_pt->finite_element_pt(0)->node_pt(0)->ndim();
475 if (use_lagrangian_coordinates)
476 {
477 SolidNode* solid_nod_pt =
478 dynamic_cast<SolidNode*>(mesh_pt->finite_element_pt(0)->node_pt(0));
479 if (solid_nod_pt != 0)
480 {
481 dim = solid_nod_pt->nlagrangian();
482 }
483 }
484 unsigned nregion = elements_to_refine.size();
485 Region_low_bound.resize(nregion);
486 Region_upp_bound.resize(nregion);
487 for (unsigned e = 0; e < nregion; e++)
488 {
489 Region_low_bound[e].resize(dim, 1.0e20);
490 Region_upp_bound[e].resize(dim, -1.0e20);
491 FiniteElement* el_pt =
492 mesh_pt->finite_element_pt(elements_to_refine[e]);
493 unsigned nnod = el_pt->nnode();
494 for (unsigned j = 0; j < nnod; j++)
495 {
496 Node* nod_pt = el_pt->node_pt(j);
497 for (unsigned i = 0; i < dim; i++)
498 {
499 double x = nod_pt->x(i);
500 if (use_lagrangian_coordinates)
501 {
502 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
503 if (solid_nod_pt != 0)
504 {
505 x = solid_nod_pt->xi(i);
506 }
507 }
508 if (x < Region_low_bound[e][i])
509 {
510 Region_low_bound[e][i] = x;
511 }
512 if (x > Region_upp_bound[e][i])
513 {
514 Region_upp_bound[e][i] = x;
515 }
516 }
517 }
518 }
519 }
520
521
522 /// Constructor. Provide vectors to "lower left" and "upper right"
523 /// vertices of refinement region
524 /// Also specify the node number of a central node
525 /// within elements -- it's used to determine if an element is
526 /// in the region where refinement is supposed to take place.
527 /// Optional boolean flag (defaulting to false) indicates that
528 /// refinement decision is based on Lagrangian coordinates -- only
529 /// applicable to solid meshes.
531 const Vector<double>& lower_left,
532 const Vector<double>& upper_right,
533 const unsigned& central_node_number,
534 const bool& use_lagrangian_coordinates = false)
535 : Use_lagrangian_coordinates(use_lagrangian_coordinates),
536 Central_node_number(central_node_number)
537 {
538#ifdef PARANOID
539 if (mesh_pt->nelement() == 0)
540 {
541 throw OomphLibError(
542 "Can't build error estimator if there are no elements in mesh\n",
543 OOMPH_CURRENT_FUNCTION,
544 OOMPH_EXCEPTION_LOCATION);
545 }
546#endif
547
548 unsigned nregion = 1;
549 Region_low_bound.resize(nregion);
550 Region_upp_bound.resize(nregion);
551 Region_low_bound[0] = lower_left;
552 Region_upp_bound[0] = upper_right;
553 }
554
555 /// Broken copy constructor
557
558 /// Broken assignment operator
559 void operator=(const DummyErrorEstimator&) = delete;
560
561 /// Empty virtual destructor
563
564
565 /// Compute the elemental error measures for a given mesh
566 /// and store them in a vector. Doc errors etc.
567 virtual void get_element_errors(Mesh*& mesh_pt,
568 Vector<double>& elemental_error,
569 DocInfo& doc_info)
570 {
571#ifdef PARANOID
572 if (doc_info.is_doc_enabled())
573 {
574 std::ostringstream warning_stream;
575 warning_stream
576 << "No output defined in DummyErrorEstimator::get_element_errors()\n"
577 << "Ignoring doc_info flag.\n";
578 OomphLibWarning(warning_stream.str(),
579 "DummyErrorEstimator::get_element_errors()",
580 OOMPH_EXCEPTION_LOCATION);
581 }
582#endif
583 unsigned nregion = Region_low_bound.size();
584 unsigned nelem = mesh_pt->nelement();
585 for (unsigned e = 0; e < nelem; e++)
586 {
587 elemental_error[e] = 0.0;
588
589 // Check if element is in the regions to be refined
590 // (based on coords of its central node)
591 Node* nod_pt =
593 for (unsigned r = 0; r < nregion; r++)
594 {
595 bool is_inside = true;
596 unsigned dim = Region_low_bound[r].size();
597 for (unsigned i = 0; i < dim; i++)
598 {
599 double x = nod_pt->x(i);
601 {
602 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
603 if (solid_nod_pt != 0)
604 {
605 x = solid_nod_pt->xi(i);
606 }
607 }
608 if (x < Region_low_bound[r][i])
609 {
610 is_inside = false;
611 break;
612 }
613 if (x > Region_upp_bound[r][i])
614 {
615 is_inside = false;
616 break;
617 }
618 }
619 if (is_inside)
620 {
621 elemental_error[e] = 1.0;
622 break;
623 }
624 }
625 }
626 }
627
628 private:
629 /// Use Lagrangian coordinates to decide which element is to be
630 /// refined
632
633 /// Number of local node that is used to identify if an element
634 /// is located in the refinement region
636
637 /// Upper bounds for the coordinates of the refinement regions
639
640 /// Lower bounds for the coordinates of the refinement regions
642 };
643
644
645} // namespace oomph
646
647#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 nvertex_node() const =0
Number of vertex nodes in the element.
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 Node * vertex_node_pt(const unsigned &j) const =0
Pointer to the j-th vertex node in the element. Needed for efficient patch assmbly.
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:1313
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
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,...
unsigned & recovery_order()
Access function for order of recovery polynomials.
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.
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.
double & reference_flux_norm()
Access function for prescribed reference flux norm.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...