Index: /issm/trunk/src/c/Bamgx/GeometricalEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/GeometricalEdge.cpp (revision 2809) +++ /issm/trunk/src/c/Bamgx/GeometricalEdge.cpp (revision 2809) @@ -0,0 +1,123 @@ +#include +#include +#include +#include +#include + +#include "Mesh2.h" +#include "QuadTree.h" +#include "SetOfE4.h" +using namespace std; + +namespace bamg { + + /*Constructor/Destructor*/ + + /*Others*/ + /*FUNCTION GeometricalEdge::R1tg{{{1*/ + Real8 GeometricalEdge::R1tg(Real8 theta,R2 & t) const // 1/R of radius of cuvature + { R2 A=v[0]->r,B=v[1]->r; + Real8 dca,dcb,dcta,dctb; + Real8 ddca,ddcb,ddcta,ddctb; + // Real8 t1 = 1 -theta; + // Real8 t1t1 = t1*t1; + Real8 tt = theta*theta; + assert( theta >=0); + assert( theta <=1); + if (TgA()) + if (TgB()) // interpolation d'hermite + { //cb = theta*theta*(3-2*theta); + dcb = 6*theta*(1-theta); + ddcb = 6*(1-2*theta); + //ca = 1-cb; + dca = -dcb; + ddca = -ddcb; + + // cta = (1-theta)*(1-theta)*theta; + dcta = (3*theta - 4)*theta + 1; + ddcta=6*theta-4; + + //ctb = (theta-1)*theta*theta ; + dctb = 3*tt - 2*theta; + ddctb = 6*theta-2; + } + else { // 1-t*t, t-t*t, t*t + Real8 t = theta; + // cb = t*t; + dcb = 2*t; + ddcb = 2; + //ca = 1-cb; + dca = -dcb; + ddca = -2; + // cta= t-cb; + dcta = 1-dcb; + ddcta = -ddcb; + // ctb =0; + dctb=0; + ddctb=0; + } + else + if (TgB()){ + Real8 t = 1-theta; + //ca = t*t; + dca = -2*t; + ddca = 2; + //cb = 1-ca; + dcb = -dca; + ddcb = -2; + //ctb= -t+ca; + dctb = 1+dca; + ddctb= ddca; + //cta=0; + dcta =0; + ddcta =0; + } + else {t=B-A;return 0;} // lagrange P1 + R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; + + R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; + Real8 d2=(d,d); + Real8 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::F{{{1*/ + R2 GeometricalEdge::F(Real8 theta) const // parametrization of the curve edge + { R2 A=v[0]->r,B=v[1]->r; + Real8 ca,cb,cta,ctb; + assert( theta >=-1e-12); + assert( 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 ; + // if(ref==4 || ref==5) + // cout << " FFF " << tg[0] << tg[1] << A << B << " => " << A*ca + B*cb + tg[0]* cta + tg[1] * ctb << endl; + } + else { // 1-t*t, t-t*t, t*t + Real8 t = theta; + cb = t*t; + ca = 1-cb; + cta= t-cb; + ctb=0; + } + else + if (TgB()){ + Real8 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*/ + +} Index: /issm/trunk/src/c/Bamgx/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Geometry.cpp (revision 2808) +++ /issm/trunk/src/c/Bamgx/Geometry.cpp (revision 2809) @@ -433,3 +433,599 @@ /*}}}1*/ + /*Other*/ + +void Geometry::EmptyGeometry() // empty geometry + { + OnDisk=0; + NbRef=0; + name =0; + quadtree=0; + curves=0; + // edgescomponante=0; + triangles=0; + edges=0; + vertices=0; + NbSubDomains=0; + // nbtf=0; +// BeginOfCurve=0; + nbiv=nbv=nbvx=0; + nbe=nbt=nbtx=0; + NbOfCurves=0; +// BeginOfCurve=0; + subdomains=0; + MaximalAngleOfCorner = 10*Pi/180; + } + + + +Geometry::Geometry(const Geometry & Gh) + { Int4 i; + *this = Gh; + NbRef =0; + quadtree=0; + name = new char[strlen(Gh.name)+4]; + strcpy(name,"cp:"); + strcat(name,Gh.name); + vertices = nbv ? new GeometricalVertex[nbv] : NULL; + triangles = nbt ? new Triangle[nbt]:NULL; + edges = nbe ? new GeometricalEdge[nbe]:NULL; + curves= NbOfCurves ? new Curve[NbOfCurves]:NULL; + subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL; + for (i=0;iAdj[0]; + else if ( (AB,BP) > 0) + on = on->Adj[1]; + else + return on; + } + return on; +} +GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const + { + Real8 save_s=s; + int NbTry=0; +retry: + s=save_s; + GeometricalEdge * on = e.on; + assert(on); + assert( e[0].on && e[1].on); + const Vertex &v0=e[0],&v1=e[1]; + V.m = Metric(1.0-s, v0,s, v1); +#define MXE__LINE __LINE__+1 + const int mxe =100; + GeometricalEdge *ge[mxe+1]; + int sensge[mxe+1]; + Real8 lge[mxe+1]; + int bge=mxe/2,tge=bge; + ge[bge] = e.on; + sensge[bge]=1; + + R2 V0 = v0,V1=v1,V01=V1-V0; + VertexOnGeom vg0= *v0.on, vg1=*v1.on; + if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl; + + // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; + GeometricalEdge * eg0=on, *eg1=on; + R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; + if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) < " << v1.r << endl + << "sg 0 = " << vg0 + << " on = " << Number(on) << ":" << Ag << Bg << "; " + << " sg 1= " << vg1 + << "--------------------------------------------" << endl; + while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0) + { + if (bge<=0) { + // int kkk; + // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ; + if(NbTry) + { + cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; + cerr << " The mesh of the Geometry is to fine: "; + cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; + cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; + cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; + cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; + MeshError(222); + } + NbTry++; + goto retry;} + GeometricalEdge* tmpge = eg0; + if(NbTry) + cout << "bug: --Edge @" << Number(tmpge) << " = "<< Number(eg0) << ":" <Adj[0]) << "," << + Number(eg0->Adj[1]) <<"," ; + ge[--bge] =eg0 = eg0->Adj[sens0]; + assert(bge>=0 && bge <= mxe); + sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); + if(NbTry) + cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S " + << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << "," + << Number(eg0->Adj[1]) <<"," << endl + <=mxe ) { + cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; + NbTry++; + if (NbTry<2) goto retry; + cerr << " The mesh of the Geometry is to fine:" ; + cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; + cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; + cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; + cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; + MeshError(223); + } + + GeometricalEdge* tmpge = eg1; + if(NbTry) + cout << "++Edge @" << tmpge << " = " << Number(eg1) <<"%" << Number(eg1->Adj[0]) << "," + << Number(eg1->Adj[1]) <<"," ; + ge[++tge] =eg1 = eg1->Adj[sens1]; + sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; + assert(tge>=0 && tge <= mxe); + if(NbTry) + cout << " Edge " << Number(eg1) << " " << sens1 << " S " + <Adj[0]) << "," << Number(eg1->Adj[1]) <<"," + <=0 && i <= mxe); + BB = (*ge[i])[sensge[i]]; + lge[i]=ll += Norme2(AA-BB); + // cout << " ll " << i << BB << ll << " " <r - vertices[i]) < eps ) + { // link v & vertices[i] + // vieille ruse pour recuperer j + GeometricalVertex * vg = (GeometricalVertex *) (void *) v; + int j = vg-v0g; + assert( v == & (Vertex &) vertices[j]); + vertices[i].link = vertices + j; + k++; + } + else quadtree.Add(vertices[i]); + } + if (k) { + cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl; + //if (verbosity>10) + { + cout << " The duplicate vertex " << endl; + for (i=0;i7) + for (i=0;ir - edges[i].v[0]->r; + Real8 lv10 = Norme2(v10); + if(lv10 == 0) { + cerr << "The length of " <9) + cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<=0) + { + ord++; + register Int4 i1= n /2; + register Int4 j1 = n % 2; + register Int4 *pn = ev + n; + float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; + n = *pn; + if (angleold > angle) // exch to have : po -> pn -> p + exch=1,*pn = *po,*po=*p,*p=n,po = pn; + else // to have : po -> p -> pn + angleold = angle, po = p,p = pn; + } + } // end while (exch) + + if (ord >= 1 ) + { /* + Int4 n = hv[i]; + while ( n >=0) + { Int4 i1 = n/2,j1 = n%2; + //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi; + n = ev[n]; + } + */ + } + if(ord == 2) { // angulare test to find a corner + Int4 n1 = hv[i]; + Int4 n2 = ev[n1]; + Int4 i1 = n1 /2, i2 = n2/2; // edge number + Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge + float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; + float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; + float da12 = Abs(angle2-angle1); + if(verbosity>9) + cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi + << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl; + + if (( da12 >= MaximalAngleOfCorner ) + && (da12 <= 2*Pi -MaximalAngleOfCorner)) { + vertices[i].SetCorner() ; + if(verbosity>7) + cout << " The vertex " << i << " is a corner (angle) " + << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;} + // if the ref a changing then is SetRequired(); + + if (edges[i1].flag != edges[i2].flag || edges[i1].Required()) + { + vertices[i].SetRequired(); + if(verbosity>7) + cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;} + + if (edges[i1].ref != edges[i2].ref) { + vertices[i].SetRequired(); + if(verbosity>7) + cout << " The vertex " << i << " is Required ref" << endl;} + } ; + + if(ord != 2) { + vertices[i].SetCorner(); + if(verbosity>7) + cout << " the vertex " << i << " is a corner ordre = " << ord << endl; + } + // close the liste around the vertex + { Int4 no=-1, ne = hv[i]; + while ( ne >=0) + ne = ev[no=ne]; + if(no>=0) + ev[no] = hv[i]; + } // now the list around the vertex is circular + + } // end for (i=0;i " + << Number(edges[i1].v[j1]) <10) + cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl; + } + + // generation of all the tangente + for (i=0;ir -edges[i].v[0]->r; + Real8 lAB = Norme2(AB); // length of current edge AB + Real8 ltg2[2]; + ltg2[0]=0;ltg2[1]=0; + for (jj=0;jj<2;jj++) { + R2 tg = edges[i].tg[jj]; + Real8 ltg = Norme2(tg); // length of tg + if(ltg == 0) {// no tg + if( ! edges[i].v[jj]->Corner()) { // not a Corner + tg = edges[i].v[1-jj]->r + - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r; + ltg = Norme2(tg); + tg = tg *(lAB/ltg),ltg=lAB; + /* + if(edges[i].ref >=4) + cout << " tg " << tg.x << " "<< tg.y << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = " + << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl; + */ + } + + //else ;// a Corner with no tangent => nothing to do + } // a tg + else + tg = tg *(lAB/ltg),ltg=lAB; + ltg2[jj] = ltg; + if ( (tg,AB) < 0) + tg = -tg; + //if(edges[i].ref >=4) cout << " tg = " << tg << endl; + edges[i].tg[jj] = tg; + } // for (jj=0;jj<2;jj++) + + if (ltg2[0]!=0) edges[i].SetTgA(); + if (ltg2[1]!=0) edges[i].SetTgB(); + } // for (i=0;i7) + for (i=0;iSetMark(); + nbgem++; + e->CurveNumber=NbOfCurves; + if(curves) { + curves[NbOfCurves].ee=e; + curves[NbOfCurves].ke=k1; + } + + GeometricalVertex *b=(*e)(k1); + if (a == b || b->Required() ) break; + k0 = e->SensAdj[k1];// vertex in next edge + e = e->Adj[k1]; // next edge + + }// for(;;) + if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl; + NbOfCurves++; + if(level) { + if(verbosity>4) + cout << " Warning: Curve "<< NbOfCurves << " without required vertex " + << "so the vertex " << Number(a) << " become required " <SetRequired(); + } + + }} + assert(nbgem && nbe); + + if(step==0) { + curves = new Curve[NbOfCurves]; + } + } + for(int i=0;ilink; + //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; + curves[i].master=true; + if(be->Equi() || be->ReverseEqui() ) + { + assert(eqbe); + int nc = eqbe->CurveNumber; + assert(i!=nc); + curves[i].next=curves[nc].next; + curves[i].master=false; + curves[nc].next=curves+i; + if(be->ReverseEqui()) + curves[i].Reverse(); + } + } + + if(verbosity>3) + cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <4) + for(int i=0;i -#include -#include -#include -#include -using namespace std; - - -#include "Mesh2.h" -#include "QuadTree.h" -#include "SetOfE4.h" -namespace bamg { -void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo - { - // if equiedges existe taille nbe - // equiedges[i]/2 == i original - // equiedges[i]/2 = j => equivalence entre i et j => meme maillage - // equiedges[i]%2 : 0 meme sens , 1 pas meme sens - // - // -------------------------- - long int verbosity=0; - if (verbosity>1) - cout << " -- construction of the geometry from the 2d mesh " << endl; - if (nbt<=0 || nbv <=0 ) { MeshError(101);} - - // construction of the edges - // Triangles * OldCurrentTh =CurrentTh; - CurrentTh=this; - // Int4 NbTold = nbt; - // generation of the integer coor - // generation of the adjacence of the triangles - if (cutoffradian>=0) - Gh.MaximalAngleOfCorner = cutoffradian; - SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv); - Int4 * st = new Int4[nbt*3]; - Int4 i,k; - int j; - if (Gh.name) delete Gh.name; - Gh.name = new char [ name ? strlen(name) + 15 : 50 ]; - Gh.name[0]=0; - strcat(Gh.name,"cons from: "); - if (name) strcat(Gh.name,name); - else strcat(Gh.name," a mesh with no name"); - for (i=0;iaddtrie(Number(edges[i][0]),Number(edges[i][1])); - } - if (nbe != edge4->nb()) - { - cerr << " Some Double edge in the mesh, the number is " << nbe - << " nbe4=" << edge4->nb() << endl; - MeshError(1002); - } - for (i=0;iaddtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), - Number(triangles[i][VerticesOfTriangularEdge[j][1]])); - Int4 invisible = triangles[i].Hidden(j); - if(st[k]==-1) - st[k]=3*i+j; - else if(st[k]>=0) { - assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); - - triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); - if (invisible) triangles[i].SetHidden(j); - if (knb(); // the total number of edges - delete edge4; - edge4 =0; - - if(verbosity>5) { - if (name) - cout << " On Mesh " << name << endl; - cout << " - The number of Vertices = " << nbv << endl; - cout << " - The number of Triangles = " << nbt << endl; - cout << " - The number of given edge = " << nbe << endl; - cout << " - The number of all edges = " << nbedges << endl; - cout << " - The Euler number = 1-Nb Of Hole = " << nbt-nbedges+nbv << endl; } - - - // check the consistant of edge[].adj and the geometrical required vertex - k=0; - kk=0; - Int4 it; - - for (i=0;i4) - cout << " Construction of the edges " << nbe << endl; - - for (i=0;i=0) // edge alone - { - it = st[i]/3; - j = (int) (st[i]%3); - add=k++; - } - - if (add>=0 && add < nbe) - { - - edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; - edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; - edges[add].on=0; - if (iref(),edges[add].v[1]->ref()); // no a good choice - } - } - assert(k==nbe); - if (edgessave) delete [] edgessave; - } - - // construction of edges[].adj - for (i=0;icolor++; - - for (i=0;icolor,j0; - if(i0<0) - edges[i ].adj[ j ]=0; // Add FH Jan 2008 - if(i0==-1) - v->color=i*2+j; - else if (i0>=0) {// i and i0 edge are adjacent by the vertex v - j0 = i0%2; - i0 = i0/2; - assert( v == edges[i0 ].v[j0]); - edges[i ].adj[ j ] =edges +i0; - edges[i0].adj[ j0] =edges +i ; - assert(edges[i0].v[j0] == v); - // if(verbosity>8) - // cout << " edges adj " << i0 << " "<< j0 << " <--> " << i << " " << j << endl; - v->color = -3;} - } - // now reconstruct the sub domain info - assert(!NbSubDomains); - NbSubDomains=0; - - { - Int4 it; - // find all the sub domain - Int4 *colorT = new Int4[nbt]; - Triangle *tt,*t; - Int4 k; - for ( it=0;it0) - if( ( j=st[level]++) <3) - { - t = &triangles[st[level-1]]; - tt=t->TriangleAdj((int)j); - - if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) - { - colorT[jt]=NbSubDomains; - st[++level]=jt; - st[++level]=0; - k++; - } - } - else - level-=2; - if (verbosity>5) - cout << " Nb of triangles " << k << " of Subdomain " - << NbSubDomains << " " << kolor << endl; - NbSubDomains++; - } - } - if (verbosity> 3) - cout << " The Number of sub domain = " << NbSubDomains << endl; - - Int4 isd; - subdomains = new SubDomain[NbSubDomains]; - for (isd=0;isdcolor != triangles[it].color) && !subdomains[isd=colorT[it]].head) - { - subdomains[isd].head = triangles+it; - subdomains[isd].ref = triangles[it].color; - subdomains[isd].sens = j; // hack - subdomains[isd].edge = 0; - k++; - } - } - assert(k== NbSubDomains); - - delete [] colorT; - - - } - delete [] st; - // now make the geometry - // 1 compress the vertices - Int4 * colorV = new Int4[nbv]; - for (i=0;i3) - cout << " Nb of vertices = " << Gh.nbv << " Nb of edges = " << Gh.nbe << endl; - NbVerticesOnGeomVertex = Gh.nbv; - VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; - NbVerticesOnGeomEdge =0; - VerticesOnGeomEdge =0; - { - Int4 j; - for (i=0;i=0) - { - - Vertex & v = Gh.vertices[j]; - v = vertices[i]; - v.color =0; - VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); - } - - } - edge4= new SetOfEdges4(nbe,nbv); - - Real4 * len = new Real4[Gh.nbv]; - for(i=0;i0); - - Real8 hmin = HUGE_VAL; - int kreq=0; - for (i=0;i=0) { - if(verbosity>9) - cout << " Edges Equi " << i << " <=> " << j << " sens = " << sens << endl; - if( sens==0) - Gh.edges[i].SetEqui(); - else - Gh.edges[i].SetReverseEqui(); - Gh.edges[i].link= & Gh.edges[j]; - //assert(sens==0);// meme sens pour l'instant - } - - } - if(requis) { // correction fevr 2009 JYU ... - Gh.edges[i].v[0]->SetRequired(); - Gh.edges[i].v[1]->SetRequired(); - Gh.edges[i].SetRequired(); // fin modif ... - } - R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; - Real8 l12=Norme2(x12); - hmin = Min(hmin,l12); - - Gh.vertices[j1].color++; - Gh.vertices[j0].color++; - - len[j0]+= l12; - len[j1] += l12; - hmin = Min(hmin,l12); - - Gh.edges[i].ref = edges[i].ref; - - k = edge4->addtrie(i0,i1); - - assert(k == i); - - } - - - for (i=0;i 0) - Gh.vertices[i].m= Metric(len[i] /(Real4) Gh.vertices[i].color); - else - Gh.vertices[i].m= Metric(hmin); - delete [] len; - for (i=0;ifindtrie(i0,i1); - if(k>=0) - { - subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1; - subdomains[i].edge = edges+k; - Gh.subdomains[i].edge = Gh.edges + k; - Gh.subdomains[i].sens = subdomains[i].sens; - Gh.subdomains[i].ref = subdomains[i].ref; - } - else - MeshError(103); - } - - delete edge4; - delete [] colorV; - // -- unset adj - for (i=0;iAdj[0]; - else if ( (AB,BP) > 0) - on = on->Adj[1]; - else - return on; - } - return on; -} -GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const - { - Real8 save_s=s; - int NbTry=0; -retry: - s=save_s; - GeometricalEdge * on = e.on; - assert(on); - assert( e[0].on && e[1].on); - const Vertex &v0=e[0],&v1=e[1]; - V.m = Metric(1.0-s, v0,s, v1); -#define MXE__LINE __LINE__+1 - const int mxe =100; - GeometricalEdge *ge[mxe+1]; - int sensge[mxe+1]; - Real8 lge[mxe+1]; - int bge=mxe/2,tge=bge; - ge[bge] = e.on; - sensge[bge]=1; - - R2 V0 = v0,V1=v1,V01=V1-V0; - VertexOnGeom vg0= *v0.on, vg1=*v1.on; - if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl; - - // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; - GeometricalEdge * eg0=on, *eg1=on; - R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; - if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) < " << v1.r << endl - << "sg 0 = " << vg0 - << " on = " << Number(on) << ":" << Ag << Bg << "; " - << " sg 1= " << vg1 - << "--------------------------------------------" << endl; - while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0) - { - if (bge<=0) { - // int kkk; - // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ; - if(NbTry) - { - cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; - cerr << " The mesh of the Geometry is to fine: "; - cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; - cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; - cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; - cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; - MeshError(222); - } - NbTry++; - goto retry;} - GeometricalEdge* tmpge = eg0; - if(NbTry) - cout << "bug: --Edge @" << Number(tmpge) << " = "<< Number(eg0) << ":" <Adj[0]) << "," << - Number(eg0->Adj[1]) <<"," ; - ge[--bge] =eg0 = eg0->Adj[sens0]; - assert(bge>=0 && bge <= mxe); - sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); - if(NbTry) - cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S " - << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << "," - << Number(eg0->Adj[1]) <<"," << endl - <=mxe ) { - cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; - NbTry++; - if (NbTry<2) goto retry; - cerr << " The mesh of the Geometry is to fine:" ; - cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; - cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; - cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; - cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; - MeshError(223); - } - - GeometricalEdge* tmpge = eg1; - if(NbTry) - cout << "++Edge @" << tmpge << " = " << Number(eg1) <<"%" << Number(eg1->Adj[0]) << "," - << Number(eg1->Adj[1]) <<"," ; - ge[++tge] =eg1 = eg1->Adj[sens1]; - sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; - assert(tge>=0 && tge <= mxe); - if(NbTry) - cout << " Edge " << Number(eg1) << " " << sens1 << " S " - <Adj[0]) << "," << Number(eg1->Adj[1]) <<"," - <=0 && i <= mxe); - BB = (*ge[i])[sensge[i]]; - lge[i]=ll += Norme2(AA-BB); - // cout << " ll " << i << BB << ll << " " <r - vertices[i]) < eps ) - { // link v & vertices[i] - // vieille ruse pour recuperer j - GeometricalVertex * vg = (GeometricalVertex *) (void *) v; - int j = vg-v0g; - assert( v == & (Vertex &) vertices[j]); - vertices[i].link = vertices + j; - k++; - } - else quadtree.Add(vertices[i]); - } - if (k) { - cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl; - //if (verbosity>10) - { - cout << " The duplicate vertex " << endl; - for (i=0;i7) - for (i=0;ir - edges[i].v[0]->r; - Real8 lv10 = Norme2(v10); - if(lv10 == 0) { - cerr << "The length of " <9) - cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<=0) - { - ord++; - register Int4 i1= n /2; - register Int4 j1 = n % 2; - register Int4 *pn = ev + n; - float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; - n = *pn; - if (angleold > angle) // exch to have : po -> pn -> p - exch=1,*pn = *po,*po=*p,*p=n,po = pn; - else // to have : po -> p -> pn - angleold = angle, po = p,p = pn; - } - } // end while (exch) - - if (ord >= 1 ) - { /* - Int4 n = hv[i]; - while ( n >=0) - { Int4 i1 = n/2,j1 = n%2; - //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi; - n = ev[n]; - } - */ - } - if(ord == 2) { // angulare test to find a corner - Int4 n1 = hv[i]; - Int4 n2 = ev[n1]; - Int4 i1 = n1 /2, i2 = n2/2; // edge number - Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge - float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; - float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; - float da12 = Abs(angle2-angle1); - if(verbosity>9) - cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi - << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl; - - if (( da12 >= MaximalAngleOfCorner ) - && (da12 <= 2*Pi -MaximalAngleOfCorner)) { - vertices[i].SetCorner() ; - if(verbosity>7) - cout << " The vertex " << i << " is a corner (angle) " - << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;} - // if the ref a changing then is SetRequired(); - - if (edges[i1].flag != edges[i2].flag || edges[i1].Required()) - { - vertices[i].SetRequired(); - if(verbosity>7) - cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;} - - if (edges[i1].ref != edges[i2].ref) { - vertices[i].SetRequired(); - if(verbosity>7) - cout << " The vertex " << i << " is Required ref" << endl;} - } ; - - if(ord != 2) { - vertices[i].SetCorner(); - if(verbosity>7) - cout << " the vertex " << i << " is a corner ordre = " << ord << endl; - } - // close the liste around the vertex - { Int4 no=-1, ne = hv[i]; - while ( ne >=0) - ne = ev[no=ne]; - if(no>=0) - ev[no] = hv[i]; - } // now the list around the vertex is circular - - } // end for (i=0;i " - << Number(edges[i1].v[j1]) <10) - cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl; - } - - // generation of all the tangente - for (i=0;ir -edges[i].v[0]->r; - Real8 lAB = Norme2(AB); // length of current edge AB - Real8 ltg2[2]; - ltg2[0]=0;ltg2[1]=0; - for (jj=0;jj<2;jj++) { - R2 tg = edges[i].tg[jj]; - Real8 ltg = Norme2(tg); // length of tg - if(ltg == 0) {// no tg - if( ! edges[i].v[jj]->Corner()) { // not a Corner - tg = edges[i].v[1-jj]->r - - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r; - ltg = Norme2(tg); - tg = tg *(lAB/ltg),ltg=lAB; - /* - if(edges[i].ref >=4) - cout << " tg " << tg.x << " "<< tg.y << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = " - << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl; - */ - } - - //else ;// a Corner with no tangent => nothing to do - } // a tg - else - tg = tg *(lAB/ltg),ltg=lAB; - ltg2[jj] = ltg; - if ( (tg,AB) < 0) - tg = -tg; - //if(edges[i].ref >=4) cout << " tg = " << tg << endl; - edges[i].tg[jj] = tg; - } // for (jj=0;jj<2;jj++) - - if (ltg2[0]!=0) edges[i].SetTgA(); - if (ltg2[1]!=0) edges[i].SetTgB(); - } // for (i=0;i7) - for (i=0;iSetMark(); - nbgem++; - e->CurveNumber=NbOfCurves; - if(curves) { - curves[NbOfCurves].ee=e; - curves[NbOfCurves].ke=k1; - } - - GeometricalVertex *b=(*e)(k1); - if (a == b || b->Required() ) break; - k0 = e->SensAdj[k1];// vertex in next edge - e = e->Adj[k1]; // next edge - - }// for(;;) - if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl; - NbOfCurves++; - if(level) { - if(verbosity>4) - cout << " Warning: Curve "<< NbOfCurves << " without required vertex " - << "so the vertex " << Number(a) << " become required " <SetRequired(); - } - - }} - assert(nbgem && nbe); - - if(step==0) { - curves = new Curve[NbOfCurves]; - } - } - for(int i=0;ilink; - //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; - curves[i].master=true; - if(be->Equi() || be->ReverseEqui() ) - { - assert(eqbe); - int nc = eqbe->CurveNumber; - assert(i!=nc); - curves[i].next=curves[nc].next; - curves[i].master=false; - curves[nc].next=curves+i; - if(be->ReverseEqui()) - curves[i].Reverse(); - } - } - - if(verbosity>3) - cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <4) - for(int i=0;ir,B=v[1]->r; - Real8 dca,dcb,dcta,dctb; - Real8 ddca,ddcb,ddcta,ddctb; - // Real8 t1 = 1 -theta; - // Real8 t1t1 = t1*t1; - Real8 tt = theta*theta; - assert( theta >=0); - assert( theta <=1); - if (TgA()) - if (TgB()) // interpolation d'hermite - { //cb = theta*theta*(3-2*theta); - dcb = 6*theta*(1-theta); - ddcb = 6*(1-2*theta); - //ca = 1-cb; - dca = -dcb; - ddca = -ddcb; - - // cta = (1-theta)*(1-theta)*theta; - dcta = (3*theta - 4)*theta + 1; - ddcta=6*theta-4; - - //ctb = (theta-1)*theta*theta ; - dctb = 3*tt - 2*theta; - ddctb = 6*theta-2; - } - else { // 1-t*t, t-t*t, t*t - Real8 t = theta; - // cb = t*t; - dcb = 2*t; - ddcb = 2; - //ca = 1-cb; - dca = -dcb; - ddca = -2; - // cta= t-cb; - dcta = 1-dcb; - ddcta = -ddcb; - // ctb =0; - dctb=0; - ddctb=0; - } - else - if (TgB()){ - Real8 t = 1-theta; - //ca = t*t; - dca = -2*t; - ddca = 2; - //cb = 1-ca; - dcb = -dca; - ddcb = -2; - //ctb= -t+ca; - dctb = 1+dca; - ddctb= ddca; - //cta=0; - dcta =0; - ddcta =0; - } - else {t=B-A;return 0;} // lagrange P1 - R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; - - R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; - Real8 d2=(d,d); - Real8 sd2 = sqrt(d2); - t=d; - if(d2>1.0e-20) {t/=sd2;return Abs(Det(d,dd))/(d2*sd2);} - else return 0; - -} - - -R2 GeometricalEdge::F(Real8 theta) const // parametrization of the curve edge -{ R2 A=v[0]->r,B=v[1]->r; - Real8 ca,cb,cta,ctb; - assert( theta >=-1e-12); - assert( 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 ; - // if(ref==4 || ref==5) - // cout << " FFF " << tg[0] << tg[1] << A << B << " => " << A*ca + B*cb + tg[0]* cta + tg[1] * ctb << endl; - } - else { // 1-t*t, t-t*t, t*t - Real8 t = theta; - cb = t*t; - ca = 1-cb; - cta= t-cb; - ctb=0; - } - else - if (TgB()){ - Real8 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; - -} - -} Index: /issm/trunk/src/c/Bamgx/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Triangles.cpp (revision 2808) +++ /issm/trunk/src/c/Bamgx/Triangles.cpp (revision 2809) @@ -461,3 +461,432 @@ /*}}}1*/ + /*Others*/ + /*FUNCTION Triangles::ConsGeometry{{{1*/ + void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo + { + // if equiedges existe taille nbe + // equiedges[i]/2 == i original + // equiedges[i]/2 = j => equivalence entre i et j => meme maillage + // equiedges[i]%2 : 0 meme sens , 1 pas meme sens + // + // -------------------------- + long int verbosity=0; + if (verbosity>1) + cout << " -- construction of the geometry from the 2d mesh " << endl; + if (nbt<=0 || nbv <=0 ) { MeshError(101);} + + // construction of the edges + // Triangles * OldCurrentTh =CurrentTh; + CurrentTh=this; + // Int4 NbTold = nbt; + // generation of the integer coor + // generation of the adjacence of the triangles + if (cutoffradian>=0) + Gh.MaximalAngleOfCorner = cutoffradian; + SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv); + Int4 * st = new Int4[nbt*3]; + Int4 i,k; + int j; + if (Gh.name) delete Gh.name; + Gh.name = new char [ name ? strlen(name) + 15 : 50 ]; + Gh.name[0]=0; + strcat(Gh.name,"cons from: "); + if (name) strcat(Gh.name,name); + else strcat(Gh.name," a mesh with no name"); + for (i=0;iaddtrie(Number(edges[i][0]),Number(edges[i][1])); + } + if (nbe != edge4->nb()) + { + cerr << " Some Double edge in the mesh, the number is " << nbe + << " nbe4=" << edge4->nb() << endl; + MeshError(1002); + } + for (i=0;iaddtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), + Number(triangles[i][VerticesOfTriangularEdge[j][1]])); + Int4 invisible = triangles[i].Hidden(j); + if(st[k]==-1) + st[k]=3*i+j; + else if(st[k]>=0) { + assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); + + triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); + if (invisible) triangles[i].SetHidden(j); + if (knb(); // the total number of edges + delete edge4; + edge4 =0; + + if(verbosity>5) { + if (name) + cout << " On Mesh " << name << endl; + cout << " - The number of Vertices = " << nbv << endl; + cout << " - The number of Triangles = " << nbt << endl; + cout << " - The number of given edge = " << nbe << endl; + cout << " - The number of all edges = " << nbedges << endl; + cout << " - The Euler number = 1-Nb Of Hole = " << nbt-nbedges+nbv << endl; } + + + // check the consistant of edge[].adj and the geometrical required vertex + k=0; + kk=0; + Int4 it; + + for (i=0;i4) + cout << " Construction of the edges " << nbe << endl; + + for (i=0;i=0) // edge alone + { + it = st[i]/3; + j = (int) (st[i]%3); + add=k++; + } + + if (add>=0 && add < nbe) + { + + edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; + edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; + edges[add].on=0; + if (iref(),edges[add].v[1]->ref()); // no a good choice + } + } + assert(k==nbe); + if (edgessave) delete [] edgessave; + } + + // construction of edges[].adj + for (i=0;icolor++; + + for (i=0;icolor,j0; + if(i0<0) + edges[i ].adj[ j ]=0; // Add FH Jan 2008 + if(i0==-1) + v->color=i*2+j; + else if (i0>=0) {// i and i0 edge are adjacent by the vertex v + j0 = i0%2; + i0 = i0/2; + assert( v == edges[i0 ].v[j0]); + edges[i ].adj[ j ] =edges +i0; + edges[i0].adj[ j0] =edges +i ; + assert(edges[i0].v[j0] == v); + // if(verbosity>8) + // cout << " edges adj " << i0 << " "<< j0 << " <--> " << i << " " << j << endl; + v->color = -3;} + } + // now reconstruct the sub domain info + assert(!NbSubDomains); + NbSubDomains=0; + + { + Int4 it; + // find all the sub domain + Int4 *colorT = new Int4[nbt]; + Triangle *tt,*t; + Int4 k; + for ( it=0;it0) + if( ( j=st[level]++) <3) + { + t = &triangles[st[level-1]]; + tt=t->TriangleAdj((int)j); + + if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) + { + colorT[jt]=NbSubDomains; + st[++level]=jt; + st[++level]=0; + k++; + } + } + else + level-=2; + if (verbosity>5) + cout << " Nb of triangles " << k << " of Subdomain " + << NbSubDomains << " " << kolor << endl; + NbSubDomains++; + } + } + if (verbosity> 3) + cout << " The Number of sub domain = " << NbSubDomains << endl; + + Int4 isd; + subdomains = new SubDomain[NbSubDomains]; + for (isd=0;isdcolor != triangles[it].color) && !subdomains[isd=colorT[it]].head) + { + subdomains[isd].head = triangles+it; + subdomains[isd].ref = triangles[it].color; + subdomains[isd].sens = j; // hack + subdomains[isd].edge = 0; + k++; + } + } + assert(k== NbSubDomains); + + delete [] colorT; + + + } + delete [] st; + // now make the geometry + // 1 compress the vertices + Int4 * colorV = new Int4[nbv]; + for (i=0;i3) + cout << " Nb of vertices = " << Gh.nbv << " Nb of edges = " << Gh.nbe << endl; + NbVerticesOnGeomVertex = Gh.nbv; + VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; + NbVerticesOnGeomEdge =0; + VerticesOnGeomEdge =0; + { + Int4 j; + for (i=0;i=0) + { + + Vertex & v = Gh.vertices[j]; + v = vertices[i]; + v.color =0; + VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); + } + + } + edge4= new SetOfEdges4(nbe,nbv); + + Real4 * len = new Real4[Gh.nbv]; + for(i=0;i0); + + Real8 hmin = HUGE_VAL; + int kreq=0; + for (i=0;i=0) { + if(verbosity>9) + cout << " Edges Equi " << i << " <=> " << j << " sens = " << sens << endl; + if( sens==0) + Gh.edges[i].SetEqui(); + else + Gh.edges[i].SetReverseEqui(); + Gh.edges[i].link= & Gh.edges[j]; + //assert(sens==0);// meme sens pour l'instant + } + + } + if(requis) { // correction fevr 2009 JYU ... + Gh.edges[i].v[0]->SetRequired(); + Gh.edges[i].v[1]->SetRequired(); + Gh.edges[i].SetRequired(); // fin modif ... + } + R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; + Real8 l12=Norme2(x12); + hmin = Min(hmin,l12); + + Gh.vertices[j1].color++; + Gh.vertices[j0].color++; + + len[j0]+= l12; + len[j1] += l12; + hmin = Min(hmin,l12); + + Gh.edges[i].ref = edges[i].ref; + + k = edge4->addtrie(i0,i1); + + assert(k == i); + + } + + + for (i=0;i 0) + Gh.vertices[i].m= Metric(len[i] /(Real4) Gh.vertices[i].color); + else + Gh.vertices[i].m= Metric(hmin); + delete [] len; + for (i=0;ifindtrie(i0,i1); + if(k>=0) + { + subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1; + subdomains[i].edge = edges+k; + Gh.subdomains[i].edge = Gh.edges + k; + Gh.subdomains[i].sens = subdomains[i].sens; + Gh.subdomains[i].ref = subdomains[i].ref; + } + else + MeshError(103); + } + + delete edge4; + delete [] colorV; + // -- unset adj + for (i=0;i