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-2022 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
30
31namespace 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);
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 *
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
////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:644
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta.
Definition: geom_objects.h:745
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary,...
Definition: nodes.h:1373
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
void wrap_into_annular_shape(const double &a, const double &h, const double &azimuthal_fraction, const double &phi)
Wrap mesh into annular shape.
const double Pi
50 digits from maple
//////////////////////////////////////////////////////////////////// ////////////////////////////////...