Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2846) +++ /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2847) @@ -220,140 +220,4 @@ } - - - - Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega) - { - register Vertex * s = this; - Vertex &vP = *s,vPsave=vP; - // if (vP.on) return 0;// Don't move boundary vertex - - register Triangle * tbegin= t , *tria = t , *ttc; - - register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc; - R2 P(s->r),PNew(0,0); - // cout << BTh.quadtree << " " << BTh.quadtree->root << endl; - // assert(BTh.quadtree && BTh.quadtree->root); - do { - k++; - - if (!tria->Hidden(j)) - { - Vertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]]; - - R2 Q = vQ,QP(P-Q); - Real8 lQP = LengthInterpole(vP,vQ,QP); - PNew += Q+QP/Max(lQP,1e-20); - kk ++; - } - ttc = tria->TriangleAdj(j); - jc = NextEdge[tria->NuEdgeTriangleAdj(j)]; - tria = ttc; - j = NextEdge[jc]; - assert(k<2000); - } while ( tbegin != tria); - if (kk<4) return 0; - PNew = PNew/(Real8)kk; - R2 Xmove((PNew-P)*omega); - PNew = P+Xmove; - Real8 delta=Norme2_2(Xmove); - - - // - Icoor2 deta[3]; - I2 IBTh = BTh.toI2(PNew); - - tstart=BTh.FindTriangleContening(IBTh,deta,tstart); - - if (tstart->det <0) - { // outside - double ba,bb; - TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ; - tstart = edge; - vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)); - } - else - { // inside - Real8 aa[3]; - Real8 s = deta[0]+deta[1]+deta[2]; - aa[0]=deta[0]/s; - aa[1]=deta[1]/s; - aa[2]=deta[2]/s; - vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]); - } - - // recompute the det of the triangle - vP.r = PNew; - - vP.i = Th.toI2(PNew); - - Vertex vPnew = vP; - - int ok=1; - int loop=1; - k=0; - while (ok) - { - ok =0; - do { - k++; - double detold = tria->det; - tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]); - if (loop) - { - Vertex *v0,*v1,*v2,*v3; - if (tria->det<0) ok =1; - else if (tria->Quadrangle(v0,v1,v2,v3)) - { - vP = vPsave; - Real8 qold =QuadQuality(*v0,*v1,*v2,*v3); - vP = vPnew; - Real8 qnew = QuadQuality(*v0,*v1,*v2,*v3); - if (qnewdet < detold/2 ) ok=1; - - } - tria->SetUnMarkUnSwap(0); - tria->SetUnMarkUnSwap(1); - tria->SetUnMarkUnSwap(2); - ttc = tria->TriangleAdj(j); - jc = NextEdge[tria->NuEdgeTriangleAdj(j)]; - tria = ttc; - j = NextEdge[jc]; - assert(k<2000); - } while ( tbegin != tria); - if (ok && loop) vP=vPsave; // no move - loop=0; - } - return delta; - } - - -double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) { - // calcul de 4 angles -- - R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d); - R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D); - // Move(A),Line(B),Line(C),Line(D),Line(A); - const Metric & Ma = a; - const Metric & Mb = b; - const Metric & Mc = c; - const Metric & Md = d; - - double lAB=Norme2(AB); - double lBC=Norme2(BC); - double lCD=Norme2(CD); - double lDA=Norme2(DA); - AB /= lAB; BC /= lBC; CD /= lCD; DA /= lDA; - // version aniso - double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB); - double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC); - double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD); - double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA); - double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA)); - if (sinmin<=0) return sinmin; - return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA))); } -} - Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2847) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 2847) @@ -0,0 +1,153 @@ +#include +#include +#include +#include + +#include "../Mesh2.h" +#include "../QuadTree.h" +#include "../SetOfE4.h" + +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +namespace bamg { + + /*Methods*/ + /*FUNCTION Vertex::Smoothing{{{1*/ + Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega){ + register Vertex * s = this; + Vertex &vP = *s,vPsave=vP; + // if (vP.on) return 0;// Don't move boundary vertex + + register Triangle * tbegin= t , *tria = t , *ttc; + + register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc; + R2 P(s->r),PNew(0,0); + // cout << BTh.quadtree << " " << BTh.quadtree->root << endl; + // assert(BTh.quadtree && BTh.quadtree->root); + do { + k++; + + if (!tria->Hidden(j)) + { + Vertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]]; + + R2 Q = vQ,QP(P-Q); + Real8 lQP = LengthInterpole(vP,vQ,QP); + PNew += Q+QP/Max(lQP,1e-20); + kk ++; + } + ttc = tria->TriangleAdj(j); + jc = NextEdge[tria->NuEdgeTriangleAdj(j)]; + tria = ttc; + j = NextEdge[jc]; + assert(k<2000); + } while ( tbegin != tria); + if (kk<4) return 0; + PNew = PNew/(Real8)kk; + R2 Xmove((PNew-P)*omega); + PNew = P+Xmove; + Real8 delta=Norme2_2(Xmove); + + + // + Icoor2 deta[3]; + I2 IBTh = BTh.toI2(PNew); + + tstart=BTh.FindTriangleContening(IBTh,deta,tstart); + + if (tstart->det <0) + { // outside + double ba,bb; + TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ; + tstart = edge; + vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)); + } + else + { // inside + Real8 aa[3]; + Real8 s = deta[0]+deta[1]+deta[2]; + aa[0]=deta[0]/s; + aa[1]=deta[1]/s; + aa[2]=deta[2]/s; + vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]); + } + + // recompute the det of the triangle + vP.r = PNew; + + vP.i = Th.toI2(PNew); + + Vertex vPnew = vP; + + int ok=1; + int loop=1; + k=0; + while (ok) + { + ok =0; + do { + k++; + double detold = tria->det; + tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]); + if (loop) + { + Vertex *v0,*v1,*v2,*v3; + if (tria->det<0) ok =1; + else if (tria->Quadrangle(v0,v1,v2,v3)) + { + vP = vPsave; + Real8 qold =QuadQuality(*v0,*v1,*v2,*v3); + vP = vPnew; + Real8 qnew = QuadQuality(*v0,*v1,*v2,*v3); + if (qnewdet < detold/2 ) ok=1; + + } + tria->SetUnMarkUnSwap(0); + tria->SetUnMarkUnSwap(1); + tria->SetUnMarkUnSwap(2); + ttc = tria->TriangleAdj(j); + jc = NextEdge[tria->NuEdgeTriangleAdj(j)]; + tria = ttc; + j = NextEdge[jc]; + assert(k<2000); + } while ( tbegin != tria); + if (ok && loop) vP=vPsave; // no move + loop=0; + } + return delta; + } + /*}}}1*/ + + /*Intermediary*/ + /*FUNCTION QuadQuality{{{{1*/ + double QuadQuality(const Vertex & a,const Vertex &b,const Vertex &c,const Vertex &d) { + // calcul de 4 angles -- + R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d); + R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D); + // Move(A),Line(B),Line(C),Line(D),Line(A); + const Metric & Ma = a; + const Metric & Mb = b; + const Metric & Mc = c; + const Metric & Md = d; + + double lAB=Norme2(AB); + double lBC=Norme2(BC); + double lCD=Norme2(CD); + double lDA=Norme2(DA); + AB /= lAB; BC /= lBC; CD /= lCD; DA /= lDA; + // version aniso + double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB); + double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC); + double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD); + double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA); + double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA)); + if (sinmin<=0) return sinmin; + return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA))); + } + /*}}}1*/ + +} Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 2846) +++ /issm/trunk/src/c/Makefile.am (revision 2847) @@ -329,4 +329,5 @@ ./Bamgx/objects/Triangles.cpp \ ./Bamgx/objects/Triangle.cpp \ + ./Bamgx/objects/Vertex.cpp \ ./Bamgx/objects/Geometry.cpp \ ./Bamgx/objects/QuadTree.cpp \ @@ -670,4 +671,5 @@ ./Bamgx/objects/Triangles.cpp \ ./Bamgx/objects/Triangle.cpp \ + ./Bamgx/objects/Vertex.cpp \ ./Bamgx/objects/Geometry.cpp \ ./Bamgx/objects/QuadTree.cpp \