algebraic_elements.cc
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-2023 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 #include "geom_objects.h"
27 #include "algebraic_elements.h"
29 
30 
31 namespace oomph
32 {
33  /// ////////////////////////////////////////////////////////////////////
34  /// ////////////////////////////////////////////////////////////////////
35  // Base class for Algebraic Elements
36  /// ////////////////////////////////////////////////////////////////////
37  /// ////////////////////////////////////////////////////////////////////
38 
39 
40  //========================================================================
41  /// Set up node update info for (newly created) algebraic node: Work out its
42  /// node update information by interpolation from its father element,
43  /// based on pointer to father element and its local coordinate in
44  /// the father element. We're creating the node update info
45  /// for update functions that are shared by all nodes in the
46  /// father element.
47  //========================================================================
49  Node*& node_pt,
50  const Vector<double>& s_father,
51  FiniteElement* father_el_pt) const
52  {
53  // Turn node into algebraic node
54  AlgebraicNode* alg_node_pt = dynamic_cast<AlgebraicNode*>(node_pt);
55 
56  // Get number of nodes in father element
57  unsigned nnod_father = father_el_pt->nnode();
58 
59  // Create map that stores the number of times a node update fct id
60  // has been encountered
61  std::map<int, unsigned> id_count;
62 
63  // Loop over all nodes in father element to extract node update ids
64  Vector<int> id;
65  for (unsigned j = 0; j < nnod_father; j++)
66  {
67  // Get vector of ids and add them to map
68  dynamic_cast<AlgebraicNode*>(father_el_pt->node_pt(j))
69  ->node_update_fct_id(id);
70  unsigned n_id = id.size();
71  for (unsigned i = 0; i < n_id; i++)
72  {
73  id_count[id[i]]++;
74  }
75  }
76 
77  /// Now loop over the ids and check if they appear in all nodes
78  Vector<int> shared_ids;
79  typedef std::map<int, unsigned>::iterator IT;
80  for (IT it = id_count.begin(); it != id_count.end(); it++)
81  {
82  if (it->second == nnod_father)
83  {
84  shared_ids.push_back(it->first);
85  }
86  }
87 
88 
89  // How many update functions we have?
90  unsigned n_update_id = shared_ids.size();
91 
92  // Loop over all node udate functions -- since it's shared by
93  // nodes we may as well read the required data from the first
94  // node in the father.
95  AlgebraicNode* father_node_pt =
96  dynamic_cast<AlgebraicNode*>(father_el_pt->node_pt(0));
97  for (unsigned i = 0; i < n_update_id; i++)
98  {
99  // Actual id:
100  int id = shared_ids[i];
101 
102  // Is this a real update fct or the default dummy one?
103  if (id >= 0)
104  {
105  // Get vector of geometric objects involved in specified node update
106  // function (create vector by copy operation)
107  Vector<GeomObject*> geom_obj_pt(
108  father_node_pt->vector_geom_object_pt(id));
109 
110  // Loop over reference values and obtain the
111  // ones for the current (son) element by interpolation from
112  // the father
113 
114  // Number of reference values for this udpate function
115  unsigned nvalue = father_node_pt->nref_value(id);
116  Vector<double> ref_value(nvalue);
117 
118  // Set up the shape functions in father element
119  Shape psi(nnod_father);
120 
121  // Get shape functions in father element
122  father_el_pt->shape(s_father, psi);
123 
124  // Initialise reference values
125  for (unsigned ivalue = 0; ivalue < nvalue; ivalue++)
126  {
127  ref_value[ivalue] = 0.0;
128  }
129 
130  // Loop over all nodes in father element for
131  // interpolation of nodes -- don't need to
132  // worry about hanging nodes here as reference values
133  // are only assigned once and then in a consistent way
134  for (unsigned j_father = 0; j_father < nnod_father; j_father++)
135  {
136  // Get reference values at node in father by copy operation
137  Vector<double> father_ref_value(
138  dynamic_cast<AlgebraicNode*>(father_el_pt->node_pt(j_father))
139  ->vector_ref_value(id));
140 
141  // Loop over reference values
142  for (unsigned ivalue = 0; ivalue < nvalue; ivalue++)
143  {
144  ref_value[ivalue] += father_ref_value[ivalue] * psi(j_father);
145  }
146  }
147 
148  // Get mesh that implements the update operation
149  AlgebraicMesh* mesh_pt = father_node_pt->mesh_pt(id);
150 
151 
152  // Setup node update info for node
153  alg_node_pt->add_node_update_info(id, // id
154  mesh_pt, // mesh
155  geom_obj_pt, // vector of geom objects
156  ref_value); // vector of ref. values
157 
158  // Update the geometric references (e.g. in FSI) if required
159  mesh_pt->update_node_update(alg_node_pt);
160 
161  } // endif for real vs. default dummy node update fct
162 
163  } // end of loop over different update fcts
164 
165  // Update the node at its current and previous positions
166  bool update_all_time_levels_for_new_node = true;
167  alg_node_pt->node_update(update_all_time_levels_for_new_node);
168  }
169 
170 
171  /// ////////////////////////////////////////////////////////////////////
172  /// ////////////////////////////////////////////////////////////////////
173  // Algebraic nodes
174  /// ////////////////////////////////////////////////////////////////////
175  /// ////////////////////////////////////////////////////////////////////
176 
177 
178  //========================================================================
179  /// Assign default value for test if different
180  /// node update functions produce the same result.
181  //========================================================================
183  1.0e-10;
184 
185 
186  //========================================================================
187  /// Default (negative!) remesh fct id for nodes for which no remesh
188  /// fct is defined
189  //========================================================================
191 
192 
193  //======================================================================
194  /// Set the dummy mesh
195  //====================================================================
197 
198  //========================================================================
199  /// Default dummy mesh to point to for nodes for which no remesh
200  /// fct is defined
201  //========================================================================
203 
204 
205  //========================================================================
206  /// Zero-sized default dummy vector of geom objects to point to for nodes
207  /// for which no remesh fct is defined
208  //========================================================================
210 
211  //========================================================================
212  /// Zero-sized default dummy vector of reference values
213  /// to point to for nodes for which no remesh fct is defined
214  //========================================================================
216 
217 
218  //========================================================================
219  /// Excute the node update function: Update the current (and if
220  /// update_all_time_levels_for_new_node==true also the previous)
221  /// nodal position. Also update the current nodal values if
222  /// an auxiliary update function is defined.
223  /// Note: updating of previous positions is only required (and should
224  /// only be performed for) newly created AlgebraicNodes
225  /// i.e. when this function is called from
226  /// AlgebraicElementBase::setup_algebraic_node_update(...). We create
227  /// the history of its nodal positions from the time-dependent
228  /// version of the specific AlgebraicMesh's algebraic_node_update(...)
229  /// function.
230  //========================================================================
232  const bool& update_all_time_levels_for_new_node)
233  {
234  // Number of time levels that need to be updated
235  unsigned ntime;
236  if (update_all_time_levels_for_new_node)
237  {
238  // Present value plus previous values
239  ntime = 1 + Position_time_stepper_pt->nprev_values();
240  }
241  else
242  {
243  // Present value only
244  ntime = 1;
245  }
246 
247  // Is it a hanging node?
248  if (is_hanging())
249  {
250  // Loop over all master nodes and update their position
251  // That's all we need to update the position of hanging nodes!
252  // (Recall that for hanging nodes Node::x(...) is not
253  // guaranteed to be kept up-to-date; the (constrained!) nodal
254  // position of hanging nodes must be determined via
255  // Node::position() which determines the position
256  // via the hanging node constraints from the position of
257  // the master nodes)
258  unsigned nmaster = hanging_pt()->nmaster();
259  for (unsigned imaster = 0; imaster < nmaster; imaster++)
260  {
261  dynamic_cast<AlgebraicNode*>(hanging_pt()->master_node_pt(imaster))
262  ->node_update();
263  }
264  }
265  // Node isn't hanging --> update it directly
266  else
267  {
268  // If no update function is defined, keep the nodal positions where
269  // they were (i.e. don't do anything), else update
270  if (nnode_update_fcts() != 0)
271  {
272  // Loop over time levels
273  for (unsigned t = 0; t < ntime; t++)
274  {
275  // Update nodal position
276  AlgebraicMesh* mesh_pt = Mesh_pt.begin()->second;
277  // Not sure why I need this intermediate variable....
278  AlgebraicNode* node_pt = this;
279  mesh_pt->algebraic_node_update(t, node_pt);
280  }
281  }
282  }
283 
284  /// Perform auxiliary update of function values?
285  if (Aux_node_update_fct_pt != 0)
286  {
288  }
289  }
290 
291 
292  //========================================================================
293  /// Perform self test: If the node has multiple update functions,
294  /// check that all update functions give the same result (with a tolerance of
295  /// AlgebraicNode::Max_allowed_difference_between_node_update_fcts
296  //========================================================================
298  {
299  // Initialise
300  bool passed = true;
301 
302  unsigned test = Node::self_test();
303  if (test != 0)
304  {
305  passed = false;
306  }
307 
308  // Loop over all update functions
309  unsigned nnode_update = nnode_update_fcts();
310 
311  // If there is just one (or no) update function, then no conflict
312  // can arise and the error is zero.
313  if (nnode_update <= 1)
314  {
315  // return 0;
316  }
317  // Multiple update funcions: check consistency
318  else
319  {
320  // Initialise error
321  double err_max = 0.0;
322 
323  // Spatial (Eulerian) position of the node
324  unsigned ndim_node = ndim();
325  Vector<double> x_0(ndim_node);
326  Vector<double> x_new(ndim_node);
327 
328  // Get vector of update fct ids
329  Vector<int> id;
330  node_update_fct_id(id);
331 
332 
333  // Quick consistency check
334 #ifdef PARANOID
335  if (id.size() != nnode_update)
336  {
337  std::ostringstream error_stream;
338  error_stream << "Inconsistency between number of node update ids:"
339  << nnode_update << " and " << id.size() << std::endl;
340  throw OomphLibError(
341  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
342  }
343 #endif
344 
345 
346  // Update with first update function
347 
348  // Set default update function
350 
351  // Update the node:
352  node_update();
353 
354  // Store coordinates
355  for (unsigned i = 0; i < ndim_node; i++)
356  {
357  x_0[i] = x(i);
358  }
359 
360 
361  // Loop over all other update functions
362  for (unsigned iupdate = 1; iupdate < nnode_update; iupdate++)
363  {
364  // Set default update function
365  set_default_node_update(id[iupdate]);
366 
367  // Update the node:
368  node_update();
369 
370  // Store coordinates
371  for (unsigned i = 0; i < ndim_node; i++)
372  {
373  x_new[i] = x(i);
374  }
375 
376  // Check error
377  double err = 0.0;
378  for (unsigned i = 0; i < ndim_node; i++)
379  {
380  err += (x_new[i] - x_0[i]) * (x_new[i] - x_0[i]);
381  }
382  err = sqrt(err);
383  if (err > err_max)
384  {
385  err_max = err;
386  }
387 
388  // Dump out
390  {
391  oomph_info << "Discrepancy in algebraic update function " << iupdate
392  << ": " << x_0[0] << " " << x_0[1] << " " << x_new[0]
393  << " " << x_new[1] << std::endl;
394 
395  passed = false;
396  }
397  }
398 
399 
400  // Update again with first update function to reset
401 
402  // Set default update function
404 
405  // Update the node:
406  node_update();
407 
408  // Return verdict
409  if (passed)
410  {
411  return 0;
412  }
413  else
414  {
415  return 1;
416  }
417  }
418  // Catch all to remove compiler warning
419  return 0;
420  }
421 
422 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
void setup_algebraic_node_update(Node *&node_pt, const Vector< double > &s_father, FiniteElement *father_el_pt) const
Set up node update info for (newly created) algebraic node: I.e. work out its node update information...
////////////////////////////////////////////////////////////////////
virtual void update_node_update(AlgebraicNode *&node_pt)=0
Update the node update info for given node, following mesh adaptation. Must be implemented for every ...
virtual void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)=0
Update the nodal position posn at time level t (t=0: present; t>0: previous). Must be implemented for...
////////////////////////////////////////////////////////////////////
std::map< int, AlgebraicMesh * > Mesh_pt
Pointer to mesh that performs the specified node update operation (Map because this node may only use...
unsigned nref_value(const int &id)
Number of reference values involved in id-th update function.
static AlgebraicMesh * Dummy_mesh_pt
Default dummy mesh to point to for nodes for which no remesh fct is defined.
static double Max_allowed_difference_between_node_update_fcts
What it says: Used in self-test to check if different node update functions produce the same result.
void set_default_node_update(const int &id)
Make id-th node update function the default.
unsigned self_test()
Perform self test: If the node has multiple node update functions, check that they all give the same ...
static Vector< GeomObject * > Dummy_geom_object_pt
Default dummy vector of geom objects to point to for nodes for which no remesh fct is defined.
void node_update(const bool &update_all_time_levels_for_new_node=false)
Broken assignment operator.
int node_update_fct_id()
Default (usually first if there are multiple ones) node update fct id.
unsigned nnode_update_fcts()
Number of node update fcts.
AlgebraicMesh * mesh_pt()
Default (usually first) mesh that implements update function.
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.
static DummyAlgebraicMesh Dummy_mesh
Static Dummy mesh to which the pointer is addressed.
void add_node_update_info(const int &id, AlgebraicMesh *mesh_pt, const Vector< GeomObject * > &geom_object_pt, const Vector< double > &ref_value, const bool &called_from_constructor=false)
Add algebraic update information for node: What's the ID of the mesh update function (typically used ...
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
static int Dummy_node_update_fct_id
Default (negative!) remesh fct id for nodes for which no remesh fct is defined.
static Vector< double > Dummy_ref_value
Default dummy vector of reference values to point to for nodes for which no remesh fct is defined.
unsigned self_test()
Self-test: Have all values been classified as pinned/unpinned? Return 0 if OK.
Definition: nodes.cc:965
////////////////////////////////////////////////////////////////////
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
AuxNodeUpdateFctPt Aux_node_update_fct_pt
Pointer to auxiliary update function – this can be used to update any nodal values following the upda...
Definition: nodes.h:967
TimeStepper * Position_time_stepper_pt
Pointer to the timestepper associated with the position data.
Definition: nodes.h:932
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An OomphLibError object which should be thrown when an run-time error is encountered....
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...