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