refineable_quad_element.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_REFINEABLE_QUAD_ELEMENT_HEADER
27#define OOMPH_REFINEABLE_QUAD_ELEMENT_HEADER
28
29// Config header generated by autoconfig
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34
35// oomph-lib headers
36#include "quadtree.h"
37#include "macro_element.h"
38#include "refineable_elements.h"
39#include "Qelements.h"
40
41namespace oomph
42{
43 // Forward definition for mesh.
44 class Mesh;
45
46 //=======================================================================
47 /// Refineable version of QElement<2,NNODE_1D>.
48 ///
49 /// Refinement is performed by quadtree procedures. When the element is
50 /// subdivided, the geometry of its sons is established by calls
51 /// to their father's
52 /// \code get_x(...) \endcode
53 /// function which refers to
54 /// - the father element's geometric FE mapping (this
55 /// is the default)
56 /// .
57 /// or
58 /// - to a MacroElement 's MacroElement::macro_map (if the pointer
59 /// to the macro element is non-NULL)
60 ///
61 /// The class provides a generic RefineableQElement<2>::build() function
62 /// which deals with generic
63 /// isoparametric QElements in which all values are associated with
64 /// nodes. The RefineableQElement<2>::further_build() function provides
65 /// an interface for any element-specific non-generic build operations.
66 ///
67 //=======================================================================
68 template<>
69 class RefineableQElement<2> : public virtual RefineableElement,
70 public virtual QuadElementBase
71 {
72 public:
73 /// Shorthand for pointer to an argument-free void member
74 /// function of the refineable element
75 typedef void (RefineableQElement<2>::*VoidMemberFctPt)();
76
77 /// Constructor: Pass refinement level (default 0 = root)
79 {
80#ifdef LEAK_CHECK
81 LeakCheckNames::RefineableQElement<2> _build += 1;
82#endif
83 }
84
85
86 /// Broken copy constructor
88
89 /// Broken assignment operator
90 // Commented out broken assignment operator because this can lead to a
91 // conflict warning when used in the virtual inheritence hierarchy.
92 // Essentially the compiler doesn't realise that two separate
93 // implementations of the broken function are the same and so, quite
94 // rightly, it shouts.
95 /*void operator=(const RefineableQElement<2>&) = delete;*/
96
97 /// Destructor
99 {
100#ifdef LEAK_CHECK
101 LeakCheckNames::RefineableQElement<2> _build -= 1;
102#endif
103 }
104
105 /// A refineable quad element has four sons
106 unsigned required_nsons() const
107 {
108 return 4;
109 }
110
111 /// If a neighbouring element has already created a node at
112 /// a position corresponding to the local fractional position within the
113 /// present element, s_fraction, return
114 /// a pointer to that node. If not, return NULL (0).
115 /// If the node is on a periodic boundary the flag is_periodic is true,
116 /// otherwise it will be false.
117 virtual Node* node_created_by_neighbour(const Vector<double>& s_fraction,
118 bool& is_periodic);
119
120 /// If a son of a neighbouring element has already created a node at
121 /// a position corresponding to the local fractional position within the
122 /// present element, s_fraction, return
123 /// a pointer to that node. If not, return NULL (0).
124 /// If the node is on a periodic boundary the flag is_periodic is true,
125 /// otherwise it will be false.
127 const Vector<double>& s_fraction, bool& is_periodic)
128 {
129 // It is impossible for this situation to arise in meshes
130 // containing elements of uniform p-order. This is here so
131 // that it can be overloaded for p-refineable elements.
132 return 0;
133 }
134
135 /// Build the element, i.e. give it nodal positions, apply BCs, etc.
136 /// Pointers to any new nodes will be returned in new_node_pt. If
137 /// it is open, the positions of the new
138 /// nodes will be written to the file stream new_nodes_file
139 virtual void build(Mesh*& mesh_pt,
140 Vector<Node*>& new_node_pt,
141 bool& was_already_built,
142 std::ofstream& new_nodes_file);
143
144 /// Check the integrity of the element: ensure that the position and
145 /// values are continuous across the element edges
146 void check_integrity(double& max_error);
147
148 /// Print corner nodes, use colour
149 void output_corners(std::ostream& outfile, const std::string& colour) const;
150
151 /// Pointer to quadtree representation of this element
153 {
154 return dynamic_cast<QuadTree*>(Tree_pt);
155 }
156
157 /// Pointer to quadtree representation of this element
159 {
160 return dynamic_cast<QuadTree*>(Tree_pt);
161 }
162
163 /// Markup all hanging nodes & document the results in
164 /// the output streams contained in the vector output_stream, if they
165 /// are open.
166 void setup_hanging_nodes(Vector<std::ofstream*>& output_stream);
167
168 /// Perform additional hanging node procedures for variables
169 /// that are not interpolated by all nodes (e.g. lower order interpolations
170 /// as for the pressure in Taylor Hood).
171 virtual void further_setup_hanging_nodes() = 0;
172
173 protected:
174 /// Coincidence between son nodal points and father boundaries:
175 /// Father_bound[node_1d](jnod_son,son_type)={SW/SE/NW/NE/S/E/N/W/OMEGA}
176 static std::map<unsigned, DenseMatrix<int>> Father_bound;
177
178 /// Setup static matrix for coincidence between son
179 /// nodal points and father boundaries
180 void setup_father_bounds();
181
182 /// Determine Vector of boundary conditions along edge (N/S/W/E)
183 void get_edge_bcs(const int& edge, Vector<int>& bound_cons) const;
184
185 public:
186 /// Determine set of (mesh) boundaries that the
187 /// element edge/vertex lives on
188 void get_boundaries(const int& edge, std::set<unsigned>& boundaries) const;
189
190 /// Determine Vector of boundary conditions along edge
191 /// (or on vertex) bound (S/W/N/E/SW/SE/NW/NE): For value ival
192 /// on this boundary, bound_cons[ival]=1 if pinned and 0 if free.
193 void get_bcs(int bound, Vector<int>& bound_cons) const;
194
195 /// Return the value of the intrinsic boundary coordinate
196 /// interpolated along the edge (S/W/N/E)
197 void interpolated_zeta_on_edge(const unsigned& boundary,
198 const int& edge,
199 const Vector<double>& s,
200 Vector<double>& zeta);
201
202 protected:
203 /// Internal helper function that is used to construct the
204 /// hanging node schemes for the value_id-th interpolated value
205 void setup_hang_for_value(const int& value_id);
206
207 /// Internal helper function that is used to construct the
208 /// hanging node schemes for the positions.
209 virtual void quad_hang_helper(const int& value_id,
210 const int& my_edge,
211 std::ofstream& output_hangfile);
212 };
213
214
215 //========================================================================
216 /// Refineable version of Solid quad elements
217 //========================================================================
218 template<>
219 class RefineableSolidQElement<2> : public virtual RefineableQElement<2>,
220 public virtual RefineableSolidElement,
221 public virtual QSolidElementBase
222 {
223 public:
224 /// Constructor, just call the constructor of the RefineableQElement<2>
227 {
228 }
229
230
231 /// Broken copy constructor
233
234 /// Broken assignment operator
235 /*void operator=(const RefineableSolidQElement<2>&) = delete;*/
236
237 /// Virtual Destructor
239
240
241 /// Final over-ride: Use version in QSolidElementBase
242 void set_macro_elem_pt(MacroElement* macro_elem_pt)
243 {
245 }
246
247 /// Final over-ride: Use version in QSolidElementBase
248 void set_macro_elem_pt(MacroElement* macro_elem_pt,
249 MacroElement* undeformed_macro_elem_pt)
250 {
252 undeformed_macro_elem_pt);
253 }
254
255 /// Use the generic finite difference routine defined in
256 /// RefineableSolidElement to calculate the Jacobian matrix
258 {
259 RefineableSolidElement::get_jacobian(residuals, jacobian);
260 }
261
262 /// Determine vector of solid (positional) boundary conditions
263 /// along edge (N/S/W/E) [Pressure does not have to be included
264 /// since it can't be subjected to bc at more than one node anyway]
265 void get_edge_solid_bcs(const int& edge,
266 Vector<int>& solid_bound_cons) const;
267
268 /// Determine vector of solid (positional) boundary conditions
269 /// along edge (or on vertex) bound (S/W/N/E/SW/SE/NW/NE): For direction i,
270 /// solid_bound_cons[i]=1 if displacement in this coordinate direction
271 /// is pinned and 0 if it's free.
272 void get_solid_bcs(int bound, Vector<int>& solid_bound_cons) const;
273
274
275 /// Build the element, i.e. give it nodal positions, apply BCs, etc.
276 /// Incl. documention into new_nodes_file
277 // NOTE: FOR SOME REASON THIS NEEDS TO LIVE IN *.H TO WORK ON INTEL
278 void build(Mesh*& mesh_pt,
279 Vector<Node*>& new_node_pt,
280 bool& was_already_built,
281 std::ofstream& new_nodes_file)
282 {
283 using namespace QuadTreeNames;
284
285 // Call the standard (non-elastic) build function
287 mesh_pt, new_node_pt, was_already_built, new_nodes_file);
288
289 // Are we done?
290 if (was_already_built) return;
291
292 // Now need to loop over the nodes again and set solid variables
293
294 // What type of son am I? Ask my quadtree representation...
295 int son_type = Tree_pt->son_type();
296
297 // Which element (!) is my father? (We must have a father
298 // since was_already_built is false...)
299 RefineableSolidQElement<2>* father_el_pt =
300 dynamic_cast<RefineableSolidQElement<2>*>(
301 Tree_pt->father_pt()->object_pt());
302
303
304#ifdef PARANOID
305 // Currently we can't handle the case of generalised coordinates
306 // since we haven't established how they should be interpolated
307 // Buffer this case:
308 if (static_cast<SolidNode*>(father_el_pt->node_pt(0))
309 ->nlagrangian_type() != 1)
310 {
311 throw OomphLibError(
312 "We can't handle generalised nodal positions (yet).\n",
313 OOMPH_CURRENT_FUNCTION,
314 OOMPH_EXCEPTION_LOCATION);
315 }
316#endif
317
318 Vector<double> s_lo(2);
319 Vector<double> s_hi(2);
320 Vector<double> s(2);
321 Vector<double> xi(2);
322 Vector<double> xi_fe(2);
323 Vector<double> x(2);
324 Vector<double> x_fe(2);
325
326 // Setup vertex coordinates in father element:
327 //--------------------------------------------
328 switch (son_type)
329 {
330 case SW:
331 s_lo[0] = -1.0;
332 s_hi[0] = 0.0;
333 s_lo[1] = -1.0;
334 s_hi[1] = 0.0;
335 break;
336
337 case SE:
338 s_lo[0] = 0.0;
339 s_hi[0] = 1.0;
340 s_lo[1] = -1.0;
341 s_hi[1] = 0.0;
342 break;
343
344 case NE:
345 s_lo[0] = 0.0;
346 s_hi[0] = 1.0;
347 s_lo[1] = 0.0;
348 s_hi[1] = 1.0;
349 break;
350
351 case NW:
352 s_lo[0] = -1.0;
353 s_hi[0] = 0.0;
354 s_lo[1] = 0.0;
355 s_hi[1] = 1.0;
356 break;
357 }
358
359 // Pass the undeformed macro element onto the son
360 // hierher why can I read this?
361 if (father_el_pt->Undeformed_macro_elem_pt != 0)
362 {
363 Undeformed_macro_elem_pt = father_el_pt->Undeformed_macro_elem_pt;
364 for (unsigned i = 0; i < 2; i++)
365 {
366 s_macro_ll(i) =
367 father_el_pt->s_macro_ll(i) +
368 0.5 * (s_lo[i] + 1.0) *
369 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
370 s_macro_ur(i) =
371 father_el_pt->s_macro_ll(i) +
372 0.5 * (s_hi[i] + 1.0) *
373 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
374 }
375 }
376
377 // Local node number
378 unsigned n = 0;
379
380 // Find number of 1D nodes in element
381 unsigned n_p = nnode_1d();
382
383 // Loop over nodes in element
384 for (unsigned i0 = 0; i0 < n_p; i0++)
385 {
386 // Local coordinate in father element
387 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * double(i0) / double(n_p - 1);
388
389 for (unsigned i1 = 0; i1 < n_p; i1++)
390 {
391 // Local coordinate in father element
392 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * double(i1) / double(n_p - 1);
393
394 // Local node number
395 n = i0 + n_p * i1;
396
397 // Get position from father element -- this uses the macro
398 // element representation(s) if appropriate. If the node
399 // turns out to be a hanging node later on, then
400 // its position gets adjusted in line with its
401 // hanging node interpolation.
402 father_el_pt->get_x_and_xi(s, x_fe, x, xi_fe, xi);
403
404 // Cast the node to an Solid node
405 SolidNode* elastic_node_pt = static_cast<SolidNode*>(node_pt(n));
406
407 for (unsigned i = 0; i < 2; i++)
408 {
409 // x_fe is the FE representation -- this is all we can
410 // work with in a solid mechanics problem. If you wish
411 // to reposition nodes on curvilinear boundaries of
412 // a domain to their exact positions on those boundaries
413 // you'll have to do this yourself! [Note: We used to
414 // use the macro-element-based representation
415 // to assign the position of pinned nodes but this is not always
416 // correct since pinning doesn't mean "pin in place" or
417 // "pin to the curvilinear boundary". For instance, we could impose
418 // the boundary displacement manually.
419 // x_fe is the FE representation
420 elastic_node_pt->x(i) = x_fe[i];
421
422 // Lagrangian coordinates can come from undeformed macro element
423
424 if (Use_undeformed_macro_element_for_new_lagrangian_coords)
425 {
426 elastic_node_pt->xi(i) = xi[i];
427 }
428 else
429 {
430 elastic_node_pt->xi(i) = xi_fe[i];
431 }
432 }
433
434
435 // Are there any history values to be dealt with?
436 TimeStepper* time_stepper_pt =
437 father_el_pt->node_pt(0)->time_stepper_pt();
438
439 // Number of history values (incl. present)
440 unsigned ntstorage = time_stepper_pt->ntstorage();
441 if (ntstorage != 1)
442 {
443 // Loop over # of history values (excluding present which has been
444 // done above)
445 for (unsigned t = 1; t < ntstorage; t++)
446 {
447 // History values can (and in the case of Newmark timestepping,
448 // the scheme most likely to be used for Solid computations, do)
449 // include non-positional values, e.g. velocities and
450 // accelerations.
451
452 // Set previous positions of the new node
453 for (unsigned i = 0; i < 2; i++)
454 {
455 elastic_node_pt->x(t, i) =
456 father_el_pt->interpolated_x(t, s, i);
457 }
458 }
459 }
460
461 } // End of vertical loop over nodes in element
462
463 } // End of horizontal loop over nodes in element
464 }
465 };
466
467} // namespace oomph
468
469#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the elemental Jacobian matrix "d equation / d variable".
Definition: elements.h:990
Base class for MacroElement s that are used during mesh refinement in domains with curvlinear and/or ...
Definition: macro_element.h:73
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
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....
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:202
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:186
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:331
virtual void set_macro_elem_pt(MacroElement *macro_elem_pt)
Broken assignment operator.
Definition: Qelements.h:349
void get_x_and_xi(const Vector< double > &s, Vector< double > &x_fe, Vector< double > &x, Vector< double > &xi_fe, Vector< double > &xi) const
Eulerian and Lagrangian coordinates as function of the local coordinates: The Eulerian position is re...
Definition: Qelements.h:379
Base class for all quad elements.
Definition: Qelements.h:821
QuadTree class: Recursively defined, generalised quadtree.
Definition: quadtree.h:104
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
Refineable version of QElement<2,NNODE_1D>.
virtual void further_setup_hanging_nodes()=0
Perform additional hanging node procedures for variables that are not interpolated by all nodes (e....
QuadTree * quadtree_pt() const
Pointer to quadtree representation of this element.
RefineableQElement()
Constructor: Pass refinement level (default 0 = root)
static std::map< unsigned, DenseMatrix< int > > Father_bound
Coincidence between son nodal points and father boundaries: Father_bound[node_1d](jnod_son,...
unsigned required_nsons() const
A refineable quad element has four sons.
virtual Node * node_created_by_son_of_neighbour(const Vector< double > &s_fraction, bool &is_periodic)
If a son of a neighbouring element has already created a node at a position corresponding to the loca...
virtual ~RefineableQElement()
Broken assignment operator.
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
RefineableQElement(const RefineableQElement< 2 > &dummy)=delete
Broken copy constructor.
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 Solid quad elements.
void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)
Build the element, i.e. give it nodal positions, apply BCs, etc. Incl. documention into new_nodes_fil...
virtual ~RefineableSolidQElement()
Broken assignment operator.
void set_macro_elem_pt(MacroElement *macro_elem_pt, MacroElement *undeformed_macro_elem_pt)
Final over-ride: Use version in QSolidElementBase.
RefineableSolidQElement()
Constructor, just call the constructor of the RefineableQElement<2>
void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Use the generic finite difference routine defined in RefineableSolidElement to calculate the Jacobian...
void set_macro_elem_pt(MacroElement *macro_elem_pt)
Final over-ride: Use version in QSolidElementBase.
RefineableSolidQElement(const RefineableSolidQElement< 2 > &dummy)=delete
Broken copy constructor.
A class that is used to template the solid refineable Q elements by dimension. It's really nothing mo...
Definition: Qelements.h:2286
MacroElement * Undeformed_macro_elem_pt
Pointer to the element's "undeformed" macro element (NULL by default)
Definition: elements.h:4076
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
unsigned nlagrangian_type() const
Number of types of Lagrangian coordinates used to interpolate the Lagrangian coordinates within the e...
Definition: nodes.h:1877
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1883
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
long RefineableQElement< 2 > _build
Vector< std::string > colour
Tecplot colours.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...