OpenLB 1.7
Loading...
Searching...
No Matches
indicatorBaseF3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2014-2016 Cyril Masquelier, Albert Mink, Mathias J. Krause, Benjamin Förster
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef INDICATOR_BASE_F_3D_HH
25#define INDICATOR_BASE_F_3D_HH
26
27
28#include<cmath>
29#include "indicatorBaseF3D.h"
31#include "math.h"
32
33#ifndef M_PI
34#define M_PI 3.14159265358979323846
35#endif
36#define M_PI2 1.57079632679489661923
37
38namespace olb {
39
40template <typename S>
43
44template <typename S>
46{
47 return _myMin;
50template <typename S>
52{
53 return _myMax;
56template <typename S>
57bool IndicatorF3D<S>::distance(S& distance, const Vector<S,3>& origin, S precision, const Vector<S,3>& direction)
58{
59 return util::distance(distance, origin, direction, precision,
60 [&](const Vector<S,3> input) {
61 return this->signedDistance(input);
62 },
63 [&](const Vector<S,3>& pos) {
64 return isInsideBox(pos);
65 });
67
68template <typename S>
69bool IndicatorF3D<S>::distance(S& distance, const Vector<S,3>& origin, const Vector<S,3>& direction, S precision, S pitch)
71 return util::distance(distance, origin, direction, precision, pitch,
72 [&](bool output[1], const S input[3]) {
73 return (*this)(output, input);
74 },
75 [&](const Vector<S,3>& pos) {
76 return isInsideBox(pos);
77 });
78}
79
80template <typename S>
81bool IndicatorF3D<S>::distance(S& distance, const Vector<S,3>& origin, const Vector<S,3>& direction, int iC)
82{
83 S precision = .0001;
84 S pitch = 0.5;
85
86 return this->distance(distance, origin, direction, precision, pitch);
87}
88
89template <typename S>
90bool IndicatorF3D<S>::distance(S& distance, const Vector<S,3>& origin)
91{
92 S input[3] = {origin[0],origin[1],origin[2]};
93 return this->distance(distance, input);
94}
95
96template <typename S>
97bool IndicatorF3D<S>::distance(S& distance, const S input[])
98{
99 distance = util::abs(signedDistance(input));
100 return true;
101}
102
103template <typename S>
104bool IndicatorF3D<S>::rotOnAxis(Vector<S,3>& vec_rot, const Vector<S,3>& vec, const Vector<S,3>& axis, S& theta)
105{
107
108 //Vector<S,3> axisN(axis*(1/(axis).norm())); // normalize rotation axis
109 Vector<S,3> axisN(axis*(1/norm(axis))); // normalize rotation axis
110
111 vec_rot = vec ;
112
113 Vector<S,3> crossProd;
114
115 crossProd[0] = axisN[1]*vec[2] - axisN[2]*vec[1];
116 crossProd[1] = axisN[2]*vec[0] - axisN[0]*vec[2];
117 crossProd[2] = axisN[0]*vec[1] - axisN[1]*vec[0];
118
119 S dotProd = axisN[0]*vec[0] + axisN[1]*vec[1] + axisN[2]*vec[2];
120
121 //v_rot = util::cos(theta)*vec + (crossProd)*util::sin(theta) + axisN*(dotProduct3D(axisN,vec))*(1 - util::cos(theta));
122 vec_rot = util::cos(theta)*vec + (crossProd)*util::sin(theta) + axisN*(dotProd)*(1 - util::cos(theta));
123
124 return true;
125
126}
127
128template <typename S>
129bool IndicatorF3D<S>::normal(Vector<S,3>& normal, const Vector<S,3>& origin, const Vector<S,3>& direction, int iC)
130{
131 //OstreamManager clout(std::cout,"IndicatorF3D");
132#ifdef OLB_DEBUG
133 std::cout << "Calculating IndicatorF3D Normal" << std::endl;
134#endif
135
136 bool originValue;
137 (*this)(&originValue, origin.data());
138 Vector<S,3> currentPoint(origin);
139
140 S precision = .0001;
141
142 S dist;
143 distance(dist, origin, direction, iC);
144
145 S dirMag = norm(direction);
146#ifdef OLB_DEBUG
147 std::cout << "magnitude = " << dirMag << std::endl;
148#endif
149
150 //Vector<S,3> POS(origin + dist*direction*(1/const_cast<Vector<S,3>&> (direction).norm())); //Point on Surface
151 //direction = direction*(1/const_cast<Vector<S,3>&> (direction).norm());
152 Vector<S,3> directionN(direction*(1/dirMag));
153 Vector<S,3> POS(origin + dist*directionN); //Point on Surface
154
156 Vector<S,3> directionPerp;
157 if ( (util::nearZero(directionN[0]) && util::nearZero(directionN[1]) && !util::nearZero(directionN[2]))
158 || (util::nearZero(directionN[0]) && !util::nearZero(directionN[1]) && util::nearZero(directionN[2])) ) {
159 directionPerp[0] = 1;
160 directionPerp[1] = 0;
161 directionPerp[2] = 0;
162 }
163 else if ( !util::nearZero(directionN[0]) && util::nearZero(directionN[1]) && util::nearZero(directionN[2]) ) {
164 directionPerp[0] = 0;
165 directionPerp[1] = 0;
166 directionPerp[2] = 1;
167 }
168 else if ( ( !util::nearZero(directionN[0]) || !util::nearZero(directionN[1]) ) && !util::nearZero(directionN[2]) ) {
169 directionPerp[0] = directionN[0];
170 directionPerp[1] = directionN[1];
171 directionPerp[2] = -(directionN[0] + directionN[1])/directionN[2];
172 }
173 else if ( ( !util::nearZero(directionN[0]) || !util::nearZero(directionN[1]) ) && util::nearZero(directionN[2]) ) {
174 directionPerp[0] = directionN[0];
175 directionPerp[1] = -(directionN[0] + directionN[2])/directionN[1];
176 directionPerp[2] = directionN[2];
177 }
178 else {
179 std::cout << "Error: unknown case for perpendicular check" << std::endl;
180 return false;
181 }
182
183 Vector<S,3> directionPerpN(directionPerp*(dirMag/norm(directionPerp)));
184
185
186 Vector<S,3> point1;
187 Vector<S,3> point2;
188 Vector<S,3> point3;
189
190 bool currentValue;
191
194
195
196 for (int n: {
197 0,120,240
198 }) {
199 S thetaMain = util::degreeToRadian(n);
200
202 Vector<S,3> perp;
203 rotOnAxis(perp, directionPerpN, directionN, thetaMain);
204 Vector<S,3> perpPoint(POS + perp);
205
206 //std::cout << "perp = [" << perp[0] << "," << perp[1] << "," << perp[2] << "]" << std::endl;
207
208 //S rotate = 90.;
209 S rotate = 179.;
210 S pitch = rotate/2.;
211
212 Vector<S,3> vec(perp);
213
214 Vector<S,3> rotAxis;
215
216 //rotAxis = cross(perp,directionN);
217 rotAxis[0] = perp[1]*directionN[2] - perp[2]*directionN[1];
218 rotAxis[1] = perp[2]*directionN[0] - perp[0]*directionN[2];
219 rotAxis[2] = perp[0]*directionN[1] - perp[1]*directionN[0];
220
221 //normal = rotAxis;
222 //return true;
223 //std::cout << "rotAxis = [" << rotAxis[0] << "," << rotAxis[1] << "," << rotAxis[2] << "]" << std::endl;
224
226 Vector<S,3> testPOS;
227 S testAngle(0.25*M_PI);
228 rotOnAxis(testPOS, perp, rotAxis, testAngle);
229 Vector<S,3> testPoint( POS + testPOS);
230 //std::cout << "testPOS = [" << testPOS[0] << "," << testPOS[1] << "," << testPOS[2] << "]" << std::endl;
231 //std::cout << "testPoint = [" << testPoint[0] << "," << testPoint[1] << "," << testPoint[2] << "]" << std::endl;
232
233 S distTestPoint = norm(testPoint - origin);
234 S distPerpPoint = norm(perpPoint - origin);
235
236 S mod = 0;
237 if (distTestPoint < distPerpPoint) { // pos. angle rotates towards
238 mod = -1;
239 }
240 else {
241 mod = 1;
242 }
243
244 while (util::abs(pitch) >= precision) {
245
246 S theta(util::degreeToRadian(pitch));
247
248 currentPoint = POS + vec;
249 (*this)(&currentValue, currentPoint.data());
250
251 S temp;
252 if (currentValue == originValue) {
253 temp = mod*theta;
254 rotOnAxis(vec, vec, rotAxis, temp);
255 }
256 else {
257 temp = -mod*theta;
258 rotOnAxis(vec, vec, rotAxis, temp);
259 }
260 pitch /= 2.;
261
262
263
264 }
265
266 if (n == 0) {
267 point1 = currentPoint;
268 }
269 else if (n == 120) {
270 point2 = currentPoint;
271 }
272 else if (n == 240) {
273 point3 = currentPoint;
274 }
275 else {
276 std::cout << "Something broke" << std::endl;
277 return false;
278 }
279
280 }
281
283 Vector<S,3> vec1 (point1 - point2);
284 Vector<S,3> vec2 (point1 - point3);
285
286 normal[0] = -(vec1[1]*vec2[2] - vec1[2]*vec2[1]);
287 normal[1] = -(vec1[2]*vec2[0] - vec1[0]*vec2[2]);
288 normal[2] = -(vec1[0]*vec2[1] - vec1[1]*vec2[0]);
289
290 normalize(normal);
291
292
293 //S dist;
294 //Vector<S,3> dist;
295 //distance(dist, origin, direction, iC);
296 //normal = Vector<S,3>(dist);
297 //normal = POS;
298 //normal = directionPerpN;
299 //normal = directionN;
300 return true;
301
302}
303
304template <typename S>
306{
307 return point >= _myMin && point <= _myMax;
308}
309
310template <typename S>
312{
313 // TODO: Add fallback algorithm here
314 assert(false);
315 return std::numeric_limits<double>::quiet_NaN();
316}
317
318template <typename S>
320{
321 return util::surfaceNormal(pos, meshSize, [&](const Vector<S,3>& pos) {
322 return this->signedDistance(pos);
323 });
324}
325
326template <typename S>
328 std::function<Vector<S,3>(const Vector<S,3>&)> transformPos)
329{
330 return this->surfaceNormal(transformPos(pos), meshSize);
331}
332
333template <typename S>
334bool IndicatorF3D<S>::operator() (bool output[1], const S input[3])
335{
336 output[0] = this->signedDistance(input) <= 0;
337 return output[0];
338}
339
340template <typename S>
342 : _f(f)
343{
344 this->_myMin = _f->getMin();
345 this->_myMax = _f->getMax();
346}
347
348template <typename S>
349bool IndicatorIdentity3D<S>::operator() (bool output[1], const S input[3])
350{
351 return (_f)->operator()(output, input);
352}
353
354
355} // namespace olb
356
357#endif
#define M_PI
GenericF is a base class, that can represent continuous as well as discrete functions.
Definition genericF.h:50
IndicatorF3D is an application from .
virtual bool operator()(bool output[1], const S input[3])
Returns true if input is inside the indicator.
virtual bool distance(S &distance, const Vector< S, 3 > &origin, S precision, const Vector< S, 3 > &direction)
virtual Vector< S, 3 > & getMax()
Vector< S, 3 > _myMin
virtual Vector< S, 3 > & getMin()
virtual bool normal(Vector< S, 3 > &normal, const Vector< S, 3 > &origin, const Vector< S, 3 > &direction, int iC=-1)
returns true and the normal if there was one found for an given origin and direction
bool isInsideBox(Vector< S, 3 > point)
Returns true if point is inside a cube with corners _myMin and _myMax
virtual S signedDistance(const Vector< S, 3 > &input)
Returns signed distance to the nearest point on the indicator surface.
virtual Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize)
Return surface normal.
virtual bool rotOnAxis(Vector< S, 3 > &vec_rot, const Vector< S, 3 > &vec, const Vector< S, 3 > &axis, S &theta)
Rotate vector around axis by angle theta.
Vector< S, 3 > _myMax
IndicatorIdentity3D(std::shared_ptr< IndicatorF3D< S > > f)
bool operator()(bool output[1], const S input[3])
Returns true if input is inside the indicator.
std::shared_ptr< IndicatorF3D< S > > _f
Plain old scalar vector.
Definition vector.h:47
constexpr const T * data() const any_platform
Definition vector.h:161
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
bool distance(S &distance, const Vector< S, D > &origin, const Vector< S, D > &direction, S precision, S pitch, F1 isInside, F2 isInsideBoundingBox)
Vector< S, D > surfaceNormal(const Vector< S, D > &pos, const S meshSize, F1 sdf)
ADf< T, DIM > sin(const ADf< T, DIM > &a)
Definition aDiff.h:569
decltype(Vector< decltype(util::sqrt(T())), D >()) degreeToRadian(const Vector< T, D > &angle)
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
Top level namespace for all of OpenLB.
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245