24#ifndef SMOOTH_INDICATOR_F_2D_HH
25#define SMOOTH_INDICATOR_F_2D_HH
34#define M_PI 3.14159265358979323846
36#define M_PI2 1.57079632679489661923
40template <
typename T,
typename S,
bool PARTICLE>
45template <
typename T,
typename S,
bool PARTICLE>
47 : _ind(xLength, yLength, center)
49 this->_epsilon = epsilon;
50 if constexpr (!PARTICLE) {
51 this->_pos = _ind.getCenter();
56 if constexpr (!PARTICLE) {
58 this->_pos[0] - this->getCircumRadius(),
59 this->_pos[1] - this->getCircumRadius()
62 this->_pos[0] + this->getCircumRadius(),
63 this->_pos[1] + this->getCircumRadius()
69template <
typename T,
typename S,
bool PARTICLE>
74template <
typename T,
typename S,
bool PARTICLE>
80template <
typename T,
typename S,
bool PARTICLE>
83 const T xLength = _ind.getxLength();
84 const T yLength = _ind.getyLength();
85 const T mass = getArea()*density;
86 const T mofi = 1./12.*mass*(xLength*xLength+yLength*yLength);
90template <
typename T,
typename S,
bool PARTICLE>
93 return _ind.surfaceNormal(pos, meshSize);
96template <
typename T,
typename S,
bool PARTICLE>
100 if constexpr(!PARTICLE) {
102 p = util::executeRotation<S,2,true>(input, this->_rotMat, this->getPos());
107 return _ind.signedDistance(p + _ind.getCenter());
110template <
typename T,
typename S,
bool PARTICLE>
114 if constexpr(!PARTICLE) {
116 p = util::executeRotation<S,2,true>(origin, this->_rotMat, this->getPos());
121 return _ind.distance(distance, p, direction, precision, pitch);
125template <
typename T,
typename S,
bool PARTICLE>
130template <
typename T,
typename S,
bool PARTICLE>
132 : _ind(center, radius)
134 this->_epsilon = epsilon;
135 if constexpr (!PARTICLE) {
136 this->_pos = _ind.getCenter();
139 this->_circumRadius = radius + 0.5*epsilon;
140 if constexpr (!PARTICLE) {
141 this->_myMin = {this->_pos[0] - this->getCircumRadius(), this->_pos[1] - this->getCircumRadius()};
142 this->_myMax = {this->_pos[0] + this->getCircumRadius(), this->_pos[1] + this->getCircumRadius()};
148template <
typename T,
typename S,
bool PARTICLE>
153template <
typename T,
typename S,
bool PARTICLE>
159template <
typename T,
typename S,
bool PARTICLE>
162 const T radius = _ind.getRadius();
163 const T mass = getArea()*density;
164 const T mofi = 0.5 * mass * radius * radius;
168template <
typename T,
typename S,
bool PARTICLE>
171 return _ind.surfaceNormal(pos, meshSize);
174template <
typename T,
typename S,
bool PARTICLE>
178 if constexpr(!PARTICLE) {
181 return _ind.signedDistance(dist + _ind.getCenter());
184template <
typename T,
typename S,
bool PARTICLE>
188 if constexpr(!PARTICLE) {
191 return _ind.distance(distance, dist, direction, precision, pitch);
195template <
typename T,
typename S,
bool PARTICLE>
200template <
typename T,
typename S,
bool PARTICLE>
202 : _ind(center, radius)
204 this->_epsilon = epsilon;
205 if constexpr (!PARTICLE) {
206 this->_pos = _ind.getCenter();
210 this->_circumRadius = radius + 0.5*epsilon;
211 if constexpr (!PARTICLE) {
212 this->_myMin = {center[0] - this->getCircumRadius(), center[1] - this->getCircumRadius()};
213 this->_myMax = {center[0] + this->getCircumRadius(), center[1] + this->getCircumRadius()};
218template <
typename T,
typename S,
bool PARTICLE>
223template <
typename T,
typename S,
bool PARTICLE>
227 const T altitude = 1.5*radius;
229 return 0.5*base*altitude;
232template <
typename T,
typename S,
bool PARTICLE>
235 const T radius = _ind.getRadius();
236 const T altitude = 1.5*radius;
238 const T mass = density*getArea();
239 const T mofi = mass*((altitude*altitude/18.)+(base*base/24.));
243template <
typename T,
typename S,
bool PARTICLE>
246 return _ind.surfaceNormal(pos, meshSize);
249template <
typename T,
typename S,
bool PARTICLE>
253 if constexpr(!PARTICLE) {
255 p = util::executeRotation<S,2,true>(input, this->_rotMat, this->getPos());
260 return _ind.signedDistance(p + _ind.getCenter());
263template <
typename T,
typename S,
bool PARTICLE>
267 if constexpr(!PARTICLE) {
269 p = util::executeRotation<S,2>(origin, this->_rotMat, this->getPos());
276 return _ind.distance(distance, p, direction, precision, pitch);
282template <
typename T,
typename S,
bool PARTICLE>
289 _latticeSpacing(latticeSpacing)
292 this->_name =
"custom2D";
293 this->_epsilon = epsilon;
294 if constexpr(!PARTICLE) {
303template <
typename T,
typename S,
bool PARTICLE>
314template <
typename T,
typename S,
bool PARTICLE>
315void SmoothIndicatorCustom2D<T,S,PARTICLE>::initBlockData(IndicatorF2D<T>& ind)
317 OstreamManager clout(std::cout,
"createIndicatorCustom2D");
320 Vector<int,2> blockDataSize;
321 Vector<int,2 >blockDataPadding;
322 for (
unsigned iD=0; iD<2; ++iD) {
324 util::ceil( (ind.getMax()[iD] - ind.getMin()[iD]) / _latticeSpacing ), 1);
326 blockDataPadding[iD] = 2*
util::ceil(0.2*blockDataSize[iD]);
327 blockDataSize[iD] += blockDataPadding[iD];
331 _cuboid = std::make_unique<Cuboid2D<T>>(PhysR<T,2>(0.), _latticeSpacing, blockDataSize);
332 this->_blockData.reset(
new BlockData<2,T,BaseType<T>>(*_cuboid));
334 for (iX[0]=0; iX[0] < this->_blockData->getNx(); ++iX[0]) {
335 for (iX[1]=0; iX[1] < this->_blockData->getNy(); ++iX[1]) {
337 for (
unsigned iD=0; iD<2; ++iD) {
338 input[iD] = (iX[iD]-blockDataPadding[iD]/2)*_latticeSpacing+ind.getMin()[iD];
340 this->_blockData->get(iX) = ind.signedDistance(input);
344 this->_cacheFunctor = std::make_unique<BlockDataF2D<T,BaseType<T>>>(*(this->_blockData));
345 this->_interpolateCache = std::make_unique<AnalyticalFfromBlockF2D<T,T>>(*(this->_cacheFunctor), *_cuboid);
348template <
typename T,
typename S,
bool PARTICLE>
349void SmoothIndicatorCustom2D<T,S,PARTICLE>::calcCenter()
354 this->_center = Vector<T,2>(0.);
355 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
356 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
357 if (regardCell(input)) {
359 const unsigned short porosity = 1;
360 this->_center[0] += porosity*this->_latticeSpacing*input[0];
361 this->_center[1] += porosity*this->_latticeSpacing*input[1];
366 this->_center *= 1./nCells;
369template <
typename T,
typename S,
bool PARTICLE>
376template <
typename T,
typename S,
bool PARTICLE>
384 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
385 const T dx =
util::abs(input[0]*_latticeSpacing - this->_center[0]);
386 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
387 if (regardCell(input)) {
388 const T dy =
util::abs(input[1]*_latticeSpacing - this->_center[1]);
389 mofi += (dx*dx+dy*dy);
394 _area = nCells *
util::pow(_latticeSpacing, 2);
396 mofi +=
util::pow(_latticeSpacing, 4)/ 6.0;
402template <
typename T,
typename S,
bool PARTICLE>
411 for (input[0] = 0; input[0] < this->_blockData->getNx(); ++input[0]) {
412 distance[0] = this->_latticeSpacing*input[0] - this->_center[0];
413 for (input[1] = 0; input[1] < this->_blockData->getNy(); ++input[1]) {
414 distance[1] = this->_latticeSpacing*input[1] - this->_center[1];
415 if (regardCell(input)) {
416 if constexpr (!PARTICLE) {
417 for (
unsigned iD=0; iD<2; ++iD) {
418 min[iD] =
util::min(distance[iD], min[iD]);
419 max[iD] =
util::max(distance[iD], max[iD]);
427 if constexpr (!PARTICLE) {
428 min -= Vector<T,2>(0.5 * this->_epsilon + _latticeSpacing);
429 max += Vector<T,2>(0.5 * this->_epsilon + _latticeSpacing);
430 this->_myMin = this->_pos +
min;
431 this->_myMax = this->_pos +
max;
434 this->_circumRadius = maxDistance + T{0.5} * (this->_epsilon +
util::sqrt(2) * _latticeSpacing);
437template <
typename T,
typename S,
bool PARTICLE>
440 return this->_center;
443template <
typename T,
typename S,
bool PARTICLE>
446 return _indPtr->surfaceNormal(pos, meshSize);
449template <
typename T,
typename S,
bool PARTICLE>
453 if constexpr (!PARTICLE) {
455 position = util::executeRotation<T,2,true>(input, this->_rotMat, this->getPos());
458 const PhysR<T,2> positionInCache = this->_center + position;
460 T signedDistance(0.);
461 if(_interpolateCache->operator()(&signedDistance, positionInCache.
data())) {
462 return signedDistance;
468 for(
unsigned iDim=0; iDim<2; ++iDim) {
469 latticePosition[iDim] =
util::round( positionInCache[iDim] / this->_latticeSpacing );
470 latticePosition[iDim] =
util::max(0, latticePosition[iDim]);
471 latticePosition[iDim] =
util::min(this->_blockData->getExtent()[iDim] - 1, latticePosition[iDim]);
473 extraDistance[iDim] =
util::abs(_latticeSpacing * latticePosition[iDim] - positionInCache[iDim]);
475 return this->_blockData->get(latticePosition.
data()) +
norm(extraDistance);
478template <
typename T,
typename S,
bool PARTICLE>
481 return this->_blockData->get(input) < std::numeric_limits<T>::epsilon();
484template <
typename T,
typename S,
bool PARTICLE>
488 if constexpr(!PARTICLE) {
489 pos = util::executeRotation<S,2,true>(pos, this->_rotMat, this->getPos());
491 if(
norm(pos) < this->_circumRadius) {
500template <
typename T,
typename S,
bool PARTICLE>
504 this->_epsilon = epsilon;
505 if constexpr (!PARTICLE) {
509 this->_circumRadius = radius + 0.5*epsilon;
510 if constexpr (!PARTICLE) {
511 this->_myMin = {center[0] - this->getCircumRadius(), center[1] - this->getCircumRadius()};
512 this->_myMax = {center[0] + this->getCircumRadius(), center[1] + this->getCircumRadius()};
517 T mass =
M_PI*radius*radius*density;
518 T mofi = 0.5 * mass * radius * radius;
521template <
typename T,
typename S,
bool PARTICLE>
524 T mass =
M_PI*_radius*_radius*density;
525 T mofi = 0.5 * mass * _radius * _radius;
530template <
typename T,
typename S,
bool PARTICLE>
533 double distToCenter2 =
util::pow(this->getPos()[0]-input[0], 2) +
534 util::pow(this->getPos()[1]-input[1], 2);
537 double d =
util::sqrt(distToCenter2) - this->_radius;
538 output[0] = T((1.-tanh(d/this->getEpsilon()))/2.);
543template <
typename T,
typename S,
bool PARTICLE>
544SmoothIndicatorFactoredCircle2D<T,S,PARTICLE>::SmoothIndicatorFactoredCircle2D(
Vector<S,2> center, S radius, S epsilon, S density,
Vector<S,2> vel, S omega, S factor)
545 : _radius(radius), _factor(factor)
547 this->_epsilon = epsilon;
548 if constexpr (!PARTICLE) {
552 this->_circumRadius = radius + 0.5*epsilon;
553 if constexpr (!PARTICLE) {
554 this->_myMin = {center[0] - this->getCircumRadius(), center[1] - this->getCircumRadius()};
555 this->_myMax = {center[0] + this->getCircumRadius(), center[1] + this->getCircumRadius()};
560template <
typename T,
typename S,
bool PARTICLE>
563 T mass =
M_PI*_radius*_radius*density;
564 T mofi = 0.5 * mass * _radius * _radius;
569template <
typename T,
typename S,
bool PARTICLE>
572 double distToCenter2 =
util::pow(this->getPos()[0]-input[0], 2) +
573 util::pow(this->getPos()[1]-input[1], 2);
589 double d =
util::sqrt(distToCenter2) - this->_radius;
590 output[0] = T(this->_factor*(1.-tanh(d/this->getEpsilon()))/2.);
595template <
typename T,
typename S,
bool PARTICLE>
596SmoothIndicatorFactoredCuboid2D<T,S,PARTICLE>::SmoothIndicatorFactoredCuboid2D(
Vector<S,2> center, S xLength, S yLength, S epsilon, S density,
Vector<S,2> vel, S omega, S factor)
597 : _xLength(xLength), _yLength(yLength), _factor(factor)
599 this->_epsilon = epsilon;
600 if constexpr (!PARTICLE) {
605 if constexpr (!PARTICLE) {
607 this->_pos[0] - this->getCircumRadius(),
608 this->_pos[1] - this->getCircumRadius()
611 this->_pos[0] + this->getCircumRadius(),
612 this->_pos[1] + this->getCircumRadius()
618template <
typename T,
typename S,
bool PARTICLE>
621 T mass = _xLength*_yLength*density;
622 T mofi = 1./12.*mass*(_xLength*_xLength+_yLength*_yLength);
627template <
typename T,
typename S,
bool PARTICLE>
632 T distToCenter2 =
util::pow(this->getPos()[1]-input[1], 2);
633 d =
util::sqrt(distToCenter2) - this->_yLength/2. ;
635 else if(_yLength == 0){
636 T distToCenter2 =
util::pow(this->getPos()[0]-input[0], 2);
637 d =
util::sqrt(distToCenter2) - this->_xLength/2. + this->getEpsilon()/2;
640 T distToCenter2x =
util::pow(this->getPos()[0]-input[0], 2);
641 T distToCenter2y =
util::pow(this->getPos()[1]-input[1], 2);
643 if(ratio >= _yLength/_xLength){
644 d =
util::sqrt(distToCenter2y) - this->_yLength/2. + this->getEpsilon()/2;
647 d =
util::sqrt(distToCenter2x) - this->_xLength/2. + this->getEpsilon()/2;
663 output[0] = T(this->_factor*(1.-tanh(d/this->getEpsilon()))/2.);
indicator function for a 2D circle
S const getRadius() const
indicator function for a 2D-cuboid, parallel to the planes x=0, y=0; theta rotates cuboid around its ...
S const getxLength() const
indicator function for a 2D equilateral triangle
S const getRadius() const
IndicatorF2D is an application from .
class for marking output with some text
implements a smooth circle in 2D with an _epsilon sector
Vector< S, 2 > surfaceNormal(const Vector< S, 2 > &pos, const S meshSize) override
IndicatorCircle2D< S > & getIndicator()
Vector< S, 2 > calcMofiAndMass(const S density) override
SmoothIndicatorCircle2D(IndicatorCircle2D< S > &indPtr, S epsilon)
const S signedDistance(const Vector< S, 2 > input) override
bool distance(S &distance, const Vector< S, 2 > &origin, const Vector< S, 2 > &direction, S precision, S pitch) override
implements a smooth cuboid in 2D with an _epsilon sector.
bool distance(S &distance, const Vector< S, 2 > &origin, const Vector< S, 2 > &direction, S precision, S pitch) override
IndicatorCuboid2D< S > & getIndicator()
Vector< S, 2 > surfaceNormal(const Vector< S, 2 > &pos, const S meshSize) override
const S signedDistance(const Vector< S, 2 > input) override
Vector< S, 2 > calcMofiAndMass(const S density) override
SmoothIndicatorCuboid2D(IndicatorCuboid2D< S > &ind, S epsilon, S theta=0)
bool regardCell(int input[2])
Vector< S, 2 > surfaceNormal(const Vector< S, 2 > &pos, const S meshSize) override
Vector< T, 2 > getLocalCenter()
bool operator()(T output[], const S input[]) override
SmoothIndicatorCustom2D(T latticeSpacing, std::shared_ptr< IndicatorF2D< T > > indPtr, Vector< T, 2 > pos, T epsilon, T theta=0.)
const S signedDistance(const Vector< S, 2 > input) override
Vector< T, 2 > calcMofiAndMass(T rhoP) override
calculates and returns mofi and mass of the particle (mofi is at index 0 and mass at index 1)
Vector< S, 2 > calcMofiAndMass(const S density) override
SmoothIndicatorFactoredCircle2D(Vector< S, 2 > center, S radius, S epsilon, S density=0, Vector< S, 2 > vel=Vector< S, 2 >(0., 0.), S omega=0, S factor=1.)
bool operator()(T output[], const S input[]) override
Vector< S, 2 > calcMofiAndMass(const S density) override
SmoothIndicatorFactoredCuboid2D(Vector< S, 2 > center, S xLength, S yLength, S epsilon, S density=0, Vector< S, 2 > vel=Vector< S, 2 >(0., 0.), S omega=0, S factor=1.)
bool operator()(T output[], const S input[]) override
Vector< S, 2 > calcMofiAndMass(const S density) override
SmoothIndicatorHTCircle2D(Vector< S, 2 > center, S radius, S epsilon, S density=0, Vector< S, 2 > vel=Vector< S, 2 >(0., 0.), S omega=0)
bool operator()(T output[], const S input[]) override
implements a smooth triangle in 2D with an _epsilon sector
Vector< S, 2 > surfaceNormal(const Vector< S, 2 > &pos, const S meshSize) override
IndicatorEquiTriangle2D< S > & getIndicator()
SmoothIndicatorTriangle2D(IndicatorEquiTriangle2D< S > &indPtr, S epsilon, S theta=0)
Vector< S, 2 > calcMofiAndMass(const S density) override
bool distance(S &distance, const Vector< S, 2 > &origin, const Vector< S, 2 > &direction, S precision, S pitch) override
const S signedDistance(const Vector< S, 2 > input) 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.