annular_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_ANNULAR_MESH_TEMPLATE_CC
27 #define OOMPH_ANNULAR_MESH_TEMPLATE_CC
28 
29 #include "annular_mesh.template.h"
30 
31 namespace oomph
32 {
33  //===================================================================
34  /// Helper function to Wrap mesh into annular shape
35  //==================================================================
36  template<class ELEMENT>
38  const double& a,
39  const double& h,
40  const double& azimuthal_fraction,
41  const double& phi)
42  {
43  // Create the hole
44  Ellipse ellipse(a, a);
45 
46  // Set all the initial positions of the nodes
47  Vector<double> xi(1);
48  Vector<double> base(2);
49  Vector<double> N(2);
50  const unsigned n_node = this->nnode();
51  for (unsigned n = 0; n < n_node; n++)
52  {
53  // Pointer to node
54  Node* nod_pt = this->node_pt(n);
55 
56  // Get the angle of the node -- rotate such that jump in angle
57  // appears at periodic boundary. Shrink domain slightly
58  // to keep angle unique
59  xi[0] = (1.0 - 1.0e-10) * (-azimuthal_fraction * nod_pt->x(0)) * 2.0 *
60  MathematicalConstants::Pi +
61  MathematicalConstants::Pi -
62  (1.0 - azimuthal_fraction) * 2.0 * MathematicalConstants::Pi;
63 
64  // Rotate
65  xi[0] += phi;
66 
67  // Get the node's fraction in the radial direction
68  double w = nod_pt->x(1);
69 
70  // Get the position on the ellipse base
71  ellipse.position(xi, base);
72 
73  // Get the unit normal, if it were a circle , by normalising the base
74  double norm = sqrt(base[0] * base[0] + base[1] * base[1]);
75  N[0] = base[0] / norm;
76  N[1] = base[1] / norm;
77 
78  // Set circular film from the ellipse
79  nod_pt->x(0) = base[0] + w * (h + a - norm) * N[0];
80  nod_pt->x(1) = base[1] + w * (h + a - norm) * N[1];
81 
82  // Set boundary coordinates
83  Vector<double> xi_bound(1);
84 
85  // Polar angle for boundary coordinate on boundary 0
86  if (nod_pt->is_on_boundary(0))
87  {
88  xi_bound[0] = atan2(nod_pt->x(1), nod_pt->x(0));
89  nod_pt->set_coordinates_on_boundary(0, xi_bound);
90  }
91 
92  // Radius for boundary coordinate on boundary 1
93  if (nod_pt->is_on_boundary(1))
94  {
95  xi_bound[0] = sqrt(pow(nod_pt->x(0), 2) + pow(nod_pt->x(1), 2));
96  nod_pt->set_coordinates_on_boundary(1, xi_bound);
97  }
98 
99  // Polar angle for boundary coordinate on boundary 2
100  if (nod_pt->is_on_boundary(2))
101  {
102  xi_bound[0] = atan2(nod_pt->x(1), nod_pt->x(0));
103  nod_pt->set_coordinates_on_boundary(2, xi_bound);
104  }
105 
106  // Radius for boundary coordinate on boundary 3
107  if (nod_pt->is_on_boundary(3))
108  {
109  xi_bound[0] = sqrt(pow(nod_pt->x(0), 2) + pow(nod_pt->x(1), 2));
110  nod_pt->set_coordinates_on_boundary(3, xi_bound);
111  }
112  }
113 
114  this->Boundary_coordinate_exists[0] = true;
115  this->Boundary_coordinate_exists[1] = true;
116  this->Boundary_coordinate_exists[2] = true;
117  this->Boundary_coordinate_exists[3] = true;
118  }
119 
120 
121 } // namespace oomph
122 
123 
124 #endif
void wrap_into_annular_shape(const double &a, const double &h, const double &azimuthal_fraction, const double &phi)
Wrap mesh into annular shape.
////////////////////////////////////////////////////////////////////// //////////////////////////////...