element_with_moving_nodes.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 base class that specified common interfaces
31 //used in elements with moving Nodes (Spine,Algebraic,MacroElementNodeUpdate)
32 
33 //Include guards to prevent multiple inclusion of the header
34 #ifndef OOMPH_ELEMENT_WITH_MOVING_NODES
35 #define OOMPH_ELEMENT_WITH_MOVING_NODES
36 
37 
38 //oomph-lib headers
39 #include "elements.h"
40 
41 namespace oomph
42 {
43 
44 
45 
46 
47 
48 
49 //////////////////////////////////////////////////////////////////////////
50 //////////////////////////////////////////////////////////////////////////
51 //////////////////////////////////////////////////////////////////////////
52 
53 
54 
55  //=======================================================================
56  /// \short A policy class that serves to establish the
57  /// common interfaces for elements that contain moving nodes. This
58  /// class provides storage for the geometric data that affect the
59  /// update of all the nodes of the element, i.e. USUALLY
60  /// all data that are using during a call
61  /// to the Element's node_update() function. In some cases
62  /// (e.g. FluidInterfaceEdge elements), node_update() is overloaded to
63  /// perform an update of the bulk element, in which case the additional
64  /// bulk geometric data become external data of the element and the
65  /// function GeneralisedElement::update_in_external_fd(i) is overloaded
66  /// to also perform the bulk node update.
67  /// The storage is populated
68  /// during the assignment of the equation numbers via the
69  /// complete_setup_of_dependencies() function and then local equations
70  /// numbers are assigned to these data, accessible via
71  /// geometric_data_local_eqn(n,i). Finally, a function is provided
72  /// that calculates the terms in the jacobian matrix by due to these
73  /// geometric data by finite differences.
74  //=======================================================================
75  class ElementWithMovingNodes : public virtual FiniteElement
76 {
77  public:
78 
79  /// Public enumeration to choose method for computing shape derivatives
83 
84  ///Constructor
90  // hierher: Anything other than the fd-based method is currently broken;
91  // at least for refineable elements -- this all needs to be checked
92  // VERY carefully again (see instructions in commit log). Until this
93  // has all been fixed/re-validated, we've frozen the default assignment
94  // (by breaking the functions that change the setting). Ultimately,
95  // we should use:
96  // Method_for_shape_derivs(Shape_derivs_by_fastest_method)
97  {}
98 
99  /// Broken copy constructor
101  {
102  BrokenCopy::broken_copy("ElementWithMovingNodes");
103  }
104 
105  /// Broken assignment operator
107  {
108  BrokenCopy::broken_assign("ElementWithMovingNodes");
109  }
110 
111  /// \short Function to describe the local dofs of the element. The ostream
112  /// specifies the output stream to which the description
113  /// is written; the string stores the currently
114  /// assembled output that is ultimately written to the
115  /// output stream by Data::describe_dofs(...); it is typically
116  /// built up incrementally as we descend through the
117  /// call hierarchy of this function when called from
118  /// Problem::describe_dofs(...)
119  void describe_local_dofs(std::ostream& out,
120  const std::string& current_string) const;
121 
122  /// Virtual destructor (clean up and allocated memory)
124  {
126  {
127  delete[] Geometric_data_local_eqn[0];
128  delete[] Geometric_data_local_eqn;
129  }
130  }
131 
132  /// Number of geometric dofs
133  unsigned ngeom_dof() const {return Ngeom_dof;}
134 
135  /// \short Return the local equation number corresponding to the i-th
136  /// value at the n-th geometric data object.
137  inline int geometric_data_local_eqn(const unsigned &n, const unsigned &i)
138  {
139 #ifdef RANGE_CHECKING
140  unsigned n_data = Geom_data_pt.size();
141  if(n >= n_data)
142  {
143  std::ostringstream error_message;
144  error_message << "Range Error: Data number " << n
145  << " is not in the range (0,"
146  << n_data-1 << ")";
147  throw OomphLibError(error_message.str(),
148  OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151  else
152  {
153  unsigned n_value = Geom_data_pt[n]->nvalue();
154  if(i >= n_value)
155  {
156  std::ostringstream error_message;
157  error_message << "Range Error: value " << i << " at data " << n
158  << " is not in the range (0,"
159  << n_value -1 << ")";
160  throw OomphLibError(error_message.str(),
161  OOMPH_CURRENT_FUNCTION,
162  OOMPH_EXCEPTION_LOCATION);
163  }
164  }
165 #endif
166 #ifdef PARANOID
167  //Check that the equations have been allocated
169  {
170  throw OomphLibError(
171  "Geometric data local equation numbers have not been allocated",
172  OOMPH_CURRENT_FUNCTION,
173  OOMPH_EXCEPTION_LOCATION);
174  }
175 #endif
176  return Geometric_data_local_eqn[n][i];
177  }
178 
179 
180  ///Return a set of all geometric data associated with the element
181  void assemble_set_of_all_geometric_data(std::set<Data*> &unique_geom_data_pt);
182 
183  /// \short Specify Data that affects the geometry of the element
184  /// by adding the element's geometric Data to the set that's passed in.
185  /// (This functionality is required in FSI problems; set is used to
186  /// avoid double counting).
187  void identify_geometric_data(std::set<Data*> &geometric_data_pt)
188  {
189  //Loop over the node update data and add to the set
190  const unsigned n_geom_data = Geom_data_pt.size();
191  for(unsigned n=0;n<n_geom_data;n++)
192  {
193  geometric_data_pt.insert(Geom_data_pt[n]);
194  }
195  }
196 
197  /// \short Return whether shape derivatives are evaluated by fd
200 
201  /// \short Insist that shape derivatives are always
202  /// evaluated by fd (using FiniteElement::get_dresidual_dnodal_coordinates())
205 
206  /// \short Insist that shape derivatives are always
207  /// evaluated using the overloaded version of this function
208  /// that may have been implemented in a derived class.
209  /// (The default behaviour will still be finite differences unless the
210  /// function has actually been overloaded
213 
214  /// Evaluate shape derivatives by direct finite differencing
216  {
217  // hierher: Harmless; simply re-sets the current (and currently only
218  // legal) default
220  }
221 
222  /// \short Evaluate shape derivatives by chain rule. Currently disabled by
223  /// default because it's broken; can re-enable use by setting optional
224  /// boolean to true.
226  const bool& i_know_what_i_am_doing=false)
227  {
228  if (!i_know_what_i_am_doing)
229  {
230  std::ostringstream error_message;
231  error_message
232  << "Evaluation of shape derivatives by chain rule is currently \n"
233  << "disabled because it's broken, at least for refineable \n"
234  << "elements. This all needs to be checked again very carefully\n"
235  << "following the instructions in the commit log\n"
236  << "If you know what you're doing and want to force this methodology\n"
237  << "call this function with the optional boolean set to true.\n";
238  throw OomphLibError(
239  error_message.str(),
240  "ElementWithMovingNodes::evaluate_shape_derivs_by_chain_rule()",
241  OOMPH_EXCEPTION_LOCATION);
242  }
244  }
245 
246  /// \short Evaluate shape derivatives by (anticipated) fastest method.
247  /// Currently disabled by default because it's broken; can re-enable
248  /// use by setting optional boolean to true.
250  const bool& i_know_what_i_am_doing=false)
251  {
252  if (!i_know_what_i_am_doing)
253  {
254  std::ostringstream error_message;
255  error_message
256  << "Evaluation of shape derivatives by fastest method is currently \n"
257  << "disabled because it's broken, at least for refineable \n"
258  << "elements. This all needs to be checked again very carefully\n"
259  << "following the instructions in the commit log\n"
260  << "If you know what you're doing and want to force this methodology\n"
261  << "call this function with the optional boolean set to true.\n";
262  throw OomphLibError(
263  error_message.str(),
264  "ElementWithMovingNodes::evaluate_shape_derivs_by_fastest_method()",
265  OOMPH_EXCEPTION_LOCATION);
266  }
268  }
269 
270  /// Access to method (enumerated flag) for determination of shape derivs
272 
273  /// Bypass the call to fill_in_jacobian_from_geometric_data
276 
277  /// Do not bypass the call to fill_in_jacobian_from_geometric_data
280 
281  /// Test whether the call to fill_in_jacobian_from_geometric_data is bypassed
284 
285  protected:
286 
287  /// \short Compute derivatives of the nodal coordinates w.r.t.
288  /// to the geometric dofs. Default implementation by FD can be overwritten
289  /// for specific elements.
290  /// dnodal_coordinates_dgeom_dofs(l,i,j) = dX_{ij} / d s_l
292  dnodal_coordinates_dgeom_dofs);
293 
294  /// Construct the vector of (unique) geometric data
296 
297  /// Assign local equation numbers for the geometric Data in the element
298  /// If the boolean argument is true then the degrees of freedom are stored
299  /// in Dof_pt
301  const bool &store_local_dof_pt);
302 
303  /// \short Calculate the contributions to the Jacobian matrix from the
304  /// geometric data. This version
305  /// assumes that the (full) residuals vector has already been calculated
306  /// and is passed in as the first argument -- needed in case
307  /// the derivatives are computed by FD.
309  Vector<double> &residuals, DenseMatrix<double> &jacobian);
310 
311  ///Calculate the contributions to the Jacobian matrix from the
312  ///geometric data. This version computes
313  ///the residuals vector before calculating the Jacobian terms
315  {
317  {
318  //Allocate storage for the residuals
319  const unsigned n_dof = ndof();
320  Vector<double> residuals(n_dof);
321 
322  //Get the residuals for the entire element
323  get_residuals(residuals);
324 
325  //Call the jacobian calculation
326  fill_in_jacobian_from_geometric_data(residuals,jacobian);
327  }
328  }
329 
330 
331  /// \short Vector that stores pointers to all Data that affect the
332  /// node update operations, i.e. the variables that can affect
333  /// the position of the node.
335 
336  public:
337 
338  /// \short Return the number of geometric data upon which the shape
339  /// of the element depends
340  unsigned ngeom_data() const {return Geom_data_pt.size();}
341 
342 private:
343 
344  /// \short Array to hold local eqn number information for the
345  /// geometric Data variables
347 
348  /// \short Number of geometric dofs (computed on the fly when
349  /// equation numbers are set up)
350  unsigned Ngeom_dof;
351 
352  /// \short Set flag to true to bypass calculation of Jacobain entries
353  /// resulting from geometric data.
355 
356  /// \short Boolean to decide if shape derivatives are to be evaluated
357  /// by fd (using FiniteElement::get_dresidual_dnodal_coordinates())
358  /// or analytically, using the overloaded version of this function
359  /// in this class.
361 
362  /// \short Choose method for evaluation of shape derivatives
363  /// (this takes one of the values in the enumeration)
365 
366 };
367 
368 
369 //===============================================================
370 /// Specific implementation of the class for specified element
371 /// and node type.
372 //==============================================================
373 template<class ELEMENT,class NODE_TYPE>
374  class ElementWithSpecificMovingNodes : public ELEMENT,
376 {
377 public:
378 
379  /// \short Function to describe the local dofs of the element. The ostream
380  /// specifies the output stream to which the description
381  /// is written; the string stores the currently
382  /// assembled output that is ultimately written to the
383  /// output stream by Data::describe_dofs(...); it is typically
384  /// built up incrementally as we descend through the
385  /// call hierarchy of this function when called from
386  /// Problem::describe_dofs(...)
387  void describe_local_dofs(std::ostream& out,
388  const std::string& current_string) const
389  {
390  ELEMENT::describe_local_dofs(out,current_string);
392  }
393 
394  /// Constructor, call the constructor of the base element
396  {}
397 
398  /// Constructor used for face elements
400  const int &face_index) :
401  ELEMENT(element_pt, face_index), ElementWithMovingNodes()
402  {}
403 
404  /// Empty Destructor,
406 
407  /// Unique final overrider for describe_dofs
408  void describe_local_dofs(std::ostream& out, std::string& curr_str)
409  {
411  ELEMENT::describe_local_dofs(out,curr_str);
412  }
413 
414 
415  /// \short Overload the node assignment routine to assign nodes of the
416  /// appropriate type.
417  Node* construct_node(const unsigned &n)
418  {
419  //Assign a node to the local node pointer
420  //The dimension and number of values are taken from internal element data
421  //The number of timesteps to be stored comes from the problem!
422  this->node_pt(n) =
423  new NODE_TYPE(this->nodal_dimension(),this->nnodal_position_type(),
424  this->required_nvalue(n));
425  //Now return a pointer to the node, so that the mesh can find it
426  return this->node_pt(n);
427  }
428 
429  /// Overloaded node allocation for unsteady problems
430  Node* construct_node(const unsigned &n,
432  {
433  //Assign a node to the local node pointer
434  //The dimension and number of values are taken from internal element data
435  //The number of timesteps to be stored comes from the problem!
436  this->node_pt(n) =
437  new NODE_TYPE(time_stepper_pt,this->nodal_dimension(),
438  this->nnodal_position_type(),
439  this->required_nvalue(n));
440  //Now return a pointer to the node, so that the mesh can find it
441  return this->node_pt(n);
442  }
443 
444  /// Overload the node assignment routine to assign boundary nodes
445  Node* construct_boundary_node(const unsigned &n)
446  {
447  //Assign a node to the local node pointer
448  //The dimension and number of values are taken from internal element data
449  //The number of timesteps to be stored comes from the problem!
450  this->node_pt(n) =
452  this->nnodal_position_type(),
453  this->required_nvalue(n));
454  //Now return a pointer to the node, so that the mesh can find it
455  return this->node_pt(n);
456  }
457 
458  /// Overloaded boundary node allocation for unsteady problems
459  Node* construct_boundary_node(const unsigned &n,
461  {
462  //Assign a node to the local node pointer
463  //The dimension and number of values are taken from internal element data
464  //The number of timesteps to be stored comes from the problem!
465  this->node_pt(n) =
467  this->nnodal_position_type(),
468  this->required_nvalue(n));
469  //Now return a pointer to the node, so that the mesh can find it
470  return this->node_pt(n);
471  }
472 
473 
474  /// \short Complete the setup of additional dependencies. Overloads
475  /// empty virtual function in GeneralisedElement to determine the "geometric
476  /// Data", i.e. the Data that affects the element's shape.
477  /// This function is called (for all elements) at the very beginning of the
478  /// equation numbering procedure to ensure that all dependencies
479  /// are accounted for.
481  {
482 
483  // Call function of underlying element
484  ELEMENT::complete_setup_of_dependencies();
485 
486  //Call function of the element with moving nodes
488  }
489 
490 
491  /// \short Assign local equation numbers for the underlying element, then
492  /// deal with the additional geometric dofs
493  void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
494  {
495  // Call the generic local equation numbering scheme of the ELEMENT
496  ELEMENT::assign_all_generic_local_eqn_numbers(store_local_dof_pt);
498  store_local_dof_pt);
499  }
500 
501  /// Compute the element's residuals vector and jacobian matrix
502  void get_jacobian(Vector<double> &residuals,
503  DenseMatrix<double> &jacobian)
504  {
505  ///Call the element's get jacobian function
506  ELEMENT::get_jacobian(residuals,jacobian);
507 
508  //Now call the additional geometric Jacobian terms
509  this->fill_in_jacobian_from_geometric_data(jacobian);
510  }
511 
512  /// Compute the element's residuals vector and jacobian matrix
514  DenseMatrix<double> &jacobian,
515  DenseMatrix<double> &mass_matrix)
516  {
517  ///Call the element's get jacobian function
518  ELEMENT::get_jacobian_and_mass_matrix(residuals,jacobian,mass_matrix);
519 
520  //Now call the additional geometric Jacobian terms
521  this->fill_in_jacobian_from_geometric_data(jacobian);
522  }
523 
524 
525  private:
526 };
527 
528 
529 }
530 #endif
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
virtual void get_residuals(Vector< double > &residuals)
Calculate the vector of residuals of the equations in the element. By default initialise the vector t...
Definition: elements.h:975
Vector< Data * > Geom_data_pt
Vector that stores pointers to all Data that affect the node update operations, i.e. the variables that can affect the position of the node.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
A policy class that serves to establish the common interfaces for elements that contain moving nodes...
bool Bypass_fill_in_jacobian_from_geometric_data
Set flag to true to bypass calculation of Jacobain entries resulting from geometric data...
void get_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residuals vector and jacobian matrix.
cstr elem_len * i
Definition: cfortran.h:607
bool Evaluate_dresidual_dnodal_coordinates_by_fd
Boolean to decide if shape derivatives are to be evaluated by fd (using FiniteElement::get_dresidual_...
Node * construct_boundary_node(const unsigned &n)
Overload the node assignment routine to assign boundary nodes.
bool is_fill_in_jacobian_from_geometric_data_bypassed() const
Test whether the call to fill_in_jacobian_from_geometric_data is bypassed.
virtual ~ElementWithMovingNodes()
Virtual destructor (clean up and allocated memory)
A general Finite Element class.
Definition: elements.h:1271
int ** Geometric_data_local_eqn
Array to hold local eqn number information for the geometric Data variables.
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2340
void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers for the underlying element, then deal with the additional geometric dof...
int Method_for_shape_derivs
Choose method for evaluation of shape derivatives (this takes one of the values in the enumeration) ...
void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residuals vector and jacobian matrix.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
unsigned ngeom_data() const
Return the number of geometric data upon which the shape of the element depends.
void evaluate_shape_derivs_by_fastest_method(const bool &i_know_what_i_am_doing=false)
Evaluate shape derivatives by (anticipated) fastest method. Currently disabled by default because it'...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
Node * construct_node(const unsigned &n)
Overload the node assignment routine to assign nodes of the appropriate type.
A Rank 3 Tensor class.
Definition: matrices.h:1337
Node * construct_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded node allocation for unsteady problems.
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh...
Definition: nodes.h:67
bool are_dresidual_dnodal_coordinates_always_evaluated_by_fd() const
Return whether shape derivatives are evaluated by fd.
Node * construct_boundary_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded boundary node allocation for unsteady problems.
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
int geometric_data_local_eqn(const unsigned &n, const unsigned &i)
Return the local equation number corresponding to the i-th value at the n-th geometric data object...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
ElementWithMovingNodes(const ElementWithMovingNodes &)
Broken copy constructor.
void disable_bypass_fill_in_jacobian_from_geometric_data()
Do not bypass the call to fill_in_jacobian_from_geometric_data.
virtual void get_dnodal_coordinates_dgeom_dofs(RankThreeTensor< double > &dnodal_coordinates_dgeom_dofs)
Compute derivatives of the nodal coordinates w.r.t. to the geometric dofs. Default implementation by ...
ElementWithSpecificMovingNodes()
Constructor, call the constructor of the base element.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values that must be stored at local node n by the element. The default is 0...
Definition: elements.h:2334
unsigned Ngeom_dof
Number of geometric dofs (computed on the fly when equation numbers are set up)
void assemble_set_of_all_geometric_data(std::set< Data * > &unique_geom_data_pt)
Return a set of all geometric data associated with the element.
void enable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated by fd (using FiniteElement::get_dresidual_dnodal_c...
void describe_local_dofs(std::ostream &out, std::string &curr_str)
Unique final overrider for describe_dofs.
void operator=(const ElementWithMovingNodes &)
Broken assignment operator.
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
void evaluate_shape_derivs_by_direct_fd()
Evaluate shape derivatives by direct finite differencing.
ElementWithSpecificMovingNodes(FiniteElement *const &element_pt, const int &face_index)
Constructor used for face elements.
void enable_bypass_fill_in_jacobian_from_geometric_data()
Bypass the call to fill_in_jacobian_from_geometric_data.
unsigned ngeom_dof() const
Number of geometric dofs.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void fill_in_jacobian_from_geometric_data(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the Jacobian matrix from the geometric data. This version assumes that...
void fill_in_jacobian_from_geometric_data(DenseMatrix< double > &jacobian)
void evaluate_shape_derivs_by_chain_rule(const bool &i_know_what_i_am_doing=false)
Evaluate shape derivatives by chain rule. Currently disabled by default because it's broken; can re-e...
void complete_setup_of_dependencies()
Construct the vector of (unique) geometric data.
void complete_setup_of_dependencies()
Complete the setup of additional dependencies. Overloads empty virtual function in GeneralisedElement...
void identify_geometric_data(std::set< Data * > &geometric_data_pt)
Specify Data that affects the geometry of the element by adding the element's geometric Data to the s...
void disable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated using the overloaded version of this function that...