refineable_line_element.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// Non-inline and non-templated functions for RefineableQElement<1> class
27
28#include <algorithm>
29
30// oomph-lib headers
31#include "mesh.h"
32#include "algebraic_elements.h"
35
36namespace oomph
37{
38 //==========================================================================
39 /// Setup static matrix for coincidence between son nodal points and father
40 /// boundaries:
41 ///
42 /// Father_bound[nnode_1d](nnode_son,son_type) = {L/R/OMEGA}
43 ///
44 /// so that node nnode_son in element of type son_type lies on boundary/
45 /// vertex Father_bound[nnode_1d](nnode_son,son_type) in its father
46 /// element. If the node doesn't lie on a boundary the value is OMEGA.
47 //==========================================================================
49 {
50 using namespace BinaryTreeNames;
51
52 // Find the number of nodes along a 1D edge (which is the number of nodes
53 // in the element for a 1D element!)
54 const unsigned n_node = nnode_1d();
55
56 // Allocate space for the boundary information
57 Father_bound[n_node].resize(n_node, 2);
58
59 // Initialise: By default points are not on the boundary
60 for (unsigned n = 0; n < n_node; n++)
61 {
62 for (unsigned ison = 0; ison < 2; ison++)
63 {
64 Father_bound[n_node](n, ison) = Tree::OMEGA;
65 }
66 }
67
68 // Left-hand son:
69 // --------------
70
71 // L node (0) is the L node of the parent
72 Father_bound[n_node](0, L) = L;
73
74 // Other boundary is in the interior
75
76 // Right-hand son:
77 // ---------------
78
79 // R node (n_node-1) is the R node of the parent
80 Father_bound[n_node](n_node - 1, R) = R;
81
82 // Other boundary is in the interior
83 }
84
85 //==========================================================================
86 /// If a neighbouring element has already created a node at a position
87 /// corresponding to the local fractional position within the present
88 /// element, s_fraction, return a pointer to that node. If not, return
89 /// NULL (0). If the node is periodic the flag is_periodic will be true.
90 //==========================================================================
92 const Vector<double>& s_fraction, bool& is_periodic)
93 {
94 using namespace BinaryTreeNames;
95
96 // Initialise edge of current element on which node lies
97 int edge_in_current = OMEGA;
98
99 // Determine the edge of the current element on which the node lies
100 if (s_fraction[0] == 0.0)
101 {
102 edge_in_current = L;
103 }
104 if (s_fraction[0] == 1.0)
105 {
106 edge_in_current = R;
107 }
108
109 // If the node does not lie on an edge then there is no neighbour:
110 // return NULL
111 if (edge_in_current == OMEGA)
112 {
113 return 0;
114 }
115
116 // Allocate storage for edge in neighbouring element
117 int edge_in_neighbour;
118
119 // Allocate storage for difference in size between current and
120 // neighbouring element
121 int diff_level;
122
123 // Allocate storage for local coordinate of node in neighbouring tree
124 Vector<double> s_in_neighbour(1);
125
126 // Allocate storage for flag indicating if the node is not in the same
127 // binary tree
128 bool in_neighbouring_tree;
129
130 // Allocate storage for the pointer to the neighbouring element
131 // (using its binary tree representation)
132 BinaryTree* neighbour_pt;
133
134 // Find pointer to neighbouring element along the edge in question and
135 // calculate the local coordinate of the node within that element
136 // s_in_neighbour
137 neighbour_pt = binary_tree_pt()->gteq_edge_neighbour(edge_in_current,
138 s_in_neighbour,
139 edge_in_neighbour,
140 diff_level,
141 in_neighbouring_tree);
142
143 // If a neighbour exists...
144 if (neighbour_pt != 0)
145 {
146 // ...check whether its nodes have been created yet
147 if (neighbour_pt->object_pt()->nodes_built())
148 {
149 // If they have, find the node in question in the neighbour
150 Node* neighbour_node_pt =
152 s_in_neighbour);
153
154 // If there is no node at this position, there is a problem, since in
155 // a 1D element (whose nodes have already been built) there should
156 // ALWAYS be a node at each edge of the element.
157 if (neighbour_node_pt == 0)
158 {
159 std::string error_message =
160 "Problem: an element claims to have had its nodes built, yet\n";
161 error_message += "it is missing (a least) a node at its edge.\n";
162 throw OomphLibError(
163 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
164 }
165 // Otherwise, carry on
166 else
167 {
168 // Now work out whether it's a periodic boundary (this is only
169 // (possible if we have moved into a neighbouring tree)
170 if (in_neighbouring_tree)
171 {
172 // Return whether the neighbour is periodic
173 is_periodic = binary_tree_pt()->root_pt()->is_neighbour_periodic(
174 edge_in_current);
175 }
176 // Return the pointer to the neighbouring node
177 return neighbour_node_pt;
178 }
179 }
180 }
181 // Node not found, return null
182 return 0;
183 }
184
185 //==========================================================================
186 /// Build the element by doing the following:
187 /// - Give it nodal positions (by establishing the pointers to its nodes).
188 /// In the process create new nodes where required (i.e. if they don't
189 /// already exist in the father element and have not already been created
190 /// whilst building new neighbouring elements). Node building involves
191 /// the following steps:
192 /// - Get the nodal position from the father element.
193 /// - Establish the time-history of the newly created nodal point
194 /// (its coordinates and previous values) consistent with the father's
195 /// history.
196 /// - Add the new node to the mesh itself.
197 /// - Doc newly created nodes in "new_nodes.dat" stored in the directory
198 /// of the DocInfo object (only if it's open!).
199 /// - NOTE: Unlike in higher-dimensional elements, in 1D it is impossible
200 /// for newly-created nodes to be on a mesh boundary, since any boundary
201 /// nodes must exist in the initial (coarse) mesh. Therefore it is not
202 /// necessary to add any nodes to the mesh's boundary node storage
203 /// schemes, and we always create normal "bulk" nodes.
204 /// - Once the element has a full complement of nodes, excute the element-
205 /// specific further_build() (empty by default -- must be overloaded for
206 /// specific elements). This deals with any build operations that are not
207 /// included in the generic process outlined above. For instance, in
208 /// Crouzeix Raviart elements we need to initialise the internal pressure
209 /// pressure values in a manner consistent with the pressure distribution
210 /// in the father element.
211 //==========================================================================
213 Vector<Node*>& new_node_pt,
214 bool& was_already_built,
215 std::ofstream& new_nodes_file)
216 {
217 using namespace BinaryTreeNames;
218
219 // Find the number of nodes along a 1D edge (which is the number of nodes
220 // in the element for a 1D element!)
221 const unsigned n_node = nnode_1d();
222
223 // Check whether static father_bound needs to be created
224 if (Father_bound[n_node].nrow() == 0)
225 {
226 setup_father_bounds();
227 }
228
229 // Pointer to the current element's father (in binary tree impersonation)
230 BinaryTree* father_pt =
231 dynamic_cast<BinaryTree*>(binary_tree_pt()->father_pt());
232
233 // What type of son is the current element?
234 // Ask its binary tree representation...
235 const int son_type = Tree_pt->son_type();
236
237 // Has somebody built the current element already?
238 // Check this by determining whether or not the first node has been built
239
240 // If this element has not already been built...
241 if (!nodes_built())
242 {
243#ifdef PARANOID
244 if (father_pt == 0)
245 {
246 std::string error_message =
247 "Something fishy here: I have no father and yet \n";
248 error_message += "I have no nodes. Who has created me then?!\n";
249
250 throw OomphLibError(
251 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
252 }
253#endif
254
255 // Indicate status
256 was_already_built = false;
257
258 // Return pointer to father element
259 RefineableQElement<1>* father_el_pt =
260 dynamic_cast<RefineableQElement<1>*>(father_pt->object_pt());
261
262 // Timestepper should be the same for all nodes in father element.
263 // Use it create timesteppers for new nodes.
264 TimeStepper* time_stepper_pt =
265 father_el_pt->node_pt(0)->time_stepper_pt();
266
267 // Determine number of history values (including present)
268 const unsigned ntstorage = time_stepper_pt->ntstorage();
269
270 // Currently we can't handle the case of generalised coordinates
271 // since we haven't established how they should be interpolated.
272 // Buffer this case:
273 if (father_el_pt->node_pt(0)->nposition_type() != 1)
274 {
275 throw OomphLibError("Can't handle generalised nodal positions (yet).",
276 OOMPH_CURRENT_FUNCTION,
277 OOMPH_EXCEPTION_LOCATION);
278 }
279
280 // Set up vertex coordinates in the father element:
281 // ------------------------------------------------
282
283 // Allocate storage for the vertex coordinates s_left and s_right of
284 // the current element as viewed by this element's father (i.e. s_left[0]
285 // stores the local coordinate within the father element at which the
286 // node on the current element's left-hand edge is located. Likewise,
287 // s_right[0] stores the local coordinate within the father element at
288 // which the node on the current element's right-hand edge is located).
289 Vector<double> s_left(1);
290 Vector<double> s_right(1);
291
292 // In order to set up the vertex coordinates we need to know which
293 // type of son the current element is
294 switch (son_type)
295 {
296 case L:
297 s_left[0] = -1.0;
298 s_right[0] = 0.0;
299 break;
300
301 case R:
302 s_left[0] = 0.0;
303 s_right[0] = 1.0;
304 break;
305 }
306
307 // Pass macro element pointer on to sons and set coordinates in
308 // macro element
309 // hierher why can I see this?
310 if (father_el_pt->Macro_elem_pt != 0)
311 {
312 set_macro_elem_pt(father_el_pt->Macro_elem_pt);
313
314 s_macro_ll(0) =
315 father_el_pt->s_macro_ll(0) +
316 0.5 * (s_left[0] + 1.0) *
317 (father_el_pt->s_macro_ur(0) - father_el_pt->s_macro_ll(0));
318 s_macro_ur(0) =
319 father_el_pt->s_macro_ll(0) +
320 0.5 * (s_right[0] + 1.0) *
321 (father_el_pt->s_macro_ur(0) - father_el_pt->s_macro_ll(0));
322 }
323
324 // If the father element hasn't been generated yet, we're stuck...
325 if (father_el_pt->node_pt(0) == 0)
326 {
327 throw OomphLibError(
328 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
329 OOMPH_CURRENT_FUNCTION,
330 OOMPH_EXCEPTION_LOCATION);
331 }
332 // Otherwise, carry on
333 else
334 {
335 // Allocate storage for the location of a node in terms of the local
336 // coordinate of the father element
337 Vector<double> s_in_father(1);
338
339 // Allocate storage for the fractional position (in the current
340 // element) of the node in the direction of s[0]
341 Vector<double> s_fraction(1);
342
343 // Loop over all nodes in the element
344 for (unsigned n = 0; n < n_node; n++)
345 {
346 // Get the fractional position (in the current element) of the node
347 // in the direction of s[0]
348 s_fraction[0] = local_one_d_fraction_of_node(n, 0);
349
350 // Evaluate the local coordinate of the node in the father element
351 s_in_father[0] = s_left[0] + (s_right[0] - s_left[0]) * s_fraction[0];
352
353 // Initialise flag: So far, this node hasn't been built or copied yet
354 bool node_done = false;
355
356 // Get the pointer to the node in the father (returns NULL if there
357 // is no node at this position)
358 Node* created_node_pt =
359 father_el_pt->get_node_at_local_coordinate(s_in_father);
360
361 // Does this node already exist in father element?
362 // -----------------------------------------------
363
364 // If it does...
365 if (created_node_pt != 0)
366 {
367 // ...copy node across
368 node_pt(n) = created_node_pt;
369
370 // Make sure that we update the values at the node so that they
371 // are consistent with the present representation. This is only
372 // needed for mixed interpolation where the value at the father
373 // could now become active.
374
375 // Loop over all history values
376 for (unsigned t = 0; t < ntstorage; t++)
377 {
378 // Get values from father element
379 // Note: get_interpolated_values() sets Vector size itself
380 Vector<double> prev_values;
381 father_el_pt->get_interpolated_values(
382 t, s_in_father, prev_values);
383
384 // Find the minimum number of values (either those stored at the
385 // node, or those returned by the function)
386 unsigned n_val_at_node = created_node_pt->nvalue();
387 unsigned n_val_from_function = prev_values.size();
388
389 // Use the ternary conditional operator here
390 unsigned n_var = n_val_at_node < n_val_from_function ?
391 n_val_at_node :
392 n_val_from_function;
393
394 // Assign the values that we can
395 for (unsigned k = 0; k < n_var; k++)
396 {
397 created_node_pt->set_value(t, k, prev_values[k]);
398 }
399 }
400
401 // Indicate that node has been created by copy
402 node_done = true;
403 }
404
405 // Node does not exist in father element but might already
406 // -------------------------------------------------------
407 // have been created by neighbouring elements
408 // ------------------------------------------
409
410 else
411 {
412 // Was the node created by one of its neighbours? Whether or not
413 // the node lies on an edge can be calculated from the fractional
414 // position.
415 bool is_periodic = false;
416 created_node_pt =
417 node_created_by_neighbour(s_fraction, is_periodic);
418
419 // If the node was so created...
420 if (created_node_pt != 0)
421 {
422 // ...assign the pointer
423 node_pt(n) = created_node_pt;
424
425 // Indicate that node has been created by copy
426 node_done = true;
427
428 // In a 1D mesh there is no way that a periodic node (which must
429 // be on a boundary) can exist without being part of the initial
430 // (coarse) mesh. Therefore issue an error message if
431 // node_created_by_neighbour(...) returns `is_periodic==true'.
432#ifdef PARANOID
433 if (is_periodic)
434 {
435 std::string error_message =
436 "node_created_by_neighbour returns a node which claims\n";
437 error_message += "to be periodic. In a 1D mesh any periodic "
438 "nodes must exist\n";
439 error_message += "in the initial (coarse) mesh.";
440
441 throw OomphLibError(error_message,
442 OOMPH_CURRENT_FUNCTION,
443 OOMPH_EXCEPTION_LOCATION);
444 }
445#endif
446 }
447 }
448
449 // Node does not exist in father element or neighbouring element
450 // -------------------------------------------------------------
451
452 // If the node has not been built anywhere ---> build it here
453 if (!node_done)
454 {
455 // In a 1D mesh any node which lies on the boundary must exist in
456 // the initial (coarse) mesh, so any newly-built nodes cannot be
457 // boundary nodes. Therefore we always create a normal "bulk" node.
458
459 // Create node and set the pointer to it from the element
460 created_node_pt = construct_node(n, time_stepper_pt);
461
462 // Add to vector of new nodes
463 new_node_pt.push_back(created_node_pt);
464
465 // Now we set the position and values at the newly created node.
466 // In the first instance use macro element or FE representation
467 // to create past and present nodal positions.
468 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC ELEMENTS AS NOT
469 // ALL OF THEM NECESSARILY IMPLEMENT NONTRIVIAL NODE UPDATE
470 // FUNCTIONS. CALLING THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL
471 // LEAVE THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
472 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL NOT ASSIGN SENSIBLE
473 // INITIAL POSITIONS!)
474
475 // Loop over history values
476 for (unsigned t = 0; t < ntstorage; t++)
477 {
478 // Allocate storage for the previous position of the node
479 Vector<double> x_prev(1);
480
481 // Get position from father element -- this uses the macro
482 // element representation if appropriate.
483 father_el_pt->get_x(t, s_in_father, x_prev);
484
485 // Set the previous position of the new node
486 created_node_pt->x(t, 0) = x_prev[0];
487
488 // Allocate storage for the previous values at the node
489 // NOTE: the size of this vector is equal to the number of values
490 // (e.g. 3 velocity components and 1 pressure, say)
491 // associated with each node and NOT the number of history values
492 // which the node stores!
493 Vector<double> prev_values;
494
495 // Get values from father element
496 // Note: get_interpolated_values() sets Vector size itself.
497 father_el_pt->get_interpolated_values(
498 t, s_in_father, prev_values);
499
500 // Determine the number of values at the new node
501 const unsigned n_value = created_node_pt->nvalue();
502
503 // Loop over all values and set the previous values
504 for (unsigned v = 0; v < n_value; v++)
505 {
506 created_node_pt->set_value(t, v, prev_values[v]);
507 }
508 } // End of loop over history values
509
510 // Add new node to mesh
511 mesh_pt->add_node_pt(created_node_pt);
512
513 } // End of case where we build the node ourselves
514
515 // Check if the element is an algebraic element
516 AlgebraicElementBase* alg_el_pt =
517 dynamic_cast<AlgebraicElementBase*>(this);
518
519 // If it is, set up node position (past and present) from algebraic
520 // node position (past and present) from algebraic node update
521 // node update function. This over-writes previous assingments that
522 // were made based on the macro-element/FE representation.
523 // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE NODE IS A
524 // MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN THE SAME NODAL
525 // POSITIONS BUT WE NEED TO ADD THE REMESH INFO FOR *ALL* ROOT
526 // ELEMENTS!
527 if (alg_el_pt != 0)
528 {
529 // Build algebraic node update info for new node. This sets up
530 // the node update data for all node update functions that are
531 // shared by all nodes in the father element.
533 node_pt(n), s_in_father, father_el_pt);
534 }
535
536 // If we have built the node and we are documenting our progress,
537 // write the (hopefully consistent position) to the outputfile
538 if ((!node_done) && (new_nodes_file.is_open()))
539 {
540 new_nodes_file << node_pt(n)->x(0) << std::endl;
541 }
542 } // End of loop over all nodes in element
543
544
545 // If the element is a MacroElementNodeUpdateElement, set the update
546 // parameters for the current element's nodes -- all this needs is
547 // the vector of (pointers to the) geometric objects that affect the
548 // MacroElement-based node update. This is the same as that in the
549 // father element
550 MacroElementNodeUpdateElementBase* father_m_el_pt =
551 dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
552 if (father_m_el_pt != 0)
553 {
554 // Get Vector of geometric objects from father (construct Vector
555 // via copy operation)
556 Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
557
558 // Cast current element to MacroElementNodeUpdateElement:
560 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
561
562#ifdef PARANOID
563 if (m_el_pt == 0)
564 {
565 std::string error_message =
566 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
567 error_message +=
568 "Strange -- if the father is a MacroElementNodeUpdateElement\n";
569 error_message += "the son should be too....\n";
570
571 throw OomphLibError(
572 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
573 }
574#endif
575 // Build update info by passing Vector of geometric objects:
576 // This sets the current element to be the update element for all
577 // of the element's nodes. This is reversed if the element is
578 // ever un-refined in the father element's rebuild_from_sons()
579 // function which overwrites this assignment to avoid nasty
580 // segmentation faults that occur when a node tries to update
581 // itself via an element that no longer exists...
582 m_el_pt->set_node_update_info(geom_object_pt);
583 }
584
585#ifdef OOMPH_HAS_MPI
586 // Pass on non-halo proc ID
587 Non_halo_proc_ID =
588 tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
589#endif
590
591 // Is the new element an ElementWithMovingNodes?
592 ElementWithMovingNodes* aux_el_pt =
593 dynamic_cast<ElementWithMovingNodes*>(this);
594
595 // Pass down the information re the method for the evaluation
596 // of the shape derivatives
597 if (aux_el_pt != 0)
598 {
599 ElementWithMovingNodes* aux_father_el_pt =
600 dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
601
602#ifdef PARANOID
603 if (aux_father_el_pt == 0)
604 {
605 std::string error_message =
606 "Failed to cast to ElementWithMovingNodes*\n";
607 error_message +=
608 "Strange -- if the son is a ElementWithMovingNodes\n";
609 error_message += "the father should be too....\n";
610
611 throw OomphLibError(
612 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
613 }
614#endif
615
616 // If evaluating the residuals by finite differences in the father
617 // continue to do so in the child
618 if (aux_father_el_pt
619 ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
620 {
621 aux_el_pt
623 }
624
625 aux_el_pt->method_for_shape_derivs() =
626 aux_father_el_pt->method_for_shape_derivs();
627
628 // If bypassing the evaluation of fill_in_jacobian_from_geometric_data
629 // continue to do so
630 if (aux_father_el_pt
631 ->is_fill_in_jacobian_from_geometric_data_bypassed())
632 {
634 }
635 }
636
637
638 // Now do further build (if any)
639 further_build();
640
641 } // Sanity check: Father element has been generated
642
643 } // End of if element has not already been built
644
645 // If element has already been built, say so
646 else
647 {
648 was_already_built = true;
649 }
650 }
651
652 //==========================================================================
653 /// Print corner nodes, use colour (default "BLACK")
654 //==========================================================================
655 void RefineableQElement<1>::output_corners(std::ostream& outfile,
656 const std::string& colour) const
657 {
658 // Allocate storage for local coordinate s
659 Vector<double> s(1);
660
661 // Allocate storage for global coordinate of element vertex
662 Vector<double> vertex(1);
663
664 outfile << "ZONE I=2,J=2, C=" << colour << std::endl;
665
666 // Left-hand vertex
667 s[0] = -1.0;
668 get_x(s, vertex);
669 outfile << vertex[0] << " " << Number << std::endl;
670
671 // Right-hand vertex
672 s[0] = 1.0;
673 get_x(s, vertex);
674 outfile << vertex[0] << " " << Number << std::endl;
675
676 outfile << "TEXT CS = GRID, X = " << vertex[0]
677 << ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\"" << Number << "\""
678 << std::endl;
679 }
680
681 //==========================================================================
682 /// Check inter-element continuity of
683 /// - nodal positions
684 /// - (nodally) interpolated function values
685 //==========================================================================
687 {
688 using namespace BinaryTreeNames;
689
690 // Calculate number of nodes
691 const unsigned n_node = nnode_1d();
692
693 // Number of timesteps (including present) for which continuity is
694 // to be checked
695 const unsigned n_time = 1;
696
697 // Initialise errors
698 max_error = 0.0;
699 Vector<double> max_error_x(1, 0.0);
700 double max_error_val = 0.0;
701
702 // Initialise vector of element edges in the current element
703 Vector<int> edge_in_current(2);
704 edge_in_current[0] = L;
705 edge_in_current[1] = R;
706
707 // Loop over both edges
708 for (unsigned edge_counter = 0; edge_counter < 2; edge_counter++)
709 {
710 // Allocate storage for the fractional position (in the current
711 // element) of the node in the direction of s[0]
712 Vector<double> s_fraction(1);
713
714 // Allocate storage for the location of a node in terms of the local
715 // coordinate of the current element
716 Vector<double> s_in_current(1);
717
718 // Allocate storage for the location of a node in terms of the local
719 // coordinate of the neighbouring element
720 Vector<double> s_in_neighbour(1);
721
722 // Allocate storage for edge in neighbouring element
723 int edge_in_neighbour;
724
725 // Allocate storage for difference in size between current and
726 // neighbouring element
727 int diff_level;
728
729 // Allocate storage for flag indicating if the node is not in the same
730 // binary tree
731 bool in_neighbouring_tree;
732
733 // Calculate the local coordinate and the local coordinate as viewed
734 // from the neighbour
735
736 // Allocate storage for the pointer to the neighbouring element
737 // (using its binary tree representation)
738 BinaryTree* neighbour_pt;
739
740 // Find pointer to neighbouring element along the edge in question and
741 // calculate the local coordinate of the node within that element,
742 // s_in_neighbour
743 neighbour_pt =
744 binary_tree_pt()->gteq_edge_neighbour(edge_in_current[edge_counter],
745 s_in_neighbour,
746 edge_in_neighbour,
747 diff_level,
748 in_neighbouring_tree);
749
750 // If neighbour exists and has existing nodes
751 if ((neighbour_pt != 0) && (neighbour_pt->object_pt()->nodes_built()))
752 {
753 // Need to exclude periodic nodes from this check.
754 // There are only periodic nodes if we are in a neighbouring tree...
755 bool is_periodic = false;
756 if (in_neighbouring_tree)
757 {
758 // Is it periodic?
759 is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(
760 edge_in_current[edge_counter]);
761 }
762
763 // Allocate storage for pointer to the local node
764 Node* local_node_pt = 0;
765
766 switch (edge_counter)
767 {
768 case 0:
769 // Local fraction of node
770 s_fraction[0] = 0.0;
771 // Get pointer to local node
772 local_node_pt = node_pt(0);
773 break;
774
775 case 1:
776 // Local fraction of node
777 s_fraction[0] = 1.0;
778 // Get pointer to local node
779 local_node_pt = node_pt(n_node - 1);
780 break;
781 }
782
783 // Evaluate the local coordinate of the node in the current element
784 s_in_current[0] = -1.0 + 2.0 * s_fraction[0];
785
786 // NOTE: We have already calculated the local coordinate of the node
787 // in the neighbouring element above when calling gteq_edge_neighbour()
788
789 // Loop over timesteps
790 for (unsigned t = 0; t < n_time; t++)
791 {
792 // Allocate storage for the nodal position in the neighbouring element
793 Vector<double> x_in_neighbour(1);
794
795 // Get the nodal position from the neighbouring element
796 neighbour_pt->object_pt()->interpolated_x(
797 t, s_in_neighbour, x_in_neighbour);
798
799 // Check error only if the node is NOT periodic
800 if (is_periodic == false)
801 {
802 // Find the spatial error
803 double err = std::fabs(local_node_pt->x(t, 0) - x_in_neighbour[0]);
804
805 // If it's bigger than our tolerance, say so
806 if (err > 1e-9)
807 {
808 oomph_info << "errx " << err << " " << t << " "
809 << local_node_pt->x(t, 0) << " " << x_in_neighbour[0]
810 << std::endl;
811
812 oomph_info << "at " << local_node_pt->x(0) << std::endl;
813 }
814
815 // If it's bigger than the previous max error, it is the
816 // new max error!
817 if (err > max_error_x[0])
818 {
819 max_error_x[0] = err;
820 }
821 }
822 // Allocate storage for the nodal values in the neighbouring element
823 Vector<double> values_in_neighbour;
824
825 // Get the values from neighbouring element. NOTE: Number of values
826 // gets set by routine (because in general we don't know how many
827 // interpolated values a certain element has.
828 neighbour_pt->object_pt()->get_interpolated_values(
829 t, s_in_neighbour, values_in_neighbour);
830
831 // Allocate storage for the nodal values in the current element
832 Vector<double> values_in_current;
833
834 // Get the values in current element
835 get_interpolated_values(t, s_in_current, values_in_current);
836
837 // Now figure out how many continuously interpolated values there are
838 const unsigned num_val =
839 neighbour_pt->object_pt()->ncont_interpolated_values();
840
841 // Check error
842 for (unsigned ival = 0; ival < num_val; ival++)
843 {
844 // Find the spatial error
845 double err =
846 std::fabs(values_in_current[ival] - values_in_neighbour[ival]);
847
848 // If it's bigger than our tolerance, say so
849 if (err > 1.0e-10)
850 {
851 oomph_info << local_node_pt->x(0) << "\n# "
852 << "erru " << err << " " << ival << " "
853 << get_node_number(local_node_pt) << " "
854 << values_in_current[ival] << " "
855 << values_in_neighbour[ival] << std::endl;
856 }
857
858 // If it's bigger than the previous max error, it is the
859 // new max error!
860 if (err > max_error_val)
861 {
862 max_error_val = err;
863 }
864 }
865 } // End of loop over timesteps
866 } // End of if neighbour exists and has existing nodes
867 } // End of loop over edges
868
869 // Update max_error if necessary
870 max_error = max_error_x[0];
871 if (max_error_val > max_error)
872 {
873 max_error = max_error_val;
874 }
875
876 // Output max_error information
877 if (max_error > 1e-9)
878 {
879 oomph_info << "\n#------------------------------------ \n#Max error ";
880 oomph_info << max_error_x[0] << " " << max_error_val << std::endl;
881 oomph_info << "#------------------------------------ \n " << std::endl;
882 }
883 }
884
885 //==========================================================================
886 /// Static matrix for coincidence between son nodal points and father
887 /// boundaries
888 //==========================================================================
889 std::map<unsigned, DenseMatrix<int>> RefineableQElement<1>::Father_bound;
890
891} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
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...
BinaryTree class: Recursively defined, generalised binary tree.
Definition: binary_tree.h:92
BinaryTree * gteq_edge_neighbour(const int &direction, Vector< double > &s_in_neighbour, int &edge, int &diff_level, bool &in_neighbouring_tree) const
Return pointer to greater or equal-sized edge neighbour in specified direction; also provide info reg...
Definition: binary_tree.cc:172
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
void enable_bypass_fill_in_jacobian_from_geometric_data()
Bypass the call to fill_in_jacobian_from_geometric_data.
int & method_for_shape_derivs()
Access to method (enumerated flag) for determination of shape derivs.
void enable_always_evaluate_dresidual_dnodal_coordinates_by_fd()
Insist that shape derivatives are always evaluated by fd (using FiniteElement::get_dresidual_dnodal_c...
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
MacroElement * Macro_elem_pt
Pointer to the element's macro element (NULL by default)
Definition: elements.h:1683
void get_x(const Vector< double > &s, Vector< double > &x) const
Global coordinates as function of local coordinates. Either via FE representation or via macro-elemen...
Definition: elements.h:1885
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
////////////////////////////////////////////////////////////////////
virtual void set_node_update_info(const Vector< GeomObject * > &geom_object_pt)=0
Set node update information: Pass the vector of (pointers to) the geometric objects that affect the n...
Vector< GeomObject * > & geom_object_pt()
Vector of (pointers to) geometric objects involved in node update function.
A general mesh class.
Definition: mesh.h:67
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
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
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
An OomphLibError object which should be thrown when an run-time error is encountered....
double & s_macro_ur(const unsigned &i)
Access fct to the i-th coordinate of the element's "upper right" vertex in the associated MacroElemen...
Definition: Qelements.h:202
double & s_macro_ll(const unsigned &i)
Access fct to the i-th coordinate of the element's "lower left" vertex in the associated MacroElement...
Definition: Qelements.h:186
virtual bool nodes_built()
Return true if all the nodes have been built, false if not.
virtual unsigned ncont_interpolated_values() const =0
Number of continuously interpolated values. Note: We assume that they are located at the beginning of...
virtual void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get all continously interpolated function values in this element as a Vector. Note: Vector sets is ow...
Refineable version of QElement<1,NNODE_1D>.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
bool is_neighbour_periodic(const int &direction)
Return whether the neighbour in the particular direction is periodic.
Definition: tree.h:364
TreeRoot *& root_pt()
Return pointer to root of the tree.
Definition: tree.h:141
int son_type() const
Return son type.
Definition: tree.h:214
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:262
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Number
The unsigned.
Vector< std::string > colour
Tecplot colours.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...