Changeset 5120
- Timestamp:
- 08/10/10 13:39:21 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 21 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Makefile.am
r5114 r5120 69 69 ./objects/Bamg/Mesh.cpp\ 70 70 ./objects/Bamg/Mesh.h\ 71 ./objects/Bamg/ MeshVertex.cpp\72 ./objects/Bamg/ MeshVertex.h\71 ./objects/Bamg/BamgVertex.cpp\ 72 ./objects/Bamg/BamgVertex.h\ 73 73 ./objects/Bamg/VertexOnEdge.h\ 74 74 ./objects/Bamg/VertexOnEdge.cpp\ … … 605 605 ./objects/Bamg/Mesh.h\ 606 606 ./objects/Bamg/Mesh.cpp\ 607 ./objects/Bamg/ MeshVertex.cpp\608 ./objects/Bamg/ MeshVertex.h\607 ./objects/Bamg/BamgVertex.cpp\ 608 ./objects/Bamg/BamgVertex.h\ 609 609 ./objects/Bamg/VertexOnEdge.h\ 610 610 ./objects/Bamg/VertexOnEdge.cpp\ -
issm/trunk/src/c/objects/Bamg/BamgVertex.cpp
r5100 r5120 9 9 10 10 /*Methods*/ 11 /*FUNCTION MeshVertex::Smoothing{{{1*/12 double MeshVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){11 /*FUNCTION BamgVertex::Smoothing{{{1*/ 12 double BamgVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){ 13 13 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/ 14 14 15 register MeshVertex* s=this;16 MeshVertex &vP = *s,vPsave=vP;15 register BamgVertex* s=this; 16 BamgVertex &vP = *s,vPsave=vP; 17 17 18 18 register Triangle* tbegin= t , *tria = t , *ttc; … … 24 24 25 25 if (!tria->Hidden(j)){ 26 MeshVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];26 BamgVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]]; 27 27 28 28 R2 Q = vQ,QP(P-Q); … … 70 70 vP.i = Th.toI2(PNew); 71 71 72 MeshVertex vPnew = vP;72 BamgVertex vPnew = vP; 73 73 74 74 int ok=1; … … 82 82 tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]); 83 83 if (loop) { 84 MeshVertex *v0,*v1,*v2,*v3;84 BamgVertex *v0,*v1,*v2,*v3; 85 85 if (tria->det<0) ok =1; 86 86 else if (tria->Quadrangle(v0,v1,v2,v3)) … … 113 113 } 114 114 /*}}}1*/ 115 /*FUNCTION MeshVertex::MetricFromHessian{{{1*/116 void MeshVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){115 /*FUNCTION BamgVertex::MetricFromHessian{{{1*/ 116 void BamgVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){ 117 117 /*Compute Metric from Hessian*/ 118 118 … … 192 192 } 193 193 /*}}}1*/ 194 /*FUNCTION MeshVertex::Echo {{{1*/195 196 void MeshVertex::Echo(void){194 /*FUNCTION BamgVertex::Echo {{{1*/ 195 196 void BamgVertex::Echo(void){ 197 197 198 198 printf("Vertex:\n"); … … 205 205 } 206 206 /*}}}*/ 207 /*FUNCTION MeshVertex::Optim {{{1*/208 long MeshVertex::Optim(int i,int koption){207 /*FUNCTION BamgVertex::Optim {{{1*/ 208 long BamgVertex::Optim(int i,int koption){ 209 209 long ret=0; 210 210 if ( t && (vint >= 0 ) && (vint <3) ){ … … 221 221 /*Intermediary*/ 222 222 /*FUNCTION QuadQuality{{{1*/ 223 double QuadQuality(const MeshVertex & a,const MeshVertex &b,const MeshVertex &c,const MeshVertex &d) {223 double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) { 224 224 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/ 225 225 -
issm/trunk/src/c/objects/Bamg/BamgVertex.h
r5100 r5120 15 15 class VertexOnEdge; 16 16 17 class MeshVertex {17 class BamgVertex { 18 18 19 19 public: 20 I2 i; // integer coordinates 21 R2 r; // real coordinates 22 Metric m; 23 long ReferenceNumber; 20 21 I2 i; // integer coordinates 22 R2 r; // real coordinates 23 Metric m; 24 long ReferenceNumber; 24 25 Direction DirOfSearch; 25 short vint; // the vertex number in triangle; varies between 0 and 2 in t 26 short vint; // the vertex number in triangle; varies between 0 and 2 in t 27 26 28 union { 27 Triangle * t;// one triangle which is containing the vertex28 long color;29 MeshVertex* to; // used in geometry MeshVertex to know the Mesh MeshVertex associated30 VertexOnGeom * onGeometry;// if vint == 8; // set with Mesh::SetVertexFieldOn()31 MeshVertex* onBackgroundVertex;// if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()32 VertexOnEdge * onBackgroundEdge;// if vint == 32 on Background edge29 Triangle *t; // one triangle which is containing the vertex 30 long color; 31 BamgVertex *to; // used in geometry BamgVertex to know the Mesh BamgVertex associated 32 VertexOnGeom *onGeometry; // if vint == 8; // set with Mesh::SetVertexFieldOn() 33 BamgVertex *onBackgroundVertex; // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh() 34 VertexOnEdge *onBackgroundEdge; // if vint == 32 on Background edge 33 35 }; 34 36 … … 47 49 48 50 //inline functions 49 inline void Set(const MeshVertex &rec,const Mesh & ,Mesh & ){*this=rec;}51 inline void Set(const BamgVertex &rec,const Mesh & ,Mesh & ){*this=rec;} 50 52 }; 51 53 52 54 //Intermediary 53 double QuadQuality(const MeshVertex &,const MeshVertex &,const MeshVertex &,const MeshVertex &);55 double QuadQuality(const BamgVertex &,const BamgVertex &,const BamgVertex &,const BamgVertex &); 54 56 } 55 57 #endif -
issm/trunk/src/c/objects/Bamg/Edge.h
r5095 r5120 3 3 4 4 #include "./include.h" 5 #include "./ MeshVertex.h"5 #include "./BamgVertex.h" 6 6 #include "../../include/include.h" 7 7 #include "../../shared/Exceptions/exceptions.h" … … 16 16 17 17 public: 18 MeshVertex* v[2];18 BamgVertex* v[2]; 19 19 long ref; 20 20 GeometricalEdge* onGeometry; … … 22 22 23 23 //Operators 24 MeshVertex & operator[](int i){return *v[i];};25 MeshVertex * operator()(int i){return v[i];};24 BamgVertex & operator[](int i){return *v[i];}; 25 BamgVertex * operator()(int i){return v[i];}; 26 26 R2 operator()(double t) const; // return the point 27 const MeshVertex & operator[](int i) const { return *v[i];};27 const BamgVertex & operator[](int i) const { return *v[i];}; 28 28 29 29 //Methods 30 void ReNumbering( MeshVertex *vb,MeshVertex *ve, long *renu){30 void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){ 31 31 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb]; 32 32 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb]; -
issm/trunk/src/c/objects/Bamg/GeometricalVertex.h
r3913 r5120 3 3 4 4 #include "./include.h" 5 #include " MeshVertex.h"5 #include "BamgVertex.h" 6 6 7 7 namespace bamg { … … 9 9 class Geometry; 10 10 11 class GeometricalVertex : public MeshVertex {11 class GeometricalVertex : public BamgVertex { 12 12 13 13 public: -
issm/trunk/src/c/objects/Bamg/Geometry.cpp
r5095 r5120 13 13 14 14 /*Constructors/Destructors*/ 15 /*FUNCTION Geometry::Geometry(const Geometry & Gh){{{1*/ 15 Geometry::Geometry(){ 16 Init(); 17 } 18 Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){ 19 Init(); 20 ReadGeometry(bamggeom,bamgopts); 21 AfterRead(); 22 } 23 /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{1*/ 16 24 Geometry::Geometry(const Geometry & Gh) { 17 25 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/ … … 22 30 quadtree=0; 23 31 vertices = nbv ? new GeometricalVertex[nbv] : NULL; 24 triangles = nbt ? new Triangle[nbt]:NULL;25 32 edges = nbe ? new GeometricalEdge[nbe]:NULL; 26 33 curves= NbOfCurves ? new Curve[NbOfCurves]:NULL; … … 34 41 for (i=0;i<NbSubDomains;i++) 35 42 subdomains[i].Set(Gh.subdomains[i],Gh,*this); 36 if (nbt); {37 ISSMERROR("nbt");38 }39 43 } 40 44 /*}}}1*/ … … 42 46 Geometry::~Geometry() { 43 47 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/ 44 45 long int verbose=0; 46 47 if (NbRef>0){ 48 ISSMERROR("NbRef>0"); 49 } 50 if(vertices) delete [] vertices;vertices=0; 51 if(edges) delete [] edges;edges=0; 52 // if(edgescomponante) delete [] edgescomponante; edgescomponante=0; 53 if(triangles) delete [] triangles;triangles=0; 54 if(quadtree) delete quadtree;quadtree=0; 55 if(curves) delete []curves;curves=0;NbOfCurves=0; 48 ISSMASSERT(NbRef<=0); 49 if(vertices) delete [] vertices; vertices=0; 50 if(edges) delete [] edges; edges=0; 51 if(triangles) delete [] triangles; triangles=0; 52 if(quadtree) delete quadtree; quadtree=0; 53 if(curves) delete []curves; curves=0;NbOfCurves=0; 56 54 if(subdomains) delete [] subdomains;subdomains=0; 57 // if(ordre) delete [] ordre; 58 EmptyGeometry(); 55 Init(); 59 56 } 60 57 /*}}}1*/ … … 65 62 66 63 int verbose; 67 nb iv=nbv=nbvx=0;68 nbe= nbt=nbtx=0;64 nbv=0; 65 nbe=0; 69 66 NbOfCurves=0; 70 67 … … 72 69 int i,j,k,n,i1,i2; 73 70 74 //initialize some variables 75 int Version=1,dim=2; 76 nbv=bamggeom->VerticesSize[0]; 77 nbe=bamggeom->EdgesSize[0]; 78 nbvx = nbv; 79 nbiv = nbv; 71 /*initialize some variables*/ 80 72 verbose=bamgopts->verbose; 73 nbv = bamggeom->VerticesSize[0]; 74 nbe = bamggeom->EdgesSize[0]; 81 75 82 76 //some checks 83 if (bamggeom->Vertices==NULL){ 84 ISSMERROR("the domain provided does not contain any vertex"); 85 } 86 if (bamggeom->Edges==NULL){ 87 ISSMERROR("the domain provided does not contain any edge"); 88 } 77 if (bamggeom->Vertices==NULL) ISSMERROR("the domain provided does not contain any vertex"); 78 if (bamggeom->Edges==NULL) ISSMERROR("the domain provided does not contain any edge"); 89 79 90 80 //Vertices … … 92 82 if(verbose>5) printf(" processing Vertices\n"); 93 83 if (bamggeom->VerticesSize[1]!=3) ISSMERROR("Vertices should have 3 columns"); 94 vertices = new GeometricalVertex[nbv x];84 vertices = new GeometricalVertex[nbv]; 95 85 for (i=0;i<nbv;i++) { 96 86 vertices[i].r.x=(double)bamggeom->Vertices[i*3+0]; … … 101 91 vertices[i].Set(); 102 92 } 103 //find domain extrema for pmin,pmax that will define the square box 104 //used for by the quadtree 93 /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/ 105 94 pmin = vertices[0].r; 106 95 pmax = vertices[0].r; … … 111 100 pmax.y = Max(pmax.y,vertices[i].r.y); 112 101 } 113 R2 DD05 = (pmax-pmin)*0.05; 114 pmin -= DD05; 115 pmax += DD05; 116 //coefIcoor is the coefficient used to have coordinates 117 //in ]0 1[^2: x'=coefIcoor*(x-pmin.x) 118 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 119 if(coefIcoor <=0){ 120 ISSMERROR("coefIcoor should be positive"); 121 } 102 /*Offset pmin and pmax to avoid round-off errors*/ 103 R2 offset = (pmax-pmin)*0.05; 104 pmin -= offset; 105 pmax += offset; 106 /*coefIcoor is the coefficient used for integer coordinates: 107 * (x-pmin.x) 108 * Icoor x = (2^30 -1) ------------ 109 * D 110 * where D is the longest side of the domain (direction x or y) 111 * so that (x-pmin.x)/D is in ]0 1[ 112 */ 113 coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 114 if(coefIcoor<=0) ISSMERROR("coefIcoor should be positive"); 122 115 } 123 116 else{ 124 ISSMERROR("No MeshVertex provided");117 ISSMERROR("No BamgVertex provided"); 125 118 } 126 119 … … 427 420 double eps=1e-20; 428 421 QuadTree quadtree; // build quadtree to find duplicates 429 MeshVertex* v0=vertices;422 BamgVertex* v0=vertices; 430 423 GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0; 431 424 … … 446 439 447 440 //find nearest vertex already present in the quadtree (NULL if empty) 448 MeshVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);441 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 449 442 450 443 //if there is a vertex found that is to close to vertices[i] -> error … … 454 447 j=vg-v0g; 455 448 //check that the clostest vertex is not itself... 456 if ( v != &( MeshVertex &) vertices[j]){457 ISSMERROR(" v != &( MeshVertex &) vertices[j]");449 if ( v != &(BamgVertex &) vertices[j]){ 450 ISSMERROR(" v != &(BamgVertex &) vertices[j]"); 458 451 } 459 452 vertices[i].link = vertices + j; … … 748 741 GeometricalEdge* on =start,* pon=0; 749 742 // walk with the cos on geometry 750 int k=0;743 int counter=0; 751 744 while(pon != on){ 752 k++; 745 counter++; 746 ISSMASSERT(counter<100); 753 747 pon = on; 754 if (k>=100){755 ISSMERROR("k>=100");756 }757 748 R2 A= (*on)[0]; 758 749 R2 B= (*on)[1]; … … 775 766 776 767 printf("Geometry:\n"); 777 printf(" nbvx (maximum number of vertices) : %i\n",nbvx);778 printf(" nbtx (maximum number of triangles): %i\n",nbtx);779 768 printf(" nbv (number of vertices) : %i\n",nbv); 780 printf(" nbt (number of triangles): %i\n",nbt);781 769 printf(" nbe (number of edges) : %i\n",nbe); 782 printf(" nbv (number of initial vertices) : %i\n",nbiv);783 770 printf(" NbSubDomains: %i\n",NbSubDomains); 784 771 printf(" NbOfCurves: %i\n",NbOfCurves); … … 797 784 } 798 785 /*}}}*/ 799 /*FUNCTION Geometry:: EmptyGeometry(){{{1*/800 void Geometry:: EmptyGeometry(){786 /*FUNCTION Geometry::Init{{{1*/ 787 void Geometry::Init(void){ 801 788 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/ 802 789 803 790 NbRef=0; 791 nbv=0; 792 nbe=0; 804 793 quadtree=0; 805 794 curves=0; 806 // edgescomponante=0;807 795 triangles=0; 808 796 edges=0; 809 797 vertices=0; 810 798 NbSubDomains=0; 811 // nbtf=0;812 // BeginOfCurve=0;813 nbiv=nbv=nbvx=0;814 nbe=nbt=nbtx=0;815 799 NbOfCurves=0; 816 // BeginOfCurve=0;817 800 subdomains=0; 818 MaxCornerAngle = 10*Pi/180; 801 MaxCornerAngle = 10*Pi/180; //default is 10 degres 819 802 } 820 803 /*}}}1*/ 821 804 /*FUNCTION Geometry::ProjectOnCurve {{{1*/ 822 GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s, MeshVertex &V,VertexOnGeom &GV) const {805 GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const { 823 806 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/ 824 807 /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/ … … 839 822 840 823 //Get the two vertices of the edge 841 const MeshVertex &v0=e[0];842 const MeshVertex &v1=e[1];824 const BamgVertex &v0=e[0]; 825 const BamgVertex &v1=e[1]; 843 826 844 827 //Get position of V0, V1 and vector v0->v1 … … 906 889 907 890 if ((*eg0)(sens0)==(GeometricalVertex*)vg0) 908 vg0=VertexOnGeom(*( MeshVertex*) vg0,*eg0,sens0); //vg0 = absisce891 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,sens0); //vg0 = absisce 909 892 910 893 if ((*eg1)(sens1)==(GeometricalVertex*)vg1) 911 vg1=VertexOnGeom(*( MeshVertex*) vg1,*eg1,sens1);894 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,sens1); 912 895 913 896 double sg; -
issm/trunk/src/c/objects/Bamg/Geometry.h
r5091 r5120 19 19 20 20 public: 21 long NbRef; // counter of ref on the this class if 0 we can delete 22 long nbvx,nbtx; // maximum number of vertices 23 long nbv,nbt,nbiv,nbe; // number of vertices 24 long NbSubDomains; // 25 long NbOfCurves; 26 GeometricalVertex* vertices; 27 Triangle* triangles; 28 GeometricalEdge* edges; 29 QuadTree* quadtree; 30 GeometricalSubDomain* subdomains; 31 Curve* curves; 32 R2 pmin,pmax; // extrema 33 double coefIcoor; // coef to integer Icoor1; 34 double MaxCornerAngle; 21 22 long NbRef; // counter of ref on the this class if 0 we can delete 23 long nbv; // number of vertices 24 long nbe; // number of edges 25 long NbSubDomains; 26 long NbOfCurves; 27 GeometricalVertex *vertices; 28 Triangle *triangles; 29 GeometricalEdge *edges; 30 QuadTree *quadtree; 31 GeometricalSubDomain *subdomains; 32 Curve *curves; 33 R2 pmin,pmax; // domain extrema coordinates 34 double coefIcoor; // coef to integer Icoor1; 35 double MaxCornerAngle; 35 36 36 37 //Constructor/Destructors 37 38 ~Geometry(); 38 Geometry(const Geometry & Gh); //Copy Operator 39 Geometry(int nbg,const Geometry** ag); // intersection operator 39 Geometry(); 40 Geometry(const Geometry & Gh); 41 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts); 40 42 41 43 //Operators … … 46 48 47 49 //Methods 48 int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); }49 50 void Echo(); 50 51 I2 toI2(const R2 & P) const { 51 52 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) 52 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} 53 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) ); 54 } 53 55 double MinimalHmin() {return 2.0/coefIcoor;} 54 56 double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 55 57 void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); 56 void EmptyGeometry(); 57 Geometry() {EmptyGeometry();}// empty Geometry 58 void Init(void); 58 59 void AfterRead(); 59 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();}60 60 long Number(const GeometricalVertex & t) const { return &t - vertices;} 61 61 long Number(const GeometricalVertex * t) const { return t - vertices;} … … 64 64 long Number(const Curve * c) const { return c - curves;} 65 65 void UnMarkEdges() {for (int i=0;i<nbe;i++) edges[i].SetUnMark();} 66 GeometricalEdge * ProjectOnCurve(const Edge & ,double, MeshVertex &,VertexOnGeom &) const ;66 GeometricalEdge * ProjectOnCurve(const Edge & ,double,BamgVertex &,VertexOnGeom &) const ; 67 67 GeometricalEdge * Containing(const R2 P, GeometricalEdge * start) const; 68 68 void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); -
issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.cpp
r5095 r5120 44 44 double ba,bb; 45 45 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); 46 MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);46 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); 47 47 NewItem(A,Metric(ba,v0,bb,v1)); 48 48 t=edge; … … 50 50 if (det(v0.i,v1.i,b)>=0) { 51 51 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); 52 MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);52 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); 53 53 NewItem(A,Metric(ba,v0,bb,v1)); 54 54 return; … … 80 80 long int verbose=2; 81 81 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); 82 MeshVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);82 BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); 83 83 NewItem(A,Metric(ba,v0,bb,v1)); 84 84 return; … … 232 232 lIntTria[Size].x = x; 233 233 Metric m0,m1,m2; 234 register MeshVertex * v;234 register BamgVertex * v; 235 235 if ((v=(*tt)(0))) m0 = v->m; 236 236 if ((v=(*tt)(1))) m1 = v->m; … … 307 307 /*}}}1*/ 308 308 /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/ 309 long ListofIntersectionTriangles::NewPoints( MeshVertex* vertices,long &nbv,long nbvx){309 long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){ 310 310 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/ 311 311 … … 354 354 355 355 si += sint; 356 if ( nbv< nbvx) {356 if ( nbv<maxnbv) { 357 357 vertices[nbv].r = C; 358 358 vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m); -
issm/trunk/src/c/objects/Bamg/ListofIntersectionTriangles.h
r5095 r5120 71 71 void SplitEdge(const Mesh & ,const R2 &,const R2 &,int nbegin=0); 72 72 double Length(); 73 long NewPoints( MeshVertex *,long & nbv,long nbvx);73 long NewPoints(BamgVertex *,long & nbv,long maxnbv); 74 74 void NewSubSeg(GeometricalEdge *e,double s0,double s1){ 75 75 long int verbosity=0; -
issm/trunk/src/c/objects/Bamg/Mesh.cpp
r5095 r5120 16 16 17 17 /*Initialize fields*/ 18 PreInit(0);18 Init(0); 19 19 20 20 /*Read Geometry if provided*/ … … 45 45 Mesh::Mesh(double* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){ 46 46 47 PreInit(0);47 Init(0); 48 48 ReadMesh(index,x,y,nods,nels); 49 49 SetIntCoor(); … … 100 100 printf(" number of triangles %i, remove = %i\n",kt,nbInT-kt); 101 101 printf(" number of New boundary edge %i\n",nbNewBedge); 102 long i nbvx=k;103 PreInit(inbvx);102 long imaxnbv =k; 103 Init(imaxnbv); 104 104 for (i=0;i<Tho.nbv;i++) 105 105 if (kk[i]>=0) … … 110 110 nbv++; 111 111 } 112 if (i nbvx!= nbv){113 ISSMERROR("i nbvx!= nbv");112 if (imaxnbv != nbv){ 113 ISSMERROR("imaxnbv != nbv"); 114 114 } 115 115 for (i=0;i<Tho.nbt;i++) … … 154 154 } 155 155 /*}}}1*/ 156 /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx) COPY{{{1*/157 Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long nbvxx)156 /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{1*/ 157 Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) 158 158 : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) { 159 159 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/ 160 160 Gh.NbRef++; 161 nbvxx = Max(nbvxx,Th.nbv);161 maxnbv_in = Max(maxnbv_in,Th.nbv); 162 162 long i; 163 163 // do all the allocation to be sure all the pointer existe 164 164 165 PreInit(nbvxx);// to make the allocation165 Init(maxnbv_in);// to make the allocation 166 166 // copy of triangles 167 167 nt=Th.nt; 168 168 nbv = Th.nbv; 169 169 nbt = Th.nbt; 170 nbiv = Th.nbiv;171 170 nbe = Th.nbe; 172 171 NbSubDomains = Th.NbSubDomains; … … 245 244 else if (BTh.NbRef==0) delete &BTh; 246 245 } 247 PreInit(0); // set all to zero246 Init(0); // set all to zero 248 247 } 249 248 /*}}}1*/ … … 260 259 261 260 nbv=nods; 262 nbvx=nbv;261 maxnbv=nbv; 263 262 nbt=nels; 264 nbiv=nbvx;265 263 266 264 //Vertices 267 265 if (verbose) printf("Reading vertices (%i)\n",nbv); 268 vertices=( MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));269 ordre=( MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));266 vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex)); 267 ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*)); 270 268 for (i=0;i<nbv;i++){ 271 269 vertices[i].r.x=x[i]; … … 276 274 vertices[i].color=0; 277 275 } 278 nbtx=2* nbvx-2; // for filling The Holes and quadrilaterals276 nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals 279 277 280 278 //Triangles … … 311 309 312 310 nbv=bamgmesh->VerticesSize[0]; 313 nbvx=nbv;311 maxnbv=nbv; 314 312 nbt=bamgmesh->TrianglesSize[0]; 315 nbiv=nbvx;316 313 317 314 //Vertices … … 319 316 if(verbose>5) printf(" processing Vertices\n"); 320 317 321 vertices=( MeshVertex*)xmalloc(nbv*sizeof(MeshVertex));322 ordre=( MeshVertex**)xmalloc(nbv*sizeof(MeshVertex*));318 vertices=(BamgVertex*)xmalloc(nbv*sizeof(BamgVertex)); 319 ordre=(BamgVertex**)xmalloc(nbv*sizeof(BamgVertex*)); 323 320 324 321 for (i=0;i<nbv;i++){ … … 330 327 vertices[i].color=0; 331 328 } 332 nbtx=2* nbvx-2; // for filling The Holes and quadrilaterals329 nbtx=2*maxnbv-2; // for filling The Holes and quadrilaterals 333 330 } 334 331 else{ … … 450 447 for (i=0;i<nbe;i++){ 451 448 for (j=0;j<2;j++) { 452 MeshVertex *v=edges[i].v[j];449 BamgVertex *v=edges[i].v[j]; 453 450 long i0=v->color,j0; 454 451 if(i0==-1){ … … 726 723 Triangle &t =triangles[i]; 727 724 Triangle* ta; 728 MeshVertex *v0,*v1,*v2,*v3;725 BamgVertex *v0,*v1,*v2,*v3; 729 726 if (reft[i]<0) continue; 730 727 if ((ta=t.Quadrangle(v0,v1,v2,v3)) !=0 && &t<ta) { … … 775 772 VertexOnGeom &v=VerticesOnGeomVertex[i]; 776 773 ISSMASSERT(v.OnGeomVertex()); 777 bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number(( MeshVertex*)v)+1; //back to Matlab indexing774 bamgmesh->VerticesOnGeometricVertex[i*2+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing 778 775 bamgmesh->VerticesOnGeometricVertex[i*2+1]=Gh.Number((GeometricalVertex*)v)+1; //back to Matlab indexing 779 776 } … … 791 788 ISSMERROR("A vertices supposed to be OnGeometricEdge is actually not"); 792 789 } 793 bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number(( MeshVertex*)v)+1; //back to Matlab indexing790 bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((BamgVertex*)v)+1; //back to Matlab indexing 794 791 bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing 795 792 bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; //absisce … … 1020 1017 for (int j=0;j<2;j++){ 1021 1018 1022 MeshVertex V;1019 BamgVertex V; 1023 1020 VertexOnGeom GV; 1024 1021 Gh.ProjectOnCurve(edges[i],ss[j],V,GV); … … 1066 1063 /*}}}1*/ 1067 1064 /*FUNCTION Mesh::AddVertex{{{1*/ 1068 void Mesh::AddVertex( MeshVertex &s,Triangle* t, Icoor2* det3) {1065 void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3) { 1069 1066 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/ 1070 1067 // ------------------------------------------- … … 1086 1083 Triangle* tt[3]; 1087 1084 //three vertices of t 1088 MeshVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];1085 BamgVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2]; 1089 1086 //three determinants 1090 1087 Icoor2 det3local[3]; … … 1435 1432 for (j=0;j<2;j++){ 1436 1433 //get current vertex 1437 MeshVertex* v=edges[i].v[j];1434 BamgVertex* v=edges[i].v[j]; 1438 1435 //get vertex color (i0) 1439 1436 long i0=v->color; … … 1575 1572 for (i=0;i<nbv;i++){ 1576 1573 if((j=colorV[i])>=0){ 1577 MeshVertex & v = Gh.vertices[j];1574 BamgVertex & v = Gh.vertices[j]; 1578 1575 v = vertices[i]; 1579 1576 v.color =0; … … 2326 2323 delete [] Edgeflags; 2327 2324 2328 //Reset MeshVertex to On2325 //Reset BamgVertex to On 2329 2326 SetVertexFieldOn(); 2330 2327 … … 2557 2554 ISSMERROR("&e==NULL"); 2558 2555 } 2559 MeshVertex * v0 = e(0),*v1 = e(1);2556 BamgVertex * v0 = e(0),*v1 = e(1); 2560 2557 Triangle *t = v0->t; 2561 2558 int sens = Gh.subdomains[i].sens; … … 2645 2642 2646 2643 /*Call NearestVertex*/ 2647 MeshVertex *a = quadtree->NearestVertex(B.x,B.y) ;2644 BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ; 2648 2645 2649 2646 /*Check output (Vertex a)*/ … … 2722 2719 } 2723 2720 /*}}}1*/ 2724 /*FUNCTION Mesh::GeomToTriangles0{{{1*/ 2725 void Mesh::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){ 2721 /*FUNCTION Mesh::Init{{{1*/ 2722 void Mesh::Init(long maxnbv_in) { 2723 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2724 2725 long int verbose=0; 2726 2727 srand(19999999); 2728 NbRef=0; 2729 NbOfTriangleSearchFind =0; 2730 NbOfSwapTriangle =0; 2731 nbv=0; 2732 maxnbv=maxnbv_in; 2733 nbt=0; 2734 NbOfQuad = 0; 2735 nbtx=2*maxnbv_in-2; 2736 NbSubDomains=0; 2737 NbVertexOnBThVertex=0; 2738 NbVertexOnBThEdge=0; 2739 VertexOnBThVertex=NULL; 2740 VertexOnBThEdge=NULL; 2741 2742 NbCrackedVertices=0; 2743 CrackedVertices =NULL; 2744 NbCrackedEdges =0; 2745 CrackedEdges =NULL; 2746 nbe = 0; 2747 2748 if (maxnbv_in) { 2749 vertices=new BamgVertex[maxnbv]; 2750 ISSMASSERT(vertices); 2751 ordre=new (BamgVertex* [maxnbv]); 2752 ISSMASSERT(ordre); 2753 triangles=new Triangle[nbtx]; 2754 ISSMASSERT(triangles); 2755 } 2756 else { 2757 vertices=NULL; 2758 ordre=NULL; 2759 triangles=NULL; 2760 nbtx=0; 2761 } 2762 2763 quadtree=NULL; 2764 edges=NULL; 2765 VerticesOnGeomVertex=NULL; 2766 VerticesOnGeomEdge=NULL; 2767 NbVerticesOnGeomVertex=0; 2768 NbVerticesOnGeomEdge=0; 2769 subdomains=NULL; 2770 NbSubDomains=0; 2771 } 2772 /*}}}1*/ 2773 /*FUNCTION Mesh::Insert{{{1*/ 2774 void Mesh::Insert() { 2775 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2776 2777 /*Insert points in the existing Geometry*/ 2778 2779 //Intermediary 2780 int i; 2781 2782 /*Get options*/ 2783 long int verbose=2; 2784 2785 //Display info 2786 if (verbose>2) printf(" Insert initial %i vertices\n",nbv); 2787 2788 //Compute integer coordinates and determinants for the existing vertices (from Geometry) 2789 SetIntCoor(); 2790 2791 /*Now we want to build a list (ordre) of the vertices in a random 2792 * order. To do so, we use the following method: 2793 * 2794 * From an initial k0 in [0 nbv[ random (vertex number) 2795 * the next k (vertex number) is computed using a big 2796 * prime number (PN>>nbv) following: 2797 * 2798 * k_{i+1} = k_i + PN [nbv] 2799 * 2800 * let's show that: 2801 * 2802 * for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j 2803 * 2804 * Let's assume that there are 2 distinct j1 and j2 such that 2805 * k_j1 = k_j2 2806 * 2807 * This means that 2808 * 2809 * k0+j1*PN = k0+j2*PN [nbv] 2810 * (j1-j2)*PN =0 [nbv] 2811 * since PN is a prime number larger than nbv, and nbv!=1 2812 * j1-j2=0 [nbv] 2813 * BUT 2814 * j1 and j2 are in [0 nbv[ which is impossible. 2815 * 2816 * We hence have built a random list of nbv elements of 2817 * [0 nbv[ all distincts*/ 2818 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ; 2819 const long PrimeNumber= BigPrimeNumber(nbv) ; 2820 int k0=rand()%nbv; 2821 for (int is3=0; is3<nbv; is3++){ 2822 ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv]; 2823 } 2824 2825 /*Modify ordre such that the first 3 vertices form a triangle*/ 2826 2827 //get first vertex i such that [0,1,i] are not aligned 2828 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){ 2829 //if i is higher than nbv, it means that all the determinants are 0, 2830 //all vertices are aligned! 2831 if ( ++i >= nbv) { 2832 ISSMERROR("all the vertices are aligned"); 2833 } 2834 } 2835 // exchange i et 2 in "ordre" so that 2836 // the first 3 vertices are not aligned (real triangle) 2837 Exchange(ordre[2], ordre[i]); 2838 2839 /*Take the first edge formed by the first two vertices and build 2840 * the initial simple mesh from this edge and 2 boundary triangles*/ 2841 2842 BamgVertex * v0=ordre[0], *v1=ordre[1]; 2843 2844 nbt = 2; 2845 triangles[0](0) = NULL; //infinite vertex 2846 triangles[0](1) = v0; 2847 triangles[0](2) = v1; 2848 triangles[1](0) = NULL;//infinite vertex 2849 triangles[1](2) = v0; 2850 triangles[1](1) = v1; 2851 2852 //Build adjacence 2853 const int e0 = OppositeEdge[0]; 2854 const int e1 = NextEdge[e0]; 2855 const int e2 = PreviousEdge[e0]; 2856 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 2857 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 2858 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 2859 2860 triangles[0].det = -1; //boundary triangle: det = -1 2861 triangles[1].det = -1; //boundary triangle: det = -1 2862 2863 triangles[0].SetTriangleContainingTheVertex(); 2864 triangles[1].SetTriangleContainingTheVertex(); 2865 2866 triangles[0].link=&triangles[1]; 2867 triangles[1].link=&triangles[0]; 2868 2869 //build quadtree 2870 if (!quadtree) quadtree = new QuadTree(this,0); 2871 quadtree->Add(*v0); 2872 quadtree->Add(*v1); 2873 2874 /*Now, add the vertices One by One*/ 2875 2876 long NbSwap=0; 2877 if (verbose>3) printf(" Begining of insertion process...\n"); 2878 2879 for (int icount=2; icount<nbv; icount++) { 2880 2881 //Get new vertex 2882 BamgVertex *newvertex=ordre[icount]; 2883 2884 //Find the triangle in which newvertex is located 2885 Icoor2 dete[3]; 2886 Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates) 2887 2888 //Add newvertex to the quadtree 2889 quadtree->Add(*newvertex); 2890 2891 //Add newvertex to the existing mesh 2892 AddVertex(*newvertex,tcvi,dete); 2893 2894 //Make the mesh Delaunay around newvertex by swaping the edges 2895 NbSwap += newvertex->Optim(1,0); 2896 } 2897 2898 //Display info 2899 if (verbose>3) { 2900 printf(" NbSwap of insertion: %i\n",NbSwap); 2901 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 2902 } 2903 2904 #ifdef NBLOOPOPTIM 2905 2906 k0 = rand()%nbv ; 2907 for (int is4=0; is4<nbv; is4++) 2908 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv]; 2909 2910 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){ 2911 long NbSwap = 0; 2912 for (int is1=0; is1<nbv; is1++) 2913 NbSwap += ordre[is1]->Optim(0,0); 2914 if (verbose>3) { 2915 printf(" Optim Loop: %i\n",Nbloop); 2916 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 2917 } 2918 if(!NbSwap) break; 2919 } 2920 ReMakeTriangleContainingTheVertex(); 2921 // because we break the TriangleContainingTheVertex 2922 #endif 2923 } 2924 /*}}}1*/ 2925 /*FUNCTION Mesh::InsertNewPoints{{{1*/ 2926 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) { 2927 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2928 2929 long int verbose=0; 2930 double seuil= 1.414/2 ;// for two close point 2931 long i; 2932 long NbSwap=0; 2933 Icoor2 dete[3]; 2934 2935 //number of new points 2936 const long nbvnew=nbv-nbvold; 2937 2938 //display info if required 2939 if (verbose>5) printf(" Try to Insert %i new points\n",nbvnew); 2940 2941 //return if no new points 2942 if (!nbvnew) return 0; 2943 2944 /*construction of a random order*/ 2945 const long PrimeNumber= BigPrimeNumber(nbv) ; 2946 //remainder of the division of rand() by nbvnew 2947 long k3 = rand()%nbvnew; 2948 //loop over the new points 2949 for (int is3=0; is3<nbvnew; is3++){ 2950 register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew); 2951 register long i=nbvold+is3; 2952 ordre[i]= vertices + j; 2953 ordre[i]->ReferenceNumber=i; 2954 } 2955 2956 // for all the new point 2957 long iv=nbvold; 2958 for (i=nbvold;i<nbv;i++){ 2959 BamgVertex &vi=*ordre[i]; 2960 vi.i=toI2(vi.r); 2961 vi.r=toR2(vi.i); 2962 double hx,hy; 2963 vi.m.Box(hx,hy); 2964 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor); 2965 if (!quadtree->ToClose(vi,seuil,hi,hj)){ 2966 // a good new point 2967 BamgVertex &vj = vertices[iv]; 2968 long j=vj.ReferenceNumber; 2969 if (&vj!=ordre[j]){ 2970 ISSMERROR("&vj!= ordre[j]"); 2971 } 2972 if(i!=j){ 2973 Exchange(vi,vj); 2974 Exchange(ordre[j],ordre[i]); 2975 } 2976 vj.ReferenceNumber=0; 2977 Triangle *tcvj=FindTriangleContaining(vj.i,dete); 2978 if (tcvj && !tcvj->link){ 2979 tcvj->Echo(); 2980 ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link); 2981 } 2982 quadtree->Add(vj); 2983 AddVertex(vj,tcvj,dete); 2984 NbSwap += vj.Optim(1); 2985 iv++; 2986 } 2987 } 2988 if (verbose>3) { 2989 printf(" number of new points: %i\n",iv); 2990 printf(" number of to close (?) points: %i\n",nbv-iv); 2991 printf(" number of swap after: %i\n",NbSwap); 2992 } 2993 nbv = iv; 2994 2995 for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1); 2996 if (verbose>3) printf(" NbSwap=%i\n",NbSwap); 2997 2998 NbTSwap += NbSwap ; 2999 return nbv-nbvold; 3000 } 3001 /*}}}1*/ 3002 /*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/ 3003 Edge** Mesh::MakeGeometricalEdgeToEdge() { 3004 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/ 3005 3006 if (!Gh.nbe){ 3007 ISSMERROR("!Gh.nbe"); 3008 } 3009 Edge **e= new (Edge* [Gh.nbe]); 3010 3011 long i; 3012 for ( i=0;i<Gh.nbe ; i++) 3013 e[i]=NULL; 3014 for ( i=0;i<nbe ; i++) 3015 { 3016 Edge * ei = edges+i; 3017 GeometricalEdge *onGeometry = ei->onGeometry; 3018 e[Gh.Number(onGeometry)] = ei; 3019 } 3020 for ( i=0;i<nbe ; i++) 3021 for (int ii=0;ii<2;ii++) { 3022 Edge * ei = edges+i; 3023 GeometricalEdge *onGeometry = ei->onGeometry; 3024 int j= ii; 3025 while (!(*onGeometry)[j].Required()) { 3026 Adj(onGeometry,j); // next geom edge 3027 j=1-j; 3028 if (e[Gh.Number(onGeometry)]) break; // optimisation 3029 e[Gh.Number(onGeometry)] = ei; 3030 } 3031 } 3032 3033 int kk=0; 3034 for ( i=0;i<Gh.nbe ; i++){ 3035 if (!e[i]){ 3036 kk++; 3037 if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i); 3038 } 3039 } 3040 if(kk) ISSMERROR("See above"); 3041 3042 return e; 3043 } 3044 /*}}}1*/ 3045 /*FUNCTION Mesh::MakeQuadrangles{{{1*/ 3046 void Mesh::MakeQuadrangles(double costheta){ 3047 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/ 3048 3049 long int verbose=0; 3050 3051 if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta); 3052 3053 if (costheta >1) { 3054 if (verbose>5) printf(" do nothing: costheta > 1\n"); 3055 } 3056 3057 3058 long nbqq = (nbt*3)/2; 3059 DoubleAndInt *qq = new DoubleAndInt[nbqq]; 3060 3061 long i,ij; 3062 int j; 3063 long k=0; 3064 for (i=0;i<nbt;i++) 3065 for (j=0;j<3;j++) 3066 if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta) 3067 qq[k++].i3j=i*3+j; 3068 // sort qq 3069 HeapSort(qq,k); 3070 3071 long kk=0; 3072 for (ij=0;ij<k;ij++) { 3073 i=qq[ij].i3j/3; 3074 j=(int) (qq[ij].i3j%3); 3075 // optisamition no float computation 3076 if (triangles[i].QualityQuad(j,0) >=costheta) 3077 triangles[i].SetHidden(j),kk++; 3078 } 3079 NbOfQuad = kk; 3080 if (verbose>2){ 3081 printf(" number of quadrilaterals = %i\n",NbOfQuad); 3082 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); 3083 printf(" number of outside triangles = %i\n",NbOutT); 3084 } 3085 delete [] qq; 3086 } 3087 /*}}}1*/ 3088 /*FUNCTION Mesh::MakeQuadTree{{{1*/ 3089 void Mesh::MakeQuadTree() { 3090 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/ 3091 3092 long int verbose=0; 3093 if ( !quadtree ) quadtree = new QuadTree(this); 3094 3095 } 3096 /*}}}1*/ 3097 /*FUNCTION Mesh::MaxSubDivision{{{1*/ 3098 void Mesh::MaxSubDivision(double maxsubdiv) { 3099 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/ 3100 3101 long int verbose=0; 3102 3103 const double maxsubdiv2 = maxsubdiv*maxsubdiv; 3104 if(verbose>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv); 3105 // for all the edges 3106 // if the len of the edge is to long 3107 long it,nbchange=0; 3108 double lmax=0; 3109 for (it=0;it<nbt;it++){ 3110 Triangle &t=triangles[it]; 3111 for (int j=0;j<3;j++){ 3112 Triangle &tt = *t.TriangleAdj(j); 3113 if ( ! &tt || it < Number(tt) && ( tt.link || t.link)){ 3114 BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 3115 BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 3116 R2 AB= (R2) v1-(R2) v0; 3117 Metric M = v0; 3118 double l = M(AB,AB); 3119 lmax = Max(lmax,l); 3120 if(l> maxsubdiv2){ 3121 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3122 double lc = M(AC,AC); 3123 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 3124 D2xD2 Rt1(Rt.inv()); 3125 D2xD2 D(maxsubdiv2,0,0,lc); 3126 D2xD2 MM = Rt1*D*Rt1.t(); 3127 v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y); 3128 nbchange++; 3129 } 3130 M = v1; 3131 l = M(AB,AB); 3132 lmax = Max(lmax,l); 3133 if(l> maxsubdiv2){ 3134 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3135 double lc = M(AC,AC); 3136 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 3137 D2xD2 Rt1(Rt.inv()); 3138 D2xD2 D(maxsubdiv2,0,0,lc); 3139 D2xD2 MM = Rt1*D*Rt1.t(); 3140 v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y); 3141 nbchange++; 3142 } 3143 } 3144 } 3145 } 3146 if(verbose>3){ 3147 printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5)); 3148 } 3149 } 3150 /*}}}1*/ 3151 /*FUNCTION Mesh::MetricAt{{{1*/ 3152 Metric Mesh::MetricAt(const R2 & A) const { 3153 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/ 3154 3155 I2 a = toI2(A); 3156 Icoor2 deta[3]; 3157 Triangle * t =FindTriangleContaining(a,deta); 3158 if (t->det <0) { // outside 3159 double ba,bb; 3160 TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ; 3161 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));} 3162 else { // inside 3163 double aa[3]; 3164 double s = deta[0]+deta[1]+deta[2]; 3165 aa[0]=deta[0]/s; 3166 aa[1]=deta[1]/s; 3167 aa[2]=deta[2]/s; 3168 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]); 3169 } 3170 } 3171 /*}}}1*/ 3172 /*FUNCTION Mesh::NearestVertex{{{1*/ 3173 BamgVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) { 3174 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/ 3175 return quadtree->NearestVertex(i,j); 3176 } 3177 /*}}}1*/ 3178 /*FUNCTION Mesh::NewPoints{{{1*/ 3179 void Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){ 3180 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/ 3181 3182 int i,j,k; 3183 long NbTSwap=0; 3184 long nbtold=nbt; 3185 long nbvold=nbv; 3186 long Headt=0; 3187 long next_t; 3188 long* first_np_or_next_t=new long[nbtx]; 3189 Triangle* t=NULL; 3190 3191 /*Recover options*/ 3192 int verbose=bamgopts->verbose; 3193 3194 /*First, insert old points if requested*/ 3195 if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){ 3196 if (verbose>5) printf(" Inserting initial mesh points\n"); 3197 for (i=0;i<Bh.nbv;i++){ 3198 BamgVertex &bv=Bh[i]; 3199 if (!bv.onGeometry){ 3200 vertices[nbv].r = bv.r; 3201 vertices[nbv++].m = bv.m; 3202 } 3203 } 3204 Bh.ReMakeTriangleContainingTheVertex(); 3205 InsertNewPoints(nbvold,NbTSwap) ; 3206 } 3207 else Bh.ReMakeTriangleContainingTheVertex(); 3208 3209 // generation of the list of next Triangle 3210 for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1); 3211 // the next traingle of i is -first_np_or_next_t[i] 3212 3213 // Big loop (most time consuming) 3214 int iter=0; 3215 if (verbose>5) printf(" Big loop\n"); 3216 do { 3217 /*Update variables*/ 3218 iter++; 3219 nbtold=nbt; 3220 nbvold=nbv; 3221 3222 /*We test all triangles*/ 3223 i=Headt; 3224 next_t=-first_np_or_next_t[i]; 3225 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){ 3226 3227 //check i 3228 if (i<0 || i>=nbt ){ 3229 ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1); 3230 } 3231 //change first_np_or_next_t[i] 3232 first_np_or_next_t[i] = iter; 3233 3234 //Loop over the edges of t 3235 for(j=0;j<3;j++){ 3236 TriangleAdjacent tj(t,j); 3237 BamgVertex &vA = *tj.EdgeVertex(0); 3238 BamgVertex &vB = *tj.EdgeVertex(1); 3239 3240 //if t is a boundary triangle, or tj locked, continue 3241 if (!t->link) continue; 3242 if (t->det <0) continue; 3243 if (t->Locked(j)) continue; 3244 3245 TriangleAdjacent tadjj = t->Adj(j); 3246 Triangle* ta=tadjj; 3247 3248 //if the adjacent triangle is a boundary triangle, continur 3249 if (ta->det<0) continue; 3250 3251 R2 A=vA; 3252 R2 B=vB; 3253 k=Number(ta); 3254 3255 //if this edge has already been done, go to next edge of triangle 3256 if(first_np_or_next_t[k]==iter) continue; 3257 3258 lIntTria.SplitEdge(Bh,A,B); 3259 lIntTria.NewPoints(vertices,nbv,maxnbv); 3260 } // end loop for each edge 3261 }// for triangle 3262 3263 if (!InsertNewPoints(nbvold,NbTSwap)) break; 3264 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3265 Headt = nbt; // empty list 3266 3267 // for all the triangle containing the vertex i 3268 for (i=nbvold;i<nbv;i++){ 3269 BamgVertex* s = vertices + i; 3270 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]); 3271 Triangle* tbegin= (Triangle*) ta; 3272 long kt; 3273 do { 3274 kt = Number((Triangle*) ta); 3275 if (first_np_or_next_t[kt]>0){ 3276 first_np_or_next_t[kt]=-Headt; 3277 Headt=kt; 3278 } 3279 if (ta.EdgeVertex(0)!=s){ 3280 ISSMERROR("ta.EdgeVertex(0)!=s"); 3281 } 3282 ta = Next(Adj(ta)); 3283 } while ( (tbegin != (Triangle*) ta)); 3284 } 3285 3286 } while (nbv!=nbvold); 3287 3288 delete [] first_np_or_next_t; 3289 3290 long NbSwapf =0,NbSwp; 3291 3292 NbSwp = NbSwapf; 3293 for (i=0;i<nbv;i++) 3294 NbSwapf += vertices[i].Optim(0); 3295 NbTSwap += NbSwapf ; 3296 } 3297 /*}}}1*/ 3298 /*FUNCTION Mesh::ProjectOnCurve{{{1*/ 3299 GeometricalEdge* Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex & vA, BamgVertex & vB, 3300 double theta,BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR) { 3301 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/ 3302 3303 void *pA=0,*pB=0; 3304 double tA=0,tB=0; 3305 R2 A=vA,B=vB; 3306 BamgVertex * pvA=&vA, * pvB=&vB; 3307 if (vA.vint == IsVertexOnVertex){ 3308 pA=vA.onBackgroundVertex; 3309 } 3310 else if (vA.vint == IsVertexOnEdge){ 3311 pA=vA.onBackgroundEdge->be; 3312 tA=vA.onBackgroundEdge->abcisse; 3313 } 3314 else { 3315 ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA)); 3316 } 3317 3318 if (vB.vint == IsVertexOnVertex){ 3319 pB=vB.onBackgroundVertex; 3320 } 3321 else if(vB.vint == IsVertexOnEdge){ 3322 pB=vB.onBackgroundEdge->be; 3323 tB=vB.onBackgroundEdge->abcisse; 3324 } 3325 else { 3326 ISSMERROR("ProjectOnCurve On BamgVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB)); 3327 } 3328 Edge * e = &BhAB; 3329 if (!pA || !pB || !e){ 3330 ISSMERROR("!pA || !pB || !e"); 3331 } 3332 // be carefull the back ground edge e is on same geom edge 3333 // of the initiale edge def by the 2 vertex A B; 3334 //check Is a background Mesh; 3335 if (e<BTh.edges || e>=BTh.edges+BTh.nbe){ 3336 ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe"); 3337 } 3338 // walk on BTh edge 3339 //not finish ProjectOnCurve with BackGround Mesh); 3340 // 1 first find a back ground edge contening the vertex A 3341 // 2 walk n back gound boundary to find the final vertex B 3342 3343 if( vA.vint == IsVertexOnEdge) 3344 { // find the start edge 3345 e = vA.onBackgroundEdge->be; 3346 3347 } 3348 else if (vB.vint == IsVertexOnEdge) 3349 { 3350 theta = 1-theta; 3351 Exchange(tA,tB); 3352 Exchange(pA,pB); 3353 Exchange(pvA,pvB); 3354 Exchange(A,B); 3355 e = vB.onBackgroundEdge->be; 3356 3357 } 3358 else{ // do the search by walking 3359 ISSMERROR("case not supported yet"); 3360 } 3361 3362 // find the direction of walking with sens of edge and pA,PB; 3363 R2 AB=B-A; 3364 3365 double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); 3366 int kkk=0; 3367 int sens = (cosE01AB>0) ? 1 : 0; 3368 3369 // double l=0; // length of the edge AB 3370 double abscisse = -1; 3371 3372 for (int cas=0;cas<2;cas++){ 3373 // 2 times algo: 3374 // 1 for computing the length l 3375 // 2 for find the vertex 3376 int iii; 3377 BamgVertex *v0=pvA,*v1; 3378 Edge *neee,*eee; 3379 double lg =0; // length of the curve 3380 double te0; 3381 // we suppose take the curve's abcisse 3382 for ( eee=e,iii=sens,te0=tA; 3383 eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; 3384 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { 3385 3386 kkk=kkk+1; 3387 if (kkk>=100){ 3388 ISSMERROR("kkk>=100"); 3389 } 3390 if (!eee){ 3391 ISSMERROR("!eee"); 3392 } 3393 double lg0 = lg; 3394 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); 3395 lg += dp; 3396 if (cas && abscisse <= lg) { // ok we find the geom edge 3397 double sss = (abscisse-lg0)/dp; 3398 double thetab = te0*(1-sss)+ sss*iii; 3399 if (thetab<0 || thetab>1){ 3400 ISSMERROR("thetab<0 || thetab>1"); 3401 } 3402 BR = VertexOnEdge(&R,eee,thetab); 3403 return Gh.ProjectOnCurve(*eee,thetab,R,GR); 3404 } 3405 } 3406 // we find the end 3407 if (v1 != pvB){ 3408 if (( void*) v1 == pB) 3409 tB = iii; 3410 3411 double lg0 = lg; 3412 if (!eee){ 3413 ISSMERROR("!eee"); 3414 } 3415 v1 = pvB; 3416 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); 3417 lg += dp; 3418 abscisse = lg*theta; 3419 if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end 3420 { // ok we find the geom edge 3421 double sss = (abscisse-lg0)/dp; 3422 double thetab = te0*(1-sss)+ sss*tB; 3423 if (thetab<0 || thetab>1){ 3424 ISSMERROR("thetab<0 || thetab>1"); 3425 } 3426 BR = VertexOnEdge(&R,eee,thetab); 3427 return Gh.ProjectOnCurve(*eee,thetab,R,GR); 3428 } 3429 } 3430 abscisse = lg*theta; 3431 3432 } 3433 ISSMERROR("Big bug..."); 3434 return 0; // just for the compiler 3435 } 3436 /*}}}1*/ 3437 /*FUNCTION Mesh::ReconstructExistingMesh{{{1*/ 3438 void Mesh::ReconstructExistingMesh(){ 3439 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/ 3440 3441 /*This routine reconstruct an existing mesh to make it CONVEX: 3442 * -all the holes are filled 3443 * -concave boundaries are filled 3444 * A convex mesh is required for a lot of operations. This is why every mesh 3445 * goes through this process. 3446 * This routine also generates mesh properties such as adjencies,... 3447 */ 3448 3449 /*Intermediary*/ 3450 int verbose=0; 3451 3452 // generation of the integer coordinate 3453 3454 // find extrema coordinates of vertices pmin,pmax 3455 long i; 3456 if(verbose>2) printf(" Reconstruct mesh of %i vertices\n",nbv); 3457 3458 //initialize ordre 3459 ISSMASSERT(ordre); 3460 for (i=0;i<nbv;i++) ordre[i]=0; 3461 3462 //Initialize NbSubDomains 3463 NbSubDomains =0; 3464 3465 /* generation of triangles adjacency*/ 3466 3467 //First add existing edges 3468 long kk =0; 3469 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 3470 for (i=0;i<nbe;i++){ 3471 kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 3472 } 3473 if (kk != nbe){ 3474 ISSMERROR("There are %i double edges in the mesh",kk-nbe); 3475 } 3476 3477 //Add edges of all triangles in existing mesh 3478 long* st = new long[nbt*3]; 3479 for (i=0;i<nbt*3;i++) st[i]=-1; 3480 for (i=0;i<nbt;i++){ 3481 for (int j=0;j<3;j++){ 3482 3483 //Add current triangle edge to edge4 3484 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 3485 3486 long invisible=triangles[i].Hidden(j); 3487 3488 //If the edge has not been added to st, add it 3489 if(st[k]==-1) st[k]=3*i+j; 3490 3491 //If the edge already exists, add adjacency 3492 else if(st[k]>=0) { 3493 ISSMASSERT(!triangles[i].TriangleAdj(j)); 3494 ISSMASSERT(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3))); 3495 3496 triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3)); 3497 if (invisible) triangles[i].SetHidden(j); 3498 if (k<nbe) triangles[i].SetLocked(j); 3499 3500 //Make st[k] negative so that it will throw an error message if it is found again 3501 st[k]=-2-st[k]; 3502 } 3503 3504 //An edge belongs to 2 triangles 3505 else { 3506 ISSMERROR("The edge (%i , %i) belongs to more than 2 triangles",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 3507 } 3508 } 3509 } 3510 3511 //Display info if required 3512 if(verbose>5) { 3513 printf(" info of Mesh:\n"); 3514 printf(" - number of vertices = %i \n",nbv); 3515 printf(" - number of triangles = %i \n",nbt); 3516 printf(" - number of given edges = %i \n",nbe); 3517 printf(" - number of all edges = %i \n" ,edge4->nb()); 3518 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 3519 } 3520 3521 //check the consistency of edge[].adj and the geometrical required vertex 3522 long k=0; 3523 for (i=0;i<edge4->nb();i++){ 3524 if (st[i]>=0){ // edge alone 3525 if (i<nbe){ 3526 long i0=edge4->i(i); 3527 ordre[i0] = vertices+i0; 3528 long i1=edge4->j(i); 3529 ordre[i1] = vertices+i1; 3530 } 3531 else { 3532 k=k+1; 3533 if (k<10) { 3534 //print only 10 edges 3535 printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i)); 3536 } 3537 else if (k==10){ 3538 printf("Other lost boundary edges not shown...\n"); 3539 } 3540 } 3541 } 3542 } 3543 if(k) { 3544 ISSMERROR("%i boundary edges (from the geometry) are not defined as mesh edges",k); 3545 } 3546 3547 /* mesh generation with boundary points*/ 3548 long nbvb=0; 3549 for (i=0;i<nbv;i++){ 3550 vertices[i].t=0; 3551 vertices[i].vint=0; 3552 if (ordre[i]) ordre[nbvb++]=ordre[i]; 3553 } 3554 3555 Triangle* savetriangles=triangles; 3556 long savenbt=nbt; 3557 long savenbtx=nbtx; 3558 SubDomain* savesubdomains=subdomains; 3559 subdomains=0; 3560 3561 long Nbtriafillhole=2*nbvb; 3562 Triangle* triafillhole=new Triangle[Nbtriafillhole]; 3563 triangles = triafillhole; 3564 3565 nbt=2; 3566 nbtx= Nbtriafillhole; 3567 3568 //Find a vertex that is not aligned with vertices 0 and 1 3569 for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;) 3570 if (++i>=nbvb) { 3571 ISSMERROR("ReconstructExistingMesh: All the vertices are aligned"); 3572 } 3573 //Move this vertex (i) to the 2d position in ordre 3574 Exchange(ordre[2], ordre[i]); 3575 3576 /*Reconstruct mesh beginning with 2 triangles*/ 3577 BamgVertex * v0=ordre[0], *v1=ordre[1]; 3578 3579 triangles[0](0) = NULL; // Infinite vertex 3580 triangles[0](1) = v0; 3581 triangles[0](2) = v1; 3582 3583 triangles[1](0) = NULL;// Infinite vertex 3584 triangles[1](2) = v0; 3585 triangles[1](1) = v1; 3586 const int e0 = OppositeEdge[0]; 3587 const int e1 = NextEdge[e0]; 3588 const int e2 = PreviousEdge[e0]; 3589 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 3590 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 3591 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 3592 3593 triangles[0].det = -1; // boundary triangles 3594 triangles[1].det = -1; // boundary triangles 3595 3596 triangles[0].SetTriangleContainingTheVertex(); 3597 triangles[1].SetTriangleContainingTheVertex(); 3598 3599 triangles[0].link=&triangles[1]; 3600 triangles[1].link=&triangles[0]; 3601 3602 if (!quadtree) delete quadtree; //ReInitialise; 3603 quadtree = new QuadTree(this,0); 3604 quadtree->Add(*v0); 3605 quadtree->Add(*v1); 3606 3607 // vertices are added one by one 3608 long NbSwap=0; 3609 for (int icount=2; icount<nbvb; icount++) { 3610 BamgVertex *vi = ordre[icount]; 3611 Icoor2 dete[3]; 3612 Triangle *tcvi = FindTriangleContaining(vi->i,dete); 3613 quadtree->Add(*vi); 3614 AddVertex(*vi,tcvi,dete); 3615 NbSwap += vi->Optim(1,1); 3616 } 3617 3618 //enforce the boundary 3619 TriangleAdjacent ta(0,0); 3620 long nbloss = 0,knbe=0; 3621 for ( i = 0; i < nbe; i++){ 3622 if (st[i] >=0){ //edge alone => on border 3623 BamgVertex &a=edges[i][0], &b=edges[i][1]; 3624 if (a.t && b.t){ 3625 knbe++; 3626 if (ForceEdge(a,b,ta)<0) nbloss++; 3627 } 3628 } 3629 } 3630 if(nbloss) { 3631 ISSMERROR("we lost %i existing edges other %i",nbloss,knbe); 3632 } 3633 3634 FindSubDomain(1); 3635 // remove all the hole 3636 // remove all the good sub domain 3637 long krm =0; 3638 for (i=0;i<nbt;i++){ 3639 if (triangles[i].link){ // remove triangles 3640 krm++; 3641 for (int j=0;j<3;j++){ 3642 TriangleAdjacent ta = triangles[i].Adj(j); 3643 Triangle &tta = *(Triangle*)ta; 3644 //if edge between remove and not remove 3645 if(! tta.link){ 3646 // change the link of ta; 3647 int ja = ta; 3648 BamgVertex *v0= ta.EdgeVertex(0); 3649 BamgVertex *v1= ta.EdgeVertex(1); 3650 long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); 3651 3652 ISSMASSERT(st[k]>=0); 3653 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); 3654 ta.SetLock(); 3655 st[k]=-2-st[k]; 3656 } 3657 } 3658 } 3659 } 3660 long NbTfillHoll =0; 3661 for (i=0;i<nbt;i++){ 3662 if (triangles[i].link) { 3663 triangles[i]=Triangle((BamgVertex *) NULL,(BamgVertex *) NULL,(BamgVertex *) NULL); 3664 triangles[i].color=-1; 3665 } 3666 else{ 3667 triangles[i].color= savenbt+ NbTfillHoll++; 3668 } 3669 } 3670 ISSMASSERT(savenbt+NbTfillHoll<=savenbtx); 3671 3672 // copy of the outside triangles in saveMesh 3673 for (i=0;i<nbt;i++){ 3674 if(triangles[i].color>=0) { 3675 savetriangles[savenbt]=triangles[i]; 3676 savetriangles[savenbt].link=0; 3677 savenbt++; 3678 } 3679 } 3680 // gestion of the adj 3681 k =0; 3682 Triangle * tmax = triangles + nbt; 3683 for (i=0;i<savenbt;i++) { 3684 Triangle & ti = savetriangles[i]; 3685 for (int j=0;j<3;j++){ 3686 Triangle * ta = ti.TriangleAdj(j); 3687 int aa = ti.NuEdgeTriangleAdj(j); 3688 int lck = ti.Locked(j); 3689 if (!ta) k++; // bug 3690 else if ( ta >= triangles && ta < tmax){ 3691 ta= savetriangles + ta->color; 3692 ti.SetAdj2(j,ta,aa); 3693 if(lck) ti.SetLocked(j); 3694 } 3695 } 3696 } 3697 3698 // restore triangles; 3699 nbt=savenbt; 3700 nbtx=savenbtx; 3701 delete [] triangles; 3702 delete [] subdomains; 3703 triangles = savetriangles; 3704 subdomains = savesubdomains; 3705 if (k) { 3706 ISSMERROR("number of triangles edges alone = %i",k); 3707 } 3708 FindSubDomain(); 3709 3710 delete edge4; 3711 delete [] st; 3712 for (i=0;i<nbv;i++) quadtree->Add(vertices[i]); 3713 3714 SetVertexFieldOn(); 3715 3716 /*Check requirements consistency*/ 3717 for (i=0;i<nbe;i++){ 3718 /*If the current mesh edge is on Geometry*/ 3719 if(edges[i].onGeometry){ 3720 for(int j=0;j<2;j++){ 3721 /*Go through the edges adjacent to current edge (if on the same curve)*/ 3722 if (!edges[i].adj[j]){ 3723 /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/ 3724 /*Check that the 2 vertices are on geometry AND required*/ 3725 if(!edges[i][j].onGeometry->IsRequiredVertex()){ 3726 printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1); 3727 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1); 3728 if (edges[i][j].onGeometry->OnGeomVertex()) 3729 printf("the vertex number %i of this edge is a geometric BamgVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->gv)+1); 3730 else if (edges[i][j].onGeometry->OnGeomEdge()) 3731 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->ge)+1); 3732 else 3733 printf("Its pointer is %p\n",edges[i][j].onGeometry); 3734 3735 printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n"); 3736 ISSMERROR("See above (might be cryptic...)"); 3737 } 3738 } 3739 } 3740 } 3741 } 3742 } 3743 /*}}}1*/ 3744 /*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/ 3745 void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){ 3746 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/ 3747 3748 long int verbose=0; 3749 long *renu= new long[nbt]; 3750 register Triangle *t0,*t,*te=triangles+nbt; 3751 register long k=0,it,i,j; 3752 3753 for ( it=0;it<nbt;it++) 3754 renu[it]=-1; // outside triangle 3755 for ( i=0;i<NbSubDomains;i++) 3756 { 3757 t=t0=subdomains[i].head; 3758 if (!t0){ // not empty sub domain 3759 ISSMERROR("!t0"); 3760 } 3761 do { 3762 long kt = Number(t); 3763 if (kt<0 || kt >= nbt ){ 3764 ISSMERROR("kt<0 || kt >= nbt"); 3765 } 3766 if (renu[kt]!=-1){ 3767 ISSMERROR("renu[kt]!=-1"); 3768 } 3769 renu[kt]=k++; 3770 } 3771 while (t0 != (t=t->link)); 3772 } 3773 // take is same numbering if possible 3774 if(justcompress) 3775 for ( k=0,it=0;it<nbt;it++) 3776 if(renu[it] >=0 ) 3777 renu[it]=k++; 3778 3779 // put the outside triangles at the end 3780 for ( it=0;it<nbt;it++){ 3781 if (renu[it]==-1) renu[it]=k++; 3782 } 3783 if (k != nbt){ 3784 ISSMERROR("k != nbt"); 3785 } 3786 // do the change on all the pointeur 3787 for ( it=0;it<nbt;it++) 3788 triangles[it].ReNumbering(triangles,te,renu); 3789 3790 for ( i=0;i<NbSubDomains;i++) 3791 subdomains[i].head=triangles+renu[Number(subdomains[i].head)]; 3792 3793 // move the Triangles without a copy of the array 3794 // be carefull not trivial code 3795 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu 3796 if (renu[it] >= 0) // a new sub cycle 3797 { 3798 i=it; 3799 Triangle ti=triangles[i],tj; 3800 while ( (j=renu[i]) >= 0) 3801 { // i is old, and j is new 3802 renu[i] = -1; // mark 3803 tj = triangles[j]; // save new 3804 triangles[j]= ti; // new <- old 3805 i=j; // next 3806 ti = tj; 3807 } 3808 } 3809 delete [] renu; 3810 nt = nbt - NbOutT; 3811 3812 } 3813 /*}}}1*/ 3814 /*FUNCTION Mesh::ReNumberingVertex{{{1*/ 3815 void Mesh::ReNumberingVertex(long * renu) { 3816 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/ 3817 3818 // warning be carfull because pointer 3819 // from on mesh to over mesh 3820 // -- so do ReNumbering at the beginning 3821 BamgVertex * ve = vertices+nbv; 3822 long it,ie,i; 3823 3824 printf("renumbering triangles\n"); 3825 for ( it=0;it<nbt;it++) 3826 triangles[it].ReNumbering(vertices,ve,renu); 3827 3828 printf("renumbering edges\n"); 3829 for ( ie=0;ie<nbe;ie++) 3830 edges[ie].ReNumbering(vertices,ve,renu); 3831 3832 printf("renumbering vertices on geom\n"); 3833 for (i=0;i< NbVerticesOnGeomVertex;i++) 3834 { 3835 BamgVertex *v = VerticesOnGeomVertex[i].mv; 3836 if (v>=vertices && v < ve) 3837 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)]; 3838 } 3839 3840 printf("renumbering vertices on edge\n"); 3841 for (i=0;i< NbVerticesOnGeomEdge;i++) 3842 { 3843 BamgVertex *v =VerticesOnGeomEdge[i].mv; 3844 if (v>=vertices && v < ve) 3845 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)]; 3846 } 3847 3848 printf("renumbering vertices on Bth vertex\n"); 3849 for (i=0;i< NbVertexOnBThVertex;i++) 3850 { 3851 BamgVertex *v=VertexOnBThVertex[i].v; 3852 if (v>=vertices && v < ve) 3853 VertexOnBThVertex[i].v=vertices+renu[Number(v)]; 3854 } 3855 3856 for (i=0;i< NbVertexOnBThEdge;i++) 3857 { 3858 BamgVertex *v=VertexOnBThEdge[i].v; 3859 if (v>=vertices && v < ve) 3860 VertexOnBThEdge[i].v=vertices+renu[Number(v)]; 3861 } 3862 3863 // move the Vertices without a copy of the array 3864 // be carefull not trivial code 3865 long j; 3866 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu 3867 if (renu[it] >= 0) // a new sub cycle 3868 { 3869 i=it; 3870 BamgVertex ti=vertices[i],tj; 3871 while ( (j=renu[i]) >= 0){ 3872 // i is old, and j is new 3873 renu[i] = -1-renu[i]; // mark 3874 tj = vertices[j]; // save new 3875 vertices[j]= ti; // new <- old 3876 i=j; // next 3877 ti = tj; 3878 } 3879 } 3880 if (quadtree){ 3881 delete quadtree; 3882 quadtree = new QuadTree(this); 3883 } 3884 3885 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1; 3886 } 3887 /*}}}1*/ 3888 /*FUNCTION Mesh::SetIntCoor{{{1*/ 3889 void Mesh::SetIntCoor(const char * strfrom) { 3890 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/ 3891 3892 /*Set integer coordinate for existing vertices*/ 3893 3894 //Get extrema coordinates of the existing vertices 3895 pmin = vertices[0].r; 3896 pmax = vertices[0].r; 3897 long i; 3898 for (i=0;i<nbv;i++) { 3899 pmin.x = Min(pmin.x,vertices[i].r.x); 3900 pmin.y = Min(pmin.y,vertices[i].r.y); 3901 pmax.x = Max(pmax.x,vertices[i].r.x); 3902 pmax.y = Max(pmax.y,vertices[i].r.y); 3903 } 3904 R2 DD = (pmax-pmin)*0.05; 3905 pmin = pmin-DD; 3906 pmax = pmax+DD; 3907 3908 //Compute coefIcoor 3909 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 3910 if (coefIcoor<=0){ 3911 ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely"); 3912 } 3913 3914 // generation of integer coord 3915 for (i=0;i<nbv;i++) { 3916 vertices[i].i = toI2(vertices[i].r); 3917 } 3918 3919 // computation of the det 3920 int number_of_errors=0; 3921 for (i=0;i<nbt;i++) { 3922 BamgVertex & v0 = triangles[i][0]; 3923 BamgVertex & v1 = triangles[i][1]; 3924 BamgVertex & v2 = triangles[i][2]; 3925 3926 //If this is not a boundary triangle 3927 if ( &v0 && &v1 && &v2 ){ 3928 3929 /*Compute determinant*/ 3930 triangles[i].det= det(v0,v1,v2); 3931 3932 /*Check that determinant is positive*/ 3933 if (triangles[i].det <=0){ 3934 3935 /*increase number_of_errors and print error only for the first 20 triangles*/ 3936 number_of_errors++; 3937 if (number_of_errors<20){ 3938 printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det); 3939 } 3940 } 3941 } 3942 3943 //else, set as -1 3944 else triangles[i].det=-1; 3945 } 3946 3947 if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above"); 3948 } 3949 /*}}}1*/ 3950 /*FUNCTION Mesh::ShowRegulaty{{{1*/ 3951 void Mesh::ShowRegulaty() const { 3952 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/ 3953 3954 const double sqrt32=sqrt(3.)*0.5; 3955 const double aireKh=sqrt32*0.5; 3956 D2 Beq(1,0),Heq(0.5,sqrt32); 3957 D2xD2 Br(D2xD2(Beq,Heq).t()); 3958 D2xD2 B1r(Br.inv()); 3959 double gammamn=1e100,hmin=1e100; 3960 double gammamx=0,hmax=0; 3961 double beta=1e100; 3962 double beta0=0; 3963 double alpha2=0; 3964 double area=0,Marea=0; 3965 // double cf= double(coefIcoor); 3966 // double cf2= 6.*cf*cf; 3967 int nt=0; 3968 for (int it=0;it<nbt;it++) 3969 if ( triangles[it].link) 3970 { 3971 nt++; 3972 Triangle &K=triangles[it]; 3973 double area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.; 3974 area+= area3; 3975 D2xD2 B_Kt(K[0],K[1],K[2]); 3976 D2xD2 B_K(B_Kt.t()); 3977 D2xD2 B1K = Br*B_K.inv(); 3978 D2xD2 BK = B_K*B1r; 3979 D2xD2 B1B1 = B1K.t()*B1K; 3980 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y); 3981 MatVVP2x2 VMK(MK); 3982 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1)); 3983 double betaK=0; 3984 3985 for (int j=0;j<3;j++) 3986 { 3987 double he= Norme2(R2(K[j],K[(j+1)%3])); 3988 hmin=Min(hmin,he); 3989 hmax=Max(hmax,he); 3990 BamgVertex & v=K[j]; 3991 D2xD2 M((Metric)v); 3992 betaK += sqrt(M.det()); 3993 3994 D2xD2 BMB = BK.t()*M*BK; 3995 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y); 3996 MatVVP2x2 VM1(M1); 3997 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2); 3998 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2); 3999 } 4000 betaK *= area3;// 1/2 (somme sqrt(det))* area(K) 4001 Marea+= betaK; 4002 beta=min(beta,betaK); 4003 beta0=max(beta0,betaK); 4004 } 4005 area*=3; 4006 gammamn=sqrt(gammamn); 4007 gammamx=sqrt(gammamx); 4008 printf(" Adaptmesh info:\n"); 4009 printf(" number of triangles = %i\n",nt); 4010 printf(" hmin = %g, hmax=%g\n",hmin,hmax); 4011 printf(" area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt)); 4012 printf(" infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx); 4013 printf(" anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5)); 4014 } 4015 /*}}}1*/ 4016 /*FUNCTION Mesh::ShowHistogram{{{1*/ 4017 void Mesh::ShowHistogram() const { 4018 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/ 4019 4020 const long kmax=10; 4021 const double llmin = 0.5,llmax=2; 4022 const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin); 4023 long histo[kmax+1]; 4024 long i,it,k, nbedges =0; 4025 for (i=0;i<=kmax;i++) histo[i]=0; 4026 for (it=0;it<nbt;it++) 4027 if ( triangles[it].link) 4028 { 4029 4030 for (int j=0;j<3;j++) 4031 { 4032 Triangle *ta = triangles[it].TriangleAdj(j); 4033 if ( !ta || !ta->link || Number(ta) >= it) 4034 { 4035 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]]; 4036 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]]; 4037 if ( !& vP || !&vQ) continue; 4038 R2 PQ = vQ.r - vP.r; 4039 double l = log(LengthInterpole(vP,vQ,PQ)); 4040 nbedges++; 4041 k = (int) ((l - lmin)*delta); 4042 k = Min(Max(k,0L),kmax); 4043 histo[k]++; 4044 } 4045 } 4046 } 4047 printf(" --- Histogram of the unit mesh, nb of edges = %i\n",nbedges); 4048 printf(" length of edge in | %% of edge | Nb of edges \n"); 4049 printf(" --------------------+-------------+-------------\n"); 4050 for (i=0;i<=kmax;i++){ 4051 if (i==0) printf(" %10i",0); 4052 else printf(" %10g",exp(lmin+i/delta)); 4053 if (i==kmax) printf(" +inf "); 4054 else printf(" %10g",exp(lmin+(i+1)/delta)); 4055 printf("| %10g |\n",((long) ((10000.0 * histo[i])/ nbedges))/100.0); 4056 printf(" %i\n",histo[i]); 4057 } 4058 printf(" --------------------+-------------+-------------\n"); 4059 } 4060 /*}}}1*/ 4061 /*FUNCTION Mesh::SmoothingVertex{{{1*/ 4062 void Mesh::SmoothingVertex(int nbiter,double omega ) { 4063 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/ 4064 4065 long int verbose=0; 4066 // if quatree exist remove it end reconstruct 4067 if (quadtree) delete quadtree; 4068 quadtree=0; 4069 ReMakeTriangleContainingTheVertex(); 4070 Triangle vide; // a triangle to mark the boundary vertex 4071 Triangle ** tstart= new Triangle* [nbv]; 4072 long i,j,k; 4073 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide 4074 if ( this == & BTh) 4075 for ( i=0;i<nbv;i++) 4076 tstart[i]=vertices[i].t; 4077 else 4078 for ( i=0;i<nbv;i++) 4079 tstart[i]=0; 4080 for ( j=0;j<NbVerticesOnGeomVertex;j++ ) 4081 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide; 4082 for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 4083 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide; 4084 if(verbose>2) printf(" SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega); 4085 for (k=0;k<nbiter;k++) 4086 { 4087 long i,NbSwap =0; 4088 double delta =0; 4089 for ( i=0;i<nbv;i++) 4090 if (tstart[i] != &vide) // not a boundary vertex 4091 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega)); 4092 if (!NbOfQuad) 4093 for ( i=0;i<nbv;i++) 4094 if (tstart[i] != &vide) // not a boundary vertex 4095 NbSwap += vertices[i].Optim(1); 4096 if (verbose>3) printf(" move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap); 4097 } 4098 4099 delete [] tstart; 4100 if (quadtree) quadtree= new QuadTree(this); 4101 } 4102 /*}}}1*/ 4103 /*FUNCTION Mesh::SmoothMetric{{{1*/ 4104 void Mesh::SmoothMetric(double raisonmax) { 4105 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/ 4106 4107 long int verbose=0; 4108 4109 if(raisonmax<1.1) return; 4110 if(verbose > 1) printf(" Mesh::SmoothMetric raisonmax = %g\n",raisonmax); 4111 ReMakeTriangleContainingTheVertex(); 4112 long i,j,kch,kk,ip; 4113 long *first_np_or_next_t0 = new long[nbv]; 4114 long *first_np_or_next_t1 = new long[nbv]; 4115 long Head0 =0,Head1=-1; 4116 double logseuil= log(raisonmax); 4117 4118 for(i=0;i<nbv-1;i++) 4119 first_np_or_next_t0[i]=i+1; 4120 first_np_or_next_t0[nbv-1]=-1;// end; 4121 for(i=0;i<nbv;i++) 4122 first_np_or_next_t1[i]=-1; 4123 kk=0; 4124 while (Head0>=0&& kk++<100) { 4125 kch=0; 4126 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) { 4127 // pour tous les triangles autour du sommet s 4128 register Triangle* t= vertices[i].t; 4129 if (!t){ 4130 ISSMERROR("!t"); 4131 } 4132 BamgVertex & vi = vertices[i]; 4133 TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]); 4134 BamgVertex *pvj0 = ta.EdgeVertex(0); 4135 while (1) { 4136 ta=Previous(Adj(ta)); 4137 if (vertices+i != ta.EdgeVertex(1)){ 4138 ISSMERROR("vertices+i != ta.EdgeVertex(1)"); 4139 } 4140 BamgVertex & vj = *(ta.EdgeVertex(0)); 4141 if ( &vj ) { 4142 j= &vj-vertices; 4143 if (j<0 || j >= nbv){ 4144 ISSMERROR("j<0 || j >= nbv"); 4145 } 4146 R2 Aij = (R2) vj - (R2) vi; 4147 double ll = Norme2(Aij); 4148 if (0) { 4149 double hi = ll/vi.m(Aij); 4150 double hj = ll/vj.m(Aij); 4151 if(hi < hj) 4152 { 4153 double dh=(hj-hi)/ll; 4154 if (dh>logseuil) { 4155 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); 4156 if(first_np_or_next_t1[j]<0) 4157 kch++,first_np_or_next_t1[j]=Head1,Head1=j; 4158 } 4159 } 4160 } 4161 else 4162 { 4163 double li = vi.m(Aij); 4164 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) 4165 if(first_np_or_next_t1[j]<0) // if the metrix change 4166 kch++,first_np_or_next_t1[j]=Head1,Head1=j; 4167 } 4168 } 4169 if ( &vj == pvj0 ) break; 4170 } 4171 } 4172 Head0 = Head1; 4173 Head1 = -1; 4174 Exchange(first_np_or_next_t0,first_np_or_next_t1); 4175 } 4176 if(verbose>2) printf(" number of iterations = %i\n",kch); 4177 delete [] first_np_or_next_t0; 4178 delete [] first_np_or_next_t1; 4179 } 4180 /*}}}1*/ 4181 /*FUNCTION Mesh::SplitElement{{{1*/ 4182 int Mesh::SplitElement(int choice){ 4183 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/ 4184 4185 long int verbose=0; 4186 4187 Direction NoDirOfSearch; 4188 const int withBackground = &BTh != this && &BTh; 4189 4190 ReNumberingTheTriangleBySubDomain(); 4191 //int nswap =0; 4192 const long nfortria( choice ? 4 : 6); 4193 if(withBackground) 4194 { 4195 BTh.SetVertexFieldOn(); 4196 SetVertexFieldOnBTh(); 4197 } 4198 else 4199 BTh.SetVertexFieldOn(); 4200 4201 long newnbt=0,newnbv=0; 4202 long * kedge = 0; 4203 long newNbOfQuad=NbOfQuad; 4204 long * ksplit= 0, * ksplitarray=0; 4205 long kkk=0; 4206 int ret =0; 4207 if (maxnbv<nbv+nbe) return 1;// 4208 // 1) create the new points by spliting the internal edges 4209 // set the 4210 long nbvold = nbv; 4211 long nbtold = nbt; 4212 long NbOutTold = NbOutT; 4213 long NbEdgeOnGeom=0; 4214 long i; 4215 4216 nbt = nbt - NbOutT; // remove all the the ouside triangles 4217 long nbtsave = nbt; 4218 Triangle * lastT = triangles + nbt; 4219 for (i=0;i<nbe;i++) 4220 if(edges[i].onGeometry) NbEdgeOnGeom++; 4221 long newnbe=nbe+nbe; 4222 // long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex; 4223 long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom; 4224 // long newNbVertexOnBThVertex=NbVertexOnBThVertex; 4225 long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0; 4226 4227 // do allocation for pointeur to the geometry and background 4228 VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge]; 4229 VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ? new VertexOnEdge[newNbVertexOnBThEdge]:0; 4230 if (NbVerticesOnGeomEdge) 4231 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge); 4232 if (NbVertexOnBThEdge) 4233 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge); 4234 Edge *newedges = new Edge [newnbe]; 4235 // memcpy(newedges,edges,sizeof(Edge)*nbe); 4236 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv); 4237 long k=nbv; 4238 long kk=0; 4239 long kvb = NbVertexOnBThEdge; 4240 long kvg = NbVerticesOnGeomEdge; 4241 long ie =0; 4242 Edge ** edgesGtoB=0; 4243 if (withBackground) 4244 edgesGtoB= BTh.MakeGeometricalEdgeToEdge(); 4245 long ferr=0; 4246 for (i=0;i<nbe;i++) 4247 newedges[ie].onGeometry=0; 4248 4249 for (i=0;i<nbe;i++) 4250 { 4251 GeometricalEdge *ong = edges[i].onGeometry; 4252 4253 newedges[ie]=edges[i]; 4254 newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ; 4255 newedges[ie].adj[1]=newedges + ie +1; 4256 R2 A = edges[i][0],B = edges[i][1]; 4257 4258 4259 kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 4260 if (ong) // a geometrical edges 4261 { 4262 if (withBackground){ 4263 // walk on back ground mesh 4264 // newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 4265 // a faire -- difficile 4266 // the first PB is to now a background edge between the 2 vertices 4267 if (!edgesGtoB){ 4268 ISSMERROR("!edgesGtoB"); 4269 } 4270 ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)], 4271 edges[i][0],edges[i][1],0.5,vertices[k], 4272 newVertexOnBThEdge[kvb], 4273 newVerticesOnGeomEdge[kvg++]); 4274 vertices[k].ReferenceNumber= edges[i].ref; 4275 vertices[k].DirOfSearch = NoDirOfSearch; 4276 ; 4277 // get the Info on background mesh 4278 double s = newVertexOnBThEdge[kvb]; 4279 BamgVertex & bv0 = newVertexOnBThEdge[kvb][0]; 4280 BamgVertex & bv1 = newVertexOnBThEdge[kvb][1]; 4281 // compute the metrix of the new points 4282 vertices[k].m = Metric(1-s,bv0,s,bv1); 4283 kvb++; 4284 } 4285 else 4286 { 4287 ong=Gh.ProjectOnCurve(edges[i], 4288 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]); 4289 // vertices[k].i = toI2( vertices[k].r); 4290 vertices[k].ReferenceNumber = edges[i].ref; 4291 vertices[k].DirOfSearch = NoDirOfSearch; 4292 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); 4293 } 4294 } 4295 else // straigth line edge --- 4296 { 4297 vertices[k].r = ((R2) edges[i][0] + (R2) edges[i][1] )*0.5; 4298 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); 4299 vertices[k].onGeometry = 0; 4300 } 4301 //vertices[k].i = toI2( vertices[k].r); 4302 R2 AB = vertices[k].r; 4303 R2 AA = (A+AB)*0.5; 4304 R2 BB = (AB+B)*0.5; 4305 vertices[k].ReferenceNumber = edges[i].ref; 4306 vertices[k].DirOfSearch = NoDirOfSearch; 4307 4308 newedges[ie].onGeometry = Gh.Containing(AA,ong); 4309 newedges[ie++].v[1]=vertices+k; 4310 4311 newedges[ie]=edges[i]; 4312 newedges[ie].adj[0]=newedges + ie -1; 4313 newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ; 4314 newedges[ie].onGeometry = Gh.Containing(BB,ong); 4315 newedges[ie++].v[0]=vertices+k; 4316 k++; 4317 } 4318 if (edgesGtoB) delete [] edgesGtoB; 4319 edgesGtoB=0; 4320 4321 newnbv=k; 4322 newNbVerticesOnGeomEdge=kvg; 4323 if (newnbv> maxnbv) goto Error;// bug 4324 4325 nbv = k; 4326 4327 4328 kedge = new long[3*nbt+1]; 4329 ksplitarray = new long[nbt+1]; 4330 ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0] 4331 4332 for (i=0;i<3*nbt;i++) 4333 kedge[i]=-1; 4334 4335 // 4336 4337 for (i=0;i<nbt;i++) { 4338 Triangle & t = triangles[i]; 4339 if (!t.link){ 4340 ISSMERROR("!t.link"); 4341 } 4342 for(int j=0;j<3;j++) 4343 { 4344 const TriangleAdjacent ta = t.Adj(j); 4345 const Triangle & tt = ta; 4346 if (&tt >= lastT) 4347 t.SetAdj2(j,0,0);// unset adj 4348 const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 4349 const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]]; 4350 long ke =edge4->SortAndFind(Number(v0),Number(v1)); 4351 if (ke>0) 4352 { 4353 long ii = Number(tt); 4354 int jj = ta; 4355 long ks = ke + nbvold; 4356 kedge[3*i+j] = ks; 4357 if (ii<nbt) // good triangle 4358 kedge[3*ii+jj] = ks; 4359 BamgVertex &A=vertices[ks]; 4360 double aa,bb,cc,dd; 4361 if ((dd=Area2(v0.r,v1.r,A.r)) >=0){ 4362 // warning PB roundoff error 4363 if (t.link && ( (aa=Area2( A.r , t[1].r , t[2].r )) < 0.0 4364 || (bb=Area2( t[0].r , A.r , t[2].r )) < 0.0 4365 || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)){ 4366 printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd); 4367 ISSMERROR("Number of triangles with P2 interpolation Problem"); 4368 } 4369 } 4370 else { 4371 if (tt.link && ( (aa=Area2( A.r , tt[1].r , tt[2].r )) < 0 4372 || (bb=Area2( tt[0].r , A.r , tt[2].r )) < 0 4373 || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)){ 4374 printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd); 4375 ISSMERROR("Number of triangles with P2 interpolation Problem"); 4376 } 4377 } 4378 } 4379 } 4380 } 4381 4382 for (i=0;i<nbt;i++){ 4383 ksplit[i]=1; // no split by default 4384 const Triangle & t = triangles[ i]; 4385 int nbsplitedge =0; 4386 int nbinvisible =0; 4387 int invisibleedge=0; 4388 int kkk[3]; 4389 for (int j=0;j<3;j++) 4390 { 4391 if (t.Hidden(j)) invisibleedge=j,nbinvisible++; 4392 4393 const TriangleAdjacent ta = t.Adj(j); 4394 const Triangle & tt = ta; 4395 4396 4397 const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 4398 const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]]; 4399 if ( kedge[3*i+j] < 0) 4400 { 4401 long ke =edge4->SortAndFind(Number(v0),Number(v1)); 4402 if (ke<0) // new 4403 { 4404 if (&tt) // internal triangles all the boundary 4405 { // new internal edges 4406 long ii = Number(tt); 4407 int jj = ta; 4408 4409 kedge[3*i+j]=k;// save the vertex number 4410 kedge[3*ii+jj]=k; 4411 if (k<maxnbv) 4412 { 4413 vertices[k].r = ((R2) v0+(R2) v1 )/2; 4414 //vertices[k].i = toI2( vertices[k].r); 4415 vertices[k].ReferenceNumber=0; 4416 vertices[k].DirOfSearch =NoDirOfSearch; 4417 vertices[k].m = Metric(0.5,v0,0.5,v1); 4418 } 4419 k++; 4420 kkk[nbsplitedge++]=j; 4421 } // tt 4422 else 4423 ISSMERROR("Bug..."); 4424 } // ke<0 4425 else 4426 { // ke >=0 4427 kedge[3*i+j]=nbvold+ke; 4428 kkk[nbsplitedge++]=j;// previously splited 4429 } 4430 } 4431 else 4432 kkk[nbsplitedge++]=j;// previously splited 4433 4434 } 4435 if (nbinvisible>=2){ 4436 ISSMERROR("nbinvisible>=2"); 4437 } 4438 switch (nbsplitedge) { 4439 case 0: ksplit[i]=10; newnbt++; break; // nosplit 4440 case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 4441 case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 4442 case 3: 4443 if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4; 4444 else ksplit[i]=10*nfortria,newnbt+=nfortria; 4445 break; 4446 } 4447 if (ksplit[i]<40){ 4448 ISSMERROR("ksplit[i]<40"); 4449 } 4450 } 4451 // now do the element split 4452 newNbOfQuad = 4*NbOfQuad; 4453 nbv = k; 4454 kkk = nbt; 4455 ksplit[-1] = nbt; 4456 // look on old true triangles 4457 4458 for (i=0;i<nbtsave;i++){ 4459 int nbmkadj=0; 4460 long mkadj [100]; 4461 mkadj[0]=i; 4462 long kk=ksplit[i]/10; 4463 int ke=(int) (ksplit[i]%10); 4464 if (kk>=7 || kk<=0){ 4465 ISSMERROR("kk>=7 || kk<=0"); 4466 } 4467 4468 // def the numbering k (edge) i vertex 4469 int k0 = ke; 4470 int k1 = NextEdge[k0]; 4471 int k2 = PreviousEdge[k0]; 4472 int i0 = OppositeVertex[k0]; 4473 int i1 = OppositeVertex[k1]; 4474 int i2 = OppositeVertex[k2]; 4475 4476 Triangle &t0=triangles[i]; 4477 BamgVertex * v0=t0(i0); 4478 BamgVertex * v1=t0(i1); 4479 BamgVertex * v2=t0(i2); 4480 4481 if (nbmkadj>=10){ 4482 ISSMERROR("nbmkadj>=10"); 4483 } 4484 // -------------------------- 4485 TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2)); 4486 // save the flag Hidden 4487 int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)}; 4488 // un set all adj -- save Hidden flag -- 4489 t0.SetAdj2(0,0,hid[0]); 4490 t0.SetAdj2(1,0,hid[1]); 4491 t0.SetAdj2(2,0,hid[2]); 4492 // -- remake 4493 switch (kk) { 4494 case 1: break;// nothing 4495 case 2: // 4496 { 4497 Triangle &t1=triangles[kkk++]; 4498 t1=t0; 4499 if (kedge[3*i+i0]<0){ 4500 ISSMERROR("kedge[3*i+i0]<0"); 4501 } 4502 BamgVertex * v3 = vertices + kedge[3*i+k0]; 4503 4504 t0(i2) = v3; 4505 t1(i1) = v3; 4506 t0.SetAllFlag(k2,0); 4507 t1.SetAllFlag(k1,0); 4508 } 4509 break; 4510 case 3: // 4511 { 4512 Triangle &t1=triangles[kkk++]; 4513 Triangle &t2=triangles[kkk++]; 4514 t2=t1=t0; 4515 if (kedge[3*i+k1]<0){ 4516 ISSMERROR("kedge[3*i+k1]<0"); 4517 } 4518 if (kedge[3*i+k2]<0){ 4519 ISSMERROR("kedge[3*i+k2]<0"); 4520 } 4521 4522 BamgVertex * v01 = vertices + kedge[3*i+k2]; 4523 BamgVertex * v02 = vertices + kedge[3*i+k1]; 4524 t0(i1) = v01; 4525 t0(i2) = v02; 4526 t1(i2) = v02; 4527 t1(i0) = v01; 4528 t2(i0) = v02; 4529 t0.SetAllFlag(k0,0); 4530 t1.SetAllFlag(k1,0); 4531 t1.SetAllFlag(k0,0); 4532 t2.SetAllFlag(k2,0); 4533 } 4534 break; 4535 case 4: // 4536 case 6: // split in 4 4537 { 4538 Triangle &t1=triangles[kkk++]; 4539 Triangle &t2=triangles[kkk++]; 4540 Triangle &t3=triangles[kkk++]; 4541 t3=t2=t1=t0; 4542 if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){ 4543 ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0"); 4544 } 4545 BamgVertex * v12 = vertices + kedge[3*i+k0]; 4546 BamgVertex * v02 = vertices + kedge[3*i+k1]; 4547 BamgVertex * v01 = vertices + kedge[3*i+k2]; 4548 t0(i1) = v01; 4549 t0(i2) = v02; 4550 t0.SetAllFlag(k0,hid[k0]); 4551 4552 t1(i0) = v01; 4553 t1(i2) = v12; 4554 t0.SetAllFlag(k1,hid[k1]); 4555 4556 t2(i0) = v02; 4557 t2(i1) = v12; 4558 t2.SetAllFlag(k2,hid[k2]); 4559 4560 t3(i0) = v12; 4561 t3(i1) = v02; 4562 t3(i2) = v01; 4563 4564 t3.SetAllFlag(0,hid[0]); 4565 t3.SetAllFlag(1,hid[1]); 4566 t3.SetAllFlag(2,hid[2]); 4567 4568 if ( kk == 6) 4569 { 4570 4571 Triangle &t4=triangles[kkk++]; 4572 Triangle &t5=triangles[kkk++]; 4573 4574 t4 = t3; 4575 t5 = t3; 4576 4577 t0.SetHidden(k0); 4578 t1.SetHidden(k1); 4579 t2.SetHidden(k2); 4580 t3.SetHidden(0); 4581 t4.SetHidden(1); 4582 t5.SetHidden(2); 4583 4584 if (nbv < maxnbv ) 4585 { 4586 vertices[nbv].r = ((R2) *v01 + (R2) *v12 + (R2) *v02 ) / 3.0; 4587 vertices[nbv].ReferenceNumber =0; 4588 vertices[nbv].DirOfSearch =NoDirOfSearch; 4589 //vertices[nbv].i = toI2(vertices[nbv].r); 4590 double a3[]={1./3.,1./3.,1./3.}; 4591 vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m); 4592 BamgVertex * vc = vertices +nbv++; 4593 t3(i0) = vc; 4594 t4(i1) = vc; 4595 t5(i2) = vc; 4596 4597 } 4598 else 4599 goto Error; 4600 } 4601 4602 } 4603 break; 4604 } 4605 4606 // save all the new triangles 4607 mkadj[nbmkadj++]=i; 4608 long jj; 4609 if (t0.link) 4610 for (jj=nbt;jj<kkk;jj++) 4611 { 4612 triangles[jj].link=t0.link; 4613 t0.link= triangles+jj; 4614 mkadj[nbmkadj++]=jj; 4615 } 4616 if (nbmkadj>13){// 13 = 6 + 4 + 4617 ISSMERROR("nbmkadj>13"); 4618 } 4619 4620 if (kk==6) newNbOfQuad+=3; 4621 for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk; 4622 ksplit[i]= nbt; // save last adresse of the new triangles 4623 kkk = nbt; 4624 } 4625 4626 for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.; 4627 4628 if(withBackground) 4629 for (i=0;i<BTh.nbv;i++) 4630 BTh.vertices[i].m = BTh.vertices[i].m*2.; 4631 4632 4633 ret = 2; 4634 if (nbt>= nbtx) goto Error; // bug 4635 if (nbv>= maxnbv) goto Error; // bug 4636 // generation of the new triangles 4637 4638 SetIntCoor("In SplitElement"); 4639 4640 ReMakeTriangleContainingTheVertex(); 4641 if(withBackground) 4642 BTh.ReMakeTriangleContainingTheVertex(); 4643 4644 delete [] edges; 4645 edges = newedges; 4646 nbe = newnbe; 4647 NbOfQuad = newNbOfQuad; 4648 4649 for (i=0;i<NbSubDomains;i++) 4650 { 4651 long k = subdomains[i].edge- edges; 4652 subdomains[i].edge = edges+2*k; // spilt all edge in 2 4653 } 4654 4655 if (ksplitarray) delete [] ksplitarray; 4656 if (kedge) delete [] kedge; 4657 if (edge4) delete edge4; 4658 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge; 4659 VerticesOnGeomEdge= newVerticesOnGeomEdge; 4660 if(VertexOnBThEdge) delete [] VertexOnBThEdge; 4661 VertexOnBThEdge = newVertexOnBThEdge; 4662 NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge; 4663 NbVertexOnBThEdge=newNbVertexOnBThEdge; 4664 // ReMakeTriangleContainingTheVertex(); 4665 4666 ReconstructExistingMesh(); 4667 4668 if (verbose>2){ 4669 printf(" number of quadrilaterals = %i\n",NbOfQuad); 4670 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); 4671 printf(" number of outside triangles = %i\n",NbOutT); 4672 } 4673 4674 return 0; //ok 4675 4676 Error: 4677 nbv = nbvold; 4678 nbt = nbtold; 4679 NbOutT = NbOutTold; 4680 // cleaning memory --- 4681 delete newedges; 4682 if (ksplitarray) delete [] ksplitarray; 4683 if (kedge) delete [] kedge; 4684 if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge; 4685 if (edge4) delete edge4; 4686 if(newVertexOnBThEdge) delete [] newVertexOnBThEdge; 4687 4688 return ret; // ok 4689 } 4690 /*}}}1*/ 4691 /*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/ 4692 long Mesh::SplitInternalEdgeWithBorderVertices(){ 4693 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/ 4694 4695 long NbSplitEdge=0; 4696 SetVertexFieldOn(); 4697 long it; 4698 long nbvold=nbv; 4699 long int verbose=2; 4700 for (it=0;it<nbt;it++){ 4701 Triangle &t=triangles[it]; 4702 if (t.link) 4703 for (int j=0;j<3;j++) 4704 if(!t.Locked(j) && !t.Hidden(j)){ 4705 Triangle &tt = *t.TriangleAdj(j); 4706 if ( &tt && tt.link && it < Number(tt)) 4707 { // an internal edge 4708 BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 4709 BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 4710 if (v0.onGeometry && v1.onGeometry){ 4711 R2 P= ((R2) v0 + (R2) v1)*0.5; 4712 if ( nbv<maxnbv) { 4713 vertices[nbv].r = P; 4714 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m); 4715 vertices[nbv].ReferenceNumber=0; 4716 vertices[nbv].DirOfSearch = NoDirOfSearch ; 4717 } 4718 NbSplitEdge++; 4719 } 4720 } 4721 } 4722 } 4723 ReMakeTriangleContainingTheVertex(); 4724 if (nbvold!=nbv){ 4725 long iv = nbvold; 4726 long NbSwap = 0; 4727 Icoor2 dete[3]; 4728 for (int i=nbvold;i<nbv;i++) {// for all the new point 4729 BamgVertex & vi = vertices[i]; 4730 vi.i = toI2(vi.r); 4731 vi.r = toR2(vi.i); 4732 4733 // a good new point 4734 vi.ReferenceNumber=0; 4735 vi.DirOfSearch =NoDirOfSearch; 4736 Triangle *tcvi = FindTriangleContaining(vi.i,dete); 4737 if (tcvi && !tcvi->link) { 4738 printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n"); 4739 } 4740 4741 quadtree->Add(vi); 4742 if (!tcvi || tcvi->det<0){// internal 4743 ISSMERROR("!tcvi || tcvi->det < 0"); 4744 } 4745 AddVertex(vi,tcvi,dete); 4746 NbSwap += vi.Optim(1); 4747 iv++; 4748 // } 4749 } 4750 if (verbose>3) { 4751 printf(" number of points: %i\n",iv); 4752 printf(" number of swap to split internal edges with border vertices: %i\n",NbSwap); 4753 nbv = iv; 4754 } 4755 } 4756 if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold)); 4757 if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge); 4758 4759 return NbSplitEdge; 4760 } 4761 /*}}}1*/ 4762 /*FUNCTION Mesh::TriangulateFromGeom0{{{1*/ 4763 void Mesh::TriangulateFromGeom0(long imaxnbv,BamgOpts* bamgopts){ 2726 4764 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/ 2727 4765 … … 2737 4775 R2 AB; 2738 4776 GeometricalVertex *a,*b; 2739 MeshVertex *va,*vb;4777 BamgVertex *va,*vb; 2740 4778 GeometricalEdge *e; 2741 4779 … … 2745 4783 //initialize this 2746 4784 if (verbose>3) printf(" Generating Boundary vertices\n"); 2747 PreInit(inbvx);4785 Init(imaxnbv); 2748 4786 nbv=0; 2749 4787 NbVerticesOnGeomVertex=0; … … 2760 4798 //initialize VerticesOnGeomVertex 2761 4799 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 2762 if( NbVerticesOnGeomVertex >= nbvx) {2763 ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex, nbvx);4800 if( NbVerticesOnGeomVertex >= maxnbv) { 4801 ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv); 2764 4802 } 2765 4803 … … 2773 4811 if (Gh[i].Required() && Gh[i].IsThe()) {//Gh vertices Required 2774 4812 2775 //Add the vertex (provided that nbv< nbvx)2776 if (nbv< nbvx){4813 //Add the vertex (provided that nbv<maxnbv) 4814 if (nbv<maxnbv){ 2777 4815 vertices[nbv]=Gh[i]; 2778 4816 } 2779 4817 else{ 2780 ISSMERROR("Maximum number of vertices ( nbvx = %i) too small",nbvx);4818 ISSMERROR("Maximum number of vertices (maxnbv = %i) too small",maxnbv); 2781 4819 } 2782 4820 … … 2860 4898 NbNewPoints=0; 2861 4899 NbEdgeCurve=0; 2862 if (nbvend>= nbvx) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase nbvx");4900 if (nbvend>=maxnbv) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase maxnbv"); 2863 4901 lcurve =0; 2864 4902 s = lstep; … … 3033 5071 } 3034 5072 /*}}}1*/ 3035 /*FUNCTION Mesh:: GeomToTriangles1{{{1*/3036 void Mesh:: GeomToTriangles1(long inbvx,BamgOpts* bamgopts,int KeepVertices){5073 /*FUNCTION Mesh::TriangulateFromGeom1{{{1*/ 5074 void Mesh::TriangulateFromGeom1(long imaxnbv,BamgOpts* bamgopts,int KeepVertices){ 3037 5075 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles1)*/ 3038 5076 … … 3064 5102 3065 5103 //Initialize new mesh 3066 this-> PreInit(inbvx);5104 this->Init(imaxnbv); 3067 5105 BTh.SetVertexFieldOn(); 3068 5106 int* bcurve = new int[Gh.NbOfCurves]; // … … 3082 5120 int i; 3083 5121 for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++; 3084 if(NbVerticesOnGeomVertex >= nbvx) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);}5122 if(NbVerticesOnGeomVertex >= maxnbv) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);} 3085 5123 3086 5124 VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex]; 3087 5125 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex]; 3088 5126 3089 //At this point there is NO vertex but vertices should have been allocated by PreInit5127 //At this point there is NO vertex but vertices should have been allocated by Init 3090 5128 ISSMASSERT(vertices); 3091 5129 for (i=0;i<Gh.nbv;i++){ … … 3103 5141 if (vog.IsRequiredVertex()){ 3104 5142 GeometricalVertex* gv=vog; 3105 MeshVertex *bv = vog;5143 BamgVertex *bv = vog; 3106 5144 ISSMASSERT(gv->to); // use of Geom -> Th 3107 5145 VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->to,bv); … … 3208 5246 long i=0;// index of new points on the curve 3209 5247 register GeometricalVertex * GA0 = *(*peequi)[k0equi].onGeometry; 3210 MeshVertex *A0;5248 BamgVertex *A0; 3211 5249 A0 = GA0->to; // the vertex in new mesh 3212 MeshVertex *A1;5250 BamgVertex *A1; 3213 5251 VertexOnGeom *GA1; 3214 5252 Edge* PreviousNewEdge = 0; … … 3228 5266 ISSMASSERT(pe && ee.onGeometry); 3229 5267 ee.onGeometry->SetMark(); 3230 MeshVertex & v0=ee[0], & v1=ee[1];5268 BamgVertex & v0=ee[0], & v1=ee[1]; 3231 5269 R2 AB=(R2)v1-(R2)v0; 3232 5270 double L0=L,LAB; … … 3241 5279 ISSMASSERT(sNew>=L0); 3242 5280 ISSMASSERT(LAB); 3243 ISSMASSERT(vertices && nbv< nbvx);5281 ISSMASSERT(vertices && nbv<maxnbv); 3244 5282 ISSMASSERT(edges && nbe<nbex); 3245 5283 ISSMASSERT(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex); … … 3328 5366 //Allocate memory 3329 5367 if(step==0){ 3330 if(nbv+NbOfNewPoints > nbvx) {3331 ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints, nbvx);5368 if(nbv+NbOfNewPoints > maxnbv) { 5369 ISSMERROR("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,maxnbv); 3332 5370 } 3333 5371 edges = new Edge[NbOfNewEdge]; … … 3365 5403 } 3366 5404 /*}}}1*/ 3367 /*FUNCTION Mesh::Insert{{{1*/3368 void Mesh::Insert() {3369 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/3370 3371 /*Insert points in the existing Geometry*/3372 3373 //Intermediary3374 int i;3375 3376 /*Get options*/3377 long int verbose=2;3378 3379 //Display info3380 if (verbose>2) printf(" Insert initial %i vertices\n",nbv);3381 3382 //Compute integer coordinates and determinants for the existing vertices (from Geometry)3383 SetIntCoor();3384 3385 /*Now we want to build a list (ordre) of the vertices in a random3386 * order. To do so, we use the following method:3387 *3388 * From an initial k0 in [0 nbv[ random (vertex number)3389 * the next k (vertex number) is computed using a big3390 * prime number (PN>>nbv) following:3391 *3392 * k_{i+1} = k_i + PN [nbv]3393 *3394 * let's show that:3395 *3396 * for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j3397 *3398 * Let's assume that there are 2 distinct j1 and j2 such that3399 * k_j1 = k_j23400 *3401 * This means that3402 *3403 * k0+j1*PN = k0+j2*PN [nbv]3404 * (j1-j2)*PN =0 [nbv]3405 * since PN is a prime number larger than nbv, and nbv!=13406 * j1-j2=0 [nbv]3407 * BUT3408 * j1 and j2 are in [0 nbv[ which is impossible.3409 *3410 * We hence have built a random list of nbv elements of3411 * [0 nbv[ all distincts*/3412 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ;3413 const long PrimeNumber= BigPrimeNumber(nbv) ;3414 int k0=rand()%nbv;3415 for (int is3=0; is3<nbv; is3++){3416 ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv];3417 }3418 3419 /*Modify ordre such that the first 3 vertices form a triangle*/3420 3421 //get first vertex i such that [0,1,i] are not aligned3422 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){3423 //if i is higher than nbv, it means that all the determinants are 0,3424 //all vertices are aligned!3425 if ( ++i >= nbv) {3426 ISSMERROR("all the vertices are aligned");3427 }3428 }3429 // exchange i et 2 in "ordre" so that3430 // the first 3 vertices are not aligned (real triangle)3431 Exchange(ordre[2], ordre[i]);3432 3433 /*Take the first edge formed by the first two vertices and build3434 * the initial simple mesh from this edge and 2 boundary triangles*/3435 3436 MeshVertex * v0=ordre[0], *v1=ordre[1];3437 3438 nbt = 2;3439 triangles[0](0) = NULL; //infinite vertex3440 triangles[0](1) = v0;3441 triangles[0](2) = v1;3442 triangles[1](0) = NULL;//infinite vertex3443 triangles[1](2) = v0;3444 triangles[1](1) = v1;3445 3446 //Build adjacence3447 const int e0 = OppositeEdge[0];3448 const int e1 = NextEdge[e0];3449 const int e2 = PreviousEdge[e0];3450 triangles[0].SetAdj2(e0, &triangles[1] ,e0);3451 triangles[0].SetAdj2(e1, &triangles[1] ,e2);3452 triangles[0].SetAdj2(e2, &triangles[1] ,e1);3453 3454 triangles[0].det = -1; //boundary triangle: det = -13455 triangles[1].det = -1; //boundary triangle: det = -13456 3457 triangles[0].SetTriangleContainingTheVertex();3458 triangles[1].SetTriangleContainingTheVertex();3459 3460 triangles[0].link=&triangles[1];3461 triangles[1].link=&triangles[0];3462 3463 //build quadtree3464 if (!quadtree) quadtree = new QuadTree(this,0);3465 quadtree->Add(*v0);3466 quadtree->Add(*v1);3467 3468 /*Now, add the vertices One by One*/3469 3470 long NbSwap=0;3471 if (verbose>3) printf(" Begining of insertion process...\n");3472 3473 for (int icount=2; icount<nbv; icount++) {3474 3475 //Get new vertex3476 MeshVertex *newvertex=ordre[icount];3477 3478 //Find the triangle in which newvertex is located3479 Icoor2 dete[3];3480 Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)3481 3482 //Add newvertex to the quadtree3483 quadtree->Add(*newvertex);3484 3485 //Add newvertex to the existing mesh3486 AddVertex(*newvertex,tcvi,dete);3487 3488 //Make the mesh Delaunay around newvertex by swaping the edges3489 NbSwap += newvertex->Optim(1,0);3490 }3491 3492 //Display info3493 if (verbose>3) {3494 printf(" NbSwap of insertion: %i\n",NbSwap);3495 printf(" NbSwap/nbv: %i\n",NbSwap/nbv);3496 }3497 3498 #ifdef NBLOOPOPTIM3499 3500 k0 = rand()%nbv ;3501 for (int is4=0; is4<nbv; is4++)3502 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv];3503 3504 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){3505 long NbSwap = 0;3506 for (int is1=0; is1<nbv; is1++)3507 NbSwap += ordre[is1]->Optim(0,0);3508 if (verbose>3) {3509 printf(" Optim Loop: %i\n",Nbloop);3510 printf(" NbSwap/nbv: %i\n",NbSwap/nbv);3511 }3512 if(!NbSwap) break;3513 }3514 ReMakeTriangleContainingTheVertex();3515 // because we break the TriangleContainingTheVertex3516 #endif3517 }3518 /*}}}1*/3519 /*FUNCTION Mesh::InsertNewPoints{{{1*/3520 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {3521 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/3522 3523 long int verbose=0;3524 double seuil= 1.414/2 ;// for two close point3525 long i;3526 long NbSwap=0;3527 Icoor2 dete[3];3528 3529 //number of new points3530 const long nbvnew=nbv-nbvold;3531 3532 //display info if required3533 if (verbose>5) printf(" Try to Insert %i new points\n",nbvnew);3534 3535 //return if no new points3536 if (!nbvnew) return 0;3537 3538 /*construction of a random order*/3539 const long PrimeNumber= BigPrimeNumber(nbv) ;3540 //remainder of the division of rand() by nbvnew3541 long k3 = rand()%nbvnew;3542 //loop over the new points3543 for (int is3=0; is3<nbvnew; is3++){3544 register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);3545 register long i=nbvold+is3;3546 ordre[i]= vertices + j;3547 ordre[i]->ReferenceNumber=i;3548 }3549 3550 // for all the new point3551 long iv=nbvold;3552 for (i=nbvold;i<nbv;i++){3553 MeshVertex &vi=*ordre[i];3554 vi.i=toI2(vi.r);3555 vi.r=toR2(vi.i);3556 double hx,hy;3557 vi.m.Box(hx,hy);3558 Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);3559 if (!quadtree->ToClose(vi,seuil,hi,hj)){3560 // a good new point3561 MeshVertex &vj = vertices[iv];3562 long j=vj.ReferenceNumber;3563 if (&vj!=ordre[j]){3564 ISSMERROR("&vj!= ordre[j]");3565 }3566 if(i!=j){3567 Exchange(vi,vj);3568 Exchange(ordre[j],ordre[i]);3569 }3570 vj.ReferenceNumber=0;3571 Triangle *tcvj=FindTriangleContaining(vj.i,dete);3572 if (tcvj && !tcvj->link){3573 tcvj->Echo();3574 ISSMERROR("problem inserting point in InsertNewPoints (tcvj=%p and tcvj->link=%i)",tcvj,tcvj->link);3575 }3576 quadtree->Add(vj);3577 AddVertex(vj,tcvj,dete);3578 NbSwap += vj.Optim(1);3579 iv++;3580 }3581 }3582 if (verbose>3) {3583 printf(" number of new points: %i\n",iv);3584 printf(" number of to close (?) points: %i\n",nbv-iv);3585 printf(" number of swap after: %i\n",NbSwap);3586 }3587 nbv = iv;3588 3589 for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);3590 if (verbose>3) printf(" NbSwap=%i\n",NbSwap);3591 3592 NbTSwap += NbSwap ;3593 return nbv-nbvold;3594 }3595 /*}}}1*/3596 /*FUNCTION Mesh::MakeGeometricalEdgeToEdge{{{1*/3597 Edge** Mesh::MakeGeometricalEdgeToEdge() {3598 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeometricalEdgeToEdge)*/3599 3600 if (!Gh.nbe){3601 ISSMERROR("!Gh.nbe");3602 }3603 Edge **e= new (Edge* [Gh.nbe]);3604 3605 long i;3606 for ( i=0;i<Gh.nbe ; i++)3607 e[i]=NULL;3608 for ( i=0;i<nbe ; i++)3609 {3610 Edge * ei = edges+i;3611 GeometricalEdge *onGeometry = ei->onGeometry;3612 e[Gh.Number(onGeometry)] = ei;3613 }3614 for ( i=0;i<nbe ; i++)3615 for (int ii=0;ii<2;ii++) {3616 Edge * ei = edges+i;3617 GeometricalEdge *onGeometry = ei->onGeometry;3618 int j= ii;3619 while (!(*onGeometry)[j].Required()) {3620 Adj(onGeometry,j); // next geom edge3621 j=1-j;3622 if (e[Gh.Number(onGeometry)]) break; // optimisation3623 e[Gh.Number(onGeometry)] = ei;3624 }3625 }3626 3627 int kk=0;3628 for ( i=0;i<Gh.nbe ; i++){3629 if (!e[i]){3630 kk++;3631 if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i);3632 }3633 }3634 if(kk) ISSMERROR("See above");3635 3636 return e;3637 }3638 /*}}}1*/3639 /*FUNCTION Mesh::MakeQuadrangles{{{1*/3640 void Mesh::MakeQuadrangles(double costheta){3641 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadrangles)*/3642 3643 long int verbose=0;3644 3645 if (verbose>2) printf("MakeQuadrangles costheta = %g\n",costheta);3646 3647 if (costheta >1) {3648 if (verbose>5) printf(" do nothing: costheta > 1\n");3649 }3650 3651 3652 long nbqq = (nbt*3)/2;3653 DoubleAndInt *qq = new DoubleAndInt[nbqq];3654 3655 long i,ij;3656 int j;3657 long k=0;3658 for (i=0;i<nbt;i++)3659 for (j=0;j<3;j++)3660 if ((qq[k].q=triangles[i].QualityQuad(j))>=costheta)3661 qq[k++].i3j=i*3+j;3662 // sort qq3663 HeapSort(qq,k);3664 3665 long kk=0;3666 for (ij=0;ij<k;ij++) {3667 i=qq[ij].i3j/3;3668 j=(int) (qq[ij].i3j%3);3669 // optisamition no float computation3670 if (triangles[i].QualityQuad(j,0) >=costheta)3671 triangles[i].SetHidden(j),kk++;3672 }3673 NbOfQuad = kk;3674 if (verbose>2){3675 printf(" number of quadrilaterals = %i\n",NbOfQuad);3676 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2);3677 printf(" number of outside triangles = %i\n",NbOutT);3678 }3679 delete [] qq;3680 }3681 /*}}}1*/3682 /*FUNCTION Mesh::MakeQuadTree{{{1*/3683 void Mesh::MakeQuadTree() {3684 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/3685 3686 long int verbose=0;3687 if ( !quadtree ) quadtree = new QuadTree(this);3688 3689 }3690 /*}}}1*/3691 /*FUNCTION Mesh::MaxSubDivision{{{1*/3692 void Mesh::MaxSubDivision(double maxsubdiv) {3693 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/3694 3695 long int verbose=0;3696 3697 const double maxsubdiv2 = maxsubdiv*maxsubdiv;3698 if(verbose>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv);3699 // for all the edges3700 // if the len of the edge is to long3701 long it,nbchange=0;3702 double lmax=0;3703 for (it=0;it<nbt;it++){3704 Triangle &t=triangles[it];3705 for (int j=0;j<3;j++){3706 Triangle &tt = *t.TriangleAdj(j);3707 if ( ! &tt || it < Number(tt) && ( tt.link || t.link)){3708 MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];3709 MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];3710 R2 AB= (R2) v1-(R2) v0;3711 Metric M = v0;3712 double l = M(AB,AB);3713 lmax = Max(lmax,l);3714 if(l> maxsubdiv2){3715 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M3716 double lc = M(AC,AC);3717 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;3718 D2xD2 Rt1(Rt.inv());3719 D2xD2 D(maxsubdiv2,0,0,lc);3720 D2xD2 MM = Rt1*D*Rt1.t();3721 v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);3722 nbchange++;3723 }3724 M = v1;3725 l = M(AB,AB);3726 lmax = Max(lmax,l);3727 if(l> maxsubdiv2){3728 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M3729 double lc = M(AC,AC);3730 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;3731 D2xD2 Rt1(Rt.inv());3732 D2xD2 D(maxsubdiv2,0,0,lc);3733 D2xD2 MM = Rt1*D*Rt1.t();3734 v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);3735 nbchange++;3736 }3737 }3738 }3739 }3740 if(verbose>3){3741 printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5));3742 }3743 }3744 /*}}}1*/3745 /*FUNCTION Mesh::MetricAt{{{1*/3746 Metric Mesh::MetricAt(const R2 & A) const {3747 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/3748 3749 I2 a = toI2(A);3750 Icoor2 deta[3];3751 Triangle * t =FindTriangleContaining(a,deta);3752 if (t->det <0) { // outside3753 double ba,bb;3754 TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ;3755 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}3756 else { // inside3757 double aa[3];3758 double s = deta[0]+deta[1]+deta[2];3759 aa[0]=deta[0]/s;3760 aa[1]=deta[1]/s;3761 aa[2]=deta[2]/s;3762 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);3763 }3764 }3765 /*}}}1*/3766 /*FUNCTION Mesh::NearestVertex{{{1*/3767 MeshVertex* Mesh::NearestVertex(Icoor1 i,Icoor1 j) {3768 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/3769 return quadtree->NearestVertex(i,j);3770 }3771 /*}}}1*/3772 /*FUNCTION Mesh::NewPoints{{{1*/3773 void Mesh::NewPoints(Mesh & Bh,BamgOpts* bamgopts,int KeepVertices){3774 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/3775 3776 int i,j,k;3777 long NbTSwap=0;3778 long nbtold=nbt;3779 long nbvold=nbv;3780 long Headt=0;3781 long next_t;3782 long* first_np_or_next_t=new long[nbtx];3783 Triangle* t=NULL;3784 3785 /*Recover options*/3786 int verbose=bamgopts->verbose;3787 3788 /*First, insert old points if requested*/3789 if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){3790 if (verbose>5) printf(" Inserting initial mesh points\n");3791 for (i=0;i<Bh.nbv;i++){3792 MeshVertex &bv=Bh[i];3793 if (!bv.onGeometry){3794 vertices[nbv].r = bv.r;3795 vertices[nbv++].m = bv.m;3796 }3797 }3798 Bh.ReMakeTriangleContainingTheVertex();3799 InsertNewPoints(nbvold,NbTSwap) ;3800 }3801 else Bh.ReMakeTriangleContainingTheVertex();3802 3803 // generation of the list of next Triangle3804 for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);3805 // the next traingle of i is -first_np_or_next_t[i]3806 3807 // Big loop (most time consuming)3808 int iter=0;3809 if (verbose>5) printf(" Big loop\n");3810 do {3811 /*Update variables*/3812 iter++;3813 nbtold=nbt;3814 nbvold=nbv;3815 3816 /*We test all triangles*/3817 i=Headt;3818 next_t=-first_np_or_next_t[i];3819 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){3820 3821 //check i3822 if (i<0 || i>=nbt ){3823 ISSMERROR("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1);3824 }3825 //change first_np_or_next_t[i]3826 first_np_or_next_t[i] = iter;3827 3828 //Loop over the edges of t3829 for(j=0;j<3;j++){3830 TriangleAdjacent tj(t,j);3831 MeshVertex &vA = *tj.EdgeVertex(0);3832 MeshVertex &vB = *tj.EdgeVertex(1);3833 3834 //if t is a boundary triangle, or tj locked, continue3835 if (!t->link) continue;3836 if (t->det <0) continue;3837 if (t->Locked(j)) continue;3838 3839 TriangleAdjacent tadjj = t->Adj(j);3840 Triangle* ta=tadjj;3841 3842 //if the adjacent triangle is a boundary triangle, continur3843 if (ta->det<0) continue;3844 3845 R2 A=vA;3846 R2 B=vB;3847 k=Number(ta);3848 3849 //if this edge has already been done, go to next edge of triangle3850 if(first_np_or_next_t[k]==iter) continue;3851 3852 lIntTria.SplitEdge(Bh,A,B);3853 lIntTria.NewPoints(vertices,nbv,nbvx);3854 } // end loop for each edge3855 }// for triangle3856 3857 if (!InsertNewPoints(nbvold,NbTSwap)) break;3858 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;3859 Headt = nbt; // empty list3860 3861 // for all the triangle containing the vertex i3862 for (i=nbvold;i<nbv;i++){3863 MeshVertex* s = vertices + i;3864 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);3865 Triangle* tbegin= (Triangle*) ta;3866 long kt;3867 do {3868 kt = Number((Triangle*) ta);3869 if (first_np_or_next_t[kt]>0){3870 first_np_or_next_t[kt]=-Headt;3871 Headt=kt;3872 }3873 if (ta.EdgeVertex(0)!=s){3874 ISSMERROR("ta.EdgeVertex(0)!=s");3875 }3876 ta = Next(Adj(ta));3877 } while ( (tbegin != (Triangle*) ta));3878 }3879 3880 } while (nbv!=nbvold);3881 3882 delete [] first_np_or_next_t;3883 3884 long NbSwapf =0,NbSwp;3885 3886 NbSwp = NbSwapf;3887 for (i=0;i<nbv;i++)3888 NbSwapf += vertices[i].Optim(0);3889 NbTSwap += NbSwapf ;3890 }3891 /*}}}1*/3892 /*FUNCTION Mesh::PreInit{{{1*/3893 void Mesh::PreInit(long inbvx) {3894 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/3895 3896 long int verbose=0;3897 3898 srand(19999999);3899 NbRef=0;3900 // allocGeometry=0;3901 NbOfTriangleSearchFind =0;3902 NbOfSwapTriangle =0;3903 nbiv=0;3904 nbv=0;3905 nbvx=inbvx;3906 nbt=0;3907 NbOfQuad = 0;3908 nbtx=2*inbvx-2;3909 NbSubDomains=0;3910 NbVertexOnBThVertex=0;3911 NbVertexOnBThEdge=0;3912 VertexOnBThVertex=NULL;3913 VertexOnBThEdge=NULL;3914 3915 NbCrackedVertices=0;3916 CrackedVertices =NULL;3917 NbCrackedEdges =0;3918 CrackedEdges =NULL;3919 nbe = 0;3920 3921 if (inbvx) {3922 vertices=new MeshVertex[nbvx];3923 if (!vertices){3924 ISSMERROR("!vertices");3925 }3926 ordre=new (MeshVertex* [nbvx]);3927 if (!ordre){3928 ISSMERROR("!ordre");3929 }3930 triangles=new Triangle[nbtx];3931 if (!triangles){3932 ISSMERROR("!triangles");3933 }3934 }3935 else {3936 vertices=NULL;3937 ordre=NULL;3938 triangles=NULL;3939 nbtx=0;3940 }3941 3942 quadtree=NULL;3943 edges=NULL;3944 VerticesOnGeomVertex=NULL;3945 VerticesOnGeomEdge=NULL;3946 NbVerticesOnGeomVertex=0;3947 NbVerticesOnGeomEdge=0;3948 subdomains=NULL;3949 NbSubDomains=0;3950 }3951 /*}}}1*/3952 /*FUNCTION Mesh::ProjectOnCurve{{{1*/3953 GeometricalEdge* Mesh::ProjectOnCurve( Edge & BhAB, MeshVertex & vA, MeshVertex & vB,3954 double theta,MeshVertex & R,VertexOnEdge & BR,VertexOnGeom & GR) {3955 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/3956 3957 void *pA=0,*pB=0;3958 double tA=0,tB=0;3959 R2 A=vA,B=vB;3960 MeshVertex * pvA=&vA, * pvB=&vB;3961 if (vA.vint == IsVertexOnVertex){3962 pA=vA.onBackgroundVertex;3963 }3964 else if (vA.vint == IsVertexOnEdge){3965 pA=vA.onBackgroundEdge->be;3966 tA=vA.onBackgroundEdge->abcisse;3967 }3968 else {3969 ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA));3970 }3971 3972 if (vB.vint == IsVertexOnVertex){3973 pB=vB.onBackgroundVertex;3974 }3975 else if(vB.vint == IsVertexOnEdge){3976 pB=vB.onBackgroundEdge->be;3977 tB=vB.onBackgroundEdge->abcisse;3978 }3979 else {3980 ISSMERROR("ProjectOnCurve On MeshVertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB));3981 }3982 Edge * e = &BhAB;3983 if (!pA || !pB || !e){3984 ISSMERROR("!pA || !pB || !e");3985 }3986 // be carefull the back ground edge e is on same geom edge3987 // of the initiale edge def by the 2 vertex A B;3988 //check Is a background Mesh;3989 if (e<BTh.edges || e>=BTh.edges+BTh.nbe){3990 ISSMERROR("e<BTh.edges || e>=BTh.edges+BTh.nbe");3991 }3992 // walk on BTh edge3993 //not finish ProjectOnCurve with BackGround Mesh);3994 // 1 first find a back ground edge contening the vertex A3995 // 2 walk n back gound boundary to find the final vertex B3996 3997 if( vA.vint == IsVertexOnEdge)3998 { // find the start edge3999 e = vA.onBackgroundEdge->be;4000 4001 }4002 else if (vB.vint == IsVertexOnEdge)4003 {4004 theta = 1-theta;4005 Exchange(tA,tB);4006 Exchange(pA,pB);4007 Exchange(pvA,pvB);4008 Exchange(A,B);4009 e = vB.onBackgroundEdge->be;4010 4011 }4012 else{ // do the search by walking4013 ISSMERROR("case not supported yet");4014 }4015 4016 // find the direction of walking with sens of edge and pA,PB;4017 R2 AB=B-A;4018 4019 double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);4020 int kkk=0;4021 int sens = (cosE01AB>0) ? 1 : 0;4022 4023 // double l=0; // length of the edge AB4024 double abscisse = -1;4025 4026 for (int cas=0;cas<2;cas++){4027 // 2 times algo:4028 // 1 for computing the length l4029 // 2 for find the vertex4030 int iii;4031 MeshVertex *v0=pvA,*v1;4032 Edge *neee,*eee;4033 double lg =0; // length of the curve4034 double te0;4035 // we suppose take the curve's abcisse4036 for ( eee=e,iii=sens,te0=tA;4037 eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;4038 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {4039 4040 kkk=kkk+1;4041 if (kkk>=100){4042 ISSMERROR("kkk>=100");4043 }4044 if (!eee){4045 ISSMERROR("!eee");4046 }4047 double lg0 = lg;4048 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);4049 lg += dp;4050 if (cas && abscisse <= lg) { // ok we find the geom edge4051 double sss = (abscisse-lg0)/dp;4052 double thetab = te0*(1-sss)+ sss*iii;4053 if (thetab<0 || thetab>1){4054 ISSMERROR("thetab<0 || thetab>1");4055 }4056 BR = VertexOnEdge(&R,eee,thetab);4057 return Gh.ProjectOnCurve(*eee,thetab,R,GR);4058 }4059 }4060 // we find the end4061 if (v1 != pvB){4062 if (( void*) v1 == pB)4063 tB = iii;4064 4065 double lg0 = lg;4066 if (!eee){4067 ISSMERROR("!eee");4068 }4069 v1 = pvB;4070 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);4071 lg += dp;4072 abscisse = lg*theta;4073 if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end4074 { // ok we find the geom edge4075 double sss = (abscisse-lg0)/dp;4076 double thetab = te0*(1-sss)+ sss*tB;4077 if (thetab<0 || thetab>1){4078 ISSMERROR("thetab<0 || thetab>1");4079 }4080 BR = VertexOnEdge(&R,eee,thetab);4081 return Gh.ProjectOnCurve(*eee,thetab,R,GR);4082 }4083 }4084 abscisse = lg*theta;4085 4086 }4087 ISSMERROR("Big bug...");4088 return 0; // just for the compiler4089 }4090 /*}}}1*/4091 /*FUNCTION Mesh::ReconstructExistingMesh{{{1*/4092 void Mesh::ReconstructExistingMesh(){4093 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/4094 4095 /*This routine reconstruct an existing mesh to make it CONVEX:4096 * -all the holes are filled4097 * -concave boundaries are filled4098 * A convex mesh is required for a lot of operations. This is why every mesh4099 * goes through this process.4100 * This routine also generates mesh properties such as adjencies,...4101 */4102 4103 /*Intermediary*/4104 int verbose=0;4105 4106 // generation of the integer coordinate4107 4108 // find extrema coordinates of vertices pmin,pmax4109 long i;4110 if(verbose>2) printf(" Reconstruct mesh of %i vertices\n",nbv);4111 4112 //initialize ordre4113 ISSMASSERT(ordre);4114 for (i=0;i<nbv;i++) ordre[i]=0;4115 4116 //Initialize NbSubDomains4117 NbSubDomains =0;4118 4119 /* generation of triangles adjacency*/4120 4121 //First add existing edges4122 long kk =0;4123 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);4124 for (i=0;i<nbe;i++){4125 kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));4126 }4127 if (kk != nbe){4128 ISSMERROR("There are %i double edges in the mesh",kk-nbe);4129 }4130 4131 //Add edges of all triangles in existing mesh4132 long* st = new long[nbt*3];4133 for (i=0;i<nbt*3;i++) st[i]=-1;4134 for (i=0;i<nbt;i++){4135 for (int j=0;j<3;j++){4136 4137 //Add current triangle edge to edge44138 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));4139 4140 long invisible=triangles[i].Hidden(j);4141 4142 //If the edge has not been added to st, add it4143 if(st[k]==-1) st[k]=3*i+j;4144 4145 //If the edge already exists, add adjacency4146 else if(st[k]>=0) {4147 ISSMASSERT(!triangles[i].TriangleAdj(j));4148 ISSMASSERT(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));4149 4150 triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));4151 if (invisible) triangles[i].SetHidden(j);4152 if (k<nbe) triangles[i].SetLocked(j);4153 4154 //Make st[k] negative so that it will throw an error message if it is found again4155 st[k]=-2-st[k];4156 }4157 4158 //An edge belongs to 2 triangles4159 else {4160 ISSMERROR("The edge (%i , %i) belongs to more than 2 triangles",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));4161 }4162 }4163 }4164 4165 //Display info if required4166 if(verbose>5) {4167 printf(" info of Mesh:\n");4168 printf(" - number of vertices = %i \n",nbv);4169 printf(" - number of triangles = %i \n",nbt);4170 printf(" - number of given edges = %i \n",nbe);4171 printf(" - number of all edges = %i \n" ,edge4->nb());4172 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv);4173 }4174 4175 //check the consistency of edge[].adj and the geometrical required vertex4176 long k=0;4177 for (i=0;i<edge4->nb();i++){4178 if (st[i]>=0){ // edge alone4179 if (i<nbe){4180 long i0=edge4->i(i);4181 ordre[i0] = vertices+i0;4182 long i1=edge4->j(i);4183 ordre[i1] = vertices+i1;4184 }4185 else {4186 k=k+1;4187 if (k<10) {4188 //print only 10 edges4189 printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));4190 }4191 else if (k==10){4192 printf("Other lost boundary edges not shown...\n");4193 }4194 }4195 }4196 }4197 if(k) {4198 ISSMERROR("%i boundary edges (from the geometry) are not defined as mesh edges",k);4199 }4200 4201 /* mesh generation with boundary points*/4202 long nbvb=0;4203 for (i=0;i<nbv;i++){4204 vertices[i].t=0;4205 vertices[i].vint=0;4206 if (ordre[i]) ordre[nbvb++]=ordre[i];4207 }4208 4209 Triangle* savetriangles=triangles;4210 long savenbt=nbt;4211 long savenbtx=nbtx;4212 SubDomain* savesubdomains=subdomains;4213 subdomains=0;4214 4215 long Nbtriafillhole=2*nbvb;4216 Triangle* triafillhole=new Triangle[Nbtriafillhole];4217 triangles = triafillhole;4218 4219 nbt=2;4220 nbtx= Nbtriafillhole;4221 4222 //Find a vertex that is not aligned with vertices 0 and 14223 for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;)4224 if (++i>=nbvb) {4225 ISSMERROR("ReconstructExistingMesh: All the vertices are aligned");4226 }4227 //Move this vertex (i) to the 2d position in ordre4228 Exchange(ordre[2], ordre[i]);4229 4230 /*Reconstruct mesh beginning with 2 triangles*/4231 MeshVertex * v0=ordre[0], *v1=ordre[1];4232 4233 triangles[0](0) = NULL; // Infinite vertex4234 triangles[0](1) = v0;4235 triangles[0](2) = v1;4236 4237 triangles[1](0) = NULL;// Infinite vertex4238 triangles[1](2) = v0;4239 triangles[1](1) = v1;4240 const int e0 = OppositeEdge[0];4241 const int e1 = NextEdge[e0];4242 const int e2 = PreviousEdge[e0];4243 triangles[0].SetAdj2(e0, &triangles[1] ,e0);4244 triangles[0].SetAdj2(e1, &triangles[1] ,e2);4245 triangles[0].SetAdj2(e2, &triangles[1] ,e1);4246 4247 triangles[0].det = -1; // boundary triangles4248 triangles[1].det = -1; // boundary triangles4249 4250 triangles[0].SetTriangleContainingTheVertex();4251 triangles[1].SetTriangleContainingTheVertex();4252 4253 triangles[0].link=&triangles[1];4254 triangles[1].link=&triangles[0];4255 4256 if (!quadtree) delete quadtree; //ReInitialise;4257 quadtree = new QuadTree(this,0);4258 quadtree->Add(*v0);4259 quadtree->Add(*v1);4260 4261 // vertices are added one by one4262 long NbSwap=0;4263 for (int icount=2; icount<nbvb; icount++) {4264 MeshVertex *vi = ordre[icount];4265 Icoor2 dete[3];4266 Triangle *tcvi = FindTriangleContaining(vi->i,dete);4267 quadtree->Add(*vi);4268 AddVertex(*vi,tcvi,dete);4269 NbSwap += vi->Optim(1,1);4270 }4271 4272 //enforce the boundary4273 TriangleAdjacent ta(0,0);4274 long nbloss = 0,knbe=0;4275 for ( i = 0; i < nbe; i++){4276 if (st[i] >=0){ //edge alone => on border4277 MeshVertex &a=edges[i][0], &b=edges[i][1];4278 if (a.t && b.t){4279 knbe++;4280 if (ForceEdge(a,b,ta)<0) nbloss++;4281 }4282 }4283 }4284 if(nbloss) {4285 ISSMERROR("we lost %i existing edges other %i",nbloss,knbe);4286 }4287 4288 FindSubDomain(1);4289 // remove all the hole4290 // remove all the good sub domain4291 long krm =0;4292 for (i=0;i<nbt;i++){4293 if (triangles[i].link){ // remove triangles4294 krm++;4295 for (int j=0;j<3;j++){4296 TriangleAdjacent ta = triangles[i].Adj(j);4297 Triangle &tta = *(Triangle*)ta;4298 //if edge between remove and not remove4299 if(! tta.link){4300 // change the link of ta;4301 int ja = ta;4302 MeshVertex *v0= ta.EdgeVertex(0);4303 MeshVertex *v1= ta.EdgeVertex(1);4304 long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);4305 4306 ISSMASSERT(st[k]>=0);4307 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));4308 ta.SetLock();4309 st[k]=-2-st[k];4310 }4311 }4312 }4313 }4314 long NbTfillHoll =0;4315 for (i=0;i<nbt;i++){4316 if (triangles[i].link) {4317 triangles[i]=Triangle((MeshVertex *) NULL,(MeshVertex *) NULL,(MeshVertex *) NULL);4318 triangles[i].color=-1;4319 }4320 else{4321 triangles[i].color= savenbt+ NbTfillHoll++;4322 }4323 }4324 ISSMASSERT(savenbt+NbTfillHoll<=savenbtx);4325 4326 // copy of the outside triangles in saveMesh4327 for (i=0;i<nbt;i++){4328 if(triangles[i].color>=0) {4329 savetriangles[savenbt]=triangles[i];4330 savetriangles[savenbt].link=0;4331 savenbt++;4332 }4333 }4334 // gestion of the adj4335 k =0;4336 Triangle * tmax = triangles + nbt;4337 for (i=0;i<savenbt;i++) {4338 Triangle & ti = savetriangles[i];4339 for (int j=0;j<3;j++){4340 Triangle * ta = ti.TriangleAdj(j);4341 int aa = ti.NuEdgeTriangleAdj(j);4342 int lck = ti.Locked(j);4343 if (!ta) k++; // bug4344 else if ( ta >= triangles && ta < tmax){4345 ta= savetriangles + ta->color;4346 ti.SetAdj2(j,ta,aa);4347 if(lck) ti.SetLocked(j);4348 }4349 }4350 }4351 4352 // restore triangles;4353 nbt=savenbt;4354 nbtx=savenbtx;4355 delete [] triangles;4356 delete [] subdomains;4357 triangles = savetriangles;4358 subdomains = savesubdomains;4359 if (k) {4360 ISSMERROR("number of triangles edges alone = %i",k);4361 }4362 FindSubDomain();4363 4364 delete edge4;4365 delete [] st;4366 for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);4367 4368 SetVertexFieldOn();4369 4370 /*Check requirements consistency*/4371 for (i=0;i<nbe;i++){4372 /*If the current mesh edge is on Geometry*/4373 if(edges[i].onGeometry){4374 for(int j=0;j<2;j++){4375 /*Go through the edges adjacent to current edge (if on the same curve)*/4376 if (!edges[i].adj[j]){4377 /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/4378 /*Check that the 2 vertices are on geometry AND required*/4379 if(!edges[i][j].onGeometry->IsRequiredVertex()){4380 printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1);4381 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1);4382 if (edges[i][j].onGeometry->OnGeomVertex())4383 printf("the vertex number %i of this edge is a geometric MeshVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->gv)+1);4384 else if (edges[i][j].onGeometry->OnGeomEdge())4385 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->ge)+1);4386 else4387 printf("Its pointer is %p\n",edges[i][j].onGeometry);4388 4389 printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");4390 ISSMERROR("See above (might be cryptic...)");4391 }4392 }4393 }4394 }4395 }4396 }4397 /*}}}1*/4398 /*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/4399 void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){4400 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/4401 4402 long int verbose=0;4403 long *renu= new long[nbt];4404 register Triangle *t0,*t,*te=triangles+nbt;4405 register long k=0,it,i,j;4406 4407 for ( it=0;it<nbt;it++)4408 renu[it]=-1; // outside triangle4409 for ( i=0;i<NbSubDomains;i++)4410 {4411 t=t0=subdomains[i].head;4412 if (!t0){ // not empty sub domain4413 ISSMERROR("!t0");4414 }4415 do {4416 long kt = Number(t);4417 if (kt<0 || kt >= nbt ){4418 ISSMERROR("kt<0 || kt >= nbt");4419 }4420 if (renu[kt]!=-1){4421 ISSMERROR("renu[kt]!=-1");4422 }4423 renu[kt]=k++;4424 }4425 while (t0 != (t=t->link));4426 }4427 // take is same numbering if possible4428 if(justcompress)4429 for ( k=0,it=0;it<nbt;it++)4430 if(renu[it] >=0 )4431 renu[it]=k++;4432 4433 // put the outside triangles at the end4434 for ( it=0;it<nbt;it++){4435 if (renu[it]==-1) renu[it]=k++;4436 }4437 if (k != nbt){4438 ISSMERROR("k != nbt");4439 }4440 // do the change on all the pointeur4441 for ( it=0;it<nbt;it++)4442 triangles[it].ReNumbering(triangles,te,renu);4443 4444 for ( i=0;i<NbSubDomains;i++)4445 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];4446 4447 // move the Triangles without a copy of the array4448 // be carefull not trivial code4449 for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu4450 if (renu[it] >= 0) // a new sub cycle4451 {4452 i=it;4453 Triangle ti=triangles[i],tj;4454 while ( (j=renu[i]) >= 0)4455 { // i is old, and j is new4456 renu[i] = -1; // mark4457 tj = triangles[j]; // save new4458 triangles[j]= ti; // new <- old4459 i=j; // next4460 ti = tj;4461 }4462 }4463 delete [] renu;4464 nt = nbt - NbOutT;4465 4466 }4467 /*}}}1*/4468 /*FUNCTION Mesh::ReNumberingVertex{{{1*/4469 void Mesh::ReNumberingVertex(long * renu) {4470 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/4471 4472 // warning be carfull because pointer4473 // from on mesh to over mesh4474 // -- so do ReNumbering at the beginning4475 MeshVertex * ve = vertices+nbv;4476 long it,ie,i;4477 4478 printf("renumbering triangles\n");4479 for ( it=0;it<nbt;it++)4480 triangles[it].ReNumbering(vertices,ve,renu);4481 4482 printf("renumbering edges\n");4483 for ( ie=0;ie<nbe;ie++)4484 edges[ie].ReNumbering(vertices,ve,renu);4485 4486 printf("renumbering vertices on geom\n");4487 for (i=0;i< NbVerticesOnGeomVertex;i++)4488 {4489 MeshVertex *v = VerticesOnGeomVertex[i].mv;4490 if (v>=vertices && v < ve)4491 VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];4492 }4493 4494 printf("renumbering vertices on edge\n");4495 for (i=0;i< NbVerticesOnGeomEdge;i++)4496 {4497 MeshVertex *v =VerticesOnGeomEdge[i].mv;4498 if (v>=vertices && v < ve)4499 VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)];4500 }4501 4502 printf("renumbering vertices on Bth vertex\n");4503 for (i=0;i< NbVertexOnBThVertex;i++)4504 {4505 MeshVertex *v=VertexOnBThVertex[i].v;4506 if (v>=vertices && v < ve)4507 VertexOnBThVertex[i].v=vertices+renu[Number(v)];4508 }4509 4510 for (i=0;i< NbVertexOnBThEdge;i++)4511 {4512 MeshVertex *v=VertexOnBThEdge[i].v;4513 if (v>=vertices && v < ve)4514 VertexOnBThEdge[i].v=vertices+renu[Number(v)];4515 }4516 4517 // move the Vertices without a copy of the array4518 // be carefull not trivial code4519 long j;4520 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu4521 if (renu[it] >= 0) // a new sub cycle4522 {4523 i=it;4524 MeshVertex ti=vertices[i],tj;4525 while ( (j=renu[i]) >= 0){4526 // i is old, and j is new4527 renu[i] = -1-renu[i]; // mark4528 tj = vertices[j]; // save new4529 vertices[j]= ti; // new <- old4530 i=j; // next4531 ti = tj;4532 }4533 }4534 if (quadtree){4535 delete quadtree;4536 quadtree = new QuadTree(this);4537 }4538 4539 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;4540 }4541 /*}}}1*/4542 /*FUNCTION Mesh::SetIntCoor{{{1*/4543 void Mesh::SetIntCoor(const char * strfrom) {4544 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/4545 4546 /*Set integer coordinate for existing vertices*/4547 4548 //Get extrema coordinates of the existing vertices4549 pmin = vertices[0].r;4550 pmax = vertices[0].r;4551 long i;4552 for (i=0;i<nbv;i++) {4553 pmin.x = Min(pmin.x,vertices[i].r.x);4554 pmin.y = Min(pmin.y,vertices[i].r.y);4555 pmax.x = Max(pmax.x,vertices[i].r.x);4556 pmax.y = Max(pmax.y,vertices[i].r.y);4557 }4558 R2 DD = (pmax-pmin)*0.05;4559 pmin = pmin-DD;4560 pmax = pmax+DD;4561 4562 //Compute coefIcoor4563 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));4564 if (coefIcoor<=0){4565 ISSMERROR("coefIcoor should be positive, a problem in the geometry is likely");4566 }4567 4568 // generation of integer coord4569 for (i=0;i<nbv;i++) {4570 vertices[i].i = toI2(vertices[i].r);4571 }4572 4573 // computation of the det4574 int number_of_errors=0;4575 for (i=0;i<nbt;i++) {4576 MeshVertex & v0 = triangles[i][0];4577 MeshVertex & v1 = triangles[i][1];4578 MeshVertex & v2 = triangles[i][2];4579 4580 //If this is not a boundary triangle4581 if ( &v0 && &v1 && &v2 ){4582 4583 /*Compute determinant*/4584 triangles[i].det= det(v0,v1,v2);4585 4586 /*Check that determinant is positive*/4587 if (triangles[i].det <=0){4588 4589 /*increase number_of_errors and print error only for the first 20 triangles*/4590 number_of_errors++;4591 if (number_of_errors<20){4592 printf("Area of Triangle %i < 0 (det=%i)\n",i+1,triangles[i].det);4593 }4594 }4595 }4596 4597 //else, set as -14598 else triangles[i].det=-1;4599 }4600 4601 if (number_of_errors) ISSMERROR("Fatal error: some triangles have negative areas, see above");4602 }4603 /*}}}1*/4604 /*FUNCTION Mesh::ShowRegulaty{{{1*/4605 void Mesh::ShowRegulaty() const {4606 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/4607 4608 const double sqrt32=sqrt(3.)*0.5;4609 const double aireKh=sqrt32*0.5;4610 D2 Beq(1,0),Heq(0.5,sqrt32);4611 D2xD2 Br(D2xD2(Beq,Heq).t());4612 D2xD2 B1r(Br.inv());4613 double gammamn=1e100,hmin=1e100;4614 double gammamx=0,hmax=0;4615 double beta=1e100;4616 double beta0=0;4617 double alpha2=0;4618 double area=0,Marea=0;4619 // double cf= double(coefIcoor);4620 // double cf2= 6.*cf*cf;4621 int nt=0;4622 for (int it=0;it<nbt;it++)4623 if ( triangles[it].link)4624 {4625 nt++;4626 Triangle &K=triangles[it];4627 double area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;4628 area+= area3;4629 D2xD2 B_Kt(K[0],K[1],K[2]);4630 D2xD2 B_K(B_Kt.t());4631 D2xD2 B1K = Br*B_K.inv();4632 D2xD2 BK = B_K*B1r;4633 D2xD2 B1B1 = B1K.t()*B1K;4634 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);4635 MatVVP2x2 VMK(MK);4636 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));4637 double betaK=0;4638 4639 for (int j=0;j<3;j++)4640 {4641 double he= Norme2(R2(K[j],K[(j+1)%3]));4642 hmin=Min(hmin,he);4643 hmax=Max(hmax,he);4644 MeshVertex & v=K[j];4645 D2xD2 M((Metric)v);4646 betaK += sqrt(M.det());4647 4648 D2xD2 BMB = BK.t()*M*BK;4649 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);4650 MatVVP2x2 VM1(M1);4651 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);4652 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);4653 }4654 betaK *= area3;// 1/2 (somme sqrt(det))* area(K)4655 Marea+= betaK;4656 beta=min(beta,betaK);4657 beta0=max(beta0,betaK);4658 }4659 area*=3;4660 gammamn=sqrt(gammamn);4661 gammamx=sqrt(gammamx);4662 printf(" Adaptmesh info:\n");4663 printf(" number of triangles = %i\n",nt);4664 printf(" hmin = %g, hmax=%g\n",hmin,hmax);4665 printf(" area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt));4666 printf(" infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx);4667 printf(" anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5));4668 }4669 /*}}}1*/4670 /*FUNCTION Mesh::ShowHistogram{{{1*/4671 void Mesh::ShowHistogram() const {4672 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/4673 4674 const long kmax=10;4675 const double llmin = 0.5,llmax=2;4676 const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);4677 long histo[kmax+1];4678 long i,it,k, nbedges =0;4679 for (i=0;i<=kmax;i++) histo[i]=0;4680 for (it=0;it<nbt;it++)4681 if ( triangles[it].link)4682 {4683 4684 for (int j=0;j<3;j++)4685 {4686 Triangle *ta = triangles[it].TriangleAdj(j);4687 if ( !ta || !ta->link || Number(ta) >= it)4688 {4689 MeshVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];4690 MeshVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];4691 if ( !& vP || !&vQ) continue;4692 R2 PQ = vQ.r - vP.r;4693 double l = log(LengthInterpole(vP,vQ,PQ));4694 nbedges++;4695 k = (int) ((l - lmin)*delta);4696 k = Min(Max(k,0L),kmax);4697 histo[k]++;4698 }4699 }4700 }4701 printf(" --- Histogram of the unit mesh, nb of edges = %i\n",nbedges);4702 printf(" length of edge in | %% of edge | Nb of edges \n");4703 printf(" --------------------+-------------+-------------\n");4704 for (i=0;i<=kmax;i++){4705 if (i==0) printf(" %10i",0);4706 else printf(" %10g",exp(lmin+i/delta));4707 if (i==kmax) printf(" +inf ");4708 else printf(" %10g",exp(lmin+(i+1)/delta));4709 printf("| %10g |\n",((long) ((10000.0 * histo[i])/ nbedges))/100.0);4710 printf(" %i\n",histo[i]);4711 }4712 printf(" --------------------+-------------+-------------\n");4713 }4714 /*}}}1*/4715 /*FUNCTION Mesh::SmoothingVertex{{{1*/4716 void Mesh::SmoothingVertex(int nbiter,double omega ) {4717 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/4718 4719 long int verbose=0;4720 // if quatree exist remove it end reconstruct4721 if (quadtree) delete quadtree;4722 quadtree=0;4723 ReMakeTriangleContainingTheVertex();4724 Triangle vide; // a triangle to mark the boundary vertex4725 Triangle ** tstart= new Triangle* [nbv];4726 long i,j,k;4727 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide4728 if ( this == & BTh)4729 for ( i=0;i<nbv;i++)4730 tstart[i]=vertices[i].t;4731 else4732 for ( i=0;i<nbv;i++)4733 tstart[i]=0;4734 for ( j=0;j<NbVerticesOnGeomVertex;j++ )4735 tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;4736 for ( k=0;k<NbVerticesOnGeomEdge;k++ )4737 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;4738 if(verbose>2) printf(" SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega);4739 for (k=0;k<nbiter;k++)4740 {4741 long i,NbSwap =0;4742 double delta =0;4743 for ( i=0;i<nbv;i++)4744 if (tstart[i] != &vide) // not a boundary vertex4745 delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));4746 if (!NbOfQuad)4747 for ( i=0;i<nbv;i++)4748 if (tstart[i] != &vide) // not a boundary vertex4749 NbSwap += vertices[i].Optim(1);4750 if (verbose>3) printf(" move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap);4751 }4752 4753 delete [] tstart;4754 if (quadtree) quadtree= new QuadTree(this);4755 }4756 /*}}}1*/4757 /*FUNCTION Mesh::SmoothMetric{{{1*/4758 void Mesh::SmoothMetric(double raisonmax) {4759 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/4760 4761 long int verbose=0;4762 4763 if(raisonmax<1.1) return;4764 if(verbose > 1) printf(" Mesh::SmoothMetric raisonmax = %g\n",raisonmax);4765 ReMakeTriangleContainingTheVertex();4766 long i,j,kch,kk,ip;4767 long *first_np_or_next_t0 = new long[nbv];4768 long *first_np_or_next_t1 = new long[nbv];4769 long Head0 =0,Head1=-1;4770 double logseuil= log(raisonmax);4771 4772 for(i=0;i<nbv-1;i++)4773 first_np_or_next_t0[i]=i+1;4774 first_np_or_next_t0[nbv-1]=-1;// end;4775 for(i=0;i<nbv;i++)4776 first_np_or_next_t1[i]=-1;4777 kk=0;4778 while (Head0>=0&& kk++<100) {4779 kch=0;4780 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {4781 // pour tous les triangles autour du sommet s4782 register Triangle* t= vertices[i].t;4783 if (!t){4784 ISSMERROR("!t");4785 }4786 MeshVertex & vi = vertices[i];4787 TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]);4788 MeshVertex *pvj0 = ta.EdgeVertex(0);4789 while (1) {4790 ta=Previous(Adj(ta));4791 if (vertices+i != ta.EdgeVertex(1)){4792 ISSMERROR("vertices+i != ta.EdgeVertex(1)");4793 }4794 MeshVertex & vj = *(ta.EdgeVertex(0));4795 if ( &vj ) {4796 j= &vj-vertices;4797 if (j<0 || j >= nbv){4798 ISSMERROR("j<0 || j >= nbv");4799 }4800 R2 Aij = (R2) vj - (R2) vi;4801 double ll = Norme2(Aij);4802 if (0) {4803 double hi = ll/vi.m(Aij);4804 double hj = ll/vj.m(Aij);4805 if(hi < hj)4806 {4807 double dh=(hj-hi)/ll;4808 if (dh>logseuil) {4809 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));4810 if(first_np_or_next_t1[j]<0)4811 kch++,first_np_or_next_t1[j]=Head1,Head1=j;4812 }4813 }4814 }4815 else4816 {4817 double li = vi.m(Aij);4818 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )4819 if(first_np_or_next_t1[j]<0) // if the metrix change4820 kch++,first_np_or_next_t1[j]=Head1,Head1=j;4821 }4822 }4823 if ( &vj == pvj0 ) break;4824 }4825 }4826 Head0 = Head1;4827 Head1 = -1;4828 Exchange(first_np_or_next_t0,first_np_or_next_t1);4829 }4830 if(verbose>2) printf(" number of iterations = %i\n",kch);4831 delete [] first_np_or_next_t0;4832 delete [] first_np_or_next_t1;4833 }4834 /*}}}1*/4835 /*FUNCTION Mesh::SplitElement{{{1*/4836 int Mesh::SplitElement(int choice){4837 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/4838 4839 long int verbose=0;4840 4841 Direction NoDirOfSearch;4842 const int withBackground = &BTh != this && &BTh;4843 4844 ReNumberingTheTriangleBySubDomain();4845 //int nswap =0;4846 const long nfortria( choice ? 4 : 6);4847 if(withBackground)4848 {4849 BTh.SetVertexFieldOn();4850 SetVertexFieldOnBTh();4851 }4852 else4853 BTh.SetVertexFieldOn();4854 4855 long newnbt=0,newnbv=0;4856 long * kedge = 0;4857 long newNbOfQuad=NbOfQuad;4858 long * ksplit= 0, * ksplitarray=0;4859 long kkk=0;4860 int ret =0;4861 if (nbvx<nbv+nbe) return 1;//4862 // 1) create the new points by spliting the internal edges4863 // set the4864 long nbvold = nbv;4865 long nbtold = nbt;4866 long NbOutTold = NbOutT;4867 long NbEdgeOnGeom=0;4868 long i;4869 4870 nbt = nbt - NbOutT; // remove all the the ouside triangles4871 long nbtsave = nbt;4872 Triangle * lastT = triangles + nbt;4873 for (i=0;i<nbe;i++)4874 if(edges[i].onGeometry) NbEdgeOnGeom++;4875 long newnbe=nbe+nbe;4876 // long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;4877 long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;4878 // long newNbVertexOnBThVertex=NbVertexOnBThVertex;4879 long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;4880 4881 // do allocation for pointeur to the geometry and background4882 VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];4883 VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ? new VertexOnEdge[newNbVertexOnBThEdge]:0;4884 if (NbVerticesOnGeomEdge)4885 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);4886 if (NbVertexOnBThEdge)4887 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);4888 Edge *newedges = new Edge [newnbe];4889 // memcpy(newedges,edges,sizeof(Edge)*nbe);4890 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);4891 long k=nbv;4892 long kk=0;4893 long kvb = NbVertexOnBThEdge;4894 long kvg = NbVerticesOnGeomEdge;4895 long ie =0;4896 Edge ** edgesGtoB=0;4897 if (withBackground)4898 edgesGtoB= BTh.MakeGeometricalEdgeToEdge();4899 long ferr=0;4900 for (i=0;i<nbe;i++)4901 newedges[ie].onGeometry=0;4902 4903 for (i=0;i<nbe;i++)4904 {4905 GeometricalEdge *ong = edges[i].onGeometry;4906 4907 newedges[ie]=edges[i];4908 newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;4909 newedges[ie].adj[1]=newedges + ie +1;4910 R2 A = edges[i][0],B = edges[i][1];4911 4912 4913 kk += (i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));4914 if (ong) // a geometrical edges4915 {4916 if (withBackground){4917 // walk on back ground mesh4918 // newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge);4919 // a faire -- difficile4920 // the first PB is to now a background edge between the 2 vertices4921 if (!edgesGtoB){4922 ISSMERROR("!edgesGtoB");4923 }4924 ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],4925 edges[i][0],edges[i][1],0.5,vertices[k],4926 newVertexOnBThEdge[kvb],4927 newVerticesOnGeomEdge[kvg++]);4928 vertices[k].ReferenceNumber= edges[i].ref;4929 vertices[k].DirOfSearch = NoDirOfSearch;4930 ;4931 // get the Info on background mesh4932 double s = newVertexOnBThEdge[kvb];4933 MeshVertex & bv0 = newVertexOnBThEdge[kvb][0];4934 MeshVertex & bv1 = newVertexOnBThEdge[kvb][1];4935 // compute the metrix of the new points4936 vertices[k].m = Metric(1-s,bv0,s,bv1);4937 kvb++;4938 }4939 else4940 {4941 ong=Gh.ProjectOnCurve(edges[i],4942 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);4943 // vertices[k].i = toI2( vertices[k].r);4944 vertices[k].ReferenceNumber = edges[i].ref;4945 vertices[k].DirOfSearch = NoDirOfSearch;4946 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]);4947 }4948 }4949 else // straigth line edge ---4950 {4951 vertices[k].r = ((R2) edges[i][0] + (R2) edges[i][1] )*0.5;4952 vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]);4953 vertices[k].onGeometry = 0;4954 }4955 //vertices[k].i = toI2( vertices[k].r);4956 R2 AB = vertices[k].r;4957 R2 AA = (A+AB)*0.5;4958 R2 BB = (AB+B)*0.5;4959 vertices[k].ReferenceNumber = edges[i].ref;4960 vertices[k].DirOfSearch = NoDirOfSearch;4961 4962 newedges[ie].onGeometry = Gh.Containing(AA,ong);4963 newedges[ie++].v[1]=vertices+k;4964 4965 newedges[ie]=edges[i];4966 newedges[ie].adj[0]=newedges + ie -1;4967 newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;4968 newedges[ie].onGeometry = Gh.Containing(BB,ong);4969 newedges[ie++].v[0]=vertices+k;4970 k++;4971 }4972 if (edgesGtoB) delete [] edgesGtoB;4973 edgesGtoB=0;4974 4975 newnbv=k;4976 newNbVerticesOnGeomEdge=kvg;4977 if (newnbv> nbvx) goto Error;// bug4978 4979 nbv = k;4980 4981 4982 kedge = new long[3*nbt+1];4983 ksplitarray = new long[nbt+1];4984 ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]4985 4986 for (i=0;i<3*nbt;i++)4987 kedge[i]=-1;4988 4989 //4990 4991 for (i=0;i<nbt;i++) {4992 Triangle & t = triangles[i];4993 if (!t.link){4994 ISSMERROR("!t.link");4995 }4996 for(int j=0;j<3;j++)4997 {4998 const TriangleAdjacent ta = t.Adj(j);4999 const Triangle & tt = ta;5000 if (&tt >= lastT)5001 t.SetAdj2(j,0,0);// unset adj5002 const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];5003 const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];5004 long ke =edge4->SortAndFind(Number(v0),Number(v1));5005 if (ke>0)5006 {5007 long ii = Number(tt);5008 int jj = ta;5009 long ks = ke + nbvold;5010 kedge[3*i+j] = ks;5011 if (ii<nbt) // good triangle5012 kedge[3*ii+jj] = ks;5013 MeshVertex &A=vertices[ks];5014 double aa,bb,cc,dd;5015 if ((dd=Area2(v0.r,v1.r,A.r)) >=0){5016 // warning PB roundoff error5017 if (t.link && ( (aa=Area2( A.r , t[1].r , t[2].r )) < 0.05018 || (bb=Area2( t[0].r , A.r , t[2].r )) < 0.05019 || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)){5020 printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd);5021 ISSMERROR("Number of triangles with P2 interpolation Problem");5022 }5023 }5024 else {5025 if (tt.link && ( (aa=Area2( A.r , tt[1].r , tt[2].r )) < 05026 || (bb=Area2( tt[0].r , A.r , tt[2].r )) < 05027 || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)){5028 printf("%i not in triangle %i In= %i %g %g %g %g\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd);5029 ISSMERROR("Number of triangles with P2 interpolation Problem");5030 }5031 }5032 }5033 }5034 }5035 5036 for (i=0;i<nbt;i++){5037 ksplit[i]=1; // no split by default5038 const Triangle & t = triangles[ i];5039 int nbsplitedge =0;5040 int nbinvisible =0;5041 int invisibleedge=0;5042 int kkk[3];5043 for (int j=0;j<3;j++)5044 {5045 if (t.Hidden(j)) invisibleedge=j,nbinvisible++;5046 5047 const TriangleAdjacent ta = t.Adj(j);5048 const Triangle & tt = ta;5049 5050 5051 const MeshVertex & v0 = t[VerticesOfTriangularEdge[j][0]];5052 const MeshVertex & v1 = t[VerticesOfTriangularEdge[j][1]];5053 if ( kedge[3*i+j] < 0)5054 {5055 long ke =edge4->SortAndFind(Number(v0),Number(v1));5056 if (ke<0) // new5057 {5058 if (&tt) // internal triangles all the boundary5059 { // new internal edges5060 long ii = Number(tt);5061 int jj = ta;5062 5063 kedge[3*i+j]=k;// save the vertex number5064 kedge[3*ii+jj]=k;5065 if (k<nbvx)5066 {5067 vertices[k].r = ((R2) v0+(R2) v1 )/2;5068 //vertices[k].i = toI2( vertices[k].r);5069 vertices[k].ReferenceNumber=0;5070 vertices[k].DirOfSearch =NoDirOfSearch;5071 vertices[k].m = Metric(0.5,v0,0.5,v1);5072 }5073 k++;5074 kkk[nbsplitedge++]=j;5075 } // tt5076 else5077 ISSMERROR("Bug...");5078 } // ke<05079 else5080 { // ke >=05081 kedge[3*i+j]=nbvold+ke;5082 kkk[nbsplitedge++]=j;// previously splited5083 }5084 }5085 else5086 kkk[nbsplitedge++]=j;// previously splited5087 5088 }5089 if (nbinvisible>=2){5090 ISSMERROR("nbinvisible>=2");5091 }5092 switch (nbsplitedge) {5093 case 0: ksplit[i]=10; newnbt++; break; // nosplit5094 case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 25095 case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 35096 case 3:5097 if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;5098 else ksplit[i]=10*nfortria,newnbt+=nfortria;5099 break;5100 }5101 if (ksplit[i]<40){5102 ISSMERROR("ksplit[i]<40");5103 }5104 }5105 // now do the element split5106 newNbOfQuad = 4*NbOfQuad;5107 nbv = k;5108 kkk = nbt;5109 ksplit[-1] = nbt;5110 // look on old true triangles5111 5112 for (i=0;i<nbtsave;i++){5113 int nbmkadj=0;5114 long mkadj [100];5115 mkadj[0]=i;5116 long kk=ksplit[i]/10;5117 int ke=(int) (ksplit[i]%10);5118 if (kk>=7 || kk<=0){5119 ISSMERROR("kk>=7 || kk<=0");5120 }5121 5122 // def the numbering k (edge) i vertex5123 int k0 = ke;5124 int k1 = NextEdge[k0];5125 int k2 = PreviousEdge[k0];5126 int i0 = OppositeVertex[k0];5127 int i1 = OppositeVertex[k1];5128 int i2 = OppositeVertex[k2];5129 5130 Triangle &t0=triangles[i];5131 MeshVertex * v0=t0(i0);5132 MeshVertex * v1=t0(i1);5133 MeshVertex * v2=t0(i2);5134 5135 if (nbmkadj>=10){5136 ISSMERROR("nbmkadj>=10");5137 }5138 // --------------------------5139 TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));5140 // save the flag Hidden5141 int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};5142 // un set all adj -- save Hidden flag --5143 t0.SetAdj2(0,0,hid[0]);5144 t0.SetAdj2(1,0,hid[1]);5145 t0.SetAdj2(2,0,hid[2]);5146 // -- remake5147 switch (kk) {5148 case 1: break;// nothing5149 case 2: //5150 {5151 Triangle &t1=triangles[kkk++];5152 t1=t0;5153 if (kedge[3*i+i0]<0){5154 ISSMERROR("kedge[3*i+i0]<0");5155 }5156 MeshVertex * v3 = vertices + kedge[3*i+k0];5157 5158 t0(i2) = v3;5159 t1(i1) = v3;5160 t0.SetAllFlag(k2,0);5161 t1.SetAllFlag(k1,0);5162 }5163 break;5164 case 3: //5165 {5166 Triangle &t1=triangles[kkk++];5167 Triangle &t2=triangles[kkk++];5168 t2=t1=t0;5169 if (kedge[3*i+k1]<0){5170 ISSMERROR("kedge[3*i+k1]<0");5171 }5172 if (kedge[3*i+k2]<0){5173 ISSMERROR("kedge[3*i+k2]<0");5174 }5175 5176 MeshVertex * v01 = vertices + kedge[3*i+k2];5177 MeshVertex * v02 = vertices + kedge[3*i+k1];5178 t0(i1) = v01;5179 t0(i2) = v02;5180 t1(i2) = v02;5181 t1(i0) = v01;5182 t2(i0) = v02;5183 t0.SetAllFlag(k0,0);5184 t1.SetAllFlag(k1,0);5185 t1.SetAllFlag(k0,0);5186 t2.SetAllFlag(k2,0);5187 }5188 break;5189 case 4: //5190 case 6: // split in 45191 {5192 Triangle &t1=triangles[kkk++];5193 Triangle &t2=triangles[kkk++];5194 Triangle &t3=triangles[kkk++];5195 t3=t2=t1=t0;5196 if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){5197 ISSMERROR("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");5198 }5199 MeshVertex * v12 = vertices + kedge[3*i+k0];5200 MeshVertex * v02 = vertices + kedge[3*i+k1];5201 MeshVertex * v01 = vertices + kedge[3*i+k2];5202 t0(i1) = v01;5203 t0(i2) = v02;5204 t0.SetAllFlag(k0,hid[k0]);5205 5206 t1(i0) = v01;5207 t1(i2) = v12;5208 t0.SetAllFlag(k1,hid[k1]);5209 5210 t2(i0) = v02;5211 t2(i1) = v12;5212 t2.SetAllFlag(k2,hid[k2]);5213 5214 t3(i0) = v12;5215 t3(i1) = v02;5216 t3(i2) = v01;5217 5218 t3.SetAllFlag(0,hid[0]);5219 t3.SetAllFlag(1,hid[1]);5220 t3.SetAllFlag(2,hid[2]);5221 5222 if ( kk == 6)5223 {5224 5225 Triangle &t4=triangles[kkk++];5226 Triangle &t5=triangles[kkk++];5227 5228 t4 = t3;5229 t5 = t3;5230 5231 t0.SetHidden(k0);5232 t1.SetHidden(k1);5233 t2.SetHidden(k2);5234 t3.SetHidden(0);5235 t4.SetHidden(1);5236 t5.SetHidden(2);5237 5238 if (nbv < nbvx )5239 {5240 vertices[nbv].r = ((R2) *v01 + (R2) *v12 + (R2) *v02 ) / 3.0;5241 vertices[nbv].ReferenceNumber =0;5242 vertices[nbv].DirOfSearch =NoDirOfSearch;5243 //vertices[nbv].i = toI2(vertices[nbv].r);5244 double a3[]={1./3.,1./3.,1./3.};5245 vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);5246 MeshVertex * vc = vertices +nbv++;5247 t3(i0) = vc;5248 t4(i1) = vc;5249 t5(i2) = vc;5250 5251 }5252 else5253 goto Error;5254 }5255 5256 }5257 break;5258 }5259 5260 // save all the new triangles5261 mkadj[nbmkadj++]=i;5262 long jj;5263 if (t0.link)5264 for (jj=nbt;jj<kkk;jj++)5265 {5266 triangles[jj].link=t0.link;5267 t0.link= triangles+jj;5268 mkadj[nbmkadj++]=jj;5269 }5270 if (nbmkadj>13){// 13 = 6 + 4 +5271 ISSMERROR("nbmkadj>13");5272 }5273 5274 if (kk==6) newNbOfQuad+=3;5275 for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;5276 ksplit[i]= nbt; // save last adresse of the new triangles5277 kkk = nbt;5278 }5279 5280 for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;5281 5282 if(withBackground)5283 for (i=0;i<BTh.nbv;i++)5284 BTh.vertices[i].m = BTh.vertices[i].m*2.;5285 5286 5287 ret = 2;5288 if (nbt>= nbtx) goto Error; // bug5289 if (nbv>= nbvx) goto Error; // bug5290 // generation of the new triangles5291 5292 SetIntCoor("In SplitElement");5293 5294 ReMakeTriangleContainingTheVertex();5295 if(withBackground)5296 BTh.ReMakeTriangleContainingTheVertex();5297 5298 delete [] edges;5299 edges = newedges;5300 nbe = newnbe;5301 NbOfQuad = newNbOfQuad;5302 5303 for (i=0;i<NbSubDomains;i++)5304 {5305 long k = subdomains[i].edge- edges;5306 subdomains[i].edge = edges+2*k; // spilt all edge in 25307 }5308 5309 if (ksplitarray) delete [] ksplitarray;5310 if (kedge) delete [] kedge;5311 if (edge4) delete edge4;5312 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;5313 VerticesOnGeomEdge= newVerticesOnGeomEdge;5314 if(VertexOnBThEdge) delete [] VertexOnBThEdge;5315 VertexOnBThEdge = newVertexOnBThEdge;5316 NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;5317 NbVertexOnBThEdge=newNbVertexOnBThEdge;5318 // ReMakeTriangleContainingTheVertex();5319 5320 ReconstructExistingMesh();5321 5322 if (verbose>2){5323 printf(" number of quadrilaterals = %i\n",NbOfQuad);5324 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2);5325 printf(" number of outside triangles = %i\n",NbOutT);5326 }5327 5328 return 0; //ok5329 5330 Error:5331 nbv = nbvold;5332 nbt = nbtold;5333 NbOutT = NbOutTold;5334 // cleaning memory ---5335 delete newedges;5336 if (ksplitarray) delete [] ksplitarray;5337 if (kedge) delete [] kedge;5338 if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;5339 if (edge4) delete edge4;5340 if(newVertexOnBThEdge) delete [] newVertexOnBThEdge;5341 5342 return ret; // ok5343 }5344 /*}}}1*/5345 /*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{1*/5346 long Mesh::SplitInternalEdgeWithBorderVertices(){5347 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/5348 5349 long NbSplitEdge=0;5350 SetVertexFieldOn();5351 long it;5352 long nbvold=nbv;5353 long int verbose=2;5354 for (it=0;it<nbt;it++){5355 Triangle &t=triangles[it];5356 if (t.link)5357 for (int j=0;j<3;j++)5358 if(!t.Locked(j) && !t.Hidden(j)){5359 Triangle &tt = *t.TriangleAdj(j);5360 if ( &tt && tt.link && it < Number(tt))5361 { // an internal edge5362 MeshVertex &v0 = t[VerticesOfTriangularEdge[j][0]];5363 MeshVertex &v1 = t[VerticesOfTriangularEdge[j][1]];5364 if (v0.onGeometry && v1.onGeometry){5365 R2 P= ((R2) v0 + (R2) v1)*0.5;5366 if ( nbv<nbvx) {5367 vertices[nbv].r = P;5368 vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);5369 vertices[nbv].ReferenceNumber=0;5370 vertices[nbv].DirOfSearch = NoDirOfSearch ;5371 }5372 NbSplitEdge++;5373 }5374 }5375 }5376 }5377 ReMakeTriangleContainingTheVertex();5378 if (nbvold!=nbv){5379 long iv = nbvold;5380 long NbSwap = 0;5381 Icoor2 dete[3];5382 for (int i=nbvold;i<nbv;i++) {// for all the new point5383 MeshVertex & vi = vertices[i];5384 vi.i = toI2(vi.r);5385 vi.r = toR2(vi.i);5386 5387 // a good new point5388 vi.ReferenceNumber=0;5389 vi.DirOfSearch =NoDirOfSearch;5390 Triangle *tcvi = FindTriangleContaining(vi.i,dete);5391 if (tcvi && !tcvi->link) {5392 printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");5393 }5394 5395 quadtree->Add(vi);5396 if (!tcvi || tcvi->det<0){// internal5397 ISSMERROR("!tcvi || tcvi->det < 0");5398 }5399 AddVertex(vi,tcvi,dete);5400 NbSwap += vi.Optim(1);5401 iv++;5402 // }5403 }5404 if (verbose>3) {5405 printf(" number of points: %i\n",iv);5406 printf(" number of swap to split internal edges with border vertices: %i\n",NbSwap);5407 nbv = iv;5408 }5409 }5410 if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold));5411 if (verbose>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge);5412 5413 return NbSplitEdge;5414 }5415 /*}}}1*/5416 5405 /*FUNCTION Mesh::TriangleReferenceList{{{1*/ 5417 5406 long Mesh::TriangleReferenceList(long* reft) const { … … 5484 5473 ISSMERROR("kkk>=1000"); 5485 5474 } 5486 MeshVertex &vI = *edge.EdgeVertex(0);5487 MeshVertex &vJ = *edge.EdgeVertex(1);5475 BamgVertex &vI = *edge.EdgeVertex(0); 5476 BamgVertex &vJ = *edge.EdgeVertex(1); 5488 5477 I2 I=vI, J=vJ, IJ= J-I; 5489 5478 IJ_IA = (IJ ,(A-I)); … … 5525 5514 // the probleme is in case of the fine and long internal hole 5526 5515 // for exemple neart the training edge of a wing 5527 MeshVertex * s=0,*s1=0, *s0=0;5516 BamgVertex * s=0,*s1=0, *s0=0; 5528 5517 Icoor2 imax = MaxICoor22; 5529 5518 Icoor2 l0 = imax,l1 = imax; … … 5612 5601 if ((link + linkp) == 1) 5613 5602 { // a boundary edge 5614 MeshVertex * st = edge.EdgeVertex(1);5603 BamgVertex * st = edge.EdgeVertex(1); 5615 5604 I2 I=*st; 5616 5605 Icoor2 ll = Norme2_2 (C-I); … … 5657 5646 /*}}}1*/ 5658 5647 /*FUNCTION ForceEdge{{{1*/ 5659 int ForceEdge( MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret) {5648 int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret) { 5660 5649 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 5661 5650 … … 5668 5657 5669 5658 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); 5670 MeshVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;5659 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 5671 5660 // we turn around a in the direct sens 5672 5661 … … 5693 5682 if((det1 < 0) && (det2 >0)) { 5694 5683 // try to force the edge 5695 MeshVertex * va = &a, *vb = &b;5684 BamgVertex * va = &a, *vb = &b; 5696 5685 tc = Previous(tc); 5697 5686 if (!v1 || !v2){ … … 5703 5692 ISSMERROR("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l); 5704 5693 } 5705 MeshVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);5694 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 5706 5695 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { 5707 5696 tc.SetLock(); … … 5733 5722 /*}}}1*/ 5734 5723 /*FUNCTION swap{{{1*/ 5735 void swap(Triangle *t1,short a1, Triangle *t2,short a2, MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2){5724 void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){ 5736 5725 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/ 5737 5726 // -------------------------------------------------------------- … … 5778 5767 /*}}}1*/ 5779 5768 /*FUNCTION SwapForForcingEdge{{{1*/ 5780 int SwapForForcingEdge( MeshVertex * & pva ,MeshVertex * & pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {5769 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) { 5781 5770 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ 5782 5771 // l'arete ta coupe l'arete pva pvb … … 5795 5784 } 5796 5785 5797 MeshVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];5798 MeshVertex & s1= (*t1)[OppositeVertex[a1]];5799 MeshVertex & s2= (*t2)[OppositeVertex[a2]];5786 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 5787 BamgVertex & s1= (*t1)[OppositeVertex[a1]]; 5788 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 5800 5789 5801 5790 -
issm/trunk/src/c/objects/Bamg/Mesh.h
r5095 r5120 22 22 23 23 class Mesh { 24 24 25 public: 25 26 26 Geometry & Gh;// Geometry27 Mesh & BTh; // Background Mesh Bth==*this =>no background28 long NbRef;// counter of ref on the this class if 0 we can delete29 long nbvx,nbtx; // nombre max de sommets , detriangles30 long nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,31 long NbOfQuad; // nb of quadrangle32 long NbSubDomains;33 long NbOutT;// Nb of oudeside triangle34 long NbOfTriangleSearchFind;35 long NbOfSwapTriangle;36 MeshVertex*vertices;37 long NbVerticesOnGeomVertex;38 VertexOnGeom *VerticesOnGeomVertex;39 long NbVerticesOnGeomEdge;40 VertexOnGeom *VerticesOnGeomEdge;41 long NbVertexOnBThVertex;42 VertexOnVertex *VertexOnBThVertex;43 long NbVertexOnBThEdge;44 VertexOnEdge *VertexOnBThEdge;45 long NbCrackedVertices;46 long *CrackedVertices;47 long NbCrackedEdges;48 CrackedEdge *CrackedEdges;49 R2 pmin,pmax;// extrema50 double coefIcoor;// coef to integer Icoor1;51 Triangle *triangles;52 Edge * edges;53 QuadTree *quadtree;54 MeshVertex**ordre;55 SubDomain *subdomains;56 ListofIntersectionTriangles lIntTria;27 Geometry & Gh; // Geometry 28 Mesh & BTh; // Background Mesh Bth== *this =>no background 29 long NbRef; // counter of ref on the this class if 0 we can delete 30 long maxnbv,nbtx; // nombre max de sommets , de triangles 31 long nt,nbv,nbt,nbe; // nb of legal triangles, nb of vertex, of triangles, of edges with reference, 32 long NbOfQuad; // nb of quadrangle 33 long NbSubDomains; 34 long NbOutT; // Nb of oudeside triangle 35 long NbOfTriangleSearchFind; 36 long NbOfSwapTriangle; 37 BamgVertex *vertices; 38 long NbVerticesOnGeomVertex; 39 VertexOnGeom *VerticesOnGeomVertex; 40 long NbVerticesOnGeomEdge; 41 VertexOnGeom *VerticesOnGeomEdge; 42 long NbVertexOnBThVertex; 43 VertexOnVertex *VertexOnBThVertex; 44 long NbVertexOnBThEdge; 45 VertexOnEdge *VertexOnBThEdge; 46 long NbCrackedVertices; 47 long *CrackedVertices; 48 long NbCrackedEdges; 49 CrackedEdge *CrackedEdges; 50 R2 pmin,pmax; // extrema 51 double coefIcoor; // coef to integer Icoor1; 52 Triangle *triangles; 53 Edge *edges; 54 QuadTree *quadtree; 55 BamgVertex **ordre; 56 SubDomain *subdomains; 57 ListofIntersectionTriangles lIntTria; 57 58 58 59 //Constructors/Destructors 59 60 Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts); 60 61 Mesh(double* index,double* x,double* y,int nods,int nels); 61 Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long nbvxx=0 ); //copy operator62 Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator 62 63 Mesh(const Mesh &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature 63 Mesh(long nbvx,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {64 try { GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}64 Mesh(long maxnbv,Mesh & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) { 65 try {TriangulateFromGeom1(maxnbv,bamgopts,keepBackVertices);} 65 66 catch(...) { this->~Mesh(); throw; } 66 67 } 67 Mesh(long nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){68 try { GeomToTriangles0(nbvx,bamgopts);}68 Mesh(long maxnbv,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){ 69 try { TriangulateFromGeom0(maxnbv,bamgopts);} 69 70 catch(...) { this->~Mesh(); throw; } 70 71 } … … 72 73 73 74 //Operators 74 const MeshVertex & operator[] (long i) const { return vertices[i];};75 MeshVertex & operator[](long i) { return vertices[i];};75 const BamgVertex & operator[] (long i) const { return vertices[i];}; 76 BamgVertex & operator[](long i) { return vertices[i];}; 76 77 const Triangle & operator() (long i) const { return triangles[i];}; 77 78 Triangle & operator()(long i) { return triangles[i];}; … … 87 88 return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y); 88 89 } 89 void AddVertex( MeshVertex & s,Triangle * t,Icoor2 * =0) ;90 void AddVertex(BamgVertex & s,Triangle * t,Icoor2 * =0) ; 90 91 void Insert(); 91 92 void ForceBoundary(); … … 110 111 void SmoothingVertex(int =3,double=0.3); 111 112 Metric MetricAt (const R2 &) const; 112 GeometricalEdge* ProjectOnCurve( Edge & AB, MeshVertex & A, MeshVertex & B,double theta, MeshVertex & R,VertexOnEdge & BR,VertexOnGeom & GR);113 GeometricalEdge* ProjectOnCurve( Edge & AB, BamgVertex & A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 113 114 long Number(const Triangle & t) const { return &t - triangles;} 114 115 long Number(const Triangle * t) const { return t - triangles;} 115 long Number(const MeshVertex & t) const { return &t - vertices;}116 long Number(const MeshVertex * t) const { return t - vertices;}116 long Number(const BamgVertex & t) const { return &t - vertices;} 117 long Number(const BamgVertex * t) const { return t - vertices;} 117 118 long Number(const Edge & t) const { return &t - edges;} 118 119 long Number(const Edge * t) const { return t - edges;} 119 120 long Number2(const Triangle * t) const { return t - triangles; } 120 MeshVertex* NearestVertex(Icoor1 i,Icoor1 j) ;121 BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ; 121 122 Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; 122 123 void ReadMesh(double* index,double* x,double* y,int nods,int nels); … … 155 156 156 157 private: 157 void GeomToTriangles1(long nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption158 void GeomToTriangles0(long nbvx,BamgOpts* bamgopts);// the real constructor mesh generator159 void PreInit(long);158 void TriangulateFromGeom1(long maxnbv,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption 159 void TriangulateFromGeom0(long maxnbv,BamgOpts* bamgopts);// the real constructor mesh generator 160 void Init(long); 160 161 161 162 }; … … 166 167 void swap(Triangle *t1,short a1, 167 168 Triangle *t2,short a2, 168 MeshVertex *s1,MeshVertex *s2,Icoor2 det1,Icoor2 det2);169 int SwapForForcingEdge( MeshVertex * & pva ,MeshVertex * & pvb ,169 BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2); 170 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , 170 171 TriangleAdjacent & tt1,Icoor2 & dets1, 171 172 Icoor2 & detsa,Icoor2 & detsb, int & nbswap); 172 int ForceEdge( MeshVertex &a, MeshVertex & b,TriangleAdjacent & taret) ;173 int ForceEdge(BamgVertex &a, BamgVertex & b,TriangleAdjacent & taret) ; 173 174 inline TriangleAdjacent Previous(const TriangleAdjacent & ta){ 174 175 return TriangleAdjacent(ta.t,PreviousEdge[ta.a]); … … 183 184 int j=i;i=on->DirAdj[i];on=on->Adj[j]; 184 185 } 185 inline double qualite(const MeshVertex &va,const MeshVertex &vb,const MeshVertex &vc){186 inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){ 186 187 double ret; 187 188 I2 ia=va,ib=vb,ic=vc; -
issm/trunk/src/c/objects/Bamg/QuadTree.cpp
r5095 r5120 24 24 /*}}}*/ 25 25 /*DOCUMENTATION What is a QuadTree? {{{1 26 * A Quadtree is a very simple way to group thevertices according27 * to their location . A square that holds all the points of the mesh28 * (or the geometry) is divided into 4 boxes. As soo mas one box26 * A Quadtree is a very simple way to group vertices according 27 * to their locations. A square that holds all the points of the mesh 28 * (or the geometry) is divided into 4 boxes. As soon as one box 29 29 * hold more than 4 vertices, it is divided into 4 new boxes, etc... 30 30 * There cannot be more than MAXDEEP (=30) subdivision. … … 67 67 * 0 1 1 1 .... 1 0 1 1 1 .... 1 in binary 68 68 * \-- 29 --/ \-- 29 --/ 69 * Using the binaries is thereforvery easy to locate a vertex in a box:69 * Using binaries is therefore very easy to locate a vertex in a box: 70 70 * we just need to look at the bits from the left to the right (See ::Add) 71 71 }}}1*/ … … 76 76 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/QuadTree)*/ 77 77 78 lenStorageQuadTreeBox(t-> nbvx/8+10),78 lenStorageQuadTreeBox(t->maxnbv/8+10), 79 79 th(t), 80 80 NbQuadTreeBox(0), … … 120 120 /*Methods*/ 121 121 /*FUNCTION QuadTree::Add{{{1*/ 122 void QuadTree::Add( MeshVertex &w){122 void QuadTree::Add(BamgVertex &w){ 123 123 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/Add)*/ 124 124 … … 163 163 164 164 //Copy the 4 vertices in the current QuadTreebox 165 MeshVertex* v4[4];165 BamgVertex* v4[4]; 166 166 v4[0]= b->v[0]; 167 167 v4[1]= b->v[1]; … … 205 205 /*}}}1*/ 206 206 /*FUNCTION QuadTree::NearestVertex{{{1*/ 207 MeshVertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) {207 BamgVertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { 208 208 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertex)*/ 209 209 … … 224 224 Icoor1 jplus( j<MaxISize?(j<0?0:j):MaxISize-1); 225 225 //initial nearest vertex pointer 226 MeshVertex* vn=NULL;226 BamgVertex* vn=NULL; 227 227 228 228 //Get initial Quadtree box (largest) … … 293 293 294 294 //if the current subbox is holding vertices, 295 if (b->n>0){ // MeshVertex QuadTreeBox not empty295 if (b->n>0){ // BamgVertex QuadTreeBox not empty 296 296 NbVerticesSearch++; 297 297 I2 i2 = b->v[k]->i; … … 344 344 /*}}}1*/ 345 345 /*FUNCTION QuadTree::NearestVertexWithNormal{{{1*/ 346 MeshVertex* QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) {346 BamgVertex* QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) { 347 347 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/NearestVertexWithNormal)*/ 348 348 … … 358 358 Icoor1 jplus( j<MaxISize?(j<0?0:j):MaxISize-1); 359 359 360 MeshVertex *vn=0;360 BamgVertex *vn=0; 361 361 362 362 // init for optimisation --- … … 411 411 int k = pi[l]; 412 412 413 if (b->n>0) // MeshVertex QuadTreeBox none empty413 if (b->n>0) // BamgVertex QuadTreeBox none empty 414 414 { 415 415 NbVerticesSearch++; … … 472 472 /*}}}1*/ 473 473 /*FUNCTION QuadTree::ToClose {{{1*/ 474 MeshVertex * QuadTree::ToClose(MeshVertex & v,double seuil,Icoor1 hx,Icoor1 hy){474 BamgVertex * QuadTree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){ 475 475 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, QuadTree.cpp/ToClose)*/ 476 476 … … 489 489 Icoor1 i0=0,j0=0; 490 490 491 // MeshVertex *vn=0;491 // BamgVertex *vn=0; 492 492 493 493 if (!root->n) … … 506 506 register int k = pi[l]; 507 507 508 if (b->n>0) // MeshVertex QuadTreeBox none empty508 if (b->n>0) // BamgVertex QuadTreeBox none empty 509 509 { 510 510 NbVerticesSearch++; -
issm/trunk/src/c/objects/Bamg/QuadTree.h
r5095 r5120 11 11 12 12 class Mesh; 13 class MeshVertex;13 class BamgVertex; 14 14 15 15 class QuadTree{ … … 18 18 public: 19 19 long n; 20 //contains only one object form the list (either MeshVertex or QuadTreeBox)21 // if n < 4 => MeshVertex else => QuadTreeBox;20 //contains only one object form the list (either BamgVertex or QuadTreeBox) 21 // if n < 4 => BamgVertex else => QuadTreeBox; 22 22 union{ 23 23 QuadTreeBox* b[4]; 24 MeshVertex* v[4];24 BamgVertex* v[4]; 25 25 }; 26 26 }; … … 50 50 long NbQuadTreeBox,NbVertices; 51 51 long NbQuadTreeBoxSearch,NbVerticesSearch; 52 MeshVertex* NearestVertex(Icoor1 i,Icoor1 j);53 MeshVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j);54 MeshVertex* ToClose(MeshVertex & ,double ,Icoor1,Icoor1);52 BamgVertex* NearestVertex(Icoor1 i,Icoor1 j); 53 BamgVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j); 54 BamgVertex* ToClose(BamgVertex & ,double ,Icoor1,Icoor1); 55 55 long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();} 56 void Add( MeshVertex & w);56 void Add( BamgVertex & w); 57 57 QuadTreeBox* NewQuadTreeBox(){ 58 58 if(! (sb->bc<sb->be)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); -
issm/trunk/src/c/objects/Bamg/R2.h
r3913 r5120 7 7 namespace bamg { 8 8 9 template <class R,class RR> class P2xP2; 10 template <class R,class RR> 11 class P2 { 9 template <class R,class RR> class P2 { 10 12 11 public: 12 13 13 //fields 14 14 R x,y; 15 15 16 //functions 16 17 P2 () :x(0),y(0) {}; … … 19 20 void Echo(){ 20 21 printf("Member of P2:\n"); 21 printf(" x: %g \n",x);22 printf(" y: %g \n",y);22 printf(" x: %g or %i\n",x,x); 23 printf(" y: %g or %i\n",y,y); 23 24 } 24 25 //operators … … 33 34 P2<R,RR> operator*=(const R r) {x *= r;y *= r;return *this;} 34 35 P2<R,RR> operator-=(const P2<R,RR> & cc) {x -= cc.x;y -= cc.y;return *this;} 36 35 37 }; 36 38 37 template <class R,class RR> 38 class P2xP2 { 39 template <class R,class RR> class P2xP2 { 40 39 41 private: 42 40 43 friend P2<R,RR> operator*(P2<R,RR> c,P2xP2<R,RR> cc){ 41 44 return P2<R,RR>(c.x*cc.x.x + c.y*cc.y.x, c.x*cc.x.y + c.y*cc.y.y); 42 45 } 43 46 public: 47 44 48 //fields 45 49 P2<R,RR> x,y; 50 46 51 //functions 47 52 P2xP2 (): x(),y() {} -
issm/trunk/src/c/objects/Bamg/Triangle.cpp
r5095 r5120 11 11 /*FUNCTION Triangle(Mesh *Th,long i,long j,long k) {{{1*/ 12 12 Triangle::Triangle(Mesh *Th,long i,long j,long k) { 13 MeshVertex *v=Th->vertices;13 BamgVertex *v=Th->vertices; 14 14 long nbv = Th->nbv; 15 15 if (i<0 || j<0 || k<0){ … … 27 27 } 28 28 /*}}}*/ 29 /*FUNCTION Triangle( MeshVertex *v0,MeshVertex *v1,MeshVertex *v2) {{{1*/30 Triangle::Triangle( MeshVertex *v0,MeshVertex *v1,MeshVertex *v2){29 /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{1*/ 30 Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){ 31 31 TriaVertices[0]=v0; 32 32 TriaVertices[1]=v1; … … 154 154 /*}}}1*/ 155 155 /*FUNCTION Triangle::Quadrangle {{{1*/ 156 Triangle* Triangle::Quadrangle( MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const{156 Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{ 157 157 // return the other triangle of the quad if a quad or 0 if not a quat 158 158 Triangle * t =0; … … 186 186 q= -1; 187 187 else if(option){ 188 const MeshVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];189 const MeshVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];190 const MeshVertex & v1 = *TriaVertices[OppositeEdge[a]];191 const MeshVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];188 const BamgVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]]; 189 const BamgVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]]; 190 const BamgVertex & v1 = *TriaVertices[OppositeEdge[a]]; 191 const BamgVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]]; 192 192 q = QuadQuality(v0,v1,v2,v3); // do the float part 193 193 } … … 220 220 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 221 221 222 register MeshVertex *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];223 register MeshVertex *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];224 register MeshVertex *s1=t1->TriaVertices[OppositeVertex[a1]];225 register MeshVertex *s2=t2->TriaVertices[OppositeVertex[a2]];222 register BamgVertex *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]]; 223 register BamgVertex *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]]; 224 register BamgVertex *s1=t1->TriaVertices[OppositeVertex[a1]]; 225 register BamgVertex *s2=t2->TriaVertices[OppositeVertex[a2]]; 226 226 227 227 Icoor2 det1=t1->det , det2=t2->det ; -
issm/trunk/src/c/objects/Bamg/Triangle.h
r5095 r5120 9 9 //classes 10 10 class Mesh; 11 class MeshVertex;11 class BamgVertex; 12 12 class Triangle; 13 13 … … 17 17 18 18 private: 19 MeshVertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer19 BamgVertex* TriaVertices[3]; // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer 20 20 Triangle* TriaAdjTriangles[3]; // 3 pointers toward the adjacent triangles 21 21 short TriaAdjSharedEdge[3]; // number of the edges in the adjacent triangles the edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1 … … 31 31 Triangle() {} 32 32 Triangle(Mesh *Th,long i,long j,long k); 33 Triangle( MeshVertex *v0,MeshVertex *v1,MeshVertex *v2);33 Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2); 34 34 35 35 //Operators 36 const MeshVertex & operator[](int i) const {return *TriaVertices[i];};37 MeshVertex & operator[](int i) {return *TriaVertices[i];};38 const MeshVertex * operator()(int i) const {return TriaVertices[i];};39 MeshVertex * & operator()(int i) {return TriaVertices[i];};36 const BamgVertex & operator[](int i) const {return *TriaVertices[i];}; 37 BamgVertex & operator[](int i) {return *TriaVertices[i];}; 38 const BamgVertex * operator()(int i) const {return TriaVertices[i];}; 39 BamgVertex * & operator()(int i) {return TriaVertices[i];}; 40 40 41 41 //Methods … … 52 52 TriangleAdjacent Adj(int i) const {return TriangleAdjacent(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);}; 53 53 Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];} 54 Triangle* Quadrangle( MeshVertex * & v0,MeshVertex * & v1,MeshVertex * & v2,MeshVertex * & v3) const ;54 Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ; 55 55 void ReNumbering(Triangle *tb,Triangle *te, long *renu){ 56 56 if (link >=tb && link <te) link = tb + renu[link -tb]; … … 59 59 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb]; 60 60 } 61 void ReNumbering( MeshVertex *vb,MeshVertex *ve, long *renu){61 void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){ 62 62 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb]; 63 63 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb]; … … 118 118 double qualite() ; 119 119 void Set(const Triangle &,const Mesh &,Mesh &); 120 int In( MeshVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}120 int In(BamgVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;} 121 121 122 122 }; -
issm/trunk/src/c/objects/Bamg/TriangleAdjacent.cpp
r5095 r5120 41 41 /*}}}*/ 42 42 /*FUNCTION TriangleAdjacent::EdgeVertex {{{1*/ 43 MeshVertex* TriangleAdjacent::EdgeVertex(const int & i) const {43 BamgVertex* TriangleAdjacent::EdgeVertex(const int & i) const { 44 44 return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; 45 45 } 46 46 /*}}}*/ 47 47 /*FUNCTION TriangleAdjacent::OppositeVertex {{{1*/ 48 MeshVertex* TriangleAdjacent::OppositeVertex() const {48 BamgVertex* TriangleAdjacent::OppositeVertex() const { 49 49 return t->TriaVertices[bamg::OppositeVertex[a]]; 50 50 } -
issm/trunk/src/c/objects/Bamg/TriangleAdjacent.h
r3913 r5120 3 3 4 4 #include "./include.h" 5 #include "./ MeshVertex.h"5 #include "./BamgVertex.h" 6 6 7 7 namespace bamg { … … 38 38 int swap(); 39 39 TriangleAdjacent Adj() const; 40 MeshVertex* EdgeVertex(const int & i) const;41 MeshVertex* OppositeVertex() const;40 BamgVertex* EdgeVertex(const int & i) const; 41 BamgVertex* OppositeVertex() const; 42 42 Icoor2& det() const; 43 43 }; -
issm/trunk/src/c/objects/Bamg/VertexOnEdge.h
r5095 r5120 9 9 //classes 10 10 class Mesh; 11 class MeshVertex;11 class BamgVertex; 12 12 13 13 class VertexOnEdge { 14 14 15 15 public: 16 MeshVertex* v;16 BamgVertex* v; 17 17 Edge* be; 18 18 double abcisse; 19 19 20 20 //Constructors 21 VertexOnEdge( MeshVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {}21 VertexOnEdge( BamgVertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {} 22 22 VertexOnEdge(){} 23 23 24 24 //Operators 25 25 operator double () const { return abcisse;} 26 operator MeshVertex* () const { return v;}26 operator BamgVertex* () const { return v;} 27 27 operator Edge* () const { return be;} 28 MeshVertex & operator[](int i) const { return (*be)[i];}28 BamgVertex & operator[](int i) const { return (*be)[i];} 29 29 30 30 //Methods -
issm/trunk/src/c/objects/Bamg/VertexOnGeom.h
r5095 r5120 9 9 //classes 10 10 class Mesh; 11 class MeshVertex;11 class BamgVertex; 12 12 class GeometricalEdge; 13 13 … … 16 16 public: 17 17 18 MeshVertex* mv;18 BamgVertex* mv; 19 19 double abscisse; 20 20 union{ … … 25 25 //Constructors/Destructors 26 26 VertexOnGeom(): mv(0),abscisse(0){gv=0;} 27 VertexOnGeom( MeshVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;}28 VertexOnGeom( MeshVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;}27 VertexOnGeom(BamgVertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} 28 VertexOnGeom(BamgVertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;} 29 29 30 30 //Operators 31 operator MeshVertex*() const {return mv;}31 operator BamgVertex*() const {return mv;} 32 32 operator GeometricalVertex * () const {return gv;} 33 33 operator GeometricalEdge * () const {return ge;} -
issm/trunk/src/c/objects/Bamg/VertexOnVertex.h
r5095 r5120 3 3 4 4 #include "./include.h" 5 #include "./ MeshVertex.h"5 #include "./BamgVertex.h" 6 6 7 7 namespace bamg { … … 13 13 14 14 public: 15 MeshVertex* v;16 MeshVertex* bv;15 BamgVertex* v; 16 BamgVertex* bv; 17 17 18 18 //Constructors 19 VertexOnVertex( MeshVertex * w,MeshVertex *bw) :v(w),bv(bw){}19 VertexOnVertex(BamgVertex * w,BamgVertex *bw) :v(w),bv(bw){} 20 20 VertexOnVertex() {}; 21 21 -
issm/trunk/src/c/objects/Bamg/typedefs.h
r3913 r5120 6 6 namespace bamg { 7 7 8 typedef int Icoor1; 9 #if LONG_BIT > 63 //64 bits or more 8 /*Integer coordinates types*/ 9 typedef int Icoor1; 10 #if LONG_BIT > 63 //64 bits or more 10 11 typedef long Icoor2; 11 #else //32 bits12 #else //32 bits 12 13 typedef double Icoor2; 13 #endif 14 #endif 15 16 /*I2 and R2*/ 14 17 typedef P2<Icoor1,Icoor2> I2; 15 18 typedef P2xP2<short,long> I2xI2; -
issm/trunk/src/c/objects/objects.h
r5095 r5120 99 99 #include "./Bamg/DoubleAndInt.h" 100 100 #include "./Bamg/Direction.h" 101 #include "./Bamg/ MeshVertex.h"101 #include "./Bamg/BamgVertex.h" 102 102 #include "./Bamg/TriangleAdjacent.h" 103 103 #include "./Bamg/Edge.h"
Note:
See TracChangeset
for help on using the changeset viewer.