#5189
Junwei Guo
Participant

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->_name = “1000”; //Ellipsoid particle code;
this->_myMin = {
};
this->_myMax = {
};
this->_epsilon = epsilon;
this->_theta = {
theta * (M_PI/180.),
theta * (M_PI/180.),
theta * (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 = 1./5.*mass*(yLength2+zLength2)/4.;
mofi = 1./5.*mass*(xLength2+zLength2)/4.;
mofi = 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<<“,”<<input<<“,”<< input<<“). Output=(“<<output<<“).”<< endl;
//cout << ” ” <<input<<“, “<<input<<“, “<< input<<“, “;
T xDist = input – this->getPos();
T yDist = input – this->getPos();
T zDist = input – this->getPos();

//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() + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;
// T y = this->getPos() + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;
// T z = this->getPos() + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;

T x = this->getPos()*.0 + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;
T y = this->getPos()*.0 + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;
T z = this->getPos()*.0 + this->getRotationMatrix()*xDist + this->getRotationMatrix()*yDist + this->getRotationMatrix()*zDist;

T xp = input;
T yp = input;
T zp = input;
T xc = this->getPos();
T yc = this->getPos();
T zc = this->getPos();
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=1.0;
return true;
}
output=.0;
return false;
}