24#ifndef INDICATOR_BASE_F_3D_HH
25#define INDICATOR_BASE_F_3D_HH
34#define M_PI 3.14159265358979323846
36#define M_PI2 1.57079632679489661923
61 return this->signedDistance(input);
64 return isInsideBox(pos);
71 return util::distance(distance, origin, direction, precision, pitch,
72 [&](
bool output[1],
const S input[3]) {
73 return (*
this)(output, input);
76 return isInsideBox(pos);
86 return this->distance(distance, origin, direction, precision, pitch);
92 S input[3] = {origin[0],origin[1],origin[2]};
93 return this->distance(distance, input);
99 distance =
util::abs(signedDistance(input));
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];
119 S dotProd = axisN[0]*vec[0] + axisN[1]*vec[1] + axisN[2]*vec[2];
133 std::cout <<
"Calculating IndicatorF3D Normal" << std::endl;
137 (*this)(&originValue, origin.
data());
143 distance(dist, origin, direction, iC);
145 S dirMag =
norm(direction);
147 std::cout <<
"magnitude = " << dirMag << std::endl;
159 directionPerp[0] = 1;
160 directionPerp[1] = 0;
161 directionPerp[2] = 0;
164 directionPerp[0] = 0;
165 directionPerp[1] = 0;
166 directionPerp[2] = 1;
169 directionPerp[0] = directionN[0];
170 directionPerp[1] = directionN[1];
171 directionPerp[2] = -(directionN[0] + directionN[1])/directionN[2];
174 directionPerp[0] = directionN[0];
175 directionPerp[1] = -(directionN[0] + directionN[2])/directionN[1];
176 directionPerp[2] = directionN[2];
179 std::cout <<
"Error: unknown case for perpendicular check" << std::endl;
183 Vector<S,3> directionPerpN(directionPerp*(dirMag/
norm(directionPerp)));
203 rotOnAxis(perp, directionPerpN, directionN, thetaMain);
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];
227 S testAngle(0.25*
M_PI);
228 rotOnAxis(testPOS, perp, rotAxis, testAngle);
233 S distTestPoint =
norm(testPoint - origin);
234 S distPerpPoint =
norm(perpPoint - origin);
237 if (distTestPoint < distPerpPoint) {
248 currentPoint = POS + vec;
249 (*this)(¤tValue, currentPoint.data());
252 if (currentValue == originValue) {
254 rotOnAxis(vec, vec, rotAxis, temp);
258 rotOnAxis(vec, vec, rotAxis, temp);
267 point1 = currentPoint;
270 point2 = currentPoint;
273 point3 = currentPoint;
276 std::cout <<
"Something broke" << std::endl;
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]);
307 return point >= _myMin && point <= _myMax;
315 return std::numeric_limits<double>::quiet_NaN();
322 return this->signedDistance(pos);
330 return this->surfaceNormal(transformPos(pos), meshSize);
336 output[0] = this->signedDistance(input) <= 0;
351 return (_f)->operator()(output, input);
GenericF is a base class, that can represent continuous as well as discrete functions.
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()
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.
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
constexpr const T * data() const any_platform
ADf< T, DIM > abs(const ADf< T, DIM > &a)
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)
decltype(Vector< decltype(util::sqrt(T())), D >()) degreeToRadian(const Vector< T, D > &angle)
bool nearZero(const ADf< T, DIM > &a)
ADf< T, DIM > cos(const ADf< T, DIM > &a)
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})