refineable_quad_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#include <algorithm>
27
28#include "mesh.h"
29#include "algebraic_elements.h"
32
33namespace oomph
34{
35 //==================================================================
36 /// Setup static matrix for coincidence between son nodal points and
37 /// father boundaries:
38 ///
39 /// Father_boundd[nnode_1d](nnode_son,son_type)={SW/SE/NW/NE/S/E/N/W/OMEGA}
40 ///
41 /// so that node nnode_son in element of type son_type lies
42 /// on boundary/vertex Father_boundd[nnode_1d](nnode_son,son_type) in its
43 /// father element. If the node doesn't lie on a boundary
44 /// the value is OMEGA.
45 //==================================================================
47 {
48 using namespace QuadTreeNames;
49
50 // Find the number of nodes along a 1D edge
51 unsigned n_p = nnode_1d();
52 // Allocate space for the boundary information
53 Father_bound[n_p].resize(n_p * n_p, 4);
54
55 // Initialise: By default points are not on the boundary
56 for (unsigned n = 0; n < n_p * n_p; n++)
57 {
58 for (unsigned ison = 0; ison < 4; ison++)
59 {
60 Father_bound[n_p](n, ison) = Tree::OMEGA;
61 }
62 }
63
64 // Southwest son
65 //--------------
66 // SW node (0) is the SW node of the parent
67 Father_bound[n_p](0, SW) = SW;
68 // Southern boundary is the southern boundary of the parent
69 for (unsigned n = 1; n < n_p; n++)
70 {
71 Father_bound[n_p](n, SW) = S;
72 }
73 // Western boundary is the western boundary of the parent
74 for (unsigned n = 1; n < n_p; n++)
75 {
76 Father_bound[n_p](n_p * n, SW) = W;
77 }
78 // Other boundaries are in the interior
79
80 // Northwest son
81 //--------------
82 // NW node (n_p*(n_p-1))is the NW node of the parent
83 Father_bound[n_p](n_p * (n_p - 1), NW) = NW;
84 // Northern boundary is the northern boundary of the parent
85 for (unsigned n = 1; n < n_p; n++)
86 {
87 Father_bound[n_p](n_p * (n_p - 1) + n, NW) = N;
88 }
89 // Western boundary is the western boundary of the parent
90 for (unsigned n = 0; n < (n_p - 1); n++)
91 {
92 Father_bound[n_p](n_p * n, NW) = W;
93 }
94 // Other boundaries are in the interior
95
96 // Northeast son
97 //--------------
98 // NE node (n_p*n_p-1) is the NE node of the parent
99 Father_bound[n_p](n_p * n_p - 1, NE) = NE;
100 // Northern boundary is the northern boundary of the parent
101 for (unsigned n = 0; n < (n_p - 1); n++)
102 {
103 Father_bound[n_p](n_p * (n_p - 1) + n, NE) = N;
104 }
105 // Eastern boundary is the eastern boundary of the parent
106 for (unsigned n = 0; n < (n_p - 1); n++)
107 {
108 Father_bound[n_p](n_p - 1 + n_p * n, NE) = E;
109 }
110 // Other boundaries are in the interior
111
112 // Southeast son
113 //--------------
114 // SE node (n_p-1) is the SE node of the parent
115 Father_bound[n_p](n_p - 1, SE) = SE;
116 // Southern boundary is the southern boundary of the parent
117 for (unsigned n = 0; n < (n_p - 1); n++)
118 {
119 Father_bound[n_p](n, SE) = S;
120 }
121 // Eastern boundary is the eastern boundary of the parent
122 for (unsigned n = 1; n < n_p; n++)
123 {
124 Father_bound[n_p](n_p - 1 + n_p * n, SE) = E;
125 }
126 }
127
128 //==================================================================
129 /// Determine Vector of boundary conditions along the element's boundary
130 /// (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
131 ///
132 /// This function assumes that the same boundary condition is applied
133 /// along the entire length of an element's edge (of course, the
134 /// vertices combine the boundary conditions of their two adjacent edges
135 /// in the most restrictive combination. Hence, if we're at a vertex,
136 /// we apply the most restrictive boundary condition of the
137 /// two adjacent edges. If we're on an edge (in its proper interior),
138 /// we apply the least restrictive boundary condition of all nodes
139 /// along the edge.
140 ///
141 /// Usual convention:
142 /// - bound_cons[ival]=0 if value ival on this boundary is free
143 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
144 //==================================================================
145 void RefineableQElement<2>::get_bcs(int bound, Vector<int>& bound_cons) const
146 {
147 using namespace QuadTreeNames;
148
149 // Max. number of nodal data values in element
150 unsigned nvalue = ncont_interpolated_values();
151 // Set up temporary vectors to hold edge boundary conditions
152 Vector<int> bound_cons1(nvalue), bound_cons2(nvalue);
153
154 // Which boundary are we on?
155 switch (bound)
156 {
157 // If on edge, just get the bcs
158 case N:
159 case S:
160 case W:
161 case E:
162 get_edge_bcs(bound, bound_cons);
163 break;
164
165 // Most restrictive boundary at SE corner
166 case SE:
167 get_edge_bcs(S, bound_cons1);
168 get_edge_bcs(E, bound_cons2);
169
170 for (unsigned k = 0; k < nvalue; k++)
171 {
172 bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
173 }
174 break;
175
176 // Most restrictive boundary at SW corner
177 case SW:
178 get_edge_bcs(S, bound_cons1);
179 get_edge_bcs(W, bound_cons2);
180
181 for (unsigned k = 0; k < nvalue; k++)
182 {
183 bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
184 }
185 break;
186
187 // Most restrictive boundary at NW corner
188 case NW:
189 get_edge_bcs(N, bound_cons1);
190 get_edge_bcs(W, bound_cons2);
191
192 for (unsigned k = 0; k < nvalue; k++)
193 {
194 bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
195 }
196 break;
197
198 // Most restrictive boundary at NE corner
199 case NE:
200 get_edge_bcs(N, bound_cons1);
201 get_edge_bcs(E, bound_cons2);
202
203 for (unsigned k = 0; k < nvalue; k++)
204 {
205 bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
206 }
207 break;
208
209 default:
210 throw OomphLibError(
211 "Wrong boundary", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
212 }
213 }
214
215 //==================================================================
216 /// Determine Vector of boundary conditions along the element's
217 /// edge (S/N/W/E) -- BC is the least restrictive combination
218 /// of all the nodes on this edge
219 ///
220 /// Usual convention:
221 /// - bound_cons[ival]=0 if value ival on this boundary is free
222 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
223 //==================================================================
225 Vector<int>& bound_cons) const
226 {
227 using namespace QuadTreeNames;
228
229 // Number of nodes along 1D edge
230 unsigned n_p = nnode_1d();
231 // Left- and Right-hand nodes
232 unsigned left_node, right_node;
233
234 // Set the left (lower) and right (upper) nodes for the edge
235 switch (edge)
236 {
237 case N:
238 left_node = n_p * (n_p - 1);
239 right_node = n_p * n_p - 1;
240 break;
241
242 case S:
243 left_node = 0;
244 right_node = n_p - 1;
245 break;
246
247 case W:
248 left_node = 0;
249 right_node = n_p * (n_p - 1);
250 break;
251
252 case E:
253 left_node = n_p - 1;
254 right_node = n_p * n_p - 1;
255 break;
256
257 default:
258 std::ostringstream error_stream;
259 error_stream << "Wrong edge " << edge << " passed to get_edge_bcs(..)"
260 << std::endl;
261
262 throw OomphLibError(
263 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
264 }
265
266 // Max. number of nodal data values in element
267 unsigned maxnvalue = ncont_interpolated_values();
268
269 // Loop over all values, the least restrictive value is
270 // the boundary condition at the left multiplied by that at the right
271 // Assuming that free is always zero and pinned is one
272 for (unsigned k = 0; k < maxnvalue; k++)
273 {
274 bound_cons[k] =
275 node_pt(left_node)->is_pinned(k) * node_pt(right_node)->is_pinned(k);
276 }
277 }
278
279 //==================================================================
280 /// Given an element edge/vertex, return a set that contains
281 /// all the (mesh-)boundary numbers that this element edge/vertex
282 /// lives on.
283 ///
284 /// For proper edges, the boundary is the one (if any) that is shared by
285 /// both vertex nodes). For vertex nodes, we just return their
286 /// boundaries.
287 //==================================================================
289 std::set<unsigned>& boundary) const
290 {
291 using namespace QuadTreeNames;
292
293 // Number of 1d nodes along an edge
294 unsigned n_p = nnode_1d();
295 // Left and right-hand nodes
296 int left_node = -1, right_node = -1;
297
298 // Set the left (lower) and right (upper) nodes for the edge
299 switch (edge)
300 {
301 case N:
302 left_node = n_p * (n_p - 1);
303 right_node = n_p * n_p - 1;
304 break;
305
306 case S:
307 left_node = 0;
308 right_node = n_p - 1;
309 break;
310
311 case W:
312 left_node = 0;
313 right_node = n_p * (n_p - 1);
314 break;
315
316 case E:
317 left_node = n_p - 1;
318 right_node = n_p * n_p - 1;
319 break;
320
321 // Vertices do not have left nodes!
322 case SE:
323 right_node = n_p - 1;
324 break;
325
326 case SW:
327 right_node = 0;
328 break;
329
330 case NE:
331 right_node = n_p * n_p - 1;
332 break;
333
334 case NW:
335 right_node = n_p * (n_p - 1);
336 break;
337
338 default:
339 std::ostringstream error_stream;
340 error_stream << "Wrong edge " << edge << " passed" << std::endl;
341
342 throw OomphLibError(
343 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
344 }
345
346 // Empty the boundary set: Edge does not live on any boundary
347 boundary.clear();
348
349 // Storage for the nodes that the right node lives on
350 std::set<unsigned>* right_boundaries_pt = 0;
351 // Get the boundaries that the right node lives on
352 node_pt(right_node)->get_boundaries_pt(right_boundaries_pt);
353
354 // If the right node lives on some boundaries
355 if (right_boundaries_pt != 0)
356 {
357 // If the node is a vertex then add the boundaries at the node
358 // into the vector boundary
359 if (left_node < 0)
360 {
361 copy(right_boundaries_pt->begin(),
362 right_boundaries_pt->end(),
363 inserter(boundary, boundary.begin()));
364 }
365 // Otherwise only add if the boundary also exists at the left hand node
366 else
367 {
368 std::set<unsigned>* left_boundaries_pt = 0;
369 node_pt(left_node)->get_boundaries_pt(left_boundaries_pt);
370 // If the left node is on some boundaries
371 if (left_boundaries_pt != 0)
372 {
373 // Use the standard algorithms to compute the boundaries in
374 // common between the left and right nodes
375 std::set_intersection(right_boundaries_pt->begin(),
376 right_boundaries_pt->end(),
377 left_boundaries_pt->begin(),
378 left_boundaries_pt->end(),
379 inserter(boundary, boundary.begin()));
380 }
381 }
382 }
383 }
384
385
386 //===================================================================
387 /// Return the value of the intrinsic boundary coordinate interpolated
388 /// along the edge (S/W/N/E)
389 //===================================================================
391 const unsigned& boundary,
392 const int& edge,
393 const Vector<double>& s,
394 Vector<double>& zeta)
395 {
396 using namespace QuadTreeNames;
397
398 // Number of 1D nodes along an edge
399 unsigned n_p = nnode_1d();
400
401 // Storage for the shape functions
402 Shape psi(n_p * n_p);
403 // Get the shape functions at the passed position
404 this->shape(s, psi);
405
406 // Unsigned data that give starts and multipliers for the loop
407 // over the nodes on the edges.
408 unsigned start = 0, multiplier = 1;
409
410 // Which edge?
411 switch (edge)
412 {
413 case S:
414#ifdef PARANOID
415 if (s[1] != -1.0)
416 {
417 std::ostringstream error_stream;
418 error_stream << "Coordinate " << s[0] << " " << s[1]
419 << " is not on South edge\n";
420
421 throw OomphLibError(error_stream.str(),
422 OOMPH_CURRENT_FUNCTION,
423 OOMPH_EXCEPTION_LOCATION);
424 }
425#endif
426 // Start is zero and multiplier is one
427 break;
428
429 case N:
430#ifdef PARANOID
431 if (s[1] != 1.0)
432 {
433 std::ostringstream error_stream;
434 error_stream << "Coordinate " << s[0] << " " << s[1]
435 << " is not on North edge\n";
436
437 throw OomphLibError(error_stream.str(),
438 OOMPH_CURRENT_FUNCTION,
439 OOMPH_EXCEPTION_LOCATION);
440 }
441#endif
442 // Start from the top left corner of the element, multiplier still one
443 start = n_p * (n_p - 1);
444 break;
445
446 case W:
447#ifdef PARANOID
448 if (s[0] != -1.0)
449 {
450 std::ostringstream error_stream;
451 error_stream << "Coordinate " << s[0] << " " << s[1]
452 << " is not on West edge\n";
453
454 throw OomphLibError(error_stream.str(),
455 OOMPH_CURRENT_FUNCTION,
456 OOMPH_EXCEPTION_LOCATION);
457 }
458#endif
459 // Loop over left-hand edge of element (start from zero)
460 multiplier = n_p;
461 break;
462
463 case E:
464#ifdef PARANOID
465 if (s[0] != 1.0)
466 {
467 std::ostringstream error_stream;
468 error_stream << "Coordinate " << s[0] << " " << s[1]
469 << " is not on East edge\n";
470
471 throw OomphLibError(error_stream.str(),
472 OOMPH_CURRENT_FUNCTION,
473 OOMPH_EXCEPTION_LOCATION);
474 }
475#endif
476 // Start from the bottom right-hand corner
477 start = n_p - 1;
478 // Loop over the right-hand edge of the element
479 multiplier = n_p;
480 break;
481
482
483 default:
484 std::ostringstream error_stream;
485 error_stream << "Wrong edge " << edge << " passed" << std::endl;
486
487 throw OomphLibError(
488 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
489 }
490
491 // Initialise the intrinsic coordinate
492 double inter_zeta = 0.0;
493 // Loop over the nodes on the edge
494 for (unsigned n = 0; n < n_p; n++)
495 {
496 // Get the node number
497 unsigned node_number = start + multiplier * n;
498 // Now get the intrinsic coordinate
499 node_pt(node_number)->get_coordinates_on_boundary(boundary, zeta);
500 // Now multiply by the shape function
501 inter_zeta += zeta[0] * psi(node_number);
502 }
503
504 // Set the value of the intrinsic coordinate
505 zeta[0] = inter_zeta;
506 }
507
508
509 //===================================================================
510 /// If a neighbouring element has already created a node at
511 /// a position corresponding to the local fractional position within the
512 /// present element, s_fraction, return
513 /// a pointer to that node. If not, return NULL (0). If the node is
514 /// periodic the flag is_periodic will be true
515 //===================================================================
517 const Vector<double>& s_fraction, bool& is_periodic)
518 {
519 using namespace QuadTreeNames;
520
521 // Calculate the edges on which the node lies
522 Vector<int> edges;
523 if (s_fraction[0] == 0.0)
524 {
525 edges.push_back(W);
526 }
527 if (s_fraction[0] == 1.0)
528 {
529 edges.push_back(E);
530 }
531 if (s_fraction[1] == 0.0)
532 {
533 edges.push_back(S);
534 }
535 if (s_fraction[1] == 1.0)
536 {
537 edges.push_back(N);
538 }
539
540 // Find the number of edges
541 unsigned n_size = edges.size();
542 // If there are no edges, then there is no neighbour, return 0
543 if (n_size == 0)
544 {
545 return 0;
546 }
547
548 Vector<unsigned> translate_s(2);
549 Vector<double> s_lo_neigh(2);
550 Vector<double> s_hi_neigh(2);
551 Vector<double> s(2);
552
553 int neigh_edge, diff_level;
554 bool in_neighbouring_tree;
555
556 // Loop over the edges
557 for (unsigned j = 0; j < n_size; j++)
558 {
559 // Find pointer to neighbouring element along edge
560 QuadTree* neigh_pt;
561 neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[j],
562 translate_s,
563 s_lo_neigh,
564 s_hi_neigh,
565 neigh_edge,
566 diff_level,
567 in_neighbouring_tree);
568
569 // Neighbour exists
570 if (neigh_pt != 0)
571 {
572 // Have its nodes been created yet?
573 // if(neigh_pt->object_pt()->node_pt(0)!=0)
574 if (neigh_pt->object_pt()->nodes_built())
575 {
576 // We now need to translate the nodal location
577 // as defined in terms of the fractional coordinates of the present
578 // element into those of its neighbour
579
580 // Calculate the local coordinate in the neighbour
581 // Note that we need to use the translation scheme in case
582 // the local coordinates are swapped in the neighbour.
583 for (unsigned i = 0; i < 2; i++)
584 {
585 s[i] = s_lo_neigh[i] +
586 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
587 }
588
589 // Find the node in the neighbour
590 Node* neighbour_node_pt =
592
593 // If there is a node, return it
594 if (neighbour_node_pt != 0)
595 {
596 // Now work out whether it's a periodic boundary
597 // only possible if we have moved into a neighbouring tree
598 if (in_neighbouring_tree)
599 {
600 // Return whether the neighbour is periodic
601 is_periodic =
602 quadtree_pt()->root_pt()->is_neighbour_periodic(edges[j]);
603 }
604 // Return the pointer to the neighbouring node
605 return neighbour_node_pt;
606 }
607 }
608 }
609 }
610 // Node not found, return null
611 return 0;
612 }
613
614 //==================================================================
615 /// Build the element by doing the following:
616 /// - Give it nodal positions (by establishing the pointers to its
617 /// nodes)
618 /// - In the process create new nodes where required (i.e. if they
619 /// don't exist in father element or have already been created
620 /// while building new neighbour elements). Node building
621 /// involves the following steps:
622 /// - Get nodal position from father element.
623 /// - Establish the time-history of the newly created nodal point
624 /// (its coordinates and the previous values) consistent with
625 /// the father's history.
626 /// - Determine the boundary conditions of the nodes (newly
627 /// created nodes can only lie on the interior of any
628 /// edges of the father element -- this makes it possible to
629 /// to figure out what their bc should be...)
630 /// - Add node to the mesh's stoarge scheme for the boundary nodes.
631 /// - Add the new node to the mesh itself
632 /// - Doc newly created nodes in "new_nodes.dat" stored in the directory
633 /// of the DocInfo object (only if it's open!)
634 /// - Finally, excute the element-specific further_build()
635 /// (empty by default -- must be overloaded for specific elements).
636 /// This deals with any build operations that are not included
637 /// in the generic process outlined above. For instance, in
638 /// Crouzeix Raviart elements we need to initialise the internal
639 /// pressure values in manner consistent with the pressure
640 /// distribution in the father element.
641 //==================================================================
643 Vector<Node*>& new_node_pt,
644 bool& was_already_built,
645 std::ofstream& new_nodes_file)
646 {
647 using namespace QuadTreeNames;
648
649 // Get the number of 1d nodes
650 unsigned n_p = nnode_1d();
651
652 // Check whether static father_bound needs to be created
653 if (Father_bound[n_p].nrow() == 0)
654 {
655 setup_father_bounds();
656 }
657
658 // Pointer to my father (in quadtree impersonation)
659 QuadTree* father_pt = dynamic_cast<QuadTree*>(quadtree_pt()->father_pt());
660
661 // What type of son am I? Ask my quadtree representation...
662 int son_type = Tree_pt->son_type();
663
664 // Has somebody build me already? (If any nodes have not been built)
665 if (!nodes_built())
666 {
667#ifdef PARANOID
668 if (father_pt == 0)
669 {
670 std::string error_message =
671 "Something fishy here: I have no father and yet \n";
672 error_message += "I have no nodes. Who has created me then?!\n";
673
674 throw OomphLibError(
675 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
676 }
677#endif
678
679 // Indicate status:
680 was_already_built = false;
681
682 // Return pointer to father element
683 RefineableQElement<2>* father_el_pt =
684 dynamic_cast<RefineableQElement<2>*>(father_pt->object_pt());
685
686 // Timestepper should be the same for all nodes in father
687 // element -- use it create timesteppers for new nodes
688 TimeStepper* time_stepper_pt =
689 father_el_pt->node_pt(0)->time_stepper_pt();
690
691 // Number of history values (incl. present)
692 unsigned ntstorage = time_stepper_pt->ntstorage();
693
694 // Currently we can't handle the case of generalised coordinates
695 // since we haven't established how they should be interpolated
696 // Buffer this case:
697 if (father_el_pt->node_pt(0)->nposition_type() != 1)
698 {
699 throw OomphLibError("Can't handle generalised nodal positions (yet).",
700 OOMPH_CURRENT_FUNCTION,
701 OOMPH_EXCEPTION_LOCATION);
702 }
703
704 Vector<double> s_lo(2);
705 Vector<double> s_hi(2);
706 Vector<double> s(2);
707 Vector<double> x(2);
708
709 // Setup vertex coordinates in father element:
710 //--------------------------------------------
711 switch (son_type)
712 {
713 case SW:
714 s_lo[0] = -1.0;
715 s_hi[0] = 0.0;
716 s_lo[1] = -1.0;
717 s_hi[1] = 0.0;
718 break;
719
720 case SE:
721 s_lo[0] = 0.0;
722 s_hi[0] = 1.0;
723 s_lo[1] = -1.0;
724 s_hi[1] = 0.0;
725 break;
726
727 case NE:
728 s_lo[0] = 0.0;
729 s_hi[0] = 1.0;
730 s_lo[1] = 0.0;
731 s_hi[1] = 1.0;
732 break;
733
734 case NW:
735 s_lo[0] = -1.0;
736 s_hi[0] = 0.0;
737 s_lo[1] = 0.0;
738 s_hi[1] = 1.0;
739 break;
740 }
741
742 // Pass macro element pointer on to sons and
743 // set coordinates in macro element
744 // hierher why can I see this?
745 if (father_el_pt->Macro_elem_pt != 0)
746 {
747 set_macro_elem_pt(father_el_pt->Macro_elem_pt);
748 for (unsigned i = 0; i < 2; i++)
749 {
750 s_macro_ll(i) =
751 father_el_pt->s_macro_ll(i) +
752 0.5 * (s_lo[i] + 1.0) *
753 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
754 s_macro_ur(i) =
755 father_el_pt->s_macro_ll(i) +
756 0.5 * (s_hi[i] + 1.0) *
757 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
758 }
759 }
760
761
762 // If the father element hasn't been generated yet, we're stuck...
763 if (father_el_pt->node_pt(0) == 0)
764 {
765 throw OomphLibError(
766 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
767 OOMPH_CURRENT_FUNCTION,
768 OOMPH_EXCEPTION_LOCATION);
769 }
770 else
771 {
772 unsigned jnod = 0;
773 Vector<double> x_small(2);
774 Vector<double> x_large(2);
775
776 Vector<double> s_fraction(2);
777 // Loop over nodes in element
778 for (unsigned i0 = 0; i0 < n_p; i0++)
779 {
780 // Get the fractional position of the node in the direction of s[0]
781 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
782 // Local coordinate in father element
783 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
784
785 for (unsigned i1 = 0; i1 < n_p; i1++)
786 {
787 // Get the fractional position of the node in the direction of s[1]
788 s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
789 // Local coordinate in father element
790 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
791
792 // Local node number
793 jnod = i0 + n_p * i1;
794
795 // Check whether the father's node is periodic if so, complain
796 /* {
797 Node* father_node_pt = father_el_pt->node_pt(jnod);
798 if((father_node_pt->is_a_copy()) ||
799 (father_node_pt->position_is_a_copy()))
800 {
801 throw OomphLibError(
802 "Can't handle periodic nodes (yet).",
803 OOMPH_CURRENT_FUNCTION,
804 OOMPH_EXCEPTION_LOCATION);
805 }
806 }*/
807
808 // Initialise flag: So far, this node hasn't been built
809 // or copied yet
810 bool node_done = false;
811
812 // Get the pointer to the node in the father, returns NULL
813 // if there is not node
814 Node* created_node_pt =
815 father_el_pt->get_node_at_local_coordinate(s);
816
817 // Does this node already exist in father element?
818 //------------------------------------------------
819 if (created_node_pt != 0)
820 {
821 // Copy node across
822 node_pt(jnod) = created_node_pt;
823
824 // Make sure that we update the values at the node so that
825 // they are consistent with the present representation.
826 // This is only need for mixed interpolation where the value
827 // at the father could now become active.
828
829 // Loop over all history values
830 for (unsigned t = 0; t < ntstorage; t++)
831 {
832 // Get values from father element
833 // Note: get_interpolated_values() sets Vector size itself.
834 Vector<double> prev_values;
835 father_el_pt->get_interpolated_values(t, s, prev_values);
836 // Find the minimum number of values
837 //(either those stored at the node, or those returned by
838 // the function)
839 unsigned n_val_at_node = created_node_pt->nvalue();
840 unsigned n_val_from_function = prev_values.size();
841 // Use the ternary conditional operator here
842 unsigned n_var = n_val_at_node < n_val_from_function ?
843 n_val_at_node :
844 n_val_from_function;
845 // Assign the values that we can
846 for (unsigned k = 0; k < n_var; k++)
847 {
848 created_node_pt->set_value(t, k, prev_values[k]);
849 }
850 }
851
852 // Node has been created by copy
853 node_done = true;
854 }
855 // Node does not exist in father element but might already
856 //--------------------------------------------------------
857 // have been created by neighbouring elements
858 //-------------------------------------------
859 else
860 {
861 // Was the node created by one of its neighbours
862 // Whether or not the node lies on an edge can be calculated
863 // by from the fractional position
864 bool is_periodic = false;
865 ;
866 created_node_pt =
867 node_created_by_neighbour(s_fraction, is_periodic);
868
869 // If the node was so created, assign the pointers
870 if (created_node_pt != 0)
871 {
872 // If the node is periodic
873 if (is_periodic)
874 {
875 // Now the node must be on a boundary, but we don't know which
876 // one
877 // The returned created_node_pt is actually the neighbouring
878 // periodic node
879 Node* neighbour_node_pt = created_node_pt;
880
881 // Determine the edge on which the new node will live
882 int father_bound = Father_bound[n_p](jnod, son_type);
883
884 // Storage for the set of Mesh boundaries on which the
885 // appropriate father edge lives.
886 // [New nodes should always be mid-edge nodes in father
887 // and therefore only live on one boundary but just to
888 // play it safe...]
889 std::set<unsigned> boundaries;
890 // Only get the boundaries if we are at the edge of
891 // an element. Nodes in the centre of an element cannot be
892 // on Mesh boundaries
893 if (father_bound != Tree::OMEGA)
894 {
895 father_el_pt->get_boundaries(father_bound, boundaries);
896 }
897
898#ifdef PARANOID
899 // Case where a new node lives on more than one boundary
900 // seems fishy enough to flag
901 if (boundaries.size() > 1)
902 {
903 throw OomphLibError(
904 "boundaries.size()!=1 seems a bit strange..\n",
905 OOMPH_CURRENT_FUNCTION,
906 OOMPH_EXCEPTION_LOCATION);
907 }
908
909 // Case when there are no boundaries, we are in big trouble
910 if (boundaries.size() == 0)
911 {
912 std::ostringstream error_stream;
913 error_stream << "Periodic node is not on a boundary...\n"
914 << "Coordinates: " << created_node_pt->x(0)
915 << " " << created_node_pt->x(1) << "\n";
916 throw OomphLibError(error_stream.str(),
917 OOMPH_CURRENT_FUNCTION,
918 OOMPH_EXCEPTION_LOCATION);
919 }
920#endif
921
922 // Create node and set the pointer to it from the element
923 created_node_pt =
924 construct_boundary_node(jnod, time_stepper_pt);
925 // Make the node periodic from the neighbour
926 created_node_pt->make_periodic(neighbour_node_pt);
927 // Add to vector of new nodes
928 new_node_pt.push_back(created_node_pt);
929
930 // Loop over # of history values
931 for (unsigned t = 0; t < ntstorage; t++)
932 {
933 // Get position from father element -- this uses the macro
934 // element representation if appropriate. If the node
935 // turns out to be a hanging node later on, then
936 // its position gets adjusted in line with its
937 // hanging node interpolation.
938 Vector<double> x_prev(2);
939 father_el_pt->get_x(t, s, x_prev);
940 // Set previous positions of the new node
941 for (unsigned i = 0; i < 2; i++)
942 {
943 created_node_pt->x(t, i) = x_prev[i];
944 }
945 }
946
947 // Next, we Update the boundary lookup schemes
948 // Loop over the boundaries stored in the set
949 for (std::set<unsigned>::iterator it = boundaries.begin();
950 it != boundaries.end();
951 ++it)
952 {
953 // Add the node to the boundary
954 mesh_pt->add_boundary_node(*it, created_node_pt);
955
956 // If we have set an intrinsic coordinate on this
957 // mesh boundary then it must also be interpolated on
958 // the new node
959 // Now interpolate the intrinsic boundary coordinate
960 if (mesh_pt->boundary_coordinate_exists(*it) == true)
961 {
962 Vector<double> zeta(1);
963 father_el_pt->interpolated_zeta_on_edge(
964 *it, father_bound, s, zeta);
965
966 created_node_pt->set_coordinates_on_boundary(*it, zeta);
967 }
968 }
969
970 // Make sure that we add the node to the mesh
971 mesh_pt->add_node_pt(created_node_pt);
972 } // End of periodic case
973 // Otherwise the node is not periodic, so just set the
974 // pointer to the neighbours node
975 else
976 {
977 node_pt(jnod) = created_node_pt;
978 }
979 // Node has been created
980 node_done = true;
981 }
982 // Node does not exist in neighbour element but might already
983 //-----------------------------------------------------------
984 // have been created by a son of a neighbouring element
985 //-----------------------------------------------------
986 else
987 {
988 // Was the node created by one of its neighbours' sons
989 // Whether or not the node lies on an edge can be calculated
990 // by from the fractional position
991 bool is_periodic = false;
992 ;
993 created_node_pt =
994 node_created_by_son_of_neighbour(s_fraction, is_periodic);
995
996 // If the node was so created, assign the pointers
997 if (created_node_pt != 0)
998 {
999 // If the node is periodic
1000 if (is_periodic)
1001 {
1002 // Now the node must be on a boundary, but we don't know
1003 // which one The returned created_node_pt is actually the
1004 // neighbouring periodic node
1005 Node* neighbour_node_pt = created_node_pt;
1006
1007 // Determine the edge on which the new node will live
1008 int father_bound = Father_bound[n_p](jnod, son_type);
1009
1010 // Storage for the set of Mesh boundaries on which the
1011 // appropriate father edge lives.
1012 // [New nodes should always be mid-edge nodes in father
1013 // and therefore only live on one boundary but just to
1014 // play it safe...]
1015 std::set<unsigned> boundaries;
1016 // Only get the boundaries if we are at the edge of
1017 // an element. Nodes in the centre of an element cannot be
1018 // on Mesh boundaries
1019 if (father_bound != Tree::OMEGA)
1020 {
1021 father_el_pt->get_boundaries(father_bound, boundaries);
1022 }
1023
1024#ifdef PARANOID
1025 // Case where a new node lives on more than one boundary
1026 // seems fishy enough to flag
1027 if (boundaries.size() > 1)
1028 {
1029 throw OomphLibError(
1030 "boundaries.size()!=1 seems a bit strange..\n",
1031 OOMPH_CURRENT_FUNCTION,
1032 OOMPH_EXCEPTION_LOCATION);
1033 }
1034
1035 // Case when there are no boundaries, we are in big trouble
1036 if (boundaries.size() == 0)
1037 {
1038 std::ostringstream error_stream;
1039 error_stream << "Periodic node is not on a boundary...\n"
1040 << "Coordinates: " << created_node_pt->x(0)
1041 << " " << created_node_pt->x(1) << "\n";
1042 throw OomphLibError(error_stream.str(),
1043 OOMPH_CURRENT_FUNCTION,
1044 OOMPH_EXCEPTION_LOCATION);
1045 }
1046#endif
1047
1048 // Create node and set the pointer to it from the element
1049 created_node_pt =
1050 construct_boundary_node(jnod, time_stepper_pt);
1051 // Make the node periodic from the neighbour
1052 created_node_pt->make_periodic(neighbour_node_pt);
1053 // Add to vector of new nodes
1054 new_node_pt.push_back(created_node_pt);
1055
1056 // Loop over # of history values
1057 for (unsigned t = 0; t < ntstorage; t++)
1058 {
1059 // Get position from father element -- this uses the macro
1060 // element representation if appropriate. If the node
1061 // turns out to be a hanging node later on, then
1062 // its position gets adjusted in line with its
1063 // hanging node interpolation.
1064 Vector<double> x_prev(2);
1065 father_el_pt->get_x(t, s, x_prev);
1066 // Set previous positions of the new node
1067 for (unsigned i = 0; i < 2; i++)
1068 {
1069 created_node_pt->x(t, i) = x_prev[i];
1070 }
1071 }
1072
1073 // Next, we Update the boundary lookup schemes
1074 // Loop over the boundaries stored in the set
1075 for (std::set<unsigned>::iterator it = boundaries.begin();
1076 it != boundaries.end();
1077 ++it)
1078 {
1079 // Add the node to the boundary
1080 mesh_pt->add_boundary_node(*it, created_node_pt);
1081
1082 // If we have set an intrinsic coordinate on this
1083 // mesh boundary then it must also be interpolated on
1084 // the new node
1085 // Now interpolate the intrinsic boundary coordinate
1086 if (mesh_pt->boundary_coordinate_exists(*it) == true)
1087 {
1088 Vector<double> zeta(1);
1089 father_el_pt->interpolated_zeta_on_edge(
1090 *it, father_bound, s, zeta);
1091
1092 created_node_pt->set_coordinates_on_boundary(*it, zeta);
1093 }
1094 }
1095
1096 // Make sure that we add the node to the mesh
1097 mesh_pt->add_node_pt(created_node_pt);
1098 } // End of periodic case
1099 // Otherwise the node is not periodic, so just set the
1100 // pointer to the neighbours node
1101 else
1102 {
1103 node_pt(jnod) = created_node_pt;
1104 }
1105 // Node has been created
1106 node_done = true;
1107 } // Node does not exist in son of neighbouring element
1108 } // Node does not exist in neighbouring element
1109 } // Node does not exist in father element
1110
1111 // Node has not been built anywhere ---> build it here
1112 if (!node_done)
1113 {
1114 // Firstly, we need to determine whether or not a node lies
1115 // on the boundary before building it, because
1116 // we actually assign a different type of node on boundaries.
1117
1118 // The node can only be on a Mesh boundary if it
1119 // lives on an edge that is shared with an edge of its
1120 // father element; i.e. it is not created inside the father
1121 // element Determine the edge on which the new node will live
1122 int father_bound = Father_bound[n_p](jnod, son_type);
1123
1124 // Storage for the set of Mesh boundaries on which the
1125 // appropriate father edge lives.
1126 // [New nodes should always be mid-edge nodes in father
1127 // and therefore only live on one boundary but just to
1128 // play it safe...]
1129 std::set<unsigned> boundaries;
1130 // Only get the boundaries if we are at the edge of
1131 // an element. Nodes in the centre of an element cannot be
1132 // on Mesh boundaries
1133 if (father_bound != Tree::OMEGA)
1134 {
1135 father_el_pt->get_boundaries(father_bound, boundaries);
1136 }
1137
1138#ifdef PARANOID
1139 // Case where a new node lives on more than one boundary
1140 // seems fishy enough to flag
1141 if (boundaries.size() > 1)
1142 {
1143 throw OomphLibError(
1144 "boundaries.size()!=1 seems a bit strange..\n",
1145 OOMPH_CURRENT_FUNCTION,
1146 OOMPH_EXCEPTION_LOCATION);
1147 }
1148#endif
1149
1150 // If the node lives on a mesh boundary,
1151 // then we need to create a boundary node
1152 if (boundaries.size() > 0)
1153 {
1154 // Create node and set the pointer to it from the element
1155 created_node_pt =
1156 construct_boundary_node(jnod, time_stepper_pt);
1157 // Add to vector of new nodes
1158 new_node_pt.push_back(created_node_pt);
1159
1160 // Now we need to work out whether to pin the values at
1161 // the new node based on the boundary conditions applied at
1162 // its Mesh boundary
1163
1164 // Get the boundary conditions from the father
1165 Vector<int> bound_cons(ncont_interpolated_values());
1166 father_el_pt->get_bcs(father_bound, bound_cons);
1167
1168 // Loop over the values and pin, if necessary
1169 unsigned n_value = created_node_pt->nvalue();
1170 for (unsigned k = 0; k < n_value; k++)
1171 {
1172 if (bound_cons[k])
1173 {
1174 created_node_pt->pin(k);
1175 }
1176 }
1177
1178 // Solid node? If so, deal with the positional boundary
1179 // conditions:
1180 SolidNode* solid_node_pt =
1181 dynamic_cast<SolidNode*>(created_node_pt);
1182 if (solid_node_pt != 0)
1183 {
1184 // Get the positional boundary conditions from the father:
1185 unsigned n_dim = created_node_pt->ndim();
1186 Vector<int> solid_bound_cons(n_dim);
1187 RefineableSolidQElement<2>* father_solid_el_pt =
1188 dynamic_cast<RefineableSolidQElement<2>*>(father_el_pt);
1189#ifdef PARANOID
1190 if (father_solid_el_pt == 0)
1191 {
1192 std::string error_message =
1193 "We have a SolidNode outside a refineable SolidElement\n";
1194 error_message +=
1195 "during mesh refinement -- this doesn't make sense";
1196
1197 throw OomphLibError(error_message,
1198 OOMPH_CURRENT_FUNCTION,
1199 OOMPH_EXCEPTION_LOCATION);
1200 }
1201#endif
1202 father_solid_el_pt->get_solid_bcs(father_bound,
1203 solid_bound_cons);
1204
1205 // Loop over the positions and pin, if necessary
1206 for (unsigned k = 0; k < n_dim; k++)
1207 {
1208 if (solid_bound_cons[k])
1209 {
1210 solid_node_pt->pin_position(k);
1211 }
1212 }
1213 } // End of if solid_node_pt
1214
1215
1216 // Next, we Update the boundary lookup schemes
1217 // Loop over the boundaries stored in the set
1218 for (std::set<unsigned>::iterator it = boundaries.begin();
1219 it != boundaries.end();
1220 ++it)
1221 {
1222 // Add the node to the boundary
1223 mesh_pt->add_boundary_node(*it, created_node_pt);
1224
1225 // If we have set an intrinsic coordinate on this
1226 // mesh boundary then it must also be interpolated on
1227 // the new node
1228 // Now interpolate the intrinsic boundary coordinate
1229 if (mesh_pt->boundary_coordinate_exists(*it) == true)
1230 {
1231 Vector<double> zeta(1);
1232 father_el_pt->interpolated_zeta_on_edge(
1233 *it, father_bound, s, zeta);
1234
1235 created_node_pt->set_coordinates_on_boundary(*it, zeta);
1236 }
1237 }
1238 }
1239 // Otherwise the node is not on a Mesh boundary and
1240 // we create a normal "bulk" node
1241 else
1242 {
1243 // Create node and set the pointer to it from the element
1244 created_node_pt = construct_node(jnod, time_stepper_pt);
1245 // Add to vector of new nodes
1246 new_node_pt.push_back(created_node_pt);
1247 }
1248
1249 // Now we set the position and values at the newly created node
1250
1251 // In the first instance use macro element or FE representation
1252 // to create past and present nodal positions.
1253 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
1254 // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
1255 // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
1256 // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
1257 // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1258 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
1259 // NOT ASSIGN SENSIBLE INITIAL POSITONS!
1260
1261 // Loop over # of history values
1262 for (unsigned t = 0; t < ntstorage; t++)
1263 {
1264 // Get position from father element -- this uses the macro
1265 // element representation if appropriate. If the node
1266 // turns out to be a hanging node later on, then
1267 // its position gets adjusted in line with its
1268 // hanging node interpolation.
1269 Vector<double> x_prev(2);
1270 father_el_pt->get_x(t, s, x_prev);
1271
1272 // Set previous positions of the new node
1273 for (unsigned i = 0; i < 2; i++)
1274 {
1275 created_node_pt->x(t, i) = x_prev[i];
1276 }
1277 }
1278
1279 // Loop over all history values
1280 for (unsigned t = 0; t < ntstorage; t++)
1281 {
1282 // Get values from father element
1283 // Note: get_interpolated_values() sets Vector size itself.
1284 Vector<double> prev_values;
1285 father_el_pt->get_interpolated_values(t, s, prev_values);
1286 // Initialise the values at the new node
1287 unsigned n_value = created_node_pt->nvalue();
1288 for (unsigned k = 0; k < n_value; k++)
1289 {
1290 created_node_pt->set_value(t, k, prev_values[k]);
1291 }
1292 }
1293
1294 // Add new node to mesh
1295 mesh_pt->add_node_pt(created_node_pt);
1296
1297 } // End of case when we build the node ourselves
1298
1299 // Check if the element is an algebraic element
1300 AlgebraicElementBase* alg_el_pt =
1301 dynamic_cast<AlgebraicElementBase*>(this);
1302
1303 // If the element is an algebraic element, setup
1304 // node position (past and present) from algebraic node update
1305 // function. This over-writes previous assingments that
1306 // were made based on the macro-element/FE representation.
1307 // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
1308 // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
1309 // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
1310 // INFO FOR *ALL* ROOT ELEMENTS!
1311 if (alg_el_pt != 0)
1312 {
1313 // Build algebraic node update info for new node
1314 // This sets up the node update data for all node update
1315 // functions that are shared by all nodes in the father
1316 // element
1317 alg_el_pt->setup_algebraic_node_update(
1318 node_pt(jnod), s, father_el_pt);
1319 }
1320
1321 // If we have built the node and we are documenting our progress
1322 // write the (hopefully consistent position) to the outputfile
1323 if ((!node_done) && (new_nodes_file.is_open()))
1324 {
1325 new_nodes_file << node_pt(jnod)->x(0) << " "
1326 << node_pt(jnod)->x(1) << std::endl;
1327 }
1328
1329 } // End of vertical loop over nodes in element
1330
1331 } // End of horizontal loop over nodes in element
1332
1333
1334 // If the element is a MacroElementNodeUpdateElement, set
1335 // the update parameters for the current element's nodes --
1336 // all this needs is the vector of (pointers to the)
1337 // geometric objects that affect the MacroElement-based
1338 // node update -- this is the same as that in the father element
1339 MacroElementNodeUpdateElementBase* father_m_el_pt =
1340 dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
1341 if (father_m_el_pt != 0)
1342 {
1343 // Get vector of geometric objects from father (construct vector
1344 // via copy operation)
1345 Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
1346
1347 // Cast current element to MacroElementNodeUpdateElement:
1349 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1350
1351#ifdef PARANOID
1352 if (m_el_pt == 0)
1353 {
1354 std::string error_message =
1355 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
1356 error_message +=
1357 "Strange -- if the father is a MacroElementNodeUpdateElement\n";
1358 error_message += "the son should be too....\n";
1359
1360 throw OomphLibError(
1361 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1362 }
1363#endif
1364 // Build update info by passing vector of geometric objects:
1365 // This sets the current element to be the update element
1366 // for all of the element's nodes -- this is reversed
1367 // if the element is ever un-refined in the father element's
1368 // rebuild_from_sons() function which overwrites this
1369 // assignment to avoid nasty segmentation faults that occur
1370 // when a node tries to update itself via an element that no
1371 // longer exists...
1372 m_el_pt->set_node_update_info(geom_object_pt);
1373 }
1374
1375#ifdef OOMPH_HAS_MPI
1376 // Pass on non-halo proc id
1377 Non_halo_proc_ID =
1378 tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
1379#endif
1380
1381 // Is it an ElementWithMovingNodes?
1382 ElementWithMovingNodes* aux_el_pt =
1383 dynamic_cast<ElementWithMovingNodes*>(this);
1384
1385 // Pass down the information re the method for the evaluation
1386 // of the shape derivatives
1387 if (aux_el_pt != 0)
1388 {
1389 ElementWithMovingNodes* aux_father_el_pt =
1390 dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
1391
1392#ifdef PARANOID
1393 if (aux_father_el_pt == 0)
1394 {
1395 std::string error_message =
1396 "Failed to cast to ElementWithMovingNodes*\n";
1397 error_message +=
1398 "Strange -- if the son is a ElementWithMovingNodes\n";
1399 error_message += "the father should be too....\n";
1400
1401 throw OomphLibError(
1402 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1403 }
1404#endif
1405
1406 // If evaluating the residuals by finite differences in the father
1407 // continue to do so in the child
1408 if (aux_father_el_pt
1409 ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
1410 {
1411 aux_el_pt
1413 }
1414
1415 aux_el_pt->method_for_shape_derivs() =
1416 aux_father_el_pt->method_for_shape_derivs();
1417
1418 // If bypassing the evaluation of fill_in_jacobian_from_geometric_data
1419 // continue to do so
1420 if (aux_father_el_pt
1421 ->is_fill_in_jacobian_from_geometric_data_bypassed())
1422 {
1424 }
1425 }
1426
1427
1428 // Now do further build (if any)
1429 further_build();
1430
1431 } // Sanity check: Father element has been generated
1432 }
1433 // Element has already been built
1434 else
1435 {
1436 was_already_built = true;
1437 }
1438 }
1439
1440 //====================================================================
1441 /// Print corner nodes, use colour (default "BLACK")
1442 //====================================================================
1443 void RefineableQElement<2>::output_corners(std::ostream& outfile,
1444 const std::string& colour) const
1445 {
1446 Vector<double> s(2);
1447 Vector<double> corner(2);
1448
1449 outfile << "ZONE I=2,J=2, C=" << colour << std::endl;
1450
1451 s[0] = -1.0;
1452 s[1] = -1.0;
1453 get_x(s, corner);
1454 outfile << corner[0] << " " << corner[1] << " " << Number << std::endl;
1455
1456 s[0] = 1.0;
1457 s[1] = -1.0;
1458 get_x(s, corner);
1459 outfile << corner[0] << " " << corner[1] << " " << Number << std::endl;
1460
1461 s[0] = -1.0;
1462 s[1] = 1.0;
1463 get_x(s, corner);
1464 outfile << corner[0] << " " << corner[1] << " " << Number << std::endl;
1465
1466 s[0] = 1.0;
1467 s[1] = 1.0;
1468 get_x(s, corner);
1469 outfile << corner[0] << " " << corner[1] << " " << Number << std::endl;
1470
1471
1472 outfile << "TEXT CS = GRID, X = " << corner[0] << ",Y = " << corner[1]
1473 << ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\"" << Number << "\""
1474 << std::endl;
1475 }
1476
1477 //====================================================================
1478 /// Set up all hanging nodes. If we are documenting the output then
1479 /// open the output files and pass the open files to the helper function
1480 //====================================================================
1482 Vector<std::ofstream*>& output_stream)
1483 {
1484#ifdef PARANOID
1485 if (output_stream.size() != 4)
1486 {
1487 throw OomphLibError("There must be four output streams",
1488 OOMPH_CURRENT_FUNCTION,
1489 OOMPH_EXCEPTION_LOCATION);
1490 }
1491#endif
1492
1493 using namespace QuadTreeNames;
1494
1495 // Setup hanging nodes on each edge of the element
1496 quad_hang_helper(-1, S, *(output_stream[0]));
1497 quad_hang_helper(-1, N, *(output_stream[1]));
1498 quad_hang_helper(-1, W, *(output_stream[2]));
1499 quad_hang_helper(-1, E, *(output_stream[3]));
1500 }
1501
1502 //================================================================
1503 /// Internal function that sets up the hanging node scheme for
1504 /// a particular continuously interpolated value
1505 //===============================================================
1507 {
1508 using namespace QuadTreeNames;
1509
1510 std::ofstream dummy_hangfile;
1511 quad_hang_helper(value_id, S, dummy_hangfile);
1512 quad_hang_helper(value_id, N, dummy_hangfile);
1513 quad_hang_helper(value_id, W, dummy_hangfile);
1514 quad_hang_helper(value_id, E, dummy_hangfile);
1515 }
1516
1517
1518 //=================================================================
1519 /// Internal function to set up the hanging nodes on a particular
1520 /// edge of the element
1521 //=================================================================
1523 const int& my_edge,
1524 std::ofstream& output_hangfile)
1525 {
1526 using namespace QuadTreeNames;
1527
1528 Vector<unsigned> translate_s(2);
1529 Vector<double> s_lo_neigh(2);
1530 Vector<double> s_hi_neigh(2);
1531 int neigh_edge, diff_level;
1532 bool in_neighbouring_tree;
1533
1534 // Find pointer to neighbour in this direction
1535 QuadTree* neigh_pt;
1536 neigh_pt = quadtree_pt()->gteq_edge_neighbour(my_edge,
1537 translate_s,
1538 s_lo_neigh,
1539 s_hi_neigh,
1540 neigh_edge,
1541 diff_level,
1542 in_neighbouring_tree);
1543
1544 // Neighbour exists and all nodes have been created
1545 if (neigh_pt != 0)
1546 {
1547 // Different sized element?
1548 if (diff_level != 0)
1549 {
1550 // Test for the periodic node case
1551 // Are we crossing a periodic boundary
1552 bool is_periodic = false;
1553 if (in_neighbouring_tree)
1554 {
1555 is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(my_edge);
1556 }
1557
1558 // If it is periodic we actually need to get the node in
1559 // the neighbour of the neighbour (which will be a parent of
1560 // the present element) so that the "fixed" coordinate is
1561 // correctly calculated.
1562 // The idea is to replace the neigh_pt and associated data
1563 // with those of the neighbour of the neighbour
1564 if (is_periodic)
1565 {
1566 // Required data for the neighbour finding routine
1567 Vector<unsigned> translate_s_in_neigh(2);
1568 Vector<double> s_lo_neigh_of_neigh(2);
1569 Vector<double> s_hi_neigh_of_neigh(2);
1570 int neigh_edge_of_neigh, diff_level_of_neigh;
1571 bool in_neighbouring_tree_of_neigh;
1572
1573 // Find pointer to neighbour of the neighbour on the edge
1574 // that we are currently considering
1575 QuadTree* neigh_of_neigh_pt;
1576 neigh_of_neigh_pt =
1577 neigh_pt->gteq_edge_neighbour(neigh_edge,
1578 translate_s_in_neigh,
1579 s_lo_neigh_of_neigh,
1580 s_hi_neigh_of_neigh,
1581 neigh_edge_of_neigh,
1582 diff_level_of_neigh,
1583 in_neighbouring_tree_of_neigh);
1584
1585 // Set the value of the NEW neighbour and edge
1586 neigh_pt = neigh_of_neigh_pt;
1587 neigh_edge = neigh_edge_of_neigh;
1588
1589 // Set the values of the translation schemes
1590 // Need to find the values of s_lo and s_hi
1591 // in the neighbour of the neighbour
1592
1593 // Get the minimum and maximum values of the coordinate
1594 // in the neighbour (don't like this, but I think it's
1595 // necessary) Note that these values are hardcoded
1596 // in the quadtrees at some point!!
1597 double s_min = neigh_pt->object_pt()->s_min();
1598 double s_max = neigh_pt->object_pt()->s_max();
1599 Vector<double> s_lo_frac(2), s_hi_frac(2);
1600 // Work out the fractional position of the low and high points
1601 // of the original element
1602 for (unsigned i = 0; i < 2; i++)
1603 {
1604 s_lo_frac[i] = (s_lo_neigh[i] - s_min) / (s_max - s_min);
1605 s_hi_frac[i] = (s_hi_neigh[i] - s_min) / (s_max - s_min);
1606 }
1607
1608 // We should now be able to construct the low and high points in
1609 // the neighbour of the neighbour
1610 for (unsigned i = 0; i < 2; i++)
1611 {
1612 s_lo_neigh[i] = s_lo_neigh_of_neigh[i] +
1613 s_lo_frac[translate_s_in_neigh[i]] *
1614 (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
1615 s_hi_neigh[i] = s_lo_neigh_of_neigh[i] +
1616 s_hi_frac[translate_s_in_neigh[i]] *
1617 (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
1618 }
1619
1620 // Finally we must sort out the translation scheme
1621 Vector<unsigned> temp_translate(2);
1622 for (unsigned i = 0; i < 2; i++)
1623 {
1624 temp_translate[i] = translate_s_in_neigh[translate_s[i]];
1625 }
1626 for (unsigned i = 0; i < 2; i++)
1627 {
1628 translate_s[i] = temp_translate[i];
1629 }
1630 } // End of special treatment for periodic hanging nodes
1631
1632 // Number of nodes in one dimension
1633 unsigned n_p = ninterpolating_node_1d(value_id);
1634 // Storage for the local nodes along the edge of the quadtree
1635 Node* local_node_pt = 0;
1636 // Loop over nodes along the edge
1637 for (unsigned i0 = 0; i0 < n_p; i0++)
1638 {
1639 // Storage for the fractional position of the node in the element
1640 Vector<double> s_fraction(2);
1641
1642 // Find the local node and the fractional position of the node
1643 // which depends on the edge, of course
1644 switch (my_edge)
1645 {
1646 case N:
1647 s_fraction[0] =
1648 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
1649 s_fraction[1] = 1.0;
1650 local_node_pt =
1651 interpolating_node_pt(i0 + n_p * (n_p - 1), value_id);
1652 break;
1653
1654 case S:
1655 s_fraction[0] =
1656 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
1657 s_fraction[1] = 0.0;
1658 local_node_pt = interpolating_node_pt(i0, value_id);
1659 break;
1660
1661 case E:
1662 s_fraction[0] = 1.0;
1663 s_fraction[1] =
1664 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
1665 local_node_pt =
1666 interpolating_node_pt(n_p - 1 + n_p * i0, value_id);
1667 break;
1668
1669 case W:
1670 s_fraction[0] = 0.0;
1671 s_fraction[1] =
1672 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
1673 local_node_pt = interpolating_node_pt(n_p * i0, value_id);
1674 break;
1675
1676 default:
1677 throw OomphLibError("my_edge not N, S, W, E\n",
1678 OOMPH_CURRENT_FUNCTION,
1679 OOMPH_EXCEPTION_LOCATION);
1680 }
1681
1682 // Calculate the local coordinates of the node in the neighbour
1683 Vector<double> s_in_neighb(2);
1684 for (unsigned i = 0; i < 2; i++)
1685 {
1686 s_in_neighb[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
1687 (s_hi_neigh[i] - s_lo_neigh[i]);
1688 }
1689
1690 // Find the Node in the neighbouring element
1691 Node* neighbouring_node_pt =
1693 s_in_neighb, value_id);
1694
1695 // If the neighbour does not have a node at this point
1696 if (0 == neighbouring_node_pt)
1697 {
1698 // Do we need to make a hanging node, we assume that we don't
1699 // initially
1700 bool make_hanging_node = false;
1701
1702 // If the node is not hanging geometrically, then we must make
1703 // it hang
1704 if (!local_node_pt->is_hanging())
1705 {
1706 make_hanging_node = true;
1707 }
1708 // Otherwise, it could be hanging geometrically, but still
1709 // require a different hanging scheme for this data value
1710 else
1711 {
1712 if (local_node_pt->hanging_pt(value_id) ==
1713 local_node_pt->hanging_pt())
1714 {
1715 make_hanging_node = true;
1716 }
1717 }
1718
1719 // If we do need to make the hanging node, then let's do it
1720 if (make_hanging_node == true)
1721 {
1722 // Cache refineable element used here
1723 RefineableElement* const obj_pt = neigh_pt->object_pt();
1724
1725 // Get shape functions in neighbour element
1726 Shape psi(obj_pt->ninterpolating_node(value_id));
1727 obj_pt->interpolating_basis(s_in_neighb, psi, value_id);
1728
1729 // Allocate the storage for the Hang pointer
1730 // which contains n_p nodes
1731 HangInfo* hang_pt = new HangInfo(n_p);
1732
1733 // Loop over nodes on edge in neighbour and mark them as nodes
1734 // that this node depends on
1735 unsigned n_neighbour;
1736
1737 // Number of nodes along edge in neighbour element
1738 for (unsigned n_edge = 0; n_edge < n_p; n_edge++)
1739 {
1740 switch (neigh_edge)
1741 {
1742 case N:
1743 n_neighbour = n_p * (n_p - 1) + n_edge;
1744 break;
1745
1746 case S:
1747 n_neighbour = n_edge;
1748 break;
1749
1750 case W:
1751 n_neighbour = n_p * n_edge;
1752 break;
1753
1754 case E:
1755 n_neighbour = n_p * n_edge + (n_p - 1);
1756 break;
1757
1758 default:
1759 throw OomphLibError("neigh_edge not N, S, W, E\n",
1760 OOMPH_CURRENT_FUNCTION,
1761 OOMPH_EXCEPTION_LOCATION);
1762 }
1763
1764 // Push back neighbouring node and weight into
1765 // Vector of (pointers to)
1766 // master nodes and weights
1767 // The weight is merely the value of the shape function
1768 // corresponding to the node in the neighbour
1769 hang_pt->set_master_node_pt(
1770 n_edge,
1771 obj_pt->interpolating_node_pt(n_neighbour, value_id),
1772 psi[n_neighbour]);
1773 }
1774
1775 // Now set the hanging data for the position
1776 // This also constrains the data values associated with the
1777 // value id
1778 local_node_pt->set_hanging_pt(hang_pt, value_id);
1779 }
1780
1781 // Dump the output if the file has been openeed
1782 if (output_hangfile.is_open())
1783 {
1784 output_hangfile << local_node_pt->x(0) << " "
1785 << local_node_pt->x(1) << std::endl;
1786 }
1787 }
1788 // Otherwise check that the nodes are the same
1789 else
1790 {
1791#ifdef PARANOID
1792 if (local_node_pt != neighbouring_node_pt)
1793 {
1794 std::ostringstream warning_stream;
1795 warning_stream << "SANITY CHECK in quad_hang_helper \n"
1796 << "Current node " << local_node_pt << " at "
1797 << "(" << local_node_pt->x(0) << ", "
1798 << local_node_pt->x(1) << ")"
1799 << " is not hanging and has " << std::endl
1800 << "Neighbour's node " << neighbouring_node_pt
1801 << " at "
1802 << "(" << neighbouring_node_pt->x(0) << ", "
1803 << neighbouring_node_pt->x(1) << ")" << std::endl
1804 << "even though the two should be "
1805 << "identical" << std::endl;
1806 OomphLibWarning(warning_stream.str(),
1807 "RefineableQElement<2>::quad_hang_helper()",
1808 OOMPH_EXCEPTION_LOCATION);
1809 }
1810#endif
1811 }
1812
1813 // If we are doing the position, then
1814 if (value_id == -1)
1815 {
1816 // Get the nodal position from neighbour element
1817 Vector<double> x_in_neighb(2);
1818 neigh_pt->object_pt()->interpolated_x(s_in_neighb, x_in_neighb);
1819
1820 // Fine adjust the coordinates (macro map will pick up boundary
1821 // accurately but will lead to different element edges)
1822 local_node_pt->x(0) = x_in_neighb[0];
1823 local_node_pt->x(1) = x_in_neighb[1];
1824 }
1825 }
1826 }
1827 }
1828 }
1829
1830
1831 //=================================================================
1832 /// Check inter-element continuity of
1833 /// - nodal positions
1834 /// - (nodally) interpolated function values
1835 //====================================================================
1836 // template<unsigned NNODE_1D>
1838 {
1839 using namespace QuadTreeNames;
1840
1841 // Number of nodes along edge
1842 unsigned n_p = nnode_1d();
1843
1844 // Number of timesteps (incl. present) for which continuity is
1845 // to be checked.
1846 unsigned n_time = 1;
1847
1848 // Initialise errors
1849 max_error = 0.0;
1850 Vector<double> max_error_x(2, 0.0);
1851 double max_error_val = 0.0;
1852
1853 Vector<int> edges(4);
1854 edges[0] = S;
1855 edges[1] = N;
1856 edges[2] = W;
1857 edges[3] = E;
1858
1859 // Loop over the edges
1860 for (unsigned edge_counter = 0; edge_counter < 4; edge_counter++)
1861 {
1862 Vector<unsigned> translate_s(2);
1863 Vector<double> s(2), s_lo_neigh(2), s_hi_neigh(2), s_fraction(2);
1864 int neigh_edge, diff_level;
1865 bool in_neighbouring_tree;
1866
1867 // Find pointer to neighbour in this direction
1868 QuadTree* neigh_pt;
1869 neigh_pt = quadtree_pt()->gteq_edge_neighbour(edges[edge_counter],
1870 translate_s,
1871 s_lo_neigh,
1872 s_hi_neigh,
1873 neigh_edge,
1874 diff_level,
1875 in_neighbouring_tree);
1876
1877 // Neighbour exists and has existing nodes
1878 if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
1879 {
1880 // Need to exclude periodic nodes from this check
1881 // There are only periodic nodes if we are in a neighbouring tree
1882 bool is_periodic = false;
1883 if (in_neighbouring_tree)
1884 {
1885 // Is it periodic
1886 is_periodic =
1887 tree_pt()->root_pt()->is_neighbour_periodic(edges[edge_counter]);
1888 }
1889
1890 // Loop over nodes along the edge
1891 for (unsigned i0 = 0; i0 < n_p; i0++)
1892 {
1893 // Storage for pointer to the local node
1894 Node* local_node_pt = 0;
1895
1896 switch (edge_counter)
1897 {
1898 case 0:
1899 // Local fraction of node
1900 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
1901 s_fraction[1] = 0.0;
1902 // Get pointer to local node
1903 local_node_pt = node_pt(i0);
1904 break;
1905
1906 case 1:
1907 // Local fraction of node
1908 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
1909 s_fraction[1] = 1.0;
1910 // Get pointer to local node
1911 local_node_pt = node_pt(i0 + n_p * (n_p - 1));
1912 break;
1913
1914 case 2:
1915 // Local fraction of node
1916 s_fraction[0] = 0.0;
1917 s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
1918 // Get pointer to local node
1919 local_node_pt = node_pt(n_p * i0);
1920 break;
1921
1922 case 3:
1923 // Local fraction of node
1924 s_fraction[0] = 1.0;
1925 s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
1926 // Get pointer to local node
1927 local_node_pt = node_pt(n_p - 1 + n_p * i0);
1928 break;
1929 }
1930
1931 // Calculate the local coordinate and the local coordinate as viewed
1932 // from the neighbour
1933 Vector<double> s_in_neighb(2);
1934 for (unsigned i = 0; i < 2; i++)
1935 {
1936 // Local coordinate in this element
1937 s[i] = -1.0 + 2.0 * s_fraction[i];
1938 // Local coordinate in the neighbour
1939 s_in_neighb[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
1940 (s_hi_neigh[i] - s_lo_neigh[i]);
1941 }
1942
1943 // Loop over timesteps
1944 for (unsigned t = 0; t < n_time; t++)
1945 {
1946 // Get the nodal position from neighbour element
1947 Vector<double> x_in_neighb(2);
1948 neigh_pt->object_pt()->interpolated_x(t, s_in_neighb, x_in_neighb);
1949
1950 // Check error only if the node is NOT periodic
1951 if (is_periodic == false)
1952 {
1953 for (int i = 0; i < 2; i++)
1954 {
1955 // Find the spatial error
1956 double err = std::fabs(local_node_pt->x(t, i) - x_in_neighb[i]);
1957
1958 // If it's bigger than our tolerance, say so
1959 if (err > 1e-9)
1960 {
1961 oomph_info << "errx " << err << " " << t << " "
1962 << local_node_pt->x(t, i) << " " << x_in_neighb[i]
1963 << std::endl;
1964
1965 oomph_info << "at " << local_node_pt->x(0) << " "
1966 << local_node_pt->x(1) << std::endl;
1967 }
1968
1969 // If it's bigger than the previous max error, it is the
1970 // new max error!
1971 if (err > max_error_x[i])
1972 {
1973 max_error_x[i] = err;
1974 }
1975 }
1976 }
1977
1978 // Get the values from neighbour element. Note: # of values
1979 // gets set by routine (because in general we don't know
1980 // how many interpolated values a certain element has
1981 Vector<double> values_in_neighb;
1983 t, s_in_neighb, values_in_neighb);
1984
1985 // Get the values in current element.
1986 Vector<double> values;
1987 get_interpolated_values(t, s, values);
1988
1989 // Now figure out how many continuously interpolated values there
1990 // are
1991 unsigned num_val =
1992 neigh_pt->object_pt()->ncont_interpolated_values();
1993
1994 // Check error
1995 for (unsigned ival = 0; ival < num_val; ival++)
1996 {
1997 double err = std::fabs(values[ival] - values_in_neighb[ival]);
1998
1999 if (err > 1.0e-10)
2000 {
2001 oomph_info << local_node_pt->x(0) << " " << local_node_pt->x(1)
2002 << " \n# "
2003 << "erru (S)" << err << " " << ival << " "
2004 << get_node_number(local_node_pt) << " "
2005 << values[ival] << " " << values_in_neighb[ival]
2006 << std::endl;
2007 }
2008
2009 if (err > max_error_val)
2010 {
2011 max_error_val = err;
2012 }
2013 }
2014 }
2015 }
2016 }
2017 }
2018
2019 max_error = max_error_x[0];
2020 if (max_error_x[1] > max_error) max_error = max_error_x[1];
2021 if (max_error_val > max_error) max_error = max_error_val;
2022
2023 if (max_error > 1e-9)
2024 {
2025 oomph_info << "\n#------------------------------------ \n#Max error ";
2026 oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
2027 << max_error_val << std::endl;
2028 oomph_info << "#------------------------------------ \n " << std::endl;
2029 }
2030 }
2031
2032
2033 //========================================================================
2034 /// Static matrix for coincidence between son nodal points and
2035 /// father boundaries
2036 ///
2037 //========================================================================
2038 std::map<unsigned, DenseMatrix<int>> RefineableQElement<2>::Father_bound;
2039
2040
2041 /// ///////////////////////////////////////////////////////////////////////
2042 /// ///////////////////////////////////////////////////////////////////////
2043 /// ///////////////////////////////////////////////////////////////////////
2044
2045
2046 //==================================================================
2047 /// Determine vector of solid (positional) boundary conditions along
2048 /// the element's boundary (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
2049 ///
2050 /// This function assumes that the same boundary condition is applied
2051 /// along the entire length of an element's edge (of course, the
2052 /// vertices combine the boundary conditions of their two adjacent edges
2053 /// in the most restrictive combination. Hence, if we're at a vertex,
2054 /// we apply the most restrictive boundary condition of the
2055 /// two adjacent edges. If we're on an edge (in its proper interior),
2056 /// we apply the least restrictive boundary condition of all nodes
2057 /// along the edge.
2058 ///
2059 /// Usual convention:
2060 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2061 /// on this boundary is free.
2062 /// - solid_bound_cons[i]=1 if it's pinned.
2063 //==================================================================
2065 int bound, Vector<int>& solid_bound_cons) const
2066 {
2067 using namespace QuadTreeNames;
2068
2069 // Spatial dimension of all nodes
2070 unsigned n_dim = this->nodal_dimension();
2071
2072 // Set up temporary vectors to hold edge boundary conditions
2073 Vector<int> bound_cons1(n_dim), bound_cons2(n_dim);
2074
2075 // Which boundary are we on?
2076 switch (bound)
2077 {
2078 // If on edge, just get the bcs
2079 case N:
2080 case S:
2081 case W:
2082 case E:
2083 get_edge_solid_bcs(bound, solid_bound_cons);
2084 break;
2085
2086 // Most restrictive boundary at SE corner
2087 case SE:
2088 get_edge_solid_bcs(S, bound_cons1);
2089 get_edge_solid_bcs(E, bound_cons2);
2090
2091 for (unsigned k = 0; k < n_dim; k++)
2092 {
2093 solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2094 }
2095 break;
2096
2097 // Most restrictive boundary at SW corner
2098 case SW:
2099 get_edge_solid_bcs(S, bound_cons1);
2100 get_edge_solid_bcs(W, bound_cons2);
2101
2102 for (unsigned k = 0; k < n_dim; k++)
2103 {
2104 solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2105 }
2106 break;
2107
2108 // Most restrictive boundary at NW corner
2109 case NW:
2110 get_edge_solid_bcs(N, bound_cons1);
2111 get_edge_solid_bcs(W, bound_cons2);
2112
2113 for (unsigned k = 0; k < n_dim; k++)
2114 {
2115 solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2116 }
2117 break;
2118
2119 // Most restrictive boundary at NE corner
2120 case NE:
2121 get_edge_solid_bcs(N, bound_cons1);
2122 get_edge_solid_bcs(E, bound_cons2);
2123
2124 for (unsigned k = 0; k < n_dim; k++)
2125 {
2126 solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2127 }
2128 break;
2129
2130 default:
2131 throw OomphLibError(
2132 "Wrong boundary", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2133 }
2134 }
2135
2136 //==================================================================
2137 /// Determine Vector of solid (positional) boundary conditions along
2138 /// the element's edge (S/N/W/E) -- BC is the least restrictive combination
2139 /// of all the nodes on this edge
2140 ///
2141 /// Usual convention:
2142 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2143 /// on this boundary is free
2144 /// - solid_bound_cons[i]=1 if it's pinned
2145 //==================================================================
2147 const int& edge, Vector<int>& solid_bound_cons) const
2148 {
2149 using namespace QuadTreeNames;
2150
2151 // Number of nodes along 1D edge
2152 unsigned n_p = nnode_1d();
2153
2154 // Left- and Right-hand nodes
2155 unsigned left_node, right_node;
2156
2157 // Set the left (lower) and right (upper) nodes for the edge
2158 switch (edge)
2159 {
2160 case N:
2161 left_node = n_p * (n_p - 1);
2162 right_node = n_p * n_p - 1;
2163 break;
2164
2165 case S:
2166 left_node = 0;
2167 right_node = n_p - 1;
2168 break;
2169
2170 case W:
2171 left_node = 0;
2172 right_node = n_p * (n_p - 1);
2173 break;
2174
2175 case E:
2176 left_node = n_p - 1;
2177 right_node = n_p * n_p - 1;
2178 break;
2179
2180 default:
2181 std::ostringstream error_stream;
2182 error_stream << "Wrong edge " << edge
2183 << " passed to get_solid_edge_bcs(..)" << std::endl;
2184
2185 throw OomphLibError(
2186 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2187 }
2188
2189
2190 // Cast to SolidNodes
2191 SolidNode* left_node_pt = dynamic_cast<SolidNode*>(node_pt(left_node));
2192 SolidNode* right_node_pt = dynamic_cast<SolidNode*>(node_pt(right_node));
2193#ifdef PARANOID
2194 if (left_node_pt == 0)
2195 {
2196 throw OomphLibError(
2197 "Left node cannot be cast to SolidNode --> something is wrong",
2198 OOMPH_CURRENT_FUNCTION,
2199 OOMPH_EXCEPTION_LOCATION);
2200 }
2201 if (right_node_pt == 0)
2202 {
2203 throw OomphLibError(
2204 "Right node cannot be cast to SolidNode --> something is wrong",
2205 OOMPH_CURRENT_FUNCTION,
2206 OOMPH_EXCEPTION_LOCATION);
2207 }
2208#endif
2209
2210
2211 // Number of coordinate directions
2212 unsigned n_dim = this->nodal_dimension();
2213
2214 // Loop over all directions, the least restrictive value is
2215 // the boundary condition at the left multiplied by that at the right
2216 // Assuming that free is always zero and pinned is one
2217 for (unsigned k = 0; k < n_dim; k++)
2218 {
2219 solid_bound_cons[k] = left_node_pt->position_is_pinned(k) *
2220 right_node_pt->position_is_pinned(k);
2221 }
2222 }
2223
2224 /// Build required templates
2225 // template class RefineableQElement<2>;
2226
2227} // namespace oomph
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
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...
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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 s_min() const
Min value of local coordinate.
Definition: elements.h:2793
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
virtual double s_max() const
Max. value of local coordinate.
Definition: elements.h:2803
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
Class that contains data for hanging nodes.
Definition: nodes.h:742
void set_master_node_pt(const unsigned &i, Node *const &master_node_pt, const double &weight)
Set the pointer to the i-th master node and its weight.
Definition: nodes.cc:1474
////////////////////////////////////////////////////////////////////
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_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:611
bool boundary_coordinate_exists(const unsigned &i) const
Indicate whether the i-th boundary has an intrinsic coordinate.
Definition: mesh.h:565
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
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
void set_hanging_pt(HangInfo *const &hang_pt, const int &i)
Set the hanging data for the i-th value. (hang_pt=0 to make non-hanging)
Definition: nodes.cc:2068
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
virtual void make_periodic(Node *const &node_pt)
Make the node periodic by copying the values from node_pt. Note that the coordinates will always rema...
Definition: nodes.cc:2257
unsigned nposition_type() const
Number of coordinate types needed in the mapping between local and global coordinates.
Definition: nodes.h:1016
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....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
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
QuadTree class: Recursively defined, generalised quadtree.
Definition: quadtree.h:104
QuadTree * gteq_edge_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, 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: quadtree.cc:413
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
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 Node * interpolating_node_pt(const unsigned &n, const int &value_id)
In mixed elements, different sets of nodes are used to interpolate different unknowns....
virtual unsigned ninterpolating_node(const int &value_id)
Return the number of nodes that are used to interpolate the value_id-th unknown. Default is to assume...
virtual void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
Return the basis functions that are used to interpolate the value_id-th unknown. By default assume is...
virtual Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
Return a pointer to the node that interpolates the value-id-th unknown at local coordinate s....
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<2,NNODE_1D>.
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Determine set of (mesh) boundaries that the element edge/vertex lives on.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along edge (or on vertex) bound (S/W/N/E/SW/SE/NW/NE): For va...
void interpolated_zeta_on_edge(const unsigned &boundary, const int &edge, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the edge (S/W/N/E)
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of Solid quad elements.
void get_solid_bcs(int bound, Vector< int > &solid_bound_cons) const
Determine vector of solid (positional) boundary conditions along edge (or on vertex) bound (S/W/N/E/S...
A class that is used to template the solid refineable Q elements by dimension. It's really nothing mo...
Definition: Qelements.h:2286
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
A Class for nodes that deform elastically (i.e. position is an unknown in the problem)....
Definition: nodes.h:1686
void pin_position(const unsigned &i)
Pin the nodal position.
Definition: nodes.h:1816
bool position_is_pinned(const unsigned &i)
Test whether the i-th coordinate is pinned, 0: false; 1: true.
Definition: nodes.h:1803
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void start(const unsigned &i)
(Re-)start i-th timer
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
unsigned Number
The unsigned.
void shape(const double &s, double *Psi)
Definition for 1D Lagrange shape functions. The value of all the shape functions at the local coordin...
Definition: shape.h:564
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...