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