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// Copyright (C) 2006-2024 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for base class that specified common interfaces
27 // used in elements with moving Nodes (Spine,Algebraic,MacroElementNodeUpdate)
28 
29 // Include guards to prevent multiple inclusion of the header
30 #ifndef OOMPH_ELEMENT_WITH_MOVING_NODES
31 #define OOMPH_ELEMENT_WITH_MOVING_NODES
32 
33 
34 // oomph-lib headers
35 #include "elements.h"
36 
37 namespace oomph
38 {
39  /// ///////////////////////////////////////////////////////////////////////
40  /// ///////////////////////////////////////////////////////////////////////
41  /// ///////////////////////////////////////////////////////////////////////
42 
43 
44  //=======================================================================
45  /// A policy class that serves to establish the
46  /// common interfaces for elements that contain moving nodes. This
47  /// class provides storage for the geometric data that affect the
48  /// update of all the nodes of the element, i.e. USUALLY
49  /// all data that are using during a call
50  /// to the Element's node_update() function. In some cases
51  /// (e.g. FluidInterfaceEdge elements), node_update() is overloaded to
52  /// perform an update of the bulk element, in which case the additional
53  /// bulk geometric data become external data of the element and the
54  /// function GeneralisedElement::update_in_external_fd(i) is overloaded
55  /// to also perform the bulk node update.
56  /// The storage is populated
57  /// during the assignment of the equation numbers via the
58  /// complete_setup_of_dependencies() function and then local equations
59  /// numbers are assigned to these data, accessible via
60  /// geometric_data_local_eqn(n,i). Finally, a function is provided
61  /// that calculates the terms in the jacobian matrix by due to these
62  /// geometric data by finite differences.
63  //=======================================================================
64  class ElementWithMovingNodes : public virtual FiniteElement
65  {
66  public:
67  /// Public enumeration to choose method for computing shape derivatives
68  enum
69  {
73  };
74 
75  /// Constructor
81  // hierher: Anything other than the fd-based method is currently broken;
82  // at least for refineable elements -- this all needs to be checked
83  // VERY carefully again (see instructions in commit log). Until this
84  // has all been fixed/re-validated, we've frozen the default assignment
85  // (by breaking the functions that change the setting). Ultimately,
86  // we should use:
87  // Method_for_shape_derivs(Shape_derivs_by_fastest_method)
88  {
89  }
90 
91  /// Broken copy constructor
93 
94  /// Broken assignment operator
95  void operator=(const ElementWithMovingNodes&) = delete;
96 
97  /// Function to describe the local dofs of the element. The ostream
98  /// specifies the output stream to which the description
99  /// is written; the string stores the currently
100  /// assembled output that is ultimately written to the
101  /// output stream by Data::describe_dofs(...); it is typically
102  /// built up incrementally as we descend through the
103  /// call hierarchy of this function when called from
104  /// Problem::describe_dofs(...)
105  void describe_local_dofs(std::ostream& out,
106  const std::string& current_string) const;
107 
108  /// Virtual destructor (clean up and allocated memory)
110  {
112  {
113  delete[] Geometric_data_local_eqn[0];
114  delete[] Geometric_data_local_eqn;
115  }
116  }
117 
118  /// Number of geometric dofs
119  unsigned ngeom_dof() const
120  {
121  return Ngeom_dof;
122  }
123 
124  /// Return the local equation number corresponding to the i-th
125  /// value at the n-th geometric data object.
126  inline int geometric_data_local_eqn(const unsigned& n, const unsigned& i)
127  {
128 #ifdef RANGE_CHECKING
129  unsigned n_data = Geom_data_pt.size();
130  if (n >= n_data)
131  {
132  std::ostringstream error_message;
133  error_message << "Range Error: Data number " << n
134  << " is not in the range (0," << n_data - 1 << ")";
135  throw OomphLibError(error_message.str(),
136  OOMPH_CURRENT_FUNCTION,
137  OOMPH_EXCEPTION_LOCATION);
138  }
139  else
140  {
141  unsigned n_value = Geom_data_pt[n]->nvalue();
142  if (i >= n_value)
143  {
144  std::ostringstream error_message;
145  error_message << "Range Error: value " << i << " at data " << n
146  << " is not in the range (0," << n_value - 1 << ")";
147  throw OomphLibError(error_message.str(),
148  OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151  }
152 #endif
153 #ifdef PARANOID
154  // Check that the equations have been allocated
155  if (Geometric_data_local_eqn == 0)
156  {
157  throw OomphLibError(
158  "Geometric data local equation numbers have not been allocated",
159  OOMPH_CURRENT_FUNCTION,
160  OOMPH_EXCEPTION_LOCATION);
161  }
162 #endif
163  return Geometric_data_local_eqn[n][i];
164  }
165 
166 
167  /// Return a set of all geometric data associated with the element
169  std::set<Data*>& unique_geom_data_pt);
170 
171  /// Specify Data that affects the geometry of the element
172  /// by adding the element's geometric Data to the set that's passed in.
173  /// (This functionality is required in FSI problems; set is used to
174  /// avoid double counting).
175  void identify_geometric_data(std::set<Data*>& geometric_data_pt)
176  {
177  // Loop over the node update data and add to the set
178  const unsigned n_geom_data = Geom_data_pt.size();
179  for (unsigned n = 0; n < n_geom_data; n++)
180  {
181  geometric_data_pt.insert(Geom_data_pt[n]);
182  }
183  }
184 
185  /// Return whether shape derivatives are evaluated by fd
187  {
189  }
190 
191  /// Insist that shape derivatives are always
192  /// evaluated by fd (using
193  /// FiniteElement::get_dresidual_dnodal_coordinates())
195  {
197  }
198 
199  /// Insist that shape derivatives are always
200  /// evaluated using the overloaded version of this function
201  /// that may have been implemented in a derived class.
202  /// (The default behaviour will still be finite differences unless the
203  /// function has actually been overloaded
205  {
207  }
208 
209  /// Evaluate shape derivatives by direct finite differencing
211  {
212  // hierher: Harmless; simply re-sets the current (and currently only
213  // legal) default
215  }
216 
217  /// Evaluate shape derivatives by chain rule. Currently disabled by
218  /// default because it's broken; can re-enable use by setting optional
219  /// boolean to true.
221  const bool& i_know_what_i_am_doing = false)
222  {
223  if (!i_know_what_i_am_doing)
224  {
225  std::ostringstream error_message;
226  error_message
227  << "Evaluation of shape derivatives by chain rule is currently \n"
228  << "disabled because it's broken, at least for refineable \n"
229  << "elements. This all needs to be checked again very carefully\n"
230  << "following the instructions in the commit log\n"
231  << "If you know what you're doing and want to force this "
232  "methodology\n"
233  << "call this function with the optional boolean set to true.\n";
234  throw OomphLibError(
235  error_message.str(),
236  "ElementWithMovingNodes::evaluate_shape_derivs_by_chain_rule()",
237  OOMPH_EXCEPTION_LOCATION);
238  }
240  }
241 
242  /// Evaluate shape derivatives by (anticipated) fastest method.
243  /// Currently disabled by default because it's broken; can re-enable
244  /// use by setting optional boolean to true.
246  const bool& i_know_what_i_am_doing = false)
247  {
248  if (!i_know_what_i_am_doing)
249  {
250  std::ostringstream error_message;
251  error_message
252  << "Evaluation of shape derivatives by fastest method is currently \n"
253  << "disabled because it's broken, at least for refineable \n"
254  << "elements. This all needs to be checked again very carefully\n"
255  << "following the instructions in the commit log\n"
256  << "If you know what you're doing and want to force this "
257  "methodology\n"
258  << "call this function with the optional boolean set to true.\n";
259  throw OomphLibError(
260  error_message.str(),
261  "ElementWithMovingNodes::evaluate_shape_derivs_by_fastest_method()",
262  OOMPH_EXCEPTION_LOCATION);
263  }
265  }
266 
267  /// Access to method (enumerated flag) for determination of shape derivs
269  {
271  }
272 
273  /// Bypass the call to fill_in_jacobian_from_geometric_data
275  {
277  }
278 
279  /// Do not bypass the call to fill_in_jacobian_from_geometric_data
281  {
283  }
284 
285  /// Test whether the call to fill_in_jacobian_from_geometric_data is
286  /// bypassed
288  {
290  }
291 
292  protected:
293  /// Compute derivatives of the nodal coordinates w.r.t.
294  /// to the geometric dofs. Default implementation by FD can be overwritten
295  /// for specific elements.
296  /// dnodal_coordinates_dgeom_dofs(l,i,j) = dX_{ij} / d s_l
298  RankThreeTensor<double>& dnodal_coordinates_dgeom_dofs);
299 
300  /// Construct the vector of (unique) geometric data
302 
303  /// Assign local equation numbers for the geometric Data in the element
304  /// If the boolean argument is true then the degrees of freedom are stored
305  /// in Dof_pt
307  const bool& store_local_dof_pt);
308 
309  /// Calculate the contributions to the Jacobian matrix from the
310  /// geometric data. This version
311  /// assumes that the (full) residuals vector has already been calculated
312  /// and is passed in as the first argument -- needed in case
313  /// the derivatives are computed by FD.
315  DenseMatrix<double>& jacobian);
316 
317  /// Calculate the contributions to the Jacobian matrix from the
318  /// geometric data. This version computes
319  /// the residuals vector before calculating the Jacobian terms
321  {
323  {
324  // Allocate storage for the residuals
325  const unsigned n_dof = ndof();
326  Vector<double> residuals(n_dof);
327 
328  // Get the residuals for the entire element
329  get_residuals(residuals);
330 
331  // Call the jacobian calculation
332  fill_in_jacobian_from_geometric_data(residuals, jacobian);
333  }
334  }
335 
336 
337  /// Vector that stores pointers to all Data that affect the
338  /// node update operations, i.e. the variables that can affect
339  /// the position of the node.
341 
342  public:
343  /// Return the number of geometric data upon which the shape
344  /// of the element depends
345  unsigned ngeom_data() const
346  {
347  return Geom_data_pt.size();
348  }
349 
350  private:
351  /// Array to hold local eqn number information for the
352  /// geometric Data variables
354 
355  /// Number of geometric dofs (computed on the fly when
356  /// equation numbers are set up)
357  unsigned Ngeom_dof;
358 
359  /// Set flag to true to bypass calculation of Jacobain entries
360  /// resulting from geometric data.
362 
363  /// Boolean to decide if shape derivatives are to be evaluated
364  /// by fd (using FiniteElement::get_dresidual_dnodal_coordinates())
365  /// or analytically, using the overloaded version of this function
366  /// in this class.
368 
369  /// Choose method for evaluation of shape derivatives
370  /// (this takes one of the values in the enumeration)
372  };
373 
374 
375  //===============================================================
376  /// Specific implementation of the class for specified element
377  /// and node type.
378  //==============================================================
379  template<class ELEMENT, class NODE_TYPE>
380  class ElementWithSpecificMovingNodes : public ELEMENT,
382  {
383  public:
384  /// Function to describe the local dofs of the element. The ostream
385  /// specifies the output stream to which the description
386  /// is written; the string stores the currently
387  /// assembled output that is ultimately written to the
388  /// output stream by Data::describe_dofs(...); it is typically
389  /// built up incrementally as we descend through the
390  /// call hierarchy of this function when called from
391  /// Problem::describe_dofs(...)
392  void describe_local_dofs(std::ostream& out,
393  const std::string& current_string) const
394  {
395  ELEMENT::describe_local_dofs(out, current_string);
396  ElementWithMovingNodes::describe_local_dofs(out, current_string);
397  }
398 
399  /// Constructor, call the constructor of the base element
401 
402  /// Constructor used for face elements
404  const int& face_index)
405  : ELEMENT(element_pt, face_index), ElementWithMovingNodes()
406  {
407  }
408 
409  /// Empty Destructor,
411 
412  /// Unique final overrider for describe_dofs
413  void describe_local_dofs(std::ostream& out, std::string& curr_str)
414  {
416  ELEMENT::describe_local_dofs(out, curr_str);
417  }
418 
419 
420  /// Overload the node assignment routine to assign nodes of the
421  /// appropriate type.
422  Node* construct_node(const unsigned& n)
423  {
424  // Assign a node to the local node pointer
425  // The dimension and number of values are taken from internal element data
426  // The number of timesteps to be stored comes from the problem!
427  this->node_pt(n) = new NODE_TYPE(this->nodal_dimension(),
428  this->nnodal_position_type(),
429  this->required_nvalue(n));
430  // Now return a pointer to the node, so that the mesh can find it
431  return this->node_pt(n);
432  }
433 
434  /// Overloaded node allocation for unsteady problems
435  Node* construct_node(const unsigned& n, TimeStepper* const& time_stepper_pt)
436  {
437  // Assign a node to the local node pointer
438  // The dimension and number of values are taken from internal element data
439  // The number of timesteps to be stored comes from the problem!
440  this->node_pt(n) = new NODE_TYPE(time_stepper_pt,
441  this->nodal_dimension(),
442  this->nnodal_position_type(),
443  this->required_nvalue(n));
444  // Now return a pointer to the node, so that the mesh can find it
445  return this->node_pt(n);
446  }
447 
448  /// Overload the node assignment routine to assign boundary nodes
449  Node* construct_boundary_node(const unsigned& n)
450  {
451  // Assign a node to the local node pointer
452  // The dimension and number of values are taken from internal element data
453  // The number of timesteps to be stored comes from the problem!
454  this->node_pt(n) =
456  this->nnodal_position_type(),
457  this->required_nvalue(n));
458  // Now return a pointer to the node, so that the mesh can find it
459  return this->node_pt(n);
460  }
461 
462  /// Overloaded boundary node allocation for unsteady problems
463  Node* construct_boundary_node(const unsigned& n,
465  {
466  // Assign a node to the local node pointer
467  // The dimension and number of values are taken from internal element data
468  // The number of timesteps to be stored comes from the problem!
469  this->node_pt(n) =
471  this->nodal_dimension(),
472  this->nnodal_position_type(),
473  this->required_nvalue(n));
474  // Now return a pointer to the node, so that the mesh can find it
475  return this->node_pt(n);
476  }
477 
478 
479  /// Complete the setup of additional dependencies. Overloads
480  /// empty virtual function in GeneralisedElement to determine the "geometric
481  /// Data", i.e. the Data that affects the element's shape.
482  /// This function is called (for all elements) at the very beginning of the
483  /// equation numbering procedure to ensure that all dependencies
484  /// are accounted for.
486  {
487  // Call function of underlying element
488  ELEMENT::complete_setup_of_dependencies();
489 
490  // Call function of the element with moving nodes
492  }
493 
494 
495  /// Assign local equation numbers for the underlying element, then
496  /// deal with the additional geometric dofs
497  void assign_all_generic_local_eqn_numbers(const bool& store_local_dof_pt)
498  {
499  // Call the generic local equation numbering scheme of the ELEMENT
500  ELEMENT::assign_all_generic_local_eqn_numbers(store_local_dof_pt);
502  store_local_dof_pt);
503  }
504 
505  /// Compute the element's residuals vector and jacobian matrix
506  void get_jacobian(Vector<double>& residuals, DenseMatrix<double>& jacobian)
507  {
508  /// Call the element's get jacobian function
509  ELEMENT::get_jacobian(residuals, jacobian);
510 
511  // Now call the additional geometric Jacobian terms
512  this->fill_in_jacobian_from_geometric_data(jacobian);
513  }
514 
515  /// Compute the element's residuals vector and jacobian matrix
517  DenseMatrix<double>& jacobian,
518  DenseMatrix<double>& mass_matrix)
519  {
520  /// Call the element's get jacobian function
521  ELEMENT::get_jacobian_and_mass_matrix(residuals, jacobian, mass_matrix);
522 
523  // Now call the additional geometric Jacobian terms
524  this->fill_in_jacobian_from_geometric_data(jacobian);
525  }
526 
527 
528  private:
529  };
530 
531 
532 } // namespace oomph
533 #endif
cstr elem_len * i
Definition: cfortran.h:603
A template Class for BoundaryNodes; that is Nodes that MAY live on the boundary of a Mesh....
Definition: nodes.h:2242
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
unsigned ngeom_dof() const
Number of geometric dofs.
Vector< Data * > Geom_data_pt
Vector that stores pointers to all Data that affect the node update operations, i....
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.
unsigned Ngeom_dof
Number of geometric dofs (computed on the fly when equation numbers are set up)
bool are_dresidual_dnodal_coordinates_always_evaluated_by_fd() const
Return whether shape derivatives are evaluated by fd.
void operator=(const ElementWithMovingNodes &)=delete
Broken assignment operator.
virtual ~ElementWithMovingNodes()
Virtual destructor (clean up and allocated memory)
bool Evaluate_dresidual_dnodal_coordinates_by_fd
Boolean to decide if shape derivatives are to be evaluated by fd (using FiniteElement::get_dresidual_...
ElementWithMovingNodes(const ElementWithMovingNodes &)=delete
Broken copy constructor.
void disable_bypass_fill_in_jacobian_from_geometric_data()
Do not bypass the call to fill_in_jacobian_from_geometric_data.
void enable_bypass_fill_in_jacobian_from_geometric_data()
Bypass the call to fill_in_jacobian_from_geometric_data.
virtual void assign_all_generic_local_eqn_numbers(const bool &store_local_dof_pt)
Assign local equation numbers for the geometric Data in the element If the boolean argument is true t...
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
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 ...
bool Bypass_fill_in_jacobian_from_geometric_data
Set flag to true to bypass calculation of Jacobain entries resulting from geometric data.
void complete_setup_of_dependencies()
Construct the vector of (unique) geometric data.
int Method_for_shape_derivs
Choose method for evaluation of shape derivatives (this takes one of the values in the enumeration)
bool is_fill_in_jacobian_from_geometric_data_bypassed() const
Test whether the call to fill_in_jacobian_from_geometric_data is bypassed.
void evaluate_shape_derivs_by_direct_fd()
Evaluate shape derivatives by direct finite differencing.
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 ...
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 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 fill_in_jacobian_from_geometric_data(DenseMatrix< double > &jacobian)
Calculate the contributions to the Jacobian matrix from the geometric data. This version computes the...
unsigned ngeom_data() const
Return the number of geometric data upon which the shape of the element depends.
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 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...
int ** Geometric_data_local_eqn
Array to hold local eqn number information for the geometric Data variables.
void disable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated using the overloaded version of this function that...
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'...
void enable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated by fd (using FiniteElement::get_dresidual_dnodal_c...
Specific implementation of the class for specified element and node type.
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.
ElementWithSpecificMovingNodes(FiniteElement *const &element_pt, const int &face_index)
Constructor used for face elements.
Node * construct_node(const unsigned &n)
Overload the node assignment routine to assign nodes of the appropriate type.
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...
void describe_local_dofs(std::ostream &out, std::string &curr_str)
Unique final overrider for describe_dofs.
Node * construct_boundary_node(const unsigned &n)
Overload the node assignment routine to assign boundary nodes.
ElementWithSpecificMovingNodes()
Constructor, call the constructor of the base element.
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 ...
void get_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residuals vector and jacobian matrix.
Node * construct_boundary_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded boundary node allocation for unsteady problems.
void complete_setup_of_dependencies()
Complete the setup of additional dependencies. Overloads empty virtual function in GeneralisedElement...
Node * construct_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded node allocation for unsteady problems.
A general Finite Element class.
Definition: elements.h:1317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2179
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2467
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:2459
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2488
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:839
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:984
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...