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-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#include "geom_objects.h"
27#include "algebraic_elements.h"
29
30
31namespace 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
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;
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;
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.
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.
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.
static DummyAlgebraicMesh Dummy_mesh
Static Dummy mesh to which the pointer is addressed.
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
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 ...
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
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 *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
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
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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
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...