StochMagnet/html/EXPR__Node_8h_source.html

 #ifndef EXPR_Node_H

 #define EXPR_Node_H


 #include "EXPR_Object.h"


 //numerics functions headers

 #include "core_numeric.h"


 class EXPR_Node : public virtual EXPR_Object {


 private:

     //attributes


     //unit bounding box

     static const std::array<tBoolean,24> BOUNDING_BOX_POINTS;


     //approximation for testing the bounding box

     static tReal EPSILON;


     //bounding box after applying the transformation

     std::array<tReal,3> mMinPoint;

     std::array<tReal,3> mMaxPoint;


     //linear transformation of the node

     //F(x)=alpha.x+beta

     //F^{-1}(y)=A.(y-beta) A = alpha^{-1}

     std::array<tReal,9> mAlpha;

     std::array<tReal,9> mA;

     std::array<tReal,3>  mBeta;


     std::array<tReal,3> mWork;


 private:

     //associations

 protected:

     //builders


     EXPR_Node();


     //deleters

     virtual ~EXPR_Node();


 public:

      //MEMORY

     //======

     virtual tMemSize getMemorySize() const override {

         return sizeof(*this)+getContentsMemorySize();

     }


     virtual tMemSize getContentsMemorySize() const override {

         tMemSize mem=CORE_Object::getContentsMemorySize();

         mem+=mBeta.size()*sizeof(tReal);

         mem+=mA.size()*sizeof(tReal);

         mem+=mAlpha.size()*sizeof(tReal);

         mem+=mMinPoint.size()*sizeof(tReal);

         mem+=mMaxPoint.size()*sizeof(tReal);

         return mem;

     }


     //methods

     virtual void copy(const EXPR_Node& node) {

         //copy the canonical bounding box


         memcpy(mMinPoint.data(),node.getBoundingBoxMinPoint().data(),sizeof(tReal)*mMinPoint.size());

         memcpy(mMaxPoint.data(),node.getBoundingBoxMaxPoint().data(),sizeof(tReal)*mMaxPoint.size());


         memcpy(mAlpha.data(),node.getAlpha().data(),sizeof(tReal)*mAlpha.size());

         memcpy(mBeta.data(),node.getBeta().data(),sizeof(tReal)*mBeta.size());


         Inverse(mAlpha,mA);

     }


     //transformation methods

     //=======================


     inline void setAlpha(std::initializer_list<tReal>&& alpha) {

         auto iAlpha=std::cbegin(alpha);

         for(auto& alpha_i:mAlpha) {

             alpha_i=(*iAlpha);

             iAlpha++;

         }


         Inverse(mAlpha,mA);

     }


     inline void setBeta(std::initializer_list<tReal>&& beta) {

         auto iBeta=std::cbegin(beta);

         for(auto& beta_i:mBeta) {

             beta_i=(*iBeta);

             iBeta++;

         }

     }

     inline void addToBeta(std::initializer_list<tReal>&& beta) {

         auto iBeta=std::cbegin(beta);

         for(auto& bi: mBeta) {

             bi+=(*iBeta);

             iBeta++;

         };

     }


     inline const std::array<tReal,3>& getBeta() const {

         return mBeta;

     }

     inline const std::array<tReal,9>& getAlpha() const {

         return mAlpha;

     }


     //affine transformation

     //=====================


     inline void leftComposition(const std::array<tReal,9>& A,

                                 const std::array<tReal,3>& B) {

         //F(X)=(A o Id(x) +B) o (alpha.Id(x)+beta)

         //F(x)=A.(alpha.x+beta)+B=A.alpha.x+A.beta+B

         //beta=A.beta+beta

         //alpha:=A.alpha


         //Y=A.beta

         tReal *iW=mWork.data();//iterator on W

         const tReal *eW=iW;eW+=mWork.size();//end iterator on W


         const tReal *Aik=A.data();//iterator on A


         tReal *iBeta=mBeta.data();//itertaor on mBeta

         const tReal *eBeta=iBeta;eBeta+=mBeta.size();//end of iterator on beta


         while (iW!=eW) {//loop on i in {0,1,2}


             //Wi=0

             (*iW)=0;


             iBeta=mBeta.data();

             while (iBeta!=eBeta) {//loop on k in {0,1,2}


                 (*iW)+=(*Aik)*(*iBeta);


                 //next k

                 Aik++;

                 iBeta++;

             }


             //next i

             iW++;

         }


         //beta=Y+B=A.beta+B

         iBeta=mBeta.data();

         iW=mWork.data();//ietartor on work array W

         const tReal *iB=B.data();

         while (iBeta!=eBeta)  {


             //beta_i=W_i+B_i

             (*iBeta)=(*iW)+(*iB);


             //next i

             iBeta++;

             iW++;

             iB++;

         }


         //Alpha:=A.Alpha

         alphaLeftComposition(A);


     }


     inline void leftComposition(const std::array<tReal,9>& A) {

         //F(X)=(A o Id(x)) o (alpha.Id(x)+beta)

         //F(x)=A.(alpha.x+beta)=A.alpha.x+A.beta

         //beta=A.beta+beta

         //alpha:=A.alpha


         //Y=A.beta


         tReal *iW=mWork.data();//iterator on W

         const tReal *eW=iW;eW+=mWork.size();//end iterator on W


         const tReal *Aik=A.data();//iterator on A


         tReal *iBeta=mBeta.data();//iterator on beta

         const tReal *eBeta=iBeta;eBeta+=mBeta.size();//end of iterator on beta


         while (iW!=eW)  {//loop on i in {0,1,2}

             (*iW)=0;


             //computes W=A.beta

             iBeta=mBeta.data();

             while (iBeta!=eBeta) {//loop on k in {0,1,2}


                 //Wi=Aik.beta_k

                 (*iW)+=(*Aik)*(*iBeta);


                 //next k

                 Aik++;

                 iBeta++;

             }


             //next i

             iW++;

         }


         //beta:=W and W=A.beta

         memcpy(mBeta.data(),mWork.data(),mWork.size()*sizeof(tReal));


         //Alpha:=A.Alpha

         alphaLeftComposition(A);

     }


     //canonical transformations

     //=========================


     inline void rotation(const tReal& theta) {

         tReal c=cos(theta);

         tReal s=sin(theta);

         if (fabs(s)<1.e-12) s=0;

         if (fabs(c)<1.e-12) c=0;

         std::array<tReal,9> rot={0,0,0,0,0,0,0,0,0};

         rot[0]=c;

         rot[1]=-s;

         rot[3]=s;

         rot[4]=c;

         rot[8]=1.;

         leftComposition(rot);

     }

     inline void rotation(const tReal& theta,const tReal& Ax,const tReal& Ay,const tReal& Az) {

          tReal c=cos(theta);

          tReal s=sin(theta);

          if (fabs(s)<1.e-12) s=0;

          if (fabs(c)<1.e-12) c=0;

          std::array<tReal,9> rot={0,0,0,0,0,0,0,0,0};

          tReal N=sqrt(Ax*Ax+Ay*Ay+Az*Az);

          tReal Ux=Ax/N;

          tReal Uy=Ay/N;

          tReal Uz=Az/N;


          rot[0]=Ux*Ux*(1-c)+c;

          rot[1]=Ux*Uy*(1-c)-Uz*s;

          rot[2]=Ux*Uz*(1-c)+Uy*s;


          rot[3]=Uy*Ux*(1-c)+Uz*s;

          rot[4]=Uy*Uy*(1-c)+c;

          rot[5]=Uy*Uz*(1-c)-Ux*s;


          rot[6]=Uz*Ux*(1-c)-Uy*s;

          rot[7]=Uz*Uy*(1-c)+Ux*s;

          rot[8]=Uz*Uz*(1-c)+c;

          leftComposition(rot);

     }

     inline void rotation(const tReal& theta,

                          const tReal& Ax,const tReal& Ay,const tReal& Az,

                          const tReal& Cx,const tReal& Cy,const tReal& Cz) {


          tReal c=cos(theta);

          tReal s=sin(theta);

          if (fabs(s)<1.e-12) s=0;

          if (fabs(c)<1.e-12) c=0;

          std::array<tReal,9> rot={0,0,0,0,0,0,0,0,0};

          tReal N=sqrt(Ax*Ax+Ay*Ay+Az*Az);


          std::array<tReal,3> U;

          tReal &Ux=U[0];

          Ux=Ax/N;

          tReal &Uy=U[1];

          Uy=Ay/N;

          tReal &Uz=U[2];

          Uz=Az/N;


          rot[0]=Ux*Ux*(1-c)+c;

          rot[1]=Ux*Uy*(1-c)-Uz*s;

          rot[2]=Ux*Uz*(1-c)+Uy*s;


          rot[3]=Uy*Ux*(1-c)+Uz*s;

          rot[4]=Uy*Uy*(1-c)+c;

          rot[5]=Uy*Uz*(1-c)-Ux*s;


          rot[6]=Uz*Ux*(1-c)-Uy*s;

          rot[7]=Uz*Uy*(1-c)+Ux*s;

          rot[8]=Uz*Uz*(1-c)+c;


          //P=Rot_{theta,U,C}(M)

          //CP=Rot_{theta,U,0)(CM)

          //P=Rot_{theta,U,0)(M)-Rot_{theta,U,0)(C)+C

          const tReal *Rij=rot.data();

          Ux=Cx;

          Ux-=(*Rij)*Cx;Rij++;

          Ux-=(*Rij)*Cy;Rij++;

          Ux-=(*Rij)*Cz;Rij++;

          Uy=Cy;

          Uy-=(*Rij)*Cx;Rij++;

          Uy-=(*Rij)*Cy;Rij++;

          Uy-=(*Rij)*Cz;Rij++;

          Uz=Cz;

          Uz-=(*Rij)*Cx;Rij++;

          Uz-=(*Rij)*Cy;Rij++;

          Uz-=(*Rij)*Cz;Rij++;


          leftComposition(rot,U);

     }


     inline void translate(const std::array<tReal,3>& T) {

         const tReal *iT=T.data();

         for(auto& Bi:mBeta) {

             Bi+=(*iT);

             iT++;

         }

     }


 private:

     inline static void Inverse(const std::array<tReal,9>& A,std::array<tReal,9>& invA) {

         const tReal *alpha=A.data();


         const tReal &A00=*alpha;alpha++;

         const tReal &A01=*alpha;alpha++;

         const tReal &A02=*alpha;alpha++;

         const tReal &A10=*alpha;alpha++;

         const tReal &A11=*alpha;alpha++;

         const tReal &A12=*alpha;alpha++;

         const tReal &A20=*alpha;alpha++;

         const tReal &A21=*alpha;alpha++;

         const tReal &A22=*alpha;


         tReal *beta=invA.data();

         tReal detA=A00*A11*A22+A01*A12*A20+A02*A10*A21-A02*A11*A20-A12*A21*A00-A22*A01*A10;

         *beta=(A11*A22-A12*A21)/detA;beta++;

         *beta=(A02*A21-A01*A22)/detA;beta++;

         *beta=(A01*A12-A02*A11)/detA;beta++;

         *beta=(A12*A20-A10*A22)/detA;beta++;

         *beta=(A00*A22-A02*A20)/detA;beta++;

         *beta=(A02*A10-A00*A12)/detA;beta++;

         *beta=(A10*A21-A11*A20)/detA;beta++;

         *beta=(A01*A20-A00*A21)/detA;beta++;

         *beta=(A00*A11-A01*A10)/detA;

     }


     inline void alphaLeftComposition(const std::array<tReal,9>& B) {

         //alpha:=B.alpha


         //mA=B.alpha;

         //=>Aij=Bik*alpha_kj


         //mAlpha:=mA


         //A=mA

         tReal *Aij=mA.data();


         const tReal *Bik,*Bi=B.data(),*eBi=Bi+B.size();

         tReal *alpha_j,*alpha_kj;

         const tReal *eAlpha=mAlpha.data()+3;

         const tReal *Bie;

         //A=B.alpha

         while (Bi!=eBi) {//loop on row i in {0,1,2}

             //alpha_.j

             alpha_j=mAlpha.data();//j-column of alpha

             while (alpha_j!=eAlpha) {//loop on colum j in {0,12}

                 (*Aij)=0;


                 Bik=Bi;//i-row of B matrix

                 Bie=Bi;Bie+=3;//end of i-row of B matrix

                 alpha_kj=alpha_j;

                 while (Bik!=Bie) {


                     //Aij=Bik.alpha_kj

                     (*Aij)+=(*Bik)*(*alpha_kj);


                     //next k

                     Bik++;

                     alpha_kj+=3;

                 }


                 //next j

                 Aij++;

                 alpha_j++;

             }

             //next i

             Bi+=3;

         }


         //alpha=mA;

         memcpy(mAlpha.data(),mA.data(),mA.size()*sizeof(tReal));


         //mA=mAlpha^{-1}

         Inverse(mAlpha,mA);

     }


 protected:

     inline void apply(std::array<tReal,3>& P) const {


         tReal Px=P[0];

         tReal Py=P[1];

         tReal Pz=P[2];


         P[0]=mAlpha[0]*(Px)+mAlpha[1]*(Py)+mAlpha[2]*(Pz)+mBeta[0];

         P[1]=mAlpha[3]*(Px)+mAlpha[4]*(Py)+mAlpha[5]*(Pz)+mBeta[1];

         P[2]=mAlpha[6]*(Px)+mAlpha[7]*(Py)+mAlpha[8]*(Pz)+mBeta[2];


     }

     inline tBoolean apply(const std::array<tReal,3>& P,std::array<tReal,3>& Q) const {


         if (P.data()==Q.data()) return false;


         //Y=alpha.P+beta


         //iterators on beta

         const tReal *iB=mBeta.data();

         const tReal *eB=iB;eB+=mBeta.size();


         //ietaror on Q=F(P);

         tReal *iQ=Q.data();


         //iterator on P

         const tReal *iP;//iterator on P

         const tReal *eP=P.data();eP+=P.size();//end iterator on P


         //iterator on alpha

         const tReal *Aki=mAlpha.data();


         while (iB!=eB) {//k loop in {0,1,2}

             //Q=B

             (*iQ)=(*iB);


             //Q=B+A.P

             iP=P.data();

             while (iP!=eP) {//i loop on {0,1,2}


                 //Q+=Aki.Pi

                 (*iQ)+=(*Aki)*(*iP);


                 //next i

                 Aki++;

                 iP++;

             }


             //next k

             iQ++;

             iB++;

         }


         return true;

     }


     inline void applyInverse(std::array<tReal,3>& P) const {


         tReal Px=P[0];

         tReal Py=P[1];

         tReal Pz=P[2];


         P[0]=mA[0]*(Px-mBeta[0])+mA[1]*(Py-mBeta[1])+mA[2]*(Pz-mBeta[2]);

         P[1]=mA[3]*(Px-mBeta[0])+mA[4]*(Py-mBeta[1])+mA[5]*(Pz-mBeta[2]);

         P[2]=mA[6]*(Px-mBeta[0])+mA[7]*(Py-mBeta[1])+mA[8]*(Pz-mBeta[2]);


     }

     inline tBoolean applyInverse(const std::array<tReal,3>& P,

                                  std::array<tReal,3>& Q) const {


         if (P.data()==Q.data()) return false;


         //iterators on beta

         const tReal *iB=mBeta.data();


         //itartor on Q=F^-1(P);

         tReal *iQ=Q.data();

         const tReal *eQ=iQ;eQ+=Q.size();


         //iterator on P

         const tReal *iP;

         const tReal *eP=P.data();eP+=P.size();


         //iterator on alpha

         const tReal *Aki=mA.data();


         //P=F(P)=alpha.Q+beta

         //Q=alpha^{-1}.(P-beta)

         //Q=A.(P-beta)


         while (iQ!=eQ)  {//loop on k in {0,1,2}


             //Q=0

             (*iQ)=0;


             iP=P.data();

             iB=mBeta.data();

             while (iP!=eP) { //loop on i in {0,1,2}


                 //Q=A.(P-beta)

                 (*iQ)+=(*Aki)*((*iP)-(*iB));


                 //next i

                 Aki++;

                 iP++;

                 iB++;

             }


             //next k

             iQ++;

         }


         return true;

     }


     //bounding box methods

     //===================

 public:

     inline static  void SetEpsilon(const tReal& eps) {EPSILON=eps;};


     inline static const tReal& GetEpsilon()  {

         return EPSILON;

     }

     virtual void adimensionize(const tReal& L,std::map<tString,tBoolean>& alreadyComputed);


     virtual void computeBoundingBox(std::map<tString,tBoolean>& alreadyComputed)=0;


     inline const std::array<tReal,3>& getBoundingBoxMinPoint() const {

        return mMinPoint;

     }

     inline const std::array<tReal,3>& getBoundingBoxMaxPoint() const {

        return mMaxPoint;

     }


     inline std::array<tReal,3>& getBoundingBoxMinPoint() {

        return mMinPoint;

     }

     inline std::array<tReal,3>& getBoundingBoxMaxPoint() {

        return mMaxPoint;

     }


 protected:


     void linearTransformBoundingBox(std::array<tReal,3>& P,std::array<tReal,3>& Q) const ;


     //is inside methods

     //=================

 public:

     inline tBoolean isInsideBoundingBox(const std::array<tReal,3>& P) const {

            //test if the point M is in the bounding box

         const tReal *mk=mMinPoint.data();

         const tReal *Mk=mMaxPoint.data();

         for(auto& Pk:P) {

             if ( (EPSILON<(*mk)-(Pk)) || ((Pk)-(*Mk)>EPSILON) ) return false;

             Mk++;

             mk++;

         }


         return true;

     }


     virtual tBoolean isInside(std::array<tReal,3> M) const=0;


 public:

     virtual tString toString() const override;

 };


 #endif


CORE_Object::getContentsMemorySize
virtual tMemSize getContentsMemorySize() const
return nthe memory size of the included associations
Definition: CORE_Object.h:278

EXPR_Node
This class is a base class of the binary tree.
Definition: EXPR_Node.h:16

EXPR_Node::adimensionize
virtual void adimensionize(const tReal &L, std::map< tString, tBoolean > &alreadyComputed)
adimensionize the node
Definition: EXPR_Node.cpp:135

EXPR_Node::getContentsMemorySize
virtual tMemSize getContentsMemorySize() const override
return nthe memory size of the included associations
Definition: EXPR_Node.h:79

EXPR_Node::addToBeta
void addToBeta(std::initializer_list< tReal > &&beta)
add to the constant by column
Definition: EXPR_Node.h:136

EXPR_Node::leftComposition
void leftComposition(const std::array< tReal, 9 > &A, const std::array< tReal, 3 > &B)
compose by the affine function
Definition: EXPR_Node.h:166

EXPR_Node::getAlpha
const std::array< tReal, 9 > & getAlpha() const
get the alpha value
Definition: EXPR_Node.h:153

EXPR_Node::SetEpsilon
static void SetEpsilon(const tReal &eps)
set the tolerance error
Definition: EXPR_Node.h:629

EXPR_Node::getBoundingBoxMaxPoint
std::array< tReal, 3 > & getBoundingBoxMaxPoint()
get the max point bounding box of the node for writing
Definition: EXPR_Node.h:671

EXPR_Node::getBoundingBoxMinPoint
const std::array< tReal, 3 > & getBoundingBoxMinPoint() const
get the min point bounding box of the node for reading
Definition: EXPR_Node.h:652

EXPR_Node::linearTransformBoundingBox
void linearTransformBoundingBox(std::array< tReal, 3 > &P, std::array< tReal, 3 > &Q) const
compute the bounding box after transformation
Definition: EXPR_Node.cpp:44

EXPR_Node::isInside
virtual tBoolean isInside(std::array< tReal, 3 > M) const =0
return true if the point is in the bounding box of the node

EXPR_Node::rotation
void rotation(const tReal &theta, const tReal &Ax, const tReal &Ay, const tReal &Az, const tReal &Cx, const tReal &Cy, const tReal &Cz)
compute the rotation of the geometry with angle theta and axis (Ax,Ay,Az)a and center (Cx,...
Definition: EXPR_Node.h:328

EXPR_Node::getMemorySize
virtual tMemSize getMemorySize() const override
return the memory size of the class
Definition: EXPR_Node.h:67

EXPR_Node::~EXPR_Node
virtual ~EXPR_Node()
delete the class
Definition: EXPR_Node.cpp:40

EXPR_Node::getBoundingBoxMinPoint
std::array< tReal, 3 > & getBoundingBoxMinPoint()
get the min point bounding box of the node for writing
Definition: EXPR_Node.h:665

EXPR_Node::leftComposition
void leftComposition(const std::array< tReal, 9 > &A)
compose by the affine function
Definition: EXPR_Node.h:225

EXPR_Node::getBeta
const std::array< tReal, 3 > & getBeta() const
get the beta value
Definition: EXPR_Node.h:147

EXPR_Node::rotation
void rotation(const tReal &theta, const tReal &Ax, const tReal &Ay, const tReal &Az)
compute the rotation of the geometry with angle theta and axis (Ax,Ay,Az)
Definition: EXPR_Node.h:295

EXPR_Node::applyInverse
tBoolean applyInverse(const std::array< tReal, 3 > &P, std::array< tReal, 3 > &Q) const
apply the inverse transformation to point P
Definition: EXPR_Node.h:567

EXPR_Node::translate
void translate(const std::array< tReal, 3 > &T)
translate the geometry
Definition: EXPR_Node.h:385

EXPR_Node::isInsideBoundingBox
tBoolean isInsideBoundingBox(const std::array< tReal, 3 > &P) const
return true if the point P is in the bounding box
Definition: EXPR_Node.h:694

EXPR_Node::setBeta
void setBeta(std::initializer_list< tReal > &&beta)
set the constant by column
Definition: EXPR_Node.h:126

EXPR_Node::setAlpha
void setAlpha(std::initializer_list< tReal > &&alpha)
set the matrix alpha by row
Definition: EXPR_Node.h:113

EXPR_Node::EXPR_Node
EXPR_Node()
create the class
Definition: EXPR_Node.cpp:18

EXPR_Node::toString
virtual tString toString() const override
return the string representation of this
Definition: EXPR_Node.cpp:145

EXPR_Node::getBoundingBoxMaxPoint
const std::array< tReal, 3 > & getBoundingBoxMaxPoint() const
get the max point bounding box of the node for reading
Definition: EXPR_Node.h:658

EXPR_Node::apply
tBoolean apply(const std::array< tReal, 3 > &P, std::array< tReal, 3 > &Q) const
apply the transformation to point P
Definition: EXPR_Node.h:502

EXPR_Node::apply
void apply(std::array< tReal, 3 > &P) const
apply the transformation to point P
Definition: EXPR_Node.h:485

EXPR_Node::GetEpsilon
static const tReal & GetEpsilon()
get the tolerance error
Definition: EXPR_Node.h:634

EXPR_Node::computeBoundingBox
virtual void computeBoundingBox(std::map< tString, tBoolean > &alreadyComputed)=0
compute the bounding box of the node

EXPR_Node::rotation
void rotation(const tReal &theta)
compute the rotation of the geometry with angle theta and k-axis
Definition: EXPR_Node.h:276

EXPR_Node::applyInverse
void applyInverse(std::array< tReal, 3 > &P) const
apply the inverse transformation to point P
Definition: EXPR_Node.h:550

EXPR_Node::copy
virtual void copy(const EXPR_Node &node)
copy
Definition: EXPR_Node.h:93

EXPR_Object
This class is the base class of all the parser package.
Definition: EXPR_Object.h:27