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