full_circle_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 
31 #ifndef OOMPH_FULL_CIRCLE_MESH_TEMPLATE_CC
32 #define OOMPH_FULL_CIRCLE_MESH_TEMPLATE_CC
33 
34 #include "full_circle_mesh.template.h"
35 
36 namespace oomph
37 {
38 
39 //====================================================================
40 /// \short Constructor for deformable quarter tube mesh class.
41 /// The domain is specified by the GeomObject that
42 /// identifies the entire volume.
43 //====================================================================
44 template <class ELEMENT>
45 FullCircleMesh<ELEMENT>::FullCircleMesh(GeomObject* area_pt,
46  const Vector<double> &theta_positions,
47  const Vector<double> &radius_box,
48  TimeStepper* time_stepper_pt) :
49  Area_pt(area_pt)
50 {
51 
52 // Check that the vectors are the correct sizes.
53 #ifdef PARANOID
54  if(radius_box.size() != 4 || theta_positions.size() != 4)
55  {
56  std::string err = "This mesh is hard coded to only work for the case when there are 5 elements: the central square and 4 surrounding elements, but you gave vectors inconsistent with this.";
57  throw OomphLibError(err, OOMPH_CURRENT_FUNCTION,
58  OOMPH_EXCEPTION_LOCATION);
59  }
60 #endif
61 
62  // Mesh can only be built with 2D Qelements.
63  MeshChecker::assert_geometric_element<QElementGeometricBase,ELEMENT>(2);
64 
65  // Build macro element-based domain
66  Domain_pt = new FullCircleDomain(area_pt,theta_positions,radius_box);
67 
68  //Set the number of boundaries
69  set_nboundary(1);
70 
71  //We have only bothered to parametrise the only boundary (boundary 0)
72  Boundary_coordinate_exists[0] = true;
73 
74  // Allocate the store for the elements
75  const unsigned nelem=5;
76  Element_pt.resize(nelem);
77 
78  // Create dummy element so we can determine the number of nodes
79  ELEMENT* dummy_el_pt=new ELEMENT;
80 
81  // Read out the number of linear points in the element
82  unsigned n_p = dummy_el_pt->nnode_1d();
83 
84  // Kill the element
85  delete dummy_el_pt;
86 
87  // Can now allocate the store for the nodes
88  unsigned nnodes_total=
89  (n_p*n_p + 4*(n_p-1)*(n_p -1));
90  Node_pt.resize(nnodes_total);
91 
92  Vector<double> s(2);
93  Vector<double> r(2);
94 
95  //Storage for the intrinsic boundary coordinate
96  Vector<double> zeta(1);
97 
98  // Loop over elements and create all nodes
99  for (unsigned ielem=0;ielem<nelem;ielem++)
100  {
101  // Create element
102  Element_pt[ielem] = new ELEMENT;
103 
104  // Loop over rows in y/s_1-direction
105  for (unsigned i1=0;i1<n_p;i1++)
106  {
107  // Loop over rows in x/s_0-direction
108  for (unsigned i0=0;i0<n_p;i0++)
109  {
110  // Local node number
111  unsigned jnod_local=i0+i1*n_p;
112 
113  // Create the node
114  Node* node_pt=finite_element_pt(ielem)->
115  construct_node(jnod_local,time_stepper_pt);
116 
117  //Set the position of the node from macro element mapping
118  s[0]=-1.0+2.0*double(i0)/double(n_p-1);
119  s[1]=-1.0+2.0*double(i1)/double(n_p-1);
120  Domain_pt->macro_element_pt(ielem)->macro_map(s,r);
121 
122  node_pt->x(0) = r[0];
123  node_pt->x(1) = r[1];
124  }
125  }
126  }
127 
128  // Initialise number of global nodes
129  unsigned node_count=0;
130 
131  // Tolerance for node killing:
132  double node_kill_tol=1.0e-12;
133 
134  // Check for error in node killing
135  bool stopit=false;
136 
137  // Define pine
138  const double pi = MathematicalConstants::Pi;
139 
140  // Loop over elements
141  for (unsigned ielem=0;ielem<nelem;ielem++)
142  {
143  // Which macro element?
144  switch(ielem)
145  {
146  // Macro element 0: Central box, create all the nodes
147  //----------------------------------------------------
148  case 0:
149 
150  // Loop over rows in y/s_1-direction
151  for (unsigned i1=0;i1<n_p;i1++)
152  {
153 
154  // Loop over rows in x/s_0-direction
155  for (unsigned i0=0;i0<n_p;i0++)
156  {
157 
158  // Local node number
159  unsigned jnod_local=i0+i1*n_p;
160 
161  //Add the node to the mesh
162  Node_pt[node_count]=finite_element_pt(ielem)->node_pt(jnod_local);
163 
164  // Increment node counter
165  node_count++;
166  }
167  }
168 
169  break;
170 
171  // Macro element 1: Bottom
172  //---------------------------------
173  case 1:
174 
175  // Loop over rows in y/s_1-direction
176  for (unsigned i1=0;i1<n_p;i1++)
177  {
178 
179  // Loop over rows in x/s_0-direction
180  for (unsigned i0=0;i0<n_p;i0++)
181  {
182 
183  // Local node number
184  unsigned jnod_local=i0+i1*n_p;
185 
186  // Has the node been killed?
187  bool killed=false;
188 
189  // Duplicate node: kill and set pointer to central element
190  if (i1==(n_p-1))
191  {
192  // Neighbour element
193  unsigned ielem_neigh=ielem-1;
194 
195  // Node in neighbour element
196  unsigned i0_neigh=i0;
197  unsigned i1_neigh=0;
198  unsigned jnod_local_neigh = i0_neigh + i1_neigh*n_p;
199 
200  // Check:
201  for (unsigned i=0;i<2;i++)
202  {
203  double error=std::fabs(
204  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
205  finite_element_pt(ielem_neigh)->
206  node_pt(jnod_local_neigh)->x(i));
207  if (error>node_kill_tol)
208  {
209  oomph_info << "Error in node killing for i "
210  << i << " " << error << std::endl;
211  stopit=true;
212  }
213  }
214 
215  // Kill node
216  delete finite_element_pt(ielem)->node_pt(jnod_local);
217  killed=true;
218 
219  // Set pointer to neighbour:
220  finite_element_pt(ielem)->node_pt(jnod_local)=
221  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
222  }
223 
224  // No duplicate node: Copy across to mesh
225  if (!killed)
226  {
227  Node_pt[node_count]=
228  finite_element_pt(ielem)->node_pt(jnod_local);
229 
230  // Set boundaries:
231 
232  // Bottom: outer boundary
233  if (i1==0)
234  {
235  this->convert_to_boundary_node(Node_pt[node_count]);
236  add_boundary_node(0,Node_pt[node_count]);
237 
238  // Get azimuthal boundary coordinate
239  zeta[0]= theta_positions[0] +
240  double(i1)/double(n_p-1)*0.5*
241  (theta_positions[1]-theta_positions[0]);
242 
243  Node_pt[node_count]->set_coordinates_on_boundary(0,zeta);
244  }
245 
246  // Increment node counter
247  node_count++;
248  }
249 
250  }
251  } //End of loop over nodes
252 
253  break;
254 
255 
256  // Macro element 2: Right element
257  //--------------------------------
258  case 2:
259 
260  // Loop over rows in y/s_1-direction
261  for (unsigned i1=0;i1<n_p;i1++)
262  {
263 
264  // Loop over rows in x/s_0-direction
265  for (unsigned i0=0;i0<n_p;i0++)
266  {
267 
268  // Local node number
269  unsigned jnod_local=i0+i1*n_p;
270 
271  // Has the node been killed?
272  bool killed=false;
273 
274  // Duplicate node: kill and set pointer to previous element
275  if (i1==0)
276  {
277 
278  // Neighbour element
279  unsigned ielem_neigh=ielem-1;
280 
281  // Node in neighbour element
282  unsigned i0_neigh=n_p-1;
283  unsigned i1_neigh= n_p-1 - i0;
284 
285  unsigned jnod_local_neigh =i0_neigh+i1_neigh*n_p;
286 
287  // Check:
288  for (unsigned i=0;i<2;i++)
289  {
290  double error=std::fabs(
291  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
292  finite_element_pt(ielem_neigh)->
293  node_pt(jnod_local_neigh)->x(i));
294  if (error>node_kill_tol)
295  {
296  oomph_info << "Error in node killing for i "
297  << i << " " << error << std::endl;
298  stopit=true;
299  }
300  }
301 
302  // Kill node
303  delete finite_element_pt(ielem)->node_pt(jnod_local);
304  killed=true;
305 
306  // Set pointer to neighbour:
307  finite_element_pt(ielem)->node_pt(jnod_local)=
308  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
309 
310  }
311 
312 
313  // Duplicate node: kill and set pointer to central element
314  if ((i0==0)&&(i1!=0))
315  {
316 
317  // Neighbour element
318  unsigned ielem_neigh=ielem-2;
319 
320  // Node in neighbour element
321  unsigned i0_neigh=n_p-1;
322  unsigned i1_neigh=i1;
323  unsigned jnod_local_neigh = i0_neigh+i1_neigh*n_p;
324 
325  // Check:
326  for (unsigned i=0;i<2;i++)
327  {
328  double error=std::fabs(
329  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
330  finite_element_pt(ielem_neigh)->
331  node_pt(jnod_local_neigh)->x(i));
332  if (error>node_kill_tol)
333  {
334  oomph_info << "Error in node killing for i "
335  << i << " " << error << std::endl;
336  stopit=true;
337  }
338  }
339 
340  // Kill node
341  delete finite_element_pt(ielem)->node_pt(jnod_local);
342  killed=true;
343 
344  // Set pointer to neighbour:
345  finite_element_pt(ielem)->node_pt(jnod_local)=
346  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
347 
348  }
349 
350  // No duplicate node: Copy across to mesh
351  if (!killed)
352  {
353  Node_pt[node_count]=finite_element_pt(ielem)->
354  node_pt(jnod_local);
355 
356  // Set boundaries:
357 
358  // FullCircle boundary
359  if (i0==n_p-1)
360  {
361  this->convert_to_boundary_node(Node_pt[node_count]);
362  add_boundary_node(0,Node_pt[node_count]);
363 
364  // Get azimuthal boundary coordinate
365  zeta[0]=theta_positions[1] +
366  double(i1)/double(n_p-1)*0.5*
367  (theta_positions[2]-theta_positions[1]);
368 
369  Node_pt[node_count]->set_coordinates_on_boundary(0,zeta);
370 
371  }
372 
373  // Increment node counter
374  node_count++;
375  }
376 
377  }
378  }
379 
380 
381  break;
382 
383  // Macro element 3: Top element
384  //--------------------------------
385  case 3:
386 
387  // Loop over rows in y/s_1-direction
388  for (unsigned i1=0;i1<n_p;i1++)
389  {
390 
391  // Loop over rows in x/s_0-direction
392  for (unsigned i0=0;i0<n_p;i0++)
393  {
394 
395  // Local node number
396  unsigned jnod_local=i0+i1*n_p;
397 
398  // Has the node been killed?
399  bool killed=false;
400 
401 
402  // Duplicate node: kill and set pointer to previous element
403  if (i0==n_p-1)
404  {
405 
406  // Neighbour element
407  unsigned ielem_neigh=ielem-1;
408 
409  // Node in neighbour element
410  unsigned i0_neigh=i1;
411  unsigned i1_neigh= n_p-1;
412  unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p;
413 
414  // Check:
415  for (unsigned i=0;i<2;i++)
416  {
417  double error=std::fabs(
418  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
419  finite_element_pt(ielem_neigh)->
420  node_pt(jnod_local_neigh)->x(i));
421  if (error>node_kill_tol)
422  {
423  oomph_info << "Error in node killing for i "
424  << i << " " << error << std::endl;
425  stopit=true;
426  }
427  }
428 
429  // Kill node
430  delete finite_element_pt(ielem)->node_pt(jnod_local);
431  killed=true;
432 
433  // Set pointer to neighbour:
434  finite_element_pt(ielem)->node_pt(jnod_local)=
435  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
436 
437  }
438 
439  // Duplicate node: kill and set pointer to central element
440  if ((i1==0)&&(i0!=n_p-1))
441  {
442 
443  // Neighbour element
444  unsigned ielem_neigh=ielem-3;
445 
446  // Node in neighbour element
447  unsigned i0_neigh=i0;
448  unsigned i1_neigh=n_p-1;
449  unsigned jnod_local_neigh= i0_neigh+i1_neigh*n_p;
450 
451  // Check:
452  for (unsigned i=0;i<2;i++)
453  {
454  double error=std::fabs(
455  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
456  finite_element_pt(ielem_neigh)->
457  node_pt(jnod_local_neigh)->x(i));
458  if (error>node_kill_tol)
459  {
460  oomph_info << "Error in node killing for i "
461  << i << " " << error << std::endl;
462  stopit=true;
463  }
464  }
465 
466  // Kill node
467  delete finite_element_pt(ielem)->node_pt(jnod_local);
468  killed=true;
469 
470  // Set pointer to neighbour:
471  finite_element_pt(ielem)->node_pt(jnod_local)=
472  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
473 
474  }
475 
476  // No duplicate node: Copy across to mesh
477  if (!killed)
478  {
479  Node_pt[node_count]=finite_element_pt(ielem)->
480  node_pt(jnod_local);
481 
482  // Set boundaries:
483 
484  // FullCircle boundary
485  if (i1==n_p-1)
486  {
487  this->convert_to_boundary_node(Node_pt[node_count]);
488  add_boundary_node(0,Node_pt[node_count]);
489 
490  // Get azimuthal boundary coordinate
491  zeta[0]=theta_positions[3] +
492  double(i0)/double(n_p-1)*0.5*
493  (theta_positions[2]-theta_positions[3]);
494 
495  Node_pt[node_count]->set_coordinates_on_boundary(0,zeta);
496 
497  }
498 
499  // Increment node counter
500  node_count++;
501  }
502 
503  }
504  }
505  break;
506 
507 
508  // Macro element 4: Left element
509  //--------------------------------
510  case 4:
511 
512  // Loop over rows in y/s_1-direction
513  for (unsigned i1=0;i1<n_p;i1++)
514  {
515 
516  // Loop over rows in x/s_0-direction
517  for (unsigned i0=0;i0<n_p;i0++)
518  {
519 
520  // Local node number
521  unsigned jnod_local=i0+i1*n_p;
522 
523  // Has the node been killed?
524  bool killed=false;
525 
526  // Duplicate node: kill and set pointer to previous element
527  if (i1==n_p-1)
528  {
529 
530  // Neighbour element
531  unsigned ielem_neigh=ielem-1;
532 
533  // Node in neighbour element
534  unsigned i0_neigh=0;
535  unsigned i1_neigh= n_p - 1 -i0;
536  unsigned jnod_local_neigh=i0_neigh+i1_neigh*n_p;
537 
538  // Check:
539  for (unsigned i=0;i<2;i++)
540  {
541  double error=std::fabs(
542  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
543  finite_element_pt(ielem_neigh)->
544  node_pt(jnod_local_neigh)->x(i));
545  if (error>node_kill_tol)
546  {
547  oomph_info << "Error in node killing for i "
548  << i << " " << error << std::endl;
549  stopit=true;
550  }
551  }
552 
553  // Kill node
554  delete finite_element_pt(ielem)->node_pt(jnod_local);
555  killed=true;
556 
557  // Set pointer to neighbour:
558  finite_element_pt(ielem)->node_pt(jnod_local)=
559  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
560 
561  }
562 
563  // Duplicate node: kill and set pointer to central element
564  if ((i0==n_p-1)&&(i1!=n_p-1))
565  {
566 
567  // Neighbour element
568  unsigned ielem_neigh=ielem-4;
569 
570  // Node in neighbour element
571  unsigned i0_neigh=0;
572  unsigned i1_neigh=i1;
573  unsigned jnod_local_neigh = i0_neigh+i1_neigh*n_p;
574 
575  // Check:
576  for (unsigned i=0;i<2;i++)
577  {
578  double error=std::fabs(
579  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
580  finite_element_pt(ielem_neigh)->
581  node_pt(jnod_local_neigh)->x(i));
582  if (error>node_kill_tol)
583  {
584  oomph_info << "Error in node killing for i "
585  << i << " " << error << std::endl;
586  stopit=true;
587  }
588  }
589 
590  // Kill node
591  delete finite_element_pt(ielem)->node_pt(jnod_local);
592  killed=true;
593 
594  // Set pointer to neighbour:
595  finite_element_pt(ielem)->node_pt(jnod_local)=
596  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
597 
598  }
599 
600  // Duplicate node: kill and set pointer to other ring element
601  if ((i1==0)&&(i0!=n_p-1))
602  {
603 
604  // Neighbour element
605  unsigned ielem_neigh=ielem-3;
606 
607  // Node in neighbour element
608  unsigned i0_neigh=0;
609  unsigned i1_neigh=i0;
610  unsigned jnod_local_neigh = i0_neigh+i1_neigh*n_p;
611 
612  // Check:
613  for (unsigned i=0;i<2;i++)
614  {
615  double error=std::fabs(
616  finite_element_pt(ielem)->node_pt(jnod_local)->x(i)-
617  finite_element_pt(ielem_neigh)->
618  node_pt(jnod_local_neigh)->x(i));
619  if (error>node_kill_tol)
620  {
621  oomph_info << "Error in node killing for i "
622  << i << " " << error << std::endl;
623  stopit=true;
624  }
625  }
626 
627  // Kill node
628  delete finite_element_pt(ielem)->node_pt(jnod_local);
629  killed=true;
630 
631  // Set pointer to neighbour:
632  finite_element_pt(ielem)->node_pt(jnod_local)=
633  finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
634 
635  }
636 
637 
638  // No duplicate node: Copy across to mesh
639  if (!killed)
640  {
641  Node_pt[node_count]=finite_element_pt(ielem)->
642  node_pt(jnod_local);
643 
644  // Set boundaries:
645 
646  // FullCircle boundary
647  if (i0==0)
648  {
649  this->convert_to_boundary_node(Node_pt[node_count]);
650  add_boundary_node(0,Node_pt[node_count]);
651 
652  // Get azimuthal boundary coordinate
653  zeta[0]=theta_positions[0] + 2.0*pi +
654  double(i1)/double(n_p-1)*0.5*
655  (theta_positions[3]-theta_positions[0] + 2.0*pi);
656 
657  Node_pt[node_count]->set_coordinates_on_boundary(0,zeta);
658  }
659 
660  // Increment node counter
661  node_count++;
662 
663  }
664  }
665  }
666  break;
667 
668  }
669  }
670 
671  // Terminate if there's been an error
672  if (stopit)
673  {
674  std::ostringstream error_stream;
675  error_stream << "Error in killing nodes\n"
676  << "The most probable cause is that the domain is not\n"
677  << "compatible with the mesh.\n"
678  << "For the FullCircleMesh, the domain must be\n"
679  << "topologically consistent with a quarter tube with a\n"
680  << "non-curved centreline.\n";
681  throw OomphLibError(error_stream.str(),
682  OOMPH_CURRENT_FUNCTION,
683  OOMPH_EXCEPTION_LOCATION);
684  }
685 
686  // Setup boundary element lookup schemes
687  setup_boundary_element_info();
688 
689 }
690 
691 }
692 #endif
oomph::FullCircleDomain
Topologically circular domain, e.g. a tube cross section. The entire domain must be defined by a Geom...
Definition: full_circle_domain.h:73
oomph::FullCircleMesh::Domain_pt
FullCircleDomain * Domain_pt
Pointer to domain.
Definition: full_circle_mesh.template.h:87
full_circle_mesh.template.h
oomph::FullCircleMesh::FullCircleMesh
FullCircleMesh(GeomObject *wall_pt, const Vector< double > &theta_positions, const Vector< double > &radius_box, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the area; values of theta at which divid...
Definition: full_circle_mesh.template.cc:45