refineable_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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Header file for classes that define refineable element objects
31 
32 //Include guard to prevent multiple inclusions of the header
33 #ifndef OOMPH_REFINEABLE_ELEMENTS_HEADER
34 #define OOMPH_REFINEABLE_ELEMENTS_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 #include "elements.h"
42 #include "tree.h"
43 
44 namespace oomph
45 {
46 
47 class Mesh;
48 
49 //=======================================================================
50 /// RefineableElements are FiniteElements that may be subdivided into
51 /// children to provide a better local approximation to the solution.
52 /// After non-uniform refinement adjacent elements need not necessarily have
53 /// nodes in common. A node that does not have a counterpart in its
54 /// neighbouring element is known as a hanging node and its position and
55 /// any data that it stores must be constrained to ensure inter-element
56 /// continuity.
57 ///
58 /// Generic data and function interfaces associated with refinement
59 /// are defined in this class.
60 ///
61 /// Additional data includes:
62 /// - a pointer to a general Tree object that is used to track the
63 /// refinement history,
64 /// - a refinement level (not necessarily the same as the level in the tree!),
65 /// - a flag indicating whether the element should be refined,
66 /// - a flag indicating whether the element should be de-refined,
67 /// - a global element number for plotting/validation purposes,
68 /// - storage for local equation numbers associated with hanging nodes.
69 ///
70 /// Additional functions perform the following generic tasks:
71 /// - provide access to additional data,
72 /// - setup local equation numbering for data associated with hanging nodes,
73 /// - generic finite-difference calculation of contributions to the
74 /// elemental jacobian from nodal data to include hanging nodes,
75 /// - split of the element into its sons,
76 /// - select and deselect the element for refinement,
77 /// - select and deselect the sons of the element for de-refinement (merging),
78 ///
79 /// In addition, there are a number of interfaces that specify element-specific
80 /// tasks. These should be overloaded in RefineableElements of particular
81 /// geometric types and perform the following tasks:
82 /// - return a pointer to the root and father elements in the tree structure,
83 /// - define the number of sons into which the element is divided,
84 /// - build the element: construct nodes, assign their positions,
85 /// values and any boundary conditions,
86 /// - recreate the element from its sons if they are merged,
87 /// - deactivate the element, perform any operations that are
88 /// required when the element is still in the tree, but no longer active
89 /// - set the number and provide access to the values interpolated
90 /// by the nodes,
91 /// - setup the hanging nodes
92 ///
93 /// In mixed element different sets of nodes are used to interpolate different
94 /// unknowns. Interfaces are provided for functions that can be used to find
95 /// the position of the nodes that interpolate the different unknowns. These
96 /// functions are used to setup hanging node information automatically in
97 /// particular elements, e.g. Taylor Hood Navier--Stokes.
98 /// The default implementation assumes that the elements are isoparameteric.
99 ///
100 //======================================================================
101 class RefineableElement : public virtual FiniteElement
102 {
103  protected:
104 
105  /// A pointer to a general tree object
107 
108  /// Refinement level
109  unsigned Refine_level;
110 
111  /// Flag for refinement
113 
114  /// Flag to indicate suppression of any refinement
116 
117  /// Flag for unrefinement
119 
120  /// Global element number -- for plotting/validation purposes
121  long Number;
122 
123  /// \short Max. allowed discrepancy in element integrity check
124  static double Max_integrity_tolerance;
125 
126  /// \short Static helper function that is used to check that the value_id
127  /// is in range
128  static void check_value_id(const int &n_continuously_interpolated_values,
129  const int &value_id);
130 
131 
132  /// \short Assemble the jacobian matrix for the mapping from local
133  /// to Eulerian coordinates, given the derivatives of the shape function
134  /// w.r.t the local coordinates.
135  /// Overload the standard version to use the hanging information for
136  /// the Eulerian coordinates.
138  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
139 
140  /// \short Assemble the the "jacobian" matrix of second derivatives of the
141  /// mapping from local to Eulerian coordinates, given
142  /// the second derivatives of the shape functions w.r.t. local coordinates.
143  /// Overload the standard version to use the hanging information for
144  /// the Eulerian coordinates.
146  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
147 
148  /// \short Assemble the covariant Eulerian base vectors, assuming that
149  /// the derivatives of the shape functions with respect to the local
150  /// coordinates have already been constructed.
151  /// Overload the standard version to account for hanging nodes.
153  const DShape &dpsids, DenseMatrix<double> &interpolated_G) const;
154 
155  /// \short Calculate the mapping from local to Eulerian coordinates given
156  /// the derivatives of the shape functions w.r.t the local coordinates.
157  /// assuming that the coordinates are aligned in the direction of the local
158  /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
159  /// This funciton returns the determinant of the jacobian, the jacobian
160  /// and the inverse jacobian. Overload the standard version to take
161  /// hanging info into account.
163  const DShape &dpsids,DenseMatrix<double> &jacobian,
164  DenseMatrix<double> &inverse_jacobian) const;
165 
166  private:
167 
168  /// \short Storage for local equation numbers of hanging node variables
169  /// (values stored at master nodes). It is
170  /// essential that these are indexed by a Node pointer because the Node
171  /// may be internal or external to the element.
172  /// local equation number = Local_hang_eqn(master_node_pt,ival)
173  std::map<Node*,int> *Local_hang_eqn;
174 
175  /// \short Lookup scheme for unique number associated with any of the nodes
176  /// that actively control the shape of the element (i.e. they are either
177  /// non-hanging nodes of this element or master nodes of hanging nodes.
178  std::map<Node*,unsigned> Shape_controlling_node_lookup;
179 
180  protected:
181 
182 
183  /// \short Assign the local equation numbers for hanging node variables
184  void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt);
185 
186  /// \short Calculate the contributions to the jacobian from the nodal
187  /// degrees of freedom using finite differences.
188  /// This version is overloaded to take hanging node information into
189  /// account
190  virtual void fill_in_jacobian_from_nodal_by_fd(Vector<double> &residuals,
191  DenseMatrix<double> &jacobian);
192 
193  public:
194 
195  /// \short Constructor, calls the FiniteElement constructor and initialises
196  /// the member data
198  To_be_refined(false), Refinement_is_enabled(true),
199  Sons_to_be_unrefined(false), Number(-1),
200  Local_hang_eqn(0) {}
201 
202  /// Destructor, delete the allocated storage for the hanging equations
203  // (The body is now in the cc file to keep the xlC compiler happy under AIX)
204  virtual ~RefineableElement();
205 
206  /// Broken copy constructor
208  {
209  BrokenCopy::broken_copy("RefineableElement");
210  }
211 
212  /// Broken assignment operator
214  {
215  BrokenCopy::broken_assign("RefineableElement");
216  }
217 
218  /// Access function: Pointer to quadtree representation of this element
219  Tree* tree_pt() {return Tree_pt;}
220 
221  /// Set pointer to quadtree representation of this element
222  void set_tree_pt(Tree* my_tree_pt) {Tree_pt=my_tree_pt;}
223 
224  /// \short Set the number of sons that can be constructed by the element
225  /// The default is none
226  virtual unsigned required_nsons() const {return 0;}
227 
228 
229  /// Flag to indicate suppression of any refinement
231 
232  /// Suppress of any refinement for this element
234 
235  /// Emnable refinement for this element
237 
238  /// \short Split the element into the number of sons to be
239  /// constructed and return a
240  /// vector of pointers to the sons. Elements are allocated, but they are
241  /// not given any properties. The refinement level of the sons is one
242  /// higher than that of the father elemern.
243  template<class ELEMENT>
244  void split(Vector<ELEMENT*>& son_pt) const
245  {
246  // Increase refinement level
247  int son_refine_level=Refine_level+1;
248 
249  //How many sons are to be constructed
250  unsigned n_sons = required_nsons();
251  //Resize the son pointer
252  son_pt.resize(n_sons);
253 
254  //Loop over the sons and construct
255  for(unsigned i=0;i<n_sons;i++)
256  {
257  son_pt[i]=new ELEMENT;
258  //Set the refinement level of the newly constructed son.
259  son_pt[i]->set_refinement_level(son_refine_level);
260  }
261  }
262 
263 
264  /// \short Access function that returns the local equation number for the
265  /// hanging node variables (values stored at master nodes). The local
266  /// equation number corresponds to the i-th unknown stored at the node
267  /// addressed by node_pt
268  inline int local_hang_eqn(Node* const &node_pt, const unsigned &i)
269  {
270 #ifdef RANGE_CHECKING
271  if(i > ncont_interpolated_values())
272  {
273  std::ostringstream error_message;
274  error_message << "Range Error: Value " << i
275  << " is not in the range (0,"
276  << ncont_interpolated_values()-1 << ")";
277  throw OomphLibError(error_message.str(),
278  OOMPH_CURRENT_FUNCTION,
279  OOMPH_EXCEPTION_LOCATION);
280  }
281 #endif
282 
283  return Local_hang_eqn[i][node_pt];
284  }
285 
286  /// \short Interface to function that builds the element: i.e. construct
287  /// the nodes, assign their positions, apply boundary conditions, etc. The
288  /// required procedures depend on the geometrical type of the element and
289  /// must be implemented in specific refineable elements. Any new nodes
290  /// created during the build process are returned in the vector
291  /// new_node_pt.
292  virtual void build(Mesh* &mesh_pt, Vector<Node*> &new_node_pt,
293  bool &was_already_built, std::ofstream &new_nodes_file)=0;
294 
295  /// Set the refinement level
296  void set_refinement_level(const int& refine_level)
297  {Refine_level=refine_level;}
298 
299  /// Return the Refinement level
300  unsigned refinement_level() const {return Refine_level;}
301 
302  /// Select the element for refinement
304 
305  /// Deselect the element for refinement
307 
308  /// Unrefinement will be performed by merging the four sons of this element
310 
311  /// No unrefinement will be performed by merging the four sons of this element
313 
314  /// Has the element been selected for refinement?
315  bool to_be_refined() {return To_be_refined;}
316 
317  /// Has the element been selected for unrefinement?
319 
320  /// \short Rebuild the element, e.g. set internal values in line with
321  /// those of the sons that have now merged
322  virtual void rebuild_from_sons(Mesh* &mesh_pt)=0;
323 
324  /// \short Unbuild the element, i.e. mark the nodes that were created
325  /// during its creation for possible deletion
326  virtual void unbuild()
327  {
328  //Get pointer to father element
329  RefineableElement* father_pt = father_element_pt();
330  //If there is no father, nothing to do
331  if(father_pt==0) {return;}
332 
333  //Loop over all the nodes
334  unsigned n_node = this->nnode();
335  for(unsigned n=0;n<n_node;n++)
336  {
337  //If any node in this element is in the father, it can't be deleted
338  if(father_pt->get_node_number(this->node_pt(n)) >= 0)
339  {
341  }
342  }
343  }
344 
345  /// \short Final operations that must be performed when the element is no
346  /// longer active in the mesh, but still resident in the QuadTree.
347  virtual void deactivate_element();
348 
349  /// Return true if all the nodes have been built, false if not
350  virtual bool nodes_built() {return node_pt(0)!=0;}
351 
352  /// Element number (for debugging/plotting)
353  long number() const {return Number;}
354 
355  /// Set element number (for debugging/plotting)
356  void set_number(const long& mynumber) {Number=mynumber;}
357 
358  /// \short Number of continuously interpolated values. Note: We assume
359  /// that they are located at the beginning of the value_pt Vector!
360  /// (Used for interpolation to son elements, for integrity check
361  /// and post-processing -- we can only expect
362  /// the continously interpolated values to be continous across
363  /// element boundaries).
364  virtual unsigned ncont_interpolated_values() const=0;
365 
366  /// \short Get all continously interpolated function values in this
367  /// element as a Vector. Note: Vector sets is own size to ensure that
368  /// that this function can be used in black-box fashion.
370  Vector<double>& values)
371  {
372  get_interpolated_values(0, s, values);
373  }
374 
375  /// \short Get all continously interpolated function values at previous
376  /// timestep in this element as a Vector. (t=0: present; t>0:
377  /// prev. timestep) Note: Vector sets is own size to ensure that that this
378  /// function can be used in black-box fashion.
379  virtual void get_interpolated_values(const unsigned& t,
380  const Vector<double>&s,
381  Vector<double>& values)=0;
382 
383  /// \short In mixed elements, different sets of nodes are used to interpolate
384  /// different unknowns. This function returns the n-th node that interpolates
385  /// the value_id-th unknown. Default implementation is that all
386  /// variables use the positional nodes, i.e. isoparametric elements. Note
387  /// that any overloaded versions of this function MUST provide a set
388  /// of nodes for the position, which always has the value_id -1.
389  virtual Node* interpolating_node_pt(const unsigned &n,
390  const int &value_id)
391 
392  {return node_pt(n);}
393 
394  /// \short Return the local one dimensional fraction of the n1d-th node
395  /// in the direction of the local coordinate s[i] that is used to interpolate
396  /// the value_id-th continuously interpolated unknown. Default assumes
397  /// isoparametric interpolation for all unknowns
399  const unsigned &n1d, const unsigned &i, const int &value_id)
400  {return local_one_d_fraction_of_node(n1d,i);}
401 
402  /// \short Return a pointer to the node that interpolates the value-id-th
403  /// unknown at local coordinate s. If there is not a node at that position,
404  /// then return 0.
405  virtual Node*
407  const int &value_id)
408 
409  {return get_node_at_local_coordinate(s);}
410 
411 
412  /// \short Return the number of nodes that are used to interpolate the
413  /// value_id-th unknown. Default is to assume isoparametric elements.
414  virtual unsigned ninterpolating_node(const int &value_id) {return nnode();}
415 
416  /// \short Return the number of nodes in a one_d direction that are
417  /// used to interpolate the value_id-th unknown. Default is to assume
418  /// an isoparametric mapping.
419  virtual unsigned ninterpolating_node_1d(const int &value_id)
420  {return nnode_1d();}
421 
422  /// \short Return the basis functions that are used to interpolate
423  /// the value_id-th unknown. By default assume isoparameteric interpolation
424  virtual void interpolating_basis(const Vector<double> &s,
425  Shape &psi,
426  const int &value_id) const
427  {shape(s,psi);}
428 
429  /// \short Check the integrity of the element: Continuity of positions
430  /// values, etc. Essentially, check that the approximation of the functions
431  /// is consistent when viewed from both sides of the element boundaries
432  /// Must be overloaded for each different geometric element
433  virtual void check_integrity(double &max_error)=0;
434 
435  /// \short Max. allowed discrepancy in element integrity check
436  static double& max_integrity_tolerance()
437  { return Max_integrity_tolerance;}
438 
439  /// \short The purpose of this function is to identify all possible
440  /// Data that can affect the fields interpolated by the FiniteElement.
441  /// This must be overloaded to include data from any hanging nodes
442  /// correctly
444  std::set<std::pair<Data*,unsigned> > &paired_field_data);
445 
446 
447  /// \short Overload the function that assigns local equation numbers
448  /// for the Data stored at the nodes so that hanging data is taken
449  /// into account
450  inline void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
451  {
453  assign_hanging_local_eqn_numbers(store_local_dof_pt);
454  }
455 
456  /// \short Pointer to the root element in refinement hierarchy (must be
457  /// implemented in specific elements that do refinement via
458  /// tree-like refinement structure. Here we provide a default
459  /// implementation that is appropriate for cases where tree-like
460  /// refinement doesn't exist or if the element doesn't have
461  /// root in that tree (i.e. if it's a root itself): We return
462  /// "this".
464  {
465  //If there is no tree -- the element is its own root
466  if(Tree_pt==0) {return this;}
467  //Otherwise it's the tree's root object
468  else {return Tree_pt->root_pt()->object_pt();}
469  }
470 
471  /// Return a pointer to the father element.
473  {
474  //If we have no tree, we have no father
475  if(Tree_pt==0) {return 0;}
476  else
477  {
478  //Otherwise get the father of the tree
479  Tree* father_pt = Tree_pt->father_pt();
480  //If the tree has no father then return null, no father
481  if(father_pt==0) {return 0;}
482  else {return father_pt->object_pt();}
483  }
484  }
485 
486  /// Return a pointer to the "father" element at the specified refinement level
488  RefineableElement*& father_at_reflevel_pt)
489  {
490  // Get the father in the tree (it shouldn't try to get a null Tree...)
491  Tree* father_pt = Tree_pt->father_pt();
492  // Get the refineable element associated with this father
493  RefineableElement* father_el_pt=dynamic_cast<RefineableElement*>
494  (father_pt->object_pt());
495  // Get the refinement level
496  unsigned level=father_el_pt->refinement_level();
497  // If the level matches the required one then return, if not call again
498  if (level==refinement_level)
499  {
500  father_at_reflevel_pt=father_el_pt;
501  }
502  else
503  {
504  // Recursive call
505  father_el_pt->get_father_at_refinement_level(refinement_level,
506  father_at_reflevel_pt);
507  }
508  }
509 
510  /// \short Initial setup of the element: e.g. set the appropriate internal
511  /// p-order. If an adopted father is specified, information from this is
512  /// used instead of using the father found from the tree.
513  virtual void initial_setup(Tree* const &adopted_father_pt=0, const unsigned &initial_p_order=0) {}
514 
515  /// \short Pre-build the element
516  virtual void pre_build(Mesh*& mesh_pt, Vector<Node*>& new_node_pt) {}
517 
518  /// \short Further build: e.g. deal with interpolation of internal values
519  virtual void further_build() {}
520 
521  /// \short Mark up any hanging nodes that arise as a result of non-uniform
522  /// refinement. Any hanging nodes will be documented in files addressed by
523  /// the streams in the vector output_stream, if the streams are open.
524  virtual void setup_hanging_nodes(Vector<std::ofstream*> &output_stream) { }
525 
526  /// \short Perform additional hanging node procedures for variables
527  /// that are not interpolated by all nodes (e.g. lower order interpolations
528  /// for the pressure in Taylor Hood).
529  virtual void further_setup_hanging_nodes() { }
530 
531  /// \short Compute derivatives of elemental residual vector with respect
532  /// to nodal coordinates. Default implementation by FD can be overwritten
533  /// for specific elements.
534  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
535  /// This version is overloaded from the version in FiniteElement
536  /// and takes hanging nodes into account -- j in the above loop
537  /// loops over all the nodes that actively control the
538  /// shape of the element (i.e. they are non-hanging or master nodes of
539  /// hanging nodes in this element).
541  dresidual_dnodal_coordinates);
542 
543 
544  /// \short Number of shape-controlling nodes = the number
545  /// of non-hanging nodes plus the number of master nodes associated
546  /// with hanging nodes.
548  {
549  return Shape_controlling_node_lookup.size();
550  }
551 
552  /// \short Return lookup scheme for unique number associated
553  /// with any of the nodes that actively control the shape of the
554  /// element (i.e. they are either non-hanging nodes of this element
555  /// or master nodes of hanging nodes.
556  std::map<Node*,unsigned> shape_controlling_node_lookup()
557  {
559  }
560 
561 
562 };
563 
564 
565 //////////////////////////////////////////////////////////////////////
566 //////////////////////////////////////////////////////////////////////
567 //////////////////////////////////////////////////////////////////////
568 
569 
570 
571 
572 //======================================================================
573 /// p-refineable version of RefineableElement
574 //======================================================================
575 class PRefineableElement : public virtual RefineableElement
576 {
577 protected:
578 
579  /// The polynomial expansion order of the elemental basis functions
580  unsigned P_order;
581 
582  /// Flag for p-refinement
584 
585  /// Flag to indicate suppression of any refinement
587 
588  /// Flag for unrefinement
590 
591 public:
592 
593  /// \short Constructor, calls the RefineableElement constructor
596  To_be_p_unrefined(false) {}
597 
598  /// Destructor, empty
599  virtual ~PRefineableElement() {}
600 
601  /// Broken copy constructor
603  {
604  BrokenCopy::broken_copy("PRefineableElement");
605  }
606 
607  /// Broken assignment operator
609  {
610  BrokenCopy::broken_assign("PRefineableElement");
611  }
612 
613 
614  /// Flag to indicate suppression of any refinement
616 
617  /// Suppress of any refinement for this element
619 
620  /// Emnable refinement for this element
622 
623  /// Access function to P_order
624  unsigned &p_order() {return P_order;}
625 
626  /// Access function to P_order (const version)
627  unsigned p_order() const {return P_order;}
628 
629  /// Get the initial P_order
630  /// This is required so that elements which are constructed with a
631  /// higher p-order initially are not un-refined past this level (e.g.
632  /// in fluid problems where elements initially use quadratic velocity
633  /// and linear pressure)
634  /// Virtual because this needs to be set in templated derived classes
635  virtual unsigned initial_p_order() const=0;
636 
637  /// Select the element for p-refinement
639 
640  /// Deselect the element for p-refinement
642 
643  /// Select the element for p-unrefinement
645 
646  /// Deselect the element for p-unrefinement
648 
649  /// Has the element been selected for refinement?
651 
652  /// Has the element been selected for p-unrefinement?
654 
655  /// p-refine the element
656  virtual void p_refine(const int &inc,
657  Mesh* const &mesh_pt,
658  GeneralisedElement* const &clone_pt)=0;
659 
660  // Overload the nodes_built function to check every node
661  bool nodes_built()
662  {
663  // Must check that EVERY node exists
664  unsigned n_node = this->nnode();
665  for(unsigned n=0; n<n_node; n++)
666  {
667  if(this->node_pt(n)==0) {return false;}
668  }
669  // If we get here then all the nodes are built
670  return true;
671  }
672 
673 };
674 
675 
676 
677 //////////////////////////////////////////////////////////////////////
678 //////////////////////////////////////////////////////////////////////
679 //////////////////////////////////////////////////////////////////////
680 
681 
682 
683 
684 
685 //=======================================================================
686 /// A base class for elements that can have hanging nodes
687 /// but are not refineable as such. This class is usually used as a
688 /// base class for FaceElements that are attached to refineable
689 /// bulk elements (and stripped out before adapting the bulk
690 /// mesh, so they don't participate in the refimenent process
691 /// itself). We therefore simply break the pure virtual functions
692 /// that don't make any sense for such elements
693 //======================================================================
695  {
696 
697  public:
698 
699 
700  /// Broken build function -- shouldn't really be needed
701  void build(Mesh* &mesh_pt, Vector<Node*> &new_node_pt,
702  bool &was_already_built, std::ofstream &new_nodes_file)
703  {
704  std::ostringstream error_message;
705  error_message << "This function is broken as it's only needed/used \n"
706  << "during \"proper\" refinement\n";
707  throw OomphLibError(error_message.str(),
708  OOMPH_CURRENT_FUNCTION,
709  OOMPH_EXCEPTION_LOCATION);
710  }
711 
712 
713 
714  /// \short Broken function -- this shouldn't really be needed.
716  Vector<double>& values)
717  {
718  std::ostringstream error_message;
719  error_message << "This function is broken as it's only needed/used \n"
720  << "during \"proper\" refinement\n";
721  throw OomphLibError(
722  error_message.str(),
723  OOMPH_CURRENT_FUNCTION,
724  OOMPH_EXCEPTION_LOCATION);
725  }
726 
727  /// \short Broken function -- this shouldn't really be needed.
728  virtual void get_interpolated_values(const unsigned& t,
729  const Vector<double>&s,
730  Vector<double>& values)
731  {
732  std::ostringstream error_message;
733  error_message << "This function is broken as it's only needed/used \n"
734  << "during \"proper\" refinement\n";
735  throw OomphLibError(
736  error_message.str(),
737  OOMPH_CURRENT_FUNCTION,
738  OOMPH_EXCEPTION_LOCATION);
739  }
740 
741 
742  /// Broken function -- this shouldn't really be needed.
743  void check_integrity(double &max_error)
744  {
745  std::ostringstream error_message;
746  error_message << "This function is broken as it's only needed/used \n"
747  << "during \"proper\" refinement\n";
748  throw OomphLibError(
749  error_message.str(),
750  OOMPH_CURRENT_FUNCTION,
751  OOMPH_EXCEPTION_LOCATION);
752  }
753 
754  /// Broken function -- this shouldn't really be needed.
755  void rebuild_from_sons(Mesh* &mesh_pt)
756  {
757  std::ostringstream error_message;
758  error_message << "This function is broken as it's only needed/used \n"
759  << "during \"proper\" refinement\n";
760  throw OomphLibError(
761  error_message.str(),
762  OOMPH_CURRENT_FUNCTION,
763  OOMPH_EXCEPTION_LOCATION);
764  }
765 
766 
767 
768  };
769 
770 
771 //////////////////////////////////////////////////////////////////////
772 //////////////////////////////////////////////////////////////////////
773 //////////////////////////////////////////////////////////////////////
774 
775 
776 
777 
778 //=======================================================================
779 /// RefineableSolidElements are SolidFiniteElements that may
780 /// be subdivided into children to provide a better local approximation
781 /// to the solution. The distinction is required to keep a clean
782 /// separation between problems that alter nodal positions and others.
783 /// A number of procedures are generic and are included in this class.
784 //======================================================================
786  public virtual SolidFiniteElement
787 {
788 
789  private:
790 
791  /// \short Storage for local equation numbers of
792  /// hanging node variables associated with nodal positions.
793  /// local position equation number =
794  /// Local_position_hang_eqn(master_node_pt,ival)
795  std::map<Node*,DenseMatrix<int> > Local_position_hang_eqn;
796 
797 
798  /// \short Assign local equation numbers to the hanging values associated
799  /// with positions or additional solid values.
800  void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt);
801 
802 
803  protected:
804 
805  /// \short Flag deciding if the Lagrangian coordinates of newly-created interior
806  /// SolidNodes are to be determined by the father element's undeformed
807  /// MacroElement representation (if it has one). Default: False as it means that,
808  /// following a refinement an element is no longer in equilbrium (as the Lagrangian
809  /// coordinate is determined differently from the Eulerian one). However, basing
810  /// the Lagrangian coordinates on the undeformed MacroElement can be
811  /// more stable numerically and for steady problems it's not a big deal
812  /// either way as the difference between the two formulations only matters
813  /// at finite resolution so we have no right to say that one is "more correct"
814  /// than the other...
816 
817  /// \short Assemble the jacobian matrix for the mapping from local
818  /// to lagrangian coordinates, given the derivatives of the shape function
819  /// Overload the standard version to use the hanging information for
820  /// the lagrangian coordinates.
822  const DShape &dpsids, DenseMatrix<double> &jacobian) const;
823 
824  /// \short Assemble the the "jacobian" matrix of second derivatives, given
825  /// the second derivatives of the shape functions w.r.t. local coordinates
826  /// Overload the standard version to use the hanging information for
827  /// the lagrangian coordinates.
829  const DShape &d2psids, DenseMatrix<double> &jacobian2) const;
830 
831  /// \short Calculate the mapping from local to Lagrangian coordinates given
832  /// the derivatives of the shape functions w.r.t the local coorindates.
833  /// assuming that the coordinates are aligned in the direction of the local
834  /// coordinates, i.e. there are no cross terms and the jacobian is diagonal.
835  /// This function returns the determinant of the jacobian, the jacobian
836  /// and the inverse jacobian.
838  const DShape &dpsids,DenseMatrix<double> &jacobian,
839  DenseMatrix<double> &inverse_jacobian) const;
840 
841  public:
842 
843  /// Constructor
846 
847  /// Virtual Destructor, delete any allocated storage
849 
850  /// \short Overload the local equation numbers for Data stored as part
851  /// of solid nodes to include hanging node data
852  void assign_solid_local_eqn_numbers(const bool &store_local_dof_pt)
853  {
855  assign_solid_hanging_local_eqn_numbers(store_local_dof_pt);
856  }
857 
858  ///\short The number of geometric data affecting a
859  /// RefineableSolidFiniteElement is the positional Data of all
860  /// non-hanging nodes plus the geometric Data of all distinct
861  /// master nodes. Recomputed on the fly.
862  unsigned ngeom_data() const;
863 
864  /// \short Return pointer to the j-th Data item that the object's
865  /// shape depends on: Positional data of non-hanging nodes and
866  /// positional data of master nodes. Recomputed on the fly.
867  Data* geom_data_pt(const unsigned& j);
868 
869  /// \short Specify Data that affects the geometry of the element
870  /// by adding the position Data to the set that's passed in.
871  /// (This functionality is required in FSI problems; set is used to
872  /// avoid double counting). Refineable version includes hanging nodes
873  void identify_geometric_data(std::set<Data*> &geometric_data_pt);
874 
875  /// \short Compute element residual Vector and element Jacobian matrix
876  /// corresponding to the solid positions. Overloaded version to take
877  /// the hanging nodes into account
879  DenseMatrix<double> &jacobian);
880 
881  /// \short Return the flag deciding if the Lagrangian coordinates of
882  /// newly-created interior SolidNodes are to be determined by the father
883  /// element's undeformed MacroElement representation (if it has one).
884  /// Default: False as it means that, following a refinement an element
885  /// is no longer in equilbrium (as the Lagrangian coordinate is
886  /// determined differently from the Eulerian one). However, basing
887  /// the Lagrangian coordinates on the undeformed MacroElement can be
888  /// more stable numerically and for steady problems it's not a big deal
889  /// either way as the difference between the two formulations only matters
890  /// at finite resolution so we have no right to say that one is "more correct"
891  /// than the other...
894 
895  /// \short Set the flag deciding if the Lagrangian coordinates of
896  /// newly-created interior SolidNodes are to be determined by the father
897  /// element's undeformed MacroElement representation (if it has one).
900 
901  /// \short Unset the flag deciding if the Lagrangian coordinates of
902  /// newly-created interior SolidNodes are to be determined by the father
903  /// element's undeformed MacroElement representation (if it has one).
906 
907  /// \short Access the local equation number of of hanging node variables
908  /// associated with nodal positions. The function returns a dense
909  /// matrix that contains all the local equation numbers corresponding to
910  /// the positional degrees of freedom.
913 
914  /// \short Further build: Pass the father's
915  /// Use_undeformed_macro_element_for_new_lagrangian_coords
916  /// flag down, then call the underlying RefineableElement's
917  /// version.
918  virtual void further_build()
919  {
921  dynamic_cast<RefineableSolidElement*>(father_element_pt())
923 
925  }
926 
927 };
928 
929 
930 //////////////////////////////////////////////////////////////////////
931 //////////////////////////////////////////////////////////////////////
932 //////////////////////////////////////////////////////////////////////
933 
934 
935 
936 
937 
938 //=======================================================================
939 /// A base class for SolidElements that can have hanging nodes
940 /// but are not refineable as such. This class is usually used as a
941 /// base class for FaceElements that are attached to refineable
942 /// bulk elements (and stripped out before adapting the bulk
943 /// mesh, so they don't participate in the refimenent process
944 /// itself). We therefore simply break the pure virtual functions
945 /// that don't make any sense for such elements
946 //======================================================================
949 public virtual RefineableSolidElement
950 {
951 
952 };
953 
954 
955 //////////////////////////////////////////////////////////////////////
956 //////////////////////////////////////////////////////////////////////
957 //////////////////////////////////////////////////////////////////////
958 
959 
960 }
961 
962 #endif
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
A Generalised Element class.
Definition: elements.h:76
virtual unsigned initial_p_order() const =0
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on: Positional data of non-hangi...
void identify_field_data_for_interactions(std::set< std::pair< Data *, unsigned > > &paired_field_data)
The purpose of this function is to identify all possible.
void select_sons_for_unrefinement()
Unrefinement will be performed by merging the four sons of this element.
void assemble_eulerian_base_vectors(const DShape &dpsids, DenseMatrix< double > &interpolated_G) const
Assemble the covariant Eulerian base vectors, assuming that the derivatives of the shape functions wi...
bool to_be_p_unrefined()
Has the element been selected for p-unrefinement?
virtual Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns. This function returns the n-th node that interpolates the value_id-th unknown. Default implementation is that all variables use the positional nodes, i.e. isoparametric elements. Note that any overloaded versions of this function MUST provide a set of nodes for the position, which always has the value_id -1.
void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)
Broken build function – shouldn't really be needed.
void assign_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for hanging node variables.
void deselect_for_p_unrefinement()
Deselect the element for p-unrefinement.
void select_for_p_refinement()
Select the element for p-refinement.
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:102
PRefineableElement()
Constructor, calls the RefineableElement constructor.
void deselect_for_refinement()
Deselect the element for refinement.
void split(Vector< ELEMENT * > &son_pt) const
Split the element into the number of sons to be constructed and return a vector of pointers to the so...
virtual void initial_setup(Tree *const &adopted_father_pt=0, const unsigned &initial_p_order=0)
Initial setup of the element: e.g. set the appropriate internal p-order. If an adopted father is spec...
Tree * Tree_pt
A pointer to a general tree object.
unsigned & p_order()
Access function to P_order.
static double Max_integrity_tolerance
Max. allowed discrepancy in element integrity check.
cstr elem_len * i
Definition: cfortran.h:607
void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Default implementation by FD can be overwritten for specific elements. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij} This version is overloaded from the version in FiniteElement and takes hanging nodes into account – j in the above loop loops over all the nodes that actively control the shape of the element (i.e. they are non-hanging or master nodes of hanging nodes in this element).
bool P_refinement_is_enabled
Flag to indicate suppression of any refinement.
static void check_value_id(const int &n_continuously_interpolated_values, const int &value_id)
Static helper function that is used to check that the value_id is in range.
void operator=(const RefineableElement &)
Broken assignment operator.
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
bool to_be_p_refined()
Has the element been selected for refinement?
void check_integrity(double &max_error)
Broken function – this shouldn't really be needed.
bool Use_undeformed_macro_element_for_new_lagrangian_coords
Flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be determined...
A general Finite Element class.
Definition: elements.h:1271
char t
Definition: cfortran.h:572
virtual ~RefineableSolidElement()
Virtual Destructor, delete any allocated storage.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
virtual Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
Return a pointer to the node that interpolates the value-id-th unknown at local coordinate s...
void assemble_local_to_lagrangian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local to lagrangian coordinates, given the derivatives of the shape function Overload the standard version to use the hanging information for the lagrangian coordinates.
void disable_p_refinement()
Suppress of any refinement for this element.
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
virtual void setup_hanging_nodes(Vector< std::ofstream * > &output_stream)
Mark up any hanging nodes that arise as a result of non-uniform refinement. Any hanging nodes will be...
virtual RefineableElement * root_element_pt()
Pointer to the root element in refinement hierarchy (must be implemented in specific elements that do...
virtual ~PRefineableElement()
Destructor, empty.
virtual void further_build()
Further build: e.g. deal with interpolation of internal values.
std::map< Node *, int > * Local_hang_eqn
Storage for local equation numbers of hanging node variables (values stored at master nodes)...
unsigned Refine_level
Refinement level.
virtual void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
Assign the local equation numbers for Data stored at the nodes Virtual so that it can be overloaded b...
Definition: elements.cc:3473
unsigned ngeom_data() const
The number of geometric data affecting a RefineableSolidFiniteElement is the positional Data of all n...
void operator=(const PRefineableElement &)
Broken assignment operator.
virtual void unbuild()
Unbuild the element, i.e. mark the nodes that were created during its creation for possible deletion...
bool To_be_refined
Flag for refinement.
bool is_undeformed_macro_element_used_for_new_lagrangian_coords() const
Return the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be...
unsigned nshape_controlling_nodes()
Number of shape-controlling nodes = the number of non-hanging nodes plus the number of master nodes a...
virtual void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes (e...
A Rank 3 Tensor class.
Definition: matrices.h:1337
virtual unsigned ninterpolating_node(const int &value_id)
Return the number of nodes that are used to interpolate the value_id-th unknown. Default is to assume...
unsigned refinement_level() const
Return the Refinement level.
virtual void assign_solid_local_eqn_numbers(const bool &store_local_dof)
Assigns local equation numbers for the generic solid local equation numbering schemes. If the boolean flag is true the the degrees of freedom are stored in Dof_pt.
Definition: elements.cc:6572
virtual void fill_in_jacobian_from_nodal_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from the nodal degrees of freedom using finite difference...
long Number
Global element number – for plotting/validation purposes.
static char t char * s
Definition: cfortran.h:572
void disable_refinement()
Suppress of any refinement for this element.
Tree * tree_pt()
Access function: Pointer to quadtree representation of this element.
bool To_be_p_unrefined
Flag for unrefinement.
unsigned p_order() const
Access function to P_order (const version)
virtual void further_build()
Further build: Pass the father's Use_undeformed_macro_element_for_new_lagrangian_coords flag down...
void deselect_for_p_refinement()
Deselect the element for p-refinement.
void set_tree_pt(Tree *my_tree_pt)
Set pointer to quadtree representation of this element.
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3804
bool Sons_to_be_unrefined
Flag for unrefinement.
RefineableElement(const RefineableElement &)
Broken copy constructor.
bool To_be_p_refined
Flag for p-refinement.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
virtual void build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt, bool &was_already_built, std::ofstream &new_nodes_file)=0
Interface to function that builds the element: i.e. construct the nodes, assign their positions...
void assemble_local_to_eulerian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives of the mapping from local to Eulerian coordi...
void deselect_sons_for_unrefinement()
No unrefinement will be performed by merging the four sons of this element.
void fill_in_jacobian_from_solid_position_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute element residual Vector and element Jacobian matrix corresponding to the solid positions...
bool p_refinement_is_enabled()
Flag to indicate suppression of any refinement.
void assign_solid_local_eqn_numbers(const bool &store_local_dof_pt)
Overload the local equation numbers for Data stored as part of solid nodes to include hanging node da...
void get_father_at_refinement_level(unsigned &refinement_level, RefineableElement *&father_at_reflevel_pt)
Return a pointer to the "father" element at the specified refinement level.
static double & max_integrity_tolerance()
Max. allowed discrepancy in element integrity check.
bool refinement_is_enabled()
Flag to indicate suppression of any refinement.
virtual unsigned required_nsons() const
Set the number of sons that can be constructed by the element The default is none.
RefineableElement()
Constructor, calls the FiniteElement constructor and initialises the member data. ...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
virtual void p_refine(const int &inc, Mesh *const &mesh_pt, GeneralisedElement *const &clone_pt)=0
p-refine the element
void disable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
Unset the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be ...
void enable_use_of_undeformed_macro_element_for_new_lagrangian_coords()
Set the flag deciding if the Lagrangian coordinates of newly-created interior SolidNodes are to be de...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
void set_refinement_level(const int &refine_level)
Set the refinement level.
void assemble_local_to_eulerian_jacobian(const DShape &dpsids, DenseMatrix< double > &jacobian) const
Assemble the jacobian matrix for the mapping from local.
std::map< Node *, unsigned > Shape_controlling_node_lookup
Lookup scheme for unique number associated with any of the nodes that actively control the shape of t...
virtual void rebuild_from_sons(Mesh *&mesh_pt)=0
Rebuild the element, e.g. set internal values in line with those of the sons that have now merged...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Broken function – this shouldn't really be needed.
void set_number(const long &mynumber)
Set element number (for debugging/plotting)
void enable_p_refinement()
Emnable refinement for this element.
void select_for_refinement()
Select the element for refinement.
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
bool Refinement_is_enabled
Flag to indicate suppression of any refinement.
int get_node_number(Node *const &node_pt) const
Return the number of the node *node_pt if this node is in the element, else return -1;...
Definition: elements.cc:3735
std::map< Node *, unsigned > shape_controlling_node_lookup()
Return lookup scheme for unique number associated with any of the nodes that actively control the sha...
bool to_be_refined()
Has the element been selected for refinement?
void assign_nodal_local_eqn_numbers(const bool &store_local_dof_pt)
Overload the function that assigns local equation numbers for the Data stored at the nodes so that ha...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
virtual void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
Return the basis functions that are used to interpolate the value_id-th unknown. By default assume is...
virtual double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
Return the local one dimensional fraction of the n1d-th node in the direction of the local coordinate...
virtual void check_integrity(double &max_error)=0
Check the integrity of the element: Continuity of positions values, etc. Essentially, check that the approximation of the functions is consistent when viewed from both sides of the element boundaries Must be overloaded for each different geometric element.
double local_to_eulerian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates given the derivatives of the shape functions...
long number() const
Element number (for debugging/plotting)
std::map< Node *, DenseMatrix< int > > Local_position_hang_eqn
Storage for local equation numbers of hanging node variables associated with nodal positions...
p-refineable version of RefineableElement
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
virtual void pre_build(Mesh *&mesh_pt, Vector< Node * > &new_node_pt)
Pre-build the element.
PRefineableElement(const PRefineableElement &)
Broken copy constructor.
virtual void deactivate_element()
Final operations that must be performed when the element is no longer active in the mesh...
bool nodes_built()
Return true if all the nodes have been built, false if not.
bool sons_to_be_unrefined()
Has the element been selected for unrefinement?
unsigned P_order
The polynomial expansion order of the elemental basis functions.
void assign_solid_hanging_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers to the hanging values associated with positions or additional solid val...
virtual ~RefineableElement()
Destructor, delete the allocated storage for the hanging equations.
double local_to_lagrangian_mapping_diagonal(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Lagrangian coordinates given the derivatives of the shape functio...
DenseMatrix< int > & local_position_hang_eqn(Node *const &node_pt)
Access the local equation number of of hanging node variables associated with nodal positions...
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:221
void set_non_obsolete()
Mark node as non-obsolete.
Definition: nodes.h:1350
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1810
virtual void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Broken function – this shouldn't really be needed.
SolidFiniteElement class.
Definition: elements.h:3320
virtual unsigned ninterpolating_node_1d(const int &value_id)
Return the number of nodes in a one_d direction that are used to interpolate the value_id-th unknown...
void select_for_p_unrefinement()
Select the element for p-unrefinement.
void assemble_local_to_lagrangian_jacobian2(const DShape &d2psids, DenseMatrix< double > &jacobian2) const
Assemble the the "jacobian" matrix of second derivatives, given the second derivatives of the shape f...
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...
void enable_refinement()
Emnable refinement for this element.
A general mesh class.
Definition: mesh.h:74
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:143
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the position Data to the set that's p...
void rebuild_from_sons(Mesh *&mesh_pt)
Broken function – this shouldn't really be needed.