refineable_poisson_elements.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// Header file for refineable QPoissonElement elements
27
28#ifndef OOMPH_REFINEABLE_POISSON_ELEMENTS_HEADER
29#define OOMPH_REFINEABLE_POISSON_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36
37// oomph-lib headers
38#include "../generic/refineable_quad_element.h"
39#include "../generic/refineable_brick_element.h"
40#include "../generic/hp_refineable_elements.h"
41#include "../generic/error_estimator.h"
42#include "poisson_elements.h"
43
44namespace oomph
45{
46 /// ////////////////////////////////////////////////////////////////////////
47 /// ////////////////////////////////////////////////////////////////////////
48
49
50 //======================================================================
51 /// Refineable version of Poisson equations
52 ///
53 ///
54 //======================================================================
55 template<unsigned DIM>
56 class RefineablePoissonEquations : public virtual PoissonEquations<DIM>,
57 public virtual RefineableElement,
58 public virtual ElementWithZ2ErrorEstimator
59 {
60 public:
61 /// Constructor, simply call other constructors
63 : PoissonEquations<DIM>(),
66 {
67 }
68
69 /// Broken copy constructor
71 delete;
72
73 /// Broken assignment operator
75
76 /// Number of 'flux' terms for Z2 error estimation
78 {
79 return DIM;
80 }
81
82 /// Get 'flux' for Z2 error recovery: Standard flux.from Poisson equations
84 {
85 this->get_flux(s, flux);
86 }
87
88 /// Get error against and norm of exact flux
90 std::ostream& outfile,
92 double& error,
93 double& norm);
94
95 /// Get the function value u in Vector.
96 /// Note: Given the generality of the interface (this function
97 /// is usually called from black-box documentation or interpolation
98 /// routines), the values Vector sets its own size in here.
100 Vector<double>& values)
101 {
102 // Set size of Vector: u
103 values.resize(1);
104
105 // Find number of nodes
106 unsigned n_node = nnode();
107
108 // Local shape function
109 Shape psi(n_node);
110
111 // Find values of shape function
112 shape(s, psi);
113
114 // Initialise value of u
115 values[0] = 0.0;
116
117 // Find the index at which the poisson unknown is stored
118 unsigned u_nodal_index = this->u_index_poisson();
119
120 // Loop over the local nodes and sum up the values
121 for (unsigned l = 0; l < n_node; l++)
122 {
123 values[0] += this->nodal_value(l, u_nodal_index) * psi[l];
124 }
125 }
126
127
128 /// Get the function value u in Vector.
129 /// Note: Given the generality of the interface (this function
130 /// is usually called from black-box documentation or interpolation
131 /// routines), the values Vector sets its own size in here.
132 void get_interpolated_values(const unsigned& t,
133 const Vector<double>& s,
134 Vector<double>& values)
135 {
136 if (t != 0)
137 {
138 std::string error_message =
139 "Time-dependent version of get_interpolated_values() ";
140 error_message += "not implemented for this element \n";
141 throw OomphLibError(
142 error_message,
143 "RefineablePoissonEquations::get_interpolated_values()",
144 OOMPH_EXCEPTION_LOCATION);
145 }
146 else
147 {
148 // Make sure that we call this particular object's steady
149 // get_interpolated_values (it could get overloaded lower down)
151 }
152 }
153
154
155 /// Further build: Copy source function pointer from father element
157 {
158 this->Source_fct_pt = dynamic_cast<RefineablePoissonEquations<DIM>*>(
159 this->father_element_pt())
160 ->source_fct_pt();
161 }
162
163
164 private:
165 /// Add element's contribution to elemental residual vector and/or
166 /// Jacobian matrix
167 /// flag=1: compute both
168 /// flag=0: compute only residual vector
170 Vector<double>& residuals,
171 DenseMatrix<double>& jacobian,
172 const unsigned& flag);
173
174 /// Compute derivatives of elemental residual vector with respect
175 /// to nodal coordinates. Overwrites default implementation in
176 /// FiniteElement base class.
177 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
179 RankThreeTensor<double>& dresidual_dnodal_coordinates);
180 };
181
182
183 //======================================================================
184 /// Refineable version of 2D QPoissonElement elements
185 ///
186 ///
187 //======================================================================
188 template<unsigned DIM, unsigned NNODE_1D>
190 : public QPoissonElement<DIM, NNODE_1D>,
191 public virtual RefineablePoissonEquations<DIM>,
192 public virtual RefineableQElement<DIM>
193 {
194 public:
195 /// Constructor, simply call the other constructors
199 RefineableQElement<DIM>(),
200 QPoissonElement<DIM, NNODE_1D>()
201 {
202 }
203
204
205 /// Broken copy constructor
207 const RefineableQPoissonElement<DIM, NNODE_1D>& dummy) = delete;
208
209 /// Broken assignment operator
211
212 /// Number of continuously interpolated values: 1
214 {
215 return 1;
216 }
217
218 /// Number of vertex nodes in the element
219 unsigned nvertex_node() const
220 {
222 }
223
224 /// Pointer to the j-th vertex node in the element
225 Node* vertex_node_pt(const unsigned& j) const
226 {
228 }
229
230 /// Rebuild from sons: empty
231 void rebuild_from_sons(Mesh*& mesh_pt) {}
232
233 /// Order of recovery shape functions for Z2 error estimation:
234 /// Same order as shape functions.
236 {
237 return (NNODE_1D - 1);
238 }
239
240 /// Perform additional hanging node procedures for variables
241 /// that are not interpolated by all nodes. Empty.
243 };
244
245
246 //======================================================================
247 /// p-refineable version of 2D QPoissonElement elements
248 //======================================================================
249 template<unsigned DIM>
251 : public QPoissonElement<DIM, 2>,
252 public virtual RefineablePoissonEquations<DIM>,
253 public virtual PRefineableQElement<DIM>
254 {
255 public:
256 /// Constructor, simply call the other constructors
260 PRefineableQElement<DIM>(),
261 QPoissonElement<DIM, 2>()
262 {
263 // Set integration scheme
264 // (To avoid memory leaks in pre-build and p-refine where new
265 // integration schemes are created)
267 }
268
269 /// Destructor (to avoid memory leaks)
271 {
272 delete this->integral_pt();
273 }
274
275
276 /// Broken copy constructor
278 delete;
279
280 /// Broken assignment operator
282
283 void further_build();
284
285 /// Number of continuously interpolated values: 1
287 {
288 return 1;
289 }
290
291 /// Number of vertex nodes in the element
292 unsigned nvertex_node() const
293 {
295 }
296
297 /// Pointer to the j-th vertex node in the element
298 Node* vertex_node_pt(const unsigned& j) const
299 {
301 }
302
303 /// Order of recovery shape functions for Z2 error estimation:
304 /// - Same order as shape functions.
305 // unsigned nrecovery_order()
306 // {
307 // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
308 // else {return 3;}
309 // }
310 /// - Constant recovery order, since recovery order of the first element
311 /// is used for the whole mesh.
313 {
314 return 3;
315 }
316
318 std::ostream& outfile,
320 double& error,
321 double& norm);
322 };
323
324
325 /// /////////////////////////////////////////////////////////////////////
326 /// /////////////////////////////////////////////////////////////////////
327 /// /////////////////////////////////////////////////////////////////////
328
329
330 //=======================================================================
331 /// Face geometry for the RefineableQuadPoissonElement elements: The spatial
332 /// dimension of the face elements is one lower than that of the
333 /// bulk element but they have the same number of points
334 /// along their 1D edges.
335 //=======================================================================
336 template<unsigned DIM, unsigned NNODE_1D>
338 : public virtual QElement<DIM - 1, NNODE_1D>
339 {
340 public:
341 /// Constructor: Call the constructor for the
342 /// appropriate lower-dimensional QElement
343 FaceGeometry() : QElement<DIM - 1, NNODE_1D>() {}
344 };
345
346} // namespace oomph
347
348#endif
static char t char * s
Definition: cfortran.h:568
char t
Definition: cfortran.h:568
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
FaceGeometry()
Constructor: Call the constructor for the appropriate lower-dimensional QElement.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
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
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
p-refineable version of 2D QPoissonElement elements
void operator=(const PRefineableQPoissonElement< DIM > &)=delete
Broken assignment operator.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
void further_build()
Further build: Copy source function pointer from father element.
unsigned nvertex_node() const
Number of vertex nodes in the element.
PRefineableQPoissonElement(const PRefineableQPoissonElement< DIM > &dummy)=delete
Broken copy constructor.
~PRefineableQPoissonElement()
Destructor (to avoid memory leaks)
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
PRefineableQPoissonElement()
Constructor, simply call the other constructors.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
A class for all isoparametric elements that solve the Poisson equations.
PoissonSourceFctPt & source_fct_pt()
Access function: Pointer to source function.
PoissonSourceFctPt Source_fct_pt
Pointer to source function:
void get_flux(const Vector< double > &s, Vector< double > &flux) const
Get flux: flux[i] = du/dx_i.
virtual unsigned u_index_poisson() const
Return the index at which the unknown value is stored. The default value, 0, is appropriate for singl...
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
////////////////////////////////////////////////////////////////////////
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void fill_in_generic_residual_contribution_poisson(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
void operator=(const RefineablePoissonEquations< DIM > &)=delete
Broken assignment operator.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Standard flux.from Poisson equations.
void compute_exact_Z2_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_flux_pt, double &error, double &norm)
Get error against and norm of exact flux.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
RefineablePoissonEquations()
Constructor, simply call other constructors.
RefineablePoissonEquations(const RefineablePoissonEquations< DIM > &dummy)=delete
Broken copy constructor.
void further_build()
Further build: Copy source function pointer from father element.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of 2D QPoissonElement elements.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 1.
RefineableQPoissonElement()
Constructor, simply call the other constructors.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void operator=(const RefineableQPoissonElement< DIM, NNODE_1D > &)=delete
Broken assignment operator.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
RefineableQPoissonElement(const RefineableQPoissonElement< DIM, NNODE_1D > &dummy)=delete
Broken copy constructor.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...