tube_mesh.template.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-2023 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 #ifndef OOMPH_TUBE_MESH_TEMPLATE_CC
27 #define OOMPH_TUBE_MESH_TEMPLATE_CC
28 
29 #include "tube_mesh.template.h"
30 
31 
32 namespace oomph
33 {
34  //====================================================================
35  /// Constructor for deformable quarter tube mesh class.
36  /// The domain is specified by the GeomObject that
37  /// identifies the entire volume.
38  //====================================================================
39  template<class ELEMENT>
40  TubeMesh<ELEMENT>::TubeMesh(GeomObject* volume_pt,
41  const Vector<double>& centreline_limits,
42  const Vector<double>& theta_positions,
43  const Vector<double>& radius_box,
44  const unsigned& nlayer,
45  TimeStepper* time_stepper_pt)
46  : Volume_pt(volume_pt)
47  {
48  // Mesh can only be built with 3D Qelements.
49  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
50 
51  // Build macro element-based domain
52  Domain_pt = new TubeDomain(
53  volume_pt, centreline_limits, theta_positions, radius_box, nlayer);
54 
55  // Set the number of boundaries
56  set_nboundary(3);
57 
58  // We have only bothered to parametrise boundary 1
59  Boundary_coordinate_exists[1] = true;
60 
61  // Allocate the store for the elements
62  unsigned nelem = 5 * nlayer;
63  Element_pt.resize(nelem);
64 
65  // Create dummy element so we can determine the number of nodes
66  ELEMENT* dummy_el_pt = new ELEMENT;
67 
68  // Read out the number of linear points in the element
69  unsigned n_p = dummy_el_pt->nnode_1d();
70 
71  // Kill the element
72  delete dummy_el_pt;
73 
74  // Can now allocate the store for the nodes
75  unsigned nnodes_total =
76  (n_p * n_p + 4 * (n_p - 1) * (n_p - 1)) * (1 + nlayer * (n_p - 1));
77  Node_pt.resize(nnodes_total);
78 
79  Vector<double> s(3);
80  Vector<double> r(3);
81 
82  // Storage for the intrinsic boundary coordinate
83  Vector<double> zeta(2);
84 
85  // Loop over elements and create all nodes
86  for (unsigned ielem = 0; ielem < nelem; ielem++)
87  {
88  // Create element
89  Element_pt[ielem] = new ELEMENT;
90 
91  // Loop over rows in z/s_2-direction
92  for (unsigned i2 = 0; i2 < n_p; i2++)
93  {
94  // Loop over rows in y/s_1-direction
95  for (unsigned i1 = 0; i1 < n_p; i1++)
96  {
97  // Loop over rows in x/s_0-direction
98  for (unsigned i0 = 0; i0 < n_p; i0++)
99  {
100  // Local node number
101  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
102 
103  // Create the node
104  Node* node_pt = finite_element_pt(ielem)->construct_node(
105  jnod_local, time_stepper_pt);
106 
107  // Set the position of the node from macro element mapping
108  s[0] = -1.0 + 2.0 * double(i0) / double(n_p - 1);
109  s[1] = -1.0 + 2.0 * double(i1) / double(n_p - 1);
110  s[2] = -1.0 + 2.0 * double(i2) / double(n_p - 1);
111  Domain_pt->macro_element_pt(ielem)->macro_map(s, r);
112 
113  node_pt->x(0) = r[0];
114  node_pt->x(1) = r[1];
115  node_pt->x(2) = r[2];
116  }
117  }
118  }
119  }
120 
121  // Initialise number of global nodes
122  unsigned node_count = 0;
123 
124  // Tolerance for node killing:
125  double node_kill_tol = 1.0e-12;
126 
127  // Check for error in node killing
128  bool stopit = false;
129 
130  // Define pine
131  const double pi = MathematicalConstants::Pi;
132 
133  // Loop over elements
134  for (unsigned ielem = 0; ielem < nelem; ielem++)
135  {
136  // Which layer?
137  unsigned ilayer = unsigned(ielem / 5);
138 
139  // Which macro element?
140  switch (ielem % 5)
141  {
142  // Macro element 0: Central box
143  //-----------------------------
144  case 0:
145 
146  // Loop over rows in z/s_2-direction
147  for (unsigned i2 = 0; i2 < n_p; i2++)
148  {
149  // Loop over rows in y/s_1-direction
150  for (unsigned i1 = 0; i1 < n_p; i1++)
151  {
152  // Loop over rows in x/s_0-direction
153  for (unsigned i0 = 0; i0 < n_p; i0++)
154  {
155  // Local node number
156  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
157 
158  // Has the node been killed?
159  bool killed = false;
160 
161  // First layer of all nodes in s_2 direction gets killed
162  // and re-directed to nodes in previous element layer
163  if ((i2 == 0) && (ilayer > 0))
164  {
165  // Neighbour element
166  unsigned ielem_neigh = ielem - 5;
167 
168  // Node in neighbour element
169  unsigned i0_neigh = i0;
170  unsigned i1_neigh = i1;
171  unsigned i2_neigh = n_p - 1;
172  unsigned jnod_local_neigh =
173  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
174 
175  // Check:
176  for (unsigned i = 0; i < 3; i++)
177  {
178  double error = std::fabs(
179  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
180  finite_element_pt(ielem_neigh)
181  ->node_pt(jnod_local_neigh)
182  ->x(i));
183  if (error > node_kill_tol)
184  {
185  oomph_info << "Error in node killing for i " << i << " "
186  << error << std::endl;
187  stopit = true;
188  }
189  }
190 
191  // Kill node
192  delete finite_element_pt(ielem)->node_pt(jnod_local);
193  killed = true;
194 
195  // Set pointer to neighbour:
196  finite_element_pt(ielem)->node_pt(jnod_local) =
197  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
198  }
199 
200  // No duplicate node: Copy across to mesh
201  if (!killed)
202  {
203  Node_pt[node_count] =
204  finite_element_pt(ielem)->node_pt(jnod_local);
205 
206  // Set boundaries:
207 
208  // Back: Boundary 0
209  if ((i2 == 0) && (ilayer == 0))
210  {
211  this->convert_to_boundary_node(Node_pt[node_count]);
212  add_boundary_node(0, Node_pt[node_count]);
213  }
214 
215  // Front: Boundary 2
216  if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
217  {
218  this->convert_to_boundary_node(Node_pt[node_count]);
219  add_boundary_node(2, Node_pt[node_count]);
220  }
221 
222  // Increment node counter
223  node_count++;
224  }
225  }
226  }
227  }
228 
229 
230  break;
231 
232  // Macro element 1: Bottom
233  //---------------------------------
234  case 1:
235 
236 
237  // Loop over rows in z/s_2-direction
238  for (unsigned i2 = 0; i2 < n_p; i2++)
239  {
240  // Loop over rows in y/s_1-direction
241  for (unsigned i1 = 0; i1 < n_p; i1++)
242  {
243  // Loop over rows in x/s_0-direction
244  for (unsigned i0 = 0; i0 < n_p; i0++)
245  {
246  // Local node number
247  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
248 
249  // Has the node been killed?
250  bool killed = false;
251 
252  // First layer of all nodes in s_2 direction gets killed
253  // and re-directed to nodes in previous element layer
254  if ((i2 == 0) && (ilayer > 0))
255  {
256  // Neighbour element
257  unsigned ielem_neigh = ielem - 5;
258 
259  // Node in neighbour element
260  unsigned i0_neigh = i0;
261  unsigned i1_neigh = i1;
262  unsigned i2_neigh = n_p - 1;
263  unsigned jnod_local_neigh =
264  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
265 
266 
267  // Check:
268  for (unsigned i = 0; i < 3; i++)
269  {
270  double error = std::fabs(
271  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
272  finite_element_pt(ielem_neigh)
273  ->node_pt(jnod_local_neigh)
274  ->x(i));
275  if (error > node_kill_tol)
276  {
277  oomph_info << "Error in node killing for i " << i << " "
278  << error << std::endl;
279  stopit = true;
280  }
281  }
282 
283  // Kill node
284  delete finite_element_pt(ielem)->node_pt(jnod_local);
285  killed = true;
286 
287  // Set pointer to neighbour:
288  finite_element_pt(ielem)->node_pt(jnod_local) =
289  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
290  }
291  // Not in first layer:
292  else
293  {
294  // Duplicate node: kill and set pointer to central element
295  if (i1 == (n_p - 1))
296  {
297  // Neighbour element
298  unsigned ielem_neigh = ielem - 1;
299 
300  // Node in neighbour element
301  unsigned i0_neigh = i0;
302  unsigned i1_neigh = 0;
303  unsigned i2_neigh = i2;
304  unsigned jnod_local_neigh =
305  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
306 
307  // Check:
308  for (unsigned i = 0; i < 3; i++)
309  {
310  double error = std::fabs(
311  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
312  finite_element_pt(ielem_neigh)
313  ->node_pt(jnod_local_neigh)
314  ->x(i));
315  if (error > node_kill_tol)
316  {
317  oomph_info << "Error in node killing for i " << i << " "
318  << error << std::endl;
319  stopit = true;
320  }
321  }
322 
323  // Kill node
324  delete finite_element_pt(ielem)->node_pt(jnod_local);
325  killed = true;
326 
327  // Set pointer to neighbour:
328  finite_element_pt(ielem)->node_pt(jnod_local) =
329  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
330  }
331  }
332 
333  // No duplicate node: Copy across to mesh
334  if (!killed)
335  {
336  Node_pt[node_count] =
337  finite_element_pt(ielem)->node_pt(jnod_local);
338 
339  // Set boundaries:
340 
341  // Back: Boundary 0
342  if ((i2 == 0) && (ilayer == 0))
343  {
344  this->convert_to_boundary_node(Node_pt[node_count]);
345  add_boundary_node(0, Node_pt[node_count]);
346  }
347 
348  // Front: Boundary 2
349  if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
350  {
351  this->convert_to_boundary_node(Node_pt[node_count]);
352  add_boundary_node(2, Node_pt[node_count]);
353  }
354 
355  // Bottom: outer boundary
356  if (i1 == 0)
357  {
358  this->convert_to_boundary_node(Node_pt[node_count]);
359  add_boundary_node(1, Node_pt[node_count]);
360 
361  // Get axial boundary coordinate
362  zeta[0] = centreline_limits[0] +
363  (double(ilayer) + double(i2) / double(n_p - 1)) *
364  (centreline_limits[1] - centreline_limits[0]) /
365  double(nlayer);
366 
367  // Get azimuthal boundary coordinate
368  zeta[1] = theta_positions[0] +
369  double(i1) / double(n_p - 1) * 0.5 *
370  (theta_positions[1] - theta_positions[0]);
371 
372  Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
373  }
374  // Increment node counter
375  node_count++;
376  }
377  }
378  }
379  } // End of loop over nodes
380 
381  break;
382 
383 
384  // Macro element 2: Right element
385  //--------------------------------
386  case 2:
387 
388  // Loop over rows in z/s_2-direction
389  for (unsigned i2 = 0; i2 < n_p; i2++)
390  {
391  // Loop over rows in y/s_1-direction
392  for (unsigned i1 = 0; i1 < n_p; i1++)
393  {
394  // Loop over rows in x/s_0-direction
395  for (unsigned i0 = 0; i0 < n_p; i0++)
396  {
397  // Local node number
398  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
399 
400  // Has the node been killed?
401  bool killed = false;
402 
403  // First layer of all nodes in s_2 direction gets killed
404  // and re-directed to nodes in previous element layer
405  if ((i2 == 0) && (ilayer > 0))
406  {
407  // Neighbour element
408  unsigned ielem_neigh = ielem - 5;
409 
410  // Node in neighbour element
411  unsigned i0_neigh = i0;
412  unsigned i1_neigh = i1;
413  unsigned i2_neigh = n_p - 1;
414  unsigned jnod_local_neigh =
415  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
416 
417  // Check:
418  for (unsigned i = 0; i < 3; i++)
419  {
420  double error = std::fabs(
421  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
422  finite_element_pt(ielem_neigh)
423  ->node_pt(jnod_local_neigh)
424  ->x(i));
425  if (error > node_kill_tol)
426  {
427  oomph_info << "Error in node killing for i " << i << " "
428  << error << std::endl;
429  stopit = true;
430  }
431  }
432 
433  // Kill node
434  delete finite_element_pt(ielem)->node_pt(jnod_local);
435  killed = true;
436 
437  // Set pointer to neighbour:
438  finite_element_pt(ielem)->node_pt(jnod_local) =
439  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
440  }
441  // Not in first layer:
442  else
443  {
444  // Duplicate node: kill and set pointer to previous element
445  if (i1 == 0)
446  {
447  // Neighbour element
448  unsigned ielem_neigh = ielem - 1;
449 
450  // Node in neighbour element
451  unsigned i0_neigh = n_p - 1;
452  unsigned i1_neigh = n_p - 1 - i0;
453  unsigned i2_neigh = i2;
454  unsigned jnod_local_neigh =
455  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
456 
457  // Check:
458  for (unsigned i = 0; i < 3; i++)
459  {
460  double error = std::fabs(
461  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
462  finite_element_pt(ielem_neigh)
463  ->node_pt(jnod_local_neigh)
464  ->x(i));
465  if (error > node_kill_tol)
466  {
467  oomph_info << "Error in node killing for i " << i << " "
468  << error << std::endl;
469  stopit = true;
470  }
471  }
472 
473  // Kill node
474  delete finite_element_pt(ielem)->node_pt(jnod_local);
475  killed = true;
476 
477  // Set pointer to neighbour:
478  finite_element_pt(ielem)->node_pt(jnod_local) =
479  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
480  }
481 
482 
483  // Duplicate node: kill and set pointer to central element
484  if ((i0 == 0) && (i1 != 0))
485  {
486  // Neighbour element
487  unsigned ielem_neigh = ielem - 2;
488 
489  // Node in neighbour element
490  unsigned i0_neigh = n_p - 1;
491  unsigned i1_neigh = i1;
492  unsigned i2_neigh = i2;
493  unsigned jnod_local_neigh =
494  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
495 
496  // Check:
497  for (unsigned i = 0; i < 3; i++)
498  {
499  double error = std::fabs(
500  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
501  finite_element_pt(ielem_neigh)
502  ->node_pt(jnod_local_neigh)
503  ->x(i));
504  if (error > node_kill_tol)
505  {
506  oomph_info << "Error in node killing for i " << i << " "
507  << error << std::endl;
508  stopit = true;
509  }
510  }
511 
512  // Kill node
513  delete finite_element_pt(ielem)->node_pt(jnod_local);
514  killed = true;
515 
516  // Set pointer to neighbour:
517  finite_element_pt(ielem)->node_pt(jnod_local) =
518  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
519  }
520 
521  // No duplicate node: Copy across to mesh
522  if (!killed)
523  {
524  Node_pt[node_count] =
525  finite_element_pt(ielem)->node_pt(jnod_local);
526 
527  // Set boundaries:
528 
529  // Back: Boundary 0
530  if ((i2 == 0) && (ilayer == 0))
531  {
532  this->convert_to_boundary_node(Node_pt[node_count]);
533  add_boundary_node(0, Node_pt[node_count]);
534  }
535 
536  // Front: Boundary 2
537  if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
538  {
539  this->convert_to_boundary_node(Node_pt[node_count]);
540  add_boundary_node(2, Node_pt[node_count]);
541  }
542 
543  // Tube boundary
544  if (i0 == n_p - 1)
545  {
546  this->convert_to_boundary_node(Node_pt[node_count]);
547  add_boundary_node(1, Node_pt[node_count]);
548 
549  // Get axial boundary coordinate
550  zeta[0] =
551  centreline_limits[0] +
552  (double(ilayer) + double(i2) / double(n_p - 1)) *
553  (centreline_limits[1] - centreline_limits[0]) /
554  double(nlayer);
555 
556  // Get azimuthal boundary coordinate
557  zeta[1] = theta_positions[1] +
558  double(i1) / double(n_p - 1) * 0.5 *
559  (theta_positions[2] - theta_positions[1]);
560 
561  Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
562  }
563 
564  // Increment node counter
565  node_count++;
566  }
567  }
568  }
569  }
570  }
571 
572  break;
573 
574  // Macro element 3: Top element
575  //--------------------------------
576  case 3:
577 
578  // Loop over rows in z/s_2-direction
579  for (unsigned i2 = 0; i2 < n_p; i2++)
580  {
581  // Loop over rows in y/s_1-direction
582  for (unsigned i1 = 0; i1 < n_p; i1++)
583  {
584  // Loop over rows in x/s_0-direction
585  for (unsigned i0 = 0; i0 < n_p; i0++)
586  {
587  // Local node number
588  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
589 
590  // Has the node been killed?
591  bool killed = false;
592 
593  // First layer of all nodes in s_2 direction gets killed
594  // and re-directed to nodes in previous element layer
595  if ((i2 == 0) && (ilayer > 0))
596  {
597  // Neighbour element
598  unsigned ielem_neigh = ielem - 5;
599 
600  // Node in neighbour element
601  unsigned i0_neigh = i0;
602  unsigned i1_neigh = i1;
603  unsigned i2_neigh = n_p - 1;
604  unsigned jnod_local_neigh =
605  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
606 
607  // Check:
608  for (unsigned i = 0; i < 3; i++)
609  {
610  double error = std::fabs(
611  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
612  finite_element_pt(ielem_neigh)
613  ->node_pt(jnod_local_neigh)
614  ->x(i));
615  if (error > node_kill_tol)
616  {
617  oomph_info << "Error in node killing for i " << i << " "
618  << error << std::endl;
619  stopit = true;
620  }
621  }
622 
623  // Kill node
624  delete finite_element_pt(ielem)->node_pt(jnod_local);
625  killed = true;
626 
627  // Set pointer to neighbour:
628  finite_element_pt(ielem)->node_pt(jnod_local) =
629  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
630  }
631  // Not in first layer:
632  else
633  {
634  // Duplicate node: kill and set pointer to previous element
635  if (i0 == n_p - 1)
636  {
637  // Neighbour element
638  unsigned ielem_neigh = ielem - 1;
639 
640  // Node in neighbour element
641  unsigned i0_neigh = i1;
642  unsigned i1_neigh = n_p - 1;
643  unsigned i2_neigh = i2;
644  unsigned jnod_local_neigh =
645  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
646 
647  // Check:
648  for (unsigned i = 0; i < 3; i++)
649  {
650  double error = std::fabs(
651  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
652  finite_element_pt(ielem_neigh)
653  ->node_pt(jnod_local_neigh)
654  ->x(i));
655  if (error > node_kill_tol)
656  {
657  oomph_info << "Error in node killing for i " << i << " "
658  << error << std::endl;
659  stopit = true;
660  }
661  }
662 
663  // Kill node
664  delete finite_element_pt(ielem)->node_pt(jnod_local);
665  killed = true;
666 
667  // Set pointer to neighbour:
668  finite_element_pt(ielem)->node_pt(jnod_local) =
669  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
670  }
671 
672  // Duplicate node: kill and set pointer to central element
673  if ((i1 == 0) && (i0 != n_p - 1))
674  {
675  // Neighbour element
676  unsigned ielem_neigh = ielem - 3;
677 
678  // Node in neighbour element
679  unsigned i0_neigh = i0;
680  unsigned i1_neigh = n_p - 1;
681  unsigned i2_neigh = i2;
682  unsigned jnod_local_neigh =
683  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
684 
685  // Check:
686  for (unsigned i = 0; i < 3; i++)
687  {
688  double error = std::fabs(
689  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
690  finite_element_pt(ielem_neigh)
691  ->node_pt(jnod_local_neigh)
692  ->x(i));
693  if (error > node_kill_tol)
694  {
695  oomph_info << "Error in node killing for i " << i << " "
696  << error << std::endl;
697  stopit = true;
698  }
699  }
700 
701  // Kill node
702  delete finite_element_pt(ielem)->node_pt(jnod_local);
703  killed = true;
704 
705  // Set pointer to neighbour:
706  finite_element_pt(ielem)->node_pt(jnod_local) =
707  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
708  }
709 
710  // No duplicate node: Copy across to mesh
711  if (!killed)
712  {
713  Node_pt[node_count] =
714  finite_element_pt(ielem)->node_pt(jnod_local);
715 
716  // Set boundaries:
717 
718  // Back: Boundary 0
719  if ((i2 == 0) && (ilayer == 0))
720  {
721  this->convert_to_boundary_node(Node_pt[node_count]);
722  add_boundary_node(0, Node_pt[node_count]);
723  }
724 
725  // Front: Boundary 2
726  if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
727  {
728  this->convert_to_boundary_node(Node_pt[node_count]);
729  add_boundary_node(2, Node_pt[node_count]);
730  }
731 
732  // Tube boundary
733  if (i1 == n_p - 1)
734  {
735  this->convert_to_boundary_node(Node_pt[node_count]);
736  add_boundary_node(1, Node_pt[node_count]);
737 
738  // Get axial boundary coordinate
739  zeta[0] =
740  centreline_limits[0] +
741  (double(ilayer) + double(i2) / double(n_p - 1)) *
742  (centreline_limits[1] - centreline_limits[0]) /
743  double(nlayer);
744 
745  // Get azimuthal boundary coordinate
746  zeta[1] = theta_positions[3] +
747  double(i0) / double(n_p - 1) * 0.5 *
748  (theta_positions[2] - theta_positions[3]);
749 
750  Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
751  }
752 
753  // Increment node counter
754  node_count++;
755  }
756  }
757  }
758  }
759  }
760  break;
761 
762 
763  // Macro element 4: Left element
764  //--------------------------------
765  case 4:
766 
767  // Loop over rows in z/s_2-direction
768  for (unsigned i2 = 0; i2 < n_p; i2++)
769  {
770  // Loop over rows in y/s_1-direction
771  for (unsigned i1 = 0; i1 < n_p; i1++)
772  {
773  // Loop over rows in x/s_0-direction
774  for (unsigned i0 = 0; i0 < n_p; i0++)
775  {
776  // Local node number
777  unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
778 
779  // Has the node been killed?
780  bool killed = false;
781 
782  // First layer of all nodes in s_2 direction gets killed
783  // and re-directed to nodes in previous element layer
784  if ((i2 == 0) && (ilayer > 0))
785  {
786  // Neighbour element
787  unsigned ielem_neigh = ielem - 5;
788 
789  // Node in neighbour element
790  unsigned i0_neigh = i0;
791  unsigned i1_neigh = i1;
792  unsigned i2_neigh = n_p - 1;
793  unsigned jnod_local_neigh =
794  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
795 
796  // Check:
797  for (unsigned i = 0; i < 3; i++)
798  {
799  double error = std::fabs(
800  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
801  finite_element_pt(ielem_neigh)
802  ->node_pt(jnod_local_neigh)
803  ->x(i));
804  if (error > node_kill_tol)
805  {
806  oomph_info << "Error in node killing for i " << i << " "
807  << error << std::endl;
808  stopit = true;
809  }
810  }
811 
812  // Kill node
813  delete finite_element_pt(ielem)->node_pt(jnod_local);
814  killed = true;
815 
816  // Set pointer to neighbour:
817  finite_element_pt(ielem)->node_pt(jnod_local) =
818  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
819  }
820  // Not in first layer:
821  else
822  {
823  // Duplicate node: kill and set pointer to previous element
824  if (i1 == n_p - 1)
825  {
826  // Neighbour element
827  unsigned ielem_neigh = ielem - 1;
828 
829  // Node in neighbour element
830  unsigned i0_neigh = 0;
831  unsigned i1_neigh = n_p - 1 - i0;
832  unsigned i2_neigh = i2;
833  unsigned jnod_local_neigh =
834  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
835 
836  // Check:
837  for (unsigned i = 0; i < 3; i++)
838  {
839  double error = std::fabs(
840  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
841  finite_element_pt(ielem_neigh)
842  ->node_pt(jnod_local_neigh)
843  ->x(i));
844  if (error > node_kill_tol)
845  {
846  oomph_info << "Error in node killing for i " << i << " "
847  << error << std::endl;
848  stopit = true;
849  }
850  }
851 
852  // Kill node
853  delete finite_element_pt(ielem)->node_pt(jnod_local);
854  killed = true;
855 
856  // Set pointer to neighbour:
857  finite_element_pt(ielem)->node_pt(jnod_local) =
858  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
859  }
860 
861  // Duplicate node: kill and set pointer to central element
862  if ((i0 == n_p - 1) && (i1 != n_p - 1))
863  {
864  // Neighbour element
865  unsigned ielem_neigh = ielem - 4;
866 
867  // Node in neighbour element
868  unsigned i0_neigh = 0;
869  unsigned i1_neigh = i1;
870  unsigned i2_neigh = i2;
871  unsigned jnod_local_neigh =
872  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
873 
874  // Check:
875  for (unsigned i = 0; i < 3; i++)
876  {
877  double error = std::fabs(
878  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
879  finite_element_pt(ielem_neigh)
880  ->node_pt(jnod_local_neigh)
881  ->x(i));
882  if (error > node_kill_tol)
883  {
884  oomph_info << "Error in node killing for i " << i << " "
885  << error << std::endl;
886  stopit = true;
887  }
888  }
889 
890  // Kill node
891  delete finite_element_pt(ielem)->node_pt(jnod_local);
892  killed = true;
893 
894  // Set pointer to neighbour:
895  finite_element_pt(ielem)->node_pt(jnod_local) =
896  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
897  }
898 
899  // Duplicate node: kill and set pointer to other ring element
900  if ((i1 == 0) && (i0 != n_p - 1))
901  {
902  // Neighbour element
903  unsigned ielem_neigh = ielem - 3;
904 
905  // Node in neighbour element
906  unsigned i0_neigh = 0;
907  unsigned i1_neigh = i0;
908  unsigned i2_neigh = i2;
909  unsigned jnod_local_neigh =
910  i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
911 
912  // Check:
913  for (unsigned i = 0; i < 3; i++)
914  {
915  double error = std::fabs(
916  finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
917  finite_element_pt(ielem_neigh)
918  ->node_pt(jnod_local_neigh)
919  ->x(i));
920  if (error > node_kill_tol)
921  {
922  oomph_info << "Error in node killing for i " << i << " "
923  << error << std::endl;
924  stopit = true;
925  }
926  }
927 
928  // Kill node
929  delete finite_element_pt(ielem)->node_pt(jnod_local);
930  killed = true;
931 
932  // Set pointer to neighbour:
933  finite_element_pt(ielem)->node_pt(jnod_local) =
934  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
935  }
936 
937 
938  // No duplicate node: Copy across to mesh
939  if (!killed)
940  {
941  Node_pt[node_count] =
942  finite_element_pt(ielem)->node_pt(jnod_local);
943 
944  // Set boundaries:
945 
946  // Back: Boundary 0
947  if ((i2 == 0) && (ilayer == 0))
948  {
949  this->convert_to_boundary_node(Node_pt[node_count]);
950  add_boundary_node(0, Node_pt[node_count]);
951  }
952 
953  // Front: Boundary 2
954  if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
955  {
956  this->convert_to_boundary_node(Node_pt[node_count]);
957  add_boundary_node(2, Node_pt[node_count]);
958  }
959 
960  // Tube boundary
961  if (i0 == 0)
962  {
963  this->convert_to_boundary_node(Node_pt[node_count]);
964  add_boundary_node(1, Node_pt[node_count]);
965 
966  // Get axial boundary coordinate
967  zeta[0] =
968  centreline_limits[0] +
969  (double(ilayer) + double(i2) / double(n_p - 1)) *
970  (centreline_limits[1] - centreline_limits[0]) /
971  double(nlayer);
972 
973  // Get azimuthal boundary coordinate
974  zeta[1] =
975  theta_positions[0] + 2.0 * pi +
976  double(i1) / double(n_p - 1) * 0.5 *
977  (theta_positions[3] - theta_positions[0] + 2.0 * pi);
978 
979  Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
980  }
981 
982  // Increment node counter
983  node_count++;
984  }
985  }
986  }
987  }
988  }
989  break;
990  }
991  }
992 
993  // Terminate if there's been an error
994  if (stopit)
995  {
996  std::ostringstream error_stream;
997  error_stream << "Error in killing nodes\n"
998  << "The most probable cause is that the domain is not\n"
999  << "compatible with the mesh.\n"
1000  << "For the TubeMesh, the domain must be\n"
1001  << "topologically consistent with a quarter tube with a\n"
1002  << "non-curved centreline.\n";
1003  throw OomphLibError(
1004  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1005  }
1006 
1007  // Setup boundary element lookup schemes
1008  setup_boundary_element_info();
1009  }
1010 
1011 } // namespace oomph
1012 #endif
oomph::TubeDomain
Tube as a domain. The entire domain must be defined by a GeomObject with the following convention: ze...
Definition: tube_domain.h:70
oomph::TubeMesh::volume_pt
GeomObject *& volume_pt()
Access function to GeomObject representing wall.
Definition: tube_mesh.template.h:80
oomph::TubeMesh::TubeMesh
TubeMesh(GeomObject *wall_pt, const Vector< double > &centreline_limits, const Vector< double > &theta_positions, const Vector< double > &radius_box, const unsigned &nlayer, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the volume, start and end coordinates fo...
Definition: tube_mesh.template.cc:40
oomph::TubeMesh::Domain_pt
TubeDomain * Domain_pt
Pointer to domain.
Definition: tube_mesh.template.h:99
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: annular_domain.h:35
tube_mesh.template.h