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