24#ifndef INDICATOR_F_3D_HH
25#define INDICATOR_F_3D_HH
41 : _translate(translate), _indicator(indicator), _myMin(indicator.getMin()), _myMax(indicator.getMax())
43 this->_myMin[0] += this->_translate[0];
44 this->_myMin[1] += this->_translate[1];
45 this->_myMin[2] += this->_translate[2];
47 this->_myMax[0] += this->_translate[0];
48 this->_myMax[1] += this->_translate[1];
49 this->_myMax[2] += this->_translate[2];
55 const S inputTranslated[3] = {
56 input[0] - _translate[0],
57 input[1] - _translate[1],
58 input[2] - _translate[2]
61 _indicator.operator()(output, inputTranslated);
69 return _indicator.signedDistance(translate);
88 : _center(center), _normal(normal), _radius2(radius*radius),
89 _cylinder(center, normal, radius, 5.0*std::numeric_limits<S>::epsilon())
92 this->
_myMin = _center - radius;
93 this->
_myMax = _center + radius;
98 S normal0, S normal1, S normal2, S radius)
107 return _cylinder(output,input);
132 : _center(center), _radius(radius), _radius2(_radius*_radius)
134 this->
_myMin = _center - radius;
135 this->
_myMax = _center + radius;
141 this->_myMin = sphere._myMin;
142 this->_myMax = sphere._myMax;
143 _center = sphere._center;
144 _radius = sphere._radius;
145 _radius2 = sphere._radius2;
172 S a = direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2];
182 S b = 2.*((origin[0] - _center[0])*direction[0] +
183 (origin[1] - _center[1])*direction[1] +
184 (origin[2] - _center[2])*direction[2])/a;
185 S c = -_radius2 + (origin[0] - _center[0])*(origin[0] - _center[0])
186 + (origin[1] - _center[1])*(origin[1] - _center[1])
187 + (origin[2] - _center[2])*(origin[2] - _center[2]);
199 if ((x1<0.) || (x2<0.)) {
221 : _indicatorF(std::move(indicatorF)), _layerSize(layerSize)
223 this->
_myMin = indicatorF->getMin() - layerSize;
224 this->
_myMax = indicatorF->getMax() + layerSize;
233 for (
int iX =- 1; iX < 2; ++iX) {
234 for (
int iY =- 1; iY < 2; ++iY) {
235 for (
int iZ =- 1; iZ < 2; ++iZ) {
236 r[0] = input[0] + iX*_layerSize;
237 r[1] = input[1] + iY*_layerSize;
238 r[2] = input[2] + iZ*_layerSize;
239 _indicatorF(output,r);
255 return sdf::rounding(_indicatorF->signedDistance(input), _layerSize);
262 : _indicatorF(indicatorF), _layerSize(layerSize)
273 _indicatorF(output,input);
279 for (
int iX =- 1; iX < 2; ++iX) {
280 for (
int iY =- 1; iY < 2; ++iY) {
281 for (
int iZ =- 1; iZ < 2; ++iZ) {
282 r[0] = input[0] + iX*_layerSize;
283 r[1] = input[1] + iY*_layerSize;
284 r[2] = input[2] + iZ*_layerSize;
285 _indicatorF(output,r);
300 : _center1(center1), _center2(center2),
301 _ba(_center2 - _center1),
302 _baba(_ba[0]*_ba[0] + _ba[1]*_ba[1] + _ba[2]*_ba[2]),
303 _radius2(radius*radius), _length(util::sqrt(_baba))
314 : _center1(center1-.5*eps*
normalize(normal)),
315 _center2(_center1 + eps*
normalize(normal)),
316 _ba(_center2 - _center1),
317 _baba(_ba[0]*_ba[0] + _ba[1]*_ba[1] + _ba[2]*_ba[2]),
318 _radius2(radius*radius), _length(util::sqrt(_baba))
328 : _center1(circleF.getCenter() - .5*eps*circleF.getNormal()),
329 _center2(_center1 + eps*circleF.getNormal()),
330 _ba(_center2 - _center1),
331 _baba(_ba[0]*_ba[0] + _ba[1]*_ba[1] + _ba[2]*_ba[2]),
332 _radius2(circleF.getRadius()*circleF.getRadius()), _length(util::sqrt(_baba))
345 S X = _I[0]*pa[0] + _I[1]*pa[1] + _I[2]*pa[2];
346 S Y = _J[0]*pa[0] + _J[1]*pa[1] + _J[2]*pa[2];
347 S Z = _K[0]*pa[0] + _K[1]*pa[1] + _K[2]*pa[2];
350 output[0] = ( Z <= _length && Z >= 0 && X*X + Y*Y <= _radius2 );
363 clout <<
"Warning: in the cylinder, the two centers have the same coordinates" << std::endl;
364 clout << _center1 << std::endl;
365 clout << _center2 << std::endl;
371 S normi =
util::sqrt (_K[1]*_K[1]+_K[0]*_K[0]);
372 _I = {-_K[1]/normi, _K[0]/normi,0};
373 _J = {_K[1]*_I[2] - _K[2]*_I[1], _K[2]*_I[0] - _K[0]*_I[2], _K[0]*_I[1] - _K[1]*_I[0]};
378 S maxx= _center1[0] +
util::sqrt(_I[0]*_I[0]+_J[0]*_J[0])*r +
util::max(_K[0]*_length, 0.);
379 S minx= _center1[0] -
util::sqrt(_I[0]*_I[0]+_J[0]*_J[0])*r +
util::min(_K[0]*_length, 0.);
381 S maxy= _center1[1] +
util::sqrt(_I[1]*_I[1]+_J[1]*_J[1])*r +
util::max(_K[1]*_length, 0.);
382 S miny= _center1[1] -
util::sqrt(_I[1]*_I[1]+_J[1]*_J[1])*r +
util::min(_K[1]*_length, 0.);
384 S maxz= _center1[2] +
util::sqrt(_I[2]*_I[2]+_J[2]*_J[2])*r +
util::max(_K[2]*_length, 0.);
385 S minz= _center1[2] -
util::sqrt(_I[2]*_I[2]+_J[2]*_J[2])*r +
util::min(_K[2]*_length, 0.);
387 this->_myMin = {minx, miny, minz};
388 this->_myMax = {maxx, maxy, maxz};
421 auto axis = _center2 - _center1;
422 auto axisPoint = _center1 + (_center2 - _center1) * randomness();
428 const S theta = randomness() * 2 *
M_PI;
430 return hyperplane.project({r *
util::cos(theta),
440 S radius1, S radius2)
441 : _center1(center1), _center2(center2),
442 _ba(center2 - center1), _baba(_ba*_ba),
443 _radius1(radius1), _radius2(radius2),
444 _length(util::sqrt((_center2 - _center1)*(_center2 - _center1)))
447 _K = (_center2 - _center1)/_length;
453 clout <<
"Warning: in the cone, the two centers have the same coordinates" << std::endl;
454 clout << _center1 << std::endl;
455 clout << _center2 << std::endl;
461 S normi =
util::sqrt(_K[1]*_K[1] + _K[0]*_K[0]);
462 _I = {-_K[1]/normi, _K[0]/normi,0};
463 _J = {_K[1]*_I[2] - _K[2]*_I[1], _K[2]*_I[0] - _K[0]*_I[2], _K[0]*_I[1] - _K[1]*_I[0]};
467 util::sqrt(_I[0]*_I[0]+_J[0]*_J[0])*_radius2 + _K[0]*_length);
469 -
util::sqrt(_I[0]*_I[0]+_J[0]*_J[0])*_radius2 + _K[0]*_length);
472 util::sqrt(_I[1]*_I[1]+_J[1]*_J[1])*_radius2 + _K[1]*_length);
474 -
util::sqrt(_I[1]*_I[1]+_J[1]*_J[1])*_radius2 + _K[1]*_length);
477 util::sqrt(_I[2]*_I[2]+_J[2]*_J[2])*_radius2 + _K[2]*_length);
479 -
util::sqrt(_I[2]*_I[2]+_J[2]*_J[2])*_radius2 + _K[2]*_length);
481 this->
_myMin = {minx, miny, minz};
482 this->
_myMax = {maxx, maxy, maxz};
491 S X = _I[0]*pa[0] + _I[1]*pa[1] + _I[2]*pa[2];
492 S Y = _J[0]*pa[0] + _J[1]*pa[1] + _J[2]*pa[2];
493 S Z = _K[0]*pa[0] + _K[1]*pa[1] + _K[2]*pa[2];
494 S radius = _radius1 + (_radius2 - _radius1)*Z / _length;
496 output[0] = ( Z <= _length && Z >= 0 && X*X + Y*Y <= radius*radius );
527 return sdf::cone(input, _center1, _ba, _baba, _radius1, _radius2);
532 : _center(center), _radius(radius)
536 clout <<
"WARNING: The indicator doesn't provide an exact signed distance function." << std::endl;
539 assert(_radius[0]>0 && _radius[1]>0 && _radius[2]>0);
541 this->
_myMin = center - radius;
542 this->
_myMax = center + radius;
566 : _center(center), _xHalfAxis(xHalfAxis), _yHalfAxis(yHalfAxis), _zHalfAxis(zHalfAxis), _exp1(exponent1), _exp2(exponent2)
568 assert(_xHalfAxis>0 && _yHalfAxis>0 && _zHalfAxis>0);
627 if ( (ab+c) <= 1. ) {
641 : _center(origin+.5*extend),_xLength(extend[0]), _yLength(extend[1]), _zLength(extend[2])
643 assert(_xLength>0 && _yLength>0 && _zLength>0);
646 this->
_myMax = origin + extend;
651 : _center(center), _xLength(xLength), _yLength(yLength), _zLength(zLength)
653 assert(_xLength>0 && _yLength>0 && _zLength>0);
655 this->
_myMin = {_center[0] - S{0.5}*_xLength, _center[1] - S{0.5}*_yLength, _center[2] - S{0.5}*_zLength};
656 this->
_myMax = {_center[0] + S{0.5}*_xLength, _center[1] + S{0.5}*_yLength, _center[2] + S{0.5}*_zLength};
688 Vector<S,3> eps = std::numeric_limits<BaseType<S>>::epsilon();
689 output[0] = ( q < eps );
696 return { S(
util::abs(S(input[0] - _center[0])) - 0.5 * _xLength),
697 S(
util::abs(S(input[1] - _center[1])) - 0.5 * _yLength),
698 S(
util::abs(S(input[2] - _center[2])) - 0.5 * _zLength)
713 auto it = std::minmax_element(v.
begin(), v.
end());
714 int min_idx = std::distance(v.
begin(), it.first);
716 axis[min_idx] = v[min_idx];
717 auto axisPoint = _center + axis * randomness();
722 int max_idx = std::distance(v.
begin(), it.second);
723 const S width = v[max_idx] * randomness();
724 int i = 3 - max_idx - min_idx;
725 const S height = v[i] * randomness();
727 return hyperplane.project({width,height});
734 :
IndicatorCuboid3D<S>(extend, origin), _theta(theta), _plane(plane), _centerRotation(centerRotation)
736 assert(plane==0 || plane==1 || plane==2);
741 :
IndicatorCuboid3D<S>(xLength, yLength, zLength, origin), _theta(theta), _plane(plane), _centerRotation(centerRotation)
743 assert(plane==0 || plane==1 || plane==2);
756 else if (_plane == 0) {
761 S x = input[i] - _centerRotation[i];
762 S y = input[j] - _centerRotation[j];
767 x =
xx + _centerRotation[i];
768 y =
yy + _centerRotation[j];
769 newInput[_plane] = input[_plane];
779 transformInput(input, newInput);
787 S newInput[3], tmp[3] = {input[0], input[1], input[2]};
788 transformInput(tmp, newInput);
806 std::stringstream xmlCenter1( params.
getAttribute(
"center") );
807 xmlCenter1 >> center[0] >> center[1] >> center[2];
808 std::stringstream xmlCenter2( params.
getAttribute(
"normal") );
809 xmlCenter2 >> normal[0] >> normal[1] >> normal[2];
810 std::stringstream xmlRadius( params.
getAttribute(
"radius") );
813 return std::make_shared<IndicatorCircle3D<S>>(center, normal, radius);
827 std::stringstream xmlCenter1( params.
getAttribute(
"center") );
828 xmlCenter1 >> center[0] >> center[1] >> center[2];
829 std::stringstream xmlRadius( params.
getAttribute(
"radius") );
832 return std::make_shared<IndicatorSphere3D<S>>(center, radius);
848 std::stringstream xmlCenter1( (params).getAttribute(
"center1") );
849 xmlCenter1 >> center1[0] >> center1[1] >> center1[2];
850 std::stringstream xmlCenter2( (params).getAttribute(
"center2") );
851 xmlCenter2 >> center2[0] >> center2[1] >> center2[2];
852 std::stringstream xmlRadius( (params).getAttribute(
"radius") );
860 return std::make_shared<IndicatorCylinder3D<S>>(center1, center2, radius);
876 std::stringstream xmlCenter1( params.
getAttribute(
"center1") );
877 xmlCenter1 >> center1[0] >> center1[1] >> center1[2];
878 std::stringstream xmlCenter2( params.
getAttribute(
"center2") );
879 xmlCenter2 >> center2[0] >> center2[1] >> center2[2];
880 std::stringstream xmlRadius1( params.
getAttribute(
"radius1") );
881 xmlRadius1 >> radius1;
882 std::stringstream xmlRadius2( params.
getAttribute(
"radius2") );
883 xmlRadius2 >> radius2;
885 return std::make_shared<IndicatorCone3D<S>>(center1, center2, radius1, radius2);
899 std::stringstream xmlOrigin( params.
getAttribute(
"origin") );
900 xmlOrigin >> origin[0] >> origin[1] >> origin[2];
901 std::stringstream xmlExtend( params.
getAttribute(
"extend") );
902 xmlExtend >> extend[0] >> extend[1] >> extend[2];
904 return std::make_shared<IndicatorCuboid3D<S>>(extend, origin);
916 std::shared_ptr<IndicatorF3D<S>> output = createIndicatorF3D<S>(**params.
begin());
917 for (
auto it = params.
begin()+1; it != params.
end(); ++it) {
918 output = output + createIndicatorF3D<S>(**it);
932 std::shared_ptr<IndicatorF3D<S>> output = createIndicatorF3D<S>(**params.
begin());
933 for (
auto it = params.
begin()+1; it != params.
end(); ++it) {
934 output = output - createIndicatorF3D<S>(**it);
948 std::shared_ptr<IndicatorF3D<S>> output = createIndicatorF3D<S>(**params.
begin());
949 for (
auto it = params.
begin()+1; it != params.
end(); ++it) {
950 output = output * createIndicatorF3D<S>(**it);
964 std::string actualName = params.
getName();
965 if ( actualName ==
"IndicatorCircle3D" ) {
966 return createIndicatorCircle3D<S>(params);
968 else if ( actualName ==
"IndicatorSphere3D" ) {
969 return createIndicatorSphere3D<S>(params);
971 else if ( actualName ==
"IndicatorCylinder3D" ) {
972 return createIndicatorCylinder3D<S>(params);
974 else if ( actualName ==
"IndicatorCone3D" ) {
975 return createIndicatorCone3D<S>(params);
977 else if ( actualName ==
"IndicatorCuboid3D" ) {
978 return createIndicatorCuboid3D<S>(params);
980 else if ( actualName ==
"IndicatorUnion3D" ) {
981 return createIndicatorUnion3D<S>(params);
983 else if ( actualName ==
"IndicatorWithout3D" ) {
984 return createIndicatorWithout3D<S>(params);
986 else if ( actualName ==
"IndicatorIntersection3D" ) {
987 return createIndicatorIntersection3D<S>(params);
990 auto firstChild = params.
begin();
991 return createIndicatorF3D<S>( **firstChild );
1000template <
typename T>
1003 output[0] = _f(input) <= 0.0;
Smart pointer for managing the various ways of passing functors around.
indicator function for a 3D circle
Vector< S, 3 > const & getCenter() const
IndicatorCircle3D(Vector< S, 3 > center, Vector< S, 3 > normal, S radius)
Vector< S, 3 > const & getNormal() const
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
Vector< S, 3 > const & getCenter1() const
IndicatorCone3D(Vector< S, 3 > center1, Vector< S, 3 > center2, S radius1, S radius2=0)
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
Vector< S, 3 > const & getCenter2() const
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
indicator function for a 3d-cuboid, parallel to the planes x=0, y=0, z=0.
Vector< S, 3 > const & getCenter() const
S const getxLength() const
IndicatorCuboid3D(Vector< S, 3 > extend, Vector< S, 3 > origin)
constructs an cuboid with x axis from origin[0] to origin[0]+extend[0], ...
S const getyLength() const
S const getzLength() const
Vector< S, 3 > getSample(const std::function< S()> &randomness) const override
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
bool operator()(bool output[], const S input[]) override
returns true if input is inside, otherwise false
indicator function for a 3d-cuboid, turned by an angle theta around an axis of rotation
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
bool operator()(bool output[], const S input[])
returns true if input is inside, otherwise false
IndicatorCuboidRotate3D(Vector< S, 3 > extend, Vector< S, 3 > origin, S theta, int plane, Vector< S, 3 > centerRotation)
indicator function for a 3d-cylinder
IndicatorCylinder3D(Vector< S, 3 > center1, Vector< S, 3 > center2, S radius)
Indicator function for a cylinder.
Vector< S, 3 > const & getCenter1() const
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
Vector< S, 3 > getSample(const std::function< S()> &randomness) const override
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
Vector< S, 3 > const & getCenter2() const
Vector< S, 3 > const & getCenter() const
Vector< S, 3 > const & getRadius() const
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
IndicatorEllipsoid3D(Vector< S, 3 > center, Vector< S, 3 > radius)
IndicatorF3D is an application from .
virtual Vector< S, 3 > & getMax()
virtual Vector< S, 3 > & getMin()
IndicatorInternal3D(IndicatorF3D< S > &indicatorF, S layerSize)
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
IndicatorLayer3D(FunctorPtr< IndicatorF3D< S > > &&indicatorF, S layerSize)
bool operator()(bool output[], const T input[]) override
IndicatorSDF3D(std::function< T(Vector< T, 3 >)> f)
indicator function for a 3D-sphere
Vector< S, 3 > const & getCenter() const
IndicatorSphere3D(Vector< S, 3 > center, S radius)
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
S const getRadius() const
bool distance(S &distance, const Vector< S, 3 > &origin, const Vector< S, 3 > &direction, int iC=-1) override
Vector< S, 3 > const & getCenter() const
IndicatorSuperEllipsoid3D(Vector< S, 3 > center, S xHalfAxis, S yHalfAxis, S zHalfAxis, S exponent1, S exponent2)
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
Vector< S, 3 > & getMin() override
bool operator()(bool output[], const S input[]) override
has to be implemented for 'every' derived class
S signedDistance(const Vector< S, 3 > &input) override
Returns signed distance to the nearest point on the indicator surface.
IndicatorTranslate3D(std::array< S, 3 > translate, IndicatorF3D< S > &indicator)
Vector< S, 3 > & getMax() override
class for marking output with some text
std::string getAttribute(const std::string &aName) const
std::vector< XMLreader * >::const_iterator begin() const
Returns an iterator.begin() of the child XMLreader This means an iterator to the next level on an XML...
std::string getName() const
return the name of the element
void setWarningsOn(bool warnings) const
switch warnings on/off
std::vector< XMLreader * >::const_iterator end() const
Returns an iterator.end() of the child XMLreader This means an iterator to the next level on an XML t...
This file contains indicator functions.
T rounding(T a, T r) any_platform
Computes a layer of a constant thickness around the surface.
T sphere(Vector< T, D > p, T r) any_platform
Exact signed distance to the surface of circle (in 2D) or sphere (in 3D).
Vector< T, D > translate(Vector< T, D > p, Vector< T, D > origin) any_platform
Translation: The output of this function is used for the calculation of the signed distance to transl...
T cone(Vector< T, 3 > p, Vector< T, 3 > a, Vector< T, 3 > ba, T baba, T ra, T rb) any_platform
Exact signed distance to the surface of three-dimensional (capped) cone.
T box(Vector< T, 2 > p, Vector< T, 2 > b) any_platform
Exact signed distance to the surface of two-dimensional cuboid.
T ellipsoid(Vector< T, 3 > p, Vector< T, 3 > r) any_platform
Calculate the NOT EXACT (!) signed distance to the surface of three-dimensional ellipsoid.
T cylinder(Vector< T, 3 > p, Vector< T, 3 > a, Vector< T, 3 > ba, T baba, T r) any_platform
Exact signed distance to the surface of three-dimensional cylinder.
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > abs(const ADf< T, DIM > &a)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
bool approxEqual(T a, U b, W epsilon)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
ADf< T, DIM > sin(const ADf< T, DIM > &a)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.
std::shared_ptr< IndicatorF3D< S > > createIndicatorF3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorIntersection3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorCylinder3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorUnion3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorSphere3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorCone3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorCircle3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorCuboid3D(XMLreader const ¶ms, bool verbose=false)
std::shared_ptr< IndicatorF3D< S > > createIndicatorWithout3D(XMLreader const ¶ms, bool verbose=false)
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition of a analytical 2D plane embedded in 3D space.
Hyperplane3D & normalTo(const Vector< T, 3 > &normal)
Calculate the spanning vectors of the hyperplane to be orthogonal to the given normal.
Hyperplane3D & originAt(const Vector< T, 3 > &origin)
Center the hyperplane at the given origin vector.