Nektar++
RefRegionCylinder.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: RefRegionCylinder.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Cylindrical surface for the refinement region.
32//
33////////////////////////////////////////////////////////////////////////////////
34
36
38{
39
40RefRegionCylinder::RefRegionCylinder(const unsigned int coordim,
41 NekDouble radius,
42 std::vector<NekDouble> coord1,
43 std::vector<NekDouble> coord2,
44 std::vector<unsigned int> numModes,
45 std::vector<unsigned int> numPoints)
46 : RefRegion(coordim, radius, coord1, coord2, numModes, numPoints)
47{
48}
49
50/**
51 * @brief Check if vertex is inside the cylinder.
52 *
53 * @param coords coordinates of the vertex
54 * @return true or false depending on if the vertex is inside
55 * or not of the surface defined by the user.
56 */
58{
59 const size_t dim = coords.size();
60 Array<OneD, NekDouble> e(dim, 0.0); // direction: rb - ra
61 Array<OneD, NekDouble> m(dim, 0.0); // momentum: ra x rb
62 NekDouble d = 0.0; // distance
63
64 // Direction
65 e[0] = m_coord2[0] - m_coord1[0];
66 e[1] = m_coord2[1] - m_coord1[1];
67 e[2] = m_coord2[2] - m_coord1[2];
68
69 // Cross product of vectors 'coord1' and 'coord2'
70 m[0] = m_coord1[1] * m_coord2[2] - m_coord1[2] * m_coord2[1];
71 m[1] = m_coord1[2] * m_coord2[0] - m_coord1[0] * m_coord2[2];
72 m[2] = m_coord1[0] * m_coord2[1] - m_coord1[1] * m_coord2[0];
73
74 // 1. Distance of P (coords) to line AB (coord1coord2) is equal or less than
75 // R
76 // d = || e x (rp - ra) || / || e ||
77
78 // rA - rP
79 Array<OneD, NekDouble> rpa(dim, 0.0);
80 rpa[0] = coords[0] - m_coord1[0];
81 rpa[1] = coords[1] - m_coord1[1];
82 rpa[2] = coords[2] - m_coord1[2];
83
84 // || e ||
85 NekDouble e_mod = sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
86
87 // || e x (rp - ra) ||
88 Array<OneD, NekDouble> exrpa(dim, 0.0);
89 exrpa[0] = e[1] * rpa[2] - e[2] * rpa[1];
90 exrpa[1] = e[2] * rpa[0] - e[0] * rpa[2];
91 exrpa[2] = e[0] * rpa[1] - e[1] * rpa[0];
92
93 NekDouble exrpa_mod =
94 sqrt(exrpa[0] * exrpa[0] + exrpa[1] * exrpa[1] + exrpa[2] * exrpa[2]);
95
96 d = exrpa_mod / e_mod;
97 if (d >= m_radius)
98 {
99 return false;
100 }
101
102 // 2. Closest point Q on line AB to P
103 // rq = rp + (e x (m + e x rp)) / ||e||^{2}
104
105 // (m + e x rp)
106 Array<OneD, NekDouble> mpexrp(dim, 0.0);
107 mpexrp[0] = m[0] + e[1] * coords[2] - e[2] * coords[1];
108 mpexrp[1] = m[1] + e[2] * coords[0] - e[0] * coords[2];
109 mpexrp[2] = m[2] + e[0] * coords[1] - e[1] * coords[0];
110
111 // e x (m + e x rp) = numerator
112 Array<OneD, NekDouble> numerator(dim, 0.0);
113 numerator[0] = e[1] * mpexrp[2] - e[2] * mpexrp[1];
114 numerator[1] = e[2] * mpexrp[0] - e[0] * mpexrp[2];
115 numerator[2] = e[0] * mpexrp[1] - e[1] * mpexrp[0];
116
117 // rq
118 Array<OneD, NekDouble> rq(dim, 0.0);
119 rq[0] = coords[0] + (numerator[0] / pow(e_mod, 2));
120 rq[1] = coords[1] + (numerator[1] / pow(e_mod, 2));
121 rq[2] = coords[2] + (numerator[2] / pow(e_mod, 2));
122
123 // 3. The baricentric coordinates of Q(wa,wb) such that rq = wa*ra+wb*rb
124 // wa = ||rq x rb||/||m||, wb = ||rq x ra||/||m||
125 NekDouble wa = 0.0;
126 NekDouble wb = 0.0;
127
128 // ||m||
129 NekDouble m_mod = sqrt(m[0] * m[0] + m[1] * m[1] + m[2] * m[2]);
130 // m_mod > 0: condition
131 ASSERTL0(m_mod, "The cylinder axis must not go through the origin");
132
133 // ||rq x rb||
134 Array<OneD, NekDouble> rqxrb(dim, 0.0);
135 rqxrb[0] = rq[1] * m_coord2[2] - rq[2] * m_coord2[1];
136 rqxrb[1] = rq[2] * m_coord2[0] - rq[0] * m_coord2[2];
137 rqxrb[2] = rq[0] * m_coord2[1] - rq[1] * m_coord2[0];
138 NekDouble rqxrb_mod =
139 sqrt(rqxrb[0] * rqxrb[0] + rqxrb[1] * rqxrb[1] + rqxrb[2] * rqxrb[2]);
140
141 // ||rq x ra||
142 Array<OneD, NekDouble> rqxra(dim, 0.0);
143 rqxra[0] = rq[1] * m_coord1[2] - rq[2] * m_coord1[1];
144 rqxra[1] = rq[2] * m_coord1[0] - rq[0] * m_coord1[2];
145 rqxra[2] = rq[0] * m_coord1[1] - rq[1] * m_coord1[0];
146 NekDouble rqxra_mod =
147 sqrt(rqxra[0] * rqxra[0] + rqxra[1] * rqxra[1] + rqxra[2] * rqxra[2]);
148
149 // wa
150 wa = rqxrb_mod / m_mod;
151 // wb
152 wb = rqxra_mod / m_mod;
153
154 // 4. Check that point Q between A and B by making sure the baricentric
155 // coordinates are between 0 and 1.
156 if ((wa >= 0) && (wa <= 1) && (wb >= 0) && (wb <= 1))
157 {
158 return true;
159 }
160 else
161 {
162 return false;
163 }
164}
165
166} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
bool v_Contains(const Array< OneD, NekDouble > &coords) override
Check if vertex is inside the surface region.
RefRegionCylinder(const unsigned int coordim, NekDouble radius, std::vector< NekDouble > coord1, std::vector< NekDouble > coord2, std::vector< unsigned int > numModes, std::vector< unsigned int > numPoints)
Constructor.
Abstract base class for the refinement surface region.
Definition: RefRegion.h:51
std::vector< NekDouble > m_coord2
Coordinate 2.
Definition: RefRegion.h:83
std::vector< NekDouble > m_coord1
Coordinate 1.
Definition: RefRegion.h:81
NekDouble m_radius
Radius of the surface region.
Definition: RefRegion.h:89
std::vector< double > d(NPUPPER *NPUPPER)
double NekDouble
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285