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