Reply To: SmoothIndicatorEllipsoid3D

[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;
}