SmoothIndicatorEllipsoid3D

[guojuw]

Dear OpenLB developers:

I am following the code “${olbRoot}/examples/particles/settlingCube3d/” to develop a ellipsoid particle sedimentation case.

1. I add these lines in “${olbRoot}/src/functors/analytical/indicator/smoothIndicatorF3D.h”:

/// implements a smooth particle Ellipsoid in 3D with an _epsilon sector.
template <typename T, typename S, bool HLBM=false>
class SmoothIndicatorEllipsoid3D final: public SmoothIndicatorF3D<T, S, HLBM> {
private:
S _xLength;
S _yLength;
S _zLength;
public:
SmoothIndicatorEllipsoid3D(Vector<S,3> center, S xLength, S yLength, S zLength, S epsilon, Vector<S,3> theta = Vector<S,3> (0.,0.,0.), S density=0, Vector<S,3> vel = Vector<S,3> (0.,0.,0.));
bool operator()(T output[],const S input[]) override;
};

2. xLength, yLength, zLenght are the diameters of the ellipsoid.

3. In “${olbRoot}/src/functors/analytical/indicator/smoothIndicatorF3D.hh”, I need to write the constructor “SmoothIndicatorEllipsoid3D<T,S,HLBM>::SmoothIndicatorEllipsoid3D(…)” and a function “bool SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator()(T output[], const S input[])”.

4. I see the constructor “SmoothIndicatorEllipsoid3D<T,S,HLBM>::SmoothIndicatorEllipsoid3D(…)” is called only once in the main program.

5. The function “bool SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator()(T output[], const S input[])” is called a lot of times during the main loop.

Would you mind to let me know what is the purpose of the function “bool SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator()(T output[], const S input[])”, and what are the variables of Input[] and output?

Regards,
Junwei Guo (PhD candidate at University of Calgary)

[Nicolas]

Dear Junwei,

before answering your questions, I want to inform you, that we are planning on releasing a new version of OpenLB (v1.4) next month. This version will fortunately include a SmothIndicatorEllipsoid3D, which you might use instead of your own implementation. Nonetheless I am going to refer to your questions in the following:

The operator is the core function used in the functor concept. It expects an input in lattice coordinates and will eventually return a value, which is to be specified in the operator function’s body and usually depends on the input itself. This enables a simple way of implementing e.g. mathematical functions, filters or masks, which can be applied to your calculation domain. For any further questions regarding the functor concept, please refer to section 2.5 in the official OpenLB user guide.

Best,

Nicolas

[guojuw]

Dear Nicolas,

Thanks again. Looking forward to seeing your new releases. I worked out the SmoothIndicatorEllipsoid3D these few days (need to adjust the epsilon at HLBM object bounds). I will compare my code with your official release.

Regards,
Junwei Guo

//Constructor: SmoothIndicatorEllipsoid3D
template <typename T, typename S, bool HLBM>
SmoothIndicatorEllipsoid3D<T,S,HLBM>::SmoothIndicatorEllipsoid3D(Vector<S,3> center, S xLength, S yLength, S zLength, S epsilon, Vector<S,3> theta, S density, Vector<S,3> vel)
: _xLength(xLength),_yLength(yLength),_zLength(zLength)
{
this->_pos = center;
T xr = xLength*.5;
T yr = yLength*.5;
T zr = zLength*.5;
T maxr = max({xr,yr,zr});
//this->_circumRadius = .5*(sqrt(xLength*xLength+yLength*yLength+zLength*zLength))+0.5*epsilon;
this->_circumRadius = maxr + 0.5*epsilon;
this->_name = “1000”; //Ellipsoid particle code;
this->_myMin = {
center[0] – this->getCircumRadius(),
center[1] – this->getCircumRadius(),
center[2] – this->getCircumRadius()
};
this->_myMax = {
center[0] + this->getCircumRadius(),
center[1] + this->getCircumRadius(),
center[2] + this->getCircumRadius()
};
this->_epsilon = epsilon;
this->_theta = {
theta[0] * (M_PI/180.),
theta[1] * (M_PI/180.),
theta[2] * (M_PI/180.)
};
//T mass = xLength*yLength*zLength*density;
T mass = 4./3.*M_PI*xLength*yLength*zLength*density/8.;

T xLength2 = xLength*xLength;
T yLength2 = yLength*yLength;
T zLength2 = zLength*zLength;
//https://en.wikipedia.org/wiki/List_of_moments_of_inertia
Vector<S,3> mofi;
mofi[0] = 1./5.*mass*(yLength2+zLength2)/4.;
mofi[1] = 1./5.*mass*(xLength2+zLength2)/4.;
mofi[2] = 1./5.*mass*(yLength2+xLength2)/4.;
this->init(this->_theta, vel, mass, mofi);
}

template <typename T, typename S, bool HLBM>
bool SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator()(T output[], const S input[])
{
//1.Calculate distance between point and center of unrotated indicator
//cout << “SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator(). Input=(” <<input[0]<<“,”<<input[1]<<“,”<< input[2]<<“). Output=(“<<output[0]<<“).”<< endl;
//cout << ” ” <<input[0]<<“, “<<input[1]<<“, “<< input[2]<<“, “;
T xDist = input[0] – this->getPos()[0];
T yDist = input[1] – this->getPos()[1];
T zDist = input[2] – this->getPos()[2];

//2.Calculate point projected to rotated indicator
// counter-clockwise rotation by _theta=-theta around center

// rotating the object BACK to normal??
// T x = this->getPos()[0] + this->getRotationMatrix()[0]*xDist + this->getRotationMatrix()[3]*yDist + this->getRotationMatrix()[6]*zDist;
// T y = this->getPos()[1] + this->getRotationMatrix()[1]*xDist + this->getRotationMatrix()[4]*yDist + this->getRotationMatrix()[7]*zDist;
// T z = this->getPos()[2] + this->getRotationMatrix()[2]*xDist + this->getRotationMatrix()[5]*yDist + this->getRotationMatrix()[8]*zDist;

T x = this->getPos()[0]*.0 + this->getRotationMatrix()[0]*xDist + this->getRotationMatrix()[3]*yDist + this->getRotationMatrix()[6]*zDist;
T y = this->getPos()[1]*.0 + this->getRotationMatrix()[1]*xDist + this->getRotationMatrix()[4]*yDist + this->getRotationMatrix()[7]*zDist;
T z = this->getPos()[2]*.0 + this->getRotationMatrix()[2]*xDist + this->getRotationMatrix()[5]*yDist + this->getRotationMatrix()[8]*zDist;

T xp = input[0];
T yp = input[1];
T zp = input[2];
T xc = this->getPos()[0];
T yc = this->getPos()[1];
T zc = this->getPos()[2];
T xr = _xLength/2.;
T yr = _yLength/2.;
T zr = _zLength/2.;
T unrotatedInside = ((xp-xc)/xr)*((xp-xc)/xr) + ((yp-yc)/yr)*((yp-yc)/yr) +((zp-zc)/zr)*((zp-zc)/zr) -1.0;

T p1 = x/xr;
T p2 = y/yr;
T p3 = z/zr;
T rotatedInside = p1*p1+p2*p2+p3*p3-1.0;

if (rotatedInside<.0)
{
output[0]=1.0;
return true;
}
output[0]=.0;
return false;
}

[Author]

Posts