#include #include #include #include "../../include/include.h" #include "../../shared/Exceptions/exceptions.h" #include "GeometricalEdge.h" #include "Geometry.h" using namespace std; namespace bamg { /*Constructor/Destructor*/ /*Methods*/ /*FUNCTION GeometricalEdge::Cracked{{{1*/ int GeometricalEdge::Cracked() const { return type &1; }/*}}}*/ /*FUNCTION GeometricalEdge::F{{{1*/ R2 GeometricalEdge::F(double theta) const{ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/F)*/ // parametrization of the curve edge R2 A=v[0]->r,B=v[1]->r; double ca,cb,cta,ctb; //Check that theta is in [0 1] ISSMASSERT(theta>-1e-12 && theta<1+1e-12); if (TgA()) if (TgB()) // interpolation d'hermite { cb = theta*theta*(3-2*theta); ca = 1-cb; cta = (1-theta)*(1-theta)*theta; ctb = (theta-1)*theta*theta ; } else { // 1-t*t, t-t*t, t*t double t = theta; cb = t*t; ca = 1-cb; cta= t-cb; ctb=0; } else if (TgB()){ double t = 1-theta; ca = t*t; cb = 1-ca; ctb= -t+ca; cta=0; } else { ca =(1-theta),cb = theta,cta=ctb=0; // lagrange P1 } return A*ca + B*cb + tg[0]* cta + tg[1] * ctb; } /*}}}1*/ /*FUNCTION GeometricalEdge::Mark{{{1*/ int GeometricalEdge::Mark() const { return type &16; }/*}}}*/ /*FUNCTION GeometricalEdge::R1tg{{{1*/ double GeometricalEdge::R1tg(double theta,R2 & t) const{ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/R1tg)*/ // 1/R of radius of cuvature R2 A=v[0]->r,B=v[1]->r; double dca,dcb,dcta,dctb; double ddca,ddcb,ddcta,ddctb; double tt = theta*theta; //check theta ISSMASSERT(theta>=0 && theta<=1); if (TgA()){ if (TgB()){ // Tangent A and B provided: // interpolation d'hermite dcb = 6*theta*(1-theta); ddcb = 6*(1-2*theta); dca = -dcb; ddca = -ddcb; dcta = (3*theta - 4)*theta + 1; ddcta=6*theta-4; dctb = 3*tt - 2*theta; ddctb = 6*theta-2; } else { //Tangent A provided but tangent B not provided // 1-t*t, t-t*t, t*t double t = theta; dcb = 2*t; ddcb = 2; dca = -dcb; ddca = -2; dcta = 1-dcb; ddcta = -ddcb; dctb=0; ddctb=0; } } else{ if (TgB()){ //Tangent B provided but tangent A not provided double t = 1-theta; dca = -2*t; ddca = 2; dcb = -dca; ddcb = -2; dctb = 1+dca; ddctb= ddca; dcta =0; ddcta =0; } else { //Neither thangent A nor tangent B provided // lagrange P1 t=B-A; return 0; } } R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; double d2=(d,d); double sd2 = sqrt(d2); t=d; if(d2>1.0e-20){ t/=sd2; return Abs(Det(d,dd))/(d2*sd2); } else return 0; } /*}}}1*/ /*FUNCTION GeometricalEdge::Required{{{1*/ int GeometricalEdge::Required() { return type &64; }/*}}}*/ /*FUNCTION GeometricalEdge::Set {{{1*/ void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew){ *this = rec; v[0] = GhNew.vertices + Gh.GetId(v[0]); v[1] = GhNew.vertices + Gh.GetId(v[1]); if (Adj[0]) Adj[0] = GhNew.edges + Gh.GetId(Adj[0]); if (Adj[1]) Adj[1] = GhNew.edges + Gh.GetId(Adj[1]); } /*}}}*/ /*FUNCTION GeometricalEdge::SetCracked{{{1*/ void GeometricalEdge::SetCracked() { type |= 1;/*=>1st digit to 1*/ }/*}}}*/ /*FUNCTION GeometricalEdge::SetTgA{{{1*/ void GeometricalEdge::SetTgA() { type |=4; /*=>2d digit to 1*/ }/*}}}*/ /*FUNCTION GeometricalEdge::SetTgB{{{1*/ void GeometricalEdge::SetTgB() { type |=8; /*=> 3d digit to 1*/ }/*}}}*/ /*FUNCTION GeometricalEdge::SetMark{{{1*/ void GeometricalEdge::SetMark() { type |=16;/*=> 4th digiy to 1*/ }/*}}}*/ /*FUNCTION GeometricalEdge::SetUnMark{{{1*/ void GeometricalEdge::SetUnMark() { type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/; }/*}}}*/ /*FUNCTION GeometricalEdge::SetRequired{{{1*/ void GeometricalEdge::SetRequired() { type |= 64;/*=>6th digit to 1*/ }/*}}}*/ /*FUNCTION GeometricalEdge::Tg{{{1*/ int GeometricalEdge::Tg(int i) const { return i==0 ? TgA() : TgB(); }/*}}}*/ /*FUNCTION GeometricalEdge::TgA{{{1*/ int GeometricalEdge::TgA() const { return type &4; }/*}}}*/ /*FUNCTION GeometricalEdge::TgB{{{1*/ int GeometricalEdge::TgB() const { return type &8; }/*}}}*/ }