24#ifndef SMOOTH_INDICATOR_F_3D_HH
25#define SMOOTH_INDICATOR_F_3D_HH
40#define M_PI 3.14159265358979323846
43#define M_PI2 1.57079632679489661923
49template <
typename T,
typename S,
bool PARTICLE>
55template <
typename T,
typename S,
bool PARTICLE>
58 :_ind(xLength, yLength, zLength, center)
60 this->_epsilon = epsilon;
61 if constexpr (!PARTICLE) {
62 this->_pos = _ind.getCenter();
66 this->_circumRadius = .5*(
util::sqrt(xLength*xLength+yLength*yLength+zLength*zLength))+0.5*epsilon;
67 if constexpr (!PARTICLE) {
69 this->_pos[0] - this->getCircumRadius(),
70 this->_pos[1] - this->getCircumRadius(),
71 this->_pos[2] - this->getCircumRadius()
75 this->_pos[0] + this->getCircumRadius(),
76 this->_pos[1] + this->getCircumRadius(),
77 this->_pos[2] + this->getCircumRadius()
83template <
typename T,
typename S,
bool PARTICLE>
88template <
typename T,
typename S,
bool PARTICLE>
91 return _ind.
getxLength()*_ind.getyLength()*_ind.getzLength();
94template <
typename T,
typename S,
bool PARTICLE>
97 T
const xLength = _ind.getxLength();
98 T
const yLength = _ind.getyLength();
99 T
const zLength = _ind.getzLength();
100 T
const mass = getVolume()*density;
101 T
const xLength2 = xLength*xLength;
102 T
const yLength2 = yLength*yLength;
103 T
const zLength2 = zLength*zLength;
105 mofi[0] = 1./12.*mass*(yLength2+zLength2);
106 mofi[1] = 1./12.*mass*(xLength2+zLength2);
107 mofi[2] = 1./12.*mass*(yLength2+xLength2);
108 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
111template <
typename T,
typename S,
bool PARTICLE>
114 return _ind.surfaceNormal(pos, meshSize);
117template <
typename T,
typename S,
bool PARTICLE>
121 if constexpr(!PARTICLE) {
123 p = util::executeRotation<S,3,true>(input, this->_rotMat, this->getPos());
128 return _ind.signedDistance(p + _ind.getCenter());
133template <
typename T,
typename S,
bool PARTICLE>
139template <
typename T,
typename S,
bool PARTICLE>
142 :_ind(center, radius)
144 this->_epsilon = epsilon;
145 if constexpr (!PARTICLE) {
151 this->_circumRadius = max_axis+0.5*epsilon;
152 if constexpr (!PARTICLE) {
154 this->_pos[0] - this->getCircumRadius(),
155 this->_pos[1] - this->getCircumRadius(),
156 this->_pos[2] - this->getCircumRadius()
159 this->_pos[0] + this->getCircumRadius(),
160 this->_pos[1] + this->getCircumRadius(),
161 this->_pos[2] + this->getCircumRadius()
167template <
typename T,
typename S,
bool PARTICLE>
172template <
typename T,
typename S,
bool PARTICLE>
176 return 4./3.*
M_PI*radius[0]*radius[1]*radius[2];
179template <
typename T,
typename S,
bool PARTICLE>
183 T
const mass = getVolume()*density;
184 T
const xHalfAxis2 = radius[0]*radius[0];
185 T
const yHalfAxis2 = radius[1]*radius[1];
186 T
const zHalfAxis2 = radius[2]*radius[2];
188 mofi[0] = 0.2*mass*(yHalfAxis2+zHalfAxis2);
189 mofi[1] = 0.2*mass*(xHalfAxis2+zHalfAxis2);
190 mofi[2] = 0.2*mass*(yHalfAxis2+xHalfAxis2);
191 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
194template <
typename T,
typename S,
bool PARTICLE>
197 return _ind.surfaceNormal(pos, meshSize);
200template <
typename T,
typename S,
bool PARTICLE>
205 if constexpr(!PARTICLE) {
207 p = util::executeRotation<S,3,true>(input, this->_rotMat, this->getPos());
212 return _ind.signedDistance(p + _ind.getCenter());
216template <
typename T,
typename S,
bool PARTICLE>
220 Vector<T,9> rotMatrix = {1, 0, 0, 0, 1, 0, 0, 0, 1};
221 if constexpr (!PARTICLE) {
222 pos = this->getPos();
223 rotMatrix = this->getRotationMatrix();
227 T xDist = input[0] - pos[0];
228 T yDist = input[1] - pos[1];
229 T zDist = input[2] - pos[2];
233 T x = pos[0] + rotMatrix[0]*xDist + rotMatrix[3]*yDist + rotMatrix[6]*zDist;
234 T y = pos[1] + rotMatrix[1]*xDist + rotMatrix[4]*yDist + rotMatrix[7]*zDist;
235 T z = pos[2] + rotMatrix[2]*xDist + rotMatrix[5]*yDist + rotMatrix[8]*zDist;
237 T a = (x - pos[0]) / (_ind.getRadius()[0] - 0.5*this->getEpsilon() );
238 T b = (y - pos[1]) / (_ind.getRadius()[1] - 0.5*this->getEpsilon() );
239 T c = (z - pos[2]) / (_ind.getRadius()[2] - 0.5*this->getEpsilon() );
240 T aEps = (x - pos[0]) / (_ind.getRadius()[0] + 0.5*this->getEpsilon() );
241 T bEps = (y - pos[1]) / (_ind.getRadius()[1] + 0.5*this->getEpsilon() );
242 T cEps = (z - pos[2]) / (_ind.getRadius()[2] + 0.5*this->getEpsilon() );
244 if ( (a*a+b*b+c*c) <= 1. ) {
248 if ( (aEps*aEps+bEps*bEps+cEps*cEps) <= 1. ) {
260template <
typename T,
typename S,
bool PARTICLE>
263 :
SmoothIndicatorSuperEllipsoid3D(ind.getCenter(), ind.getXHalfAxis(), ind.getYHalfAxis(), ind.getZHalfAxis(), ind.getExponent1(), ind.getExponent2(), epsilon, theta)
266template <
typename T,
typename S,
bool PARTICLE>
268 S xHalfAxis, S yHalfAxis, S zHalfAxis, S exponent1, S exponent2,
270 :_ind(center, xHalfAxis, yHalfAxis, zHalfAxis, exponent1, exponent2)
272 this->_epsilon = epsilon;
273 if constexpr (!PARTICLE) {
274 this->_pos = _ind.getCenter();
279 this->_circumRadius =
util::sqrt(2.)*max_axis+0.5*epsilon;
280 if constexpr (!PARTICLE) {
282 this->_pos[0] - this->getCircumRadius(),
283 this->_pos[1] - this->getCircumRadius(),
284 this->_pos[2] - this->getCircumRadius()
287 this->_pos[0] + this->getCircumRadius(),
288 this->_pos[1] + this->getCircumRadius(),
289 this->_pos[2] + this->getCircumRadius()
295template <
typename T,
typename S,
bool PARTICLE>
300template <
typename T,
typename S,
bool PARTICLE>
303 return moments(0., 0., 0.);
306template <
typename T,
typename S,
bool PARTICLE>
309 T
const mass = getVolume() * density;
311 mofi[0] = ( moments(0., 2., 0.) + moments(0., 0., 2.) ) * density;
312 mofi[1] = ( moments(2., 0., 0.) + moments(0., 0., 2.) ) * density;
313 mofi[2] = ( moments(0., 2., 0.) + moments(2., 0., 0.) ) * density;
314 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
317template <
typename T,
typename S,
bool PARTICLE>
320 return (std::tgamma(arg1)*std::tgamma(arg2)) / std::tgamma(arg1+arg2);
323template <
typename T,
typename S,
bool PARTICLE>
326 S ex1 = 2./_ind.getExponent1();
327 S ex2 = 2./_ind.getExponent2();
328 S tmp1 = 2./(p+q+2.);
329 S tmp2 =
util::pow(_ind.getXHalfAxis(), p+1.) *
util::pow(_ind.getYHalfAxis(), q+1.) *
util::pow(_ind.getZHalfAxis(), r+1.) * ex1 * ex2;
330 S tmp3 = beta( (r+1.)*(ex1/2.), (p+q+2.)*(ex2/2.)+1. ) * beta( (q+1.)*(ex2/2.), (p+1.)*(ex2/2.) );
331 return tmp1 * tmp2 * tmp3;
334template <
typename T,
typename S,
bool PARTICLE>
337 return _ind.surfaceNormal(pos, meshSize);
340template <
typename T,
typename S,
bool PARTICLE>
345 if constexpr(!PARTICLE) {
347 p = util::executeRotation<S,3,true>(input, this->_rotMat, this->getPos());
352 return _ind.signedDistance(p + _ind.getCenter());
355template <
typename T,
typename S,
bool PARTICLE>
359 Vector<T,9> rotMatrix = {1, 0, 0, 0, 1, 0, 0, 0, 1};
360 if constexpr (!PARTICLE) {
361 pos = this->getPos();
362 rotMatrix = this->getRotationMatrix();
366 T xDist = input[0] - pos[0];
367 T yDist = input[1] - pos[1];
368 T zDist = input[2] - pos[2];
372 T x = pos[0] + rotMatrix[0]*xDist + rotMatrix[3]*yDist + rotMatrix[6]*zDist;
373 T y = pos[1] + rotMatrix[1]*xDist + rotMatrix[4]*yDist + rotMatrix[7]*zDist;
374 T z = pos[2] + rotMatrix[2]*xDist + rotMatrix[5]*yDist + rotMatrix[8]*zDist;
376 T a =
util::pow (
util::abs( (x-pos[0]) / _ind.getXHalfAxis() ), _ind.getExponent1() );
377 T b =
util::pow (
util::abs( (y-pos[1]) / _ind.getYHalfAxis() ), _ind.getExponent1() );
378 T c =
util::pow (
util::abs( (z-pos[2]) / _ind.getZHalfAxis() ), _ind.getExponent2() );
379 T ab =
util::pow( a+b, _ind.getExponent2()/_ind.getExponent1() );
381 if ( (ab+c) <= 1. ) {
392template <
typename T,
typename S,
bool PARTICLE>
397template <
typename T,
typename S,
bool PARTICLE>
399 : _ind(center, radius)
401 this->_epsilon = epsilon;
402 if constexpr (!PARTICLE) {
406 this->_circumRadius = radius + 0.5*epsilon;
407 if constexpr (!PARTICLE) {
409 this->_pos[0] - this->getCircumRadius(),
410 this->_pos[1] - this->getCircumRadius(),
411 this->_pos[2] - this->getCircumRadius()
414 this->_pos[0] + this->getCircumRadius(),
415 this->_pos[1] + this->getCircumRadius(),
416 this->_pos[2] + this->getCircumRadius()
423template <
typename T,
typename S,
bool PARTICLE>
428template <
typename T,
typename S,
bool PARTICLE>
434template <
typename T,
typename S,
bool PARTICLE>
437 T
const radius = _ind.getRadius();
438 T
const radius2 = radius*radius;
439 T
const mass = getVolume()*density;
441 mofi[0] = 2./5.*mass*radius2;
442 mofi[1] = 2./5.*mass*radius2;
443 mofi[2] = 2./5.*mass*radius2;
444 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
447template <
typename T,
typename S,
bool PARTICLE>
450 return _ind.surfaceNormal(pos, meshSize);
453template <
typename T,
typename S,
bool PARTICLE>
458 if constexpr (!PARTICLE) {
461 return _ind.signedDistance(dist);
465template <
typename T,
typename S,
bool PARTICLE>
468 this->_circumRadius =
util::sqrt(_ind.getRadius()*_ind.getRadius()+(0.5*length)*(0.5*length))+0.5*this->getEpsilon();
469 if constexpr (!PARTICLE) {
471 this->_pos[0] - this->getCircumRadius(),
472 this->_pos[1] - this->getCircumRadius(),
473 this->_pos[2] - this->getCircumRadius()
476 this->_pos[0] + this->getCircumRadius(),
477 this->_pos[1] + this->getCircumRadius(),
478 this->_pos[2] + this->getCircumRadius()
484template <
typename T,
typename S,
bool PARTICLE>
489template <
typename T,
typename S,
bool PARTICLE>
494template <
typename T,
typename S,
bool PARTICLE>
496 : _ind(center1, center2, radius)
498 this->_epsilon = epsilon;
499 if constexpr (!PARTICLE) {
500 this->_pos = 0.5 * (_ind.getCenter1() + _ind.getCenter2());
503 Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
504 S
const length =
norm(dist);
505 if constexpr (!PARTICLE) {
506 initIndicatorCylinder3D(this->_theta, length);
514template <
typename T,
typename S,
bool PARTICLE>
515SmoothIndicatorCylinder3D<T,S,PARTICLE>::SmoothIndicatorCylinder3D(
Vector<S,3> center,
Vector<S,3> normal, S radius, S height, S epsilon,
Vector<S,3> theta)
519template <
typename T,
typename S,
bool PARTICLE>
522 const Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
523 S
const length =
norm(dist);
524 return M_PI*_ind.getRadius()*_ind.getRadius()*length;
527template <
typename T,
typename S,
bool PARTICLE>
530 Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
531 S
const length =
norm(dist);
532 T
const radius2 = _ind.getRadius() * _ind.getRadius();
533 T
const mass = getVolume()*density;
535 mofi[0] = 0.5*mass*radius2;
536 mofi[1] = 1/12.*mass*(length*length+3.*radius2);
537 mofi[2] = 1/12.*mass*(length*length+3.*radius2);
538 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
541template <
typename T,
typename S,
bool PARTICLE>
544 return _ind.surfaceNormal(pos, meshSize);
547template <
typename T,
typename S,
bool PARTICLE>
551 if constexpr(!PARTICLE) {
553 p = util::executeRotation<S,3,true>(input, this->_rotMat, this->getPos());
560 return _ind.signedDistance(p + 0.5 * (_ind.getCenter2() + _ind.getCenter1()));
564template <
typename T,
typename S,
bool PARTICLE>
566 :
SmoothIndicatorCone3D(ind.getCenter1(), ind.getCenter2(), ind.getRadius1(), ind.getRadius2(), epsilon, theta)
569template <
typename T,
typename S,
bool PARTICLE>
570SmoothIndicatorCone3D<T,S,PARTICLE>::SmoothIndicatorCone3D(
Vector<S,3> center1,
Vector<S,3> center2, S radius1, S radius2, S epsilon,
Vector<S,3> theta)
571 :_ind(center1, center2, radius1, radius2)
573 this->_epsilon = epsilon;
576 if constexpr (!PARTICLE) {
577 this->_pos = _startPos;
581 Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
582 S
const length =
norm(dist);
583 initIndicatorCone3D(theta, length);
586template <
typename T,
typename S,
bool PARTICLE>
592template <
typename T,
typename S,
bool PARTICLE>
595 const T _radiusA = _ind.getRadius1();
596 const T _radiusB = _ind.getRadius2();
598 if (_radiusA >= _radiusB) {
599 this->_circumRadius =
util::sqrt(_radiusA*_radiusA+(0.5*length)*(0.5*length))+0.5*this->getEpsilon();
602 this->_circumRadius =
util::sqrt(_radiusB*_radiusB+(0.5*length)*(0.5*length))+0.5*this->getEpsilon();
605 if constexpr (!PARTICLE) {
607 this->_pos[0] - this->getCircumRadius(),
608 this->_pos[1] - this->getCircumRadius(),
609 this->_pos[2] - this->getCircumRadius()
612 this->_pos[0] + this->getCircumRadius(),
613 this->_pos[1] + this->getCircumRadius(),
614 this->_pos[2] + this->getCircumRadius()
621template <
typename T,
typename S,
bool PARTICLE>
624 const Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
625 S
const length =
norm(dist);
626 return (1/3.)*
M_PI*length* (_ind.getRadius1()*_ind.getRadius1()
627 +_ind.getRadius1()*_ind.getRadius2()
628 +_ind.getRadius2()*_ind.getRadius2());
631template <
typename T,
typename S,
bool PARTICLE>
634 const T _radiusA = _ind.getRadius1();
635 const T _radiusB = _ind.getRadius2();
637 Vector<T,3> dist = _ind.getCenter2() - _ind.getCenter1();
638 S
const length =
norm(dist);
640 const T mass = getVolume()*density;
646 T topRadius = _radiusA;
647 T baseRadius = _radiusB;
648 if (_radiusA > _radiusB) {
649 topRadius = _radiusB;
650 baseRadius = _radiusA;
653 mofi[1] = 3/20.*mass*( (
util::pow(baseRadius-topRadius, 2) + 4*length*length) * (
util::pow(baseRadius,5) -
util::pow(topRadius,5)) )
654 / (
util::pow(baseRadius-topRadius,3) * (baseRadius*baseRadius+baseRadius*topRadius+topRadius*topRadius));
656 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
659template <
typename T,
typename S,
bool PARTICLE>
662 T topRadius, baseRadius;
665 if (_ind.getRadius1() < _ind.getRadius2()) {
666 topRadius = _ind.getRadius1();
667 baseRadius = _ind.getRadius2();
668 centerTop = _ind.getCenter1();
669 centerBase = _ind.getCenter2();
671 else if (_ind.getRadius1() > _ind.getRadius2()) {
672 topRadius = _ind.getRadius2();
673 baseRadius = _ind.getRadius1();
674 centerTop = _ind.getCenter2();
675 centerBase = _ind.getCenter1();
678 std::cerr <<
"Error calculating a cone's center of mass." << std::endl;
682 const T height =
norm(centerTop - centerBase);
683 const T centerOfMassHeight = 0.25*height * (baseRadius*baseRadius + 2*baseRadius*topRadius + 3*topRadius*topRadius)
684 / (baseRadius*baseRadius + baseRadius*topRadius + topRadius*topRadius);
685 return centerBase +
normalize(centerTop-centerBase) * centerOfMassHeight;
688template <
typename T,
typename S,
bool PARTICLE>
691 return _ind.surfaceNormal(pos, meshSize);
694template <
typename T,
typename S,
bool PARTICLE>
699 if constexpr(!PARTICLE) {
701 p = util::executeRotation<S,3,true>(input, this->_rotMat, this->getPos());
708 return _ind.signedDistance(p + _startPos);
714template <
typename T,
typename S,
bool PARTICLE>
721 _latticeSpacing(latticeSpacing),
725 this->_name =
"custom3D";
726 this->_epsilon = epsilon;
727 if constexpr (!PARTICLE) {
733 initBlockData(*indPtr);
741template <
typename T,
typename S,
bool PARTICLE>
749 for (
unsigned iD=0; iD<3; ++iD) {
753 blockDataPadding[iD] = 2*
util::ceil(0.2*blockDataSize[iD]);
754 blockDataSize[iD] += blockDataPadding[iD];
758 _cuboid = Cuboid3D<T>(PhysR<T,3>(0.), _latticeSpacing, blockDataSize);
759 this->_blockData.reset(
new BlockData<3,T,BaseType<T>>(_cuboid));
761 for (iX[0]=0; iX[0] < this->_blockData->getNx(); ++iX[0]) {
762 for (iX[1]=0; iX[1] < this->_blockData->getNy(); ++iX[1]) {
763 for (iX[2]=0; iX[2] < this->_blockData->getNz(); ++iX[2]) {
765 for (
unsigned iD=0; iD<3; ++iD) {
766 input[iD] = (iX[iD]-blockDataPadding[iD]/2)*_latticeSpacing+ind.
getMin()[iD];
773 this->_cacheFunctor = std::make_unique<BlockDataF3D<T,BaseType<T>>>(*(this->_blockData));
774 this->_interpolateCache = std::make_unique<AnalyticalFfromBlockF3D<T,T>>(*(this->_cacheFunctor), _cuboid);
777template <
typename T,
typename S,
bool PARTICLE>
778void SmoothIndicatorCustom3D<T,S,PARTICLE>::calcCenter()
783 this->_center = PhysR<T,3>(0.);
784 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
785 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
786 for (input[2] = 0; input[2] < this->_blockData->getNz(); ++input[2]) {
787 if (regardCell(input)) {
789 this->_center[0] += this->_latticeSpacing*input[0];
790 this->_center[1] += this->_latticeSpacing*input[1];
791 this->_center[2] += this->_latticeSpacing*input[2];
797 this->_center *= 1./nCells;
800template <
typename T,
typename S,
bool PARTICLE>
806template <
typename T,
typename S,
bool PARTICLE>
810 T cuboidMofi =
util::pow(this->_latticeSpacing, 2)/ 6.0;
815 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
816 dx =
util::abs(this->_latticeSpacing*input[0] - this->_center[0]);
817 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
818 dy =
util::abs(this->_latticeSpacing*input[1] - this->_center[1]);
819 for (input[2] = 0; input[2] < this->_blockData->getNz(); ++input[2]) {
820 if (regardCell(input)) {
821 dz =
util::abs(this->_latticeSpacing*input[2] - this->_center[2]);
822 mofi[0] += (dy*dy+dz*dz+cuboidMofi);
823 mofi[1] += (dx*dx+dz*dz+cuboidMofi);
824 mofi[2] += (dx*dx+dy*dy+cuboidMofi);
830 _volume = nCells *
util::pow(_latticeSpacing,3);
831 const T mass = rhoP * _volume;
832 const T cuboidMass = mass/nCells;
835 return Vector<T,4>(mofi[0], mofi[1], mofi[2], mass);
839template <
typename T,
typename S,
bool PARTICLE>
848 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
849 distance[0] = this->_latticeSpacing * input[0] - this->_center[0];
850 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
851 distance[1] = this->_latticeSpacing * input[1] - this->_center[1];
852 for (input[2] = 0; input[2] < this->_blockData->getNz(); ++input[2]) {
853 distance[2] = this->_latticeSpacing * input[2] - this->_center[2];
854 if (regardCell(input)) {
855 if constexpr (!PARTICLE) {
856 for (
unsigned iD=0; iD<3; ++iD) {
857 min[iD] =
util::min(distance[iD], min[iD]);
858 max[iD] =
util::max(distance[iD], max[iD]);
867 if constexpr (!PARTICLE) {
868 min -= Vector<T,3>(0.5 * this->_epsilon + _latticeSpacing);
869 max += Vector<T,3>(0.5 * this->_epsilon + _latticeSpacing);
870 this->_myMin = this->_pos +
min;
871 this->_myMax = this->_pos +
max;
874 this->_circumRadius = maxDistance + T{0.5} * (this->_epsilon +
util::sqrt(3) * _latticeSpacing);
877template <
typename T,
typename S,
bool PARTICLE>
880 return this->_center;
883template <
typename T,
typename S,
bool PARTICLE>
886 return _indPtr->surfaceNormal(pos, meshSize);
889template <
typename T,
typename S,
bool PARTICLE>
893 return _indPtr->surfaceNormal( pos, meshSize, transformPos );
896template <
typename T,
typename S,
bool PARTICLE>
900 if constexpr (!PARTICLE) {
902 position = util::executeRotation<T,3,true>(input, this->_rotMat, this->getPos());
905 const PhysR<T,3> positionInCache = this->_center + position;
907 T signedDistance(0.);
908 if(_interpolateCache->operator()(&signedDistance, positionInCache.
data())) {
909 return signedDistance;
915 for(
unsigned iDim=0; iDim<3; ++iDim) {
916 latticePosition[iDim] =
util::round( positionInCache[iDim] / this->_latticeSpacing );
917 latticePosition[iDim] =
util::max(0, latticePosition[iDim]);
918 latticePosition[iDim] =
util::min(this->_blockData->getExtent()[iDim] - 1, latticePosition[iDim]);
920 extraDistance[iDim] =
util::abs(_latticeSpacing * latticePosition[iDim] - positionInCache[iDim]);
922 return this->_blockData->get(latticePosition.
data()) +
norm(extraDistance);
925template <
typename T,
typename S,
bool PARTICLE>
928 return this->_blockData->get(input) < std::numeric_limits<T>::epsilon();
931template <
typename T,
typename S,
bool PARTICLE>
935 if constexpr(!PARTICLE) {
936 if(
norm(this->getPos() - pos) > this->_circumRadius) {
940 pos = util::executeRotation<S,3,true>(pos, this->_rotMat, this->getPos());
942 if(
norm(pos) < this->_circumRadius) {
950template <
typename T,
typename S,
bool PARTICLE>
951SmoothIndicatorFactoredCircle3D<T,S,PARTICLE>::SmoothIndicatorFactoredCircle3D(
Vector<S,3> center, S radius, S epsilon, S density,
Vector<S,3> vel, S omega, S factor)
952 : _radius(radius), _factor(factor)
954 this->_epsilon = epsilon;
955 if constexpr (!PARTICLE) {
959 this->_circumRadius = radius + 0.5*epsilon;
960 if constexpr (!PARTICLE) {
961 this->_myMin = {center[0] - this->getCircumRadius(), center[1] - this->getCircumRadius(), center[2] - this->getCircumRadius()};
962 this->_myMax = {center[0] + this->getCircumRadius(), center[1] + this->getCircumRadius(), center[2] + this->getCircumRadius()};
967template <
typename T,
typename S,
bool PARTICLE>
970 T mass = 4/3*
M_PI*_radius*_radius*_radius*density;
972 mofi[0] = 2./5.*mass*_radius*_radius;
973 mofi[1] = 2./5.*mass*_radius*_radius;
974 mofi[2] = 2./5.*mass*_radius*_radius;
975 return Vector<S,4>(mofi[0], mofi[1], mofi[2], mass);
979template <
typename T,
typename S,
bool PARTICLE>
982 double distToCenter2 =
util::pow(this->getPos()[0]-input[0], 2) +
983 util::pow(this->getPos()[1]-input[1], 2) +
984 util::pow(this->getPos()[2]-input[2], 2);
1000 double d =
util::sqrt(distToCenter2) - this->_radius;
1001 output[0] = T(this->_factor*(1.-tanh(d/this->getEpsilon()))/2.);
indicator function for a 3d frustum
indicator function for a 3d-cuboid, parallel to the planes x=0, y=0, z=0.
S const getxLength() const
indicator function for a 3d-cylinder
indicator function for an ellipsoid
IndicatorF3D is an application from .
virtual Vector< S, 3 > & getMax()
virtual Vector< S, 3 > & getMin()
virtual S signedDistance(const Vector< S, 3 > &input)
Returns signed distance to the nearest point on the indicator surface.
indicator function for a 3D-sphere
indicator function for a super ellipsoid
class for marking output with some text
implements a smooth particle cone in 3D with an _epsilon sector
const S signedDistance(const PhysR< S, 3 > input) override
SmoothIndicatorCone3D(IndicatorCone3D< S > &indPtr, S epsilon, Vector< S, 3 > theta=Vector< S, 3 >(0., 0., 0.))
IndicatorCone3D< S > & getIndicator()
Vector< S, 4 > calcMofiAndMass(const S density) override
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
Vector< S, 3 > calcCenterOfMass() override
implements a smooth particle cuboid in 3D with an _epsilon sector.
Vector< S, 4 > calcMofiAndMass(const S density) override
SmoothIndicatorCuboid3D(IndicatorCuboid3D< S > &ind, S epsilon, Vector< S, 3 > theta=Vector< S, 3 >(0., 0., 0.))
const S signedDistance(const PhysR< S, 3 > input) override
IndicatorCuboid3D< S > & getIndicator()
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
SmoothIndicatorCustom3D(T latticeSpacing, std::shared_ptr< IndicatorF3D< T > > indPtr, Vector< T, 3 > pos, T epsilon, Vector< T, 3 > theta=Vector< T, 3 >(0.))
const S signedDistance(const PhysR< S, 3 > input) override
bool operator()(T output[], const S input[]) override
Vector< T, 3 > getLocalCenter()
Vector< T, 4 > calcMofiAndMass(T rhoP) override
bool regardCell(int input[3])
implements a smooth particle cylinder in 3D with an _epsilon sector.
const S signedDistance(const PhysR< S, 3 > input) override
SmoothIndicatorCylinder3D(IndicatorCylinder3D< S > &ind, S epsilon, Vector< S, 3 > theta=Vector< S, 3 >(0., 0., 0.))
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
Vector< S, 4 > calcMofiAndMass(const S density) override
IndicatorCylinder3D< S > & getIndicator()
implements a smooth particle ellipsoid in 3D with an _epsilon sector.
Vector< S, 4 > calcMofiAndMass(const S density) override
const S signedDistance(const PhysR< S, 3 > input) override
IndicatorEllipsoid3D< S > & getIndicator()
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
bool operator()(T output[], const S input[]) override
SmoothIndicatorEllipsoid3D(IndicatorEllipsoid3D< S > &ind, S epsilon, Vector< S, 3 > theta=Vector< S, 3 >(0., 0., 0.))
SmoothIndicatorFactoredCircle3D(Vector< S, 3 > center, S radius, S epsilon, S density=0, Vector< S, 3 > vel=Vector< S, 3 >(0., 0., 0.), S omega=0, S factor=1.)
bool operator()(T output[], const S input[]) override
Vector< S, 4 > calcMofiAndMass(const S density) override
implements a smooth sphere in 3D with an _epsilon sector
SmoothIndicatorSphere3D(IndicatorSphere3D< S > &ind, S epsilon)
const S signedDistance(const PhysR< S, 3 > input) override
IndicatorSphere3D< S > & getIndicator()
Vector< S, 4 > calcMofiAndMass(const S density) override
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
implements a smooth particle super-ellipsoid in 3D. The epsilon sector is currently missing.
SmoothIndicatorSuperEllipsoid3D(IndicatorSuperEllipsoid3D< S > &ind, S epsilon, Vector< S, 3 > theta=Vector< S, 3 >(0., 0., 0.))
const S signedDistance(const PhysR< S, 3 > input) override
IndicatorSuperEllipsoid3D< S > & getIndicator()
Vector< S, 3 > surfaceNormal(const Vector< S, 3 > &pos, const S meshSize) override
bool operator()(T output[], const S input[]) override
Vector< S, 4 > calcMofiAndMass(const S density) override
constexpr const T * data() const any_platform
Pack< T > min(Pack< T > rhs, Pack< T > lhs)
Pack< T > max(Pack< T > rhs, Pack< T > lhs)
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > abs(const ADf< T, DIM > &a)
ADf< T, DIM > ceil(const ADf< T, DIM > &a)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
ADf< T, DIM > round(const ADf< T, DIM > &a)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
decltype(Vector< decltype(util::sqrt(T())), D >()) degreeToRadian(const Vector< T, D > &angle)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Top level namespace for all of OpenLB.
typename util::BaseTypeHelper< T >::type BaseType
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})