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