refineable_brick_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
33
34namespace oomph
35{
36 //========================================================================
37 /// Print corner nodes, use colour (default "BLACK")
38 /// in right order so that tecplot can draw a cube without crossed lines
39 //========================================================================
40 void RefineableQElement<3>::output_corners(std::ostream& outfile,
41 const std::string& colour) const
42 {
44 Vector<double> corner(3);
45
46 outfile << "ZONE I=2, J=2, K=2 C=" << colour << std::endl;
47
48 s[0] = -1.0;
49 s[1] = -1.0;
50 s[2] = -1.0;
51 get_x(s, corner);
52 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
53 << Number << std::endl;
54
55 s[0] = -1.0;
56 s[1] = -1.0;
57 s[2] = 1.0;
58 get_x(s, corner);
59 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
60 << Number << std::endl;
61
62 s[0] = -1.0;
63 s[1] = 1.0;
64 s[2] = -1.0;
65 get_x(s, corner);
66 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
67 << Number << std::endl;
68
69 s[0] = -1.0;
70 s[1] = 1.0;
71 s[2] = 1.0;
72 get_x(s, corner);
73 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
74 << Number << std::endl;
75
76 // next face
77
78
79 s[0] = 1.0;
80 s[1] = -1.0;
81 s[2] = -1.0;
82 get_x(s, corner);
83 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
84 << Number << std::endl;
85
86 s[0] = 1.0;
87 s[1] = -1.0;
88 s[2] = 1.0;
89 get_x(s, corner);
90 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
91 << Number << std::endl;
92
93 s[0] = 1.0;
94 s[1] = 1.0;
95 s[2] = -1.0;
96 get_x(s, corner);
97 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
98 << Number << std::endl;
99
100 s[0] = 1.0;
101 s[1] = 1.0;
102 s[2] = 1.0;
103 get_x(s, corner);
104 outfile << corner[0] << " " << corner[1] << " " << corner[2] << " "
105 << Number << std::endl;
106
107
108 // outfile << "TEXT CS = GRID, X = " << corner[0] <<
109 // ",Y = " << corner[1] << ",Z = " << corner[2] <<
110 // ", HU = GRID, H = 0.01, AN = MIDCENTER, T =\""
111 // << Number << "\"" << std::endl;
112 }
113
114
115 //==================================================================
116 /// Setup static matrix for coincidence between son nodal points and
117 /// father boundaries:
118 ///
119 /// Father_bound[nnode_1d](nnode_son,son_type)={RU/RF/RD/RB/.../ OMEGA}
120 ///
121 /// so that node nnode_son in element of type son_type lies
122 /// on boundary/vertex Father_bound[nnode_1d](nnode_son,son_type) in its
123 /// father element. If the node doesn't lie on a boundary
124 /// the value is OMEGA.
125 //==================================================================
127 {
128 using namespace OcTreeNames;
129
130 // Find the number of nodes along a 1D edge
131 unsigned n_p = nnode_1d();
132
133 // Allocate space for the boundary information
134 Father_bound[n_p].resize(n_p * n_p * n_p, 8);
135
136 // Initialise: By default points are not on the boundary
137 for (unsigned n = 0; n < n_p * n_p * n_p; n++)
138 {
139 for (unsigned ison = 0; ison < 8; ison++)
140 {
141 Father_bound[n_p](n, ison) = Tree::OMEGA;
142 }
143 }
144
145 for (int i_son = LDB; i_son <= RUF; i_son++)
146 {
147 // vector representing the son
148 Vector<int> vect_son(3);
149 // vector representing (at the end) the boundaries
150 Vector<int> vect_bound(3);
151 vect_son = OcTree::Direction_to_vector[i_son];
152 for (unsigned i0 = 0; i0 < n_p; i0++)
153 {
154 for (unsigned i1 = 0; i1 < n_p; i1++)
155 {
156 for (unsigned i2 = 0; i2 < n_p; i2++)
157 {
158 // Initialisation to make it work
159 for (unsigned i = 0; i < 3; i++)
160 {
161 vect_bound[i] = -vect_son[i];
162 }
163 // Seting up the boundaries coordinates as if the coordinates
164 // of vect_bound had been initialised to 0 if it were so,
165 // vect_bound would be the vector of the boundaries in the son
166 // itself.
167
168 if (i0 == 0)
169 {
170 vect_bound[0] = -1;
171 }
172 if (i0 == n_p - 1)
173 {
174 vect_bound[0] = 1;
175 }
176 if (i1 == 0)
177 {
178 vect_bound[1] = -1;
179 }
180 if (i1 == n_p - 1)
181 {
182 vect_bound[1] = 1;
183 }
184 if (i2 == 0)
185 {
186 vect_bound[2] = -1;
187 }
188 if (i2 == n_p - 1)
189 {
190 vect_bound[2] = 1;
191 }
192
193 // The effect of this is to filter the boundaries to keep only the
194 // ones which are effectively father boundaries.
195 // -- if the node is not on a "i0 boundary", we still
196 // have vect_bound[0]=-vect_son[0]
197 // and the result is vect_bound[0]=0
198 // (he is not on this boundary for his father)
199 // -- if he is on a son's boundary which is not one of
200 // the father -> same thing
201 // -- if he is on a boundary which is one of his father's,
202 // vect_bound[i]=vect_son[i]
203 // and the new vect_bound[i] is the same as the old one
204 for (int i = 0; i < 3; i++)
205 {
206 vect_bound[i] = (vect_bound[i] + vect_son[i]) / 2;
207 }
208
209 // Return the result as {U,R,D,...RDB,LUF,OMEGA}
210 Father_bound[n_p](i0 + n_p * i1 + n_p * n_p * i2, i_son) =
211 OcTree::Vector_to_direction[vect_bound];
212
213 } // Loop over i2
214 } // Loop over i1
215 } // Loop over i0
216 } // Loop over i_son
217
218 } // setup_father_bounds()
219
220
221 //==================================================================
222 /// Determine Vector of boundary conditions along the element's boundary
223 /// bound.
224 ///
225 /// This function assumes that the same boundary condition is applied
226 /// along the entire area of an element's face (of course, the
227 /// vertices combine the boundary conditions of their two adjacent edges
228 /// in the most restrictive combination. Hence, if we're at a vertex,
229 /// we apply the most restrictive boundary condition of the
230 /// two adjacent edges. If we're on an edge (in its proper interior),
231 /// we apply the least restrictive boundary condition of all nodes
232 /// along the edge.
233 ///
234 /// Usual convention:
235 /// - bound_cons[ival]=0 if value ival on this boundary is free
236 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
237 //==================================================================
238 void RefineableQElement<3>::get_bcs(int bound, Vector<int>& bound_cons) const
239 {
240 using namespace OcTreeNames;
241
242 // Max. number of nodal data values in element
243 unsigned nvalue = ncont_interpolated_values();
244 // Set up temporary vectors to hold edge boundary conditions
245 Vector<int> bound_cons1(nvalue), bound_cons2(nvalue);
246 Vector<int> bound_cons3(nvalue);
247
248 Vector<int> vect1(3), vect2(3), vect3(3);
249 Vector<int> vect_elem;
250 Vector<int> notzero;
251 int n = 0;
252
253 vect_elem = OcTree::Direction_to_vector[bound];
254
255 // Just to see if bound is a face, an edge, or a vertex, n stores
256 // the number of non-zero values in the vector reprensenting the bound
257 // and the vector notzero stores the position of these values
258 for (int i = 0; i < 3; i++)
259 {
260 if (vect_elem[i] != 0)
261 {
262 n++;
263 notzero.push_back(i);
264 }
265 }
266
267 switch (n)
268 {
269 // If there is only one non-zero value, bound is a face
270 case 1:
271 get_face_bcs(bound, bound_cons);
272 break;
273
274 // If there are two non-zero values, bound is an edge
275 case 2:
276
277 for (unsigned i = 0; i < 3; i++)
278 {
279 vect1[i] = 0;
280 vect2[i] = 0;
281 }
282 // vect1 and vect2 are the vector of the two faces adjacent to bound
283 vect1[notzero[0]] = vect_elem[notzero[0]];
284 vect2[notzero[1]] = vect_elem[notzero[1]];
285
286 get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
287 get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
288 // get the most restrictive bc
289 for (unsigned k = 0; k < nvalue; k++)
290 {
291 bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
292 }
293 break;
294
295 // If there are three non-zero value, bound is a vertex
296 case 3:
297
298 for (unsigned i = 0; i < 3; i++)
299 {
300 vect1[i] = 0;
301 vect2[i] = 0;
302 vect3[i] = 0;
303 }
304 // vectors to the three adjacent faces of the vertex
305 vect1[0] = vect_elem[0];
306 vect2[1] = vect_elem[1];
307 vect3[2] = vect_elem[2];
308
309 get_face_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
310 get_face_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
311 get_face_bcs(OcTree::Vector_to_direction[vect3], bound_cons3);
312
313
314 // set the bcs to the most restrictive ones
315 for (unsigned k = 0; k < nvalue; k++)
316 {
317 bound_cons[k] = (bound_cons1[k] || bound_cons2[k] || bound_cons3[k]);
318 }
319 break;
320
321 default:
322 throw OomphLibError("Make sure you are not giving OMEGA as bound",
323 OOMPH_CURRENT_FUNCTION,
324 OOMPH_EXCEPTION_LOCATION);
325 }
326 }
327
328 //==================================================================
329 /// Determine Vector of boundary conditions along the element's
330 /// face (R/L/U/D/B/F) -- BC is the least restrictive combination
331 /// of all the nodes on this face
332 ///
333 /// Usual convention:
334 /// - bound_cons[ival]=0 if value ival on this boundary is free
335 /// - bound_cons[ival]=1 if value ival on this boundary is pinned
336 //==================================================================
338 Vector<int>& bound_cons) const
339 {
340 using namespace OcTreeNames;
341
342 // Number of nodes along 1D edge
343 unsigned n_p = nnode_1d();
344 // the four corner nodes on the boundary
345 unsigned node1, node2, node3, node4;
346
347 // Set the four corner nodes for the face
348 switch (face)
349 {
350 case U:
351 node1 = n_p * n_p * n_p - 1;
352 node2 = n_p * n_p - 1;
353 node3 = n_p * (n_p - 1);
354 node4 = n_p * (n_p * n_p - 1);
355 break;
356
357 case D:
358 node1 = 0;
359 node2 = n_p - 1;
360 node3 = (n_p * n_p + 1) * (n_p - 1);
361 node4 = n_p * n_p * (n_p - 1);
362 break;
363
364 case R:
365 node1 = n_p - 1;
366 node2 = (n_p * n_p + 1) * (n_p - 1);
367 node3 = n_p * n_p * n_p - 1;
368 node4 = n_p * n_p - 1;
369 break;
370
371 case L:
372 node1 = 0;
373 node2 = n_p * (n_p - 1);
374 node3 = n_p * (n_p * n_p - 1);
375 node4 = n_p * n_p * (n_p - 1);
376 break;
377
378 case B:
379 node1 = 0;
380 node2 = n_p - 1;
381 node3 = n_p * n_p - 1;
382 node4 = n_p * (n_p - 1);
383 break;
384
385 case F:
386 node1 = n_p * n_p * n_p - 1;
387 node2 = n_p * (n_p * n_p - 1);
388 node3 = n_p * n_p * (n_p - 1);
389 node4 = (n_p - 1) * (n_p * n_p + 1);
390 break;
391
392 default:
393 std::ostringstream error_stream;
394 error_stream << "Wrong edge " << face << " passed\n";
395
396 throw OomphLibError(
397 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
398 }
399
400 // Max. number of nodal data values in element
401 unsigned maxnvalue = ncont_interpolated_values();
402
403 // Loop over all values, the least restrictive value is
404 // the multiplication of the boundary conditions at the 4 nodes
405 // Assuming that free is always zero and pinned is one
406 for (unsigned k = 0; k < maxnvalue; k++)
407 {
408 bound_cons[k] =
409 node_pt(node1)->is_pinned(k) * node_pt(node2)->is_pinned(k) *
410 node_pt(node3)->is_pinned(k) * node_pt(node4)->is_pinned(k);
411 }
412 }
413
414
415 //==================================================================
416 /// Given an element edge/vertex, return a Vector which contains
417 /// all the (mesh-)boundary numbers that this element edge/vertex
418 /// lives on.
419 ///
420 /// For proper edges, the boundary is the one (if any) that is shared by
421 /// both vertex nodes). For vertex nodes, we just return their
422 /// boundaries.
423 //==================================================================
425 std::set<unsigned>& boundary) const
426 {
427 using namespace OcTreeNames;
428
429 // Number of 1d nodes along an edge
430 unsigned n_p = nnode_1d();
431 // Left and right-hand nodes
432 int node[4];
433 int a = 0, b = 0;
434 int n = 0;
435 Vector<int> vect_face(3);
436 vect_face = OcTree::Direction_to_vector[element];
437 // Set the left (lower) and right (upper) nodes for the edge
438
439 // this is to know what is the type of element (face, edge, vertex)
440 // we just need to count the number of values equal to 0 in the
441 // vector representing this element
442
443 // Local storage for the node-numbers in given directions
444 // Initialise to zero (assume at LH end of each domain)
445 int one_d_node_number[3] = {0, 0, 0};
446 // n stores the number of values equal to 0, a is the position of the
447 // last 0-value ;b is the position of the last non0-value
448 for (int i = 0; i < 3; i++)
449 {
450 if (vect_face[i] == 0)
451 {
452 a = i;
453 n++;
454 }
455 else
456 {
457 b = i;
458 // If we are at the extreme (RH) end of the face,
459 // set the node number accordingly
460 if (vect_face[i] == 1)
461 {
462 one_d_node_number[i] = n_p - 1;
463 }
464 }
465 }
466
467 switch (n)
468 {
469 // if n=0 element is a vertex, and need to look at only one node
470 case 0:
471 node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
472 n_p * n_p * one_d_node_number[2];
473 node[1] = node[0];
474 node[2] = node[0];
475 node[3] = node[0];
476 break;
477
478 // if n=1 element is an edge, and we need to look at two nodes
479 case 1:
480 if (a == 0)
481 {
482 node[0] = (n_p - 1) + n_p * one_d_node_number[1] +
483 n_p * n_p * one_d_node_number[2];
484 node[1] =
485 n_p * one_d_node_number[1] + n_p * n_p * one_d_node_number[2];
486 }
487 else if (a == 1)
488 {
489 node[0] = n_p * (n_p - 1) + one_d_node_number[0] +
490 n_p * n_p * one_d_node_number[2];
491 node[1] = one_d_node_number[0] + n_p * n_p * one_d_node_number[2];
492 }
493 else if (a == 2)
494 {
495 node[0] = one_d_node_number[0] + n_p * one_d_node_number[1] +
496 n_p * n_p * (n_p - 1);
497 node[1] = one_d_node_number[0] + n_p * one_d_node_number[1];
498 }
499 node[2] = node[1];
500 node[3] = node[1];
501 break;
502
503 // if n=2 element is a face, and we need to look at its 4 nodes
504 case 2:
505 if (b == 0)
506 {
507 node[0] =
508 one_d_node_number[0] + n_p * n_p * (n_p - 1) + n_p * (n_p - 1);
509 node[1] = one_d_node_number[0] + n_p * (n_p - 1);
510 node[2] = one_d_node_number[0] + n_p * n_p * (n_p - 1);
511 node[3] = one_d_node_number[0];
512 }
513 else if (b == 1)
514 {
515 node[0] =
516 n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1) + (n_p - 1);
517 node[1] = n_p * one_d_node_number[1] + (n_p - 1);
518 node[2] = n_p * one_d_node_number[1] + n_p * n_p * (n_p - 1);
519 node[3] = n_p * one_d_node_number[1];
520 }
521 else if (b == 2)
522 {
523 node[0] =
524 n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1) + (n_p - 1);
525 node[1] = n_p * n_p * one_d_node_number[2] + (n_p - 1);
526 node[2] = n_p * n_p * one_d_node_number[2] + n_p * (n_p - 1);
527 node[3] = n_p * n_p * one_d_node_number[2];
528 }
529 break;
530 default:
531 throw OomphLibError("Make sure you are not giving OMEGA as boundary",
532 OOMPH_CURRENT_FUNCTION,
533 OOMPH_EXCEPTION_LOCATION);
534 }
535
536
537 // Empty boundary set: Edge does not live on any boundary
538 boundary.clear();
539
540 // Storage for the boundaries at the four nodes
541 Vector<std::set<unsigned>*> node_boundaries_pt(4, 0);
542
543 // Loop over the four nodes and get the boundary information
544 for (unsigned i = 0; i < 4; i++)
545 {
546 node_pt(node[i])->get_boundaries_pt(node_boundaries_pt[i]);
547 }
548
549
550 // Now work out the intersections
551 Vector<std::set<unsigned>> boundary_aux(2);
552
553 for (unsigned i = 0; i < 2; i++)
554 {
555 // If the two nodes both lie on boundaries
556 if ((node_boundaries_pt[2 * i] != 0) &&
557 (node_boundaries_pt[2 * i + 1] != 0))
558 {
559 // Find the intersection (the common boundaries) of these nodes
560 std::set_intersection(
561 node_boundaries_pt[2 * i]->begin(),
562 node_boundaries_pt[2 * i]->end(),
563 node_boundaries_pt[2 * i + 1]->begin(),
564 node_boundaries_pt[2 * i + 1]->end(),
565 inserter(boundary_aux[i], boundary_aux[i].begin()));
566 }
567 }
568
569 // Now calculate the total intersection
570 set_intersection(boundary_aux[0].begin(),
571 boundary_aux[0].end(),
572 boundary_aux[1].begin(),
573 boundary_aux[1].end(),
574 inserter(boundary, boundary.begin()));
575 }
576
577
578 //===================================================================
579 /// Return the value of the intrinsic boundary coordinate interpolated
580 /// along the face
581 //===================================================================
583 const unsigned& boundary,
584 const int& face,
585 const Vector<double>& s,
586 Vector<double>& zeta)
587 {
588 using namespace OcTreeNames;
589
590 // Number of nodes along an edge
591 unsigned nnodes_1d = nnode_1d();
592
593 // Number of nodes on a face
594 unsigned nnodes_2d = nnodes_1d * nnodes_1d;
595
596 // Total number of nodes
597 unsigned nnodes_3d = nnode();
598
599 // Storage for the shape functions
600 Shape psi(nnodes_3d);
601
602 // Get the shape functions at the passed position
603 this->shape(s, psi);
604
605 // Unsigned data that give starts and increments for the loop
606 // over the nodes on the faces.
607 unsigned start = 0, increment1 = 1, increment2 = 1;
608
609 // Flag to record if actually on a face or an edge
610 bool on_edge = true;
611
612 // Which face?
613 switch (face)
614 {
615 case L:
616#ifdef PARANOID
617 if (s[0] != -1.0)
618 {
619 std::ostringstream error_stream;
620 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
621 << " is not on Left face\n";
622
623 throw OomphLibError(error_stream.str(),
624 OOMPH_CURRENT_FUNCTION,
625 OOMPH_EXCEPTION_LOCATION);
626 }
627#endif
628 // Start is zero (bottom-left-back corner)
629 increment1 = nnodes_1d;
630 increment2 = 0;
631 on_edge = false;
632 break;
633
634 case R:
635#ifdef PARANOID
636 if (s[0] != 1.0)
637 {
638 std::ostringstream error_stream;
639 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
640 << " is not on Right face\n";
641
642 throw OomphLibError(error_stream.str(),
643 OOMPH_CURRENT_FUNCTION,
644 OOMPH_EXCEPTION_LOCATION);
645 }
646#endif
647 // Start is bottom-right-back corner
648 start = nnodes_1d - 1;
649 increment1 = nnodes_1d;
650 increment2 = 0;
651 on_edge = false;
652 break;
653
654 case D:
655#ifdef PARANOID
656 if (s[1] != -1.0)
657 {
658 std::ostringstream error_stream;
659 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
660 << " is not on Bottom face\n";
661
662 throw OomphLibError(error_stream.str(),
663 OOMPH_CURRENT_FUNCTION,
664 OOMPH_EXCEPTION_LOCATION);
665 }
666#endif
667 // Start is zero and increments2 is nnode_2d-nnode_1d
668 increment2 = nnodes_2d - nnodes_1d;
669 on_edge = false;
670 break;
671
672 case U:
673#ifdef PARANOID
674 if (s[1] != 1.0)
675 {
676 std::ostringstream error_stream;
677 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
678 << " is not on Upper face\n";
679
680 throw OomphLibError(error_stream.str(),
681 OOMPH_CURRENT_FUNCTION,
682 OOMPH_EXCEPTION_LOCATION);
683 }
684#endif
685 // Start is top-left-back corner and increments2 is nnode_2d-nnode_1d
686 start = nnodes_2d - nnodes_1d;
687 increment2 = start;
688 on_edge = false;
689 break;
690
691 case B:
692#ifdef PARANOID
693 if (s[2] != -1.0)
694 {
695 std::ostringstream error_stream;
696 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
697 << " is not on Back face\n";
698
699 throw OomphLibError(error_stream.str(),
700 OOMPH_CURRENT_FUNCTION,
701 OOMPH_EXCEPTION_LOCATION);
702 }
703#endif
704 // Start is zero and increments are 1 and 0
705 increment2 = 0;
706 on_edge = false;
707 break;
708
709 case F:
710#ifdef PARANOID
711 if (s[2] != 1.0)
712 {
713 std::ostringstream error_stream;
714 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
715 << " is not on Front face\n";
716
717 throw OomphLibError(error_stream.str(),
718 OOMPH_CURRENT_FUNCTION,
719 OOMPH_EXCEPTION_LOCATION);
720 }
721#endif
722 // Start is bottom-left-front corner
723 start = nnodes_3d - nnodes_2d;
724 increment2 = 0;
725 on_edge = false;
726 break;
727
728 case LF:
729#ifdef PARANOID
730 if ((s[0] != -1.0) || (s[2] != 1.0))
731 {
732 std::ostringstream error_stream;
733 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
734 << " is not on Front-Left edge\n";
735
736 throw OomphLibError(error_stream.str(),
737 OOMPH_CURRENT_FUNCTION,
738 OOMPH_EXCEPTION_LOCATION);
739 }
740#endif
741 start = nnodes_3d - nnodes_2d;
742 increment1 = nnodes_1d;
743 break;
744
745 case LD:
746#ifdef PARANOID
747 if ((s[0] != -1.0) || (s[1] != -1.0))
748 {
749 std::ostringstream error_stream;
750 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
751 << " is not on Bottom-Left edge\n";
752
753 throw OomphLibError(error_stream.str(),
754 OOMPH_CURRENT_FUNCTION,
755 OOMPH_EXCEPTION_LOCATION);
756 }
757#endif
758 increment1 = nnodes_2d;
759 break;
760
761 case LU:
762#ifdef PARANOID
763 if ((s[0] != -1.0) || (s[1] != 1.0))
764 {
765 std::ostringstream error_stream;
766 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
767 << " is not on Upper-Left edge\n";
768
769 throw OomphLibError(error_stream.str(),
770 OOMPH_CURRENT_FUNCTION,
771 OOMPH_EXCEPTION_LOCATION);
772 }
773#endif
774 start = nnodes_2d - nnodes_1d;
775 increment1 = nnodes_2d;
776 break;
777
778 case LB:
779#ifdef PARANOID
780 if ((s[0] != -1.0) || (s[2] != -1.0))
781 {
782 std::ostringstream error_stream;
783 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
784 << " is not on Back-Left edge\n";
785
786 throw OomphLibError(error_stream.str(),
787 OOMPH_CURRENT_FUNCTION,
788 OOMPH_EXCEPTION_LOCATION);
789 }
790#endif
791 increment1 = nnodes_1d;
792 break;
793
794 case RF:
795#ifdef PARANOID
796 if ((s[0] != 1.0) || (s[2] != 1.0))
797 {
798 std::ostringstream error_stream;
799 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
800 << " is not on Front-Right edge\n";
801
802 throw OomphLibError(error_stream.str(),
803 OOMPH_CURRENT_FUNCTION,
804 OOMPH_EXCEPTION_LOCATION);
805 }
806#endif
807 start = nnodes_3d - 1;
808 increment1 = -nnodes_1d;
809 break;
810
811 case RD:
812#ifdef PARANOID
813 if ((s[0] != 1.0) || (s[1] != -1.0))
814 {
815 std::ostringstream error_stream;
816 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
817 << " is not on Bottom-Right edge\n";
818
819 throw OomphLibError(error_stream.str(),
820 OOMPH_CURRENT_FUNCTION,
821 OOMPH_EXCEPTION_LOCATION);
822 }
823#endif
824 start = nnodes_1d - 1;
825 increment1 = nnodes_2d;
826 break;
827
828 case RU:
829#ifdef PARANOID
830 if ((s[0] != 1.0) || (s[1] != 1.0))
831 {
832 std::ostringstream error_stream;
833 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
834 << " is not on Upper-Right edge\n";
835
836 throw OomphLibError(error_stream.str(),
837 OOMPH_CURRENT_FUNCTION,
838 OOMPH_EXCEPTION_LOCATION);
839 }
840#endif
841 start = nnodes_2d - 1;
842 increment1 = nnodes_2d;
843 break;
844
845 case RB:
846#ifdef PARANOID
847 if ((s[0] != 1.0) || (s[2] != -1.0))
848 {
849 std::ostringstream error_stream;
850 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
851 << " is not on Back-Right edge\n";
852
853 throw OomphLibError(error_stream.str(),
854 OOMPH_CURRENT_FUNCTION,
855 OOMPH_EXCEPTION_LOCATION);
856 }
857#endif
858 start = nnodes_1d - 1;
859 increment1 = nnodes_1d;
860 break;
861
862 case DB:
863#ifdef PARANOID
864 if ((s[1] != -1.0) || (s[2] != -1.0))
865 {
866 std::ostringstream error_stream;
867 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
868 << " is not on Back-Bottom edge\n";
869
870 throw OomphLibError(error_stream.str(),
871 OOMPH_CURRENT_FUNCTION,
872 OOMPH_EXCEPTION_LOCATION);
873 }
874#endif
875 break;
876
877 case DF:
878#ifdef PARANOID
879 if ((s[1] != -1.0) || (s[2] != 1.0))
880 {
881 std::ostringstream error_stream;
882 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
883 << " is not on Front-Bottom edge\n";
884
885 throw OomphLibError(error_stream.str(),
886 OOMPH_CURRENT_FUNCTION,
887 OOMPH_EXCEPTION_LOCATION);
888 }
889#endif
890 start = nnodes_3d - nnodes_2d;
891 break;
892
893 case UB:
894#ifdef PARANOID
895 if ((s[1] != 1.0) || (s[2] != -1.0))
896 {
897 std::ostringstream error_stream;
898 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
899 << " is not on Back-Upper edge\n";
900
901 throw OomphLibError(error_stream.str(),
902 OOMPH_CURRENT_FUNCTION,
903 OOMPH_EXCEPTION_LOCATION);
904 }
905#endif
906 start = nnodes_2d - nnodes_1d;
907 break;
908
909 case UF:
910#ifdef PARANOID
911 if ((s[1] != 1.0) || (s[2] != 1.0))
912 {
913 std::ostringstream error_stream;
914 error_stream << "Coordinate " << s[0] << " " << s[1] << " " << s[2]
915 << " is not on Upper-Front edge\n";
916
917 throw OomphLibError(error_stream.str(),
918 OOMPH_CURRENT_FUNCTION,
919 OOMPH_EXCEPTION_LOCATION);
920 }
921#endif
922 start = nnodes_3d - nnodes_1d;
923 break;
924
925 default:
926
927 std::ostringstream error_stream;
928 error_stream << "Wrong face " << OcTree::Direct_string[face]
929 << " passed" << std::endl;
930 error_stream << "Trouble at : s= [" << s[0] << " " << s[1] << " "
931 << s[2] << "]\n";
932 Vector<double> x(3);
933 this->interpolated_x(s, x);
934 error_stream << "corresponding to : x= [" << x[0] << " " << x[1] << " "
935 << x[2] << "]\n";
936 throw OomphLibError(
937 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
938 }
939
940 // Initialise the intrinsic coordinate
941 zeta[0] = 0.0;
942 zeta[1] = 0.0;
943
944 // Set the start node number
945 unsigned node = start;
946
947 if (on_edge)
948 {
949 for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
950 {
951 // Get the intrinsic coordinate
952 Vector<double> node_zeta(2);
953 node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
954
955 // Now multiply by the shape function
956 zeta[0] += node_zeta[0] * psi(node);
957 zeta[1] += node_zeta[1] * psi(node);
958
959 // Update node
960 node += increment1;
961 }
962 }
963 else
964 {
965 for (unsigned l2 = 0; l2 < nnodes_1d; l2++)
966 {
967 for (unsigned l1 = 0; l1 < nnodes_1d; l1++)
968 {
969 // Get the intrinsic coordinate
970 Vector<double> node_zeta(2);
971 node_pt(node)->get_coordinates_on_boundary(boundary, node_zeta);
972
973 // Now multiply by the shape function
974 zeta[0] += node_zeta[0] * psi(node);
975 zeta[1] += node_zeta[1] * psi(node);
976
977 // Update node
978 node += increment1;
979 }
980 // Update node
981 node += increment2;
982 }
983 }
984 }
985
986 //===================================================================
987 /// If a neighbouring element has already created a node at a
988 /// position corresponding to the local fractional position within the
989 /// present element, s_fraction, return
990 /// a pointer to that node. If not, return NULL (0).
991 //===================================================================
993 const Vector<double>& s_fraction, bool& is_periodic)
994 {
995 using namespace OcTreeNames;
996
997 // Calculate the faces/edges on which the node lies
998 Vector<int> faces;
999 Vector<int> edges;
1000
1001 if (s_fraction[0] == 0.0)
1002 {
1003 faces.push_back(L);
1004 if (s_fraction[1] == 0.0)
1005 {
1006 edges.push_back(LD);
1007 }
1008 if (s_fraction[2] == 0.0)
1009 {
1010 edges.push_back(LB);
1011 }
1012 if (s_fraction[1] == 1.0)
1013 {
1014 edges.push_back(LU);
1015 }
1016 if (s_fraction[2] == 1.0)
1017 {
1018 edges.push_back(LF);
1019 }
1020 }
1021
1022 if (s_fraction[0] == 1.0)
1023 {
1024 faces.push_back(R);
1025 if (s_fraction[1] == 0.0)
1026 {
1027 edges.push_back(RD);
1028 }
1029 if (s_fraction[2] == 0.0)
1030 {
1031 edges.push_back(RB);
1032 }
1033 if (s_fraction[1] == 1.0)
1034 {
1035 edges.push_back(RU);
1036 }
1037 if (s_fraction[2] == 1.0)
1038 {
1039 edges.push_back(RF);
1040 }
1041 }
1042
1043 if (s_fraction[1] == 0.0)
1044 {
1045 faces.push_back(D);
1046 if (s_fraction[2] == 0.0)
1047 {
1048 edges.push_back(DB);
1049 }
1050 if (s_fraction[2] == 1.0)
1051 {
1052 edges.push_back(DF);
1053 }
1054 }
1055
1056 if (s_fraction[1] == 1.0)
1057 {
1058 faces.push_back(U);
1059 if (s_fraction[2] == 0.0)
1060 {
1061 edges.push_back(UB);
1062 }
1063 if (s_fraction[2] == 1.0)
1064 {
1065 edges.push_back(UF);
1066 }
1067 }
1068
1069 if (s_fraction[2] == 0.0)
1070 {
1071 faces.push_back(B);
1072 }
1073
1074 if (s_fraction[2] == 1.0)
1075 {
1076 faces.push_back(F);
1077 }
1078
1079 // Find the number of faces
1080 unsigned n_face = faces.size();
1081
1082 // Find the number of edges
1083 unsigned n_edge = edges.size();
1084
1085 Vector<unsigned> translate_s(3);
1086 Vector<double> s_lo_neigh(3);
1087 Vector<double> s_hi_neigh(3);
1088 Vector<double> s(3);
1089
1090 int neigh_face, diff_level;
1091
1092 // Loop over the faces on which the node lies
1093 //-------------------------------------------
1094 for (unsigned j = 0; j < n_face; j++)
1095 {
1096 // Boolean to indicate whether or not the neighbour has a
1097 // different Tree root
1098 bool in_neighbouring_tree;
1099
1100 // Find pointer to neighbouring element along face
1101 OcTree* neigh_pt;
1102 neigh_pt = octree_pt()->gteq_face_neighbour(faces[j],
1103 translate_s,
1104 s_lo_neigh,
1105 s_hi_neigh,
1106 neigh_face,
1107 diff_level,
1108 in_neighbouring_tree);
1109
1110 // Neighbour exists
1111 if (neigh_pt != 0)
1112 {
1113 // Have its nodes been created yet?
1114 if (neigh_pt->object_pt()->nodes_built())
1115 {
1116 // We now need to translate the nodal location, defined in terms
1117 // of the fractional coordinates of the present element into
1118 // those of its neighbour. For this we use the information returned
1119 // to use from the octree function.
1120
1121 // Calculate the local coordinate in the neighbour
1122 // Note that we need to use the translation scheme in case
1123 // the local coordinates are swapped in the neighbour.
1124 for (unsigned i = 0; i < 3; i++)
1125 {
1126 s[i] = s_lo_neigh[i] +
1127 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
1128 }
1129
1130 // Find the node in the neighbour
1131 Node* neighbour_node_pt =
1133
1134 // If there is a node, return it
1135 if (neighbour_node_pt != 0)
1136 {
1137 // Now work out whether it's a periodic boundary. This is
1138 // only possible if we have moved into a neighbouring tree
1139 if (in_neighbouring_tree)
1140 {
1141 // Return whether the neighbour is periodic
1142 is_periodic =
1143 octree_pt()->root_pt()->is_neighbour_periodic(faces[j]);
1144 }
1145
1146 // Return the neighbour node pointer
1147 return neighbour_node_pt;
1148 }
1149 } // if (neigh_pt->object_pt()->nodes_built())
1150 } // if (neigh_pt!=0)
1151 } // for (unsigned j=0;j<n_face;j++)
1152
1153 // Loop over the edges on which the node lies
1154 //------------------------------------------
1155 for (unsigned j = 0; j < n_edge; j++)
1156 {
1157 // Even if we restrict ourselves to true edge neighbours (i.e.
1158 // elements that are not also face neighbours) there may be multiple
1159 // edge neighbours across the edges between multiple root octrees.
1160 // When making the first call to OcTree::gteq_true_edge_neighbour(...)
1161 // we simply return the first one of these multiple edge neighbours
1162 // (if there are any at all, of course) and also return the total number
1163 // of true edge neighbours. If the node in question already exists
1164 // on the first edge neighbour we're done. If it doesn't it may exist
1165 // on other edge neighbours so we repeat the process over all
1166 // other edge neighbours (bailing out if a node is found, of course).
1167
1168 // Initially return the zero-th true edge neighbour
1169 unsigned i_root_edge_neighbour = 0;
1170
1171 // Initialise the total number of true edge neighbours
1172 unsigned nroot_edge_neighbour = 0;
1173
1174 // Keep searching until we've found the node or until we've checked
1175 // all available edge neighbours
1176 bool keep_searching = true;
1177 while (keep_searching)
1178 {
1179 // Find pointer to neighbouring element along edge
1180 OcTree* neigh_pt;
1181 neigh_pt = octree_pt()->gteq_true_edge_neighbour(edges[j],
1182 i_root_edge_neighbour,
1183 nroot_edge_neighbour,
1184 translate_s,
1185 s_lo_neigh,
1186 s_hi_neigh,
1187 neigh_face,
1188 diff_level);
1189
1190 // Neighbour exists
1191 if (neigh_pt != 0)
1192 {
1193 // Have its nodes been created yet?
1194 if (neigh_pt->object_pt()->nodes_built())
1195 {
1196 // We now need to translate the nodal location, defined in terms
1197 // of the fractional coordinates of the present element into
1198 // those of its neighbour. For this we use the information returned
1199 // to use from the octree function.
1200
1201 // Calculate the local coordinate in the neighbour
1202 // Note that we need to use the translation scheme in case
1203 // the local coordinates are swapped in the neighbour.
1204 for (unsigned i = 0; i < 3; i++)
1205 {
1206 s[i] = s_lo_neigh[i] + s_fraction[translate_s[i]] *
1207 (s_hi_neigh[i] - s_lo_neigh[i]);
1208 }
1209
1210 // Find the node in the neighbour
1211 Node* neighbour_node_pt =
1213
1214 // If there is a node, return it
1215 if (neighbour_node_pt != 0)
1216 {
1217 // Get the faces on which the edge lies
1218 Vector<int> faces_attached_to_edge =
1220
1221 // Get the number of entries in the vector
1222 unsigned n_faces_attached_to_edge = faces_attached_to_edge.size();
1223
1224 // Loop over the faces
1225 for (unsigned i_face = 0; i_face < n_faces_attached_to_edge;
1226 i_face++)
1227 {
1228 // Is the node periodic in the face direction?
1229 is_periodic = octree_pt()->root_pt()->is_neighbour_periodic(
1230 faces_attached_to_edge[i_face]);
1231
1232 // Check if the edge is periodic in the i_face-th face direction
1233 if (is_periodic)
1234 {
1235 // We're done!
1236 break;
1237 }
1238 } // for (unsigned
1239 // i_face=0;i_face<n_faces_attached_to_edge;i_face++)
1240
1241 // Return the neighbour node pointer
1242 return neighbour_node_pt;
1243 } // if (neighbour_node_pt!=0)
1244 } // if (neigh_pt->object_pt()->nodes_built())
1245 } // if (neigh_pt!=0)
1246
1247 // Keep searching, but only if there are further edge neighbours
1248 // Try next root edge neighbour
1249 i_root_edge_neighbour++;
1250
1251 // Have we exhausted the search?
1252 if (i_root_edge_neighbour >= nroot_edge_neighbour)
1253 {
1254 // Stop searching
1255 keep_searching = false;
1256 }
1257 } // End of while keep searching over all true edge neighbours
1258 } // End of loop over edges
1259
1260 // Node not found, return null
1261 return 0;
1262 }
1263
1264
1265 //==================================================================
1266 /// Build the element by doing the following:
1267 /// - Give it nodal positions (by establishing the pointers to its
1268 /// nodes)
1269 /// - In the process create new nodes where required (i.e. if they
1270 /// don't exist in father element or have already been created
1271 /// while building new neighbour elements). Node building
1272 /// involves the following steps:
1273 /// - Get nodal position from father element.
1274 /// - Establish the time-history of the newly created nodal point
1275 /// (its coordinates and the previous values) consistent with
1276 /// the father's history.
1277 /// - Determine the boundary conditions of the nodes (newly
1278 /// created nodes can only lie on the interior of any
1279 /// edges of the father element -- this makes it possible to
1280 /// to figure out what their bc should be...)
1281 /// - Add node to the mesh's stoarge scheme for the boundary nodes.
1282 /// - Add the new node to the mesh itself
1283 /// - Doc newly created nodes in file "new_nodes.dat" in the directory
1284 /// / of the DocInfo object (only if it's open!)
1285 /// - Finally, excute the element-specific further_build()
1286 /// (empty by default -- must be overloaded for specific elements).
1287 /// This deals with any build operations that are not included
1288 /// in the generic process outlined above. For instance, in
1289 /// Crouzeix Raviart elements we need to initialise the internal
1290 /// pressure values in manner consistent with the pressure
1291 /// distribution in the father element.
1292 //==================================================================
1294 Vector<Node*>& new_node_pt,
1295 bool& was_already_built,
1296 std::ofstream& new_nodes_file)
1297 {
1298 using namespace OcTreeNames;
1299
1300 // Number of dimensions
1301 unsigned n_dim = 3;
1302
1303 // Get the number of 1d nodes
1304 unsigned n_p = nnode_1d();
1305
1306 // Check whether static father_bound needs to be created
1307 if (Father_bound[n_p].nrow() == 0)
1308 {
1309 setup_father_bounds();
1310 }
1311
1312 // Pointer to my father (in octree impersonation)
1313 OcTree* father_pt = dynamic_cast<OcTree*>(octree_pt()->father_pt());
1314
1315 // What type of son am I? Ask my octree representation...
1316 int son_type = octree_pt()->son_type();
1317
1318 // Has somebody build me already? (If any nodes have not been built)
1319 if (!nodes_built())
1320 {
1321#ifdef PARANOID
1322 if (father_pt == 0)
1323 {
1324 std::string error_message =
1325 "Something fishy here: I have no father and yet \n";
1326 error_message += "I have no nodes. Who has created me then?!\n";
1327
1328 throw OomphLibError(
1329 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1330 }
1331#endif
1332
1333 // Indicate status:
1334 was_already_built = false;
1335
1336 // Return pointer to father element
1337 RefineableQElement<3>* father_el_pt =
1338 dynamic_cast<RefineableQElement<3>*>(father_pt->object_pt());
1339
1340 // Timestepper should be the same for all nodes in father
1341 // element -- use it create timesteppers for new nodes
1342 TimeStepper* time_stepper_pt =
1343 father_el_pt->node_pt(0)->time_stepper_pt();
1344
1345 // Number of history values (incl. present)
1346 unsigned ntstorage = time_stepper_pt->ntstorage();
1347
1348 // Currently we can't handle the case of generalised coordinates
1349 // since we haven't established how they should be interpolated
1350 // Buffer this case:
1351 if (father_el_pt->node_pt(0)->nposition_type() != 1)
1352 {
1353 throw OomphLibError("Can't handle generalised nodal positions (yet).",
1354 OOMPH_CURRENT_FUNCTION,
1355 OOMPH_EXCEPTION_LOCATION);
1356 }
1357
1358 Vector<int> s_lo(n_dim);
1359 Vector<int> s_hi(n_dim);
1360 Vector<double> s(n_dim);
1361 Vector<double> x(n_dim);
1362
1363 // Setup vertex coordinates in father element:
1364 //--------------------------------------------
1365 // find the s_lo coordinates
1366 s_lo = octree_pt()->Direction_to_vector[son_type];
1367
1368 // Just scale them, because the Direction_to_vector
1369 // doesn't really gives s_lo;
1370 for (unsigned i = 0; i < n_dim; i++)
1371 {
1372 s_lo[i] = (s_lo[i] + 1) / 2 - 1;
1373 }
1374
1375 // setup s_hi (Actually s_hi[i]=s_lo[i]+1)
1376 for (unsigned i = 0; i < n_dim; i++)
1377 {
1378 s_hi[i] = s_lo[i] + 1;
1379 }
1380
1381 // Pass macro element pointer on to sons and
1382 // set coordinates in macro element
1383 if (father_el_pt->macro_elem_pt() != 0)
1384 {
1385 set_macro_elem_pt(father_el_pt->macro_elem_pt());
1386 for (unsigned i = 0; i < n_dim; i++)
1387 {
1388 s_macro_ll(i) =
1389 father_el_pt->s_macro_ll(i) +
1390 0.5 * (s_lo[i] + 1.0) *
1391 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1392 s_macro_ur(i) =
1393 father_el_pt->s_macro_ll(i) +
1394 0.5 * (s_hi[i] + 1.0) *
1395 (father_el_pt->s_macro_ur(i) - father_el_pt->s_macro_ll(i));
1396 }
1397 }
1398
1399
1400 // If the father element hasn't been generated yet, we're stuck...
1401 if (father_el_pt->node_pt(0) == 0)
1402 {
1403 throw OomphLibError(
1404 "Trouble: father_el_pt->node_pt(0)==0\n Can't build son element!\n",
1405 OOMPH_CURRENT_FUNCTION,
1406 OOMPH_EXCEPTION_LOCATION);
1407 }
1408 else
1409 {
1410 unsigned jnod = 0;
1411 Vector<double> s_fraction(n_dim);
1412
1413 // Loop over nodes in element
1414 for (unsigned i0 = 0; i0 < n_p; i0++)
1415 {
1416 // Get the fractional position of the node in the direction of s[0]
1417 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
1418 // Local coordinate in father element
1419 s[0] = s_lo[0] + (s_hi[0] - s_lo[0]) * s_fraction[0];
1420
1421 for (unsigned i1 = 0; i1 < n_p; i1++)
1422 {
1423 // Get the fractional position of the node in the direction of s[1]
1424 s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
1425 // Local coordinate in father element
1426 s[1] = s_lo[1] + (s_hi[1] - s_lo[1]) * s_fraction[1];
1427
1428 for (unsigned i2 = 0; i2 < n_p; i2++)
1429 {
1430 // Get the fractional position of the node in the direction of
1431 // s[2]
1432 s_fraction[2] = local_one_d_fraction_of_node(i2, 2);
1433
1434 // Local coordinate in father element
1435 s[2] = s_lo[2] + (s_hi[2] - s_lo[2]) * s_fraction[2];
1436
1437 // Local node number
1438 jnod = i0 + n_p * i1 + n_p * n_p * i2;
1439
1440 // Initialise flag: So far, this node hasn't been built
1441 // or copied yet
1442 bool node_done = false;
1443
1444 // Get the pointer to the node in the father; returns NULL
1445 // if there is not a node
1446 Node* created_node_pt =
1447 father_el_pt->get_node_at_local_coordinate(s);
1448
1449 // Does this node already exist in father element?
1450 //------------------------------------------------
1451 bool node_exists_in_father = false;
1452 if (created_node_pt != 0)
1453 {
1454 // Remember this!
1455 node_exists_in_father = true;
1456
1457 // Copy node across
1458 node_pt(jnod) = created_node_pt;
1459
1460 // Make sure that we update the values at the node so that
1461 // they are consistent with the present representation.
1462 // This is only need for mixed interpolation where the value
1463 // at the father could now become active.
1464
1465 // Loop over all history values
1466 for (unsigned t = 0; t < ntstorage; t++)
1467 {
1468 // Get values from father element
1469 // Note: get_interpolated_values() sets Vector size itself.
1470 Vector<double> prev_values;
1471 father_el_pt->get_interpolated_values(t, s, prev_values);
1472 // Find the minimum number of values
1473 //(either those stored at the node, or those returned by
1474 // the function)
1475 unsigned n_val_at_node = created_node_pt->nvalue();
1476 unsigned n_val_from_function = prev_values.size();
1477 // Use the ternary conditional operator here
1478 unsigned n_var = n_val_at_node < n_val_from_function ?
1479 n_val_at_node :
1480 n_val_from_function;
1481 // Assign the values that we can
1482 for (unsigned k = 0; k < n_var; k++)
1483 {
1484 created_node_pt->set_value(t, k, prev_values[k]);
1485 }
1486 }
1487
1488 // Node has been created by copy
1489 node_done = true;
1490 }
1491 // Node does not exist in father element but might already
1492 //--------------------------------------------------------
1493 // have been created by neighbouring elements
1494 //-------------------------------------------
1495 else
1496 {
1497 // Boolean to check if the node is periodic
1498 bool is_periodic = false;
1499
1500 // Was the node created by one of its neighbours
1501 // Whether or not the node lies on an edge can be determined
1502 // from the fractional position
1503 created_node_pt =
1504 node_created_by_neighbour(s_fraction, is_periodic);
1505
1506 // If so, then copy the pointer across
1507 if (created_node_pt != 0)
1508 {
1509 // Now the node must be on a boundary, but we don't know which
1510 // one. The returned created_node_pt is actually the
1511 // neighbouring periodic node
1512 Node* neighbour_node_pt = created_node_pt;
1513
1514 // Determine the edge on which the new node will live
1515 int father_bound = Father_bound[n_p](jnod, son_type);
1516
1517 // Storage for the set of Mesh boundaries on which the
1518 // appropriate father edge lives.
1519 std::set<unsigned> boundaries;
1520
1521 // Only get the boundaries if we are at the edge of
1522 // an element. Nodes in the centre of an element cannot be
1523 // on Mesh boundaries
1524 if (father_bound != Tree::OMEGA)
1525 {
1526 father_el_pt->get_boundaries(father_bound, boundaries);
1527 }
1528
1529#ifdef PARANOID
1530 // Case where a new node lives on more than one boundary
1531 // seems fishy enough to flag
1532 if (boundaries.size() > 2)
1533 {
1534 throw OomphLibError(
1535 "boundaries.size()>2 seems a bit strange..\n",
1536 OOMPH_CURRENT_FUNCTION,
1537 OOMPH_EXCEPTION_LOCATION);
1538 }
1539#endif
1540
1541 // If the node is periodic and is definitely a boundary node.
1542 // NOTE: the reason for this is that confusion can arise when
1543 // a node is created on an edge that joins a periodic face.
1544 if ((is_periodic) && (boundaries.size() > 0))
1545 {
1546 // Create node and set the pointer to it from the element
1547 created_node_pt =
1548 construct_boundary_node(jnod, time_stepper_pt);
1549
1550 // Make the node periodic from the neighbour
1551 created_node_pt->make_periodic(neighbour_node_pt);
1552
1553 // Add to vector of new nodes
1554 new_node_pt.push_back(created_node_pt);
1555
1556 // Loop over # of history values
1557 for (unsigned t = 0; t < ntstorage; t++)
1558 {
1559 // Get position from father element -- this uses the macro
1560 // element representation if appropriate. If the node
1561 // turns out to be a hanging node later on, then
1562 // its position gets adjusted in line with its
1563 // hanging node interpolation.
1564 Vector<double> x_prev(n_dim);
1565 father_el_pt->get_x(t, s, x_prev);
1566 // Set previous positions of the new node
1567 for (unsigned i = 0; i < n_dim; i++)
1568 {
1569 created_node_pt->x(t, i) = x_prev[i];
1570 }
1571 }
1572
1573 // Next, we Update the boundary lookup schemes
1574 // Loop over the boundaries stored in the set
1575 for (std::set<unsigned>::iterator it = boundaries.begin();
1576 it != boundaries.end();
1577 ++it)
1578 {
1579 // Add the node to the boundary
1580 mesh_pt->add_boundary_node(*it, created_node_pt);
1581
1582 // If we have set an intrinsic coordinate on this
1583 // mesh boundary then it must also be interpolated on
1584 // the new node
1585 // Now interpolate the intrinsic boundary coordinate
1586 if (mesh_pt->boundary_coordinate_exists(*it) == true)
1587 {
1588 Vector<double> zeta(2, 0.0);
1589 father_el_pt->interpolated_zeta_on_face(
1590 *it, father_bound, s, zeta);
1591 created_node_pt->set_coordinates_on_boundary(*it, zeta);
1592 }
1593 }
1594
1595 // Make sure that we add the node to the mesh
1596 mesh_pt->add_node_pt(created_node_pt);
1597 } // End of periodic case
1598 // Otherwise the node is not periodic, so just set the
1599 // pointer to the neighbours node
1600 else
1601 {
1602 node_pt(jnod) = created_node_pt;
1603 }
1604 node_done = true;
1605 }
1606 // Node does not exist in neighbour element but might already
1607 //-----------------------------------------------------------
1608 // have been created by a son of a neighbouring element
1609 //-----------------------------------------------------
1610 else
1611 {
1612 // Was the node created by one of its neighbours' sons
1613 // Whether or not the node lies on an edge can be calculated
1614 // by from the fractional position
1615 bool is_periodic = false;
1616 created_node_pt =
1617 node_created_by_son_of_neighbour(s_fraction, is_periodic);
1618
1619 // If the node was so created, assign the pointers
1620 if (created_node_pt != 0)
1621 {
1622 // If the node is periodic
1623 if (is_periodic)
1624 {
1625 // Now the node must be on a boundary, but we don't know
1626 // which one The returned created_node_pt is actually the
1627 // neighbouring periodic node
1628 Node* neighbour_node_pt = created_node_pt;
1629
1630 // Determine the edge on which the new node will live
1631 int father_bound = Father_bound[n_p](jnod, son_type);
1632
1633 // Storage for the set of Mesh boundaries on which the
1634 // appropriate father edge lives.
1635 std::set<unsigned> boundaries;
1636
1637 // Only get the boundaries if we are at the edge of
1638 // an element. Nodes in the centre of an element cannot be
1639 // on Mesh boundaries
1640 if (father_bound != Tree::OMEGA)
1641 {
1642 father_el_pt->get_boundaries(father_bound, boundaries);
1643 }
1644
1645#ifdef PARANOID
1646 // Case where a new node lives on more than one boundary
1647 // seems fishy enough to flag
1648 if (boundaries.size() > 2)
1649 {
1650 throw OomphLibError(
1651 "boundaries.size()>2 seems a bit strange..\n",
1652 OOMPH_CURRENT_FUNCTION,
1653 OOMPH_EXCEPTION_LOCATION);
1654 }
1655
1656 // Case when there are no boundaries, we are in big
1657 // trouble
1658 if (boundaries.size() == 0)
1659 {
1660 std::ostringstream error_stream;
1661 error_stream
1662 << "Periodic node is not on a boundary...\n"
1663 << "Coordinates: " << created_node_pt->x(0) << " "
1664 << created_node_pt->x(1) << "\n";
1665 throw OomphLibError(error_stream.str(),
1666 OOMPH_CURRENT_FUNCTION,
1667 OOMPH_EXCEPTION_LOCATION);
1668 }
1669#endif
1670
1671 // Create node and set the pointer to it from the element
1672 created_node_pt =
1673 construct_boundary_node(jnod, time_stepper_pt);
1674
1675 // Make the node periodic from the neighbour
1676 created_node_pt->make_periodic(neighbour_node_pt);
1677
1678 // Add to vector of new nodes
1679 new_node_pt.push_back(created_node_pt);
1680
1681 // Loop over # of history values
1682 for (unsigned t = 0; t < ntstorage; t++)
1683 {
1684 // Get position from father element -- this uses the
1685 // macro element representation if appropriate. If the
1686 // node turns out to be a hanging node later on, then
1687 // its position gets adjusted in line with its
1688 // hanging node interpolation.
1689 Vector<double> x_prev(n_dim, 0.0);
1690 father_el_pt->get_x(t, s, x_prev);
1691 // Set previous positions of the new node
1692 for (unsigned i = 0; i < n_dim; i++)
1693 {
1694 created_node_pt->x(t, i) = x_prev[i];
1695 }
1696 }
1697
1698 // Next, we Update the boundary lookup schemes
1699 // Loop over the boundaries stored in the set
1700 for (std::set<unsigned>::iterator it = boundaries.begin();
1701 it != boundaries.end();
1702 ++it)
1703 {
1704 // Add the node to the boundary
1705 mesh_pt->add_boundary_node(*it, created_node_pt);
1706
1707 // If we have set an intrinsic coordinate on this
1708 // mesh boundary then it must also be interpolated on
1709 // the new node
1710 // Now interpolate the intrinsic boundary coordinate
1711 if (mesh_pt->boundary_coordinate_exists(*it) == true)
1712 {
1713 Vector<double> zeta(2, 0.0);
1714 father_el_pt->interpolated_zeta_on_face(
1715 *it, father_bound, s, zeta);
1716 created_node_pt->set_coordinates_on_boundary(*it,
1717 zeta);
1718 }
1719 }
1720
1721 // Make sure that we add the node to the mesh
1722 mesh_pt->add_node_pt(created_node_pt);
1723 } // End of periodic case
1724 // Otherwise the node is not periodic, so just set the
1725 // pointer to the neighbours node
1726 else
1727 {
1728 node_pt(jnod) = created_node_pt;
1729 }
1730 // Node has been created
1731 node_done = true;
1732 } // Node does not exist in son of neighbouring element
1733 } // Node does not exist in neighbouring element
1734 } // Node does not exist in father element
1735
1736 // Node already exists in father: No need to do anything else!
1737 // otherwise deal with boundary information etc.
1738 if (!node_exists_in_father)
1739 {
1740 // Check boundary status
1741 //----------------------
1742
1743 // Firstly, we need to determine whether or not a node lies
1744 // on the boundary before building it, because
1745 // we actually assign a different type of node on boundaries.
1746
1747 // If the new node lives on a face that is
1748 // shared with a face of its father element,
1749 // it needs to inherit the bounday conditions
1750 // from the father face
1751 int father_bound = Father_bound[n_p](jnod, son_type);
1752
1753 // Storage for the set of Mesh boundaries on which the
1754 // appropriate father face lives.
1755 // [New nodes should always be mid-edge nodes in father
1756 // and therefore only live on one boundary but just to
1757 // play it safe...]
1758 std::set<unsigned> boundaries;
1759
1760 // Only get the boundaries if we are at the edge of
1761 // an element. Nodes in the centre of an element cannot be
1762 // on Mesh boundaries
1763 if (father_bound != Tree::OMEGA)
1764 {
1765 father_el_pt->get_boundaries(father_bound, boundaries);
1766 }
1767
1768#ifdef PARANOID
1769 // Case where a new node lives on more than two boundaries
1770 // seems fishy enough to flag
1771 if (boundaries.size() > 2)
1772 {
1773 throw OomphLibError(
1774 "boundaries.size()>2 seems a bit strange..\n",
1775 OOMPH_CURRENT_FUNCTION,
1776 OOMPH_EXCEPTION_LOCATION);
1777 }
1778#endif
1779
1780 // If the node lives on a mesh boundary,
1781 // then we need to create a boundary node
1782 if (boundaries.size() > 0)
1783 {
1784 // Do we need a new node?
1785 if (!node_done)
1786 {
1787 // Create node and set the internal pointer
1788 created_node_pt =
1789 construct_boundary_node(jnod, time_stepper_pt);
1790 // Add to vector of new nodes
1791 new_node_pt.push_back(created_node_pt);
1792 }
1793
1794 // Now we need to work out whether to pin the values at
1795 // the new node based on the boundary conditions applied at
1796 // its Mesh boundary
1797
1798 // Get the boundary conditions from the father.
1799 // Note: We can only deal with the values that are
1800 // continuously interpolated in the bulk element.
1801 unsigned n_cont = ncont_interpolated_values();
1802 Vector<int> bound_cons(n_cont);
1803 father_el_pt->get_bcs(father_bound, bound_cons);
1804
1805 // Loop over the continuously interpolated values and pin,
1806 // if necessary
1807 for (unsigned k = 0; k < n_cont; k++)
1808 {
1809 if (bound_cons[k])
1810 {
1811 created_node_pt->pin(k);
1812 }
1813 }
1814
1815 // Solid node? If so, deal with the positional boundary
1816 // conditions:
1817 SolidNode* solid_node_pt =
1818 dynamic_cast<SolidNode*>(created_node_pt);
1819 if (solid_node_pt != 0)
1820 {
1821 // Get the positional boundary conditions from the father:
1822 unsigned n_dim = created_node_pt->ndim();
1823 Vector<int> solid_bound_cons(n_dim);
1824 RefineableSolidQElement<3>* father_solid_el_pt =
1825 dynamic_cast<RefineableSolidQElement<3>*>(father_el_pt);
1826#ifdef PARANOID
1827 if (father_solid_el_pt == 0)
1828 {
1829 std::string error_message = "We have a SolidNode outside "
1830 "a refineable SolidElement\n";
1831 error_message +=
1832 "during mesh refinement -- this doesn't make sense";
1833
1834 throw OomphLibError(error_message,
1835 OOMPH_CURRENT_FUNCTION,
1836 OOMPH_EXCEPTION_LOCATION);
1837 }
1838#endif
1839 father_solid_el_pt->get_solid_bcs(father_bound,
1840 solid_bound_cons);
1841
1842 // Loop over the positions and pin, if necessary
1843 for (unsigned k = 0; k < n_dim; k++)
1844 {
1845 if (solid_bound_cons[k])
1846 {
1847 solid_node_pt->pin_position(k);
1848 }
1849 }
1850 } // End of if solid_node_pt
1851
1852 // Next update the boundary look-up schemes
1853
1854 // Loop over the boundaries in the set
1855 for (std::set<unsigned>::iterator it = boundaries.begin();
1856 it != boundaries.end();
1857 ++it)
1858 {
1859 // Add the node to the bounadry
1860 mesh_pt->add_boundary_node(*it, created_node_pt);
1861
1862 // If we have set an intrinsic coordinate on this
1863 // mesh boundary then it must also be interpolated on
1864 // the new node
1865 // Now interpolate the intrinsic boundary coordinate
1866 if (mesh_pt->boundary_coordinate_exists(*it) == true)
1867 {
1868 // Usually there will be two coordinates
1869 Vector<double> zeta(2);
1870 father_el_pt->interpolated_zeta_on_face(
1871 *it, father_bound, s, zeta);
1872
1873 created_node_pt->set_coordinates_on_boundary(*it, zeta);
1874 }
1875 }
1876 }
1877 // Otherwise the node is not on a Mesh boundary and
1878 // we create a normal "bulk" node
1879 else
1880 {
1881 // Do we need a new node?
1882 if (!node_done)
1883 {
1884 // Create node and set the pointer to it from the element
1885 created_node_pt = construct_node(jnod, time_stepper_pt);
1886 // Add to vector of new nodes
1887 new_node_pt.push_back(created_node_pt);
1888 }
1889 }
1890
1891 // In the first instance use macro element or FE representation
1892 // to create past and present nodal positions.
1893 // (THIS STEP SHOULD NOT BE SKIPPED FOR ALGEBRAIC
1894 // ELEMENTS AS NOT ALL OF THEM NECESSARILY IMPLEMENT
1895 // NONTRIVIAL NODE UPDATE FUNCTIONS. CALLING
1896 // THE NODE UPDATE FOR SUCH ELEMENTS/NODES WILL LEAVE
1897 // THEIR NODAL POSITIONS WHERE THEY WERE (THIS IS APPROPRIATE
1898 // ONCE THEY HAVE BEEN GIVEN POSITIONS) BUT WILL
1899 // NOT ASSIGN SENSIBLE INITIAL POSITONS!
1900
1901
1902 // Have we created a new node?
1903 if (!node_done)
1904 {
1905 // Loop over # of history values
1906 for (unsigned t = 0; t < ntstorage; t++)
1907 {
1908 // Get position from father element -- this uses the macro
1909 // element representation if appropriate. If the node
1910 // turns out to be a hanging node later on, then
1911 // its position gets adjusted in line with its
1912 // hanging node interpolation.
1913 Vector<double> x_prev(n_dim);
1914 father_el_pt->get_x(t, s, x_prev);
1915
1916 // Set previous positions of the new node
1917 for (unsigned i = 0; i < n_dim; i++)
1918 {
1919 created_node_pt->x(t, i) = x_prev[i];
1920 }
1921 }
1922
1923 // Now set the values
1924 // Loop over all history values
1925 for (unsigned t = 0; t < ntstorage; t++)
1926 {
1927 // Get values from father element
1928 // Note: get_interpolated_values() sets Vector size itself.
1929 Vector<double> prev_values;
1930 father_el_pt->get_interpolated_values(t, s, prev_values);
1931
1932 // Initialise the values at the new node
1933 unsigned n_value = created_node_pt->nvalue();
1934 for (unsigned k = 0; k < n_value; k++)
1935 {
1936 created_node_pt->set_value(t, k, prev_values[k]);
1937 }
1938 }
1939
1940 // Add new node to mesh
1941 mesh_pt->add_node_pt(created_node_pt);
1942
1943 } // End of whether the node has been created by us or not
1944
1945 } // End of if node already existed in father in which case
1946 // everything above gets bypassed.
1947
1948 // Check if the element is an algebraic element
1949 AlgebraicElementBase* alg_el_pt =
1950 dynamic_cast<AlgebraicElementBase*>(this);
1951
1952 // If the element is an algebraic element, setup
1953 // node position (past and present) from algebraic update
1954 // function.
1955 // NOTE: YES, THIS NEEDS TO BE CALLED REPEATEDLY IF THE
1956 // NODE IS MEMBER OF MULTIPLE ELEMENTS: THEY ALL ASSIGN
1957 // THE SAME NODAL POSITIONS BUT WE NEED TO ADD THE REMESH
1958 // INFO FOR *ALL* ROOT ELEMENTS!
1959 if (alg_el_pt != 0)
1960 {
1961 // Build algebraic node update info for new node
1962 // This sets up the node update data for all node update
1963 // functions that are shared by all nodes in the father
1964 // element
1965 alg_el_pt->setup_algebraic_node_update(
1966 node_pt(jnod), s, father_el_pt);
1967 }
1968
1969 // We have built the node and we are documenting
1970 if ((!node_done) && (new_nodes_file.is_open()))
1971 {
1972 new_nodes_file << node_pt(jnod)->x(0) << " "
1973 << node_pt(jnod)->x(1) << " "
1974 << node_pt(jnod)->x(2) << std::endl;
1975 }
1976
1977 } // End of Z loop over nodes in element
1978
1979 } // End of vertical loop over nodes in element
1980
1981 } // End of horizontal loop over nodes in element
1982
1983 // If the element is a MacroElementNodeUpdateElement, set
1984 // the update parameters for the current element's nodes --
1985 // all this needs is the vector of (pointers to the)
1986 // geometric objects that affect the MacroElement-based
1987 // node update -- this is the same as that in the father element
1988 MacroElementNodeUpdateElementBase* father_m_el_pt =
1989 dynamic_cast<MacroElementNodeUpdateElementBase*>(father_el_pt);
1990 if (father_m_el_pt != 0)
1991 {
1992 // Get vector of geometric objects from father (construct vector
1993 // via copy operation)
1994 Vector<GeomObject*> geom_object_pt(father_m_el_pt->geom_object_pt());
1995
1996 // Cast current element to MacroElementNodeUpdateElement:
1998 dynamic_cast<MacroElementNodeUpdateElementBase*>(this);
1999
2000#ifdef PARANOID
2001 if (m_el_pt == 0)
2002 {
2003 std::string error_message =
2004 "Failed to cast to MacroElementNodeUpdateElementBase*\n";
2005 error_message +=
2006 "Strange -- if the father is a MacroElementNodeUpdateElement\n";
2007 error_message += "the son should be too....\n";
2008
2009 throw OomphLibError(
2010 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2011 }
2012#endif
2013 // Build update info by passing vector of geometric objects:
2014 // This sets the current element to be the update element
2015 // for all of the element's nodes -- this is reversed
2016 // if the element is ever un-refined in the father element's
2017 // rebuild_from_sons() function which overwrites this
2018 // assignment to avoid nasty segmentation faults that occur
2019 // when a node tries to update itself via an element that no
2020 // longer exists...
2021 m_el_pt->set_node_update_info(geom_object_pt);
2022 }
2023
2024#ifdef OOMPH_HAS_MPI
2025 // Is the new element a halo element?
2026 if (tree_pt()->father_pt()->object_pt()->is_halo())
2027 {
2028 Non_halo_proc_ID =
2029 tree_pt()->father_pt()->object_pt()->non_halo_proc_ID();
2030 }
2031#endif
2032
2033 // Is it an ElementWithMovingNodes?
2034 ElementWithMovingNodes* aux_el_pt =
2035 dynamic_cast<ElementWithMovingNodes*>(this);
2036
2037 // Pass down the information re the method for the evaluation
2038 // of the shape derivatives
2039 if (aux_el_pt != 0)
2040 {
2041 ElementWithMovingNodes* aux_father_el_pt =
2042 dynamic_cast<ElementWithMovingNodes*>(father_el_pt);
2043
2044#ifdef PARANOID
2045 if (aux_father_el_pt == 0)
2046 {
2047 std::string error_message =
2048 "Failed to cast to ElementWithMovingNodes*\n";
2049 error_message +=
2050 "Strange -- if the son is a ElementWithMovingNodes\n";
2051 error_message += "the father should be too....\n";
2052
2053 throw OomphLibError(
2054 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2055 }
2056#endif
2057
2058 // If evaluating the residuals by finite differences in the father
2059 // continue to do so in the child
2060 if (aux_father_el_pt
2061 ->are_dresidual_dnodal_coordinates_always_evaluated_by_fd())
2062 {
2063 aux_el_pt
2065 }
2066
2067
2068 aux_el_pt->method_for_shape_derivs() =
2069 aux_father_el_pt->method_for_shape_derivs();
2070
2071 // If bypassing the evaluation of fill_in_jacobian_from_geometric_data
2072 // continue to do so
2073 if (aux_father_el_pt
2074 ->is_fill_in_jacobian_from_geometric_data_bypassed())
2075 {
2077 }
2078 }
2079
2080 // Now do further build (if any)
2081 further_build();
2082
2083 } // Sanity check: Father element has been generated
2084
2085 } // End for element has not been built yet
2086 else
2087 {
2088 was_already_built = true;
2089 }
2090 }
2091
2092 //====================================================================
2093 /// Set up all hanging nodes.
2094 ///
2095 /// (Wrapper to avoid specification of the unwanted output file).
2096 //====================================================================
2098 Vector<std::ofstream*>& output_stream)
2099 {
2100#ifdef PARANOID
2101 if (output_stream.size() != 6)
2102 {
2103 throw OomphLibError("There must be six output streams",
2104 OOMPH_CURRENT_FUNCTION,
2105 OOMPH_EXCEPTION_LOCATION);
2106 }
2107#endif
2108
2109 using namespace OcTreeNames;
2110
2111 // Setup hanging nodes on each edge of the element
2112 oc_hang_helper(-1, U, *(output_stream[0]));
2113 oc_hang_helper(-1, D, *(output_stream[1]));
2114 oc_hang_helper(-1, L, *(output_stream[2]));
2115 oc_hang_helper(-1, R, *(output_stream[3]));
2116 oc_hang_helper(-1, B, *(output_stream[4]));
2117 oc_hang_helper(-1, F, *(output_stream[5]));
2118 }
2119
2120 //================================================================
2121 /// Internal function that sets up the hanging node scheme for a
2122 /// particular value
2123 //=================================================================
2125 {
2126 using namespace OcTreeNames;
2127
2128 std::ofstream dummy_hangfile;
2129 // Setup hanging nodes on each edge of the element
2130 oc_hang_helper(value_id, U, dummy_hangfile);
2131 oc_hang_helper(value_id, D, dummy_hangfile);
2132 oc_hang_helper(value_id, R, dummy_hangfile);
2133 oc_hang_helper(value_id, L, dummy_hangfile);
2134 oc_hang_helper(value_id, B, dummy_hangfile);
2135 oc_hang_helper(value_id, F, dummy_hangfile);
2136 }
2137
2138 //=================================================================
2139 /// Internal function to set up the hanging nodes on a particular
2140 /// face of the element
2141 //=================================================================
2143 const int& my_face,
2144 std::ofstream& output_hangfile)
2145 {
2146 using namespace OcTreeNames;
2147
2148 // Number of dimensions
2149 unsigned n_dim = 3;
2150
2151 Vector<unsigned> translate_s(n_dim);
2152 Vector<double> s_lo_neigh(n_dim);
2153 Vector<double> s_hi_neigh(n_dim);
2154 int neigh_face;
2155 int diff_level;
2156 bool in_neighbouring_tree;
2157
2158 // Find pointer to neighbour in this direction
2159 OcTree* neigh_pt;
2160 neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2161 translate_s,
2162 s_lo_neigh,
2163 s_hi_neigh,
2164 neigh_face,
2165 diff_level,
2166 in_neighbouring_tree);
2167
2168 // Neighbour exists
2169 if (neigh_pt != 0)
2170 {
2171 // Different sized element?
2172 if (diff_level != 0)
2173 {
2174 // Test for the periodic node case
2175 // Are we crossing a periodic boundary
2176 bool is_periodic = false;
2177 if (in_neighbouring_tree)
2178 {
2179 is_periodic = tree_pt()->root_pt()->is_neighbour_periodic(my_face);
2180 }
2181
2182 // If it is periodic we actually need to get the node in
2183 // the neighbour of the neighbour (which will be a parent of
2184 // the present element) so that the "fixed" coordinate is
2185 // correctly calculated.
2186 // The idea is to replace the neigh_pt and associated data
2187 // with those of the neighbour of the neighbour
2188 if (is_periodic)
2189 {
2190 // Required data for the neighbour finding routine
2191 Vector<unsigned> translate_s_in_neigh(n_dim, 0.0);
2192 Vector<double> s_lo_neigh_of_neigh(n_dim, 0.0);
2193 Vector<double> s_hi_neigh_of_neigh(n_dim, 0.0);
2194 int neigh_face_of_neigh;
2195 int diff_level_of_neigh;
2196 bool in_neighbouring_tree_of_neigh;
2197
2198 // Find pointer to neighbour of the neighbour on the edge
2199 // that we are currently considering
2200 OcTree* neigh_of_neigh_pt;
2201 neigh_of_neigh_pt =
2202 neigh_pt->gteq_face_neighbour(neigh_face,
2203 translate_s_in_neigh,
2204 s_lo_neigh_of_neigh,
2205 s_hi_neigh_of_neigh,
2206 neigh_face_of_neigh,
2207 diff_level_of_neigh,
2208 in_neighbouring_tree_of_neigh);
2209
2210 // Set the value of the NEW neighbour and edge
2211 neigh_pt = neigh_of_neigh_pt;
2212 neigh_face = neigh_face_of_neigh;
2213
2214 // Set the values of the translation schemes
2215 // Need to find the values of s_lo and s_hi
2216 // in the neighbour of the neighbour
2217
2218 // Get the minimum and maximum values of the coordinate
2219 // in the neighbour (don't like this, but I think it's
2220 // necessary) Note that these values are hardcoded
2221 // in the quadtrees at some point!!
2222 double s_min = neigh_pt->object_pt()->s_min();
2223 double s_max = neigh_pt->object_pt()->s_max();
2224 Vector<double> s_lo_frac(n_dim), s_hi_frac(n_dim);
2225 // Work out the fractional position of the low and high points
2226 // of the original element
2227 for (unsigned i = 0; i < n_dim; i++)
2228 {
2229 s_lo_frac[i] = (s_lo_neigh[i] - s_min) / (s_max - s_min);
2230 s_hi_frac[i] = (s_hi_neigh[i] - s_min) / (s_max - s_min);
2231 }
2232
2233 // We should now be able to construct the low and high points in
2234 // the neighbour of the neighbour
2235 for (unsigned i = 0; i < n_dim; i++)
2236 {
2237 s_lo_neigh[i] = s_lo_neigh_of_neigh[i] +
2238 s_lo_frac[translate_s_in_neigh[i]] *
2239 (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2240 s_hi_neigh[i] = s_lo_neigh_of_neigh[i] +
2241 s_hi_frac[translate_s_in_neigh[i]] *
2242 (s_hi_neigh_of_neigh[i] - s_lo_neigh_of_neigh[i]);
2243 }
2244
2245 // Finally we must sort out the translation scheme
2246 Vector<unsigned> temp_translate(n_dim, 0.0);
2247 for (unsigned i = 0; i < n_dim; i++)
2248 {
2249 temp_translate[i] = translate_s_in_neigh[translate_s[i]];
2250 }
2251 for (unsigned i = 0; i < n_dim; i++)
2252 {
2253 translate_s[i] = temp_translate[i];
2254 }
2255 } // End of special treatment for periodic hanging nodes
2256
2257 // Number of nodes in one dimension
2258 unsigned n_p = ninterpolating_node_1d(value_id);
2259 // Storage for the local nodes along the face of the element
2260 Node* local_node_pt = 0;
2261
2262 // Loop over nodes along the face
2263 for (unsigned i0 = 0; i0 < n_p; i0++)
2264 {
2265 for (unsigned i1 = 0; i1 < n_p; i1++)
2266 {
2267 // Storage for the fractional position of the node in the element
2268 Vector<double> s_fraction(n_dim);
2269
2270 // Local node number
2271 switch (my_face)
2272 {
2273 case U:
2274 s_fraction[0] =
2275 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2276 s_fraction[1] = 1.0;
2277 s_fraction[2] =
2278 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2279 local_node_pt = interpolating_node_pt(
2280 i0 + (n_p - 1) * n_p + n_p * n_p * i1, value_id);
2281 break;
2282
2283 case D:
2284 s_fraction[0] =
2285 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2286 s_fraction[1] = 0.0;
2287 s_fraction[2] =
2288 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2289 local_node_pt =
2290 interpolating_node_pt(i0 + n_p * n_p * i1, value_id);
2291 break;
2292
2293 case R:
2294 s_fraction[0] = 1.0;
2295 s_fraction[1] =
2296 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2297 s_fraction[2] =
2298 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2299 local_node_pt = interpolating_node_pt(
2300 n_p - 1 + i0 * n_p + i1 * n_p * n_p, value_id);
2301 break;
2302
2303 case L:
2304 s_fraction[0] = 0.0;
2305 s_fraction[1] =
2306 local_one_d_fraction_of_interpolating_node(i0, 1, value_id);
2307 s_fraction[2] =
2308 local_one_d_fraction_of_interpolating_node(i1, 2, value_id);
2309 local_node_pt =
2310 interpolating_node_pt(n_p * i0 + i1 * n_p * n_p, value_id);
2311 break;
2312
2313 case B:
2314 s_fraction[0] =
2315 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2316 s_fraction[1] =
2317 local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2318 s_fraction[2] = 0.0;
2319 local_node_pt = interpolating_node_pt(i0 + i1 * n_p, value_id);
2320 break;
2321
2322 case F:
2323 s_fraction[0] =
2324 local_one_d_fraction_of_interpolating_node(i0, 0, value_id);
2325 s_fraction[1] =
2326 local_one_d_fraction_of_interpolating_node(i1, 1, value_id);
2327 s_fraction[2] = 1.0;
2328 local_node_pt = interpolating_node_pt(
2329 i0 + i1 * n_p + (n_p - 1) * n_p * n_p, value_id);
2330 break;
2331
2332 default:
2333 throw OomphLibError("my_face not U, D, L, R, B, F\n",
2334 OOMPH_CURRENT_FUNCTION,
2335 OOMPH_EXCEPTION_LOCATION);
2336 }
2337
2338 // Set up the coordinate in the neighbour
2339 Vector<double> s_in_neighb(n_dim);
2340 for (unsigned i = 0; i < n_dim; i++)
2341 {
2342 s_in_neighb[i] =
2343 s_lo_neigh[i] +
2344 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2345 }
2346
2347 // Find the Node in the neighbouring element
2348 Node* neighbouring_node_pt =
2350 s_in_neighb, value_id);
2351
2352 // If the neighbour does not have a node at this point
2353 if (0 == neighbouring_node_pt)
2354 {
2355 // Do we need to make a hanging node, we assume that we don't
2356 // initially
2357 bool make_hanging_node = false;
2358
2359 // If the node is not hanging geometrically, then we must make it
2360 // hang
2361 if (!local_node_pt->is_hanging())
2362 {
2363 make_hanging_node = true;
2364 }
2365 // Otherwise is could be hanging geometrically, but still require
2366 // a different hanging scheme for this data value
2367 else
2368 {
2369 if ((value_id != -1) && (local_node_pt->hanging_pt(value_id) ==
2370 local_node_pt->hanging_pt()))
2371 {
2372 make_hanging_node = true;
2373 }
2374 }
2375
2376 if (make_hanging_node == true)
2377 {
2378 // Cache refineable element used here
2379 RefineableElement* const obj_pt = neigh_pt->object_pt();
2380
2381 // Get shape functions in neighbour element
2382 Shape psi(obj_pt->ninterpolating_node(value_id));
2383 obj_pt->interpolating_basis(s_in_neighb, psi, value_id);
2384
2385 // Allocate the storage for the Hang pointer
2386 // We know that it will hold n_p*n_p nodes
2387 HangInfo* hang_pt = new HangInfo(n_p * n_p);
2388
2389 // Loop over nodes on edge in neighbour and mark them as nodes
2390 // that this node depends on
2391 unsigned n_neighbour;
2392
2393 // Number of nodes along edge in neighbour element
2394 for (unsigned ii0 = 0; ii0 < n_p; ii0++)
2395 {
2396 for (unsigned ii1 = 0; ii1 < n_p; ii1++)
2397 {
2398 switch (neigh_face)
2399 {
2400 case U:
2401 n_neighbour = ii0 + n_p * (n_p - 1) + ii1 * n_p * n_p;
2402 break;
2403
2404 case D:
2405 n_neighbour = ii0 + ii1 * n_p * n_p;
2406 break;
2407
2408 case L:
2409 n_neighbour = ii0 * n_p + ii1 * n_p * n_p;
2410 break;
2411
2412 case R:
2413 n_neighbour = (n_p - 1) + ii0 * n_p + ii1 * n_p * n_p;
2414 break;
2415
2416 case B:
2417 n_neighbour = ii0 + ii1 * n_p;
2418 break;
2419
2420 case F:
2421 n_neighbour = ii0 + ii1 * n_p + n_p * n_p * (n_p - 1);
2422 break;
2423 default:
2424 throw OomphLibError("neigh_face not U, L, R, B, F\n",
2425 OOMPH_CURRENT_FUNCTION,
2426 OOMPH_EXCEPTION_LOCATION);
2427 }
2428
2429 // Push back neighbouring node and weight into
2430 // Vector of (pointers to)
2431 // master nodes and weights
2432 // The weight is merely the value of the shape function
2433 // corresponding to the node in the neighbour
2434 hang_pt->set_master_node_pt(
2435 ii0 * n_p + ii1,
2436 obj_pt->interpolating_node_pt(n_neighbour, value_id),
2437 psi[n_neighbour]);
2438 }
2439 }
2440 // Now set the hanging data for the position
2441 // This also constrains the data values associated with the
2442 // value id
2443 local_node_pt->set_hanging_pt(hang_pt, value_id);
2444 }
2445
2446 if (output_hangfile.is_open())
2447 {
2448 // output_hangfile
2449 output_hangfile << local_node_pt->x(0) << " "
2450 << local_node_pt->x(1) << " "
2451 << local_node_pt->x(2) << std::endl;
2452 }
2453 }
2454 else
2455 {
2456#ifdef PARANOID
2457 if (local_node_pt != neighbouring_node_pt)
2458 {
2459 std::ofstream reportage("dodgy.dat", std::ios_base::app);
2460 reportage << local_node_pt->x(0) << " " << local_node_pt->x(1)
2461 << " " << local_node_pt->x(2) << std::endl;
2462 reportage.close();
2463
2464 std::ostringstream warning_stream;
2465 warning_stream
2466 << "SANITY CHECK in oc_hang_helper \n"
2467 << "Current node " << local_node_pt << " at "
2468 << "(" << local_node_pt->x(0) << ", " << local_node_pt->x(1)
2469 << ", " << local_node_pt->x(2) << ")" << std::endl
2470 << " is not hanging and has " << std::endl
2471 << "Neighbour's node " << neighbouring_node_pt << " at "
2472 << "(" << neighbouring_node_pt->x(0) << ", "
2473 << neighbouring_node_pt->x(1) << ", "
2474 << neighbouring_node_pt->x(2) << ")" << std::endl
2475 << "even though the two should be "
2476 << "identical" << std::endl;
2477 OomphLibWarning(warning_stream.str(),
2478 OOMPH_CURRENT_FUNCTION,
2479 OOMPH_EXCEPTION_LOCATION);
2480 }
2481#endif
2482 }
2483
2484 // If we are doing the position then
2485 if (value_id == -1)
2486 {
2487 // Get the nodal position from neighbour element
2488 Vector<double> x_in_neighb(n_dim);
2489 neigh_pt->object_pt()->interpolated_x(s_in_neighb, x_in_neighb);
2490
2491 // Fine adjust the coordinates (macro map will pick up boundary
2492 // accurately but will lead to different element edges)
2493 local_node_pt->x(0) = x_in_neighb[0];
2494 local_node_pt->x(1) = x_in_neighb[1];
2495 local_node_pt->x(2) = x_in_neighb[2];
2496 }
2497 }
2498 }
2499 }
2500 }
2501 }
2502
2503
2504 //=================================================================
2505 /// Check inter-element continuity of
2506 /// - nodal positions
2507 /// - (nodally) interpolated function values
2508 //====================================================================
2510 {
2511 using namespace OcTreeNames;
2512
2513 // std::ofstream error_file("errors.dat");
2514
2515 // Number of nodes along edge
2516 unsigned n_p = nnode_1d();
2517
2518 // Number of timesteps (incl. present) for which continuity is
2519 // to be checked.
2520 unsigned n_time = 1;
2521
2522 // Initialise errors
2523 max_error = 0.0;
2524 Vector<double> max_error_x(3, 0.0);
2525 double max_error_val = 0.0;
2526
2527 // Set the faces
2528 Vector<int> faces(6);
2529 faces[0] = D;
2530 faces[1] = U;
2531 faces[2] = L;
2532 faces[3] = R;
2533 faces[4] = B;
2534 faces[5] = F;
2535
2536 // Loop over the edges
2537 for (unsigned face_counter = 0; face_counter < 6; face_counter++)
2538 {
2539 Vector<double> s(3), s_lo_neigh(3), s_hi_neigh(3);
2540 Vector<double> s_fraction(3);
2541 Vector<unsigned> translate_s(3);
2542 int neigh_face, diff_level;
2543 int my_face;
2544 bool in_neighbouring_tree;
2545
2546 my_face = faces[face_counter];
2547
2548 // Find pointer to neighbour in this direction
2549 OcTree* neigh_pt;
2550 neigh_pt = octree_pt()->gteq_face_neighbour(my_face,
2551 translate_s,
2552 s_lo_neigh,
2553 s_hi_neigh,
2554 neigh_face,
2555 diff_level,
2556 in_neighbouring_tree);
2557
2558 // Neighbour exists and has existing nodes
2559 if ((neigh_pt != 0) && (neigh_pt->object_pt()->nodes_built()))
2560 {
2561 // if (diff_level!=0)
2562 {
2563 // Need to exclude periodic nodes from this check
2564 // There are only periodic nodes if we are in a neighbouring tree
2565 bool is_periodic = false;
2566 if (in_neighbouring_tree)
2567 {
2568 // Is it periodic
2569 is_periodic =
2570 tree_pt()->root_pt()->is_neighbour_periodic(faces[face_counter]);
2571 }
2572
2573 // Loop over nodes along the edge
2574 for (unsigned i0 = 0; i0 < n_p; i0++)
2575 {
2576 for (unsigned i1 = 0; i1 < n_p; i1++)
2577 {
2578 // Local node number
2579 unsigned n = 0;
2580 switch (face_counter)
2581 {
2582 case 0:
2583 // Fractions
2584 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2585 s_fraction[1] = 0.0;
2586 s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2587 // Set local node number
2588 n = i0 + i1 * n_p * n_p;
2589 break;
2590
2591 case 1:
2592 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2593 s_fraction[1] = 1.0;
2594 s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2595 // Set local node number
2596 n = i0 + n_p * (n_p - 1) + i1 * n_p * n_p;
2597 break;
2598
2599 case 2:
2600 s_fraction[0] = 0.0;
2601 s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2602 s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2603 // Set local node number
2604 n = n_p * i0 + i1 * n_p * n_p;
2605 break;
2606
2607 case 3:
2608 s_fraction[0] = 1.0;
2609 s_fraction[1] = local_one_d_fraction_of_node(i0, 1);
2610 s_fraction[2] = local_one_d_fraction_of_node(i1, 2);
2611 // Set local node number
2612 n = n_p - 1 + n_p * i0 + n_p * n_p * i1;
2613 break;
2614
2615 case 4:
2616 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2617 s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2618 s_fraction[2] = 0.0;
2619 // Set local node number
2620 n = i0 + n_p * i1;
2621 break;
2622
2623 case 5:
2624 s_fraction[0] = local_one_d_fraction_of_node(i0, 0);
2625 s_fraction[1] = local_one_d_fraction_of_node(i1, 1);
2626 s_fraction[2] = 1.0;
2627 // Set local node number
2628 n = i0 + n_p * i1 + n_p * n_p * (n_p - 1);
2629 break;
2630 }
2631
2632
2633 // We have to check if the hi and lo directions along the
2634 // face are inverted or not
2635 Vector<double> s_in_neighb(3);
2636 for (unsigned i = 0; i < 3; i++)
2637 {
2638 s[i] = -1.0 + 2.0 * s_fraction[i];
2639 s_in_neighb[i] =
2640 s_lo_neigh[i] +
2641 s_fraction[translate_s[i]] * (s_hi_neigh[i] - s_lo_neigh[i]);
2642 }
2643
2644
2645 // Loop over timesteps
2646 for (unsigned t = 0; t < n_time; t++)
2647 {
2648 // Get the nodal position from neighbour element
2649 Vector<double> x_in_neighb(3);
2650 neigh_pt->object_pt()->interpolated_x(
2651 t, s_in_neighb, x_in_neighb);
2652
2653 // Check error only if the node is NOT periodic
2654 if (is_periodic == false)
2655 {
2656 // Check error
2657 for (unsigned i = 0; i < 3; i++)
2658 {
2659 // Find the spatial error
2660 double err =
2661 std::fabs(node_pt(n)->x(t, i) - x_in_neighb[i]);
2662
2663 // If it's bigger than our tolerance, say so
2664 if (err > 1e-9)
2665 {
2666 oomph_info << "errx[" << i << "], t x, x_neigh: " << err
2667 << " " << t << " " << node_pt(n)->x(t, i)
2668 << " " << x_in_neighb[i] << std::endl;
2669 oomph_info << "at " << node_pt(n)->x(0) << " "
2670 << node_pt(n)->x(1) << " " << node_pt(n)->x(2)
2671 << " " << std::endl;
2672 }
2673
2674 // If it's bigger than the previous max error, it is the
2675 // new max error!
2676 if (err > max_error_x[i])
2677 {
2678 max_error_x[i] = err;
2679 }
2680 }
2681 }
2682
2683 // Get the values from neighbour element. Note: # of values
2684 // gets set by routine (because in general we don't know
2685 // how many interpolated values a certain element has
2686 Vector<double> values_in_neighb;
2688 t, s_in_neighb, values_in_neighb);
2689
2690 // Get the values in current element.
2691 Vector<double> values;
2692 get_interpolated_values(t, s, values);
2693
2694 // Now figure out how many continuously interpolated
2695 // values there are
2696 unsigned num_val =
2697 neigh_pt->object_pt()->ncont_interpolated_values();
2698
2699 // Check error
2700 for (unsigned ival = 0; ival < num_val; ival++)
2701 {
2702 double err = std::fabs(values[ival] - values_in_neighb[ival]);
2703
2704 if (err > 1.0e-10)
2705 {
2706 oomph_info << node_pt(n)->x(0) << " " << node_pt(n)->x(1)
2707 << " " << node_pt(n)->x(2) << " \n# "
2708 << "erru (S)" << err << " " << ival << " " << n
2709 << " " << values[ival] << " "
2710 << values_in_neighb[ival] << std::endl;
2711 // error_file<<"ZONE"<<std::endl
2712 // << node_pt(n)->x(0) << " "
2713 // << node_pt(n)->x(1) << " "
2714 // << node_pt(n)->x(2) << std::endl;
2715 }
2716
2717 if (err > max_error_val)
2718 {
2719 max_error_val = err;
2720 }
2721 }
2722 }
2723 }
2724 }
2725 }
2726 }
2727 }
2728
2729 max_error = max_error_x[0];
2730 if (max_error_x[1] > max_error) max_error = max_error_x[1];
2731 if (max_error_x[2] > max_error) max_error = max_error_x[2];
2732 if (max_error_val > max_error) max_error = max_error_val;
2733
2734 if (max_error > 1e-9)
2735 {
2736 oomph_info << "\n#------------------------------------ \n#Max error ";
2737 oomph_info << max_error_x[0] << " " << max_error_x[1] << " "
2738 << max_error_x[2] << " " << max_error_val << std::endl;
2739 oomph_info << "#------------------------------------ \n " << std::endl;
2740 }
2741
2742 // error_file.close();
2743 }
2744
2745
2746 //==================================================================
2747 /// Determine vector of solid (positional) boundary conditions along
2748 /// the element's boundary (or vertex) bound (S/W/N/E/SW/SE/NW/NE).
2749 ///
2750 /// This function assumes that the same boundary condition is applied
2751 /// along the entire length of an element's edge (of course, the
2752 /// vertices combine the boundary conditions of their two adjacent edges
2753 /// in the most restrictive combination. Hence, if we're at a vertex,
2754 /// we apply the most restrictive boundary condition of the
2755 /// two adjacent edges. If we're on an edge (in its proper interior),
2756 /// we apply the least restrictive boundary condition of all nodes
2757 /// along the edge.
2758 ///
2759 /// Usual convention:
2760 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2761 /// on this boundary is free.
2762 /// - solid_bound_cons[i]=1 if it's pinned.
2763 //==================================================================
2765 int bound, Vector<int>& solid_bound_cons) const
2766 {
2767 using namespace OcTreeNames;
2768
2769 // Spatial dimension of all nodes
2770 unsigned n_dim = this->nodal_dimension();
2771 // Set up temporary vectors to hold edge boundary conditions
2772 Vector<int> bound_cons1(n_dim), bound_cons2(n_dim);
2773 Vector<int> bound_cons3(n_dim);
2774
2775 Vector<int> vect1(3), vect2(3), vect3(3);
2776 Vector<int> vect_elem;
2777 Vector<int> notzero;
2778 int n = 0;
2779
2780 vect_elem = OcTree::Direction_to_vector[bound];
2781
2782 // Just to see if bound is a face, an edge, or a vertex, n stores
2783 // the number of non-zero values in the vector reprensenting the bound
2784 // and the vector notzero stores the position of these values
2785 for (int i = 0; i < 3; i++)
2786 {
2787 if (vect_elem[i] != 0)
2788 {
2789 n++;
2790 notzero.push_back(i);
2791 }
2792 }
2793
2794 switch (n)
2795 {
2796 // If there is only one non-zero value, bound is a face
2797 case 1:
2798 get_face_solid_bcs(bound, solid_bound_cons);
2799 break;
2800
2801 // If there are two non-zero values, bound is an edge
2802 case 2:
2803
2804 for (unsigned i = 0; i < 3; i++)
2805 {
2806 vect1[i] = 0;
2807 vect2[i] = 0;
2808 }
2809 // vect1 and vect2 are the vector of the two faces adjacent to bound
2810 vect1[notzero[0]] = vect_elem[notzero[0]];
2811 vect2[notzero[1]] = vect_elem[notzero[1]];
2812
2813 get_face_solid_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
2814 get_face_solid_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
2815 // get the most restrictive bc
2816 for (unsigned k = 0; k < n_dim; k++)
2817 {
2818 solid_bound_cons[k] = (bound_cons1[k] || bound_cons2[k]);
2819 }
2820 break;
2821
2822 // If there are three non-zero value, bound is a vertex
2823 case 3:
2824
2825 for (unsigned i = 0; i < 3; i++)
2826 {
2827 vect1[i] = 0;
2828 vect2[i] = 0;
2829 vect3[i] = 0;
2830 }
2831 // vectors to the three adjacent faces of the vertex
2832 vect1[0] = vect_elem[0];
2833 vect2[1] = vect_elem[1];
2834 vect3[2] = vect_elem[2];
2835
2836 get_face_solid_bcs(OcTree::Vector_to_direction[vect1], bound_cons1);
2837 get_face_solid_bcs(OcTree::Vector_to_direction[vect2], bound_cons2);
2838 get_face_solid_bcs(OcTree::Vector_to_direction[vect3], bound_cons3);
2839
2840
2841 // set the bcs to the most restrictive ones
2842 for (unsigned k = 0; k < n_dim; k++)
2843 {
2844 solid_bound_cons[k] =
2845 (bound_cons1[k] || bound_cons2[k] || bound_cons3[k]);
2846 }
2847 break;
2848
2849 default:
2850 throw OomphLibError("Make sure you are not giving OMEGA as bound",
2851 OOMPH_CURRENT_FUNCTION,
2852 OOMPH_EXCEPTION_LOCATION);
2853 }
2854 }
2855
2856 //==================================================================
2857 /// Determine Vector of solid (positional) boundary conditions along
2858 /// the element's edge (S/N/W/E) -- BC is the least restrictive combination
2859 /// of all the nodes on this edge
2860 ///
2861 /// Usual convention:
2862 /// - solid_bound_cons[i]=0 if displacement in coordinate direction i
2863 /// on this boundary is free
2864 /// - solid_bound_cons[i]=1 if it's pinned
2865 //==================================================================
2867 const int& face, Vector<int>& solid_bound_cons) const
2868 {
2869 using namespace OcTreeNames;
2870
2871 // Number of nodes along 1D edge
2872 unsigned n_p = nnode_1d();
2873 // the four corner nodes on the boundary
2874 unsigned node1, node2, node3, node4;
2875
2876 // Set the four corner nodes for the face
2877 switch (face)
2878 {
2879 case U:
2880 node1 = n_p * n_p * n_p - 1;
2881 node2 = n_p * n_p - 1;
2882 node3 = n_p * (n_p - 1);
2883 node4 = n_p * (n_p * n_p - 1);
2884 break;
2885
2886 case D:
2887 node1 = 0;
2888 node2 = n_p - 1;
2889 node3 = (n_p * n_p + 1) * (n_p - 1);
2890 node4 = n_p * n_p * (n_p - 1);
2891 break;
2892
2893 case R:
2894 node1 = n_p - 1;
2895 node2 = (n_p * n_p + 1) * (n_p - 1);
2896 node3 = n_p * n_p * n_p - 1;
2897 node4 = n_p * n_p - 1;
2898 break;
2899
2900 case L:
2901 node1 = 0;
2902 node2 = n_p * (n_p - 1);
2903 node3 = n_p * (n_p * n_p - 1);
2904 node4 = n_p * n_p * (n_p - 1);
2905 break;
2906
2907 case B:
2908 node1 = 0;
2909 node2 = n_p - 1;
2910 node3 = n_p * n_p - 1;
2911 node4 = n_p * (n_p - 1);
2912 break;
2913
2914 case F:
2915 node1 = n_p * n_p * n_p - 1;
2916 node2 = n_p * (n_p * n_p - 1);
2917 node3 = n_p * n_p * (n_p - 1);
2918 node4 = (n_p - 1) * (n_p * n_p + 1);
2919 break;
2920
2921 default:
2922 std::ostringstream error_stream;
2923 error_stream << "Wrong edge " << face << " passed\n";
2924
2925 throw OomphLibError(
2926 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2927 }
2928
2929 // Cast to solid nodes
2930 SolidNode* solid_node1_pt = dynamic_cast<SolidNode*>(node_pt(node1));
2931 SolidNode* solid_node2_pt = dynamic_cast<SolidNode*>(node_pt(node2));
2932 SolidNode* solid_node3_pt = dynamic_cast<SolidNode*>(node_pt(node3));
2933 SolidNode* solid_node4_pt = dynamic_cast<SolidNode*>(node_pt(node4));
2934
2935 //#ifdef PARANOID
2936 if (solid_node1_pt == 0)
2937 {
2938 throw OomphLibError(
2939 "Corner node 1 cannot be cast to SolidNode --> something is wrong",
2940 OOMPH_CURRENT_FUNCTION,
2941 OOMPH_EXCEPTION_LOCATION);
2942 }
2943
2944 if (solid_node2_pt == 0)
2945 {
2946 throw OomphLibError(
2947 "Corner node 2 cannot be cast to SolidNode --> something is wrong",
2948 OOMPH_CURRENT_FUNCTION,
2949 OOMPH_EXCEPTION_LOCATION);
2950 }
2951
2952 if (solid_node3_pt == 0)
2953 {
2954 throw OomphLibError(
2955 "Corner node 3 cannot be cast to SolidNode --> something is wrong",
2956 OOMPH_CURRENT_FUNCTION,
2957 OOMPH_EXCEPTION_LOCATION);
2958 }
2959
2960 if (solid_node4_pt == 0)
2961 {
2962 throw OomphLibError(
2963 "Corner node 4 cannot be cast to SolidNode --> something is wrong",
2964 OOMPH_CURRENT_FUNCTION,
2965 OOMPH_EXCEPTION_LOCATION);
2966 }
2967
2968 //#endif
2969
2970 // Number of coordinates
2971 unsigned n_dim = this->nodal_dimension();
2972
2973 // Loop over all directions, the least restrictive value is
2974 // the multiplication of the boundary conditions at the 4 nodes
2975 // Assuming that free is always zero and pinned is one
2976 for (unsigned k = 0; k < n_dim; k++)
2977 {
2978 solid_bound_cons[k] = solid_node1_pt->position_is_pinned(k) *
2979 solid_node2_pt->position_is_pinned(k) *
2980 solid_node3_pt->position_is_pinned(k) *
2981 solid_node4_pt->position_is_pinned(k);
2982 }
2983 }
2984
2985
2986 //========================================================================
2987 /// Static matrix for coincidence between son nodal points and
2988 /// father boundaries
2989 ///
2990 //========================================================================
2991 std::map<unsigned, DenseMatrix<int>> RefineableQElement<3>::Father_bound;
2992
2993
2994} // 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
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
MacroElement * macro_elem_pt()
Access function to pointer to macro element.
Definition: elements.h:1878
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
OcTree class: Recursively defined, generalised octree.
Definition: octree.h:114
static Vector< int > faces_of_common_edge(const int &edge)
Function that, given an edge, returns the two faces on which it.
Definition: octree.cc:268
static Vector< std::string > Direct_string
Translate (enumerated) directions into strings.
Definition: octree.h:329
OcTree * gteq_face_neighbour(const int &direction, Vector< unsigned > &translate_s, Vector< double > &s_sw, Vector< double > &s_ne, int &face, int &diff_level, bool &in_neighbouring_tree) const
Find (pointer to) ‘greater-or-equal-sized face neighbour’ in given direction (L/R/U/D/F/B)....
Definition: octree.cc:3373
OcTree * gteq_true_edge_neighbour(const int &direction, const unsigned &i_root_edge_neighbour, unsigned &nroot_edge_neighbour, Vector< unsigned > &translate_s, Vector< double > &s_lo, Vector< double > &s_hi, int &edge, int &diff_level) const
Find (pointer to) ‘greater-or-equal-sized true edge neighbour’ in the given direction (LB,...
Definition: octree.cc:3618
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:353
static std::map< Vector< int >, int > Vector_to_direction
Each vector representing a direction can be translated into a direction, either a son type (vertex),...
Definition: octree.h:348
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
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<3,NNODE_1D>.
void interpolated_zeta_on_face(const unsigned &boundary, const int &face, const Vector< double > &s, Vector< double > &zeta)
Return the value of the intrinsic boundary coordinate interpolated along the face.
void get_bcs(int bound, Vector< int > &bound_cons) const
Determine Vector of boundary conditions along the element's boundary (or vertex) bound (S/W/N/E/SW/SE...
void get_boundaries(const int &edge, std::set< unsigned > &boundaries) const
Given an element edge/vertex, return a Vector which contains all the (mesh-)boundary numbers that thi...
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 brick 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
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...