Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2740) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2740) @@ -0,0 +1,375 @@ +/*!\file Bamgx + * \brief: use Bamg capabilities. + */ + +#include "./Bamgx.h" + +#undef __FUNCT__ +#define __FUNCT__ "Bamgx" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + + +/*Bamg: */ +#include +#include +#include +#include +#include +#include +#include +#include "Meshio.h" +#include +#include "Mesh2.h" +#include "QuadTree.h" +#include + +using namespace bamg; +using namespace std; +#define NBVMAX 50000 + + +int Bamgx(double** pelementsout,double** pxout,double** pyout, int* pnelout,int* pnodsout, + double* elements, double* x, double* y, int nel, int nods, double* metric, double gradation, int splitpbedge, int nbv){ + + + + int noerr=1; + + /*output: */ + double* elementsout=NULL; + double* xout=NULL; + double* yout=NULL; + int nelout, nodsout; + + /*Bamg variables: */ + double time0, time1,time2,time3,time4; + // 2 way for uses --- + // 1 first mesh + // or adaptation + + // double *sol=0; + // int adaptation = 0 ; + Int4 i; + hinterpole=1; // by def interpolation a h + // int Err=0; + int fileout=0; + int nbvx = NBVMAX; + int iso =0,AbsError=0,nbjacoby=1,allquad=0; + int NoMeshReconstruction=0; + int Rescaling = 1; + double costheta=2; + Real8 cutoffradian=-1; + double anisomax = 1e6; + double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,raison=0,CutOff=1e-5; + int KeepBackVertices=1; + double hminaniso=1e-100; + const double boundmaxsubdiv = 10; + double maxsubdiv=boundmaxsubdiv; + double omega=1.8; + int NbSmooth=3; + double *solMbb =0,*solMBB=0; + int *typesolsBB =0; + Int4 nbsolbb=0,lsolbb=0; + Int4 nbsolBB=0,lsolBB=0; + int SplitEdgeWith2Boundary=0; + int rbbeqMbb=0,rBBeqMBB=0; + int ChoiseHessien = 0; + double power=1; + Triangles *Thr=0,*Thb=0; + + char * fgeom=0,*fmeshback=0,*fmeshout=0,*fmeshr=0,*fmetrix=0,*famfmt=0,*fmsh=0, + *fnopo=0,*fftq=0,*fam=0,*famdba=0,*rbb=0,*rBB=0,*wbb=0,*wBB=0, + *fMbb=0,*foM=0,*fMBB=0; + + long int verbosity =10; + + char *datargv[128] ; + int datargc=1; + datargv[0]= datargv[1]=0;// for create a error if no parameter + + //Recover command line options: */ + nbvx=nbv; + if(splitpbedge)SplitEdgeWith2Boundary=1; + + /*Rest of the code taken directly from bamg.cpp, with our own mods: */ + + // some verification + NoMeshReconstruction= fmeshr !=0; + if (!fmeshback) fmeshback=fmeshr; + + if ( maxsubdiv > boundmaxsubdiv || maxsubdiv <= 1.0) + { + cerr << " -maxsubdiv " << maxsubdiv << " is not in ]1,"<< boundmaxsubdiv << "]" << endl; + exit(3); + } + if (iso) + anisomax=1; + if (!(fmetrix||fMbb)) + NbSmooth=0; // no metric -> no smoothing + if( ! rbb) // to set the rbb filename by default + rbb=fMbb; + if( ! rBB) // to set the rbb filename by default + rBB=fMBB; + + if (fMbb && rbb) + rbbeqMbb = !strcmp(rbb,fMbb); + if (fMBB && rBB) + rBBeqMBB = !strcmp(rBB,fMBB); + + time0 = CPUtime(); + + /*Read background mesh from simple delaunay triangulation: */ + Triangles BTh(elements,x,y,nel,nods,nbvx,cutoffradian); // read the background mesh + + hmin = Max(hmin,BTh.MinimalHmin()); + hmax = Min(hmax,BTh.MaximalHmax()); + + throw ErrorException(__FUNCT__," done"); + + BTh.MakeQuadTree(); + time1 = CPUtime(); + if (fmetrix) + BTh.ReadMetric(fmetrix,hmin,hmax,coef); + else + { // init with hmax + Metric Mhmax(hmax); + for (Int4 iv=0;iv2) + cout << " -- write Metric file " << foM <Gh,Thb,nbvx) // copy the mesh + free space to modification + : new Triangles(nbvx,BTh,KeepBackVertices) // construct a new mesh + )); + if (Thr != &BTh) delete Thr; + + if(costheta<=1.0) + Th.MakeQuadrangles(costheta); + if (allquad) + Th.SplitElement(allquad==2); + if(SplitEdgeWith2Boundary) + Th.SplitInternalEdgeWithBorderVertices(); + Th.ReNumberingTheTriangleBySubDomain(); + time3 = CPUtime(); + + if(verbosity>0){ + cout << " Cpu time for meshing : " << time3-time2 << "s Nb Vertices/s = " << (Th.nbv) /(time3-time2) << endl; + } + if(verbosity>3) Th.ShowHistogram(); + if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega); + if(verbosity>3 && NbSmooth>0) Th.ShowHistogram(); + + Th.BTh.ReMakeTriangleContainingTheVertex(); + + if (fmeshout) Th.Write(fmeshout ,Triangles::BDMesh); + if (famfmt) Th.Write(famfmt ,Triangles::am_fmtMesh); + if (fam) Th.Write(fam ,Triangles::amMesh); + if (famdba) Th.Write(famdba ,Triangles::amdbaMesh); + if (fftq) Th.Write(fftq ,Triangles::ftqMesh); + if (fmsh) Th.Write(fmsh ,Triangles::mshMesh); + if (fnopo) Th.Write(fnopo ,Triangles::NOPOMesh); + + if ( ( rbb && wbb) ||( rBB && wBB)){ // the code for interpolation + + if(verbosity >1) + { + if (rbb ) + cout << " -- interpolation P1 files : " << rbb << " -> " << wbb < " << wBB <precision(15); + //if(fbb) fbb->precision(15); + for(i=0;i0) { // internal point + a[0]= (Real8) dete[0]/ tb.det; + a[1]= (Real8) dete[1] / tb.det; + a[2] = (Real8) dete[2] / tb.det; + i0=Th.BTh.Number(tb[0]); + i1=Th.BTh.Number(tb[1]); + i2=Th.BTh.Number(tb[2]);} + else { + double aa,bb; + + TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj(); + + int k = ta; + Triangle & tc = *(Triangle *)ta; + i0=Th.BTh.Number(tc[0]); + i1=Th.BTh.Number(tc[1]); + i2=Th.BTh.Number(tc[2]); + a[VerticesOfTriangularEdge[k][1]] =aa; + a[VerticesOfTriangularEdge[k][0]] = bb; + a[OppositeVertex[k]] = 1- aa -bb;} + + Int4 ibb0 = nbfieldbb*i0; + Int4 ibb1 = nbfieldbb*i1; + Int4 ibb2 = nbfieldbb*i2; + + Int4 iBB0 = nbfieldBB*i0; + Int4 iBB1 = nbfieldBB*i1; + Int4 iBB2 = nbfieldBB*i2; + Int4 j=0; + + for ( j=0;j0) + cout << " Cpu time for interpolation " << time04-time00 << " s" <0) + { + cout << " Cpu time for meshing with io : " << time4-time0 + << "s Nb Vertices/s = " << Th.nbv/(time4-time0) + << endl + << " Nb vertices = " << Th.nbv; + if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2) + cout << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2); + if (Th.NbOfQuad ) + cout << " Nb Quadrilaterals = " << Th.NbOfQuad ; + cout << endl; + } + + + // cout << "delete &Th " ; + delete &Th; + //cout << " end " << endl; + } + + throw ErrorException(__FUNCT__," almost there\n"); + + + /*Assign output pointers:*/ + *pelementsout=elementsout; + *pxout=xout; + *pyout=yout; + *pnelout=nelout; + *pnodsout=nodsout; + + return noerr; +} Index: /issm/trunk/src/c/Bamgx/Bamgx.h =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.h (revision 2740) +++ /issm/trunk/src/c/Bamgx/Bamgx.h (revision 2740) @@ -0,0 +1,16 @@ +/*!\file: Bamgx.h + * \brief header file for Bamg module + */ + +#ifndef _BAMGX_H +#define _BAMGX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +int Bamgx(double** pelementsout,double** pxout,double** pyout, int* pnelout,int* pnodsout, + double* elements, double* x, double* y, int nel, int nods, double* metric, double gradation, int splitpbedge, int nbv); + +#endif /* _BAMGX_H */ + + Index: /issm/trunk/src/c/Bamgx/Makefile.am =================================================================== --- /issm/trunk/src/c/Bamgx/Makefile.am (revision 2740) +++ /issm/trunk/src/c/Bamgx/Makefile.am (revision 2740) @@ -0,0 +1,12 @@ +# Makefile using Automake + Autoconf +# ---------------------------------- +# $Id: Makefile.am,v 1.4 2008/11/24 13:22:26 hecht Exp $ +CXXLD=$(STATICTOOL) $(CXX) + +bin_PROGRAMS=$(BAMGPROG) +EXTRA_PROGRAMS=bamg cvmsh2 drawbdmesh +bamg_SOURCES=bamg.cpp +bamg_DEPENDENCIES=../fflib/libff.a + +LDADD= ../fflib/libff.a +INCLUDES=-I../fflib -I../bamglib -I../Graphics Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2740) +++ /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2740) @@ -0,0 +1,5189 @@ +// -*- Mode : c++ -*- +// +// SUMMARY : +// USAGE : +// ORG : +// AUTHOR : Frederic Hecht +// E-MAIL : hecht@ann.jussieu.fr +// + +/* + + This file is part of Freefem++ + + Freefem++ is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + Freefem++ is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with Freefem++; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + */ +#ifdef __MWERKS__ +#ifdef __INTEL__ +//#pragma global_optimizer off +//#pragma inline_depth(0) +//#pragma optimization_level 2 +#endif +//#pragma inline_depth 0 +#endif +extern bool withrgraphique; +#include +#include +#include +#include +#include +using namespace std; + +#include "Mesh2.h" +#include "QuadTree.h" +#include "SetOfE4.h" + +namespace bamg { + + +#ifdef DEBUG1 + extern int SHOW ; // for debugging + int SHOW = 0; // for debugging + +#endif + + int Triangles::counter = 0; + + Triangles * CurrentTh =0; + + int hinterpole=1; + + + long NbUnSwap =0; + int ForDebugging = 0; + const Direction NoDirOfSearch = Direction(); +#ifndef NDEBUG + inline void MyAssert(int i,char*ex,char * file,long line) + { + if( i) { + cerr << "Error Assert:" << ex << " in " << file << " line: " << line << endl; +#ifdef NOTFREEFEM + exit(1); +#else + throw(ErrorExec("exit",1000)); +#endif + } + } +#endif + + Int4 AGoodNumberPrimeWith(Int4 n) + { + const Int4 BigPrimeNumber[] ={ 567890359L, + 567890431L, 567890437L, 567890461L, 567890471L, + 567890483L, 567890489L, 567890497L, 567890507L, + 567890591L, 567890599L, 567890621L, 567890629L , 0}; + + Int4 o = 0; + Int4 pi = BigPrimeNumber[1]; + for (int i=0; BigPrimeNumber[i]; i++) { + Int4 r = BigPrimeNumber[i] % n; + Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); + if ( o < oo) + o=oo,pi=BigPrimeNumber[i];} + // cout << " AGoodNumberPrimeWith " << n << " " <Number(ta) << "::" + << CurrentTh->Number(ta.ns[0]) << "," + << CurrentTh->Number(ta.ns[1]) << "," + << CurrentTh->Number(ta.ns[2]) << "," + << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} " + << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} " + << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} " + << "]" ; + else + f << "[" + << ta.ns[0] << "," + << ta.ns[1] << "," + << ta.ns[2] << "," + << "{" << ta.at[0] << " " << ta.aa[0] << "} " + << "{" << ta.at[1] << " " << ta.aa[1] << "} " + << "{" << ta.at[2] << " " << ta.aa[2] << "} " + << "]" ; + return f;} + + void swap(Triangle *t1,Int1 a1, + Triangle *t2,Int1 a2, + Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2) + { // swap + // -------------------------------------------------------------- + // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles + // + // sb sb + // / | \ / \ ! + // as1/ | \ /a2 \ ! + // / | \ / t2 \ ! + // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 ! + // \ a1|a2 / \ as1 / + // \ | / \ t1 / + // \ | / as2 \ a1/ + // \ | / \ / + // sa sa + // ------------------------------------------------------------- + int as1 = NextEdge[a1]; + int as2 = NextEdge[a2]; + int ap1 = PreviousEdge[a1]; + int ap2 = PreviousEdge[a2]; +#ifdef DRAWING1 + couleur(0); + t1->Draw(); + t2->Draw(); +#endif +#ifdef DEBUG1 + t1->check(); + t2->check(); +#endif + (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb + (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa + // mise a jour des 2 adjacences externes + TriangleAdjacent taas1 = t1->Adj(as1), + taas2 = t2->Adj(as2), + tas1(t1,as1), tas2(t2,as2), + ta1(t1,a1),ta2(t2,a2); +#ifdef DEBUG + assert( ! ta1.Locked()); + assert( ! ta2.Locked()); +#endif + // externe haut gauche + taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); + // externe bas droite + taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); + // remove the Mark UnMarkSwap + t1->SetUnMarkUnSwap(ap1); + t2->SetUnMarkUnSwap(ap2); + // interne + tas1.SetAdj2(tas2); + + t1->det = det1; + t2->det = det2; + + t1->SetTriangleContainingTheVertex(); + t2->SetTriangleContainingTheVertex(); +#ifdef DEBUG1 + t1->check(); + t2->check(); +#endif +#ifdef DRAWING1 + couleur(1); + t1->Draw(); + t2->Draw(); +#endif +#ifdef DRAWING1 + if( CurrentTh) + CurrentTh->inquire(); +#endif + + } // end swap + + + + + + Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside) + { + CurrentTh=&Th; + assert(&Th); + I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); + Icoor2 dete[3]; + Triangle & tb = *Th.FindTriangleContening(I,dete); + + if (tb.link) + { // internal point in a true triangles + a[0]= (Real8) dete[0]/ tb.det; + a[1]= (Real8) dete[1] / tb.det; + a[2] = (Real8) dete[2] / tb.det; + inside = 1; + return Th.Number(tb); + } + else + { + inside = 0; + double aa,bb; + TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb); + int k = ta; + Triangle * tc = ta; + if (!tc->link) + { ta = ta.Adj(); + tc=ta; + k = ta; + Exchange(aa,bb); + assert(tc->link); + } + a[VerticesOfTriangularEdge[k][0]] = aa; + a[VerticesOfTriangularEdge[k][1]] = bb; + a[OppositeVertex[k]] = 1- aa -bb; + return Th.Number(tc); + } + } + + + TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { + // + // cout << " - "; + int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + int dir=0; + assert(k>=0); + int kkk=0; + Icoor2 IJ_IA,IJ_AJ; + TriangleAdjacent edge(t,OppositeEdge[k]); + for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) + { + + assert(kkk++<1000); + Vertex &vI = *edge.EdgeVertex(0); + Vertex &vJ = *edge.EdgeVertex(1); + I2 I=vI, J=vJ, IJ= J-I; + IJ_IA = (IJ ,(A-I)); + // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <0) {a=1;b=0;return edge;}// change of signe => I + else {dir=-1; + continue;}};// go in direction i + IJ_AJ = (IJ ,(J-A)); + if (IJ_AJ<0) { + if(dir<0) {a=0;b=1;return edge;} + else {dir = 1; + continue;}}// go in direction j + double IJ2 = IJ_IA + IJ_AJ; + assert(IJ2); + a= IJ_AJ/IJ2; + b= IJ_IA/IJ2; + // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl; + return edge; + } + } + + TriangleAdjacent Triangle::FindBoundaryEdge(int i) const + { + // turn around the vertex ns[i] also call s +#ifdef DEBUG + register Vertex * s = ns[i]; +#endif + Triangle *t = (Triangle *) this , *ttc; + int k=0,j = EdgesVertexTriangle[i][0],jc; + int exterieur = !link ; + + do + { + int exterieurp = exterieur; + k++; +#ifdef DEBUG + assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] ); +#endif + ttc = t->at[j]; + exterieur = !ttc->link; + if (exterieur+exterieurp == 1) + return TriangleAdjacent(t,j); + jc = NextEdge[t->aa[j]&3]; + t = ttc; + j = NextEdge[jc]; + assert(k<2000); + } while ( (this!= t)); + return TriangleAdjacent(0,0); + + } + + + TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) + { + // walk around the vertex + // version 2 for remove the probleme if we fill the hole + //int bug=1; + // Triangle *torigine = t; + // restart: + // int dir=0; + assert(t->link == 0); + // to have a starting edges + // try the 3 edge bourna-- in case of internal hole + // and choice the best + // + // + // the probleme is in case of the fine and long internal hole + // for exemple neart the training edge of a wing + // + Vertex * s=0,*s1=0, *s0=0; + Icoor2 imax = MaxICoor22; + Icoor2 l0 = imax,l1 = imax; + double dd2 = imax;// infinity + TriangleAdjacent er; + int cas=-2; + for (int j=0;j<3;j++) + { + TriangleAdjacent ta=t->FindBoundaryEdge(j); + if (! (Triangle *) ta) continue; + s0 = ta.EdgeVertex(0); + s1 = ta.EdgeVertex(1); + I2 A = * s0; + I2 B = *ta.EdgeVertex(1); + I2 AB = B-A,AC=C-A,BC=B-C; + Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); + Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 + Icoor2 ABAC = (AB,AC); // AB.AC| + + double d2; + if ( ABAC < 0 ) // DIST A + { + if ( (d2=(double) ACAC) < dd2) + { + // cout << " A " << d2 << " " << dd2; + er = ta; + l0 = ACAC; + l1 = BCBC; + cas = 0; + s = s0; + } + } + else if (ABAC > AB2) // DIST B + { + if ( (d2=(double) BCBC) < dd2) + { + // cout << " B " << d2 << " " << dd2; + dd2 = d2; + er = Adj(ta); // other direction + l0 = BCBC; + l1 = ACAC; + cas = 1; + s = s1; + } + } + else // DIST AB + { + + double det_2 = (double) Det(AB,AC); + det_2 *= det_2; // square of area*2 of triangle ABC + d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC + // cout << " AB " << d2 << " " << dd2 + // << " " << CurrentTh->Number(ta.EdgeVertex(0)) + // << " " << CurrentTh->Number(ta.EdgeVertex(1)) << " " ; + + if (d2 < dd2) + { + dd2 = d2; + er = ta; + l0 = (AC,AC); + l1 = (BC,BC); + s = 0; + cas = -1; + // cout << " ABAC " << ABAC << " ABAC " << ABAC + // << " AB2 " << AB2 << endl; + b = ((double) ABAC/(double) AB2); + a = 1 - b; + } + } + } + assert(cas !=-2); + // l1 = ||C s1|| , l0 = ||C s0|| + // where s0,s1 are the vertex of the edge er + + if ( s) + { + t=er; + TriangleAdjacent edge(er); + + int kkk=0; + int linkp = t->link == 0; + + Triangle * tt=t=edge=Adj(Previous(edge)); + // cout << CurrentTh->Number(t) << " " << linkp << endl; + do { // loop around vertex s + + assert(edge.EdgeVertex(0)==s && kkk++<10000); + + int link = tt->link == 0; + // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s) + // << " " << CurrentTh->Number(er.EdgeVertex(0)) + // << " " << CurrentTh->Number(er.EdgeVertex(1)) + // << " " << CurrentTh->Number(edge.EdgeVertex(0)) + // << " " << CurrentTh->Number(edge.EdgeVertex(1)) + // << endl; + if ((link + linkp) == 1) + { // a boundary edge + Vertex * st = edge.EdgeVertex(1); + I2 I=*st; + Icoor2 ll = Norme2_2 (C-I); + if (ll < l1) { // the other vertex is neart + s1=st; + l1=ll; + er = edge; + if(llNumber(er.EdgeVertex(0)) + // << " " << CurrentTh->Number(er.EdgeVertex(1)) ; + if (aa<0) a=1,b=0; + else if(bb<0) a=0,b=1; + else + { + a = bb/(aa+bb); + b = aa/(aa+bb); + } + } + + // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " " + // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a + // << " " << b <<" " << l0 << " " <det <0) { // outside + double ba,bb; + TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ; + return 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; + return Metric(aa,(*t)[0],(*t)[1],(*t)[2]); + } + } + + + void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh, + const R2 &A,const R2 &B,int nbegin) + { // SplitEdge + // if(SHOW) cout << " splitedge " << A << B << " " << nbegin << endl; + Triangle *tbegin, *t; + + long int verbosity=2; + Icoor2 deta[3], deti,detj; + Real8 ba[3]; + int nbt =0,ifirst=-1,ilast; + int i0,i1,i2; + int ocut,i,j,k=-1; + // int OnAVertices =0; + Icoor2 dt[3]; + I2 a = Bh.toI2(A) ,b= Bh.toI2(B);// compute the Icoor a,b + I2 vi,vj; + int iedge =-1;// not a edge + + if(nbegin) {// optimisation + // we suppose knowing the starting triangle + t=tbegin=lIntTria[ilast=(Size-1)].t; + if (tbegin->det>=0) + ifirst = ilast;} + else {// not optimisation + init(); + t=tbegin = Bh.FindTriangleContening(a,deta); + // if(SHOW) cout <det<< " " << Real8(deta[1])/t->det + // << " " << Real8(deta[2])/t->det << endl; + if( t->det>=0) + ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det); + else + {// find the nearest boundary edge of the vertex A + // find a edge or such normal projection a the edge IJ is on the edge + // <=> IJ.IA >=0 && IJ.AJ >=0 + ilast=ifirst; + double ba,bb; + TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); + Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); + NewItem(A,Metric(ba,v0,bb,v1)); + t=edge; + // test if the point b is in the same side + if (det(v0.i,v1.i,b)>=0) { + //cout << " All the edge " << A << B << endl; + TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); + Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); + NewItem(A,Metric(ba,v0,bb,v1)); + return; + } + } // find the nearest boundary edge of the vertex A + } // end not optimisation + if (t->det<0) { // outside departure + while (t->det <0) { // intersection boundary edge and a,b, + k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + assert(k>=0); + ocut = OppositeEdge[k]; + i=VerticesOfTriangularEdge[ocut][0]; + j=VerticesOfTriangularEdge[ocut][1]; + vi=(*t)[i]; + vj=(*t)[j]; + deti = bamg::det(a,b,vi); + detj = bamg::det(a,b,vj); + // if(SHOW) { penthickness(3); + // Move(vi);Line(vj);CurrentTh->inquire();penthickness(1); + // cout << Bh.Number(tbegin) << " " << Bh.Number(t) << " i= " << i <<" j= " << j << " k=" << k + // << " deti= " << deti << " detj= " << detj + // << " v = " << Bh.Number((*t)[i]) << (*t)[i].r << " " << Bh.Number((*t)[j]) << (*t)[j].r << endl;} + if (deti>0) // go to i direction on gamma + ocut = PreviousEdge[ocut]; + else if (detj<=0) // go to j direction on gamma + ocut = NextEdge[ocut]; + TriangleAdjacent tadj =t->Adj(ocut); + t = tadj; + iedge= tadj; + if (t == tbegin) { // + double ba,bb; + long int verbosity=2; + if (verbosity>7) + cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b) + << " deti= " << deti << " detj=" <Draw(); + if (iedge < 0) { + i0 =0;i1=1;i2=2; + dt[0] =bamg::det(a,b,(*t)[0]); + dt[1] =bamg::det(a,b,(*t)[1]); + dt[2] =bamg::det(a,b,(*t)[2]);} + else { + i2 = iedge; + i0 = NextEdge[i2]; + i1 = NextEdge[i0]; + dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because + dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side + dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);} + + // so we have just to see the transition from - to + of the det0..2 on edge of t + // because we are going from a to b + if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) + ocut =i2; + else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) + ocut =i2; + else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] == 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] == 0)) + ocut =i2; + else { // On a edge (2 zero) + k =0; + if (dt[0]) ocut=0,k++; + if (dt[1]) ocut=1,k++; + if (dt[2]) ocut=2,k++; + if(k == 1) { + if (dt[ocut] >0) // triangle upper AB + ocut = NextEdge[ocut]; + i= VerticesOfTriangularEdge[ocut][0]; + j= VerticesOfTriangularEdge[ocut][1]; + } + else { + cerr << " Bug Split Edge " << endl; + cerr << " dt[0]= " << dt[0] + << " dt[1]= " << dt[1] + << " dt[2]= "<< dt[2] << endl; + cerr << i0 << " " << i1 << " " << i2 << endl; + cerr << " A = " << A << " B= " << B << endl; + cerr << " Triangle t = " << *t << endl; + cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl; + cerr << " nbt = " << nbt << endl; + MeshError(100);}} + + k = OppositeVertex[ocut]; + + Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b); + + + if (detbij >= 0) { //we find the triangle contening b + dt[0]=bamg::det((*t)[1],(*t)[2],b); + dt[1]=bamg::det((*t)[2],(*t)[0],b); + dt[2]=bamg::det((*t)[0],(*t)[1],b); +#ifdef DEBUG + assert(dt[0] >= 0); + assert(dt[1] >= 0); + assert(dt[2] >= 0); +#endif + Real8 dd = t->det; + NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd); + return ;} + else { // next triangle by adjacent by edge ocut + deti = dt[i]; + detj = dt[j]; + Real4 dij = detj-deti; + ba[i] = detj/dij; + ba[j] = -deti/dij; + ba[3-i-j ] = 0; + ilast=NewItem(t, ba[0],ba[1],ba[2]); + + TriangleAdjacent ta =t->Adj(ocut); + t = ta; + iedge= ta; + if (t->det <= 0) { + double ba,bb; + TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); + NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); + // cout << " return " << ba << " " << bb << endl; + // ajoute le 03 frev 1997 par F. hecht + return; + } + }// we go outside of omega + } // for(;;) + + + } // routine SplitEdge + + + int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) { + register int n; + R2 x(0,0); + if ( d0) x = (*tt)[0].r * d0; + if ( d1) x = x + (*tt)[1].r * d1; + if ( d2) x = x + (*tt)[2].r * d2; + // newer add same point + if(!Size || Norme2_2(lIntTria[Size-1].x-x)) { + if (Size==MaxSize) ReShape(); + lIntTria[Size].t=tt; + lIntTria[Size].bary[0]=d0; + lIntTria[Size].bary[1]=d1; + lIntTria[Size].bary[2]=d2; + lIntTria[Size].x = x; + Metric m0,m1,m2; + register Vertex * v; + if ((v=(*tt)(0))) m0 = v->m; + if ((v=(*tt)(1))) m1 = v->m; + if ((v=(*tt)(2))) m2 = v->m; + lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2); +#ifdef DEBUG1 + if(SHOW) { cout << "SHOW ++ NewItem =" << Size << x ; + cout << " " << d0 << " " << d1 << " " << d2 <0); + // computation of the length + R2 C; + Metric Mx,My; + int ii,jj; + R2 x,y,xy; + + SegInterpolation *SegI=lSegsI; + SegI=lSegsI; + lSegsI[NbSeg].last=Size;// improvement + + int EndSeg=Size; + + y = lIntTria[0].x; + Real8 sxy, s = 0; + lIntTria[0].s =0; + SegI->lBegin=s; + + for (jj=0,ii=1;iilEnd=s, + SegI++, + EndSeg=SegI->last, + SegI->lBegin=s; + + // cout << ii << " " << jj << x<< y <lEnd=s; + + // cout << " len= " << s << endl; + return s; + } + + Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx) + { + long int verbosity=0; + + + const Int4 nbvold = nbv; + Real8 s = Length(); + if (s < 1.5 ) return 0; + ////////////////////// + int ii = 1 ; + R2 y,x; + Metric My,Mx ; + Real8 sx =0,sy; + int nbi = Max(2,(int) (s+0.5)); + Real8 sint = s/nbi; + Real8 si = sint; + + int EndSeg=Size; + SegInterpolation *SegI=0; + if (NbSeg) + SegI=lSegsI,EndSeg=SegI->last; + + for (int k=1;klast; + + int ii1=ii-1; + x =lIntTria[ii1].x; + sx =lIntTria[ii1].s; + Metric Mx=lIntTria[ii1].m; +#ifdef DEBUG + double lx = lIntTria[ii-1].sn; +#endif + y =lIntTria[ii].x; + sy =lIntTria[ii].s; + Metric My=lIntTria[ii].m; +#ifdef DEBUG + double ly =lIntTria[ii].sp; + assert( sx <= si); + assert( si <= sy); + assert( sy != sx); +#endif + + Real8 lxy = sy-sx; + Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy); + + R2 C; + Real8 cx = 1-cy; + C = SegI ? SegI->F(si): x * cx + y *cy; + + si += sint; + if ( nbv1) { + R2 AB = vertices[nbv-2].r - vertices[nbv-1].r ; + Real8 dp = LengthInterpole(vertices[nbv-2].m,vertices[nbv-1].m,AB); + if (dp > 1.6) { + cerr << "PB calcul new Int. points trop loin l=" << dp << " v=" << nbv-1 << " " << nbv-2 <0 et detsb <0 + // attention la routine echange pva et pvb + + if(tt1.Locked()) return 0; // frontiere croise + + TriangleAdjacent tt2 = Adj(tt1); + Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent + Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles + assert ( a1 >= 0 && a1 < 3 ); + + Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; + Vertex & s1= (*t1)[OppositeVertex[a1]]; + Vertex & s2= (*t2)[OppositeVertex[a2]]; + + + Icoor2 dets2 = det(*pva,*pvb,s2); + +#ifdef DEBUG + Vertex & sb= (*t1)[VerticesOfTriangularEdge[a1][1]]; + Icoor2 wdets1 = det(*pva,*pvb,s1); + Icoor2 wdetsa = det(*pva,*pvb,sa); + Icoor2 wdetsb = det(*pva,*pvb,sb); + assert(wdets1 == dets1); + assert(wdetsa == detsa); + assert(wdetsb == detsb); +#endif + + Icoor2 det1=t1->det , det2=t2->det ; +#ifdef DEBUG + assert(det1>0 && det2 >0); + Icoor2 ddet1 = det((*t1)[0],(*t1)[1],(*t1)[2]); + Icoor2 ddet2 = det((*t2)[0],(*t2)[1],(*t2)[2]); + if ((det1 != ddet1) || (det2 != ddet2) ) + { + assert(det1 == ddet1); + assert(det2 == ddet2); + } + Icoor2 detvasasb = det(*pva,sa,sb); + Icoor2 detvbsasb = det(*pvb,sa,sb); + if ( CurrentTh && ! ( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0)))) + { + cout << " detvasasb =" << detvasasb << "detvbsasb = " << detvbsasb + << " " << pva << " " << pvb << " " <Draw(); + penthickness(10); + pva->MoveTo();pvb->LineTo(); + penthickness(1); + CurrentTh->inquire(); +#endif + } + assert( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0))); +#endif + + Icoor2 detT = det1+det2; + assert((det1>0 ) && (det2 > 0)); + assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb + Icoor2 ndet1 = bamg::det(s1,sa,s2); + Icoor2 ndet2 = detT - ndet1; + + int ToSwap =0; //pas de swap + if ((ndet1 >0) && (ndet2 >0)) + { // on peut swaper + if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) + ToSwap =1; + else // swap alleatoire + if (BinaryRand()) + ToSwap =2; + } +#ifdef DEBUG + if (ForDebugging) { + cerr << "swap = " << ToSwap << " ndet1 " << ndet1 << ", ndet2 " << ndet2 << "det1 " << det1 << " det2 " << det2 + << " if1 = " << ((ndet1 >0) && (ndet2 >0)) + << " if2 = " << ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) << endl; +#ifdef DRAWING + couleur(0); + t1->Draw(); + t2->Draw(); +#endif + } +#endif + if (ToSwap) NbSwap++, + bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); + +#ifdef DEBUG + if (ForDebugging) { +#ifdef DRAWING + couleur(4); + t1->Draw(); + t2->Draw(); + rattente(1); +#endif + } +#endif + int ret=1; + + if (dets2 < 0) {// haut + dets1 = ToSwap ? dets1 : detsa ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = ToSwap ? dets1 : detsb ; + detsb = dets2; + //xxxx tt1 = ToSwap ? tt1 : Next(tt2); + if(!ToSwap) tt1 = Next(tt2); + } + else { // changement de sens + if (ForDebugging) cout << "changement de sens" << endl; + ret = -1; + Exchange(pva,pvb); + Exchange(detsa,detsb); + Exchange(dets1,dets2); + Exchange(tt1,tt2); + dets1=-dets1; + dets2=-dets2; + detsa=-detsa; + detsb=-detsb; + + if (ToSwap) + if (dets2 < 0) {// haut + dets1 = (ToSwap ? dets1 : detsa) ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = (ToSwap ? dets1 : detsb) ; + detsb = dets2; + if(!ToSwap) tt1 = Next(tt2); + } + else {// on a fin ??? + tt1 = Next(tt2); + ret =0;} + + } + return ret; + } + + int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) + { +#ifdef DEBUG +restart: // for debug +#endif + int NbSwap =0; + assert(a.t && b.t); // the 2 vertex is in a mesh + int k=0; + taret=TriangleAdjacent(0,0); // erreur + + TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); + Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; + // we turn around a in the direct sens + + Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; + if(v2) // normal case + det2 = det(*v2,a,b); + else { // no chance infini vertex try the next + tta= Previous(Adj(tta)); + v2 = tta.EdgeVertex(0); + vbegin =v2; + assert(v2); + det2 = det(*v2,a,b); + // cout << " No Change try the next" << endl; + } + +#ifdef DRAWING1 + a.MoveTo();b.LineTo(); +#endif + + while (v2 != &b) { + TriangleAdjacent tc = Previous(Adj(tta)); + v1 = v2; + v2 = tc.EdgeVertex(0); + det1 = det2; +#ifdef DEBUG + assert( v1 == tta.EdgeVertex(0)); + assert( &a == tc.EdgeVertex(1) ); +#endif + det2 = v2 ? det(*v2,a,b): det2; + + if((det1 < 0) && (det2 >0)) { + // try to force the edge + Vertex * va = &a, *vb = &b; + tc = Previous(tc); + assert ( v1 && v2); + Icoor2 detss = 0,l=0,ks; + // cout << "Real ForcingEdge " << *va << *vb << detss << endl; +#ifdef DEBUG + Icoor2 dettt1 = det(*v1,a,b); + Icoor2 dettt2 = det(*v2,a,b); + + if (!(dettt1==det1 && dettt2==det2)) + { + assert(ForDebugging==0); + ForDebugging=1; + goto restart; + } + +#endif + while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) + if(l++ > 10000000) { + cerr << " Loop in forcing Egde AB" + <<"\n vertex A " << a + <<"\n vertex B " << b + <<"\n nb de swap " << NbSwap + <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl; + + if ( CurrentTh ) + cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; +#ifdef DEBUG + ForDebugging = 1; +#endif +#ifdef DRAWING1 + if ( CurrentTh ) { + reffecran(); + couleur(6); + CurrentTh->Draw(); + couleur(1); + penthickness(10); + a.MoveTo();b.LineTo(); + penthickness(1); + CurrentTh->inquire(); + couleur(6); + l=0; + reffecran(); + while (ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap) && (l++ < 1000)) + cerr << " " << CurrentTh->Number(tc.EdgeVertex(0))<<" " <Number(tc.EdgeVertex(1)) << " "; + } +#endif + MeshError(990); + } + Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); + if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { + tc.SetLock(); + a.Optim(1,0); + b.Optim(1,0); + taret = tc; + return NbSwap; + } + else + { + taret = tc; + return -2; // error boundary is crossing + /* cerr << "Fatal Error boundary is crossing "; + if(CurrentTh) + { + cerr << " edge: [" << CurrentTh->Number(a) << ", " << CurrentTh->Number(b) << " ] and [ "; + cerr << CurrentTh->Number(aa) << " " << CurrentTh->Number(bb) << " ] " << endl; + } + MeshError(991); + */ + } + } + tta = tc; + assert(k++<2000); + if ( vbegin == v2 ) return -1;// error + } + + tta.SetLock(); + taret=tta; + a.Optim(1,0); + b.Optim(1,0); + return NbSwap; + } + + + int Triangle::swap(Int2 a,int koption){ +#ifdef DEBUG + if(a &4 ) return 0;// arete lock + int munswap1 = a/4; + a &=3; +#else + if(a/4 !=0) return 0;// arete lock or MarkUnSwap +#endif + + register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent + register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles +#ifdef DEBUG + if(a2 & 4) return 0; // arete lock + int munswap2 = a2/4; + a2 &= 3; +#else + if(a2/4 !=0) return 0; // arete lock or MarkUnSwap +#endif + + register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; + register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; + register Vertex *s1=t1->ns[OppositeVertex[a1]]; + register Vertex *s2=t2->ns[OppositeVertex[a2]]; + +#ifdef DEBUG + assert ( a >= 0 && a < 3 ); +#endif + + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + Icoor2 detA = Abs(det1) + Abs(det2); + Icoor2 detMin = Min(det1,det2); + + int OnSwap = 0; + // si 2 triangle infini (bord) => detT = -2; + if (sa == 0) {// les deux triangles sont frontieres + det2=bamg::det(s2->i,sb->i,s1->i); + OnSwap = det2 >0;} + else if (sb == 0) { // les deux triangles sont frontieres + det1=bamg::det(s1->i,sa->i,s2->i); + OnSwap = det1 >0;} + else if(( s1 != 0) && (s2 != 0) ) { + det1 = bamg::det(s1->i,sa->i,s2->i); + det2 = detT - det1; + OnSwap = (Abs(det1) + Abs(det2)) < detA; + + Icoor2 detMinNew=Min(det1,det2); + // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test + if (! OnSwap &&(detMinNew>0)) { + OnSwap = detMin ==0; + if (! OnSwap) { + int kopt = koption; + while (1) + if(kopt) { + // critere de Delaunay pure isotrope + register Icoor2 xb1 = sb->i.x - s1->i.x, + x21 = s2->i.x - s1->i.x, + yb1 = sb->i.y - s1->i.y, + y21 = s2->i.y - s1->i.y, + xba = sb->i.x - sa->i.x, + x2a = s2->i.x - sa->i.x, + yba = sb->i.y - sa->i.y, + y2a = s2->i.y - sa->i.y; + register double + cosb12 = double(xb1*x21 + yb1*y21), + cosba2 = double(xba*x2a + yba*y2a) , + sinb12 = double(det2), + sinba2 = double(t2->det); + + + // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) + OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); + // if(CurrentTh) + // cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ; + // cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12 + // << " Onswap = " << OnSwap << endl; + break; + } + else + { + // critere de Delaunay anisotrope + Real8 som; + I2 AB=(I2) *sb - (I2) *sa; + I2 MAB2=((I2) *sb + (I2) *sa); + R2 MAB(MAB2.x*0.5,MAB2.y*0.5); + I2 A1=(I2) *s1 - (I2) *sa; + I2 D = (I2) * s1 - (I2) * sb ; + R2 S2(s2->i.x,s2->i.y); + R2 S1(s1->i.x,s1->i.y); + { + Metric M=s1->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if (Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som = M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + + } + { + Metric M=s2->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if(Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + OnSwap = som < 2; + break; + } + + } // OnSwap + } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) + } +#ifdef DEBUG1 + if (OnSwap && ( munswap1 || munswap2)) { + cout << " erreur Mark unswap T " << CurrentTh->Number(t1) << " " << CurrentTh->Number(t2) << endl + << *t1 << endl + << *t2 << endl; + return 0; + } +#endif + if( OnSwap ) + bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); + else { + NbUnSwap ++; + t1->SetMarkUnSwap(a1); + } + return OnSwap; + } + + Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega) + { +#ifdef DEBUG + register Int4 NbSwap =0; +#endif + 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++; + +#ifdef DEBUG + assert( s == & (*tria)[VerticesOfTriangularEdge[j][1]] ); + assert( tria->det >0); +#endif + 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; + } + + + + + void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3) + { + // ------------------------------------------- + // s2 + // ! + // /|\ ! + // / | \ ! + // / | \ ! + // tt1 / | \ tt0 ! + // / |s \ ! + // / . \ ! + // / . ` \ ! + // / . ` \ ! + // ---------------- ! + // s0 tt2 s1 + //-------------------------------------------- + + Triangle * tt[3]; // the 3 new Triangles + Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2]; + Icoor2 det3local[3]; + int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; + // infv = ordre of the infini vertex (null) + register int nbd0 =0; // number of zero det3 + register int izerodet=-1,iedge; // izerodet = egde contening the vertex s + Icoor2 detOld = t->det; + + if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ) + { + cerr << " infv " << infv << " det = " << detOld << endl; + cerr << Number(s) << " "<< Number(s0) << " " + << Number(s1) << " " << Number(s2) << endl; + MeshError(3); + } + + // if det3 do not exist then constuct det3 + if (!det3) { + det3 = det3local; // alloc + if ( infv<0 ) { + det3[0]=bamg::det(s ,s1,s2); + det3[1]=bamg::det(s0,s ,s2); + det3[2]=bamg::det(s0,s1,s );} + else { + // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk + det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; + det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; + det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}} + + + if (!det3[0]) izerodet=0,nbd0++; + if (!det3[1]) izerodet=1,nbd0++; + if (!det3[2]) izerodet=2,nbd0++; + + if (nbd0 >0 ) // point s on a egde or on a vertex + if (nbd0 == 1) { + iedge = OppositeEdge[izerodet]; + TriangleAdjacent ta = t->Adj(iedge); + +#ifdef DEBUG1 + cout << " the point " << Number(s) << " is the edge " << izerodet + << " of " << Number(t) << " det3 = " + << det3[0] << " " << det3[1] << " " << det3[2] << " " << endl; + cout << " ta = " << ta << "ta->det =" << ((Triangle*) ta)->det + << " "<< t->det<< endl; +#endif + + // the point is on the edge + // if the point is one the boundary + // add the point in outside part + if ( t->det >=0) { // inside triangle + if ((( Triangle *) ta)->det < 0 ) { + // add in outside triangle + Add(s,( Triangle *) ta); + return;} + }} + else { + cerr << " bug " << nbd0 <SetUnMarkUnSwap(0); + t->SetUnMarkUnSwap(1); + t->SetUnMarkUnSwap(2); + + tt[0]= t; + tt[1]= &triangles[nbt++]; + tt[2]= &triangles[nbt++]; + + if (nbt>nbtx) { + cerr << " No enougth triangles " << endl; + MeshError(999,this); + } + + *tt[1]= *tt[2]= *t; + // gestion of the link + tt[0]->link=tt[1]; + tt[1]->link=tt[2]; + + (* tt[0])(OppositeVertex[0])=&s; + (* tt[1])(OppositeVertex[1])=&s; + (* tt[2])(OppositeVertex[2])=&s; + + tt[0]->det=det3[0]; + tt[1]->det=det3[1]; + tt[2]->det=det3[2]; + + // update adj des triangles externe + tt[0]->SetAdjAdj(0); + tt[1]->SetAdjAdj(1); + tt[2]->SetAdjAdj(2); + // update des adj des 3 triangle interne + const int i0 = 0; + const int i1= NextEdge[i0]; + const int i2 = PreviousEdge[i0]; + + tt[i0]->SetAdj2(i2,tt[i2],i0); + tt[i1]->SetAdj2(i0,tt[i0],i1); + tt[i2]->SetAdj2(i1,tt[i1],i2); + + tt[0]->SetTriangleContainingTheVertex(); + tt[1]->SetTriangleContainingTheVertex(); + tt[2]->SetTriangleContainingTheVertex(); + + + // swap if the point s is on a edge + if(izerodet>=0) { + // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; + int rswap =tt[izerodet]->swap(iedge); + + if (!rswap) + { + cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; +#ifdef DRAWING + if( CurrentTh && withrgraphique) + { + reffecran(); + + DrawMark(s.r); + CurrentTh->inquire(); + DrawMark(s.r); + rattente(1); + } +#endif + } + assert(rswap); + } + +#ifdef DEBUG + tt[0]->check(); + tt[1]->check(); + tt[2]->check(); +#endif +#ifdef DRAWING1 + tt[0]->Draw(); + tt[1]->Draw(); + tt[2]->Draw(); +#endif + + } + + + Int4 Triangles::SplitInternalEdgeWithBorderVertices() + { + Int4 NbSplitEdge=0; + SetVertexFieldOn(); + Int4 it; + Int4 nbvold=nbv; + long int verbosity=2; + for (it=0;it7) + cout <<" Internal edge with two vertices on boundary" + << Number(v0) << " " << Number(v1) << " by " << endl; + } + } + } + } + ReMakeTriangleContainingTheVertex(); + if (nbvold!=nbv) + { + Int4 iv = nbvold; + Int4 NbSwap = 0; + Icoor2 dete[3]; + for (Int4 i=nbvold;iToClose(vi,seuil,hi,hj)) { + // a good new point + vi.ReferenceNumber=0; + vi.DirOfSearch =NoDirOfSearch; + // cout << " Add " << Number(vi) << " " << vi + // << " " << Number(vi) << " <--> " << Number(vi) <link) { + cout << i << " PB insert point " << Number(vi) << vi << Number(vi) + << " tcvi = " << tcvi << " " << tcvi->link << endl; + cout << (*tcvi)[1] << (*tcvi)[2] << endl; + tcvi = FindTriangleContening(vi.i,dete); + cout << (*tcvi)[1] << (*tcvi)[2] << endl; +#ifdef DRAWING1 + inquire(); + penthickness(5); + DrawMark(vi.r); + penthickness(1); + inquire(); +#endif + + MeshError(1001,this); + } + + + quadtree->Add(vi); +#ifdef DRAWING1 + DrawMark(vi.r); +#endif + assert (tcvi && tcvi->det >= 0) ;// internal + Add(vi,tcvi,dete); + NbSwap += vi.Optim(1); + iv++; + // } + } + if (verbosity>3) + { + cout << " Nb Of New Point " << iv ; + cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;} + + nbv = iv; + } + if (NbSplitEdge > nbv-nbvold) + cout << " Warning not enough vertices to split all internal edges " << endl + << " we lost " << NbSplitEdge - ( nbv-nbvold) << " Edges Sorry " << endl; + if (verbosity>2) + cout << "SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << endl; + return NbSplitEdge; + } + Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) + { + long int verbosity=2; + Real8 seuil= 1.414/2 ;// for two close point + Int4 i; + // insertion part --- + + const Int4 nbvnew = nbv-nbvold; + if (verbosity>5) + cout << " Try to Insert the " <ReferenceNumber=i; + } + // be carefull + Int4 iv = nbvold; + for (i=nbvold;iToClose(vi,seuil,hi,hj)) + { + // a good new point + Vertex & vj = vertices[iv]; + Int4 j = vj.ReferenceNumber; + assert( &vj== ordre[j]); + if(i!=j) + { // for valgring + Exchange(vi,vj); + Exchange(ordre[j],ordre[i]); + } + vj.ReferenceNumber=0; + // cout << " Add " << Number(vj) << " " << vj + // << " " << Number(vi) << " <--> " << Number(vj) <link) + { + cerr << i << " PB insert point " << Number(vj) << vj << Number(vi) + << " tcvj = " << tcvj << " " << tcvj->link << endl; + cerr << (*tcvj)[1] << (*tcvj)[2] << endl; + tcvj = FindTriangleContening(vj.i,dete); + cout << (*tcvj)[1] << (*tcvj)[2] << endl; +#ifdef DRAWING1 + inquire(); + penthickness(5); + DrawMark(vj.r); + penthickness(1); + + inquire(); +#endif + + MeshError(1001,this); + } + + + quadtree->Add(vj); +#ifdef DRAWING1 + DrawMark(vj.r); +#endif + assert (tcvj && tcvj->det >= 0) ;// internal + Add(vj,tcvj,dete); + NbSwap += vj.Optim(1); + iv++; + } + } + if (verbosity>3) { + cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ; + cout << " Nb swap = " << NbSwap << " after " ;} + + nbv = iv; + } + +#ifdef DRAWING1 + inquire(); +#endif + + for (i=nbvold;i3) + cout << " NbSwap = " << NbSwap << endl; + + + NbTSwap += NbSwap ; +#ifdef DEBUG + { + Int4 NbErr=0; + Int4 i; + for (i=0;i2) + cout << " -- Triangles::NewPoints "; + if (verbosity>3)cout << " nbv (in) on Boundary = " << nbv <2) + cout << " (Nb of Points from background mesh = " + << nbv-nbvold << " / " << nbv1-nbvold << ")" << endl; + } + else + Bh.ReMakeTriangleContainingTheVertex(); + + Triangle *t; + // generation of the list of next Triangle + // at 1 time we test all the triangles + Int4 Headt =0,next_t; + for(i=0;i= nbt + // the list of test triangle is + // the next traingle on i is -first_np_or_next_t[i] + int iter=0; + // Big loop + do { + iter++; + nbtold = nbt; + nbvold = nbv; +#ifdef DRAWING1 + inquire(); +#endif + + // default size of IntersectionTriangle + + i=Headt; + next_t=-first_np_or_next_t[i]; + for(t=&triangles[i];it, EdgesVertexTriangle[s->vint][1]); + Triangle * tbegin= (Triangle*) ta; + Int4 kt; + do { + kt = Number((Triangle*) ta); + if (first_np_or_next_t[kt]>0) + first_np_or_next_t[kt]=-Headt,Headt=kt; + assert( ta.EdgeVertex(0) == s); + ta = Next(Adj(ta)); + } while ( (tbegin != (Triangle*) ta)); + } + + } while (nbv!=nbvold); + + delete [] first_np_or_next_t; + + Int4 NbSwapf =0,NbSwp; + + // bofbof + + + NbSwp = NbSwapf; + for (i=0;i3) cout << " " ; + if (verbosity>2) + cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT + << " NbSwap final = " << NbSwapf << " Nb Total Of Swap = " << NbTSwap << endl; + + + } + + void Triangles::NewPointsOld(Triangles & Bh) + { // Triangles::NewPointsOld + long int verbosity=2; + Real8 seuil= 0.7 ;// for two neart point + if (verbosity>1) + cout << " begin : Triangles::NewPointsOld " << endl; + Int4 i,k; + int j; + Int4 BeginNewPoint[3]; + Int4 EndNewPoint[3]; +#ifdef TRACETRIANGLE + Int4 trace=0; +#endif + int step[3]; + Int4 *first_np_or_next_t = new Int4[nbtx]; + Int4 ColorEdge[3]; + Int4 color=-1; + Triangle *t; + // generation of the list of next Triangle + // at 1 time we test all the triangles + Int4 Headt =0,next_t; + for(i=0;i= nbt + // the list of test triangle is + // the next Triangle on i is -first_np_or_next_t[i] + Int4 nbtold,nbvold; + + // Big loop + do { + nbtold = nbt; + nbvold = nbv; +#ifdef DRAWING1 + inquire(); +#endif + + // default size of IntersectionTriangle + + i=Headt; + next_t=-first_np_or_next_t[i]; + for(t=&triangles[i];iLocked(j) ; + } +#endif + if (!t->link) continue;// boundary + if (t->det <0) continue; + if (t->Locked(j)) continue; + + TriangleAdjacent tadjj = t->Adj(j); + Triangle * ta= tadjj; + + if (ta->det <0) continue; + + R2 A = vA; + R2 B = vB; + + k=Number(ta); + // the 2 opposite vertices + const Vertex & vC1 = *tj.OppositeVertex(); + const Vertex & vC2 = *tadjj.OppositeVertex(); + +#ifdef TRACETRIANGLE + trace = trace || k == TRACETRIANGLE; + if(trace) { + cout << "Test Arete " << i << " AB = " << A << B + << "i " <link << " ta=" << Number( ta) + << " det " <det ; + cout << " hA = " <Locked(j) << endl; + } +#endif + + if(first_np_or_next_t[k]>0) { // this edge is done before + // find the color of the edge and begin , end of newpoint + register int kk = t->NuEdgeTriangleAdj(j); + assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1])); + assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0])); + register Int4 kolor =3*k + kk; + ColorEdge[j]=kolor; + register Int4 kkk= 1; + step[j]=-1;// other sens + BeginNewPoint[j]=0; + EndNewPoint[j]=-1; // empty list + for (Int4 iv=first_np_or_next_t[k];iv kolor) + break; // the color is passed + else if (vertices[iv].color == kolor) { + EndNewPoint[j]=iv; + if (kkk) // one time test + kkk=0,BeginNewPoint[j]=iv;} + continue; // next edge of the triangle + } // end if( k < i) + + +#ifdef DRAWING1 + penthickness(2); Move(A);Line(B); penthickness(1); +#endif + const Int4 NbvOld = nbv; + lIntTria.SplitEdge(Bh,A,B); + // Int4 NbvNp = + lIntTria.NewPoints(vertices,nbv,nbvx); + Int4 nbvNew=nbv; + nbv = NbvOld; + for (Int4 iv=NbvOld;iv seuil) + && ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) ) + nbv++; + } + + EndNewPoint[j] = nbv-1; + } // end loop for each edge + +#ifdef TRACETRIANGLE + if(trace) { + // verification des point cree + cout << "\n ------------ " << t->link << " " << t->det + << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1] + << " " << BeginNewPoint[2] << " " + << " e " << EndNewPoint[0] << " " << EndNewPoint[1] + << " " << EndNewPoint[2] << " " + << " s " << step[0] << " " << step[1] << " " << step[2] << " " + << endl; + } +#endif + + if (!t->link) continue;// boundary + if (t->det<=0) continue;// outside + // continue; + for(int j0=0;j0<3;j0++) + for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++) + { + // find the neart point one the opposite edge + // to compute i1 + Vertex & vi0 = vertices[i0]; + int kstack = 0; + Int4 stack[10]; + // Int4 savRef[10]; + int j1 = j0; + while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge + // computation of the intersection of edge j1 and DOrto + // take the good sens + + if (BeginNewPoint[j1]> EndNewPoint[j1]) + continue; // + else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) { + for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++) + stack[kstack++] = ii1; + continue;} + + + int k0,k1; + if (step[j1]<0) k0=1,k1=0; // reverse + else k0=0,k1=1; + R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]]; + R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]]; + R2 D = V11-V10; + Real8 c0 = vi0.m(D,(R2) vi0); + + Real8 c10 = vi0.m(D,V10); + Real8 c11 = vi0.m(D,V11); + + Real8 s; + //cout << " --i0 = " << i0 << D << V10 << V11 << endl ; + //cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl; + if (( c10 < c0 ) && (c0 < c11)) + s = (c11-c0)/(c11-c10); + else if (( c11 < c0 ) && (c0 < c10)) + s = (c11-c0) /(c11-c10); + else break; + R2 VP = V10*s + V11*(1-s); + int sss = (c11-c10) >0 ? 1 : -1; +#ifdef DRAWING1 + penthickness(2); + Move((R2) vi0); + Line(VP); + penthickness(1); +#endif + // find the 2 point by dichotomie + //cout << " t =" << Number(t) << " c0 " << c0 ; + Int4 ii0 = BeginNewPoint[j1]; + Int4 ii1 = EndNewPoint[j1]; + Real8 ciii=-1,cii0=-1,cii1=-1 ; + if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 ) + stack[kstack++] = ii0;//cout << " add+0 " << ii0; + else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 ) + stack[kstack++] = ii1;//cout << " add+1 " << ii1; + else { + while ((ii1-ii0)> 1) { + Int4 iii = (ii0+ii1)/2; + ciii = vi0.m( D ,(R2) vertices[iii]); + //cout << " (iii = " << iii << " " << ciii << ") "; + if ( sss * (ciii - c0) <0 ) ii0 = iii; + else ii1 = iii;} + stack[kstack++] = ii0;// cout << " add0 " << ii0; + if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1; + } +#ifdef DEBUG2 + cout << "ii1 = " << ii1 + << " ii0 = " << ii0 << endl; + cout << " cccc = " << cii0 << " " << ciii + << " " << cii1 << " sss=" << sss << endl; +#endif + if (kstack >5) // bug ? + cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl; + } + + stack[kstack++] = -1; // to stop + Int4 i1; + kstack =0; + while( (i1=stack[kstack++]) >= 0) + { // the two parameter is i0 and i1 + assert(i1 < nbv && i1 >= 0); + assert(i0 < nbv && i0 >= 0); + assert(i1 != i0); + R2 v01 = (R2) vertices[i1]- (R2) vertices[i0]; + Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01)); + + +#ifdef DRAWING1 + Move(vertices[i0].r); + Line(vertices[i1].r); +#endif +#ifdef TRACETRIANGLE + if(trace) { + cout << "\n test j" << j <<" " << i0 + << " " << i1 << " d01=" << d01 <= nbvold); + assert (i1 >= nbvold); + assert(i0 != i1); + if (d01 == 0) + break; + if ( d01 < seuil) + if (i1=0) {// good points + // cout <<" i = " << i ; + for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp) + vertices[ip].ReferenceNumber = -1;// mark remove + vertices[ip].ReferenceNumber = -1;// mark remove + // cout << i << " ---> " << kkk << endl; + vertices[kkk] = vertices[i]; + vertices[kkk].i = toI2(vertices[kkk].r); +#ifdef DRAWING1 + DrawMark(vertices[kkk]); +#endif + vertices[kkk++].ReferenceNumber = 0; + + } +#ifdef DRAWING1 + penthickness(1); +#endif + + // insertion part --- + + const Int4 nbvnew = kkk-nbvold; + + cout << " Remove " << nbv - kkk << " to close vertex " ; + cout << " and Insert the " <i,dete); + // Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y); + // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is " + // << Number(nv) << nv->i << endl; + // Int4 kt = Number(tcvi); + // + + quadtree->Add(*vi); // +#ifdef DRAWING1 + DrawMark(vi->r); +#endif + assert (tcvi->det >= 0) ;// internal + Add(*vi,tcvi,dete),NbSwap += vi->Optim(1); + } + } + cout << " Nb swap = " << NbSwap << " after " ; +#ifdef DRAWING1 + inquire(); +#endif + + for (i=nbvold;it, EdgesVertexTriangle[s->vint][1]); + Triangle * tbegin= (Triangle*) ta; + Int4 kt; + do { + kt = Number((Triangle*) ta); + if (first_np_or_next_t[kt]>0) + first_np_or_next_t[kt]=-Headt,Headt=kt; + assert( ta.EdgeVertex(0) == s); + ta = Next(Adj(ta)); + } while ( (tbegin != (Triangle*) ta)); + } + + + } while (nbv!=nbvold); + delete [] first_np_or_next_t; +#ifdef DEBUG + int nberr=0; + for (int it=0;iti, ordre[1]->i, ordre[i]->i ) == 0;) + if ( ++i >= nbv) { + cerr << " All the vertices are aline " << endl; + MeshError(998,this); } + + // echange i et 2 dans ordre afin + // que les 3 premiers ne soit pas aligne + Exchange( ordre[2], ordre[i]); + + // on ajoute un point a l'infini pour construire le maillage + // afin d'avoir une definition simple des aretes frontieres + nbt = 2; + + + // on construit un maillage trivale forme + // d'une arete et de 2 triangles + // construit avec le 2 aretes orientes et + Vertex * v0=ordre[0], *v1=ordre[1]; + + triangles[0](0) = 0; // sommet pour infini + triangles[0](1) = v0; + triangles[0](2) = v1; + + triangles[1](0) = 0;// sommet pour infini + triangles[1](2) = v0; + triangles[1](1) = v1; + const int e0 = OppositeEdge[0]; + const int e1 = NextEdge[e0]; + const int e2 = PreviousEdge[e0]; + triangles[0].SetAdj2(e0, &triangles[1] ,e0); + triangles[0].SetAdj2(e1, &triangles[1] ,e2); + triangles[0].SetAdj2(e2, &triangles[1] ,e1); + + triangles[0].det = -1; // faux triangles + triangles[1].det = -1; // faux triangles + + triangles[0].SetTriangleContainingTheVertex(); + triangles[1].SetTriangleContainingTheVertex(); + + triangles[0].link=&triangles[1]; + triangles[1].link=&triangles[0]; + +#ifdef DEBUG + triangles[0].check(); + triangles[1].check(); +#endif + // nbtf = 2; + if ( !quadtree ) quadtree = new QuadTree(this,0); + quadtree->Add(*v0); + quadtree->Add(*v1); + + // on ajoute les sommets un Ò un + Int4 NbSwap=0; + + time1=CPUtime(); + + if (verbosity>3) cout << " -- Begin of insertion process " << endl; + + for (Int4 icount=2; icounti,dete); + quadtree->Add(*vi); + Add(*vi,tcvi,dete); + NbSwap += vi->Optim(1,0); +#ifdef DRAWING1 + inquire(); +#endif + + }// fin de boucle en icount + time2=CPUtime(); + if (verbosity>3) + cout << " NbSwap of insertion " << NbSwap + << " NbSwap/Nbv " << (float) NbSwap / (float) nbv + << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " + << (float)NbUnSwap /(float) nbv + <Optim(0,0); + timeloop = CPUtime(); + if (verbosity>3) + cout << " Optim Loop "<4) + cout << " init " << time1 - time0 << " initialisation, " + << time2 - time1 << "s, insert point " + << time3 -time2 << "s, optim " << endl + << " Init Total Cpu Time = " << time3 - time0 << "s " << endl; + +#ifdef DRAWING1 + inquire(); +#endif + CurrentTh=OldCurrentTh; + } + + + void Triangles::ForceBoundary() + { + long int verbosity=2; + if (verbosity > 2) + cout << " -- ForceBoundary nb of edge " << nbe << endl; + int k=0; + Int4 nbfe=0,nbswp=0,Nbswap=0; + for (Int4 t = 0; t < nbt; t++) + if (!triangles[t].det) + k++,cerr << " det T" << t << " = " << 0 << endl; + if (k!=0) { + cerr << " ther is " << k << " triangles of mes = 0 " << endl; + MeshError(11,this);} + + TriangleAdjacent ta(0,0); + for (Int4 i = 0; i < nbe; i++) + { + nbswp = ForceEdge(edges[i][0],edges[i][1],ta); + + if ( nbswp < 0) k++; + else Nbswap += nbswp; + if (nbswp) nbfe++; + if ( nbswp < 0 && k < 5) + { + cerr << " Missing Edge " << i << " v0 = " << Number(edges[i][0]) << edges[i][0].r + <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ; + if(ta.t) + { + Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1); + cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n"; + } + else cerr << endl; + + } + if ( nbswp >=0 && edges[i].on->Cracked()) + ta.SetCracked(); + } + + + if (k!=0) { + cerr << " they is " << k << " lost edges " << endl; + cerr << " The boundary is crossing may be!" << endl; + MeshError(10,this); + } + for (Int4 j=0;j 3) + cout << " Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl; + + } + + void Triangles::FindSubDomain(int OutSide=0) + { + long int verbosity=2; + //#define DRAWING1 + + if (verbosity >2) + { + if (OutSide) + cout << " -- Find all external sub-domain "; + else + cout << " -- Find all internal sub-domain "; + if(verbosity>99) + { + + for(int i=0;i 4) cout << " OutSide=" << OutSide << endl; + short * HeapArete = new short[nbt]; + Triangle ** HeapTriangle = new Triangle* [nbt]; + Triangle *t,*t1; + Int4 k,it; + + for (Int4 itt=0;ittlink = t ; // sd forme d'un triangle cicular link +#ifdef DRAWING1 + t->Draw(NbSubDomTot-1); +#endif + + + HeapTriangle[i] =t ; + HeapArete[i] = 3; + + while (i >= 0) // boucle sur la pile + { while ( HeapArete[i]--) // boucle sur les 3 aretes + { + int na = HeapArete[i]; + Triangle * tc = HeapTriangle[i]; // triangle courant + if( ! tc->Locked(na)) // arete non frontiere + { + Triangle * ta = tc->TriangleAdj(na) ; // næ triangle adjacent + if (ta->link == 0 ) // non deja chainer => on enpile + { + i++; +#ifdef DRAWING1 + ta->Draw(NbSubDomTot-1); +#endif + ta->link = t->link ; // on chaine les triangles + t->link = ta ; // d'un meme sous domaine + HeapArete[i] = 3; // pour les 3 triangles adjacents + HeapTriangle[i] = ta; + }} + } // deplie fin de boucle sur les 3 adjacences + i--; + } + } + } + + // supression de tous les sous domaine infini <=> contient le sommet NULL + it =0; + NbOutT = 0; + while (itlink; + t1->link=0;}//while (t) + } + } + it++;} // end while (it All triangles are outside " << endl; + MeshError(888,this); + } + + delete [] HeapArete; + delete [] HeapTriangle; + + + if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains ) + { // No geom sub domain + Int4 i; + if (subdomains) delete [] subdomains; + subdomains = new SubDomain[ NbSubDomTot]; + NbSubDomains= NbSubDomTot; + for ( i=0;ilink; + mark[it]=k; +#ifdef DRAWING1 + t1->Draw(k); +#endif + subdomains[k].head = t1; + // cout << " New -- " << Number(t1) << " " << it << endl; + do {// cout << " k " << k << " " << Number(t) << endl; + mark[Number(t)]=k; +#ifdef DRAWING1 + t->Draw(k); +#endif + t=t->link; + } while (t!=t1); +#ifdef DRAWING1 + t1->Draw(k); +#endif + mark[it]=k++;} + // else if(mark[it] == -2 ) triangles[it].Draw(999); + it++;} // end white (it=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0) + nbk--,subdomains[kr].ref=subdomains[kl].ref-1; + if (kr >=0 && subdomains[kr].ref <0 && kl >=0 && subdomains[kl].ref>=0) + nbk--,subdomains[kl].ref=subdomains[kr].ref-1; + if(kr<0 && kl >=0 && subdomains[kl].ref>=0) + nbk--,subdomains[kl].ref=-1; + if(kl<0 && kr >=0 && subdomains[kr].ref>=0) + nbk--,subdomains[kr].ref=-1; + // cout << " after \t " + // << kl << subdomains[kl].ref << " rr " << kr + // << subdomains[kr].ref << endl; + } + } + // cout << subdomains[0].ref << subdomains[1].ref << endl; + Int4 j=0; + for ( i=0;ilink; + t1->link=0;}//while (t) + } + + if(verbosity>4) + cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl; + NbSubDomains=j; + } + + delete [] mark; + + } + else + { // find the head for all sub domaine + if (Gh.NbSubDomains != NbSubDomains && subdomains) + delete [] subdomains, subdomains=0; + if (! subdomains ) + subdomains = new SubDomain[ Gh.NbSubDomains]; + NbSubDomains =Gh.NbSubDomains; + if(verbosity>4) + cout << " find the " << NbSubDomains << " sub domain " << endl; + Int4 err=0; + ReMakeTriangleContainingTheVertex(); + Int4 * mark = new Int4[nbt]; + Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge(); + + for (it=0;itt; + int sens = Gh.subdomains[i].sens; + // test if ge and e is in the same sens + // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref + // << " ref edge =" << eg.ref << " sens " << sens ; + if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) + sens = -sens ; + subdomains[i].sens = sens; + subdomains[i].edge = &e; + // cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl; + // cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl; + assert(t && sens); + + TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges + + while (1) + { + assert( v0 == ta.EdgeVertex(1) ); + // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl; + if (ta.EdgeVertex(0) == v1) { // ok we find the edge + if (sens>0) + subdomains[i].head=t=Adj(ta); + else + subdomains[i].head=t=ta; + //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl; + if(t= triangles+nbt || t->det < 0 + || t->link == 0) // Ajoute aout 200 + { + cerr << " Error in the def of sub domain "<=0) { + if(verbosity>10) + cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref + << " is previouly defined with " <Draw(i); +#endif + mark[Number(tt)]=i; + tt=tt->link; + } while (tt!=t); + if(verbosity>7) + cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl; + break;} + ta = Previous(Adj(ta)); + if(t == (Triangle *) ta) { + err++; + cerr << " Error in the def of sub domain " << i + << " edge=" << Gh.Number(eg) << " " << sens << endl; + break;} + // cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl; + + } + + } + if (err) MeshError(777,this); + + if (inew < NbSubDomains) { + if (verbosity>5) + cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl; + NbSubDomains=inew;} + + + for (it=0;it 4) + cout << " " ; + if (verbosity> 2) + cout << " Nb of Sub borned Domain = " << NbSubDomTot << " NbOutTriangles = " << NbOutT <=vertices && v < ve) + VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)]; + } + + for (i=0;i< NbVerticesOnGeomEdge;i++) + { + Vertex *v =VerticesOnGeomEdge[i].mv; + if (v>=vertices && v < ve) + VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)]; + } + + for (i=0;i< NbVertexOnBThVertex;i++) + { + Vertex *v=VertexOnBThVertex[i].v; + if (v>=vertices && v < ve) + VertexOnBThVertex[i].v=vertices+renu[Number(v)]; + } + + for (i=0;i< NbVertexOnBThEdge;i++) + { + Vertex *v=VertexOnBThEdge[i].v; + if (v>=vertices && v < ve) + VertexOnBThEdge[i].v=vertices+renu[Number(v)]; + } + + // move the Vertices without a copy of the array + // be carefull not trivial code + Int4 j; + for ( it=0;it= 0) // a new sub cycle + { + i=it; + Vertex ti=vertices[i],tj; + while ( (j=renu[i]) >= 0) + { // i is old, and j is new + renu[i] = -1-renu[i]; // mark + tj = vertices[j]; // save new + vertices[j]= ti; // new <- old + i=j; // next + ti = tj; + } + } + if (quadtree) + { delete quadtree; + quadtree = new QuadTree(this); + } + for ( it=0;it=0 && kt < nbt ); + assert(renu[kt]==-1); + renu[kt]=k++; + } + while (t0 != (t=t->link)); + } + if (verbosity>9) + cout << " number of inside triangles " << k << " nbt = " << nbt << endl; + // take is same numbering if possible + if(justcompress) + for ( k=0,it=0;it=0 ) + renu[it]=k++; + + // put the outside triangles at the end + for ( it=0;it= 0) // a new sub cycle + { + i=it; + Triangle ti=triangles[i],tj; + while ( (j=renu[i]) >= 0) + { // i is old, and j is new + renu[i] = -1; // mark + tj = triangles[j]; // save new + triangles[j]= ti; // new <- old + i=j; // next + ti = tj; + } + } + delete [] renu; + nt = nbt - NbOutT; + +#ifdef DEBUG + // verif + for ( it=0;it=0 &&num < nbt); + reft[num]=i; + // cout << Number(t0) << " " <link)); + } + // NbOutT = nbt - k; + if (verbosity>5) + cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k + << " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl; + + return k; + + } + + /* + void Triangles::ConsLinkTriangle() + { + for (Int4 i=0;ilink) + { + register Int4 color = t->color-1; + assert(color=0); + if (hst=subdomains[color].head) { + t->link=hst->link; + hst->link=t; } + else { + subdomains[color].head = t; + t->link=t;}// circular link + } + } + { + for (Int4 i=0;iNearestVertex(i,j); } + + + void Triangles::PreInit(Int4 inbvx,char *fname) + { + long int verbosity=0; + + srand(19999999); + OnDisk =0; + NbRef=0; + // allocGeometry=0; + identity=0; + NbOfTriangleSearchFind =0; + NbOfSwapTriangle =0; + nbiv=0; + nbv=0; + nbvx=inbvx; + nbt=0; + NbOfQuad = 0; + nbtx=2*inbvx-2; + NbSubDomains=0; + NbVertexOnBThVertex=0; + NbVertexOnBThEdge=0; + VertexOnBThVertex=0; + VertexOnBThEdge=0; + + NbCrackedVertices=0; + NbCrackedEdges =0; + CrackedEdges =0; + nbe = 0; + name = fname ; + + if (inbvx) { + vertices=new Vertex[nbvx]; + assert(vertices); + ordre=new (Vertex* [nbvx]); + assert(ordre); + triangles=new Triangle[nbtx]; + assert(triangles);} + else { + vertices=0; + ordre=0; + triangles=0; + nbtx=0; + } + if ( name || inbvx) + { + time_t timer =time(0); + char buf[70]; + strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer)); + counter++; + char countbuf[30]; + sprintf(countbuf,"%d",counter); + int lg =0 ; + if (&BTh != this && BTh.name) + lg = strlen(BTh.name)+4; + identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)]; + identity[0]=0; + if (lg) + strcat(strcat(strcat(identity,"B="),BTh.name),", "); + + if (Gh.name) + strcat(strcat(identity,"G="),Gh.name); + strcat(strcat(identity,";"),countbuf); + strcat(identity,buf); + // cout << "New MAILLAGE "<< identity << endl; + } + + quadtree=0; + // edgescomponante=0; + edges=0; + VerticesOnGeomVertex=0; + VerticesOnGeomEdge=0; + NbVerticesOnGeomVertex=0; + NbVerticesOnGeomEdge=0; + // nbMaxIntersectionTriangles=0; + // lIntTria; + subdomains=0; + NbSubDomains=0; + // Meshbegin = vertices; + // Meshend = vertices + nbvx; + if (verbosity>98) + cout << "Triangles::PreInit() " << nbvx << " " << nbtx + << " " << vertices + << " " << ordre << " " << triangles <= nbvx) + { + cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; + MeshError(1,this); + } + assert(vertices); + for (i=0;i Th + VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); + // cout << "--------- " <to);// use of Geom -> Th + VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv); + gv->to->m = bv->m; // for taking the metrix of the background mesh + ;} + } + assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex); + // new stuff FH with curve + // find the begin of the curve in BTh + { + Gh.UnMarkEdges(); + int bfind=0; + /* + cout << " nb curves = " << Gh.NbOfCurves << endl; + for(int i=0;iMark() && ei[jedge].on->IsRequiredVertex() ) + { + */ + // new code FH 2004 + Real8 L=0; + for (int icurve=0;icurveMark() can be change in + // loop for(jedge=0;jedge<2;jedge++) + // new curve + // good the find a starting edge + Real8 Lstep=0,Lcurve=0;// step between two points (phase==1) + Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) + + + // cout.precision(16); + for(int phase=0;phase<=step;phase++) + { + + for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) + { + + int icurveequi= Gh.Number(curve); + + if( phase == 0 && icurveequi != icurve) continue; + + int k0=jedge,k1; + Edge * pe= BTh.edges+iedge; + //GeometricalEdge *ong = ei.on; + int iedgeequi=bcurve[icurveequi]/2; + int jedgeequi=bcurve[icurveequi]%2; + + int k0equi=jedgeequi,k1equi; + Edge * peequi= BTh.edges+iedgeequi; + GeometricalEdge *ongequi = peequi->on; + + Real8 sNew=Lstep;// abcisse of the new points (phase==1) + L=0;// length of the curve + Int4 i=0;// index of new points on the curve + register GeometricalVertex * GA0 = *(*peequi)[k0equi].on; + Vertex *A0; + A0 = GA0->to; // the vertex in new mesh + Vertex *A1; + VertexOnGeom *GA1; + Edge * PreviousNewEdge = 0; + // cout << " --------------New Curve phase " << phase + // << "---------- A0=" << *A0 << ei[k0] <=0 && A0-vertices Required() ) + { + GeometricalVertex *GA1 = *(*peequi)[1-k0equi].on; + A1 = GA1->to; // + } + else + for(;;) + { + // assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe); + Edge &ee=*pe; + Edge &eeequi=*peequi; + k1 = 1-k0; // next vertex of the edge + k1equi= 1 - k0equi; + + assert(pe && ee.on); + ee.on->SetMark(); + Vertex & v0=ee[0], & v1=ee[1]; + R2 AB= (R2) v1 - (R2) v0; + Real8 L0=L,LAB; + LAB = LengthInterpole(v0.m,v1.m,AB); + L+= LAB; + if (phase) {// computation of the new points + while ((i!=NbCreatePointOnCurve) && sNew <= L) { + // cout << " L0= " << L0 << " L " << L << " sN=" + // << sNew << " LAB=" << LAB << " NBPC =" <= L0); + assert(LAB); + + + assert(vertices && nbv=0 && se < 1.000000001); +#ifdef DEBUG + se = abscisseInterpole(v0.m,v1.m,AB,se); // because code \ref(xxx) +#else + se = abscisseInterpole(v0.m,v1.m,AB,se,1); +#endif + assert(se>=0 && se <= 1); + //((k1==1) != (k1==k1equi)) + se = k1 ? se : 1. - se; + se = k1==k1equi ? se : 1. - se; + VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save + ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); + A1->ReferenceNumber = eeequi.ref; + A1->DirOfSearch =NoDirOfSearch; + //cout << icurveequi << " " << i << " " << *A1 << endl; + e->on = ongequi; + e->v[0]= A0; + e->v[1]= A1; + if(verbosity>99) + cout << i << "+ New P "<< nbv-1 << " " <r <r - A0->r; + Real8 dp = LengthInterpole(A0->m,A1->m,A1A0); + if (dp > 1.4) { + cerr << " PB new Points "<< nbv-1 ; + cerr << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= " ; + cerr << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <CurveNumber==ei.on->CurveNumber); + if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl; + if ( ee[k1].on->IsRequiredVertex()) { + assert(eeequi[k1equi].on->IsRequiredVertex()); + register GeometricalVertex * GA1 = *eeequi[k1equi].on; + A1=GA1->to;// the vertex in new mesh + assert (A1-vertices>=0 && A1-vertices " << icurveequi <<"-----" << + NbCreatePointOnCurve << " == " <on = ongequi; + e->v[0]= A0; + e->v[1]= A1; + e->ref = peequi->ref; + e->adj[0]=PreviousNewEdge; + e->adj[1]=0; + if (PreviousNewEdge) + PreviousNewEdge->adj[1] = e; + PreviousNewEdge = e; + // cout << "Last new edge " << nbe << " " << " on " << Gh.Number(pe->on) + // << " of curve =" <on->CurveNumber <Draw(); + A1->Draw(); + A0->Draw(); + // inquire(); +#endif + assert(i==NbCreatePointOnCurve); + + } + } // end loop on equi curve + + if (!phase) { // + Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg + Lstep = L/NbSegOnCurve; + Lcurve = L; + NbCreatePointOnCurve = NbSegOnCurve-1; + + for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) + { + NbOfNewEdge += NbSegOnCurve; + NbOfNewPoints += NbCreatePointOnCurve; + } + if(verbosity>5) + cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep=" + << Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve < nbvx) + { + cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl; + MeshError(3,this); + } + //cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl; + edges = new Edge[NbOfNewEdge]; + nbex = NbOfNewEdge; + if(NbOfNewPoints) { // + VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints]; + NbVertexOnBThEdge =NbOfNewPoints; + VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints]; + NbVerticesOnGeomEdgex = NbOfNewPoints; } + NbOfNewPoints =0; + NbOfNewEdge = 0; + } + } // for(step;;) + assert(nbe); + + delete [] bcurve; + + +#ifdef DRAWING1 + reffecran(); + InitDraw(); + Draw(); + inquire(); +#endif + + Insert(); + ForceBoundary(); + FindSubDomain(); + +#ifdef DRAWING1 + reffecran(); + Draw(); + inquire(); +#endif + // NewPointsOld(*this) ; + // BTh.ReMakeTriangleContainingTheVertex(); // FH change => put in NewPoints + // for (Int4 iv=0;iv= nbvx) + { + cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; + MeshError(1,this); + } + for (i=0;i Th + VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); + // cout << "--------- " <The(); + AB = b->r - a->r; + Metric MA = background ? BTh.MetricAt(a->r) :a->m ; + Metric MB = background ? BTh.MetricAt(b->r) :b->m ; + Real8 ledge = (MA(AB) + MB(AB))/2; + // + const int MaxSubEdge = 10; + int NbSubEdge = 1; + Real8 lSubEdge[MaxSubEdge]; + R2 A,B; + if (ledge < 1.5) + lSubEdge[0] = ledge; + else { + NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); + A= a->r; + Metric MAs =MA,MBs; + // cout << " lSubEdge old=" << ledge + // << " new " << A << MA << endl; + ledge = 0; + Real8 x =0, xstep= 1. / NbSubEdge; + for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { + x += xstep; + B = e->F(k ? x : 1-x); + MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); + AB = A-B; + lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); + // cout << " " << lSubEdge[kk] << " x " << x + // << " " << A << B << MA << MB<< endl ; + } + // cout << endl; + } + + Real8 lcurveb = lcurve+ ledge ; + while (lcurve<=s && s <= lcurveb && nbv < nbvend) + { + // New points + + // Real8 aa=(lcurveb-s)/ledge; + // Real8 bb=(s-lcurve)/ledge; + + Real8 ss = s-lcurve; + // 1) find the SubEdge containing ss by dichotomie + int kk0=-1,kk1=NbSubEdge-1,kkk; + Real8 ll0=0,ll1=ledge,llk; + while (kk1-kk0>1) + { + if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) + kk1=kkk,ll1=llk; + else + kk0=kkk,ll0=llk;} + assert(kk1 != kk0); + + Real8 sbb = (ss-ll0 )/(ll1-ll0); + Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb; + + // new vertex on edge + vb = &vertices[nbv++]; + vb->m = Metric(aa,a->m,bb,b->m); + vb->ReferenceNumber = e->ref; + vb->DirOfSearch =NoDirOfSearch; + Real8 abcisse = k ? bb : aa; + vb->r = e->F( abcisse ); + VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); + + // to take in account the sens of the edge + + s += lstep; + edges[nbe].v[0]=va; + edges[nbe].v[1]=vb; + edges[nbe].ref = e->ref; + edges[nbe].on = e; + edges[nbe].adj[0] = PreviousNewEdge; + if(PreviousNewEdge) + PreviousNewEdge->adj[1] = &edges[nbe]; +#ifdef DRAWING1 + vb->Draw(); + edges[nbe].Draw(); +#endif + PreviousNewEdge = edges + nbe; + nbe++; +#ifdef DEBUG1 + cout << " new points " << nbv-1 << " " << vb->r ; + cout << " new edge " << nbe-1 << " " ; + cout << va << vb << " kk0 = " << kk0 + << " " << kk1 << " ss=" << ss ; + cout << " " << sbb << endl; + cout << " " << aa << va->r << bb << vb->r + <<" length=" << Norme(va->r-vb->r) << endl; + cout << " s " << s << " lstep= " << lstep + << " ledge= " << ledge + << " lcurve= " << lcurve << endl; +#endif + va = vb; + } + lcurve = lcurveb; + e->SetMark(); + // cout << e-Gh.edges << ", " << k << " " + // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" " + // << lcurve<< ";; " ; + a=b; + if (b->Required() ) break; + int kprev=k; + k = e->SensAdj[kprev];// next vertices + e = e->Adj[kprev]; + assert(e); + }// for(;;) + vb = b->to; + // cout << endl; + NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1); + NbNewPoints = NbEdgeCurve-1; + if(!kstep) + { NbVerticesOnGeomEdge0 += NbNewPoints; + NbOfCurves++;} + + nbvend=nbv+NbNewPoints; + + lstep = lcurve / NbEdgeCurve; + // cout <<"lstep " << lstep << " lcurve " + // << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " < " ; + Adj(on,j); // next geom edge + j=1-j; + // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)] + // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl; + if (e[Gh.Number(on)]) break; // optimisation + e[Gh.Number(on)] = ei; + } + } + + int kk=0; + for ( i=0;i10) + cout << " ~Triangles "<< this <<" "<< identity << endl; + if(vertices) delete [] vertices; + if(edges) delete [] edges; + if(triangles) delete [] triangles; + if(quadtree) delete quadtree; + if(ordre) delete [] ordre; + if( subdomains) delete [] subdomains; + if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; + if (VerticesOnGeomVertex) delete [] VerticesOnGeomVertex; + if (name) delete [] name; + if (identity) delete [] identity; + if (VertexOnBThVertex) delete [] VertexOnBThVertex; + if (VertexOnBThEdge) delete [] VertexOnBThEdge; + + if (&Gh) + { + if (Gh.NbRef>0) Gh.NbRef--; + else if (Gh.NbRef==0) delete &Gh; + } + if (&BTh && (&BTh != this)) + { + if (BTh.NbRef>0) BTh.NbRef--; + else if (BTh.NbRef==0) delete &BTh; + } + PreInit(0); // set all to zero + + } + + void Triangles::SetIntCoor(const char * strfrom) + { + pmin = vertices[0].r; + pmax = vertices[0].r; + + // recherche des extrema des vertices pmin,pmax + Int4 i; + for (i=0;i0); + + // generation of integer coord + for (i=0;i2) + cout << " -- FillHoleInMesh: Nb of vertices =" << nbv + << " Pmin = "<< pmin << " Pmax = "<< pmax << endl; + + assert(ordre); + for (i=0;iaddtrie(Number(edges[i][0]),Number(edges[i][1]))); + if (kk != nbe) + { + cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl; + MeshError(1002,this); + } + 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 (k5) { + 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 = " << edge4->nb() << endl; + cout << " - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; + } + + + // check the consistant of edge[].adj and the geometrical required vertex + Int4 k=0; + for (i=0;inb();i++) + if (st[i] >=0) // edge alone + if (i < nbe) + { + Int4 i0=edge4->i(i);ordre[i0] = vertices+i0; + Int4 i1=edge4->j(i);ordre[i1] = vertices+i1; + } + else { + k++; + if (verbosity>20 && k <20) + { + Int4 i0=edge4->i(i); + Int4 i1=edge4->j(i); + cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl; + } + } + + if(k != 0) { + if (verbosity>20) + { + cout << " The given edge are " << endl; + for (int i=0;i< nbe;i++) + cout << " Edge " << i << " : " << Number(edges[i][0]) << " " << Number(edges[i][1]) + << " " << edges[i].ref << endl; + } + cerr << k << " boundary edges are not defined as edges " << endl; + MeshError(9998,this); + } + printf("salope1\n"); + + // generation of the mesh with boundary points + Int4 nbvb = 0; + for (i=0;i9) + cout << " Nbtriafillhole triafillhole*" << triafillhole << endl; + triangles = triafillhole; + + nbt=2; + nbtx= Nbtriafillhole; + + for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){ + if ( ++i >= nbvb) { + cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl; + MeshError(998,this); + } + } + Exchange( ordre[2], ordre[i]); + + Vertex * v0=ordre[0], *v1=ordre[1]; + + + triangles[0](0) = 0; // sommet pour infini + triangles[0](1) = v0; + triangles[0](2) = v1; + + triangles[1](0) = 0;// sommet pour infini + triangles[1](2) = v0; + triangles[1](1) = v1; + const int e0 = OppositeEdge[0]; + const int e1 = NextEdge[e0]; + const int e2 = PreviousEdge[e0]; + triangles[0].SetAdj2(e0, &triangles[1] ,e0); + triangles[0].SetAdj2(e1, &triangles[1] ,e2); + triangles[0].SetAdj2(e2, &triangles[1] ,e1); + + triangles[0].det = -1; // faux triangles + triangles[1].det = -1; // faux triangles + + triangles[0].SetTriangleContainingTheVertex(); + triangles[1].SetTriangleContainingTheVertex(); + + triangles[0].link=&triangles[1]; + triangles[1].link=&triangles[0]; + + printf("salope2\n"); + // nbtf = 2; + if ( !quadtree ) delete quadtree; // ->ReInitialise(); + + quadtree = new QuadTree(this,0); + quadtree->Add(*v0); + quadtree->Add(*v1); + + // on ajoute les sommets un a un + Int4 NbSwap=0; + for (Int4 icount=2; icounti,dete); + quadtree->Add(*vi); + Add(*vi,tcvi,dete); + NbSwap += vi->Optim(1,1); + + }// end loop on icount + printf("salope3\n"); + + //Int4 nbtfillhole = nbt; + // inforce the boundary + TriangleAdjacent ta(0,0); + Int4 nbloss = 0,knbe=0; + for ( i = 0; i < nbe; i++) + if (st[i] >=0) // edge alone => on border ... FH oct 2009 + { + Vertex & a=edges[i][0], & b = edges[i][1]; + if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. + { + knbe++; + if (ForceEdge(a,b,ta)<0) + nbloss++; + } + } + if(nbloss) + { + cerr << " we loss some " << nbloss << " " << " edges other " << knbe << endl; + MeshError(1100,this); + } + FindSubDomain(1); + // remove all the hole + // remove all the good sub domain + Int4 krm =0; + for (i=0;iaddtrie(v0?Number(v0):nbv,v1? Number(v1):nbv); + assert(st[k] >=0); + tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); + ta.SetLock(); + st[k]=-2-st[k]; + } + } + } + Int4 NbTfillHoll =0; + for (i=0;i=0) + { + savetriangles[savenbt]=triangles[i]; + savetriangles[savenbt].link=0; + savenbt++; + } + // gestion of the adj + k =0; + Triangle * tmax = triangles + nbt; + for (i=0;i= triangles && ta < tmax) + { + ta= savetriangles + ta->color; + ti.SetAdj2(j,ta,aa); + if(lck) ti.SetLocked(j); + } + } + } + // OutSidesTriangles = triangles; + // Int4 NbOutSidesTriangles = nbt; + + // restore triangles; + nbt=savenbt; + nbtx=savenbtx; + delete [] triangles; + delete [] subdomains; + triangles = savetriangles; + subdomains = savesubdomains; + // cout << triangles << " <> " << OutSidesTriangles << endl; + /* k=0; + for (i=0;iAdd(vertices[i]); + } + printf("salope4c\n"); + + SetVertexFieldOn(); + printf("salope5a\n"); + + for (i=0;iIsRequiredVertex()) { + cerr << " Erreur adj et sommet requis edges [" << i << "][ " << j << "]= " + << Number(edges[i][j]) << " : " << " on = " << Gh.Number(edges[i].on) ; + if (edges[i][j].on->OnGeomVertex()) + cerr << " vertex " << Gh.Number(edges[i][j].on->gv); + else if (edges[i][j].on->OnGeomEdge()) + cerr << "Edges " << Gh.Number(edges[i][j].on->ge); + else + cerr << " = " << edges[i][j].on ; + cerr << endl; + } + } + } + } + } + + printf("salope5\n"); + CurrentTh=OldCurrentTh; + } + + Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx) // COPY OPERATOR + : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) + { + Gh.NbRef++; + nbvxx = Max(nbvxx,Th.nbv); + Int4 i; + // do all the allocation to be sure all the pointer existe + + char * cname = 0; + if (Th.name) + { + cname = new char[strlen(Th.name)+1]; + strcpy(cname,Th.name); + } + PreInit(nbvxx,cname);// to make the allocation + // copy of triangles + nt=Th.nt; + nbv = Th.nbv; + nbt = Th.nbt; + nbiv = Th.nbiv; + nbe = Th.nbe; + NbSubDomains = Th.NbSubDomains; + NbOutT = Th.NbOutT; + NbOfQuad = Th.NbOfQuad ; + NbOfSwapTriangle =0; + NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex; + if(NbVerticesOnGeomVertex) + VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; + NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge; + if (NbVerticesOnGeomEdge) + VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ; + if (& BTh == & Th.BTh) // same back ground + { + BTh.NbRef++; + NbVertexOnBThVertex = Th.NbVertexOnBThVertex; + if(NbVertexOnBThVertex) + VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex]; + NbVertexOnBThEdge = Th.NbVertexOnBThEdge; + if(NbVertexOnBThEdge) + VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge]; + } + else + { // no add on back ground mesh + BTh.NbRef++; + NbVertexOnBThVertex=0; + VertexOnBThVertex=0; + NbVertexOnBThEdge=0; + VertexOnBThEdge=0; + // assert (& BTh == this); // --- a voir + + } + + + if(nbe) + edges = new Edge[nbe]; + if(NbSubDomains) + subdomains = new SubDomain[NbSubDomains]; + pmin = Th.pmin; + pmax = Th.pmax; + coefIcoor = Th.coefIcoor; + for(i=0;iDraw(); +DrawMark( s->r); +#endif +ttc = t->at[j]; +jc = NextEdge[t->aa[j]&3]; +cout << *t << " " << VerticesOfTriangularEdge[j][1] << "\n\t try swap " << * ttc << " " << jc ; +while ( ttc->swap(jc,koption)) { +NbSwap++,assert(k++<20000); +ttc = t->at[j]; +jc = NextEdge[t->aa[j]&3]; +cout << "\n\t s " << *ttc << " " << jc << endl; +} +cout << endl; +t = ttc; +j = NextEdge[jc]; +assert(k<20000); +} while ( (tbegin != t)); + +return NbSwap; +} +*/ +Int4 Triangle::Optim(Int2 i,int koption) +{ + // turne around in positif sens + Int4 NbSwap =0; +#ifdef DEBUG + Vertex * s = ns[i]; +#endif + Triangle *t = this; + int k=0,j =OppositeEdge[i]; + int jp = PreviousEdge[j]; + // initialise tp, jp the previous triangle & edge + Triangle *tp= at[jp]; + jp = aa[jp]&3; +#ifdef DEBUG + assert(tp->at[jp] == this); +#endif + do { +#ifdef DEBUG + assert(k++<20000); + assert( s == & (*t)[OppositeVertex[j]] ); +#endif + // cout << *t << " " << j << "\n\t try swap " ; + while (t->swap(j,koption)) + { + NbSwap++; + assert(k++<20000); + t= tp->at[jp]; // set unchange t qnd j for previous triangles + j= NextEdge[tp->aa[jp]&3]; + // cout << "\n\t s " << *t << " " << j << endl; +#ifdef DEBUG + assert( s == & (*t)[OppositeVertex[j]] ); +#endif + } + // end on this Triangle + tp = t; + jp = NextEdge[j]; + + t= tp->at[jp]; // set unchange t qnd j for previous triangles + j= NextEdge[tp->aa[jp]&3]; + + } while( t != this); + return NbSwap; +} + +void Triangles::SmoothingVertex(int nbiter,Real8 omega ) +{ + long int verbosity=0; + // if quatree exist remove it end reconstruct + if (quadtree) delete quadtree; + quadtree=0; + ReMakeTriangleContainingTheVertex(); + Triangle vide; // a triangle to mark the boundary vertex + Triangle ** tstart= new Triangle* [nbv]; + Int4 i,j,k; + // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide + if ( this == & BTh) + for ( i=0;i2) + cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl; + for (k=0;k3) + cout << " Move max = " << sqrt(delta) << " iteration = " + << k << " Nb of Swap = " << NbSwap << endl; + } + + delete [] tstart; + if (quadtree) quadtree= new QuadTree(this); +} +void Triangles::MakeQuadTree() +{ + long int verbosity=2; + if(verbosity>8) + cout << " MakeQuadTree" << endl; + if ( !quadtree ) quadtree = new QuadTree(this); + + +#ifdef DRAWING1 + quadtree->Draw(); + rattente(1); + reffecran(); + quadtree->Draw(); + rattente(1); +#endif + +} +void Triangles::ShowRegulaty() const// Add FH avril 2007 +{ + const Real8 sqrt32=sqrt(3.)*0.5; + const Real8 aireKh=sqrt32*0.5; + D2 Beq(1,0),Heq(0.5,sqrt32); + D2xD2 Br(D2xD2(Beq,Heq).t()); + D2xD2 B1r(Br.inv()); + /* D2xD2 BB = Br.t()*Br; + cout << " BB = " << BB << " " << Br*B1r << endl; + MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y); + MatVVP2x2 VMM(MMM); + cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl; + */ + double gammamn=1e100,hmin=1e100; + double gammamx=0,hmax=0; + double beta=1e100; + double beta0=0; + double alpha2=0; + double area=0,Marea=0; + // Real8 cf= Real8(coefIcoor); + // Real8 cf2= 6.*cf*cf; + int nt=0; + for (int it=0;itlink || Number(ta) >= it) + { + Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]]; + Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]]; + if ( !& vP || !&vQ) continue; + R2 PQ = vQ.r - vP.r; + Real8 l = log(LengthInterpole(vP,vQ,PQ)); +#ifdef DRAWING2 + if (l>1.4) { + penthickness(3); + vP.MoveTo(),vQ.LineTo(); + penthickness(1); + cout << " l = " << l << Number(vP) << " edge = " << Number(vQ) << endl; + } +#endif + nbedges++; + k = (int) ((l - lmin)*delta); + k = Min(Max(k,0L),kmax); + histo[k]++; + } + } + } + cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <0); + for (int i=0;i0); + for (int i=0;iCracked()) k++; + if( k==0) return 0; + CurrentTh = this; + cout << " Nb of Cracked Edges = " << k << endl; + NbCrackedEdges =k; + CrackedEdges = new CrackedEdge[k]; + // new edge + Edge * e = new Edge[ nbe + k]; + + // copy + for (i=0;iCracked()) + { + e[nbe] = e[i]; + // return the edge + e[nbe].v[0] = e[i].v[1]; + e[nbe].v[1] = e[i].v[0]; + e[nbe].on = e[i].on->link ; // fqux + CrackedEdges[k++]=CrackedEdge(edges,i,nbe); + nbe++; + } + ReMakeTriangleContainingTheVertex() ; + // + int nbcrakev =0; + Vertex *vlast = vertices + nbv; + Vertex *vend = vertices + nbvx; // end of array + for (int iv=0;iv= 0 ) && (i <3)); + // turn around the vertex v + TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge + int k=0; + do { + int kv = VerticesOfTriangularEdge[ta][1]; + k++; + Triangle * tt (ta); + if ( ta.Cracked() ) + { + TriangleAdjacent tta=(ta.Adj()); + assert(tta.Cracked()); + if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart + if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use + kk++;// number of cracked edge view + } + if ( tt->link ) { // if good triangles store the value + int it = Number(tt); + assert(it < nt); + (*tt)(kv)= vv; // Change the vertex of triangle + if(vvReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref + // tt->SetTriangleContainingTheVertex(); + kkk++; + } else if (kk) { // crack + boundary + if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use + } + + ta = Next(ta).Adj(); + } while ( (tbegin != ta)); + assert(k); + if (kc) nbcrakev++; + } + + if ( nbcrakev ) + for (int iec =0;iec < NbCrackedEdges; iec ++) + CrackedEdges[iec].Set(); + + // set the ref + cout << " set the ref " << endl ; + NbCrackedVertices = nbcrakev; + // int nbvo = nbv; + nbv = vlast - vertices; + int nbnewv = nbv - nbv; // nb of new vrtices + if (nbcrakev && verbosity > 1 ) + cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl; + // all the new vertices are on geometry + // BOFBO-- A VOIR + if (nbnewv) + { // + Int4 n = nbnewv+NbVerticesOnGeomVertex; + Int4 i,j,k; + VertexOnGeom * vog = new VertexOnGeom[n]; + for ( i =0; i= LastOld) + { // a new vertex + Int4 old = v->ReferenceNumber ; // the old same vertex + Int4 i = ( v - LastOld); + // if the old is on vertex => warning + // else the old is on edge => ok + vog[i] = vog[old]; + // vog[i].mv = v; + //g[i].ge = ; + //og[i].abcisse = ; + } + + } + } + + NbVerticesOnGeomVertex = n; + } + SetVertexFieldOn(); + + + if (vlast >= vend) + { + cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl; + MeshError(555,this); + } + cout << " NbCrackedVertices " << NbCrackedVertices << endl; + CurrentTh = CurrentThOld; + return NbCrackedVertices; +} + + Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb) +: Gh(*(new Geometry())), BTh(*this) +{ // truncature + // + + char cname[] = "trunc"; + + int i,k,itadj; + int kt=0; + int * kk = new int [Tho.nbv]; + Int4 * reft = new Int4[Tho.nbt]; + Int4 nbInT = Tho.ConsRefTriangle(reft); + Int4 * refv = new Int4[Tho.nbv]; + + for (i=0;i=0 && flag[i]) + { + const Triangle & t = Tho.triangles[i]; + kt++; + kk[Tho.Number(t[0])]=1; + kk[Tho.Number(t[1])]=1; + kk[Tho.Number(t[2])]=1; + itadj=Tho.Number(t.TriangleAdj(0)); + if ( reft[itadj] >=0 && !flag[itadj]) + { nbNewBedge++; + refv[Tho.Number(t[VerticesOfTriangularEdge[0][0]])]=bb[i]; + refv[Tho.Number(t[VerticesOfTriangularEdge[0][1]])]=bb[i]; + } + itadj=Tho.Number(t.TriangleAdj(1)); + if ( reft[itadj] >=0 && !flag[itadj]) + { nbNewBedge++; + refv[Tho.Number(t[VerticesOfTriangularEdge[1][0]])]=bb[i]; + refv[Tho.Number(t[VerticesOfTriangularEdge[1][1]])]=bb[i];} + itadj=Tho.Number(t.TriangleAdj(2)); + if ( reft[itadj] >=0 && !flag[itadj]) + { nbNewBedge++; + refv[Tho.Number(t[VerticesOfTriangularEdge[2][0]])]=bb[i]; + refv[Tho.Number(t[VerticesOfTriangularEdge[2][1]])]=bb[i];} + } + k=0; + for (i=0;i=0) + kk[i]=k++; + cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl; + cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl; + cout << " number of New boundary edge " << nbNewBedge << endl; + Int4 inbvx =k; + PreInit(inbvx,cname); + for (i=0;i=0) + { + vertices[nbv] = Tho.vertices[i]; + if (!vertices[nbv].ref()) + vertices[nbv].ReferenceNumber = refv[i]; + nbv++; + } + assert(inbvx == nbv); + for (i=0;i=0 && flag[i]) + { + const Triangle & t = Tho.triangles[i]; + int i0 = Tho.Number(t[0]); + int i1 = Tho.Number(t[1]); + int i2 = Tho.Number(t[2]); + assert(i0>=0 && i1 >= 0 && i2 >= 0); + assert(i0link); + +} + +Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { // in: B + // out: t + // out : dete[3] + // t the triangle and s0,s1,s2 the 3 vertices of t + // in dete[3] = { det(B,s1,s2) , det(s0,B,s2), det(s0,s1,B)} + // with det(a,b,c ) = -1 if one of 3 vertices a,b,c is NULL + Triangle * t=0; + int j,jp,jn,jj; + if (tstart) + t=tstart; + else + { + assert(quadtree); + Vertex *a = quadtree->NearestVertex(B.x,B.y) ; + + if (! a || !a->t ) { + if (a) + {cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl; + cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;} + cerr << " Pb with " << B << toR2(B) << endl; + MeshError(7777); + } + assert(a>= vertices && a < vertices+nbv); +#ifdef DRAWING1 + a->Draw(); +#endif + // int k=0; + t = a->t; + assert(t>= triangles && t < triangles+nbt); + + } + Icoor2 detop ; + int kkkk =0; // number of test triangle + + while ( t->det < 0) + { // the initial triangles is outside + int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + assert(k0>=0); // k0 the NULL vertex + int k1=NextVertex[k0],k2=PreviousVertex[k0]; + dete[k0]=det(B,(*t)[k1],(*t)[k2]); + dete[k1]=dete[k2]=-1; + if (dete[k0] > 0) // outside B + return t; + t = t->TriangleAdj(OppositeEdge[k0]); + assert(kkkk++ < 2); + } + + jj=0; + detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B); + + while(t->det > 0 ) + { + assert( kkkk++ < 2000 ); + j= OppositeVertex[jj]; + +#ifdef DRAWING1 + t->Draw(); +#endif + dete[j] = detop; //det(*b,*s1,*s2); + jn = NextVertex[j]; + jp = PreviousVertex[j]; + dete[jp]= det(*(*t)(j),*(*t)(jn),B); + dete[jn] = t->det-dete[j] -dete[jp]; + +#ifdef DEBUG + const Vertex * s0 = (*t)(0); + const Vertex * s1 = (*t)(1); + const Vertex * s2 = (*t)(2); + assert(dete[0] == det(B ,*s1,*s2)); + assert(dete[1] == det(*s0,B ,*s2)); + assert(dete[2] == det(*s0,*s1,B )); + assert(t->det== (dete[0] + dete[1] +dete[2])); +#endif + // count the number k of dete <0 + int k=0,ii[3]; + if (dete[0] < 0 ) ii[k++]=0; + if (dete[1] < 0 ) ii[k++]=1; + if (dete[2] < 0 ) ii[k++]=2; + // 0 => ok + // 1 => go in way 1 + // 2 => two way go in way 1 or 2 randomly + + if (k==0) + break; + if (k == 2 && BinaryRand()) + Exchange(ii[0],ii[1]); + assert ( k < 3); + TriangleAdjacent t1 = t->Adj(jj=ii[0]); + if ((t1.det() < 0 ) && (k == 2)) + t1 = t->Adj(jj=ii[1]); + t=t1; + j=t1;// for optimisation we now the -det[OppositeVertex[j]]; + detop = -dete[OppositeVertex[jj]]; + jj = j; + } + + if (t->det<0) // outside triangle + dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; + // NbOfTriangleSearchFind += kkkk; + return t; +} + +} + Index: /issm/trunk/src/c/Bamgx/Mesh2.h =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2740) +++ /issm/trunk/src/c/Bamgx/Mesh2.h (revision 2740) @@ -0,0 +1,1503 @@ +// -*- Mode : c++ -*- +// +// SUMMARY : +// USAGE : +// ORG : +// AUTHOR : Frederic Hecht +// E-MAIL : hecht@ann.jussieu.fr +// + +/* + + This file is part of Freefem++ + + Freefem++ is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + Freefem++ is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with Freefem++; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + */ +#ifndef _MESH2_H_ +#define _MESH2_H_ + +#include +#include +#include +#include +#if (defined(unix) || defined(__unix)) && !defined(__AIX) +#define SYSTIMES +#include +#include +#endif +#ifdef DRAWING +#include "rgraph.hpp" +#endif + +extern int SHOW; +#include "meshtype.h" + +#include "error.hpp" + + +#include "R2.h" + +namespace bamg { + + + +const double Pi = 3.14159265358979323846264338328; +const float fPi = 3.14159265358979323846264338328; + + +class MeshIstream; +class OFortranUnFormattedFile; +class IFortranUnFormattedFile; + +extern int hinterpole; + + +typedef P2 I2; + +inline int BinaryRand(){ +#ifdef RAND_MAX + const long HalfRandMax = RAND_MAX/2; + return rand() R2; +typedef P2xP2 I2xI2; +typedef P2 R2xR2; + +} + +#include "Metric.h" + +namespace bamg { +inline float OppositeAngle(float a) + {return a<0 ? fPi + a :a - fPi ;} +inline double OppositeAngle(double a) + {return a<0 ? Pi + a :a - Pi ;} + +#ifdef DRAWING +extern Real4 xGrafCoef,yGrafCoef,xGrafOffSet,yGrafOffSet; +extern R2 GrafPMin,GrafPMax; +extern Real8 Grafh; +#endif + +Icoor2 inline det(const I2 &a,const I2 & b,const I2 &c) +{ + register Icoor2 bax = b.x - a.x ,bay = b.y - a.y; + register Icoor2 cax = c.x - a.x ,cay = c.y - a.y; + return bax*cay - bay*cax;} + + + +// def de numerotation dans un triangles +static const Int2 VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}}; +static const Int2 EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}}; +static const Int2 OppositeVertex[3] = {0,1,2}; +static const Int2 OppositeEdge[3] = {0,1,2}; +static const Int2 NextEdge[3] = {1,2,0}; +static const Int2 PreviousEdge[3] = {2,0,1}; +static const Int2 NextVertex[3] = {1,2,0}; +static const Int2 PreviousVertex[3] = {2,0,1}; + +Int4 AGoodNumberPrimeWith(Int4 n); + +// remark all the angle are in radian beetwen [-Pi,Pi] + + +class Geometry; +//static Geometry *NULLGeometry=0; +class Triangles; +class Triangle; +class QuadTree; +class GeometricalEdge; +class VertexOnGeom; +class VertexOnEdge; +///////////////////////////////////////////////////////////////////////////////////// +const int IsVertexOnGeom = 8; +const int IsVertexOnVertex = 16; +const int IsVertexOnEdge = 32; +///////////////////////////////////////////////////////////////////////////////////// +#ifndef NOTFREEFEM +class ErrorMesh : public Error +{ +public: + Triangles *Th; + ErrorMesh(const char * Text,int l,Triangles * TTh=0, const char *t2="") : + Error(MESH_ERROR,"Meshing error: ",Text,"\n number : ",l,", ",t2),Th(TTh) {} +}; +#endif + +class Direction { // + private: + Icoor1 dir; + public: + Direction(): dir(MaxICoor){}; // no direction set + Direction(Icoor1 i,Icoor1 j) { Icoor2 n2 = 2*(Abs(i)+Abs(j)); + Icoor2 r = MaxICoor* (Icoor2) i; + Icoor1 r1 = (Icoor1) (2*(r/ n2)); // odd number + dir = (j>0) ? r1 : r1+1; // odd -> j>0 even -> j<0 + } + int sens( Icoor1 i,Icoor1 j) { int r =1; + if (dir!= MaxICoor) { + Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1); + r = (x*i + y*j) >=0;} + return r;} +#ifdef DRAWING + void Draw() { + if (dir!= MaxICoor) { + Icoor2 x(dir/2),y1(MaxICoor/2-Abs(x)),y(dir%2?-y1:y1); + R2 D(x,y); + double eps = Grafh/Norme2(D)/20; + D = D*eps; + rmoveto(D.x,D.y); + } + } +#endif + + + + +}; +///////////////////////////////////////////////////////////////////////////////////// +class Vertex {public: + I2 i; // allow to use i.x, and i.y in long int (beware of scale and centering) + R2 r; // allow to use r.x, and r.y in double + Metric m; + Int4 ReferenceNumber; + Direction DirOfSearch; + union { + Triangle * t; // one triangle which contained the vertex + Int4 color; + Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed + VertexOnGeom * on; // if vint 8; // set with Triangles::SetVertexFieldOn() + Vertex * onbv; // if vint == 16 on Background vertex Triangles::SetVertexFieldOnBTh() + VertexOnEdge * onbe; // if vint == 32 on Background edge + }; + Int1 vint; // the vertex number in triangle; varies between 0 and 2 in t + operator I2 () const {return i;} // operator de cast + operator const R2 & () const {return r;}// operator de cast +// operator R2 & () {return r;}// operator de cast + Real8 operator()(R2 x) const { return m(x);} + operator Metric () const {return m;}// operator de cast + Int4 Optim(int = 1,int =0); + // Vertex(){} + // ~Vertex(){} + Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8 =1); + int ref() const { return ReferenceNumber;} + + friend std::ostream& operator <<(std::ostream& f, const Vertex & v) + {f << "(" << v.i << "," << v.r << MatVVP2x2(v.m) << ")" ; return f;} + inline void Set(const Vertex & rec,const Triangles &,Triangles &); + +#ifdef DRAWING + void Draw(Int4 =-1) const ; + void MoveTo() const { rmoveto(r.x,r.y); } + void LineTo() const { rlineto(r.x,r.y); } +#endif +}; + +double QuadQuality(const Vertex &,const Vertex &,const Vertex &,const Vertex &); + +// extern Vertex *Meshend , *Meshbegin; + +///////////////////////////////////////////////////////////////////////////////////// +class TriangleAdjacent { + friend std::ostream& operator <<(std::ostream& f, const TriangleAdjacent & ta) + {f << "{" << ta.t << "," << ((int) ta.a) << "}" ; + return f;} + +public: + Triangle * t; // le triangle + int a; // le numero de l arete + + TriangleAdjacent(Triangle * tt,int aa): t(tt),a(aa &3) {}; + TriangleAdjacent() {}; + + operator Triangle * () const {return t;} + operator Triangle & () const {return *t;} + operator int() const {return a;} + TriangleAdjacent & operator++() + { + a= NextEdge[a]; + return *this;} + TriangleAdjacent operator--() + { + a= PreviousEdge[a]; + return *this;} + inline TriangleAdjacent Adj() const ; + int swap(); + inline void SetAdj2(const TriangleAdjacent& , int =0); + inline Vertex * EdgeVertex(const int &) const ; + inline Vertex * OppositeVertex() const ; + inline Icoor2 & det() const; + inline int Locked() const ; + inline int GetAllFlag_UnSwap() const ; + inline void SetLock(); + inline int MarkUnSwap() const; + inline void SetMarkUnSwap(); + inline void SetCracked(); + inline int Cracked() const ; +};// end of Vertex class + + +///////////////////////////////////////////////////////////////////////////////////// +class Edge { public: + Vertex * v[2]; + Int4 ref; + GeometricalEdge * on; + Vertex & operator[](int i){return *v[i];}; + Vertex * operator()(int i){return v[i];}; + + void ReNumbering(Vertex *vb,Vertex *ve, Int4 *renu) + { + if (v[0] >=vb && v[0] =vb && v[1] no continuite + GeometricalEdge * Adj [2]; + int SensAdj[2]; +// private: + int flag ; + public: + GeometricalEdge * link; // if Cracked() or Equi() + +// end of data + + GeometricalVertex & operator[](int i){return *v[i];}; + const GeometricalVertex & operator[](int i) const { return *v[i];}; + GeometricalVertex * operator()(int i){return v[i];}; + // inline void Set(const Geometry &,Int4,Geometry &); + + R2 F(Real8 theta) const ; // parametrization of the curve edge + Real8 R1tg(Real8 theta,R2 &t) const ; // 1/radius of curvature + tangente + int Cracked() const {return flag & 1;} + int Dup() const { return flag & 32;} + int Equi()const {return flag & 2;} + int ReverseEqui()const {return flag & 128;} + int TgA()const {return flag &4;} + int TgB()const {return flag &8;} + int Tg(int i) const{return i==0 ? TgA() : TgB();} + int Mark()const {return flag &16;} + int Required() { return flag & 64;} + void SetCracked() { flag |= 1;} + void SetDup() { flag |= 32;} // not a real edge + void SetEqui() { flag |= 2;} + void SetTgA() { flag|=4;} + void SetTgB() { flag|=8;} + void SetMark() { flag|=16;} + void SetUnMark() { flag &= 1007 /* 1023-16*/;} + void SetRequired() { flag |= 64;} + void SetReverseEqui() {flag |= 128;} + + inline void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); + +#ifdef DRAWING + void Draw(Int4 =-1); +#endif + +}; + +class Curve {public: + GeometricalEdge * be,*ee; // begin et end edge + int kb,ke; // begin vetex and end vertex + Curve *next; // next curve equi to this + bool master; // true => of equi curve point on this curve + inline void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); + Curve() : be(0),ee(0),kb(0),ke(0),next(0),master(true) {} + void Reverse() { Exchange(be,ee); Exchange(kb,ke);} // revese the sens of the curse +}; + + + +///////////////////////////////////////////////////////////////////////////////////// +class Triangle { + friend class TriangleAdjacent; + friend std::ostream& operator <<(std::ostream& f, const Triangle & ta); + + + private: // les arete sont opposes a un sommet + Vertex * ns [3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer + Triangle * at [3]; // nu triangle adjacent + Int1 aa[3]; // les nu des arete dans le triangles (mod 4) + public: + Icoor2 det; // determinant du triangle (2 fois l aire des vertices entieres) + union { + Triangle * link ; + Int4 color; + }; + void SetDet() { + if(ns[0] && ns[1] && ns[2]) det = bamg::det(*ns[0],*ns[1],*ns[2]); + else det = -1; } + Triangle() {} + Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k); + Triangle(Vertex *v0,Vertex *v1,Vertex *v2); + inline void Set(const Triangle &,const Triangles &,Triangles &); + inline int In(Vertex *v) const { return ns[0]==v || ns[1]==v || ns[2]==v ;} + TriangleAdjacent FindBoundaryEdge(int ) const; + + void ReNumbering(Triangle *tb,Triangle *te, Int4 *renu) + { + if (link >=tb && link =tb && at[0] =tb && at[1] =tb && at[2] =vb && ns[0] =vb && ns[1] =vb && ns[2] at[aatt]=this; + tt->aa[aatt]=a + (aa[a] & 60 ) ;}// Copy all the mark + } + + void SetAdj2(Int1 a,Triangle *t,Int1 aat) + { at[a]=t;aa[a]=aat; + if(t) {t->at[aat]=this;t->aa[aat]=a;} + } + + void SetTriangleContainingTheVertex() + { + if (ns[0]) (ns[0]->t=this,ns[0]->vint=0); + if (ns[1]) (ns[1]->t=this,ns[1]->vint=1); + if (ns[2]) (ns[2]->t=this,ns[2]->vint=2); + } + + int swap(Int2 a1,int=0); + Int4 Optim(Int2 a,int =0); + + int Locked(int a)const { return aa[a]&4;} + int Hidden(int a)const { return aa[a]&16;} + int Cracked(int a) const { return aa[a] & 32;} + // for optimisation + int GetAllflag(int a){return aa[a] & 1020;} + void SetAllFlag(int a,int f){aa[a] = (aa[a] &3) + (1020 & f);} + + void SetHidden(int a){ + register Triangle * t = at[a]; + if(t) t->aa[aa[a] & 3] |=16; + aa[a] |= 16;} + void SetCracked(int a){ + register Triangle * t = at[a]; + if(t) t->aa[aa[a] & 3] |=32; + aa[a] |= 32;} + + double QualityQuad(int a,int option=1) const; + Triangle * Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; + + void SetLocked(int a){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] |=4; + aa[a] |= 4;} + + void SetMarkUnSwap(int a){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] |=8; + aa[a] |=8 ;} + + + void SetUnMarkUnSwap(int a){ + register Triangle * t = at[a]; + t->aa[aa[a] & 3] &=55; // 23 + 32 + aa[a] &=55 ;} + + + +#ifdef DEBUG + void inline checka(Int1 a); + void inline check(); +#endif + +#ifdef DRAWING + void Draw(Int4 i=-1) const; + int swapDRAW(Int2 a1); + +#endif + +}; // end of Triangle class + + + + +class ListofIntersectionTriangles { +///////////////////////////////////////////////////////////////////////////////////// +class IntersectionTriangles { +public: + Triangle *t; + Real8 bary[3]; // use if t != 0 + R2 x; + Metric m; + Real8 s;// abscisse curviline + Real8 sp; // len of the previous seg in m + Real8 sn;// len of the next seg in m +}; +///////////////////////////////////////////////////////////////////////////////////// +class SegInterpolation { + public: + GeometricalEdge * e; + Real8 sBegin,sEnd; // abscisse of the seg on edge parameter + Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length + int last;// last index in ListofIntersectionTriangles for this Sub seg of edge + R2 F(Real8 s){ + Real8 c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; + assert(lBegin<= s && s <=lEnd); + return e->F(sBegin*c0+sEnd*c1);} +}; + + int MaxSize; // + int Size; // + Real8 len; // + int state; + IntersectionTriangles * lIntTria; + int NbSeg; + int MaxNbSeg; + SegInterpolation * lSegsI; + public: + IntersectionTriangles & operator[](int i) {return lIntTria[i];} + operator int&() {return Size;} + ListofIntersectionTriangles(int n=256,int m=16) + : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , + NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]) + { + long int verbosity=0; + + if (verbosity>9) + std::cout << " construct ListofIntersectionTriangles" + << MaxSize << " " << MaxNbSeg<< std::endl;}; + ~ListofIntersectionTriangles(){ + if (lIntTria) delete [] lIntTria,lIntTria=0; + if (lSegsI) delete [] lSegsI,lSegsI=0;} + void init(){state=0;len=0;Size=0;} + + int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2); + int NewItem(R2,const Metric & ); + void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1) + { + long int verbosity=0; + + if (NbSeg>=MaxNbSeg) { + int mneo= MaxNbSeg; + MaxNbSeg *= 2; + if (verbosity>3) + std::cout <<" reshape lSegsI from " << mneo << " to " + << MaxNbSeg <3) + std::cout << " ListofIntersectionTriangles ReShape MaxSize " + << MaxSize << " -> " + << newsize << std::endl; + MaxSize = newsize; + delete [] lIntTria;// remove old + lIntTria = nw; // copy pointer + } + + void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); + Real8 Length(); + Int4 NewPoints(Vertex *,Int4 & nbv,Int4 nbvx); +}; + + +///////////////////////////////////////////////////////////////////////////////////// +class GeometricalSubDomain { +public: + GeometricalEdge *edge; + int sens; // -1 or 1 + Int4 ref; + inline void Set(const GeometricalSubDomain &,const Geometry & ,const Geometry &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class SubDomain { +public: + Triangle * head; + Int4 ref; + int sens; // -1 or 1 + Edge * edge; // to geometrical + inline void Set(const Triangles &,Int4,Triangles &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class VertexOnGeom { public: + + Vertex * mv; + Real8 abscisse; + union{ + GeometricalVertex * gv; // if abscisse <0; + GeometricalEdge * ge; // if abscisse in [0..1] + }; + inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); + int OnGeomVertex()const {return this? abscisse <0 :0;} + int OnGeomEdge() const {return this? abscisse >=0 :0;} + VertexOnGeom(): mv(0),abscisse(0){gv=0;} + VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} + // std::cout << " mv = " <Required():0):0 )) : 0;} + void SetOn(){mv->on=this;mv->vint=IsVertexOnGeom;} + friend std::ostream& operator <<(std::ostream& f, const VertexOnGeom & vog){ + f << vog.abscisse << " " << vog.mv << " " << vog.gv << " ; "; + if (vog.abscisse < 0) f << *vog.gv << " ;; " ; + // else f << *vog.ge << " ;; " ; + return f;} + inline void Set(const Triangles &,Int4,Triangles &); + +}; +///////////////////////////////////////////////////////////////////////////////////// +class VertexOnVertex {public: + Vertex * v, *bv; + VertexOnVertex(Vertex * w,Vertex *bw) :v(w),bv(bw){} + VertexOnVertex() {}; + inline void Set(const Triangles &,Int4,Triangles &); + void SetOnBTh(){v->onbv=bv;v->vint=IsVertexOnVertex;} +}; +///////////////////////////////////////////////////////////////////////////////////// +class VertexOnEdge {public: + Vertex * v; + Edge * be; + Real8 abcisse; + VertexOnEdge( Vertex * w, Edge *bw,Real8 s) :v(w),be(bw),abcisse(s) {} + VertexOnEdge(){} + inline void Set(const Triangles &,Int4,Triangles &); + void SetOnBTh(){v->onbe=this;v->vint=IsVertexOnEdge;} + Vertex & operator[](int i) const { return (*be)[i];} + operator Real8 () const { return abcisse;} + operator Vertex * () const { return v;} + operator Edge * () const { return be;} +}; + + inline TriangleAdjacent FindTriangleAdjacent(Edge &E); + inline Vertex * TheVertex(Vertex * a); // for remove crak in mesh +///////////////////////////////////////////////////////////////////////////////////// + +class CrackedEdge { // a small class to store on crack an uncrack information + friend class Triangles; + friend std::ostream& operator <<(std::ostream& f, const Triangles & Th) ; + class CrackedTriangle { + friend class Triangles; + friend class CrackedEdge; + friend std::ostream& operator <<(std::ostream& f, const Triangles & Th) ; + Triangle * t; // edge of triangle t + int i; // edge number of in triangle + Edge *edge; // the 2 edge + Vertex *New[2]; // new vertex number + CrackedTriangle() : t(0),i(0),edge(0) { New[0]=New[1]=0;} + CrackedTriangle(Edge * a) : t(0),i(0),edge(a) { New[0]=New[1]=0;} + void Crack(){ + Triangle & T(*t); + int i0=VerticesOfTriangularEdge[i][0]; + int i1=VerticesOfTriangularEdge[i][0]; + assert(New[0] && New[1]); + T(i0) = New[0]; + T(i1) = New[1];} + void UnCrack(){ + Triangle & T(*t); + int i0=VerticesOfTriangularEdge[i][0]; + int i1=VerticesOfTriangularEdge[i][0]; + assert(New[0] && New[1]); + T(i0) = TheVertex(T(i0)); + T(i1) = TheVertex(T(i1));} + void Set() { + TriangleAdjacent ta ( FindTriangleAdjacent(*edge)); + t = ta; + i = ta; + + New[0] = ta.EdgeVertex(0); + New[1] = ta.EdgeVertex(1); + // warning the ref + + } + + }; // end of class CrackedTriangle + public: + CrackedTriangle a,b; + CrackedEdge() :a(),b() {} + CrackedEdge(Edge * start, Int4 i,Int4 j) : a(start+i),b(start+j) {}; + CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; + + void Crack() { a.Crack(); b.Crack();} + void UnCrack() { a.UnCrack(); b.UnCrack();} + void Set() { a.Set(), b.Set();} +}; + +///////////////////////////////////////////////////////////////////////////////////// +class Triangles { +public: + + enum TypeFileMesh { + AutoMesh=0,BDMesh=1,NOPOMesh=2,amMesh=3,am_fmtMesh=4,amdbaMesh=5, + ftqMesh=6,mshMesh=7}; + + int static counter; // to kown the number of mesh in memory + int OnDisk; // true if on disk + Geometry & Gh; // Geometry + Triangles & BTh; // Background Mesh Bth==*this =>no background + + Int4 NbRef; // counter of ref on the this class if 0 we can delete + Int4 nbvx,nbtx; // nombre max de sommets , de triangles + + Int4 nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, + // of initial vertices, of edges with reference, + Int4 NbOfQuad; // nb of quadrangle + + Int4 NbSubDomains; // + Int4 NbOutT; // Nb of oudeside triangle + Int4 NbOfTriangleSearchFind; + Int4 NbOfSwapTriangle; + char * name, *identity; + Vertex * vertices; // data of vertices des sommets + + Int4 NbVerticesOnGeomVertex; + VertexOnGeom * VerticesOnGeomVertex; + + Int4 NbVerticesOnGeomEdge; + VertexOnGeom * VerticesOnGeomEdge; + + Int4 NbVertexOnBThVertex; + VertexOnVertex *VertexOnBThVertex; + + Int4 NbVertexOnBThEdge; + VertexOnEdge *VertexOnBThEdge; + + + Int4 NbCrackedVertices; + + + Int4 NbCrackedEdges; + CrackedEdge *CrackedEdges; + + + R2 pmin,pmax; // extrema + Real8 coefIcoor; // coef to integer Icoor1; + + Triangle * triangles; + Edge * edges; + + QuadTree *quadtree; + Vertex ** ordre; + SubDomain * subdomains; + ListofIntersectionTriangles lIntTria; +// end of variable + + Triangles(Int4 i);//:BTh(*this),Gh(*new Geometry()){PreInit(i);} + + + ~Triangles(); + Triangles(const char * ,Real8=-1) ; + Triangles(double* elements, double* x, double* y, int nel, int nods,int NBV, Real8=-1); + + Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) + :Gh(BT.Gh),BTh(BT) { + try {GeomToTriangles1(nbvx,keepBackVertices);} + catch(...) { this->~Triangles(); throw; } } + + Triangles(Int4 nbvx,Geometry & G) + :Gh(G),BTh(*this){ + try { GeomToTriangles0(nbvx);} + catch(...) { this->~Triangles(); throw; } } + Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR + // Triangles(Triangles &){ std::cerr << " BUG call copy opretor of Triangles" << std::endl;MeshError(111);} + Triangles(const Triangles &,const int *flag,const int *bb); // truncature + + int gibbsv (long int* ptvoi, long int* vois,long int* lvois,long int* w,long int* v); + int gibbs(void); + + void SetIntCoor(const char * from =0); + + // void RandomInit(); + // void CubeInit(int ,int); + + Real8 MinimalHmin() {return 2.0/coefIcoor;} + Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} + const Vertex & operator[] (Int4 i) const { return vertices[i];}; + Vertex & operator[](Int4 i) { return vertices[i];}; + const Triangle & operator() (Int4 i) const { return triangles[i];}; + Triangle & operator()(Int4 i) { return triangles[i];}; + I2 toI2(const R2 & P) const { + return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) + ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} + R2 toR2(const I2 & P) const { + return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);} + void Add( Vertex & s,Triangle * t,Icoor2 * =0) ; + void Insert(); + // void InsertOld(); + void ForceBoundary(); + void Heap(); + void FindSubDomain(int ); + Int4 ConsRefTriangle(Int4 *) const; + void ShowHistogram() const; + void ShowRegulaty() const; // Add FH avril 2007 +// void ConsLinkTriangle(); + + void ReMakeTriangleContainingTheVertex(); + void UnMarkUnSwapTriangle(); + void SmoothMetric(Real8 raisonmax) ; + void BoundAnisotropy(Real8 anisomax,double hminaniso= 1e-100) ; + void MaxSubDivision(Real8 maxsubdiv); + void WriteMetric(std::ostream &,int iso) ; + Edge** MakeGeometricalEdgeToEdge(); + void SetVertexFieldOn(); + void SetVertexFieldOnBTh(); + Int4 SplitInternalEdgeWithBorderVertices(); + void MakeQuadrangles(double costheta); + int SplitElement(int choice); + void MakeQuadTree(); + void NewPoints( Triangles &,int KeepBackVertex =1 ); + Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; + void NewPointsOld( Triangles & ); + void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);} + void ReNumberingTheTriangleBySubDomain(bool justcompress=false); + void ReNumberingVertex(Int4 * renu); + void SmoothingVertex(int =3,Real8=0.3); + Metric MetricAt (const R2 &) const; + GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, + Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); + + + void WriteElements(std::ostream& f,Int4 * reft ,Int4 nbInT) const; + + + Int4 Number(const Triangle & t) const { return &t - triangles;} + Int4 Number(const Triangle * t) const { return t - triangles;} + Int4 Number(const Vertex & t) const { return &t - vertices;} + Int4 Number(const Vertex * t) const { return t - vertices;} + Int4 Number(const Edge & t) const { return &t - edges;} + Int4 Number(const Edge * t) const { return t - edges;} + Int4 Number2(const Triangle * t) const { + // if(t>= triangles && t < triangles + nbt ) + return t - triangles; + // else return t - OutSidesTriangles; + } + + Vertex * NearestVertex(Icoor1 i,Icoor1 j) ; + Triangle * FindTriangleContening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; + void Write(const char * filename,const TypeFileMesh type = AutoMesh); + void Write_am_fmt(std::ostream &) const ; + void Write_am(std::ostream &) const ; + void Write_ftq(std::ostream &) const ; + void Write_nopo(std::ostream &) const ; + void Write_msh(std::ostream &) const ; + void Write_amdba(std::ostream &) const ; + + void Read(MeshIstream &,int version,Real8 cutoffradian); + void ReadFromMatlabMesh(double* elements,double* x,double* y,int nel,int nods,int inbvx,Real8 cutoffradian); + void Read_am_fmt(MeshIstream &); + void Read_amdba(MeshIstream &); + void Read_am(MeshIstream &); + void Read_nopo(MeshIstream &); + void Read_ftq(MeshIstream &); + void Read_msh(MeshIstream &); + + void ReadMetric(const char * fmetrix,const Real8 hmin,const Real8 hmax,const Real8 coef); + void IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols, + const Real8 hmin,const Real8 hmax, const Real8 coef, + const Real8 anisomax,const Real8 CutOff=1.e-4,const int NbJacobi=1, + const int DoNormalisation=1, + const double power=1.0, + const int choise=0); + void IntersectGeomMetric(const Real8 err,const int iso); + + + int isCracked() const {return NbCrackedVertices != 0;} + int Crack(); + int UnCrack(); + +#ifdef DEBUG + void inline Check(); +#endif +#ifdef DRAWING + void Draw() const ; + void InitDraw() const ; + void inquire() ; +#endif + friend std::ostream& operator <<(std::ostream& f, const Triangles & Th); + void Write(const char * filename); + void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo + void FillHoleInMesh() ; + int CrackMesh(); + private: + void GeomToTriangles1(Int4 nbvx,int KeepBackVertices=1);// the real constructor mesh adaption + void GeomToTriangles0(Int4 nbvx);// the real constructor mesh generator + void PreInit(Int4,char * =0 ); + // + void Write_nop5(OFortranUnFormattedFile * f, + Int4 &lnop5,Int4 &nef,Int4 &lgpdn,Int4 ndsr) const ; + + +}; // End Class Triangles + + +///////////////////////////////////////////////////////////////////////////////////// +class Geometry { +public: + int OnDisk; + Int4 NbRef; // counter of ref on the this class if 0 we can delete + + char *name; + Int4 nbvx,nbtx; // nombre max de sommets , de Geometry + Int4 nbv,nbt,nbiv,nbe; // nombre de sommets, de Geometry, de sommets initiaux, + Int4 NbSubDomains; // + Int4 NbEquiEdges; + Int4 NbCrackedEdges; +// Int4 nbtf;// de triangle frontiere + Int4 NbOfCurves; + int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } + GeometricalVertex * vertices; // data of vertices des sommets + Triangle * triangles; + GeometricalEdge * edges; + QuadTree *quadtree; + GeometricalSubDomain *subdomains; + Curve *curves; + ~Geometry(); + Geometry(const Geometry & Gh); //Copy Operator + Geometry(int nbg,const Geometry ** ag); // intersection operator + + R2 pmin,pmax; // extrema + Real8 coefIcoor; // coef to integer Icoor1; + Real8 MaximalAngleOfCorner; + +// end of data + + + I2 toI2(const R2 & P) const { + return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) + ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} + + Real8 MinimalHmin() {return 2.0/coefIcoor;} + Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} + void ReadGeometry(const char * ) ; + void ReadGeometry(MeshIstream & ,const char *) ; + + void EmptyGeometry(); + Geometry() {EmptyGeometry();}// empty Geometry + void AfterRead(); + Geometry(const char * filename) {EmptyGeometry();OnDisk=1;ReadGeometry(filename);AfterRead();} + + void ReadMetric(const char *,Real8 hmin,Real8 hmax,Real8 coef); + const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; + GeometricalVertex & operator[](Int4 i) { return vertices[i];}; + const GeometricalEdge & operator() (Int4 i) const { return edges[i];}; + GeometricalEdge & operator()(Int4 i) { return edges[i];}; + Int4 Number(const GeometricalVertex & t) const { return &t - vertices;} + Int4 Number(const GeometricalVertex * t) const { return t - vertices;} + Int4 Number(const GeometricalEdge & t) const { return &t - edges;} + Int4 Number(const GeometricalEdge * t) const { return t - edges;} + Int4 Number(const Curve * c) const { return c - curves;} + + void UnMarkEdges() { + for (Int4 i=0;i=0) { + t = at[a]; + // if (t-this<0) return 0; + v2 = ns[VerticesOfTriangularEdge[a][0]]; + v0 = ns[VerticesOfTriangularEdge[a][1]]; + v1 = ns[OppositeEdge[a]]; + v3 = t->ns[OppositeEdge[aa[a]&3]]; + } + } + return t; +} + +inline double Triangle::QualityQuad(int a,int option) const +{ // first do the logique part + double q; + if (!link || aa[a] &4) + q= -1; + else { + Triangle * t = at[a]; + if (t-this<0) q= -1;// because we do 2 times + else if (!t->link ) q= -1; + else if (aa[0] & 16 || aa[1] & 16 || aa[2] & 16 || t->aa[0] & 16 || t->aa[1] & 16 || t->aa[2] & 16 ) + q= -1; + else if(option) + { + const Vertex & v2 = *ns[VerticesOfTriangularEdge[a][0]]; + const Vertex & v0 = *ns[VerticesOfTriangularEdge[a][1]]; + const Vertex & v1 = *ns[OppositeEdge[a]]; + const Vertex & v3 = * t->ns[OppositeEdge[aa[a]&3]]; + q = QuadQuality(v0,v1,v2,v3); // do the float part + } + else q= 1; + } + return q; +} + + +inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ) + { + *this = rec; + } +inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ) + { + *this = rec; + } +inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) + { + *this = Th.edges[i]; + v[0] = ThNew.vertices + Th.Number(v[0]); + v[1] = ThNew.vertices + Th.Number(v[1]); + if (on) + on = ThNew.Gh.edges+Th.Gh.Number(on); + if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]); + if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]); + + } +inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew) + { + *this = rec; + v[0] = GhNew.vertices + Gh.Number(v[0]); + v[1] = GhNew.vertices + Gh.Number(v[1]); + if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]); + if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]); + } + +inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew) +{ + *this = rec; + be = GhNew.edges + Gh.Number(be); + ee = GhNew.edges + Gh.Number(ee); + if(next) next= GhNew.curves + Gh.Number(next); +} + +inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew) + { + *this = rec; + if ( ns[0] ) ns[0] = ThNew.vertices + Th.Number(ns[0]); + if ( ns[1] ) ns[1] = ThNew.vertices + Th.Number(ns[1]); + if ( ns[2] ) ns[2] = ThNew.vertices + Th.Number(ns[2]); + if(at[0]) at[0] = ThNew.triangles + Th.Number(at[0]); + if(at[1]) at[1] = ThNew.triangles + Th.Number(at[1]); + if(at[2]) at[2] = ThNew.triangles + Th.Number(at[2]); + if (link >= Th.triangles && link < Th.triangles + Th.nbt) + link = ThNew.triangles + Th.Number(link); + } +inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.VertexOnBThVertex[i]; + v = ThNew.vertices + Th.Number(v); + +} +inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.subdomains[i]; + assert( head - Th.triangles >=0 && head - Th.triangles < Th.nbt); + head = ThNew.triangles + Th.Number(head) ; + assert(edge - Th.edges >=0 && edge - Th.edges < Th.nbe); + edge = ThNew.edges+ Th.Number(edge); +} + inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew) +{ + *this = rec; + edge = Gh.Number(edge) + GhNew.edges; +} +inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) +{ + *this = Th.VertexOnBThEdge[i]; + v = ThNew.vertices + Th.Number(v); +} + +inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew) +{ + *this = rec; + mv = ThNew.vertices + Th.Number(mv); + if (gv) + if (abscisse < 0 ) + gv = ThNew.Gh.vertices + Th.Gh.Number(gv); + else + ge = ThNew.Gh.edges + Th.Gh.Number(ge); + +} +inline Real8 Edge::MetricLength() const + { + return LengthInterpole(v[0]->m,v[1]->m,v[1]->r - v[0]->r) ; + } + +inline void Triangles::ReMakeTriangleContainingTheVertex() + { + register Int4 i; + for ( i=0;iat[a]=ta.t; + t->aa[a]=ta.a|l;} + if(ta.t) { + ta.t->at[ta.a] = t ; + ta.t->aa[ta.a] = a| l ; + } +} + + +inline int TriangleAdjacent::Locked() const +{ return t->aa[a] &4;} +inline int TriangleAdjacent::Cracked() const +{ return t->aa[a] &32;} +inline int TriangleAdjacent::GetAllFlag_UnSwap() const +{ return t->aa[a] & 1012;} // take all flag except MarkUnSwap + +inline int TriangleAdjacent::MarkUnSwap() const +{ return t->aa[a] &8;} + +inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);} + +inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);} + +inline TriangleAdjacent TriangleAdjacent::Adj() const +{ return t->Adj(a);} + +inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const + {return t->ns[VerticesOfTriangularEdge[a][i]]; } +inline Vertex * TriangleAdjacent::OppositeVertex() const +{return t->ns[bamg::OppositeVertex[a]]; } +inline Icoor2 & TriangleAdjacent::det() const +{ return t->det;} +inline TriangleAdjacent Adj(const TriangleAdjacent & a) +{ return a.Adj();} + +inline TriangleAdjacent Next(const TriangleAdjacent & ta) +{ return TriangleAdjacent(ta.t,NextEdge[ta.a]);} + +inline TriangleAdjacent Previous(const TriangleAdjacent & ta) +{ return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);} + +inline void Adj(GeometricalEdge * & on,int &i) + {int j=i;i=on->SensAdj[i];on=on->Adj[j];} + +inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) +{ + Real4 ret; + I2 ia=va,ib=vb,ic=vc; + I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; + Icoor2 deta=Det(ab,ac); + if (deta <=0) ret = -1; + else { + Real8 a = sqrt((Real8) (ac,ac)), + b = sqrt((Real8) (bc,bc)), + c = sqrt((Real8) (ab,ab)), + p = a+b+c; + Real8 h= Max(Max(a,b),c),ro=deta/p; + ret = ro/h;} + return ret; +} + + +inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) { + Vertex *v=Th->vertices; + Int4 nbv = Th->nbv; + assert(i >=0 && j >=0 && k >=0); + assert(i < nbv && j < nbv && k < nbv); + ns[0]=v+i; + ns[1]=v+j; + ns[2]=v+k; + at[0]=at[1]=at[2]=0; + aa[0]=aa[1]=aa[2]=0; + det=0; +} + +inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){ + ns[0]=v0; + ns[1]=v1; + ns[2]=v2; + at[0]=at[1]=at[2]=0; + aa[0]=aa[1]=aa[2]=0; + if (v0) det=0; + else { + det=-1; + link=NULL;}; +} + +inline Real4 Triangle::qualite() +{ + return det < 0 ? -1 : bamg::qualite(*ns[0],*ns[1],*ns[2]); +} + +Int4 inline Vertex::Optim(int i,int koption) +{ + Int4 ret=0; + if ( t && (vint >= 0 ) && (vint <3) ) + { + ret = t->Optim(vint,koption); + if(!i) + { + t =0; // for no future optime + vint= 0; } + } + return ret; +} + +Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c) +{ + register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y; + register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y; + return bax*cay - bay*cax;} + + +void swap(Triangle *t1,Int1 a1, + Triangle *t2,Int1 a2, + Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); + + + +int inline TriangleAdjacent::swap() +{ return t->swap(a);} + + + +int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , + TriangleAdjacent & tt1,Icoor2 & dets1, + Icoor2 & detsa,Icoor2 & detsb, int & nbswap); + +int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; + +// inline bofbof FH +inline TriangleAdjacent FindTriangleAdjacent(Edge &E) + { + Vertex * a = E.v[0]; + Vertex * b = E.v[1]; + + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + assert(t && i>=0 && i < 3); + assert( a == (*t)(i)); + int k=0; + do { // turn around vertex in direct sens (trigo) + k++;assert(k< 20000); + // in no crack => ta.EdgeVertex(1) == a otherwise ??? + if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find + ta = ta.Adj(); + if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find + --ta; + } while (t != (Triangle *)ta); + assert(0); + return TriangleAdjacent(0,0);// error + } + +inline Vertex * TheVertex(Vertex * a) // give a unique vertex with smallest number +{ // in case on crack in mesh + Vertex * r(a), *rr; + Triangle * t = a->t; + int i = a->vint; + TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge + assert(t && i>=0 && i < 3); + assert( a == (*t)(i)); + int k=0; + do { // turn around vertex in direct sens (trigo) + k++;assert(k< 20000); + // in no crack => ta.EdgeVertex(1) == a + if ((rr=ta.EdgeVertex(0)) < r) r = rr; + ta = ta.Adj(); + if ((rr=ta.EdgeVertex(1)) < r) r =rr; + --ta; + } while (t != (Triangle*) ta); + return r; +} + +inline double CPUtime(){ +#ifdef SYSTIMES + struct tms buf; + if (times(&buf)!=-1) + return ((double)buf.tms_utime+(double)buf.tms_stime)/(long) sysconf(_SC_CLK_TCK); + else +#endif + return ((double) clock())/CLOCKS_PER_SEC; +} + +#ifdef DEBUG +void inline Triangle::checka(Int1 a) { + // verif de la coherence des adjacences de l arete a + a = a%4; + assert(a < 3 && a >= 0 ); + Triangle *t1=this,*t2=at[a]; + Int2 a1=a,a2=aa[a]%4; + + assert(a2 < 3 && a2 >= 0 ); + if (t2 && ( ((*t1).ns[VerticesOfTriangularEdge[a1][0]] != (*t2).ns[VerticesOfTriangularEdge[a2][1]]) + || ((*t1).ns[VerticesOfTriangularEdge[a1][1]] != (*t2).ns[VerticesOfTriangularEdge[a2][0]]))) + { + if (CurrentTh) std::cerr << " In Triangles beetween Triangle " << CurrentTh->Number(t1) << " and " + << CurrentTh->Number(t2) << std::endl; + std::cerr << "---- t1="<< t1 << " " << a1 <<", t2="<< t2 << " " << a2 << std::endl; + std::cerr <<"t1="<< t1 << " " << a1 << " " << t1->ns[VerticesOfTriangularEdge[a1][0]] + << " " << t1->ns[VerticesOfTriangularEdge[a1][1]] <Number(t1->ns[VerticesOfTriangularEdge[a1][0]]) + << " " << CurrentTh->Number(t1->ns[VerticesOfTriangularEdge[a1][1]]) <ns[VerticesOfTriangularEdge[a2][0]] + << " " << t2->ns[VerticesOfTriangularEdge[a2][1]] <Number(t2->ns[VerticesOfTriangularEdge[a2][0]]) + << " " << CurrentTh->Number(t2->ns[VerticesOfTriangularEdge[a2][1]]) <aa[a1]/4 == t2->aa[a2]/4); // lock compatibite +} + + +void inline Triangle::check() { + Icoor2 det2=0; + // std::cout << " check " << this << std::endl; + int infv=ns[0] ? (( ns[1] ? ( ns[2] ? -1 : 2) : 1 )) : 0; + if (det<0) { + if (infv<0 ) + { if (CurrentTh) std::cerr << " In Triangles " << CurrentTh->Number(this) << std::endl; + std::cerr << " det = " << det << " and " << infv << std::endl; + MeshError(5); + }} + else if (infv>=0 ) + { if (CurrentTh) std::cerr << " In Triangles " << CurrentTh->Number(this) << std::endl; + std::cerr << " det = " << det << " and " << infv << std::endl; + MeshError(5); + } + + if (det >=0) + if( det != (det2=bamg::det(*ns[0],*ns[1],*ns[2]))) + { // penthickness(4);Draw(); + if (CurrentTh) std::cerr << " In Triangles" << CurrentTh->Number(this) + << std::endl; + std::cerr << *ns[0] << *ns[1] << " " << *ns[2] << " " << std::endl; + std::cerr << " Bug in triangle " << this + << ":" << det << " != " << det2 << std::endl; + MeshError(5); + } + checka(0); + checka(1); + checka(2); +// if (ns[0]) assert( ns[0] - Meshbegin >= 0 ); +// if (ns[0]) assert( Meshend - ns[0] >= 0 ); +// if (ns[1]) assert( ns[1] - Meshbegin >= 0 ); +// if (ns[1]) assert( Meshend - ns[1] >= 0 ); +// if (ns[2]) assert( ns[2] - Meshbegin >= 0 ); +// if (ns[2]) assert( Meshend - ns[2] >= 0 ); + assert(ns[0] != ns[2]); + assert(ns[1] != ns[2]); + assert(ns[0] != ns[1]); +} + + +#endif + + + + +#ifdef DRAWING +extern Real4 xGrafCoef,yGrafCoef,xGrafOffSet,yGrafOffSet; // R2 -> I2 transform +extern R2 Gpmin,Gpmax; +//extern Real8 Gh; +// cf routine ILineTo IMoveto + +extern void IMoveTo(long i,long j); +extern void ILineTo(long i,long j); +extern char Getxyc(long &i,long &j); +extern void Draw(float ,float ); +extern void Draw(long ,long ); +extern void DrawMark(R2 r); +//inline void DrawMark(D2 r) {DrawMark(R2(r.x,r.y));} +inline void Move(I2 x) {IMoveTo(x.x,x.y);} +inline void Move(R2 x) {rmoveto(x.x,x.y);} +//inline void Move(D2 x) {rmoveto(x.x,x.y);} +inline void Line(I2 x){ILineTo(x.x,x.y);} +inline void Line(R2 x) {rlineto(x.x,x.y);} +//inline void Line(D2 x) {rlineto(x.x,x.y);} +#endif + +} + + +#endif + Index: /issm/trunk/src/c/Bamgx/MeshDraw.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/MeshDraw.cpp (revision 2740) +++ /issm/trunk/src/c/Bamgx/MeshDraw.cpp (revision 2740) @@ -0,0 +1,808 @@ +// -*- Mode : c++ -*- +// +// SUMMARY : +// USAGE : +// ORG : +// AUTHOR : Frederic Hecht +// E-MAIL : hecht@ann.jussieu.fr +// + +/* + + This file is part of Freefem++ + + Freefem++ is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + Freefem++ is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with Freefem++; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + */ +#ifdef DRAWING + +#include +#include +#include +#include +#include +using namespace std; + +#include "Mesh2.h" +#include "QuadTree.h" +#include "SetOfE4.h" +extern bool withrgraphique; + + + +namespace bamg { +Real4 xGrafCoef,yGrafCoef,xGrafOffSet,yGrafOffSet; +R2 GrafPMin,GrafPMax; +Real8 Grafh=0; + + + +void IMoveTo(long i,long j) +{ if (!withrgraphique) return; + rmoveto(float(i)/xGrafCoef+xGrafOffSet ,float(j)/yGrafCoef+yGrafOffSet );} + +void ILineTo(long i,long j) +{ if (!withrgraphique) return; + rlineto(float(i)/xGrafCoef+xGrafOffSet ,float(j)/yGrafCoef+yGrafOffSet );} + + +void Triangles::Draw( ) const { + if (!withrgraphique) return; + showgraphic(); +// if (!init) InitDraw(); + Int4 i; +// for (i=0;ir.x,v[1]->r.x), + Max(v[0]->r.y,v[1]->r.y), + Min(v[0]->r.x,v[1]->r.x), + Max(v[0]->r.y,v[1]->r.y))) { + v[0]->MoveTo(); + v[1]->LineTo(); + R2 M= ((R2) *v[0] + (R2) * v[1])*0.5; + Move(M); + char VertexDraw_i10[20]; + if (i<0) + sprintf(VertexDraw_i10,"%p",this); + else + sprintf(VertexDraw_i10,"%ld",i); + if (i>=0) + plotstring(&VertexDraw_i10[0]); + + } +} +void Vertex::Draw(Int4 i) const +{ + if (!withrgraphique) return; + if (CurrentTh && i<0&& i != -2 ) + { + if (CurrentTh->vertices <= this && this < CurrentTh->vertices + CurrentTh->nbv) + i = CurrentTh->Number(this); + } + if (InPtScreen(r.x,r.y)) { + char VertexDraw_i10[20]; + if (i<0) + sprintf(VertexDraw_i10,"%p",this); + else + sprintf(VertexDraw_i10,"%ld",i); + + showgraphic(); + // float eps = (MaxICoor/yGrafCoef)/100; + DrawMark(r); + + if (i>=0) + plotstring(&VertexDraw_i10[0]); } +} + +void DrawMark(R2 r) +{ + if (!withrgraphique) return; + if (InPtScreen(r.x,r.y)) { + float eps = Grafh/100; + rmoveto(r.x+eps,r.y); + rlineto(r.x,r.y+eps); + rlineto(r.x-eps,r.y); + rlineto(r.x,r.y-eps); + rlineto(r.x+eps,r.y);} + +} +void Triangle::Draw(Int4 i ) const +{ + if (!withrgraphique) return; + // int cc=LaCouleur(); + if (CurrentTh && i<0 && i != -2) + { + if (CurrentTh->triangles <= this && this < CurrentTh->triangles + CurrentTh->nbt) + i = CurrentTh->Number(this); + } + char i10[20]; + if (i<0) sprintf(i10,"%p",this); + else sprintf(i10,"%ld",i); + showgraphic(); + + if (ns[0] == 0) { + if (InRecScreen(Min(ns[1]->r.x,ns[2]->r.x), + Min(ns[1]->r.y,ns[2]->r.y), + Max(ns[1]->r.x,ns[2]->r.x), + Max(ns[1]->r.y,ns[2]->r.y))) + { rmoveto(ns[1]->r.x,ns[1]->r.y); + rlineto(ns[2]->r.x,ns[2]->r.y); + rmoveto( (ns[1]->r.x + ns[2]->r.x)/2.0 + (ns[1]->r.y - ns[2]->r.y)*0.1, + (ns[1]->r.y + ns[2]->r.y)/2.0 - (ns[1]->r.x - ns[2]->r.x)*0.1 ); + if(i>=0) plotstring(&i10[0]); + }} + else + if (InRecScreen(Min3(ns[0]->r.x,ns[1]->r.x,ns[2]->r.x), + Min3(ns[0]->r.y,ns[1]->r.y,ns[2]->r.y), + Max3(ns[0]->r.x,ns[1]->r.x,ns[2]->r.x), + Max3(ns[0]->r.y,ns[1]->r.y,ns[2]->r.y))) { + { + const int i0=0 , j01=EdgesVertexTriangle[i0][1];// j01==2 + const int i1=VerticesOfTriangularEdge[j01][1], j12=EdgesVertexTriangle[i1][1];//i1=1,j12=0 + const int i2=VerticesOfTriangularEdge[j12][1], j20=EdgesVertexTriangle[i2][1];//i1=2,j20=1 + + rmoveto(ns[i0]->r.x,ns[i0]->r.y); + if(Hidden(j01)) rmoveto(ns[i1]->r.x,ns[i1]->r.y); else rlineto(ns[i1]->r.x,ns[i1]->r.y); + if(Hidden(j12)) rmoveto(ns[i2]->r.x,ns[i2]->r.y); else rlineto(ns[i2]->r.x,ns[i2]->r.y); + if(Hidden(j20)) rmoveto(ns[i0]->r.x,ns[i0]->r.y); else rlineto(ns[i0]->r.x,ns[i0]->r.y); + + rmoveto( (ns[0]->r.x + ns[1]->r.x + ns[2]->r.x)/3.0 , + (ns[0]->r.y + ns[1]->r.y + ns[2]->r.y)/3.0); + + } + + if(i>=0) plotstring(&i10[0]); + } + // LaCouleur(cc); +} + +void Triangles::InitDraw() const +{ + if (!withrgraphique) return; + couleur(1); + GrafPMin = vertices[0].r; + GrafPMax = vertices[0].r; + // recherche des extrema des vertices GrafPMin,GrafPMax + Int4 i; + for (i=0;i2) + cout << "inquire: Nb de Triangle reel = " < draw metric o=> optim around (debug)",i++); + Show("t: find triangle contening the point with brute force",i++); + Show("e: find the nearest edge of triangle contening the point",i++); + Show("C: construct the in-circle in anistrope way",i++); + Show("n: show normal for find ",i++); + } + if (c=='p') openPS(0),PS=1; + if (c=='P') closePS(),PS=0; + if (c=='B') BTh.inquire(); + if ( c=='D') { + for(int iter=0;iter < 50;iter++) + { int k = 0; + for (Int4 icount=0; icount=0) + { couleur(2+(reft[i])%6); + triangles[i].Draw(reft[i]); + couleur(1); + } + } + + if (c=='b') { + Int4 i; + reffecran(); + for (i=0;iDraw(Number(tb)); + penthickness(1); + if(tb->link==0) { + double a,b; + TriangleAdjacent ta=CloseBoundaryEdgeV2(I,tb,a,b); + R2 A = *ta.EdgeVertex(0); + R2 B = *ta.EdgeVertex(1); + //Triangle * tt = ta; + // tt->Draw(Number(tt)); + penthickness(5); + // ta.EdgeVertex(0)->MoveTo(); + // ta.EdgeVertex(1)->LineTo(); + DrawMark(A*a+B*b); + penthickness(1); + } + } + if (c=='v' || c=='o' || c=='s' ) + { + couleur(3); + if(!setv) ReMakeTriangleContainingTheVertex(); + setv=1; + + R2 XY(x,y),P; + Real8 d; + Int4 i,j=0; + P = XY - vertices[0].r; + d = (P,P); + for (i=0;iMoveTo(); + tta.EdgeVertex(1)->LineTo(); + cout << " Edge " << Number(tta.EdgeVertex(0)) << " " << Number(tta.EdgeVertex(1)) << endl; + for (Int4 k=0;kr.x,s1->r.y); + D2 SA(sa->r.x,sa->r.y); + D2 SB(sb->r.x,sb->r.y); + D2 AB= SB-SA; + D2 MAB2=SB + SA; + D2 M1A=(S1 + SA)*0.5; + D2 MAB(MAB2.x*0.5,MAB2.y*0.5); + D2 A1=S1-SA; + D2 D = S1 - SB ; + { + Metric M=s1->m; + D2 ABo = M.Orthogonal(AB); + D2 A1o = M.Orthogonal(A1); + penthickness(1); + Move(MAB); + Line(MAB+ABo); + Line(MAB-ABo); + Move(M1A); + Line(M1A+A1o); + Line(M1A-A1o); + penthickness(3); + + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + cout << " d = " << d << " dd= " << dd << endl; + if (Abs(d) > dd*1.e-10) { + D2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + cout << C << s1->r <r <r << endl; + DrawMark(C); + cout << M << " l = "<< M(C-S1) << " lnew = 1 == " ; + //M.Draw(R2(C.x,C.y)); + M=M/M(C-S1); + cout << M(C-S1) << M << endl; + M.Draw(R2(C.x,C.y)); + } + + } + } + break; + } + cout << " fin recherche triangle " << endl; + couleur(1); + } + if (c=='='||c==249) { + rd=1; + float hx = (GrafPMax.x-GrafPMin.x); + float hy = (GrafPMax.y-GrafPMin.y); + Grafh = Max(hx,hy)*0.55; + cadreortho((GrafPMin.x+GrafPMax.x)*0.5,(GrafPMin.y+GrafPMax.y)*0.55,Grafh); + } + penthickness(1); + if (rd) + reffecran(),Draw(); + + } + // cout << endl; + delete [] reft; + CurrentTh = OldCurrentTh; + } +void Draw(float x,float y) +{ + if (!withrgraphique) return; + + if (InPtScreen(x,y)) { + float eps = Max(GrafPMax.x-GrafPMin.x,GrafPMax.y-GrafPMin.y)/100; + rmoveto(x+eps,y); + rlineto(x,y+eps); + rlineto(x-eps,y); + rlineto(x,y-eps); + rlineto(x+eps,y);} +} +char Getxyc(long &i,long &j) +{ + if (!withrgraphique) return 0; + + float x,y; + char c=::Getxyc(x,y); + i = (long)( (x - xGrafOffSet)*xGrafCoef); + j = (long)((y - yGrafOffSet)*yGrafCoef); + return c; +} + + +void Draw(long i,long j) +{ + if (!withrgraphique) return; + Draw(((float) i)/xGrafCoef+xGrafOffSet ,((float) j)/yGrafCoef+yGrafOffSet ); +} + +int Triangle::swapDRAW(Int2 a){ + int NbUnSwap=0; + if(a/4 !=0) {cout << a << "arete lock"<< a <ns[VerticesOfTriangularEdge[a1][0]]; + register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; + register Vertex *s1=t1->ns[OppositeVertex[a1]]; + register Vertex *s2=t2->ns[OppositeVertex[a2]]; + +#ifdef DEBUG + assert ( a >= 0 && a < 3 ); +#endif + + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + Icoor2 detA = Abs(det1) + Abs(det2); + Icoor2 detMin = Min(det1,det2); + + int OnSwap = 0; + // si 2 triangle infini (bord) => detT = -2; + if (sa == 0) {// les deux triangles sont frontieres + det2=bamg::det(s2->i,sb->i,s1->i); + OnSwap = det2 >0;} + else if (sb == 0) { // les deux triangles sont frontieres + det1=bamg::det(s1->i,sa->i,s2->i); + OnSwap = det1 >0;} + else if(( s1 != 0) && (s2 != 0) ) { + det1 = bamg::det(s1->i,sa->i,s2->i); + det2 = detT - det1; + OnSwap = (Abs(det1) + Abs(det2)) < detA; + Icoor2 detMinNew=Min(det1,det2); + cout << " detMin = " << detMin << " detMinNew " << detMinNew << endl; + if (! OnSwap &&(detMinNew>0)) { + OnSwap = detMin ==0; + if (! OnSwap) { + int kopt = 1; + while (1) + if(kopt) { + // critere de Delaunay pure + register Real8 xb1 = sb->i.x - s1->i.x, + x21 = s2->i.x - s1->i.x, + yb1 = sb->i.y - s1->i.y, + y21 = s2->i.y - s1->i.y, + xba = sb->i.x - sa->i.x, + x2a = s2->i.x - sa->i.x, + yba = sb->i.y - sa->i.y, + y2a = s2->i.y - sa->i.y, + cosb12 = xb1*x21 + yb1*y21 , + cosba2 = xba*x2a + yba*y2a , + sinb12 = det2, + sinba2 = t2->det, + //zsinba2 = t2->det; + + // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) + OnSwap = (cosb12 * sinba2) < (cosba2 * sinb12); + if(CurrentTh) + cout << "swap s1=" << CurrentTh->Number(sa) << " s2=" << CurrentTh->Number(sb) + << " t1= " << CurrentTh->Number(t1) << " t2=" << CurrentTh->Number(t2) << " " ; + + + cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12 + << " Onswap = " << OnSwap << endl; + + break; + } + else + { + // critere de Delaunay pure + Real8 som; + I2 AB=(I2) *sb - (I2) *sa; + I2 MAB2=((I2) *sb + (I2) *sa); + D2 MAB(MAB2.x*0.5,MAB2.y*0.5); + I2 A1=(I2) *s1 - (I2) *sa; + I2 D = (I2) * s1 - (I2) * sb ; + D2 S2(s2->i.x,s2->i.y); + D2 S1(s1->i.x,s1->i.y); + D2 SA(sa->i.x,sa->i.y); + D2 SB(sb->i.x,sb->i.y); + DrawMark(s2->r); + DrawMark(s1->r); + DrawMark(sa->r); + DrawMark(sb->r); + { + Metric M=s1->m; + D2 ABo = M.Orthogonal(SB-SA); + D2 A1o = M.Orthogonal(S1-SA); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if (Abs(d) > dd*1.e-3) { + D2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + cout << "M1 r2 =" << M(C - S2) << " r1 = " << M(C - S1) + << "ra = " << M(C - SA) << " rb = " << M(C-SB) ; + som = M(C - S2)/M(C - S1); + cout << " r1/r2 = " << som << endl; + } else + {kopt=1;continue;} + + } + { + Metric M=s2->m; + D2 ABo = M.Orthogonal(SB-SA); + D2 A1o = M.Orthogonal(S1-SA); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + cout << " d = " << Abs(d) << " dd " << dd << endl; + if(Abs(d) > dd*1.e-3) { + D2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + cout << "M2 r1 =" << M(C - S2) << " r2 = " << M(C - S1) + << "ra = " << M(C - SA) << " rb = " << M(C-SB) + << " r1/r2 = " << M(C - S2)/M(C - S1) << endl; + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + { + Metric M=sa->m; + D2 ABo = M.Orthogonal(SB-SA); + D2 A1o = M.Orthogonal(S1-SA); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + cout << " d = " << Abs(d) << " dd " << dd << endl; + if(Abs(d) > dd*1.e-3) { + D2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + cout << "M2 r1 =" << M(C - S2) << " r2 = " << M(C - S1) + << "ra = " << M(C - SA) << " rb = " << M(C-SB) + << " r1/r2 = " << M(C - S2)/M(C - S1) << endl; + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + { + Metric M=sb->m; + D2 ABo = M.Orthogonal(SB-SA); + D2 A1o = M.Orthogonal(S1-SA); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + cout << " d = " << Abs(d) << " dd " << dd << endl; + if(Abs(d) > dd*1.e-3) { + D2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + cout << "M2 r1 =" << M(C - S2) << " r2 = " << M(C - S1) + << "ra = " << M(C - SA) << " rb = " << M(C-SB) + << " r1/r2 = " << M(C - S2)/M(C - S1) << endl; + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + cout << som << endl; + OnSwap = som < 4; + break; + } + + } // OnSwap + } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) + } + cout << OnSwap << endl; + if( OnSwap ) { + couleur(0); + t1->Draw(); + t2->Draw(); + bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); + couleur(1); + t1->Draw(); + t2->Draw(); + } + else { + NbUnSwap ++; + t1->SetMarkUnSwap(a1); + } + return OnSwap; +} + + + + +void GeometricalEdge::Draw(Int4 i) +{ + if (!withrgraphique) return; + + if (CurrentTh && i<0 && i != -2) + { + if (CurrentTh->Gh.edges <= this && this < CurrentTh->Gh.edges + CurrentTh->Gh.nbe) + i = CurrentTh->Gh.Number((this)); + } + + v[0]->MoveTo(); + R2 x,x50; + int k=0,k50=0; + for (int ii = 0;ii<100;ii++) { + x= F( Real4(ii)/100.0); + if(ii==50) x50=x,k50=1; + if (InPtScreen(x.x,x.y) ) + if(k) rlineto(x.x,x.y); + else k=1,rmoveto(x.x,x.y); + } + if (InPtScreen(v[1]->r.x,v[1]->r.y)) + v[1]->LineTo(); + + char VertexDraw_i10[20]; + if( k50) { + if (i<0) + sprintf(VertexDraw_i10,"Eg%p",this); + else + sprintf(VertexDraw_i10,"Eg%ld",i); + rmoveto(x50.x,x50.y); + plotstring(&VertexDraw_i10[0]); } + +} + +void Geometry::InitDraw() const +{ + GrafPMin = vertices[0].r; + GrafPMax = vertices[0].r; + // recherche des extrema des vertices GrafPMin,GrafPMax + Int4 i; + for (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) < V0Ag.V0V1 > 0 || V0Bg.V0V1 >0 + // => V1Ag.V0V1 < 0 || V0Bg.V0V1 <0 + R2 V0V1 = V1-V0; + Real8 cos0A = ((Ag-V0),V0V1); + Real8 cos0B = ((Bg-V0),V0V1); + Real8 cos1A = ((Ag-V1),V0V1); + Real8 cos1B = ((Bg-V1),V0V1); + if ( (cos0A < 0 && cos0B <0) || (cos1A> 0 && cos1B >0)) + { + cerr << " Bug on pointer edge [" << V0 << " , " << V1 << " ] " + << " on geometrical edge " << Number(on) << " = [" << Ag << " , " << Bg << " ] " << endl; + cerr << cos0A << "> 0 || " << cos0B << " > 0 and "; + cerr << cos1A << "< 0 || " << cos1B << " < 0 " << endl; + + exit (1); + } + +#endif + if(NbTry) cout << "bug: edge = " << v0.r << " -> " << 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/MeshQuad.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/MeshQuad.cpp (revision 2740) +++ /issm/trunk/src/c/Bamgx/MeshQuad.cpp (revision 2740) @@ -0,0 +1,947 @@ +// -*- Mode : c++ -*- +// +// SUMMARY : +// USAGE : +// ORG : +// AUTHOR : Frederic Hecht +// E-MAIL : hecht@ann.jussieu.fr +// +// ********** DO NOT REMOVE THIS BANNER ********** +/* + + This file is part of Freefem++ + + Freefem++ is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. + + Freefem++ is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with Freefem++; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + */ + +// +// SUMMARY: Bamg: Bidimensional Anisotrope Mesh Generator +// RELEASE: 0 +// AUTHOR: F. Hecht, +// ORG : UMPC +// E-MAIL : Frederic.Hecht@Inria.fr +// +// ORIG-DATE: frev 98 +//--------------------------------------------------------- +// to make quad +// ------------------- + +#include +#include +#include +#include + +#include "Meshio.h" +#include "Mesh2.h" +#include "QuadTree.h" +#include "SetOfE4.h" + +namespace bamg { + +static const Direction NoDirOfSearch=Direction(); + +#ifdef __MWERKS__ +#pragma global_optimizer on +#pragma optimization_level 1 +//#pragma inline_depth 0 +#endif + +class DoubleAndInt4 { +public: double q; Int4 i3j; + int operator<(DoubleAndInt4 a){return q > a.q;} +};// to sort by decroissant + + + +template inline void HeapSort(T *c,long n) +{ + long l,j,r,i; + T crit; + c--; // on decale de 1 pour que le tableau commence a 1 + if( n <= 1) return; + l = n/2 + 1; + r = n; + while (1) { // label 2 + if(l <= 1 ) { // label 20 + crit = c[r]; + c[r--] = c[1]; + if ( r == 1 ) { c[1]=crit; return;} + } else crit = c[--l]; + j=l; + while (1) {// label 4 + i=j; + j=2*j; + if (j>r) {c[i]=crit;break;} // L8 -> G2 + if ((j G2 + } + } +} + +Triangle * swapTest(Triangle *t1,Int2 a); +// + + +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)); + // cout << A << B << C << D ; + // cout << " sinmin " << sinmin << " " << cosDAB << " " << cosABC << " " << cosBCD << " " << cosCDA << endl; + // rattente(1); + if (sinmin<=0) return sinmin; + return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA))); +} + + GeometricalEdge * Triangles::ProjectOnCurve( Edge & BhAB, Vertex & vA, Vertex & vB, + Real8 theta, + Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR) + +{ + void *pA=0,*pB=0; + Real8 tA=0,tB=0; + R2 A=vA,B=vB; + Vertex * pvA=&vA, * pvB=&vB; + if (vA.vint == IsVertexOnVertex) + { + // cout << " Debut vertex = " << BTh.Number(vA.onbv) ; + pA=vA.onbv; + } + else if (vA.vint == IsVertexOnEdge) + { + pA=vA.onbe->be; + tA=vA.onbe->abcisse; + // cout << " Debut edge = " << BTh.Number(vA.onbv) << " " << tA ; + + } + else + {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vA) << " " << endl; + cerr << " forget call to SetVertexFieldOnBTh" << endl; + MeshError(-1); + } + + if (vB.vint == IsVertexOnVertex) + { + // cout << " Fin vertex = " << BTh.Number(vB.onbv) << endl; + pB=vB.onbv; + } + else if(vB.vint == IsVertexOnEdge) + { + pB=vB.onbe->be; + tB=vB.onbe->abcisse; + // cout << " Fin edge = " << BTh.Number(vB.onbe->be) << " " << tB ; + + } + else + {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vB) << " " << endl; + cerr << " forget call to SetVertexFieldOnBTh" << endl; + MeshError(-1); + } + Edge * e = &BhAB; + assert( pA && pB && e); + // be carefull the back ground edge e is on same geom edge + // of the initiale edge def by the 2 vertex A B; + assert(e>=BTh.edges && e !IsRequiredVertex + if (vA.vint == IsVertexOnEdge) + if (tA<=0) + assert(! (*vA.onbe->be)[0].on->IsRequiredVertex()); + else if (tA>=1) + assert(!(*vA.onbe->be)[1].on->IsRequiredVertex()); +#endif + + if( vA.vint == IsVertexOnEdge) + { // find the start edge + e = vA.onbe->be; + + } + else if (vB.vint == IsVertexOnEdge) + { + theta = 1-theta; + Exchange(tA,tB); + Exchange(pA,pB); + Exchange(pvA,pvB); + Exchange(A,B); + e = vB.onbe->be; + + // cout << " EXCHANGE A et B) " << endl; + } + else + { // do the search by walking + assert(0 /* A FAIRE */); + } + + // find the direction of walking with sens of edge and pA,PB; + R2 AB=B-A; + + Real8 cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); + int kkk=0; + int sens = (cosE01AB>0) ? 1 : 0; + + // Real8 l=0; // length of the edge AB + Real8 abscisse = -1; + + for (int cas=0;cas<2;cas++) + {// 2 times algo: + // 1 for computing the length l + // 2 for find the vertex + int iii; + Vertex *v0=pvA,*v1; + Edge *neee,*eee; + Real8 lg =0; // length of the curve + Real8 te0; + // we suppose take the curve's abcisse + // cout << kkk << " e = " << BTh.Number(e) << " v0= " + // << BTh.Number(v0) << " v1 = " << BTh.Number((*e)[sens]) << endl; + for ( eee=e,iii=sens,te0=tA; + eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; + neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) + { + // cout << kkk << " eee = " << BTh.Number(eee) << " v0= " + // << BTh.Number(v0) << " v1 = " << BTh.Number(v1) << endl; + + assert(kkk++<100); + assert(eee); + Real8 lg0 = lg; + Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); + lg += dp; + if (cas && abscisse <= lg) + { // ok we find the geom edge + Real8 sss = (abscisse-lg0)/dp; + Real8 thetab = te0*(1-sss)+ sss*iii; + assert(thetab>=0 && thetab<=1); + BR = VertexOnEdge(&R,eee,thetab); + + // cout << Number(R) << " = " << thetab << " on " << BTh.Number(eee) + // << " = " << R << endl; + + return Gh.ProjectOnCurve(*eee,thetab,R,GR); + + } + } + // we find the end + if (v1 != pvB) + { + if (( void*) v1 == pB) + tB = iii; + + Real8 lg0 = lg; + assert(eee); + v1 = pvB; + Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); + lg += dp; + abscisse = lg*theta; + if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end + { // ok we find the geom edge + Real8 sss = (abscisse-lg0)/dp; + Real8 thetab = te0*(1-sss)+ sss*tB; + assert(thetab>=0 && thetab<=1); + BR = VertexOnEdge(&R,eee,thetab); + + // cout << kkk << " eee = " << BTh.Number(eee) << " v0= " + // << BTh.Number(v0) << " " << te0 + // << " v1 = " << BTh.Number(v1) << " " << tB << endl; + + //out << Number(R) << " Opt = " << thetab << " on " << BTh.Number(eee) + // << " = " << R << endl; + + return Gh.ProjectOnCurve(*eee,thetab,R,GR); + } + } + abscisse = lg*theta; + + } + cerr << " Big Bug" << endl; + MeshError(678); + return 0; // just for the compiler + +} + + +void Triangles::MakeQuadrangles(double costheta) +{ + long int verbosity=0; + + if (verbosity>2) + cout << " -- MakeQuadrangles costheta = " << costheta << endl; + if (verbosity>5) + cout << " (in) Nb of Quadrilaterals = " << NbOfQuad + << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 + << " Nb of outside triangles = " << NbOutT << endl; + + if (costheta >1) { + if (verbosity>5) + cout << " do nothing costheta >1 "<< endl; + return;} + + + Int4 nbqq = (nbt*3)/2; + DoubleAndInt4 *qq = new DoubleAndInt4[nbqq]; + + Int4 i,ij; + int j; + Int4 k=0; + for (i=0;i=costheta) + qq[k++].i3j=i*3+j; +// sort qq + HeapSort(qq,k); + + Int4 kk=0; + for (ij=0;ij=costheta) + triangles[i].SetHidden(j),kk++; + } + NbOfQuad = kk; + if (verbosity>2) + { + cout << " (out) Nb of Quadrilaterals = " << NbOfQuad + << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 + << " Nb of outside triangles = " << NbOutT << endl; + } + delete [] qq; +#ifdef DRAWING2 + Draw(); + inquire(); +#endif + +} +/* +Triangles::BThBoundary(Edge e,Real 8) const +{ + // pointeur of the background must be on + // + Edge be = e.on; +} +*/ +int Triangles::SplitElement(int choice) +{ + long int verbosity=0; + + Direction NoDirOfSearch; + const int withBackground = &BTh != this && &BTh; + if (verbosity>2) + cout << " -- SplitElement " << (choice? " Q->4Q and T->4T " : " Q->4Q or T->3Q " ) << endl;; + if (verbosity>5) + cout << endl << " (in) Nb of Quadrilaterals = " << NbOfQuad + << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 + << " Nb of outside triangles = " << NbOutT << endl; + + ReNumberingTheTriangleBySubDomain(); +#ifdef DRAWING2 + Draw(); + inquire(); +#endif + //int nswap =0; + const Int4 nfortria( choice ? 4 : 6); + if(withBackground) + { + BTh.SetVertexFieldOn(); + SetVertexFieldOnBTh(); + } + else + BTh.SetVertexFieldOn(); + + Int4 newnbt=0,newnbv=0; + Int4 * kedge = 0; + Int4 newNbOfQuad=NbOfQuad; + Int4 * ksplit= 0, * ksplitarray=0; + Int4 kkk=0; + int ret =0; + if (nbvxfindtrie(Number(v0),Number(v1)); + // cout << ":" << ke << "," << !!t.link << " " << &tt ; + if (ke<0) // new + { + if (&tt) // internal triangles all the boundary + { // new internal edges + Int4 ii = Number(tt); + int jj = ta; + + kedge[3*i+j]=k;// save the vertex number + kedge[3*ii+jj]=k; + if (k=40); + } + // now do the element split + newNbOfQuad = 4*NbOfQuad; + nbv = k; +#ifdef DRAWING2 + inquire(); +#endif +// cout << " Nbv = " << nbv << endl; + kkk = nbt; + ksplit[-1] = nbt; + // look on old true triangles + + for (i=0;i0); + + // def the numbering k (edge) i vertex + int k0 = ke; + int k1 = NextEdge[k0]; + int k2 = PreviousEdge[k0]; + int i0 = OppositeVertex[k0]; + int i1 = OppositeVertex[k1]; + int i2 = OppositeVertex[k2]; + + Triangle &t0=triangles[i]; + Vertex * v0=t0(i0); + Vertex * v1=t0(i1); + Vertex * v2=t0(i2); + + // cout << "nbmkadj " << nbmkadj << " it=" << i <