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-2022 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
37namespace 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 //=======================================================================
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];
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
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);
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
507 {
508 /// Call the element's get jacobian function
509 ELEMENT::get_jacobian(residuals, jacobian);
510
511 // Now call the additional geometric Jacobian terms
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
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...
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.
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.
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.
Node * construct_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded node allocation for unsteady problems.
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.
Node * construct_boundary_node(const unsigned &n)
Overload the node assignment routine to assign boundary nodes.
Node * construct_boundary_node(const unsigned &n, TimeStepper *const &time_stepper_pt)
Overloaded boundary node allocation for unsteady problems.
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.
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.
void complete_setup_of_dependencies()
Complete the setup of additional dependencies. Overloads empty virtual function in GeneralisedElement...
A general Finite Element class.
Definition: elements.h:1313
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
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:2455
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
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:980
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...