Ice Sheet System Model  4.18
Code documentation
Public Member Functions | Data Fields | Private Member Functions
bamg::Mesh Class Reference

#include <Mesh.h>

Public Member Functions

 Mesh (BamgGeom *bamggeom, BamgMesh *bamgmesh, BamgOpts *bamgopts)
 
 Mesh (int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts)
 
 Mesh (double *x, double *y, int nods, BamgOpts *bamgopts)
 
 Mesh (Mesh &, Geometry *pGh=0, Mesh *pBTh=0, long maxnbv_in=0)
 
 Mesh (const Mesh &, const int *flag, const int *bb, BamgOpts *bamgopts)
 
 Mesh (long maxnbv, Mesh &BT, BamgOpts *bamgopts, int keepBackVertices=1)
 
 Mesh (long maxnbv, Geometry &G, BamgOpts *bamgopts)
 
 ~Mesh ()
 
const BamgVertexoperator[] (long i) const
 
BamgVertexoperator[] (long i)
 
const Triangleoperator() (long i) const
 
Triangleoperator() (long i)
 
void SetIntCoor (const char *from=0)
 
double MinimalHmin ()
 
double MaximalHmax ()
 
I2 R2ToI2 (const R2 &P) const
 
R2 I2ToR2 (const I2 &P) const
 
void AddVertex (BamgVertex &s, Triangle *t, long long *=0)
 
void Insert (BamgOpts *bamgopts)
 
void Echo (void)
 
void ForceBoundary (BamgOpts *bamgopts)
 
void FindSubDomain (BamgOpts *bamgopts, int OutSide=0)
 
long TriangleReferenceList (long *) const
 
void TriangleIntNumbering (long *renumbering)
 
void CrackMesh (BamgOpts *bamgopts)
 
void SmoothMetric (BamgOpts *bamgopts, double raisonmax)
 
void BoundAnisotropy (BamgOpts *bamgopts, double anisomax, double hminaniso=1e-100)
 
Edge ** MakeGeomEdgeToEdge ()
 
long SplitInternalEdgeWithBorderVertices ()
 
void MakeBamgQuadtree ()
 
void MaxSubDivision (BamgOpts *bamgopts, double maxsubdiv)
 
void NewPoints (Mesh &, BamgOpts *bamgopts, int KeepVertices=1)
 
long InsertNewPoints (long nbvold, long &NbTSwap, BamgOpts *bamgopts)
 
void TrianglesRenumberBySubDomain (bool justcompress=false)
 
void SmoothingVertex (BamgOpts *bamgopts, int=3, double=0.3)
 
Metric MetricAt (const R2 &)
 
GeomEdgeProjectOnCurve (Edge &AB, BamgVertex &A, BamgVertex &B, double theta, BamgVertex &R, VertexOnEdge &BR, VertexOnGeom &GR)
 
long GetId (const Triangle &t) const
 
long GetId (const Triangle *t) const
 
long GetId (const BamgVertex &t) const
 
long GetId (const BamgVertex *t) const
 
long GetId (const Edge &t) const
 
long GetId (const Edge *t) const
 
BamgVertexNearestVertex (int i, int j)
 
TriangleTriangleFindFromCoord (const I2 &, long long[3], Triangle *tstart=0)
 
void ReadMesh (int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts)
 
void ReadMesh (BamgMesh *bamgmesh, BamgOpts *bamgopts)
 
void WriteMesh (BamgMesh *bamgmesh, BamgOpts *bamgopts)
 
void ReadMetric (const BamgOpts *bamgopts)
 
void WriteMetric (BamgOpts *bamgopts)
 
void WriteIndex (int **pindex, int *pnels)
 
void AddMetric (BamgOpts *bamgopts)
 
void BuildMetric0 (BamgOpts *bamgopts)
 
void BuildMetric1 (BamgOpts *bamgopts)
 
void BuildGeometryFromMesh (BamgOpts *bamgopts=NULL)
 
long RandomNumber (long max)
 
void ReconstructExistingMesh (BamgOpts *bamgopts)
 
void CreateSingleVertexToTriangleConnectivity ()
 
void UnMarkUnSwapTriangle ()
 
void SetVertexFieldOn ()
 
void SetVertexFieldOnBTh ()
 

Data Fields

GeometryGh
 
MeshBTh
 
BamgVertexvertices
 
Triangletriangles
 
Edgeedges
 
BamgQuadtreequadtree
 
BamgVertex ** orderedvertices
 
SubDomainsubdomains
 
long NbRef
 
long maxnbv
 
long maxnbt
 
long nbv
 
long nbt
 
long nbe
 
long nbsubdomains
 
long nbtout
 
R2 pmin
 
R2 pmax
 
double coefIcoor
 
ListofIntersectionTriangles lIntTria
 
long randomseed
 
long NbVerticesOnGeomVertex
 
VertexOnGeomVerticesOnGeomVertex
 
long NbVerticesOnGeomEdge
 
VertexOnGeomVerticesOnGeomEdge
 
long NbVertexOnBThVertex
 
VertexOnVertexVertexOnBThVertex
 
long NbVertexOnBThEdge
 
VertexOnEdgeVertexOnBThEdge
 
long NbCrackedVertices
 
long * CrackedVertices
 
long NbCrackedEdges
 
CrackedEdgeCrackedEdges
 

Private Member Functions

void TriangulateFromGeom1 (BamgOpts *bamgopts, int KeepVertices=1)
 
void TriangulateFromGeom0 (BamgOpts *bamgopts)
 
void Triangulate (double *x, double *y, int nods, BamgOpts *bamgopts)
 
void Init (long)
 
int ForceEdge (BamgVertex &a, BamgVertex &b, AdjacentTriangle &taret)
 
int SwapForForcingEdge (BamgVertex *&pva, BamgVertex *&pvb, AdjacentTriangle &tt1, long long &dets1, long long &detsa, long long &detsb, int &nbswap)
 

Detailed Description

Definition at line 21 of file Mesh.h.

Constructor & Destructor Documentation

◆ Mesh() [1/7]

bamg::Mesh::Mesh ( BamgGeom bamggeom,
BamgMesh bamgmesh,
BamgOpts bamgopts 
)

Definition at line 13 of file Mesh.cpp.

13  :Gh(*(new Geometry())),BTh(*this){ /*{{{*/
14 
15  /*Initialize fields*/
16  Init(0);
17 
18  /*Read Geometry if provided*/
19  if(bamggeom->Edges) {
20  Gh.ReadGeometry(bamggeom,bamgopts);
21  Gh.PostRead();
22  }
23 
24  /*Read background mesh*/
25  ReadMesh(bamgmesh,bamgopts);
26 
27  /*Build Geometry if not provided*/
28  if(bamggeom->Edges==NULL) {
29  /*Recreate geometry if needed*/
30  _printf_("WARNING: mesh present but no geometry found. Reconstructing...\n");
31  BuildGeometryFromMesh(bamgopts);
32  Gh.PostRead(true);
33  }
34 
35  /*Set integer coordinates*/
36  SetIntCoor();
37 
38  /*Fill holes and generate mesh properties*/
39  ReconstructExistingMesh(bamgopts);
40  }

◆ Mesh() [2/7]

bamg::Mesh::Mesh ( int *  index,
double *  x,
double *  y,
int  nods,
int  nels,
BamgOpts bamgopts 
)

Definition at line 42 of file Mesh.cpp.

42  :Gh(*(new Geometry())),BTh(*this){/*{{{*/
43 
44  Init(0);
45  ReadMesh(index,x,y,nods,nels,bamgopts);
46  SetIntCoor();
47  ReconstructExistingMesh(bamgopts);
48  }

◆ Mesh() [3/7]

bamg::Mesh::Mesh ( double *  x,
double *  y,
int  nods,
BamgOpts bamgopts 
)

Definition at line 50 of file Mesh.cpp.

50  :Gh(*(new Geometry())),BTh(*this){/*{{{*/
51  Triangulate(x,y,nods,bamgopts);
52  }

◆ Mesh() [4/7]

bamg::Mesh::Mesh ( Mesh Th,
Geometry pGh = 0,
Mesh pBTh = 0,
long  maxnbv_in = 0 
)

Definition at line 153 of file Mesh.cpp.

154  : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) {
155  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/
156  Gh.NbRef++;
157  maxnbv_in = Max(maxnbv_in,Th.nbv);
158  long i;
159  // do all the allocation to be sure all the pointer existe
160 
161  Init(maxnbv_in);// to make the allocation
162  // copy of triangles
163  nbv = Th.nbv;
164  nbt = Th.nbt;
165  nbe = Th.nbe;
166  nbsubdomains = Th.nbsubdomains;
167  nbtout = Th.nbtout;
168  NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex;
170  VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
171  NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge;
173  VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ;
174  if (& BTh == & Th.BTh){ // same background
175  BTh.NbRef++;
176  NbVertexOnBThVertex = Th.NbVertexOnBThVertex;
178  VertexOnBThVertex = new VertexOnVertex[NbVertexOnBThVertex];
179  NbVertexOnBThEdge = Th.NbVertexOnBThEdge;
181  VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge];
182  }
183  else { // no add on background mesh
184  BTh.NbRef++;
188  VertexOnBThEdge=0;
189  }
190 
191  if(nbe)
192  edges = new Edge[nbe];
193  if(nbsubdomains)
194  subdomains = new SubDomain[nbsubdomains];
195  pmin = Th.pmin;
196  pmax = Th.pmax;
197  coefIcoor = Th.coefIcoor;
198  for(i=0;i<nbt;i++)
199  triangles[i].Set(Th.triangles[i],Th,*this);
200  for(i=0;i<nbe;i++)
201  edges[i].Set(Th,i,*this);
202  for(i=0;i<nbv;i++)
203  vertices[i].Set(Th.vertices[i],Th,*this);
204  for(i=0;i<nbsubdomains;i++)
205  subdomains[i].Set(Th,i,*this);
206  for (i=0;i<NbVerticesOnGeomVertex;i++)
207  VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this);
208  for (i=0;i<NbVerticesOnGeomEdge;i++)
209  VerticesOnGeomEdge[i].Set(Th.VerticesOnGeomEdge[i],Th,*this);
210  quadtree=0;
211 
212  }

◆ Mesh() [5/7]

bamg::Mesh::Mesh ( const Mesh Tho,
const int *  flag,
const int *  bb,
BamgOpts bamgopts 
)

Definition at line 54 of file Mesh.cpp.

54  : Gh(*(new Geometry())), BTh(*this) {/*{{{*/
55  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/
56 
57  int i,k,itadj;
58  int kt=0;
59  int * kk = new int [Tho.nbv];
60  long * reft = new long[Tho.nbt];
61  long nbInT = Tho.TriangleReferenceList(reft);
62  long * refv = new long[Tho.nbv];
63 
64  for (i=0;i<Tho.nbv;i++)
65  kk[i]=-1;
66  for (i=0;i<Tho.nbv;i++)
67  refv[i]=0;
68  int nbNewBedge =0;
69  // int nbOldBedge =0;
70  for (i=0;i<Tho.nbt;i++)
71  if( reft[i] >=0 && flag[i])
72  {
73  const Triangle & t = Tho.triangles[i];
74  kt++;
75  kk[Tho.GetId(t[0])]=1;
76  kk[Tho.GetId(t[1])]=1;
77  kk[Tho.GetId(t[2])]=1;
78  itadj=Tho.GetId(t.TriangleAdj(0));
79  if ( reft[itadj] >=0 && !flag[itadj])
80  { nbNewBedge++;
81  refv[Tho.GetId(t[VerticesOfTriangularEdge[0][0]])]=bb[i];
82  refv[Tho.GetId(t[VerticesOfTriangularEdge[0][1]])]=bb[i];
83  }
84  itadj=Tho.GetId(t.TriangleAdj(1));
85  if ( reft[itadj] >=0 && !flag[itadj])
86  { nbNewBedge++;
87  refv[Tho.GetId(t[VerticesOfTriangularEdge[1][0]])]=bb[i];
88  refv[Tho.GetId(t[VerticesOfTriangularEdge[1][1]])]=bb[i];}
89  itadj=Tho.GetId(t.TriangleAdj(2));
90  if ( reft[itadj] >=0 && !flag[itadj])
91  { nbNewBedge++;
92  refv[Tho.GetId(t[VerticesOfTriangularEdge[2][0]])]=bb[i];
93  refv[Tho.GetId(t[VerticesOfTriangularEdge[2][1]])]=bb[i];}
94  }
95  k=0;
96  for (i=0;i<Tho.nbv;i++){
97  if (kk[i]>=0) kk[i]=k++;
98  }
99  _printf_(" number of vertices " << k << ", remove = " << Tho.nbv - k << "\n");
100  _printf_(" number of triangles " << kt << ", remove = " << nbInT-kt << "\n");
101  _printf_(" number of New boundary edge " << nbNewBedge << "\n");
102  long imaxnbv =k;
103  Init(imaxnbv);
104  for (i=0;i<Tho.nbv;i++)
105  if (kk[i]>=0)
106  {
107  vertices[nbv] = Tho.vertices[i];
108  if (!vertices[nbv].GetReferenceNumber())
109  vertices[nbv].ReferenceNumber = refv[i];
110  nbv++;
111  }
112  if (imaxnbv != nbv){
113  delete [] kk;
114  delete [] refv;
115  _error_("imaxnbv != nbv");
116  }
117  for (i=0;i<Tho.nbt;i++)
118  if(reft[i] >=0 && flag[i]){
119  const Triangle & t = Tho.triangles[i];
120  int i0 = Tho.GetId(t[0]);
121  int i1 = Tho.GetId(t[1]);
122  int i2 = Tho.GetId(t[2]);
123  if(i0<0 || i1<0 || i2<0){
124  delete [] refv;
125  _error_("i0<0 || i1<0 || i2< 0");
126  }
127  if(i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
128  delete [] refv;
129  _error_("i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv");
130  }
131  triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
132  triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
133  nbt++;
134  }
135  if (nbt==0 && nbv==0){
136  delete [] refv;
137  _error_("All triangles have been removed");
138  }
139  delete [] kk;
140  delete [] reft;
141  delete [] refv;
142  BuildGeometryFromMesh(bamgopts);
143  Gh.PostRead();
144  SetIntCoor();
145  ReconstructExistingMesh(bamgopts);
146 
147  /*Final checks*/
148  _assert_(kt==nbt);
150  _assert_(subdomains[0].head && subdomains[0].head->link);
151  }

◆ Mesh() [6/7]

bamg::Mesh::Mesh ( long  maxnbv,
Mesh BT,
BamgOpts bamgopts,
int  keepBackVertices = 1 
)

Definition at line 214 of file Mesh.cpp.

214  :Gh(BT.Gh),BTh(BT) {/*{{{*/
215  this->Init(imaxnbv);
216  TriangulateFromGeom1(bamgopts,keepBackVertices);
217  }

◆ Mesh() [7/7]

bamg::Mesh::Mesh ( long  maxnbv,
Geometry G,
BamgOpts bamgopts 
)

Definition at line 219 of file Mesh.cpp.

219  :Gh(G),BTh(*this){/*{{{*/
220  Init(imaxnbv);
221  TriangulateFromGeom0(bamgopts);
222  }

◆ ~Mesh()

bamg::Mesh::~Mesh ( )

Definition at line 224 of file Mesh.cpp.

224  {/*{{{*/
225  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Triangles)*/
226 
227  if(vertices) delete [] vertices;
228  if(edges) delete [] edges;
229  if(triangles) delete [] triangles;
230  if(quadtree) delete quadtree;
231  if(orderedvertices) delete [] orderedvertices;
232  if(subdomains) delete [] subdomains;
235  if(VertexOnBThVertex) delete [] VertexOnBThVertex;
236  if(VertexOnBThEdge) delete [] VertexOnBThEdge;
237 
238  if (Gh.NbRef>0) Gh.NbRef--;
239  else if (Gh.NbRef==0) delete &Gh;
240  if(&BTh != this){
241  if (BTh.NbRef>0) BTh.NbRef--;
242  else if (BTh.NbRef==0) delete &BTh;
243  }
244  Init(0); // set all to zero
245  }

Member Function Documentation

◆ operator[]() [1/2]

const BamgVertex& bamg::Mesh::operator[] ( long  i) const
inline

Definition at line 68 of file Mesh.h.

68 { return vertices[i]; };

◆ operator[]() [2/2]

BamgVertex& bamg::Mesh::operator[] ( long  i)
inline

Definition at line 69 of file Mesh.h.

69 { return vertices[i]; };

◆ operator()() [1/2]

const Triangle& bamg::Mesh::operator() ( long  i) const
inline

Definition at line 70 of file Mesh.h.

70 { return triangles[i]; };

◆ operator()() [2/2]

Triangle& bamg::Mesh::operator() ( long  i)
inline

Definition at line 71 of file Mesh.h.

71 { return triangles[i]; };

◆ SetIntCoor()

void bamg::Mesh::SetIntCoor ( const char *  from = 0)

Definition at line 3678 of file Mesh.cpp.

3678  {/*{{{*/
3679  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SetIntCoor)*/
3680 
3681  /*Set integer coordinate for existing vertices*/
3682 
3683  //Get extrema coordinates of the existing vertices
3684  pmin = vertices[0].r;
3685  pmax = vertices[0].r;
3686  long i;
3687  for (i=0;i<nbv;i++) {
3688  pmin.x = Min(pmin.x,vertices[i].r.x);
3689  pmin.y = Min(pmin.y,vertices[i].r.y);
3690  pmax.x = Max(pmax.x,vertices[i].r.x);
3691  pmax.y = Max(pmax.y,vertices[i].r.y);
3692  }
3693  R2 DD = (pmax-pmin)*0.05;
3694  pmin = pmin-DD;
3695  pmax = pmax+DD;
3696 
3697  //Compute coefIcoor
3698  long MaxICoord = 1073741823; //2^30 - 1 = =111...111 (29 times one)
3699  coefIcoor= (MaxICoord)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
3700  if (coefIcoor<=0){
3701  _error_("coefIcoor should be positive, a problem in the geometry is likely");
3702  }
3703 
3704  // generation of integer coord
3705  for (i=0;i<nbv;i++) {
3706  vertices[i].i = R2ToI2(vertices[i].r);
3707  }
3708 
3709  // computation of the det
3710  int number_of_errors=0;
3711  for (i=0;i<nbt;i++) {
3712  BamgVertex* v0 = triangles[i](0);
3713  BamgVertex* v1 = triangles[i](1);
3714  BamgVertex* v2 = triangles[i](2);
3715 
3716  //If this is not a boundary triangle
3717  if (v0 && v1 && v2){
3718 
3719  /*Compute determinant*/
3720  triangles[i].det= det(v0->GetIntegerCoordinates(),v1->GetIntegerCoordinates(),v2->GetIntegerCoordinates());
3721 
3722  /*Check that determinant is positive*/
3723  if (triangles[i].det <=0){
3724 
3725  /*increase number_of_errors and print error only for the first 20 triangles*/
3726  number_of_errors++;
3727  if (number_of_errors<20){
3728  _printf_("Area of Triangle " << i+1 << " < 0 (det=" << triangles[i].det << ")\n");
3729  }
3730  }
3731  }
3732 
3733  //else, set as -1
3734  else triangles[i].det=-1;
3735  }
3736 
3737  if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
3738  }

◆ MinimalHmin()

double bamg::Mesh::MinimalHmin ( )

Definition at line 3031 of file Mesh.cpp.

3031  {/*{{{*/
3032  return 2.0/coefIcoor;
3033  }

◆ MaximalHmax()

double bamg::Mesh::MaximalHmax ( )

Definition at line 2949 of file Mesh.cpp.

2949  {/*{{{*/
2950  return Max(pmax.x-pmin.x,pmax.y-pmin.y);
2951  }

◆ R2ToI2()

I2 bamg::Mesh::R2ToI2 ( const R2 P) const

Definition at line 3937 of file Mesh.cpp.

3937  {/*{{{*/
3938  return I2( (int) (coefIcoor*(P.x-pmin.x)),(int) (coefIcoor*(P.y-pmin.y)) );
3939  }

◆ I2ToR2()

R2 bamg::Mesh::I2ToR2 ( const I2 P) const

Definition at line 3941 of file Mesh.cpp.

3941  {/*{{{*/
3942  return R2( (double) P.x/coefIcoor+pmin.x, (double) P.y/coefIcoor+pmin.y);
3943  }

◆ AddVertex()

void bamg::Mesh::AddVertex ( BamgVertex s,
Triangle t,
long long *  det3 = 0 
)

Definition at line 1020 of file Mesh.cpp.

1020  {/*{{{*/
1021  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/
1022  // -------------------------------
1023  // s2
1024  // !
1025  // /|\ !
1026  // / | \ !
1027  // / | \ !
1028  // tt1 / | \ tt0 !
1029  // / |s \ !
1030  // / . \ !
1031  // / . ` \ !
1032  // / . ` \ !
1033  // ---------------- !
1034  // s0 tt2 s1
1035  //-------------------------------
1036 
1037  /*Intermediaries*/
1038  Triangle* tt[3]; //the three triangles
1039  long long det3local[3]; //three determinants (integer)
1040  int nbzerodet =0; //number of zeros in det3
1041  int izerodet=-1; //egde containing the vertex s
1042  int iedge;
1043 
1044  /*three vertices of t*/
1045  BamgVertex* s0=t->GetVertex(0);
1046  BamgVertex* s1=t->GetVertex(1);
1047  BamgVertex* s2=t->GetVertex(2);
1048 
1049  //determinant of t
1050  long long detOld=t->det;
1051 
1052  /* infvertexindex = index of the infinite vertex (NULL)
1053  if no infinite vertex (NULL) infvertexindex=-1
1054  else if v_i is infinite, infvertexindex=i*/
1055  int infvertexindex = s0 ? ((s1? (s2?-1:2):1)):0;
1056 
1057  //some checks
1058  if(((infvertexindex <0 ) && (detOld <0)) || ((infvertexindex >=0) && (detOld >0)) ){
1059  _error_("inconsistent configuration (Contact ISSM developers)");
1060  }
1061 
1062  // if det3 does not exist, build it
1063  if (!det3){
1064  //allocate
1065  det3 = det3local;
1066  //if no infinite vertex
1067  if (infvertexindex<0 ) {
1068  det3[0]=bamg::det(s .GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates());
1069  det3[1]=bamg::det(s0->GetIntegerCoordinates(),s .GetIntegerCoordinates(),s2->GetIntegerCoordinates());
1070  det3[2]=bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates());}
1071  else {
1072  // one of &s1 &s2 &s0 is NULL
1073  det3[0]= s0 ? -1 : bamg::det(s.GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
1074  det3[1]= s1 ? -1 : bamg::det(s0->GetIntegerCoordinates(),s.GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
1075  det3[2]= s2 ? -1 : bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates()) ;
1076  }
1077  }
1078 
1079  if (!det3[0]) izerodet=0,nbzerodet++;
1080  if (!det3[1]) izerodet=1,nbzerodet++;
1081  if (!det3[2]) izerodet=2,nbzerodet++;
1082 
1083  //if nbzerodet>0, point s is on an egde or on a vertex
1084  if (nbzerodet>0){
1085  /*s is on an edge*/
1086  if (nbzerodet==1) {
1087  iedge = OppositeEdge[izerodet];
1088  AdjacentTriangle ta = t->Adj(iedge);
1089 
1090  /*if the point is one the boundary
1091  add the point in outside part */
1092  if (t->det>=0){ // inside triangle
1093  if (((Triangle*)ta)->det<0 ) {
1094  // add in outside triangle
1095  AddVertex(s,( Triangle *)ta);
1096  return;
1097  }
1098  }
1099  }
1100  else{
1101  _error_("Cannot add a vertex more than once. Check duplicates");
1102  }
1103  }
1104 
1105  // remove de MarkUnSwap edge
1106  t->SetUnMarkUnSwap(0);
1107  t->SetUnMarkUnSwap(1);
1108  t->SetUnMarkUnSwap(2);
1109 
1110  tt[0]= t;
1111  tt[1]= &triangles[nbt++];
1112  tt[2]= &triangles[nbt++];
1113 
1114  if (nbt>maxnbt) _error_("Not enough triangles");
1115 
1116  *tt[1]=*tt[2]=*t;
1117  tt[0]->link=tt[1];
1118  tt[1]->link=tt[2];
1119 
1120  (*tt[0])(OppositeVertex[0])=&s;
1121  (*tt[1])(OppositeVertex[1])=&s;
1122  (*tt[2])(OppositeVertex[2])=&s;
1123 
1124  tt[0]->det=det3[0];
1125  tt[1]->det=det3[1];
1126  tt[2]->det=det3[2];
1127 
1128  // update adj des triangles externe
1129  tt[0]->SetAdjAdj(0);
1130  tt[1]->SetAdjAdj(1);
1131  tt[2]->SetAdjAdj(2);
1132  // update des adj des 3 triangle interne
1133  const int i0 = 0;
1134  const int i1= NextEdge[i0];
1135  const int i2 = PreviousEdge[i0];
1136 
1137  tt[i0]->SetAdj2(i2,tt[i2],i0);
1138  tt[i1]->SetAdj2(i0,tt[i0],i1);
1139  tt[i2]->SetAdj2(i1,tt[i1],i2);
1140 
1141  tt[0]->SetSingleVertexToTriangleConnectivity();
1142  tt[1]->SetSingleVertexToTriangleConnectivity();
1143  tt[2]->SetSingleVertexToTriangleConnectivity();
1144 
1145  // swap if the point s is on a edge
1146  if(izerodet>=0) {
1147  int rswap=tt[izerodet]->swap(iedge);
1148 
1149  if (!rswap) {
1150  _error_("swap the point s is on a edge");
1151  }
1152  }
1153 
1154  }/*}}}*/

◆ Insert()

void bamg::Mesh::Insert ( BamgOpts bamgopts)

Definition at line 2679 of file Mesh.cpp.

2679  {/*{{{*/
2680  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/
2681 
2682  /*Insert points in the existing Geometry*/
2683 
2684  //Intermediary
2685  long i;
2686  int verbose=0;
2687 
2688  /*Get options*/
2689  if(bamgopts) verbose=bamgopts->verbose;
2690 
2691  //Display info
2692  if (verbose>2) _printf_(" Insert initial " << nbv << " vertices\n");
2693 
2694  //Compute integer coordinates for the existing vertices
2695  SetIntCoor();
2696 
2697  /*Now we want to build a list (orderedvertices) of the vertices in a random
2698  * order. To do so, we use the following method:
2699  *
2700  * From an initial k0 in [0 nbv[ random (vertex number)
2701  * the next k (vertex number) is computed using a big
2702  * prime number (PN>>nbv) following:
2703  *
2704  * k_{i+1} = k_i + PN [nbv]
2705  *
2706  * let's show that:
2707  *
2708  * for all j in [0 nbv[, there is a unique i in [0 nbv[ such that k_i=j
2709  *
2710  * Let's assume that there are 2 distinct j1 and j2 such that
2711  * k_j1 = k_j2
2712  *
2713  * This means that
2714  *
2715  * k0+j1*PN = k0+j2*PN [nbv]
2716  * (j1-j2)*PN =0 [nbv]
2717  * since PN is a prime number larger than nbv, and nbv!=1
2718  * j1-j2=0 [nbv]
2719  * BUT
2720  * j1 and j2 are in [0 nbv[ which is impossible.
2721  *
2722  * We hence have built a random list of nbv elements of
2723  * [0 nbv[ all distincts*/
2724 
2725  //Get Prime number
2726  const long PrimeNumber= BigPrimeNumber(nbv);
2727  long k0=this->RandomNumber(nbv);
2728  if (verbose>4) _printf_(" Prime Number = "<<PrimeNumber<<"\n");
2729  if (verbose>4) _printf_(" k0 = "<<k0<<"\n");
2730 
2731  //Build orderedvertices
2732  for (i=0; i<nbv; i++){
2733  orderedvertices[i]=&vertices[k0=(k0+PrimeNumber)%nbv];
2734  }
2735 
2736  /*Modify orderedvertices such that the first 3 vertices form a triangle*/
2737 
2738  //get first vertex i such that [0,1,i] are not aligned
2739  for (i=2; det(orderedvertices[0]->i,orderedvertices[1]->i,orderedvertices[i]->i)==0;){
2740  //if i is higher than nbv, it means that all the determinants are 0,
2741  //all vertices are aligned!
2742  if (++i>=nbv) _error_("all the vertices are aligned");
2743  }
2744  if (verbose>4) _printf_(" i = "<<i<<"\n");
2745  // exchange i et 2 in "orderedvertices" so that
2746  // the first 3 vertices are not aligned (real triangle)
2748 
2749  /*Take the first edge formed by the first two vertices and build
2750  * the initial simple mesh from this edge and 2 boundary triangles*/
2751 
2752  BamgVertex *v0=orderedvertices[0], *v1=orderedvertices[1];
2753 
2754  nbt = 2;
2755  triangles[0](0) = NULL;//infinite vertex
2756  triangles[0](1) = v0;
2757  triangles[0](2) = v1;
2758  triangles[1](0) = NULL;//infinite vertex
2759  triangles[1](2) = v0;
2760  triangles[1](1) = v1;
2761 
2762  //Build adjacence
2763  const int e0 = OppositeEdge[0];
2764  const int e1 = NextEdge[e0];
2765  const int e2 = PreviousEdge[e0];
2766  triangles[0].SetAdj2(e0, &triangles[1] ,e0);
2767  triangles[0].SetAdj2(e1, &triangles[1] ,e2);
2768  triangles[0].SetAdj2(e2, &triangles[1] ,e1);
2769 
2770  triangles[0].det = -1; //boundary triangle: det = -1
2771  triangles[1].det = -1; //boundary triangle: det = -1
2772 
2775 
2776  triangles[0].link=&triangles[1];
2777  triangles[1].link=&triangles[0];
2778 
2779  //build quadtree
2780  if (!quadtree) quadtree = new BamgQuadtree(this,0);
2781  quadtree->Add(*v0);
2782  quadtree->Add(*v1);
2783 
2784  /*Now, add the vertices One by One*/
2785  long NbSwap=0;
2786  if (verbose>3) _printf_(" Begining of insertion process...\n");
2787  if (verbose>4) _printf_(" nbv = "<<nbv<<"\n");
2788 
2789  for (long icount=2; icount<nbv; icount++) {
2790 
2791  //Get new vertex
2792  BamgVertex *newvertex=orderedvertices[icount];
2793 
2794  //Find the triangle in which newvertex is located
2795  long long det3[3];
2796  Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates)
2797 
2798  //Add newvertex to the quadtree
2799  quadtree->Add(*newvertex);
2800 
2801  //Add newvertex to the existing mesh
2802  AddVertex(*newvertex,tcvi,det3);
2803 
2804  //Make the mesh Delaunay around newvertex by swaping the edges
2805  NbSwap += newvertex->Optim(1,0);
2806  }
2807 
2808  //Display info
2809  if (verbose>3) {
2810  _printf_(" NbSwap of insertion: " << NbSwap << "\n");
2811  _printf_(" NbSwap/nbv: " << NbSwap/nbv << "\n");
2812  }
2813  }

◆ Echo()

void bamg::Mesh::Echo ( void  )

Definition at line 2296 of file Mesh.cpp.

2296  {/*{{{*/
2297 
2298  int i;
2299 
2300  _printf_("Mesh Echo:\n");
2301  _printf_(" nbv = " << nbv << "\n");
2302  _printf_(" nbt = " << nbt << "\n");
2303  _printf_(" nbe = " << nbe << "\n");
2304  _printf_(" index:\n");
2305  for (i=0;i<nbt;i++){
2306  _printf_(" " << setw(4) << i+1 << ": ["
2307  << setw(4) << (((BamgVertex*)triangles[i](0))?GetId(triangles[i][0])+1:0) << " "
2308  << setw(4) << (((BamgVertex*)triangles[i](1))?GetId(triangles[i][1])+1:0) << " "
2309  << setw(4) << (((BamgVertex*)triangles[i](2))?GetId(triangles[i][2])+1:0) << "]\n");
2310  }
2311  _printf_(" coordinates:\n");
2312  for (i=0;i<nbv;i++){
2313  _printf_(" " << setw(4) << i+1 << ": [" << vertices[i].r.x << " " << vertices[i].r.y << "] and [" << vertices[i].i.x << " " << vertices[i].i.y << "]\n");
2314  }
2315 
2316  }

◆ ForceBoundary()

void bamg::Mesh::ForceBoundary ( BamgOpts bamgopts)

Definition at line 2318 of file Mesh.cpp.

2318  {/*{{{*/
2319  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceBoundary)*/
2320 
2321  int verbose=0;
2322  int k=0;
2323  int nbfe=0,nbswp=0,Nbswap=0;
2324 
2325  /*Get options*/
2326  if(bamgopts) verbose=bamgopts->verbose;
2327 
2328  //display
2329  if (verbose > 2) _printf_(" ForceBoundary nb of edge: " << nbe << "\n");
2330 
2331  //check that there is no triangle with 0 determinant
2332  for (int t = 0; t < nbt; t++){
2333  if (!triangles[t].det) k++;
2334  }
2335  if (k!=0) {
2336  _error_("there is " << k << " triangles of mes = 0");
2337  }
2338 
2339  //Force Edges
2340  AdjacentTriangle ta(0,0);
2341  for (int i = 0; i < nbe; i++){
2342 
2343  //Force edge i
2344  nbswp = ForceEdge(edges[i][0],edges[i][1],ta);
2345  if (nbswp<0) k++;
2346  else Nbswap += nbswp;
2347 
2348  if (nbswp) nbfe++;
2349  if ( nbswp < 0 && k < 5){
2350  _error_("Missing Edge " << i << ", v0=" << GetId(edges[i][0]) << ",v1=" << GetId(edges[i][1]));
2351  }
2352  }
2353 
2354  if (k!=0) {
2355  _error_("There are " << k << " lost edges, the boundary might be crossing");
2356  }
2357  for (int j=0;j<nbv;j++){
2358  Nbswap += vertices[j].Optim(1,0);
2359  }
2360  if (verbose > 3) _printf_(" number of inforced edge = " << nbfe << ", number of swap= " << Nbswap << "\n");
2361  }

◆ FindSubDomain()

void bamg::Mesh::FindSubDomain ( BamgOpts bamgopts,
int  OutSide = 0 
)

Definition at line 2363 of file Mesh.cpp.

2363  {/*{{{*/
2364  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindSubDomain)*/
2365 
2366  int verbose=0;
2367 
2368  /*Get options*/
2369  if(bamgopts) verbose=bamgopts->verbose;
2370 
2371  if (verbose >2){
2372  if (OutSide) _printf_(" Find all external sub-domain\n");
2373  else _printf_(" Find all internal sub-domain\n");
2374  }
2375  short * HeapArete = new short[nbt];
2376  Triangle ** HeapTriangle = new Triangle* [nbt];
2377  Triangle *t,*t1;
2378  long k,it;
2379 
2380  /*No color by default*/
2381  for(int itt=0;itt<nbt;itt++) triangles[itt].link=0;
2382 
2383  long NbSubDomTot =0;
2384  for(it=0;it<nbt;it++) {
2385  if( !triangles[it].link ){
2386  t = triangles + it;
2387  NbSubDomTot++;; // new composante connexe
2388  long i = 0; // niveau de la pile
2389  t->link = t ; // sd forme d'un triangle cicular link
2390 
2391  HeapTriangle[i] =t ;
2392  HeapArete[i] = 3;
2393 
2394  while(i >= 0) // boucle sur la pile
2395  { while (HeapArete[i]--) // boucle sur les 3 aretes
2396  {
2397  int na = HeapArete[i];
2398  Triangle * tc = HeapTriangle[i]; // triangle courant
2399  if( ! tc->Locked(na)) // arete non frontiere
2400  {
2401  Triangle * ta = tc->TriangleAdj(na) ; // næ triangle adjacent
2402  if (ta->link == 0 ) // non deja chainer => on enpile
2403  {
2404  i++;
2405  ta->link = t->link ; // on chaine les triangles
2406  t->link = ta ; // d'un meme sous domaine
2407  HeapArete[i] = 3; // pour les 3 triangles adjacents
2408  HeapTriangle[i] = ta;
2409  }}
2410  } // deplie fin de boucle sur les 3 adjacences
2411  i--;
2412  }
2413  }
2414  }
2415 
2416  // supression de tous les sous domaine infini <=> contient le sommet NULL
2417  it =0;
2418  nbtout = 0;
2419  while (it<nbt) {
2420  if (triangles[it].link)
2421  {
2422  if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))
2423  {
2424  // infini triangle
2425  NbSubDomTot --;
2426  t=&triangles[it];
2427  nbtout--; // on fait un coup de trop.
2428  while (t){
2429  nbtout++;
2430  t1=t;
2431  t=t->link;
2432  t1->link=0;
2433  }
2434  }
2435  }
2436  it++;} // end while (it<nbt)
2437  if (nbt == nbtout || !NbSubDomTot) {
2438  delete [] HeapArete;
2439  delete [] HeapTriangle;
2440  _error_("The boundary is not close: all triangles are outside");
2441  }
2442 
2443  delete [] HeapArete;
2444  delete [] HeapTriangle;
2445 
2446  if (OutSide|| !Gh.subdomains || !Gh.nbsubdomains )
2447  { // No geom sub domain
2448  long i;
2449  if (subdomains) delete [] subdomains;
2450  subdomains = new SubDomain[ NbSubDomTot];
2451  nbsubdomains= NbSubDomTot;
2452  for ( i=0;i<nbsubdomains;i++) {
2453  subdomains[i].head=NULL;
2454  subdomains[i].ReferenceNumber=i+1;
2455  }
2456  long * mark = new long[nbt];
2457  for (it=0;it<nbt;it++)
2458  mark[it]=triangles[it].link ? -1 : -2;
2459 
2460  it =0;
2461  k = 0;
2462  while (it<nbt) {
2463  if (mark[it] == -1) {
2464  t1 = & triangles[it];
2465  t = t1->link;
2466  mark[it]=k;
2467  subdomains[k].head = t1;
2468  do {
2469  mark[GetId(t)]=k;
2470  t=t->link;
2471  } while (t!=t1);
2472  mark[it]=k++;}
2473  // else if(mark[it] == -2 ) triangles[it].Draw(999);
2474  it++;} // end white (it<nbt)
2475  if (k!=nbsubdomains){
2476  delete [] mark;
2477  _error_("k!=nbsubdomains");
2478  }
2479  if(OutSide){
2480  // to remove all the sub domain by parity adjacents
2481  // because in this case we have only the true boundary edge
2482  // so the boundary is manifold
2483  long nbk = nbsubdomains;
2484  while (nbk)
2485  for (it=0;it<nbt && nbk ;it++)
2486  for (int na=0;na<3 && nbk ;na++)
2487  {
2488  Triangle *ta = triangles[it].TriangleAdj(na);
2489  long kl = ta ? mark[GetId(ta)] : -2;
2490  long kr = mark[it];
2491  if(kr !=kl) {
2492  if (kl >=0 && subdomains[kl].ReferenceNumber <0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
2494  if (kr >=0 && subdomains[kr].ReferenceNumber <0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
2496  if(kr<0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
2497  nbk--,subdomains[kl].ReferenceNumber=-1;
2498  if(kl<0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
2499  nbk--,subdomains[kr].ReferenceNumber=-1;
2500  }
2501  }
2502  long j=0;
2503  for ( i=0;i<nbsubdomains;i++)
2504  if((-subdomains[i].ReferenceNumber) %2) { // good
2505  if(i != j)
2507  j++;}
2508  else{
2509  t= subdomains[i].head;
2510  while (t){
2511  nbtout++;
2512  t1=t;
2513  t=t->link;
2514  t1->link=0;
2515  }//while (t)
2516  }
2517  if(verbose>4) _printf_(" Number of removes subdomains (OutSideMesh) = " << nbsubdomains-j << "\n");
2518  nbsubdomains=j;
2519  }
2520 
2521  delete [] mark;
2522 
2523  }
2524  else{ // find the head for all subdomains
2526  delete [] subdomains, subdomains=0;
2527  if (! subdomains )
2528  subdomains = new SubDomain[ Gh.nbsubdomains];
2531  long * mark = new long[nbt];
2532  Edge **GeomEdgetoEdge = MakeGeomEdgeToEdge();
2533 
2534  for (it=0;it<nbt;it++)
2535  mark[it]=triangles[it].link ? -1 : -2;
2536  long inew =0;
2537  for (int i=0;i<nbsubdomains;i++) {
2538  GeomEdge &eg = *Gh.subdomains[i].edge;
2540  // by carefull is not easy to find a edge create from a GeomEdge
2541  // see routine MakeGeomEdgeToEdge
2542  Edge &e = *GeomEdgetoEdge[Gh.GetId(eg)];
2543  _assert_(&e);
2544  BamgVertex * v0 = e(0),*v1 = e(1);
2545  Triangle *t = v0->t;
2546  int direction = Gh.subdomains[i].direction;
2547  // test if ge and e is in the same direction
2548  if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) direction = -direction ;
2549  subdomains[i].direction = direction;
2550  subdomains[i].edge = &e;
2551  _assert_(t && direction);
2552 
2553  AdjacentTriangle ta(t,EdgesVertexTriangle[v0->IndexInTriangle][0]);// previous edges
2554 
2555  while (1) {
2556  _assert_(v0==ta.EdgeVertex(1));
2557  if (ta.EdgeVertex(0) == v1) { // ok we find the edge
2558  if (direction>0)
2559  subdomains[i].head=t=Adj(ta);
2560  else
2561  subdomains[i].head=t=ta;
2562  if(t<triangles || t >= triangles+nbt || t->det < 0 || t->link == 0) {
2563  _error_("bad definition of SubSomain " << i);
2564  }
2565  long it = GetId(t);
2566  if (mark[it] >=0) {
2567  break;
2568  }
2569  if(i != inew)
2570  Exchange(subdomains[i],subdomains[inew]);
2571  inew++;
2572  Triangle *tt=t;
2573  long kkk=0;
2574  do
2575  {
2576  kkk++;
2577  if (mark[GetId(tt)]>=0){
2578  _error_("mark[GetId(tt)]>=0");
2579  }
2580  mark[GetId(tt)]=i;
2581  tt=tt->link;
2582  } while (tt!=t);
2583  break;
2584  }
2585  ta = Previous(Adj(ta));
2586  if(t == (Triangle *) ta) {
2587  _error_("bad definition of SubSomain " << i);
2588  }
2589  }
2590  }
2591 
2592  if (inew < nbsubdomains) {
2593  if (verbose>5) _printf_("WARNING: " << nbsubdomains-inew << " SubDomains are being removed\n");
2594  nbsubdomains=inew;}
2595 
2596  for (it=0;it<nbt;it++)
2597  if ( mark[it] ==-1 )
2598  nbtout++,triangles[it].link =0;
2599  delete [] GeomEdgetoEdge;
2600  delete [] mark;
2601 
2602  }
2603  nbtout=0;
2604  for (it=0;it<nbt;it++)
2605  if(!triangles[it].link) nbtout++;
2606  }

◆ TriangleReferenceList()

long bamg::Mesh::TriangleReferenceList ( long *  reft) const

Definition at line 4046 of file Mesh.cpp.

4046  {/*{{{*/
4047  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
4048 
4049  Triangle *t0,*t;
4050  long k=0, num;
4051 
4052  //initialize all triangles as -1 (outside)
4053  for (int it=0;it<nbt;it++) reft[it]=-1;
4054 
4055  //loop over all subdomains
4056  for (int i=0;i<nbsubdomains;i++){
4057 
4058  //first triangle of the subdomain i
4059  t=t0=subdomains[i].head;
4060 
4061  //check that the subdomain is not empty
4062  if (!t0){ _error_("At least one subdomain is empty");}
4063 
4064  //loop
4065  do{
4066  k++;
4067 
4068  //get current triangle number
4069  num = GetId(t);
4070 
4071  //check that num is in [0 nbt[
4072  _assert_(num>=0 && num<nbt);
4073 
4074  //reft of this triangle is the subdomain number
4075  reft[num]=i;
4076 
4077  } while (t0 != (t=t->link));
4078  //stop when all triangles of subdomains have been tagged
4079 
4080  }
4081  return k;
4082  }

◆ TriangleIntNumbering()

void bamg::Mesh::TriangleIntNumbering ( long *  renumbering)

Definition at line 4036 of file Mesh.cpp.

4036  {/*{{{*/
4037 
4038  long num=0;
4039  for (int i=0;i<nbt;i++){
4040  if (triangles[i].det>0) renumbering[i]=num++;
4041  else renumbering[i]=-1;
4042  }
4043  return;
4044  }

◆ CrackMesh()

void bamg::Mesh::CrackMesh ( BamgOpts bamgopts)

Definition at line 2148 of file Mesh.cpp.

2148  {/*{{{*/
2149  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/CrackMesh)*/
2150 
2151  /*Intermediary*/
2152  int i,j,k,num,count;
2153  int i1,i2;
2154  int j1,j2;
2155  int verbose=0;
2156 
2157  /*Options*/
2158  if(bamgopts) verbose=bamgopts->verbose;
2159 
2160  // computed the number of cracked edge
2161  for (k=i=0;i<nbe;i++){
2162  if(edges[i].GeomEdgeHook->Cracked()) k++;
2163  }
2164 
2165  //Return if no edge is cracked
2166  if(k==0) return;
2167  if (verbose>4) _printf_(" number of Cracked Edges = " << k << "\n");
2168 
2169  //Initialize Cracked edge
2170  NbCrackedEdges=k;
2171  CrackedEdges=new CrackedEdge[k];
2172 
2173  //Compute number of Cracked Vertices
2174  k=0;
2176 
2177  int* splitvertex=new int[nbv];
2178  for (i=0;i<nbv;i++) splitvertex[i]=0;
2179 
2180  for (i=0;i<nbe;i++){
2181  if(edges[i].GeomEdgeHook->Cracked()){
2182 
2183  //Fill edges fields of CrackedEdges
2184  CrackedEdges[k ].E =edges[i].GeomEdgeHook;
2185  CrackedEdges[k++].e1=&edges[i];
2186 
2187  //Get number of the two vertices on the edge
2188  i1=GetId(edges[i][0]);
2189  i2=GetId(edges[i][1]);
2190  _assert_(i1>=0 && i1<nbv && i2>=0 && i2<nbv);
2191  splitvertex[i1]++;
2192  splitvertex[i2]++;
2193 
2194  //If the vertex has already been flagged once, it is a cracked vertex (tip otherwise)
2195  if (splitvertex[i1]==2) NbCrackedVertices++;
2196  if (splitvertex[i2]==2) NbCrackedVertices++;
2197 
2198  //The vertex cannot be marked more than twice
2199  if (splitvertex[i1]==3 || splitvertex[i2]==3){
2200  delete [] splitvertex;
2201  _error_("Crossing rifts not supported yet");
2202  }
2203  }
2204  }
2206 
2207  //Add new vertices
2208  if (verbose>4) _printf_(" number of Cracked Vertices = " << NbCrackedVertices << "\n");
2209  if (NbCrackedVertices){
2210  CrackedVertices=xNew<long>(2*NbCrackedVertices);
2211  num=0;
2212  for (i=0;i<nbv;i++){
2213  if (splitvertex[i]==2){
2214  CrackedVertices[num*2+0]=i; //index of first vertex
2215  CrackedVertices[num*2+1]=nbv+num;//index of new vertex
2216  num++;
2217  }
2218  }
2220  }
2221  delete [] splitvertex;
2222 
2223  //Now, find the triangles that hold a cracked edge
2225 
2226  long* Edgeflags=new long[NbCrackedEdges];
2227  for(i=0;i<NbCrackedEdges;i++) Edgeflags[i]=0;
2228 
2229  for(i=0;i<NbCrackedEdges;i++){
2230  //Get the numbers of the 2 vertices of the crren cracked edge
2231  i1=GetId((*CrackedEdges[i].e1)[0]);
2232  i2=GetId((*CrackedEdges[i].e1)[1]);
2233 
2234  //find a triangle holding the vertex i1 (first vertex of the ith cracked edge)
2235  Triangle* tbegin=vertices[i1].t;
2236  k=vertices[i1].IndexInTriangle;//local number of i in triangle tbegin
2237  _assert_(GetId((*tbegin)[k])==GetId(vertices[i1]));
2238 
2239  //Now, we are going to go through the adjacent triangle that hold i1 till
2240  //we find one that has the cracked edge
2241  AdjacentTriangle ta(tbegin,EdgesVertexTriangle[k][0]);
2242  count=0;
2243  do {
2244  for(j=0;j<3;j++){
2245  //Find the position of i1 in the triangle index
2246  if (GetId((*ta.t)[j])==i1){
2247  j1=j;
2248  break;
2249  }
2250  }
2251  for(j=0;j<3;j++){
2252  //Check wether i2 is also in the triangle index
2253  if (GetId((*ta.t)[j])==i2){
2254  j2=j;
2255  //Invert j1 and j2 if necessary
2256  if ((j1+1)%3==j2){
2257  int j3=j1;
2258  j1=j2;
2259  j2=j3;
2260  }
2261  if (Edgeflags[i]==0){
2262  //first element
2263  CrackedEdges[i].a=ta.t;
2264  CrackedEdges[i].length=Norme2((*ta.t)[j1].r-(*ta.t)[j2].r);
2265  CrackedEdges[i].normal=Orthogonal((*ta.t)[j1].r-(*ta.t)[j2].r);
2266  }
2267  else{
2268  //Second element -> to renumber
2269  CrackedEdges[i].b=ta.t;
2270  CrackedEdges[i].length=Norme2((*ta.t)[j1].r-(*ta.t)[j2].r);
2271  CrackedEdges[i].normal=Orthogonal((*ta.t)[j1].r-(*ta.t)[j2].r);
2272  }
2273  Edgeflags[i]++;
2274  break;
2275  }
2276  }
2277  //_printf_(element_renu[GetId(ta.t)] << " -> " << GetId((*ta.t)[0])+1 << " " << GetId((*ta.t)[1])+1 << " " << GetId((*ta.t)[2])+1 << ", edge [" << i1 << "->" << j1 << " " << i2 << "->" << j2 << "]\n");
2278  ta = Next(ta).Adj();
2279  if (count++>50) _error_("Maximum number of iteration exceeded");
2280  }while ((tbegin != ta));
2281  }
2282 
2283  //Check EdgeFlag
2284  for(i=0;i<NbCrackedEdges;i++){
2285  if (Edgeflags[i]!=2){
2286  _error_("A problem occured: at least one crack edge (number " << i+1 << ") does not belong to 2 elements");
2287  }
2288  }
2289  delete [] Edgeflags;
2290 
2291  //Reset BamgVertex to On
2292  SetVertexFieldOn();
2293 
2294  }

◆ SmoothMetric()

void bamg::Mesh::SmoothMetric ( BamgOpts bamgopts,
double  raisonmax 
)

Definition at line 3785 of file Mesh.cpp.

3785  { /*{{{*/
3786  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
3787 
3788  /*Intermediaries*/
3789  int verbose=0;
3790 
3791  /*Get options*/
3792  if(bamgopts) verbose=bamgopts->verbose;
3793 
3794  if(raisonmax<1.1) return;
3795  if(verbose > 1) _printf_(" Mesh::SmoothMetric raisonmax = " << raisonmax << "\n");
3797  long i,j,kch,kk,ip;
3798  long *first_np_or_next_t0 = new long[nbv];
3799  long *first_np_or_next_t1 = new long[nbv];
3800  long Head0 =0,Head1=-1;
3801  double logseuil= log(raisonmax);
3802 
3803  for(i=0;i<nbv-1;i++)
3804  first_np_or_next_t0[i]=i+1;
3805  first_np_or_next_t0[nbv-1]=-1;// end;
3806  for(i=0;i<nbv;i++)
3807  first_np_or_next_t1[i]=-1;
3808  kk=0;
3809  while(Head0>=0&& kk++<100){
3810  kch=0;
3811  for(i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
3812  // pour tous les triangles autour du sommet s
3813  Triangle* t= vertices[i].t;
3814  if (!t){
3815  _error_("!t");
3816  }
3817  BamgVertex & vi = vertices[i];
3818  AdjacentTriangle ta(t,EdgesVertexTriangle[vertices[i].IndexInTriangle][0]);
3819  BamgVertex *pvj0 = ta.EdgeVertex(0);
3820  while (1) {
3821  ta=Previous(Adj(ta));
3822  if (vertices+i != ta.EdgeVertex(1)){
3823  _error_("vertices+i != ta.EdgeVertex(1)");
3824  }
3825  BamgVertex *pvj = (ta.EdgeVertex(0));
3826  BamgVertex & vj = *pvj;
3827  if(pvj){
3828  j= &vj-vertices;
3829  if (j<0 || j >= nbv){
3830  _error_("j<0 || j >= nbv");
3831  }
3832  R2 Aij = (R2) vj - (R2) vi;
3833  double ll = Norme2(Aij);
3834  if (0) {
3835  double hi = ll/vi.m.Length(Aij.x,Aij.y);
3836  double hj = ll/vj.m.Length(Aij.x,Aij.y);
3837  if(hi < hj)
3838  {
3839  double dh=(hj-hi)/ll;
3840  if (dh>logseuil) {
3841  vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
3842  if(first_np_or_next_t1[j]<0)
3843  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
3844  }
3845  }
3846  }
3847  else
3848  {
3849  double li = vi.m.Length(Aij.x,Aij.y);
3850  if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
3851  if(first_np_or_next_t1[j]<0) // if the metrix change
3852  kch++,first_np_or_next_t1[j]=Head1,Head1=j;
3853  }
3854  }
3855  if ( &vj == pvj0 ) break;
3856  }
3857  }
3858  Head0 = Head1;
3859  Head1 = -1;
3860  Exchange(first_np_or_next_t0,first_np_or_next_t1);
3861  }
3862  if(verbose>2) _printf_(" number of iterations = " << kch << "\n");
3863  delete [] first_np_or_next_t0;
3864  delete [] first_np_or_next_t1;
3865  }

◆ BoundAnisotropy()

void bamg::Mesh::BoundAnisotropy ( BamgOpts bamgopts,
double  anisomax,
double  hminaniso = 1e-100 
)

Definition at line 1155 of file Mesh.cpp.

1155  {/*{{{*/
1156  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/BoundAnisotropy)*/
1157 
1158  int verbose=0;
1159  double lminaniso = 1./ (Max(hminaniso*hminaniso,1e-100));
1160 
1161  //Get options
1162  if(bamgopts) verbose=bamgopts->verbose;
1163 
1164  //display info
1165  if (verbose > 1) _printf_(" BoundAnisotropy by " << anisomax << "\n");
1166 
1167  double h1=1.e30,h2=1e-30;
1168  double coef = 1./(anisomax*anisomax);
1169  double hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0;
1170 
1171  //loop over all vertices
1172  for (int i=0;i<nbv;i++){
1173  EigenMetric Vp(vertices[i]);
1174  double lmax=Vp.lmax();
1175  Vp*=Min(lminaniso,lmax)/lmax;
1176  Vp.BoundAniso2(coef);
1177  vertices[i].m = Vp;
1178 
1179  //info to be displayed
1180  if (verbose>2){
1181  h1 =Min(h1,Vp.lmin());
1182  h2 =Max(h2,Vp.lmax());
1183  hn1=Min(hn1,Vp.lmin());
1184  hn2=Max(hn2,Vp.lmax());
1185  rx =Max(rx,Vp.Aniso2());
1186  rnx= Max(rnx,Vp.Aniso2());
1187  }
1188  }
1189 
1190  //display info
1191  if (verbose>2){
1192  _printf_(" input: Hmin = " << pow(h2,-0.5) << ", Hmax = " << pow(h1,-0.5) << ", factor of anisotropy max = " << pow(rx,0.5) << "\n");
1193  _printf_(" output: Hmin = " << pow(hn2,-0.5) << ", Hmax = " << pow(hn1,-0.5)<< ", factor of anisotropy max = " <<pow(rnx,0.5) << "\n");
1194  }
1195  }

◆ MakeGeomEdgeToEdge()

Edge ** bamg::Mesh::MakeGeomEdgeToEdge ( )

Definition at line 2899 of file Mesh.cpp.

2899  {/*{{{*/
2900  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeGeomEdgeToEdge)*/
2901 
2902  if (!Gh.nbe){
2903  _error_("!Gh.nbe");
2904  }
2905  Edge **e= new Edge* [Gh.nbe];
2906 
2907  long i;
2908  for ( i=0;i<Gh.nbe ; i++)
2909  e[i]=NULL;
2910  for ( i=0;i<nbe ; i++)
2911  {
2912  Edge * ei = edges+i;
2913  GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
2914  e[Gh.GetId(GeomEdgeHook)] = ei;
2915  }
2916  for ( i=0;i<nbe ; i++)
2917  for (int ii=0;ii<2;ii++) {
2918  Edge * ei = edges+i;
2919  GeomEdge *GeomEdgeHook = ei->GeomEdgeHook;
2920  int j= ii;
2921  while (!(*GeomEdgeHook)[j].Required()) {
2922  Adj(GeomEdgeHook,j); // next geom edge
2923  j=1-j;
2924  if (e[Gh.GetId(GeomEdgeHook)]) break; // optimisation
2925  e[Gh.GetId(GeomEdgeHook)] = ei;
2926  }
2927  }
2928 
2929  int kk=0;
2930  for ( i=0;i<Gh.nbe ; i++){
2931  if (!e[i]){
2932  kk++;
2933  if(kk<10) _printf_("BUG: the geometrical edge " << i << " is on no edge curve\n");
2934  }
2935  }
2936  if(kk){
2937  delete [] e;
2938  _error_("See above");
2939  }
2940 
2941  return e;
2942  }

◆ SplitInternalEdgeWithBorderVertices()

long bamg::Mesh::SplitInternalEdgeWithBorderVertices ( )

Definition at line 3867 of file Mesh.cpp.

3867  {/*{{{*/
3868  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SplitInternalEdgeWithBorderVertices)*/
3869 
3870  long NbSplitEdge=0;
3871  SetVertexFieldOn();
3872  long it;
3873  long nbvold=nbv;
3874  long int verbose=2;
3875  for (it=0;it<nbt;it++){
3876  Triangle &t=triangles[it];
3877  if (t.link)
3878  for (int j=0;j<3;j++)
3879  if(!t.Locked(j) && !t.Hidden(j)){
3880 
3881  Triangle *ptt = t.TriangleAdj(j);
3882  Triangle &tt = *ptt;
3883 
3884  if( ptt && tt.link && it < GetId(tt))
3885  { // an internal edge
3886  BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
3887  BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
3888  if (v0.GeomEdgeHook && v1.GeomEdgeHook){
3889  R2 P= ((R2) v0 + (R2) v1)*0.5;
3890  if ( nbv<maxnbv) {
3891  vertices[nbv].r = P;
3892  vertices[nbv++].m = Metric(0.5,v0.m,0.5,v1.m);
3894  }
3895  NbSplitEdge++;
3896  }
3897  }
3898  }
3899  }
3901  if (nbvold!=nbv){
3902  long iv = nbvold;
3903  long NbSwap = 0;
3904  long long det3[3];
3905  for (int i=nbvold;i<nbv;i++) {// for all the new point
3906  BamgVertex & vi = vertices[i];
3907  vi.i = R2ToI2(vi.r);
3908  vi.r = I2ToR2(vi.i);
3909 
3910  // a good new point
3911  vi.ReferenceNumber=0;
3912  Triangle *tcvi = TriangleFindFromCoord(vi.i,det3);
3913  if (tcvi && !tcvi->link) {
3914  _printf_("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
3915  }
3916 
3917  quadtree->Add(vi);
3918  if (!tcvi || tcvi->det<0){// internal
3919  _error_("!tcvi || tcvi->det < 0");
3920  }
3921  AddVertex(vi,tcvi,det3);
3922  NbSwap += vi.Optim(1);
3923  iv++;
3924  }
3925  if (verbose>3) {
3926  _printf_(" number of points: " << iv << "\n");
3927  _printf_(" number of swap to split internal edges with border vertices: " << NbSwap << "\n");
3928  nbv = iv;
3929  }
3930  }
3931  if (NbSplitEdge>nbv-nbvold) _printf_("WARNING: not enough vertices to split all internal edges, we lost " << NbSplitEdge - ( nbv-nbvold) << " edges...\n");
3932  if (verbose>2) _printf_("SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << "\n");
3933 
3934  return NbSplitEdge;
3935  }

◆ MakeBamgQuadtree()

void bamg::Mesh::MakeBamgQuadtree ( )

Definition at line 2944 of file Mesh.cpp.

2944  { /*{{{*/
2945  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/
2946  if(!quadtree) quadtree = new BamgQuadtree(this);
2947  }

◆ MaxSubDivision()

void bamg::Mesh::MaxSubDivision ( BamgOpts bamgopts,
double  maxsubdiv 
)

Definition at line 2953 of file Mesh.cpp.

2953  {/*{{{*/
2954  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
2955 
2956  /*Intermediaries*/
2957  int verbose=0;
2958 
2959  /*Get options*/
2960  if(bamgopts) verbose=bamgopts->verbose;
2961 
2962  const double maxsubdiv2 = maxsubdiv*maxsubdiv;
2963  if(verbose>1) _printf_(" Limit the subdivision of a edges in the new mesh by " << maxsubdiv << "\n");
2964  // for all the edges
2965  // if the len of the edge is to long
2966  long it,nbchange=0;
2967  double lmax=0;
2968  for (it=0;it<nbt;it++){
2969  Triangle &t=triangles[it];
2970  for (int j=0;j<3;j++){
2971  Triangle *ptt=t.TriangleAdj(j);
2972  Triangle &tt = *ptt;
2973  if ( (!ptt || it < GetId(tt)) && ( tt.link || t.link)){
2974  BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
2975  BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
2976  R2 AB= (R2) v1-(R2) v0;
2977  Metric M = v0;
2978  double l = M(AB,AB);
2979  lmax = Max(lmax,l);
2980  if(l> maxsubdiv2){
2981  R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
2982  double lc = M(AC,AC);
2983  D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
2984  D2xD2 Rt1(Rt.inv());
2985  D2xD2 D(maxsubdiv2,0,0,lc);
2986  D2xD2 MM = Rt1*D*Rt1.t();
2987  v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
2988  nbchange++;
2989  }
2990  M = v1;
2991  l = M(AB,AB);
2992  lmax = Max(lmax,l);
2993  if(l> maxsubdiv2){
2994  R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M
2995  double lc = M(AC,AC);
2996  D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC;
2997  D2xD2 Rt1(Rt.inv());
2998  D2xD2 D(maxsubdiv2,0,0,lc);
2999  D2xD2 MM = Rt1*D*Rt1.t();
3000  v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
3001  nbchange++;
3002  }
3003  }
3004  }
3005  }
3006  if(verbose>3){
3007  _printf_(" number of metric changes = " << nbchange << ", maximum number of subdivision of a edges before change = " << pow(lmax,0.5) << "\n");
3008  }
3009  }

◆ NewPoints()

void bamg::Mesh::NewPoints ( Mesh Bh,
BamgOpts bamgopts,
int  KeepVertices = 1 
)

Definition at line 3040 of file Mesh.cpp.

3040  {/*{{{*/
3041  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
3042 
3043  int i,j,k;
3044  int verbose=0;
3045  long NbTSwap=0;
3046  long nbtold=nbt;
3047  long nbvold=nbv;
3048  long Headt=0;
3049  long next_t;
3050  long* first_np_or_next_t=new long[maxnbt];
3051  Triangle* t=NULL;
3052 
3053  /*Recover options*/
3054  if(bamgopts) verbose=bamgopts->verbose;
3055 
3056  /*First, insert old points if requested*/
3057  if(KeepVertices && (&Bh != this) && (nbv+Bh.nbv< maxnbv)){
3058  if (verbose>5) _printf_(" Inserting initial mesh points\n");
3059  bool pointsoutside = false;
3060  for(i=0;i<Bh.nbv;i++){
3061  BamgVertex &bv=Bh[i];
3062  /*Do not insert if the point is outside*/
3063  long long det3[3];
3064  Triangle* tcvj=TriangleFindFromCoord(bv.i,det3);
3065  if(tcvj->det<0 || !tcvj->link){
3066  pointsoutside = true;
3067  continue;
3068  }
3069  IssmDouble area_1=((bv.r.x -(*tcvj)(2)->r.x)*((*tcvj)(1)->r.y-(*tcvj)(2)->r.y)
3070  - (bv.r.y -(*tcvj)(2)->r.y)*((*tcvj)(1)->r.x-(*tcvj)(2)->r.x));
3071  IssmDouble area_2=(((*tcvj)(0)->r.x -(*tcvj)(2)->r.x)*(bv.r.y -(*tcvj)(2)->r.y)
3072  - ((*tcvj)(0)->r.y -(*tcvj)(2)->r.y)*(bv.r.x -(*tcvj)(2)->r.x));
3073  IssmDouble area_3 =((bv.r.x -(*tcvj)(1)->r.x)*((*tcvj)(0)->r.y-(*tcvj)(1)->r.y)
3074  - (bv.r.y -(*tcvj)(1)->r.y)*((*tcvj)(0)->r.x-(*tcvj)(1)->r.x));
3075  if(area_1<0 || area_2<0 || area_3<0){
3076  pointsoutside = true;
3077  continue;
3078  }
3079  if(!bv.GeomEdgeHook){
3080  vertices[nbv].r = bv.r;
3081  vertices[nbv].PreviousNumber = i+1;
3082  vertices[nbv++].m = bv.m;
3083  }
3084  }
3085  //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");
3086  Bh.CreateSingleVertexToTriangleConnectivity();
3087  InsertNewPoints(nbvold,NbTSwap,bamgopts);
3088  }
3089  else Bh.CreateSingleVertexToTriangleConnectivity();
3090 
3091  // generation of the list of next Triangle
3092  for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
3093  // the next traingle of i is -first_np_or_next_t[i]
3094 
3095  // Big loop (most time consuming)
3096  int iter=0;
3097  if (verbose>5) _printf_(" Big loop\n");
3098  do {
3099  /*Update variables*/
3100  iter++;
3101  nbtold=nbt;
3102  nbvold=nbv;
3103 
3104  /*We test all triangles*/
3105  i=Headt;
3106  next_t=-first_np_or_next_t[i];
3107  for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
3108 
3109  //check i
3110  if (i<0 || i>=nbt ){
3111  _error_("Index problem in NewPoints (i=" << i << " not in [0 " << nbt-1 << "])");
3112  }
3113  //change first_np_or_next_t[i]
3114  first_np_or_next_t[i] = iter;
3115 
3116  //Loop over the edges of t
3117  for(j=0;j<3;j++){
3118  AdjacentTriangle tj(t,j);
3119  BamgVertex &vA = *tj.EdgeVertex(0);
3120  BamgVertex &vB = *tj.EdgeVertex(1);
3121 
3122  //if t is a boundary triangle, or tj locked, continue
3123  if (!t->link) continue;
3124  if (t->det <0) continue;
3125  if (t->Locked(j)) continue;
3126 
3127  AdjacentTriangle tadjj = t->Adj(j);
3128  Triangle* ta=tadjj;
3129 
3130  //if the adjacent triangle is a boundary triangle, continue
3131  if (ta->det<0) continue;
3132 
3133  R2 A=vA;
3134  R2 B=vB;
3135  k=GetId(ta);
3136 
3137  //if this edge has already been done, go to next edge of triangle
3138  if(first_np_or_next_t[k]==iter) continue;
3139 
3140  lIntTria.SplitEdge(Bh,A,B);
3142  } // end loop for each edge
3143  }// for triangle
3144 
3145  if (!InsertNewPoints(nbvold,NbTSwap,bamgopts)) break;
3146  for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
3147  Headt = nbt; // empty list
3148 
3149  // for all the triangle containing the vertex i
3150  for (i=nbvold;i<nbv;i++){
3151  BamgVertex* s = vertices + i;
3152  AdjacentTriangle ta(s->t, EdgesVertexTriangle[s->IndexInTriangle][1]);
3153  Triangle* tbegin= (Triangle*) ta;
3154  long kt;
3155  do {
3156  kt = GetId((Triangle*) ta);
3157  if (first_np_or_next_t[kt]>0){
3158  first_np_or_next_t[kt]=-Headt;
3159  Headt=kt;
3160  }
3161  if (ta.EdgeVertex(0)!=s){
3162  _error_("ta.EdgeVertex(0)!=s");
3163  }
3164  ta = Next(Adj(ta));
3165  } while ( (tbegin != (Triangle*) ta));
3166  }
3167 
3168  }while(nbv!=nbvold);
3169  delete [] first_np_or_next_t;
3170 
3171  long NbSwapf =0;
3172  for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
3173  }/*}}}*/

◆ InsertNewPoints()

long bamg::Mesh::InsertNewPoints ( long  nbvold,
long &  NbTSwap,
BamgOpts bamgopts 
)

Definition at line 2815 of file Mesh.cpp.

2815  {/*{{{*/
2816  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/
2817 
2818  int verbose=0;
2819  double seuil= 1.414/2.;// for two close point
2820  long i;
2821  long NbSwap=0;
2822  long long det3[3];
2823 
2824  /*Get options*/
2825  if(bamgopts) verbose=bamgopts->verbose;
2826 
2827  //number of new points
2828  const long nbvnew=nbv-nbvold;
2829 
2830  //display info if required
2831  if (verbose>5) _printf_(" Try to Insert " << nbvnew << " new points\n");
2832 
2833  //return if no new points
2834  if (!nbvnew) return 0;
2835 
2836  /*construction of a random order*/
2837  const long PrimeNumber= BigPrimeNumber(nbv) ;
2838  long k3 = this->RandomNumber(nbvnew);
2839  //loop over the new points
2840  for (int is3=0; is3<nbvnew; is3++){
2841  long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
2842  long i=nbvold+is3;
2843  orderedvertices[i]= vertices + j;
2845  }
2846 
2847  // for all the new point
2848  long iv=nbvold;
2849  for(i=nbvold;i<nbv;i++){
2850  BamgVertex &vi=*orderedvertices[i];
2851  vi.i=R2ToI2(vi.r);
2852  vi.r=I2ToR2(vi.i);
2853  double hx,hy;
2854  vi.m.Box(hx,hy);
2855  int hi=(int) (hx*coefIcoor),hj=(int) (hy*coefIcoor);
2856  if(!quadtree->TooClose(&vi,seuil,hi,hj)){
2857  // a good new point
2858  BamgVertex &vj = vertices[iv];
2859  long j=vj.ReferenceNumber;
2860  if (&vj!=orderedvertices[j]){
2861  _error_("&vj!= orderedvertices[j]");
2862  }
2863  if(i!=j){
2864  Exchange(vi,vj);
2866  }
2867  vj.ReferenceNumber=0;
2868  Triangle *tcvj=TriangleFindFromCoord(vj.i,det3);
2869  if (tcvj && !tcvj->link){
2870  _printf_("While trying to add the following point:\n");
2871  vj.Echo();
2872  _printf_("BAMG determined that it was inside the following triangle, which probably lies outside of the geometric domain\n");
2873  tcvj->Echo();
2874  _error_("problem inserting point in InsertNewPoints (tcvj=" << tcvj << " and tcvj->link=" << tcvj->link << ")");
2875  }
2876  quadtree->Add(vj);
2877  AddVertex(vj,tcvj,det3);
2878  NbSwap += vj.Optim(1);
2879  iv++;
2880  }
2881  else{
2882  vi.PreviousNumber = 0;
2883  }
2884  }
2885  if (verbose>3) {
2886  _printf_(" number of new points: " << iv << "\n");
2887  _printf_(" number of to close (?) points: " << nbv-iv << "\n");
2888  _printf_(" number of swap after: " << NbSwap << "\n");
2889  }
2890  nbv = iv;
2891 
2892  for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1);
2893  if (verbose>3) _printf_(" NbSwap=" << NbSwap << "\n");
2894 
2895  NbTSwap += NbSwap ;
2896  return nbv-nbvold;
2897  }

◆ TrianglesRenumberBySubDomain()

void bamg::Mesh::TrianglesRenumberBySubDomain ( bool  justcompress = false)

Definition at line 3611 of file Mesh.cpp.

3611  {/*{{{*/
3612  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
3613 
3614  long *renu= new long[nbt];
3615  Triangle *t0,*t,*te=triangles+nbt;
3616  long k=0,it,i,j;
3617 
3618  for ( it=0;it<nbt;it++)
3619  renu[it]=-1; // outside triangle
3620  for ( i=0;i<nbsubdomains;i++)
3621  {
3622  t=t0=subdomains[i].head;
3623  if (!t0){ // not empty sub domain
3624  _error_("!t0");
3625  }
3626  do {
3627  long kt = GetId(t);
3628  if (kt<0 || kt >= nbt ){
3629  _error_("kt<0 || kt >= nbt");
3630  }
3631  if (renu[kt]!=-1){
3632  _error_("renu[kt]!=-1");
3633  }
3634  renu[kt]=k++;
3635  }
3636  while (t0 != (t=t->link));
3637  }
3638  // take is same numbering if possible
3639  if(justcompress)
3640  for ( k=0,it=0;it<nbt;it++)
3641  if(renu[it] >=0 )
3642  renu[it]=k++;
3643 
3644  // put the outside triangles at the end
3645  for ( it=0;it<nbt;it++){
3646  if (renu[it]==-1) renu[it]=k++;
3647  }
3648  if (k != nbt){
3649  _error_("k != nbt");
3650  }
3651  // do the change on all the pointeur
3652  for ( it=0;it<nbt;it++)
3653  triangles[it].Renumbering(triangles,te,renu);
3654 
3655  for ( i=0;i<nbsubdomains;i++)
3656  subdomains[i].head=triangles+renu[GetId(subdomains[i].head)];
3657 
3658  // move the Triangles without a copy of the array
3659  // be carefull not trivial code
3660  for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
3661  if (renu[it] >= 0) // a new sub cycle
3662  {
3663  i=it;
3664  Triangle ti=triangles[i],tj;
3665  while ( (j=renu[i]) >= 0)
3666  { // i is old, and j is new
3667  renu[i] = -1; // mark
3668  tj = triangles[j]; // save new
3669  triangles[j]= ti; // new <- old
3670  i=j; // next
3671  ti = tj;
3672  }
3673  }
3674  delete [] renu;
3675 
3676  }

◆ SmoothingVertex()

void bamg::Mesh::SmoothingVertex ( BamgOpts bamgopts,
int  nbiter = 3,
double  omega = 0.3 
)

Definition at line 3740 of file Mesh.cpp.

3740  { /*{{{*/
3741  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
3742 
3743  /*Intermediaries*/
3744  int verbose=0;
3745 
3746  /*Get options*/
3747  if(bamgopts) verbose=bamgopts->verbose;
3748 
3749  // if quatree exist remove it end reconstruct
3750  if (quadtree) delete quadtree;
3751  quadtree=0;
3753  Triangle vide; // a triangle to mark the boundary vertex
3754  Triangle ** tstart= new Triangle* [nbv];
3755  long i,j,k;
3756  // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
3757  if ( this == & BTh)
3758  for ( i=0;i<nbv;i++)
3759  tstart[i]=vertices[i].t;
3760  else
3761  for ( i=0;i<nbv;i++)
3762  tstart[i]=0;
3763  for ( j=0;j<NbVerticesOnGeomVertex;j++ )
3764  tstart[ GetId(VerticesOnGeomVertex[j].meshvertex)]=&vide;
3765  for ( k=0;k<NbVerticesOnGeomEdge;k++ )
3766  tstart[ GetId(VerticesOnGeomEdge[k].meshvertex)]=&vide;
3767  if(verbose>2) _printf_(" SmoothingVertex: nb Iteration = " << nbiter << ", Omega=" << omega << "\n");
3768  for (k=0;k<nbiter;k++)
3769  {
3770  long i,NbSwap =0;
3771  double delta =0;
3772  for ( i=0;i<nbv;i++)
3773  if (tstart[i] != &vide) // not a boundary vertex
3774  delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
3775  for ( i=0;i<nbv;i++)
3776  if (tstart[i] != &vide) // not a boundary vertex
3777  NbSwap += vertices[i].Optim(1);
3778  if (verbose>3) _printf_(" move max = " << pow(delta,0.5) << ", iteration = " << k << ", nb of swap = " << NbSwap << "\n");
3779  }
3780 
3781  delete [] tstart;
3782  if (quadtree) quadtree= new BamgQuadtree(this);
3783  }

◆ MetricAt()

Metric bamg::Mesh::MetricAt ( const R2 A)

Definition at line 3011 of file Mesh.cpp.

3011  { /*{{{*/
3012  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/
3013 
3014  I2 a = R2ToI2(A);
3015  long long deta[3];
3016  Triangle * t =TriangleFindFromCoord(a,deta);
3017  if (t->det <0) { // outside
3018  double ba,bb;
3019  AdjacentTriangle edge= CloseBoundaryEdge(a,t,ba,bb) ;
3020  return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
3021  else { // inside
3022  double aa[3];
3023  double s = deta[0]+deta[1]+deta[2];
3024  aa[0]=deta[0]/s;
3025  aa[1]=deta[1]/s;
3026  aa[2]=deta[2]/s;
3027  return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
3028  }
3029  }

◆ ProjectOnCurve()

GeomEdge * bamg::Mesh::ProjectOnCurve ( Edge AB,
BamgVertex A,
BamgVertex B,
double  theta,
BamgVertex R,
VertexOnEdge BR,
VertexOnGeom GR 
)

Definition at line 3174 of file Mesh.cpp.

3175  {
3176  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/ProjectOnCurve)*/
3177 
3178  void *pA=0,*pB=0;
3179  double tA=0,tB=0;
3180  R2 A=vA,B=vB;
3181  BamgVertex * pvA=&vA, * pvB=&vB;
3182  if (vA.IndexInTriangle == IsVertexOnVertex){
3183  pA=vA.BackgroundVertexHook;
3184  }
3185  else if (vA.IndexInTriangle == IsVertexOnEdge){
3186  pA=vA.BackgroundEdgeHook->be;
3187  tA=vA.BackgroundEdgeHook->abcisse;
3188  }
3189  else {
3190  _error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vA) << " forget call to SetVertexFieldOnBTh");
3191  }
3192 
3193  if (vB.IndexInTriangle == IsVertexOnVertex){
3194  pB=vB.BackgroundVertexHook;
3195  }
3196  else if(vB.IndexInTriangle == IsVertexOnEdge){
3197  pB=vB.BackgroundEdgeHook->be;
3198  tB=vB.BackgroundEdgeHook->abcisse;
3199  }
3200  else {
3201  _error_("ProjectOnCurve On BamgVertex " << BTh.GetId(vB) << " forget call to SetVertexFieldOnBTh");
3202  }
3203  Edge * e = &BhAB;
3204  if (!pA || !pB || !e){
3205  _error_("!pA || !pB || !e");
3206  }
3207  // be carefull the back ground edge e is on same geom edge
3208  // of the initiale edge def by the 2 vertex A B;
3209  //check Is a background Mesh;
3210  if (e<BTh.edges || e>=BTh.edges+BTh.nbe){
3211  _error_("e<BTh.edges || e>=BTh.edges+BTh.nbe");
3212  }
3213  // walk on BTh edge
3214  //not finish ProjectOnCurve with BackGround Mesh);
3215  // 1 first find a back ground edge contening the vertex A
3216  // 2 walk n back gound boundary to find the final vertex B
3217 
3218  if( vA.IndexInTriangle == IsVertexOnEdge)
3219  { // find the start edge
3220  e = vA.BackgroundEdgeHook->be;
3221 
3222  }
3223  else if (vB.IndexInTriangle == IsVertexOnEdge)
3224  {
3225  theta = 1-theta;
3226  Exchange(tA,tB);
3227  Exchange(pA,pB);
3228  Exchange(pvA,pvB);
3229  Exchange(A,B);
3230  e = vB.BackgroundEdgeHook->be;
3231 
3232  }
3233  else{ // do the search by walking
3234  _error_("case not supported yet");
3235  }
3236 
3237  // find the direction of walking with direction of edge and pA,PB;
3238  R2 AB=B-A;
3239 
3240  double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
3241  int kkk=0;
3242  int direction = (cosE01AB>0) ? 1 : 0;
3243 
3244  // double l=0; // length of the edge AB
3245  double abscisse = -1;
3246 
3247  for (int step=0;step<2;step++){
3248  // 2 times algo:
3249  // 1 for computing the length l
3250  // 2 for find the vertex
3251  int iii;
3252  BamgVertex *v0=pvA,*v1;
3253  Edge *neee,*eee;
3254  double lg =0; // length of the curve
3255  double te0;
3256  // we suppose take the curve's abcisse
3257  for ( eee=e,iii=direction,te0=tA;
3258  eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
3259  neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
3260 
3261  kkk=kkk+1;
3262  _assert_(kkk<100);
3263  _assert_(eee);
3264  double lg0 = lg;
3265  double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
3266  lg += dp;
3267  if (step && abscisse <= lg) { // ok we find the geom edge
3268  double sss = (abscisse-lg0)/dp;
3269  double thetab = te0*(1-sss)+ sss*iii;
3270  _assert_(thetab>=0 && thetab<=1);
3271  BR = VertexOnEdge(&R,eee,thetab);
3272  return Gh.ProjectOnCurve(*eee,thetab,R,GR);
3273  }
3274  }
3275  // we find the end
3276  if (v1 != pvB){
3277  if (( void*) v1 == pB)
3278  tB = iii;
3279 
3280  double lg0 = lg;
3281  _assert_(eee);
3282  v1 = pvB;
3283  double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
3284  lg += dp;
3285  abscisse = lg*theta;
3286  if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end
3287  { // ok we find the geom edge
3288  double sss = (abscisse-lg0)/dp;
3289  double thetab = te0*(1-sss)+ sss*tB;
3290  _assert_(thetab>=0 && thetab<=1);
3291  BR = VertexOnEdge(&R,eee,thetab);
3292  return Gh.ProjectOnCurve(*eee,thetab,R,GR);
3293  }
3294  }
3295  abscisse = lg*theta;
3296 
3297  }
3298  _error_("Big bug...");
3299  return 0; // just for the compiler
3300  }

◆ GetId() [1/6]

long bamg::Mesh::GetId ( const Triangle t) const

Definition at line 2608 of file Mesh.cpp.

2608  { /*{{{*/
2609  return &t - triangles;
2610  }

◆ GetId() [2/6]

long bamg::Mesh::GetId ( const Triangle t) const

Definition at line 2612 of file Mesh.cpp.

2612  { /*{{{*/
2613  return t - triangles;
2614  }

◆ GetId() [3/6]

long bamg::Mesh::GetId ( const BamgVertex t) const

Definition at line 2616 of file Mesh.cpp.

2616  { /*{{{*/
2617  return &t - vertices;
2618  }

◆ GetId() [4/6]

long bamg::Mesh::GetId ( const BamgVertex t) const

Definition at line 2620 of file Mesh.cpp.

2620  { /*{{{*/
2621  return t - vertices;
2622  }

◆ GetId() [5/6]

long bamg::Mesh::GetId ( const Edge t) const

Definition at line 2624 of file Mesh.cpp.

2624  { /*{{{*/
2625  return &t - edges;
2626  }

◆ GetId() [6/6]

long bamg::Mesh::GetId ( const Edge t) const

Definition at line 2628 of file Mesh.cpp.

2628  { /*{{{*/
2629  return t - edges;
2630  }

◆ NearestVertex()

BamgVertex * bamg::Mesh::NearestVertex ( int  i,
int  j 
)

Definition at line 3035 of file Mesh.cpp.

3035  {/*{{{*/
3036  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NearestVertex)*/
3037  return quadtree->NearestVertex(i,j);
3038  }

◆ TriangleFindFromCoord()

Triangle * bamg::Mesh::TriangleFindFromCoord ( const I2 B,
long long  det3[3],
Triangle tstart = 0 
)

Definition at line 3945 of file Mesh.cpp.

3945  {/*{{{*/
3946  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
3947 
3948  Triangle * t=0;
3949  int j,jp,jn,jj;
3950  int counter;
3951 
3952  /*Get starting triangle. Take tsart if provided*/
3953  if (tstart) t=tstart;
3954 
3955  /*Else find the closest Triangle using the quadtree*/
3956  else {
3957 
3958  /*Check that the quadtree does exist*/
3959  if (!quadtree) _error_("no starting triangle provided and no quadtree available");
3960 
3961  /*Call NearestVertex*/
3962  BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
3963 
3964  /*Check output (Vertex a)*/
3965  if (!a) _error_("problem while trying to find nearest vertex from a given point. No output found");
3966  if (!a->t) _error_("no triangle is associated to vertex number " << GetId(a)+1 << " (orphan?)");
3967  _assert_(a>=vertices && a<vertices+nbv);
3968 
3969  /*Get starting triangle*/
3970  t = a->t;
3971  _assert_(t>=triangles && t<triangles+nbt);
3972  }
3973 
3974  long long detop ;
3975 
3976  /*initialize number of test triangle*/
3977  counter=0;
3978 
3979  /*The initial triangle might be outside*/
3980  while (t->det < 0){
3981 
3982  /*Get a real vertex from this triangle (k0)*/
3983  int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
3984  _assert_(k0>=0);// k0 the NULL vertex
3985  int k1=NextVertex[k0],k2=PreviousVertex[k0];
3986  det3[k0]=det(B,(*t)[k1],(*t)[k2]);
3987  det3[k1]=det3[k2]=-1;
3988  if (det3[k0] > 0) // outside B
3989  return t;
3990  t = t->TriangleAdj(OppositeEdge[k0]);
3991  counter++;
3992  _assert_(counter<2);
3993  }
3994 
3995  jj=0;
3996  detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
3997 
3998  while(t->det>0){
3999 
4000  /*Increase counter*/
4001  if (++counter>=10000) _error_("Maximum number of iteration reached (threshold = " << counter << ").");
4002 
4003  j= OppositeVertex[jj];
4004  det3[j] = detop; //det(*b,*s1,*s2);
4005  jn = NextVertex[j];
4006  jp = PreviousVertex[j];
4007  det3[jp]= det(*(*t)(j),*(*t)(jn),B);
4008  det3[jn] = t->det-det3[j] -det3[jp];
4009 
4010  // count the number k of det3 <0
4011  int k=0,ii[3];
4012  if (det3[0] < 0 ) ii[k++]=0;
4013  if (det3[1] < 0 ) ii[k++]=1;
4014  if (det3[2] < 0 ) ii[k++]=2;
4015  // 0 => ok
4016  // 1 => go in way 1
4017  // 2 => two way go in way 1 or 2 randomly
4018 
4019  if (k==0) break;
4020  if (k==2 && this->RandomNumber(1)) Exchange(ii[0],ii[1]);
4021  _assert_(k<3);
4022  AdjacentTriangle t1 = t->Adj(jj=ii[0]);
4023  if ((t1.det() < 0 ) && (k == 2))
4024  t1 = t->Adj(jj=ii[1]);
4025  t=t1;
4026  j=t1;// for optimisation we now the -det[OppositeVertex[j]];
4027  detop = -det3[OppositeVertex[jj]];
4028  jj = j;
4029  }
4030 
4031  if (t->det<0) // outside triangle
4032  det3[0]=det3[1]=det3[2]=-1,det3[OppositeVertex[jj]]=detop;
4033  return t;
4034  }

◆ ReadMesh() [1/2]

void bamg::Mesh::ReadMesh ( int *  index,
double *  x,
double *  y,
int  nods,
int  nels,
BamgOpts bamgopts 
)

Definition at line 249 of file Mesh.cpp.

249  {/*{{{*/
250 
251  long i1,i2,i3;
252  long i;
253  Metric M1(1);
254  int verbose=0;
255  bool* nodeflags=NULL;
256 
257  nbv=nods;
258  maxnbv=nbv;
259  nbt=nels;
260  if(bamgopts) verbose=bamgopts->verbose;
261 
262  //Vertices
263  if (verbose) _printf_("Reading vertices (" << nbv << ")\n");
264  vertices=xNew<BamgVertex>(nbv);
265  orderedvertices=xNew<BamgVertex*>(nbv);
266  for (i=0;i<nbv;i++){
267  vertices[i].r.x=x[i];
268  vertices[i].r.y=y[i];
270  vertices[i].m=M1;
271  vertices[i].color=0;
272  }
273  maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
274 
275  //Triangles
276  if (verbose) _printf_("Reading triangles (" << nbt << ")\n");
277  triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
278  nodeflags=xNew<bool>(nbv);
279  for(i=0;i<nbv;i++) nodeflags[i]=false;
280  //other triangles will be added for each edge
281  for (i=0;i<nbt;i++){
282  Triangle & t = triangles[i];
283  i1=(long)index[i*3+0]-1; //for C indexing
284  i2=(long)index[i*3+1]-1; //for C indexing
285  i3=(long)index[i*3+2]-1; //for C indexing
286  t=Triangle(this,i1,i2,i3);
287  t.color=1;
288  nodeflags[i1]=nodeflags[i2]=nodeflags[i3]=true;
289  }
290 
291  /*Recreate geometry: */
292  if (verbose) _printf_("Building Geometry\n");
294  if (verbose) _printf_("Completing geometry\n");
295  Gh.PostRead();
296 
297  /*Check that there is no orphan*/
298  bool isorphan=false;
299  for(i=0;i<nbv;i++){
300  if(!nodeflags[i]){
301  _printf_("Vertex " << i+1 << " does not belong to any element\n");
302  isorphan=true;
303  }
304  }
305  if(isorphan) _error_("Orphan found in mesh, see ids above");
306 
307  /*Clean up*/
308  xDelete<bool>(nodeflags);
309  }

◆ ReadMesh() [2/2]

void bamg::Mesh::ReadMesh ( BamgMesh bamgmesh,
BamgOpts bamgopts 
)

Definition at line 311 of file Mesh.cpp.

311  {/*{{{*/
312 
313  int verbose=0;
314  double Hmin = HUGE_VAL; // the infinie value
315  long i1,i2,i3;
316  long i,j;
317  Metric M1(1);
318 
319  /*Check needed pointer*/
320  _assert_(bamgmesh);
321 
322  if(bamgopts) verbose=bamgopts->verbose;
323 
324  nbv=bamgmesh->VerticesSize[0];
325  maxnbv=nbv;
326  nbt=bamgmesh->TrianglesSize[0];
327 
328  //Vertices
329  if(bamgmesh->Vertices){
330  if(verbose>5) _printf_(" processing Vertices\n");
331 
332  vertices=xNew<BamgVertex>(nbv);
333  orderedvertices=xNew<BamgVertex*>(nbv);
334 
335  for (i=0;i<nbv;i++){
336  vertices[i].r.x=bamgmesh->Vertices[i*3+0];
337  vertices[i].r.y=bamgmesh->Vertices[i*3+1];
338  vertices[i].ReferenceNumber=(long)bamgmesh->Vertices[i*3+2];
339  vertices[i].m=M1;
340  vertices[i].color=0;
341  }
342  maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
343  }
344  else{
345  if(verbose>5) _error_("no Vertices found in the initial mesh");
346  }
347 
348  //Triangles
349  if(bamgmesh->Triangles){
350  if(verbose>5) _printf_(" processing Triangles\n");
351  triangles =new Triangle[maxnbt]; //we cannot allocate only nbt triangles since
352  //other triangles will be added for each edge
353  for (i=0;i<nbt;i++){
354  Triangle &t=triangles[i];
355  i1=(long)bamgmesh->Triangles[i*4+0]-1; //for C indexing
356  i2=(long)bamgmesh->Triangles[i*4+1]-1; //for C indexing
357  i3=(long)bamgmesh->Triangles[i*4+2]-1; //for C indexing
358  t=Triangle(this,i1,i2,i3);
359  t.color=(long)bamgmesh->Triangles[i*4+3];
360  }
361  }
362  else{
363  if(verbose>5) _error_("no Triangles found in the initial mesh");
364  }
365 
366  //VerticesOnGeomEdge
367  if(bamgmesh->VerticesOnGeomEdge){
368  if(verbose>5) _printf_(" processing VerticesOnGeomEdge\n");
370  VerticesOnGeomEdge= new VertexOnGeom[NbVerticesOnGeomEdge] ;
371  for (i=0;i<NbVerticesOnGeomEdge;i++){
372  long i1,i2;
373  double s;
374  i1=(long) bamgmesh->VerticesOnGeomEdge[i*3+0]-1; //for C indexing
375  i2=(long) bamgmesh->VerticesOnGeomEdge[i*3+1]-1; //for C indexing
376  s =(double)bamgmesh->VerticesOnGeomEdge[i*3+2];
377  VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
378  }
379  }
380 
381  //VerticesOnGeomVertex
382  if(bamgmesh->VerticesOnGeomVertexSize[0]){
383  if(verbose>5) _printf_(" processing VerticesOnGeomVertex\n");
385  VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex] ;
386  for (i=0;i<NbVerticesOnGeomVertex;i++){
387  long i1,i2;
388  i1=(long)bamgmesh->VerticesOnGeomVertex[i*2+0]-1; //for C indexing
389  i2=(long)bamgmesh->VerticesOnGeomVertex[i*2+1]-1; //for C indexing
390  VerticesOnGeomVertex[i]=VertexOnGeom(vertices[i1],Gh.vertices[i2]);
391  }
392  }
393 
394  //Edges
395  if (bamgmesh->Edges){
396  int i1,i2;
397  double* len=NULL;
398 
399  if(verbose>5) _printf_(" processing Edges\n");
400  nbe=bamgmesh->EdgesSize[0];
401  edges= new Edge[nbe];
402  //initialize length of each edge (used to provided metric)
403  len= new double[nbv];
404  for(i=0;i<nbv;i++) len[i]=0;
405 
406  for (i=0;i<nbe;i++){
407  i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
408  i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
409  edges[i].ReferenceNumber=(long)bamgmesh->Edges[i*3+2];
410  edges[i].v[0]= vertices +i1;
411  edges[i].v[1]= vertices +i2;
412  edges[i].adj[0]=NULL;
413  edges[i].adj[1]=NULL;
414  R2 x12=vertices[i2].r-vertices[i1].r;
415  double l12=sqrt((x12,x12));
416 
417  //prepare metric
418  vertices[i1].color++;
419  vertices[i2].color++;
420  len[i1]+=l12;
421  len[i2]+=l12;
422  Hmin = Min(Hmin,l12);
423  }
424 
425  // definition the default of the given mesh size
426  for (i=0;i<nbv;i++){
427  if (vertices[i].color>0)
428  vertices[i].m=Metric(len[i]/(double)vertices[i].color);
429  else
430  vertices[i].m=Metric(Hmin);
431  }
432  delete [] len;
433 
434  // construction of edges[].adj
435  for (i=0;i<nbv;i++){
436  vertices[i].color=(vertices[i].color ==2) ?-1:-2;
437  }
438  for (i=0;i<nbe;i++){
439  for (j=0;j<2;j++) {
440  BamgVertex *v=edges[i].v[j];
441  long i0=v->color,j0;
442  if(i0==-1){
443  v->color=i*2+j;
444  }
445  else if (i0>=0) {// i and i0 edge are adjacent by the vertex v
446  j0 = i0%2;
447  i0 = i0/2;
448  _assert_(v==edges[i0 ].v[j0]);
449  edges[i ].adj[j ] =edges +i0;
450  edges[i0].adj[j0] =edges +i ;
451  v->color = -3;
452  }
453  }
454  }
455  }
456 
457  //EdgeOnGeomEdge
458  if(bamgmesh->EdgesOnGeomEdge){
459  if(verbose>5) _printf_(" processing EdgesOnGeomEdge\n");
460  int i1,i2,i,j;
461  i2=bamgmesh->EdgesOnGeomEdgeSize[0];
462  for (i1=0;i1<i2;i1++) {
463  i=(int)bamgmesh->EdgesOnGeomEdge[i1*2+0]-1; //C indexing
464  j=(int)bamgmesh->EdgesOnGeomEdge[i1*2+1]-1; //C indexing
465  //Check value
466  if(!(i>=0 && j>=0 && i<nbe && j<Gh.nbe)) {
467  _error_("ReadMesh error: EdgesOnGeomEdge edge provided (line " << i1+1 << ": [" << i+1 << " " << j+1 << "]) is incorrect (must be positive, [0<i<nbe=" << nbe << " 0<j<Gh.nbe=" << Gh.nbe << "]");
468  }
469  edges[i].GeomEdgeHook=Gh.edges+j;
470  }
471  }
472 
473  //SubDomain
474  if(bamgmesh->SubDomains){
475  long i3,head,direction;
476  if(verbose>5) _printf_(" processing SubDomains\n");
477  nbsubdomains=bamgmesh->SubDomainsSize[0];
478  subdomains = new SubDomain [ nbsubdomains ];
479  for (i=0;i<nbsubdomains;i++) {
480  i3 =(int)bamgmesh->SubDomains[i*4+0];
481  head=(int)bamgmesh->SubDomains[i*4+1]-1;//C indexing
482  direction=(int)bamgmesh->SubDomains[i*4+2];
483  if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
484  if (head<0 || head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
485  subdomains[i].head = triangles+head;
486  subdomains[i].direction = direction;
487  subdomains[i].ReferenceNumber = i3;
488  }
489  }
490 
491  }

◆ WriteMesh()

void bamg::Mesh::WriteMesh ( BamgMesh bamgmesh,
BamgOpts bamgopts 
)

Definition at line 493 of file Mesh.cpp.

493  {/*{{{*/
494 
495  /*Intermediary*/
496  int i,j,k,num,i1,i2;
497  long n;
498  int* head_1=NULL;
499  int* next_1=NULL;
500  int* connectivitysize_1=NULL;
501  int connectivitymax_1=0;
502  int verbose=0;
503 
504  /*Check needed pointer*/
505  _assert_(bamgmesh);
506 
507  /*Get options*/
508  if(bamgopts) verbose=bamgopts->verbose;
509 
510  /*Build reft that holds the number the subdomain number of each triangle, and the real numbering of the elements*/
511  long* reft = new long[nbt];
512  long* numt = new long[nbt];
513  long nbInT = TriangleReferenceList(reft);
514  TriangleIntNumbering(numt);
515 
516  /*Chaining algorithm used to generate connectivity tables and other outputs*/
517 
518  //Memory Allocation
519  head_1=xNew<int>(nbv);
520  next_1=xNew<int>(3*nbt);
521  connectivitysize_1=xNew<int>(nbv);
522 
523  //Initialization
524  for (i=0;i<nbv;i++) head_1[i]=-1;
525  for (i=0;i<nbv;i++) connectivitysize_1[i]=0;
526  k=0;
527  //Chains generation
528  for (i=0;i<nbt;i++) {
529  //Do not take into account outside triangles (reft<0)
530  if (reft[i]>=0){
531  for (j=0;j<3;j++){
532  int v=GetId(triangles[i][j]); //jth vertex of the ith triangle
533  if (k>3*nbt-1 || k<0) _error_("k = " << k << ", nbt = " << nbt);
534  next_1[k]=head_1[v];
535  if (v>nbv-1 || v<0) _error_("v = " << v << ", nbv = " << nbv);
536  head_1[v]=k++;
537  connectivitysize_1[v]+=1;
538  }
539  }
540  }
541  //Get maximum connectivity
542  connectivitymax_1=0;
543  for (i=0;i<nbv;i++){
544  if (connectivitysize_1[i]>connectivitymax_1) connectivitymax_1=connectivitysize_1[i];
545  }
546 
547  /*OK, now build outputs*/
548 
549  /*Vertices*/
550  if(verbose>5) _printf_(" writing Vertices\n");
551  bamgmesh->VerticesSize[0]=nbv;
552  bamgmesh->VerticesSize[1]=3;
553  if(nbv){
554  bamgmesh->Vertices=xNew<double>(3*nbv);
555  bamgmesh->PreviousNumbering=xNew<double>(nbv);
556  for (i=0;i<nbv;i++){
557  bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
558  bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
559  bamgmesh->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
560  bamgmesh->PreviousNumbering[i]=vertices[i].PreviousNumber;
561  }
562  }
563 
564  /*Edges*/
565  if(verbose>5) _printf_(" writing Edges\n");
566  bamgmesh->EdgesSize[0]=nbe;
567  bamgmesh->EdgesSize[1]=3;
568  int NumIssmSegments=0;
569  if (nbe){
570  bamgmesh->Edges=xNew<double>(3*nbe);
571  for (i=0;i<nbe;i++){
572  bamgmesh->Edges[i*3+0]=GetId(edges[i][0])+1; //back to M indexing
573  bamgmesh->Edges[i*3+1]=GetId(edges[i][1])+1; //back to M indexing
574  bamgmesh->Edges[i*3+2]=edges[i].ReferenceNumber;
575  if(edges[i].GeomEdgeHook){
576  NumIssmSegments++;
577  }
578  }
579  }
580 
581  /*Element edges*/
582  if(verbose>5) _printf_(" writing element edges\n");
583  SetOfEdges4* edge4=new SetOfEdges4(nbt*3,nbv);
584  double* elemedge=NULL;
585  elemedge=xNew<double>(3*nbt);
586  for (i=0;i<3*nbt;i++) elemedge[i]=-2.;//will become -1
587  k=0;
588  for (i=0;i<nbt;i++){
589  //Do not take into account outside triangles (reft<0)
590  if (reft[i]>=0){
591  for (j=0;j<3;j++) {
594  n =edge4->SortAndFind(i1,i2);
595  if (n==-1){
596  //first time
597  n=edge4->SortAndAdd(i1,i2);
598  elemedge[n*2+0]=double(k);
599  }
600  else{
601  //second time
602  elemedge[n*2+1]=double(k);
603  }
604  }
605  k++;
606  }
607  }
608  bamgmesh->IssmEdgesSize[0]=edge4->nb();
609  bamgmesh->IssmEdgesSize[1]=4;
610  bamgmesh->IssmEdges=xNew<double>(4*edge4->nb());
611  for (i=0;i<edge4->nb();i++){
612  /*Invert first two vertices if necessary*/
613  bool found=false;
614  for (j=0;j<3;j++){
615  if (triangles[(int)elemedge[2*i+0]](j)==vertices+edge4->i(i)){
616  if (triangles[(int)elemedge[2*i+0]]((j+1)%3)==vertices+edge4->j(i)){
617  //trigonometric direction
618  bamgmesh->IssmEdges[i*4+0]=edge4->i(i)+1;// back to M indexing
619  bamgmesh->IssmEdges[i*4+1]=edge4->j(i)+1;// back to M indexing
620  }
621  else{
622  bamgmesh->IssmEdges[i*4+0]=edge4->j(i)+1;// back to M indexing
623  bamgmesh->IssmEdges[i*4+1]=edge4->i(i)+1;// back to M indexing
624  }
625  found=true;
626  break;
627  }
628  }
629  _assert_(found);
630  bamgmesh->IssmEdges[i*4+2]=elemedge[2*i+0]+1; // back to M indexing
631  bamgmesh->IssmEdges[i*4+3]=elemedge[2*i+1]+1; // back to M indexing
632  }
633  //clean up
634  delete edge4;
635  xDelete<double>(elemedge);
636 
637  /*IssmSegments*/
638  if(verbose>5) _printf_(" writing IssmSegments\n");
639  bamgmesh->IssmSegmentsSize[0]=NumIssmSegments;
640  bamgmesh->IssmSegmentsSize[1]=4;
641  bamgmesh->IssmSegments=xNew<double>(4*NumIssmSegments);
642  num=0;
643  for (i=0;i<nbe;i++){
644  if(edges[i].GeomEdgeHook){
645  //build segment
646  int i1=GetId(edges[i][0]);
647  int i2=GetId(edges[i][1]);
648  bool stop=false;
649  for(j=head_1[i1];j!=-1;j=next_1[j]){
650  for(k=0;k<3;k++){
651  if (GetId(triangles[(int)j/3][k])==i1){
652  if (GetId(triangles[(int)j/3][(int)((k+1)%3)])==i2){
653  bamgmesh->IssmSegments[num*4+0]=GetId(edges[i][0])+1; //back to M indexing
654  bamgmesh->IssmSegments[num*4+1]=GetId(edges[i][1])+1; //back to M indexing
655  bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
656  bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
657  num+=1;
658  stop=true;
659  break;
660  }
661  if (GetId(triangles[(int)j/3][(int)((k+2)%3)])==i2){
662  bamgmesh->IssmSegments[num*4+0]=GetId(edges[i][1])+1; //back to M indexing
663  bamgmesh->IssmSegments[num*4+1]=GetId(edges[i][0])+1; //back to M indexing
664  bamgmesh->IssmSegments[num*4+2]=(int)j/3+1; //back to M indexing
665  bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
666  num+=1;
667  stop=true;
668  break;
669  }
670  }
671  }
672  if(stop) break;
673  }
674  if (!stop){
675  _error_("Element holding segment [" << i1+1 << " " << i2+1 << "] not found...");
676  }
677  }
678  }
679 
680  /*Triangles*/
681  if(verbose>5) _printf_(" writing Triangles\n");
682  k=nbInT;
683  num=0;
684  bamgmesh->TrianglesSize[0]=k;
685  bamgmesh->TrianglesSize[1]=4;
686  if (k){
687  bamgmesh->Triangles=xNew<double>(4*k);
688  for (i=0;i<nbt;i++){
689  Triangle &t=triangles[i];
690  //reft[i]=-1 for outside triangle
691  if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
692  bamgmesh->Triangles[num*4+0]=GetId(t[0])+1; //back to M indexing
693  bamgmesh->Triangles[num*4+1]=GetId(t[1])+1; //back to M indexing
694  bamgmesh->Triangles[num*4+2]=GetId(t[2])+1; //back to M indexing
695  bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber;
696  num=num+1;
697  }
698  }
699  }
700 
701  /*SubDomains*/
702  if(verbose>5) _printf_(" writing SubDomains\n");
703  bamgmesh->SubDomainsSize[0]=nbsubdomains;
704  bamgmesh->SubDomainsSize[1]=4;
705  if (nbsubdomains){
706  bamgmesh->SubDomains=xNew<double>(4*nbsubdomains);
707  for (i=0;i<nbsubdomains;i++){
708  bamgmesh->SubDomains[i*4+0]=3;
709  bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing
710  bamgmesh->SubDomains[i*4+2]=1;
711  bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber;
712  }
713  }
714 
715  /*SubDomainsFromGeom*/
716  if(verbose>5) _printf_(" writing SubDomainsFromGeom\n");
718  bamgmesh->SubDomainsFromGeomSize[1]=4;
719  if (Gh.nbsubdomains){
720  bamgmesh->SubDomainsFromGeom=xNew<double>(4*Gh.nbsubdomains);
721  for (i=0;i<Gh.nbsubdomains;i++){
722  bamgmesh->SubDomainsFromGeom[i*4+0]=2;
723  bamgmesh->SubDomainsFromGeom[i*4+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
724  bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].direction;
725  bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ReferenceNumber;
726  }
727  }
728 
729  /*VerticesOnGeomVertex*/
730  if(verbose>5) _printf_(" writing VerticesOnGeomVertex\n");
732  bamgmesh->VerticesOnGeomVertexSize[1]=2;
734  bamgmesh->VerticesOnGeomVertex=xNew<double>(2*NbVerticesOnGeomVertex);
735  for (i=0;i<NbVerticesOnGeomVertex;i++){
736  VertexOnGeom &v=VerticesOnGeomVertex[i];
737  _assert_(v.OnGeomVertex());
738  bamgmesh->VerticesOnGeomVertex[i*2+0]=GetId((BamgVertex*)v)+1; //back to Matlab indexing
739  bamgmesh->VerticesOnGeomVertex[i*2+1]=Gh.GetId((GeomVertex*)v)+1; //back to Matlab indexing
740  }
741  }
742 
743  /*VertexOnGeomEdge*/
744  if(verbose>5) _printf_(" writing VerticesOnGeomEdge\n");
746  bamgmesh->VerticesOnGeomEdgeSize[1]=3;
748  bamgmesh->VerticesOnGeomEdge=xNew<double>(3*NbVerticesOnGeomEdge);
749  for (i=0;i<NbVerticesOnGeomEdge;i++){
750  const VertexOnGeom &v=VerticesOnGeomEdge[i];
751  if (!v.OnGeomEdge()){
752  _error_("A vertices supposed to be OnGeomEdge is actually not");
753  }
754  bamgmesh->VerticesOnGeomEdge[i*3+0]=GetId((BamgVertex*)v)+1; //back to Matlab indexing
755  bamgmesh->VerticesOnGeomEdge[i*3+1]=Gh.GetId((const GeomEdge*)v)+1; //back to Matlab indexing
756  bamgmesh->VerticesOnGeomEdge[i*3+2]=(double)v; //absisce
757  }
758  }
759 
760  /*EdgesOnGeomEdge*/
761  if(verbose>5) _printf_(" writing EdgesOnGeomEdge\n");
762  k=0;
763  for (i=0;i<nbe;i++){
764  if (edges[i].GeomEdgeHook) k=k+1;
765  }
766  bamgmesh->EdgesOnGeomEdgeSize[0]=k;
767  bamgmesh->EdgesOnGeomEdgeSize[1]=2;
768  if (k){
769  bamgmesh->EdgesOnGeomEdge=xNew<double>(2*(int)k);
770  int count=0;
771  for (i=0;i<nbe;i++){
772  if (edges[i].GeomEdgeHook){
773  bamgmesh->EdgesOnGeomEdge[count*2+0]=(double)i+1; //back to Matlab indexing
774  bamgmesh->EdgesOnGeomEdge[count*2+1]=(double)Gh.GetId(edges[i].GeomEdgeHook)+1; //back to Matlab indexing
775  count=count+1;
776  }
777  }
778  }
779 
780  /*Element Connectivity*/
781  if(verbose>5) _printf_(" writing Element connectivity\n");
782  bamgmesh->ElementConnectivitySize[0]=nbt-nbtout;
783  bamgmesh->ElementConnectivitySize[1]=3;
784  bamgmesh->ElementConnectivity=xNew<double>(3*(nbt-nbtout));
785  for (i=0;i<3*(nbt-nbtout);i++) bamgmesh->ElementConnectivity[i]=NAN;
786  num=0;
787  for (i=0;i<nbt;i++){
788  if (reft[i]>=0){
789  for (j=0;j<3;j++){
790  k=GetId(triangles[i].TriangleAdj(j));
791  if (reft[k]>=0){
792  _assert_(3*num+j<3*(nbt-nbtout));
793  bamgmesh->ElementConnectivity[3*num+j]=k+1; // back to Matlab indexing
794  }
795  }
796  num+=1;
797  }
798  }
799 
800  /*ElementNodal Connectivity*/
801  if(verbose>5) _printf_(" writing Nodal element connectivity\n");
802  bamgmesh->NodalElementConnectivitySize[0]=nbv;
803  bamgmesh->NodalElementConnectivitySize[1]=connectivitymax_1;
804  bamgmesh->NodalElementConnectivity=xNew<double>(connectivitymax_1*nbv);
805  for (i=0;i<connectivitymax_1*nbv;i++) bamgmesh->NodalElementConnectivity[i]=NAN;
806  for (i=0;i<nbv;i++){
807  k=0;
808  for(j=head_1[i];j!=-1;j=next_1[j]){
809  _assert_(connectivitymax_1*i+k < connectivitymax_1*nbv);
810  bamgmesh->NodalElementConnectivity[connectivitymax_1*i+k]=floor((double)j/3)+1;
811  k++;
812  }
813  }
814 
815  /*Nodal Connectivity*/
816  if(verbose>5) _printf_(" writing Nodal connectivity\n");
817  //chaining algorithm (again...)
818  int* head_2=NULL;
819  int* next_2=NULL;
820  int* connectivitysize_2=NULL;
821  int connectivitymax_2=0;
822  i1=bamgmesh->IssmEdgesSize[0];
823  i2=bamgmesh->IssmEdgesSize[1];
824  head_2=xNew<int>(nbv);
825  next_2=xNew<int>(2*i1);
826  connectivitysize_2=xNew<int>(nbv);
827  //Initialization
828  for (i=0;i<nbv;i++) head_2[i]=-1;
829  for (i=0;i<nbv;i++) connectivitysize_2[i]=0;
830  k=0;
831  //Chains generation
832  for (i=0;i<i1;i++) {
833  for (j=0;j<2;j++){
834  int v=(int)bamgmesh->IssmEdges[i*i2+j]-1; //back to C indexing
835  if (k>2*i1-1 || k<0) _error_("Index exceed matrix dimensions (k=" << k << " not in [0 " << 2*i1-1 << "]");
836  next_2[k]=head_2[v];
837  if (v>nbv-1 || v<0) _error_("Index exceed matrix dimensions (v=" << v << " not in [0 " << nbv-1 << "])");
838  head_2[v]=k++;
839  connectivitysize_2[v]+=1;
840  }
841  }
842  //Get maximum connectivity
843  for (i=0;i<nbv;i++){
844  if (connectivitysize_2[i]>connectivitymax_2) connectivitymax_2=connectivitysize_2[i];
845  }
846  //Build output
847  connectivitymax_2++;//add last column for size
848  bamgmesh->NodalConnectivitySize[0]=nbv;
849  bamgmesh->NodalConnectivitySize[1]=connectivitymax_2;
850  bamgmesh->NodalConnectivity=xNew<double>(connectivitymax_2*nbv);
851  for (i=0;i<connectivitymax_2*nbv;i++) bamgmesh->NodalConnectivity[i]=0;
852  for (i=0;i<nbv;i++){
853  k=0;
854  for(j=head_2[i];j!=-1;j=next_2[j]){
855  _assert_(connectivitymax_2*i+k < connectivitymax_2*nbv);
856  num=(int)bamgmesh->IssmEdges[int(j/2)*i2+0];
857  if (i+1==num){ //carefull, ElementEdge is in M indexing
858  //i is the first vertex of the edge, it is therefore connected to the second vertex
859  bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=bamgmesh->IssmEdges[int(j/2)*i2+1];
860  }
861  else{
862  bamgmesh->NodalConnectivity[connectivitymax_2*i+k]=num;
863  }
864  k++;
865  }
866  bamgmesh->NodalConnectivity[connectivitymax_2*(i+1)-1]=k;
867  }
868 
869  /*Cracked vertices*/
870  if(verbose>5) _printf_(" writing Cracked vertices\n");
872  bamgmesh->CrackedVerticesSize[1]=2;
873  if (NbCrackedVertices){
874  bamgmesh->CrackedVertices=xNew<double>(2*NbCrackedVertices);
875  for (i=0;i<NbCrackedVertices;i++){
876  bamgmesh->CrackedVertices[i*2+0]=CrackedVertices[i*2+0]+1; //M indexing
877  bamgmesh->CrackedVertices[i*2+1]=CrackedVertices[i*2+1]+1; //M indexing
878  }
879  }
880 
881  /*Cracked vertices*/
882  if(verbose>5) _printf_(" writing Cracked vertices\n");
883  bamgmesh->CrackedEdgesSize[0]=NbCrackedEdges;
884  bamgmesh->CrackedEdgesSize[1]=4;
885  if (NbCrackedEdges){
886  bamgmesh->CrackedEdges=xNew<double>(2*NbCrackedEdges);
887  for (i=0;i<NbCrackedEdges;i++){
888  bamgmesh->CrackedEdges[i*2+0]=0;//CrackedEdges[i]->+1; //M indexing
889  bamgmesh->CrackedEdges[i*2+1]=0;//CrackedEdges[i]-]->+1; //M indexing
890  }
891  }
892 
893  //clean up
894  xDelete<int>(connectivitysize_1);
895  xDelete<int>(head_1);
896  xDelete<int>(next_1);
897  xDelete<int>(connectivitysize_2);
898  xDelete<int>(head_2);
899  xDelete<int>(next_2);
900  delete [] reft;
901  delete [] numt;
902  }

◆ ReadMetric()

void bamg::Mesh::ReadMetric ( const BamgOpts bamgopts)

Definition at line 904 of file Mesh.cpp.

904  {/*{{{*/
905 
906  /*Intermediary*/
907  int i,j;
908  int verbose=0;
909 
910  /*Check pointer*/
911  _assert_(bamgopts);
912 
913  /*Get options*/
914  verbose=bamgopts->verbose;
915 
916  if(verbose>3) _printf_(" processing metric\n");
917  double hmin = Max(bamgopts->hmin,MinimalHmin());
918  double hmax = Min(bamgopts->hmax,MaximalHmax());
919  double coef = bamgopts->coeff;
920 
921  //for now we only use j==3
922  j=3;
923 
924  for (i=0;i<nbv;i++){
925  double h;
926  if (j == 1){
927  h=bamgopts->metric[i];
928  vertices[i].m=Metric(Max(hmin,Min(hmax, h*coef)));
929  }
930  else if (j==3){
931  //do not erase metric computed by hVertices
932  if (vertices[i].m.a11==1 && vertices[i].m.a21==0 && vertices[i].m.a22==1){
933  double a,b,c;
934  a=bamgopts->metric[i*3+0];
935  b=bamgopts->metric[i*3+1];
936  c=bamgopts->metric[i*3+2];
937  Metric M(a,b,c);
938  EigenMetric Vp(M/coef);
939 
940  Vp.Maxh(hmax);
941  Vp.Minh(hmin);
942  vertices[i].m = Vp;
943  }
944  }
945  }
946  }

◆ WriteMetric()

void bamg::Mesh::WriteMetric ( BamgOpts bamgopts)

Definition at line 948 of file Mesh.cpp.

948  {/*{{{*/
949  int i;
950  _assert_(bamgopts);
951  xDelete<double>(bamgopts->metric);
952  bamgopts->metric=xNew<double>(3*nbv);
953  for (i=0;i<nbv;i++){
954  bamgopts->metric[i*3+0]=vertices[i].m.a11;
955  bamgopts->metric[i*3+1]=vertices[i].m.a21;
956  bamgopts->metric[i*3+2]=vertices[i].m.a22;
957  }
958  }

◆ WriteIndex()

void bamg::Mesh::WriteIndex ( int **  pindex,
int *  pnels 
)

Definition at line 960 of file Mesh.cpp.

960  {/*{{{*/
961 
962  /*Intermediary*/
963  int i,k;
964 
965  /*output*/
966  int* index=NULL;
967  int num=0;
968 
969  /*Get number of triangles*/
970  k=0;
971  for (i=0;i<nbt;i++){
972  Triangle &t=triangles[i];
973  if(t.det>0) k++;
974  }
975 
976  if (k){
977  index=xNew<int>(3*k);
978  for (i=0;i<nbt;i++){
979  Triangle &t=triangles[i];
980  if (t.det>0 && !(t.Hidden(0)||t.Hidden(1) || t.Hidden(2) )){
981  /*Remove triangles that have a bad aspect ratio*/
982  //if(t.Anisotropy()<2 & t.Length()<1.e+5){
983  index[num*3+0]=GetId(t[0])+1; //back to M indexing
984  index[num*3+1]=GetId(t[1])+1; //back to M indexing
985  index[num*3+2]=GetId(t[2])+1; //back to M indexing
986  num=num+1;
987  //}
988  }
989  }
990  }
991 
992  /*Assign output pointers*/
993  *pindex=index;
994  *pnels=num;
995  }

◆ AddMetric()

void bamg::Mesh::AddMetric ( BamgOpts bamgopts)

Definition at line 999 of file Mesh.cpp.

999  {/*{{{*/
1000  // Hessiantype = 0 => H is computed using double L2 projection
1001  // Hessiantype = 1 => H is computed with green formula
1002 
1003  /*Check pointer*/
1004  _assert_(bamgopts);
1005 
1006  /*Options*/
1007  int Hessiantype=bamgopts->Hessiantype;
1008 
1009  if (Hessiantype==0){
1010  BuildMetric0(bamgopts);
1011  }
1012  else if (Hessiantype==1){
1013  BuildMetric1(bamgopts);
1014  }
1015  else{
1016  _error_("Hessiantype " << Hessiantype << " not supported yet (1->use Green formula, 0-> double L2 projection)");
1017  }
1018  }

◆ BuildMetric0()

void bamg::Mesh::BuildMetric0 ( BamgOpts bamgopts)

Definition at line 1638 of file Mesh.cpp.

1638  {/*{{{*/
1639 
1640  /*Options*/
1641  double* s=NULL;
1642  long nbsol;
1643  int verbose=0;
1644 
1645  int i,j,k,iA,iB,iC;
1646  int iv;
1647 
1648  /*Check pointer*/
1649  _assert_(bamgopts);
1650 
1651  /*Recover options*/
1652  verbose=bamgopts->verbose;
1653 
1654  /*Get and process fields*/
1655  s=bamgopts->field;
1656  nbsol=bamgopts->fieldSize[1];
1657 
1658  /*Check size*/
1659  if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
1660 
1661  //initialization of some variables
1662  double* ss=(double*)s;
1663  double sA,sB,sC;
1664  double* detT = new double[nbt];
1665  double* sumareas = new double[nbv];
1666  double* alpha= new double[nbt*3];
1667  double* beta = new double[nbt*3];
1668  double* dx_elem = new double[nbt];
1669  double* dy_elem = new double[nbt];
1670  double* dx_vertex = new double[nbv];
1671  double* dy_vertex = new double[nbv];
1672  double* dxdx_elem = new double[nbt];
1673  double* dxdy_elem = new double[nbt];
1674  double* dydy_elem = new double[nbt];
1675  double* dxdx_vertex= new double[nbv];
1676  double* dxdy_vertex= new double[nbv];
1677  double* dydy_vertex= new double[nbv];
1678 
1679  //display infos
1680  if(verbose>1) {
1681  _printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
1682  }
1683 
1684  //first, build the chains that will be used for the Hessian computation, as weel as the area of each element
1685  int* head_s=NULL;
1686  head_s=xNew<int>(nbv);
1687  int* next_p=NULL;
1688  next_p=xNew<int>(3*nbt);
1689  int p=0;
1690  //initialization
1691  for(i=0;i<nbv;i++){
1692  sumareas[i]=0;
1693  head_s[i]=-1;
1694  }
1695  for(i=0;i<nbt;i++){
1696 
1697  //lopp over the real triangles (no boundary elements)
1698  if(triangles[i].link){
1699 
1700  //get current triangle t
1701  const Triangle &t=triangles[i];
1702 
1703  // coor of 3 vertices
1704  R2 A=t[0];
1705  R2 B=t[1];
1706  R2 C=t[2];
1707 
1708  //compute triangle determinant (2*Area)
1709  double dett = bamg::Area2(A,B,C);
1710  detT[i]=dett;
1711 
1712  /*The nodal functions are such that for a vertex A:
1713  * N_A(x,y)=alphaA x + beta_A y +gamma_A
1714  * N_A(A) = 1, N_A(B) = 0, N_A(C) = 0
1715  * solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
1716  * leads to:
1717  * N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
1718  * and this gives:
1719  * alpha_A = (yB-yC)/(2*Area(T))*/
1720  alpha[i*3+0]=(B.y-C.y)/dett;
1721  alpha[i*3+1]=(C.y-A.y)/dett;
1722  alpha[i*3+2]=(A.y-B.y)/dett;
1723  beta[ i*3+0]=(C.x-B.x)/dett;
1724  beta[ i*3+1]=(A.x-C.x)/dett;
1725  beta[ i*3+2]=(B.x-A.x)/dett;
1726 
1727  //compute chains
1728  for(j=0;j<3;j++){
1729  k=GetId(triangles[i][j]);
1730  next_p[p]=head_s[k];
1731  head_s[k]=p++;
1732 
1733  //add area to sumareas
1734  sumareas[k]+=dett;
1735  }
1736 
1737  }
1738  }
1739 
1740  //for all Solutions
1741  for (int nusol=0;nusol<nbsol;nusol++) {
1742  double smin=ss[nusol],smax=ss[nusol];
1743 
1744  //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
1745  for ( iv=0,k=0; iv<nbv; iv++){
1746  smin=Min(smin,ss[iv*nbsol+nusol]);
1747  smax=Max(smax,ss[iv*nbsol+nusol]);
1748  }
1749  double sdelta=smax-smin;
1750  double absmax=Max(Abs(smin),Abs(smax));
1751 
1752  //display info
1753  if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << "\n");
1754 
1755  //skip constant field
1756  if (sdelta < 1.0e-10*Max(absmax,1e-20)){
1757  _printf_(" Solution " << nusol << " is constant, skipping...\n");
1758  continue;
1759  }
1760 
1761  //initialize the hessian and gradient matrices
1762  for ( iv=0,k=0; iv<nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
1763 
1764  //1: Compute gradient for each element (exact)
1765  for (i=0;i<nbt;i++){
1766  if(triangles[i].link){
1767  // number of the 3 vertices
1768  iA = GetId(triangles[i][0]);
1769  iB = GetId(triangles[i][1]);
1770  iC = GetId(triangles[i][2]);
1771 
1772  // value of the P1 fonction on 3 vertices
1773  sA = ss[iA*nbsol+nusol];
1774  sB = ss[iB*nbsol+nusol];
1775  sC = ss[iC*nbsol+nusol];
1776 
1777  //gradient = (sum alpha_i s_i, sum_i beta_i s_i)
1778  dx_elem[i]=sA*alpha[3*i+0]+sB*alpha[3*i+1]+sC*alpha[3*i+2];
1779  dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
1780  }
1781  }
1782 
1783  //2: then compute a gradient for each vertex using a P2 projection
1784  for(i=0;i<nbv;i++){
1785  for(p=head_s[i];p!=-1;p=next_p[p]){
1786  //Get triangle number
1787  k=(long)(p/3);
1788  dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
1789  dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
1790  }
1791  }
1792 
1793  //3: compute Hessian matrix on each element
1794  for (i=0;i<nbt;i++){
1795  if(triangles[i].link){
1796  // number of the 3 vertices
1797  iA = GetId(triangles[i][0]);
1798  iB = GetId(triangles[i][1]);
1799  iC = GetId(triangles[i][2]);
1800 
1801  //Hessian
1802  dxdx_elem[i]=dx_vertex[iA]*alpha[3*i+0]+dx_vertex[iB]*alpha[3*i+1]+dx_vertex[iC]*alpha[3*i+2];
1803  dxdy_elem[i]=dy_vertex[iA]*alpha[3*i+0]+dy_vertex[iB]*alpha[3*i+1]+dy_vertex[iC]*alpha[3*i+2];
1804  dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
1805  }
1806  }
1807 
1808  //4: finaly compute Hessian on each vertex using the second P2 projection
1809  for(i=0;i<nbv;i++){
1810  for(p=head_s[i];p!=-1;p=next_p[p]){
1811  //Get triangle number
1812  k=(long)(p/3);
1813  dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
1814  dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
1815  dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
1816  }
1817  }
1818 
1819  /*Compute Metric from Hessian*/
1820  for ( iv=0;iv<nbv;iv++){
1821  vertices[iv].MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
1822  }
1823 
1824  }//for all solutions
1825 
1826  //clean up
1827  xDelete<int>(head_s);
1828  xDelete<int>(next_p);
1829  delete [] detT;
1830  delete [] alpha;
1831  delete [] beta;
1832  delete [] sumareas;
1833  delete [] dx_elem;
1834  delete [] dy_elem;
1835  delete [] dx_vertex;
1836  delete [] dy_vertex;
1837  delete [] dxdx_elem;
1838  delete [] dxdy_elem;
1839  delete [] dydy_elem;
1840  delete [] dxdx_vertex;
1841  delete [] dxdy_vertex;
1842  delete [] dydy_vertex;
1843  }

◆ BuildMetric1()

void bamg::Mesh::BuildMetric1 ( BamgOpts bamgopts)

Definition at line 1845 of file Mesh.cpp.

1845  {/*{{{*/
1846  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectConsMetric)*/
1847 
1848  /*Options*/
1849  double* s=NULL;
1850  long nbsol;
1851  int NbJacobi;
1852  int verbose;
1853 
1854  /*Check pointer*/
1855  _assert_(bamgopts);
1856 
1857  /*Recover options*/
1858  verbose=bamgopts->verbose;
1859  NbJacobi=bamgopts->nbjacobi;
1860 
1861  /*Get and process fields*/
1862  s=bamgopts->field;
1863  nbsol=bamgopts->fieldSize[1];
1864 
1865  /*Check size*/
1866  if (bamgopts->fieldSize[0] != nbv) _error_("'field' should have " << nbv << " rows");
1867 
1868  //initialization of some variables
1869  long i,k,iA,iB,iC,iv;
1870  R2 O(0,0);
1871  double* ss=(double*)s;
1872  double sA,sB,sC;
1873  double* detT = new double[nbt];
1874  double* Mmass= new double[nbv];
1875  double* Mmassxx= new double[nbv];
1876  double* dxdx= new double[nbv];
1877  double* dxdy= new double[nbv];
1878  double* dydy= new double[nbv];
1879  double* workT= new double[nbt];
1880  double* workV= new double[nbv];
1881  int* OnBoundary = new int[nbv];
1882 
1883  //display infos
1884  if(verbose>1) {
1885  _printf_(" Construction of Metric: number of field: " << nbsol << " (nbt=" << nbt << ", nbv=" << nbv << ")\n");
1886  }
1887 
1888  //initialize Mmass, OnBoundary and Massxx by zero
1889  for (iv=0;iv<nbv;iv++){
1890  Mmass[iv]=0;
1891  OnBoundary[iv]=0;
1892  Mmassxx[iv]=0;
1893  }
1894 
1895  //Build detT Mmas Mmassxx workT and OnBoundary
1896  for (i=0;i<nbt;i++){
1897 
1898  //lopp over the real triangles (no boundary elements)
1899  if(triangles[i].link){
1900 
1901  //get current triangle t
1902  const Triangle &t=triangles[i];
1903 
1904  // coor of 3 vertices
1905  R2 A=t[0];
1906  R2 B=t[1];
1907  R2 C=t[2];
1908 
1909  // number of the 3 vertices
1910  iA = GetId(t[0]);
1911  iB = GetId(t[1]);
1912  iC = GetId(t[2]);
1913 
1914  //compute triangle determinant (2*Area)
1915  double dett = bamg::Area2(A,B,C);
1916  detT[i]=dett;
1917  dett /= 6;
1918 
1919  // construction of OnBoundary (flag=1 if on boundary, else 0)
1920  int nbb=0;
1921  for(int j=0;j<3;j++){
1922  //get adjacent triangle
1923  Triangle *ta=t.Adj(j);
1924  //if there is no adjacent triangle, the edge of the triangle t is on boundary
1925  if ( !ta || !ta->link){
1926  //mark the two vertices of the edge as OnBoundary
1927  OnBoundary[GetId(t[VerticesOfTriangularEdge[j][0]])]=1;
1928  OnBoundary[GetId(t[VerticesOfTriangularEdge[j][1]])]=1;
1929  nbb++;
1930  }
1931  }
1932 
1933  //number of vertices on boundary for current triangle t
1934  workT[i] = nbb;
1935 
1936  //Build Mmass Mmass[i] = Mmass[i] + Area/3
1937  Mmass[iA] += dett;
1938  Mmass[iB] += dett;
1939  Mmass[iC] += dett;
1940 
1941  //Build Massxx = Mmass
1942  Mmassxx[iA] += dett;
1943  Mmassxx[iB] += dett;
1944  Mmassxx[iC] += dett;
1945  }
1946 
1947  //else: the triangle is a boundary triangle -> workT=-1
1948  else workT[i]=-1;
1949  }
1950 
1951  //for all Solution
1952  for (int nusol=0;nusol<nbsol;nusol++) {
1953 
1954  double smin=ss[nusol],smax=ss[nusol];
1955  double h1=1.e30,h2=1e-30,rx=0;
1956  double hn1=1.e30,hn2=1e-30,rnx =1.e-30;
1957 
1958  //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy)
1959  for ( iv=0,k=0; iv<nbv; iv++ ){
1960  dxdx[iv]=dxdy[iv]=dydy[iv]=0;
1961  smin=Min(smin,ss[iv*nbsol+nusol]);
1962  smax=Max(smax,ss[iv*nbsol+nusol]);
1963  }
1964  double sdelta=smax-smin;
1965  double absmax=Max(Abs(smin),Abs(smax));
1966 
1967  //display info
1968  if(verbose>2) _printf_(" Solution " << nusol << ", Min = " << smin << ", Max = " << smax << ", Delta = " << sdelta << ", number of fields = " << nbsol << "\n");
1969 
1970  //skip constant field
1971  if (sdelta < 1.0e-10*Max(absmax,1e-20) ){
1972  if (verbose>2) _printf_(" Solution " << nusol << " is constant, skipping...\n");
1973  continue;
1974  }
1975 
1976  //pointer toward ss that is also a pointer toward s (solutions)
1977  double* sf=ss;
1978 
1979  //initialize the hessian matrix
1980  for ( iv=0,k=0; iv<nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
1981 
1982  //loop over the triangles
1983  for (i=0;i<nbt;i++){
1984 
1985  //for real all triangles
1986  if(triangles[i].link){
1987 
1988  // coor of 3 vertices
1989  R2 A=triangles[i][0];
1990  R2 B=triangles[i][1];
1991  R2 C=triangles[i][2];
1992 
1993  //warning: the normal is internal and the size is the length of the edge
1994  R2 nAB = Orthogonal(B-A);
1995  R2 nBC = Orthogonal(C-B);
1996  R2 nCA = Orthogonal(A-C);
1997  //note that : nAB + nBC + nCA == 0
1998 
1999  // number of the 3 vertices
2000  iA = GetId(triangles[i][0]);
2001  iB = GetId(triangles[i][1]);
2002  iC = GetId(triangles[i][2]);
2003 
2004  // for the test of boundary edge
2005  // the 3 adj triangles
2006  Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
2007  Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
2008  Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);
2009 
2010  // value of the P1 fonction on 3 vertices
2011  sA = ss[iA*nbsol+nusol];
2012  sB = ss[iB*nbsol+nusol];
2013  sC = ss[iC*nbsol+nusol];
2014 
2015  /*The nodal functions are such that for a vertex A:
2016  N_A(x,y)=alphaA x + beta_A y +gamma_A
2017  N_A(A) = 1, N_A(B) = 0, N_A(C) = 0
2018  solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined)
2019  leads to:
2020  N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T))
2021  and this gives:
2022  alpha_A = (yB-yC)/(2*Area(T))
2023  beta_A = (xC-xB)/(2*Area(T))
2024  and therefore:
2025  grad N_A = nA / detT
2026  for an interpolation of a solution s:
2027  grad(s) = s * sum_{i=A,B,C} grad(N_i) */
2028 
2029  R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];
2030 
2031  //Use Green to compute Hessian Matrix
2032 
2033  // if edge on boundary no contribution => normal = 0
2034  if ( !tBC || !tBC->link ) nBC=O;
2035  if ( !tCA || !tCA->link ) nCA=O;
2036  if ( !tAB || !tAB->link ) nAB=O;
2037 
2038  // remark we forgot a 1/2 because
2039  // int_{edge} w_i = 1/2 if i is in edge
2040  // 0 if not
2041  // if we don't take the boundary
2042  dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
2043  dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
2044  dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
2045 
2046  //warning optimization (1) the division by 2 is done on the metric construction
2047  dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
2048  dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
2049  dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ;
2050 
2051  dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
2052  dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
2053  dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
2054 
2055  } // for real all triangles
2056  }
2057 
2058  long kk=0;
2059  for ( iv=0,k=0 ; iv<nbv; iv++){
2060  if(Mmassxx[iv]>0){
2061  dxdx[iv] /= 2*Mmassxx[iv];
2062  // warning optimization (1) on term dxdy[iv]*ci/2
2063  dxdy[iv] /= 4*Mmassxx[iv];
2064  dydy[iv] /= 2*Mmassxx[iv];
2065  // Compute the matrix with abs(eigen value)
2066  Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
2067  EigenMetric Vp(M);
2068  Vp.Abs();
2069  M = Vp;
2070  dxdx[iv] = M.a11;
2071  dxdy[iv] = M.a21;
2072  dydy[iv] = M.a22;
2073  }
2074  else kk++;
2075  }
2076 
2077  // correction of second derivative
2078  // by a laplacien
2079  double* dd;
2080  for (int xy = 0;xy<3;xy++) {
2081  if (xy==0) dd=dxdx;
2082  else if (xy==1) dd=dxdy;
2083  else if (xy==2) dd=dydy;
2084  else{
2085  delete [] detT;
2086  delete [] Mmass;
2087  delete [] workT;
2088  delete [] workV;
2089  delete [] Mmassxx;
2090  delete [] OnBoundary;
2091  _error_("not supported yet");
2092  }
2093  // do leat 2 iteration for boundary problem
2094  for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
2095  for (i=0;i<nbt;i++)
2096  if(triangles[i].link){// the real triangles
2097  // number of the 3 vertices
2098  iA = GetId(triangles[i][0]);
2099  iB = GetId(triangles[i][1]);
2100  iC = GetId(triangles[i][2]);
2101  double cc=3;
2102  if(ijacobi==0)
2103  cc = Max((double) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
2104  workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
2105  }
2106  for (iv=0;iv<nbv;iv++) workV[iv]=0;
2107 
2108  for (i=0;i<nbt;i++){
2109  if(triangles[i].link){ // the real triangles
2110  // number of the 3 vertices
2111  iA = GetId(triangles[i][0]);
2112  iB = GetId(triangles[i][1]);
2113  iC = GetId(triangles[i][2]);
2114  double cc = workT[i]*detT[i];
2115  workV[iA] += cc;
2116  workV[iB] += cc;
2117  workV[iC] += cc;
2118  }
2119  }
2120 
2121  for (iv=0;iv<nbv;iv++){
2122  if( ijacobi<NbJacobi || OnBoundary[iv]){
2123  dd[iv] = workV[iv]/(Mmass[iv]*6);
2124  }
2125  }
2126  }
2127  }
2128 
2129  /*Compute Metric from Hessian*/
2130  for ( iv=0;iv<nbv;iv++){
2131  vertices[iv].MetricFromHessian(dxdx[iv],dxdy[iv],dydy[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->err[nusol],bamgopts);
2132  }
2133 
2134  }// end for all solution
2135 
2136  delete [] detT;
2137  delete [] Mmass;
2138  delete [] dxdx;
2139  delete [] dxdy;
2140  delete [] dydy;
2141  delete [] workT;
2142  delete [] workV;
2143  delete [] Mmassxx;
2144  delete [] OnBoundary;
2145 
2146  }

◆ BuildGeometryFromMesh()

void bamg::Mesh::BuildGeometryFromMesh ( BamgOpts bamgopts = NULL)

Definition at line 1197 of file Mesh.cpp.

1197  {/*{{{*/
1198  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ConsGeometry)*/
1199 
1200  /*Reconstruct Geometry from Mesh*/
1201 
1202  /*Intermediary*/
1203  int i,j,k,kk,it,jt;
1204  int verbose=0;
1205 
1206  /*Recover options*/
1207  if(bamgopts) verbose=bamgopts->verbose;
1208 
1209  //display info
1210  if (verbose>1) _printf_(" construction of the geometry from the 2d mesh\n");
1211 
1212  //check that the mesh is not empty
1213  if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)");
1214 
1215  /*Construction of the edges*/
1216 
1217  /*initialize st and edge4*/
1218  SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
1219  long* st = new long[nbt*3];
1220 
1221  /*initialize st as -1 (chaining algorithm)*/
1222  for (i=0;i<nbt*3;i++) st[i]=-1;
1223 
1224  /*build edge4 (chain)*/
1225  for (i=0;i<nbe;i++){
1226  edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1]));
1227  }
1228  /*check that there is no double edge*/
1229  if (nbe != edge4->nb()){
1230  delete [] st;
1231  _error_("Some Double edge in the mesh, the number is " << nbe << ", nbe4=" << edge4->nb());
1232  }
1233  //keep nbe in nbeold
1234  long nbeold = nbe;
1235 
1236  //Go through the triangles and add the edges in edge4 if they are not there yet
1237  for (i=0;i<nbt;i++){
1238  //3 edges per triangle
1239  for (j=0;j<3;j++) {
1240  //Add Edge to edge4 (k=numberofedges in edge4)
1241  long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]), GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
1242  long invisible = triangles[i].Hidden(j);
1243 
1244  //if st[k] has not been changed yet, add 3*i+j (= vertex position in the index)
1245  if(st[k]==-1) st[k]=3*i+j;
1246 
1247  //else st[k]>=0 -> the edge already exist, check
1248  else if(st[k]>=0) {
1249  //check that it is not an edge on boundary (should not already exist)
1250  if (triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))){
1251  _error_("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)");
1252  }
1253  /*OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list*/
1254  triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
1255  if (invisible) triangles[i].SetHidden(j);
1256  /* if k < nbe mark the edge as on Boundary (Locked)*/
1257  if (k<nbe) {
1258  triangles[i].SetLocked(j);
1259  }
1260  /*set st[k] as negative so that the edge should not be called again*/
1261  st[k]=-2-st[k];
1262  }
1263  else {
1264  /*else (see 3 lines above), the edge has been called more than twice: return error*/
1265  _printf_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")\n");
1266  _printf_("Edge " << j << " of triangle " << i << "\n");
1267  _printf_("Edge " << (-st[k]+2)%3 << " of triangle " << (-st[k]+2)/3 << "\n");
1268  _printf_("Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)) << " of triangle " << GetId(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << "\n");
1269  _error_("An edge belongs to more than 2 triangles");
1270  }
1271  }
1272  }
1273 
1274  //delete edge4
1275  long nbedges = edge4->nb(); // the total number of edges
1276  delete edge4; edge4=NULL;
1277 
1278  //display info
1279  if(verbose>5) {
1280  _printf_(" info on Mesh:\n");
1281  _printf_(" - number of vertices = " << nbv << "\n");
1282  _printf_(" - number of triangles = " << nbt << "\n");
1283  _printf_(" - number of given edges = " << nbe << "\n");
1284  _printf_(" - number of all edges = " << nbedges << "\n");
1285  _printf_(" - Euler number 1 - nb of holes = " << nbt-nbedges+nbv << "\n");
1286  }
1287 
1288  /*check consistency of edge[].adj and geometrical required vertices*/
1289  k=0; kk=0;
1290  for (i=0;i<nbedges;i++){
1291  //internal edge
1292  if (st[i] <-1) {
1293  //get triangle number back
1294  it = (-2-st[i])/3;
1295  //get edge position back
1296  j = (int) ((-2-st[i])%3);
1297  Triangle &tt=*triangles[it].TriangleAdj(j);
1298  if (triangles[it].color != tt.color|| i < nbeold) k++;
1299  }
1300  //boundary edge (alone)
1301  else if (st[i] >=0)
1302  kk++;
1303  }
1304 
1305  /*Constructions of edges*/
1306  k += kk;
1307  kk=0;
1308  if (k) {
1309  nbe = k;
1310  Edge* edgessave=edges;
1311  edges = new Edge[nbe];
1312  k =0;
1313 
1314  //display info
1315  if(verbose>4) _printf_(" Construction of the edges " << nbe << "\n");
1316 
1317  for (i=0;i<nbedges;i++){
1318  long add= -1;
1319 
1320  //internal edge (belongs to two triangles)
1321  if (st[i] <-1){
1322  it = (-2-st[i])/3;
1323  j = (int) ((-2-st[i])%3);
1324  Triangle & tt = * triangles[it].TriangleAdj(j);
1325  if (triangles[it].color != tt.color || i < nbeold) add=k++;
1326  }
1327  //boundary edge
1328  else if (st[i] >=0){
1329  it = st[i]/3;
1330  j = (int) (st[i]%3);
1331  add=k++;
1332  }
1333  if (add>=0 && add < nbe){
1334  edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]];
1335  edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]];
1336  edges[add].GeomEdgeHook=NULL;
1337  //if already existed
1338  if (i<nbeold){
1339  edges[add].ReferenceNumber=edgessave[i].ReferenceNumber;
1340  edges[add].GeomEdgeHook=edgessave[i].GeomEdgeHook; // HACK to get required edges
1341  _printf_("oh no...\n");
1342  }
1343  else
1344  edges[add].ReferenceNumber=Min(edges[add].v[0]->GetReferenceNumber(),edges[add].v[1]->GetReferenceNumber());
1345  }
1346  }
1347 
1348  //check that we have been through all edges
1349  if (k!=nbe){
1350  _error_("problem in edge construction process: k!=nbe (should not happen)");
1351  }
1352  //delete edgessave
1353  if (edgessave) delete [] edgessave;
1354  }
1355 
1356  /*Color the vertices*/
1357 
1358  //initialize color of all vertices as 0
1359  for (i=0;i<nbv;i++) vertices[i].color =0;
1360 
1361  //go through the edges and add a color to corresponding vertices
1362  //(A vertex in 4 edges will have a color 4)
1363  for (i=0;i<nbe;i++){
1364  for (j=0;j<2;j++) edges[i].v[j]->color++;
1365  }
1366 
1367  //change the color: if a vertex belongs to 2 edges -1, else -2
1368  for (i=0;i<nbv;i++) {
1369  vertices[i].color=(vertices[i].color ==2)? -1 : -2;
1370  }
1371 
1372  /*Build edges[i].adj: adjacency of each edge (if on the same curve)*/
1373  for (i=0;i<nbe;i++){
1374  for (j=0;j<2;j++){
1375  //get current vertex
1376  BamgVertex* v=edges[i].v[j];
1377  //get vertex color (i0)
1378  long i0=v->color;
1379  long j0;
1380 
1381  //if color<0 (first time), no adjacent edge
1382  if(i0<0) edges[i].adj[j]=NULL;
1383 
1384  //if color=-1 (corner),change the vertex color as 2*i+j (position of the vertex in edges)
1385  if(i0==-1) v->color=i*2+j;
1386 
1387  //if color>=0 (i and i0 edge are adjacent by the vertex v)
1388  else if (i0>=0) {
1389  //get position of v in edges back
1390  j0 = i0%2; //column in edges
1391  i0 = i0/2; //line in edges
1392 
1393  //check that we have the correct vertex
1394  if (v!=edges[i0 ].v[j0]){
1395  _error_("v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge");
1396  }
1397 
1398  //Add adjacence
1399  edges[i ].adj[j ]=edges +i0;
1400  edges[i0].adj[j0]=edges +i ;
1401 
1402  //change color to -3
1403  v->color = -3;
1404  }
1405  }
1406  }
1407 
1408  /*Reconstruct subdomains info*/
1409 
1410  //check that nbsubdomains is empty
1411  if(nbsubdomains) _error_("nbsubdomains should be 0");
1412  nbsubdomains=0;
1413 
1414  //color the subdomains
1415  long* colorT= new long[nbt];
1416  Triangle *tt,*t;
1417 
1418  //initialize the color of each triangle as -1
1419  for (it=0;it<nbt;it++) colorT[it]=-1;
1420 
1421  //loop over the triangles
1422  for (it=0;it<nbt;it++){
1423 
1424  //if the triangle has not been colored yet:
1425  if (colorT[it]<0){
1426 
1427  //color = number of subdomains
1428  colorT[it]=nbsubdomains;
1429 
1430  //color all the adjacent triangles of T that share a non marked edge
1431  int level =1;
1432  int kolor=triangles[it].color;
1433  st[0]=it; // stack
1434  st[1]=0;
1435  k=1;
1436  while (level>0){
1437  if( (j=st[level]++)<3 ){
1438  t = &triangles[st[level-1]];
1439  tt=t->TriangleAdj((int)j);
1440 
1441  //color the adjacent triangle
1442  if ( ! t->Locked(j) && tt && (colorT[jt = GetId(tt)] == -1) && ( tt->color==kolor)) {
1443  colorT[jt]=nbsubdomains;
1444  st[++level]=jt;
1445  st[++level]=0;
1446  k++;
1447  }
1448  }
1449  else level-=2;
1450  }
1451  nbsubdomains++;
1452  }
1453  }
1454  if (verbose> 3) _printf_(" The Number of sub domain = " << nbsubdomains << "\n");
1455 
1456  //build subdomains
1457  long isd;
1458  subdomains = new SubDomain[nbsubdomains];
1459 
1460  //initialize subdomains[isd].head as 0
1461  for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0;
1462 
1463  k=0;
1464  for (it=0;it<nbt;it++){
1465  for (int j=0;j<3;j++){
1466  tt=triangles[it].TriangleAdj(j);
1467  if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){
1468  subdomains[isd].head = triangles+it;
1470  subdomains[isd].direction = j; // hack
1471  subdomains[isd].edge = 0;
1472  k++;
1473  }
1474  }
1475  }
1476  //check that we have been through all subdomains
1477  if (k!= nbsubdomains){
1478  delete [] colorT;
1479  _error_("k!= nbsubdomains");
1480  }
1481  //delete colorT and st
1482  delete [] colorT;
1483  delete [] st;
1484 
1485  /*Reconstruct Geometry Gh*/
1486 
1487  //build colorV -1 for all vertex and 0 for the vertices belonging to edges
1488  long* colorV = new long[nbv];
1489  for (i=0;i<nbv;i++) colorV[i]=-1;
1490  for (i=0;i<nbe;i++){
1491  for ( j=0;j<2;j++) colorV[GetId(edges[i][j])]=0;
1492  }
1493  //number the vertices belonging to edges
1494  k=0;
1495  for (i=0;i<nbv;i++){
1496  if(!colorV[i]) colorV[i]=k++;
1497  }
1498 
1499  //Build Gh
1500  Gh.nbv=k;
1501  Gh.nbe = nbe;
1502  Gh.vertices = new GeomVertex[k];
1503  Gh.edges = new GeomEdge[nbe];
1505  Gh.subdomains = new GeomSubDomain[nbsubdomains];
1506  if (verbose>3) _printf_(" number of vertices = " << Gh.nbv << "\n number of edges = " << Gh.nbe << "\n");
1508  VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
1510  VerticesOnGeomEdge =0;
1511 
1512  //Build VertexOnGeom
1513  for (i=0;i<nbv;i++){
1514  if((j=colorV[i])>=0){
1515  BamgVertex & v = Gh.vertices[j];
1516  v = vertices[i];
1517  v.color =0;
1518  VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]);
1519  }
1520  }
1521 
1522  //Buid pmin and pmax of Gh (extrema coordinates)
1523  Gh.pmin = Gh.vertices[0].r;
1524  Gh.pmax = Gh.vertices[0].r;
1525  // recherche des extrema des vertices pmin,pmax
1526  for (i=0;i<Gh.nbv;i++) {
1527  Gh.pmin.x = Min(Gh.pmin.x,Gh.vertices[i].r.x);
1528  Gh.pmin.y = Min(Gh.pmin.y,Gh.vertices[i].r.y);
1529  Gh.pmax.x = Max(Gh.pmax.x,Gh.vertices[i].r.x);
1530  Gh.pmax.y = Max(Gh.pmax.y,Gh.vertices[i].r.y);
1531  }
1532  R2 DD05 = (Gh.pmax-Gh.pmin)*0.05;
1533  Gh.pmin -= DD05;
1534  Gh.pmax += DD05;
1535 
1536  //Build Gh.coefIcoor
1537  long MaxICoord = 1073741823; //2^30 - 1 = =111...111 (29 times one)
1538  Gh.coefIcoor= (MaxICoord)/(Max(Gh.pmax.x-Gh.pmin.x,Gh.pmax.y-Gh.pmin.y));
1539  if (Gh.coefIcoor<=0){
1540  delete [] colorV;
1541  _error_("Gh.coefIcoor<=0 in infered Geometry (this should not happen)");
1542  }
1543 
1544  /*Build Gh.edges*/
1545 
1546  //initialize len as 0
1547  double * len = new double[Gh.nbv];
1548  for(i=0;i<Gh.nbv;i++) len[i]=0;
1549 
1550  //initialize edge4 again
1551  edge4= new SetOfEdges4(nbe,nbv);
1552  double hmin = HUGE_VAL;
1553  int kreq=0;
1554  for (i=0;i<nbe;i++){
1555 
1556  long i0 = GetId(edges[i][0]);
1557  long i1 = GetId(edges[i][1]);
1558  long j0 = colorV[i0];
1559  long j1 = colorV[i1];
1560 
1561  Gh.edges[i].v[0] = Gh.vertices + j0;
1562  Gh.edges[i].v[1] = Gh.vertices + j1;
1563 
1564  Gh.edges[i].type = 0;
1565 
1566  Gh.edges[i].tg[0]=R2();
1567  Gh.edges[i].tg[1]=R2();
1568 
1569  bool required= edges[i].GeomEdgeHook;
1570  if(required) kreq++;
1571  edges[i].GeomEdgeHook = Gh.edges + i;
1572  if(required){
1573  Gh.edges[i].v[0]->SetRequired();
1574  Gh.edges[i].v[1]->SetRequired();
1575  Gh.edges[i].SetRequired();
1576  }
1577 
1578  R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r;
1579  double l12=Norme2(x12);
1580  hmin = Min(hmin,l12);
1581 
1582  Gh.vertices[j1].color++;
1583  Gh.vertices[j0].color++;
1584 
1585  len[j0]+= l12;
1586  len[j1] += l12;
1587  hmin = Min(hmin,l12);
1589 
1590  k = edge4->SortAndAdd(i0,i1);
1591  if (k != i){
1592  delete [] len;
1593  delete [] colorV;
1594  _error_("problem in Edge4 construction: k != i");
1595  }
1596  }
1597 
1598  //Build metric for all vertices of Gh
1599  for (i=0;i<Gh.nbv;i++){
1600  if (Gh.vertices[i].color > 0)
1601  Gh.vertices[i].m= Metric(len[i] /(double) Gh.vertices[i].color);
1602  else
1603  Gh.vertices[i].m= Metric(hmin);
1604  }
1605  //delete len
1606  delete [] len;
1607 
1608  //Build Gh.subdomains
1609  for (i=0;i<nbsubdomains;i++){
1610  it = GetId(subdomains[i].head);
1611  j = subdomains[i].direction;
1612  long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
1613  long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
1614  k = edge4->SortAndFind(i0,i1);
1615  if(k<0){
1616  delete [] colorV;
1617  _error_("This should not happen");
1618  }
1619  subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
1620  subdomains[i].edge = edges+k;
1621  Gh.subdomains[i].edge = Gh.edges + k;
1624  }
1625 
1626  delete edge4;
1627  delete [] colorV;
1628 
1629  //unset adj
1630  for (i=0;i<nbt;i++){
1631  for (j=0;j<3;j++){
1632  triangles[i].SetAdj2(j,0,triangles[i].GetAllflag(j));
1633  }
1634  }
1635 
1636  }

◆ RandomNumber()

long bamg::Mesh::RandomNumber ( long  max)

Definition at line 4756 of file Mesh.cpp.

4756  {/*{{{*/
4757  /* Generate a random number from 0 to max-1 using linear congruential
4758  * random number generator*/
4759 
4760  this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l;
4761  return this->randomseed / (714025l / max + 1);
4762  } /*}}}*/

◆ ReconstructExistingMesh()

void bamg::Mesh::ReconstructExistingMesh ( BamgOpts bamgopts)

Definition at line 3302 of file Mesh.cpp.

3302  {/*{{{*/
3303  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/
3304 
3305  /*This routine reconstruct an existing mesh to make it CONVEX:
3306  * -all the holes are filled
3307  * -concave boundaries are filled
3308  * A convex mesh is required for a lot of operations. This is why every mesh
3309  * goes through this process.
3310  * This routine also generates mesh properties such as adjencies,...
3311  */
3312 
3313  /*Intermediary*/
3314  int verbose=0;
3315 
3316  /*Get options*/
3317  if(bamgopts) verbose=bamgopts->verbose;
3318 
3319  // generation of the integer coordinate
3320 
3321  // find extrema coordinates of vertices pmin,pmax
3322  long i;
3323  if(verbose>2) _printf_(" Reconstruct mesh of " << nbv << " vertices\n");
3324 
3325  //initialize orderedvertices
3327  for (i=0;i<nbv;i++) orderedvertices[i]=0;
3328 
3329  //Initialize nbsubdomains
3330  nbsubdomains =0;
3331 
3332  /* generation of triangles adjacency*/
3333 
3334  //First add existing edges
3335  long kk =0;
3336  SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);
3337  for (i=0;i<nbe;i++){
3338  kk=kk+(i==edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])));
3339  }
3340  if (kk != nbe){
3341  _error_("There are " << kk-nbe << " double edges in the mesh");
3342  }
3343 
3344  //Add edges of all triangles in existing mesh
3345  long* st = new long[nbt*3];
3346  for (i=0;i<nbt*3;i++) st[i]=-1;
3347  for (i=0;i<nbt;i++){
3348  for (int j=0;j<3;j++){
3349 
3350  //Add current triangle edge to edge4
3351  long k =edge4->SortAndAdd(GetId(triangles[i][VerticesOfTriangularEdge[j][0]]),GetId(triangles[i][VerticesOfTriangularEdge[j][1]]));
3352 
3353  long invisible=triangles[i].Hidden(j);
3354 
3355  //If the edge has not been added to st, add it
3356  if(st[k]==-1) st[k]=3*i+j;
3357 
3358  //If the edge already exists, add adjacency
3359  else if(st[k]>=0) {
3360  _assert_(!triangles[i].TriangleAdj(j));
3361  _assert_(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));
3362 
3363  triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));
3364  if (invisible) triangles[i].SetHidden(j);
3365  if (k<nbe) triangles[i].SetLocked(j);
3366 
3367  //Make st[k] negative so that it will throw an error message if it is found again
3368  st[k]=-2-st[k];
3369  }
3370 
3371  //An edge belongs to 2 triangles
3372  else {
3373  _error_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << " , " << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles");
3374  }
3375  }
3376  }
3377 
3378  //Display info if required
3379  if(verbose>5) {
3380  _printf_(" info of Mesh:\n");
3381  _printf_(" - number of vertices = " << nbv << " \n");
3382  _printf_(" - number of triangles = " << nbt << " \n");
3383  _printf_(" - number of given edges = " << nbe << " \n");
3384  _printf_(" - number of all edges = " << edge4->nb() << "\n");
3385  _printf_(" - Euler number 1 - nb of holes = " << nbt-edge4->nb()+nbv << "\n");
3386  }
3387 
3388  //check the consistency of edge[].adj and the geometrical required vertex
3389  long k=0;
3390  for (i=0;i<edge4->nb();i++){
3391  if (st[i]>=0){ // edge alone
3392  if (i<nbe){
3393  long i0=edge4->i(i);
3394  orderedvertices[i0] = vertices+i0;
3395  long i1=edge4->j(i);
3396  orderedvertices[i1] = vertices+i1;
3397  }
3398  else {
3399  k=k+1;
3400  if (k<10) {
3401  //print only 10 edges
3402  _printf_("Lost boundary edges " << i << " : " << edge4->i(i) << " " << edge4->j(i) << "\n");
3403  }
3404  else if (k==10){
3405  _printf_("Other lost boundary edges not shown...\n");
3406  }
3407  }
3408  }
3409  }
3410  if(k) {
3411  _error_(k << " boundary edges (from the geometry) are not defined as mesh edges");
3412  }
3413 
3414  /* mesh generation with boundary points*/
3415  long nbvb=0;
3416  for (i=0;i<nbv;i++){
3417  vertices[i].t=0;
3419  if (orderedvertices[i]) orderedvertices[nbvb++]=orderedvertices[i];
3420  }
3421 
3422  Triangle* savetriangles=triangles;
3423  long savenbt=nbt;
3424  long savemaxnbt=maxnbt;
3425  SubDomain* savesubdomains=subdomains;
3426  subdomains=0;
3427 
3428  long Nbtriafillhole=2*nbvb;
3429  Triangle* triafillhole=new Triangle[Nbtriafillhole];
3430  triangles = triafillhole;
3431 
3432  nbt=2;
3433  maxnbt= Nbtriafillhole;
3434 
3435  //Find a vertex that is not aligned with vertices 0 and 1
3436  for (i=2;det(orderedvertices[0]->i,orderedvertices[1]->i,orderedvertices[i]->i)==0;)
3437  if (++i>=nbvb) {
3438  _error_("ReconstructExistingMesh: All the vertices are aligned");
3439  }
3440  //Move this vertex (i) to the 2d position in orderedvertices
3442 
3443  /*Reconstruct mesh beginning with 2 triangles*/
3444  BamgVertex * v0=orderedvertices[0], *v1=orderedvertices[1];
3445 
3446  triangles[0](0) = NULL; // Infinite vertex
3447  triangles[0](1) = v0;
3448  triangles[0](2) = v1;
3449 
3450  triangles[1](0) = NULL;// Infinite vertex
3451  triangles[1](2) = v0;
3452  triangles[1](1) = v1;
3453  const int e0 = OppositeEdge[0];
3454  const int e1 = NextEdge[e0];
3455  const int e2 = PreviousEdge[e0];
3456  triangles[0].SetAdj2(e0, &triangles[1] ,e0);
3457  triangles[0].SetAdj2(e1, &triangles[1] ,e2);
3458  triangles[0].SetAdj2(e2, &triangles[1] ,e1);
3459 
3460  triangles[0].det = -1; // boundary triangles
3461  triangles[1].det = -1; // boundary triangles
3462 
3465 
3466  triangles[0].link=&triangles[1];
3467  triangles[1].link=&triangles[0];
3468 
3469  if (!quadtree) delete quadtree; //ReInitialise;
3470  quadtree = new BamgQuadtree(this,0);
3471  quadtree->Add(*v0);
3472  quadtree->Add(*v1);
3473 
3474  // vertices are added one by one
3475  long NbSwap=0;
3476  for (int icount=2; icount<nbvb; icount++) {
3477  BamgVertex *vi = orderedvertices[icount];
3478  long long det3[3];
3479  Triangle *tcvi = TriangleFindFromCoord(vi->i,det3);
3480  quadtree->Add(*vi);
3481  AddVertex(*vi,tcvi,det3);
3482  NbSwap += vi->Optim(1,1);
3483  }
3484 
3485  //enforce the boundary
3486  AdjacentTriangle ta(0,0);
3487  long nbloss = 0,knbe=0;
3488  for ( i = 0; i < nbe; i++){
3489  if (st[i] >=0){ //edge alone => on border
3490  BamgVertex &a=edges[i][0], &b=edges[i][1];
3491  if (a.t && b.t){
3492  knbe++;
3493  if (ForceEdge(a,b,ta)<0) nbloss++;
3494  }
3495  }
3496  }
3497  if(nbloss) {
3498  _error_("we lost " << nbloss << " existing edges other " << knbe);
3499  }
3500 
3501  FindSubDomain(bamgopts,1);
3502  // remove all the hole
3503  // remove all the good sub domain
3504  long krm =0;
3505  for (i=0;i<nbt;i++){
3506  if (triangles[i].link){ // remove triangles
3507  krm++;
3508  for (int j=0;j<3;j++){
3509  AdjacentTriangle ta = triangles[i].Adj(j);
3510  Triangle &tta = *(Triangle*)ta;
3511  //if edge between remove and not remove
3512  if(! tta.link){
3513  // change the link of ta;
3514  int ja = ta;
3515  BamgVertex *v0= ta.EdgeVertex(0);
3516  BamgVertex *v1= ta.EdgeVertex(1);
3517  long k =edge4->SortAndAdd(v0?GetId(v0):nbv,v1? GetId(v1):nbv);
3518 
3519  _assert_(st[k]>=0);
3520  tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
3521  ta.SetLock();
3522  st[k]=-2-st[k];
3523  }
3524  }
3525  }
3526  }
3527  long NbTfillHoll =0;
3528  for (i=0;i<nbt;i++){
3529  if (triangles[i].link) {
3530  triangles[i]=Triangle((BamgVertex *) NULL,(BamgVertex *) NULL,(BamgVertex *) NULL);
3531  triangles[i].color=-1;
3532  }
3533  else{
3534  triangles[i].color= savenbt+ NbTfillHoll++;
3535  }
3536  }
3537  _assert_(savenbt+NbTfillHoll<=savemaxnbt);
3538 
3539  // copy of the outside triangles in saveMesh
3540  for (i=0;i<nbt;i++){
3541  if(triangles[i].color>=0) {
3542  savetriangles[savenbt]=triangles[i];
3543  savetriangles[savenbt].link=0;
3544  savenbt++;
3545  }
3546  }
3547  // gestion of the adj
3548  k =0;
3549  Triangle * tmax = triangles + nbt;
3550  for (i=0;i<savenbt;i++) {
3551  Triangle & ti = savetriangles[i];
3552  for (int j=0;j<3;j++){
3553  Triangle * ta = ti.TriangleAdj(j);
3554  int aa = ti.NuEdgeTriangleAdj(j);
3555  int lck = ti.Locked(j);
3556  if (!ta) k++; // bug
3557  else if ( ta >= triangles && ta < tmax){
3558  ta= savetriangles + ta->color;
3559  ti.SetAdj2(j,ta,aa);
3560  if(lck) ti.SetLocked(j);
3561  }
3562  }
3563  }
3564 
3565  // restore triangles;
3566  nbt=savenbt;
3567  maxnbt=savemaxnbt;
3568  delete [] triangles;
3569  delete [] subdomains;
3570  triangles = savetriangles;
3571  subdomains = savesubdomains;
3572  if (k) {
3573  _error_("number of triangles edges alone = " << k);
3574  }
3575  FindSubDomain(bamgopts);
3576 
3577  delete edge4;
3578  delete [] st;
3579  for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);
3580 
3581  SetVertexFieldOn();
3582 
3583  /*Check requirements consistency*/
3584  for (i=0;i<nbe;i++){
3585  /*If the current mesh edge is on Geometry*/
3586  if(edges[i].GeomEdgeHook){
3587  for(int j=0;j<2;j++){
3588  /*Go through the edges adjacent to current edge (if on the same curve)*/
3589  if (!edges[i].adj[j]){
3590  /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
3591  /*Check that the 2 vertices are on geometry AND required*/
3592  if(!edges[i][j].GeomEdgeHook->IsRequiredVertex()){
3593  _printf_("ReconstructExistingMesh error message: problem with the edge number " << i+1 << ": [" << GetId(edges[i][0])+1 << " " << GetId(edges[i][1])+1 << "]\n");
3594  _printf_("This edge is on geometrical edge number " << Gh.GetId(edges[i].GeomEdgeHook)+1 << "\n");
3595  if (edges[i][j].GeomEdgeHook->OnGeomVertex())
3596  _printf_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric BamgVertex number " << Gh.GetId(edges[i][j].GeomEdgeHook->gv)+1 << "\n");
3597  else if (edges[i][j].GeomEdgeHook->OnGeomEdge())
3598  _printf_("the vertex number " << GetId(edges[i][j])+1 << " of this edge is a geometric Edge number " << Gh.GetId(edges[i][j].GeomEdgeHook->ge)+1 << "\n");
3599  else
3600  _printf_("Its pointer is " << edges[i][j].GeomEdgeHook << "\n");
3601 
3602  _printf_("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
3603  _error_("See above (might be cryptic...)");
3604  }
3605  }
3606  }
3607  }
3608  }
3609  }

◆ CreateSingleVertexToTriangleConnectivity()

void bamg::Mesh::CreateSingleVertexToTriangleConnectivity ( )
inline

Definition at line 121 of file Mesh.h.

121  {
122  for (int i=0;i<nbv;i++) vertices[i].IndexInTriangle=0, vertices[i].t=NULL;
123  for (int i=0;i<nbt;i++) triangles[i].SetSingleVertexToTriangleConnectivity();
124  }

◆ UnMarkUnSwapTriangle()

void bamg::Mesh::UnMarkUnSwapTriangle ( )
inline

Definition at line 125 of file Mesh.h.

125  {
126  for (int i=0;i<nbt;i++)
127  for(int j=0;j<3;j++)
128  triangles[i].SetUnMarkUnSwap(j);
129  }

◆ SetVertexFieldOn()

void bamg::Mesh::SetVertexFieldOn ( )
inline

Definition at line 130 of file Mesh.h.

130  {
131  for (int i=0;i<nbv;i++) vertices[i].GeomEdgeHook=NULL;
132  for (int j=0;j<NbVerticesOnGeomVertex;j++) VerticesOnGeomVertex[j].SetOn();
133  for (int k=0;k<NbVerticesOnGeomEdge;k++ ) VerticesOnGeomEdge[k].SetOn();
134  }

◆ SetVertexFieldOnBTh()

void bamg::Mesh::SetVertexFieldOnBTh ( )
inline

Definition at line 135 of file Mesh.h.

135  {
136  for (int i=0;i<nbv;i++) vertices[i].GeomEdgeHook=NULL;
137  for (int j=0;j<NbVertexOnBThVertex;j++) VertexOnBThVertex[j].SetOnBTh();
138  for (int k=0;k<NbVertexOnBThEdge;k++ ) VertexOnBThEdge[k].SetOnBTh();
139  }

◆ TriangulateFromGeom1()

void bamg::Mesh::TriangulateFromGeom1 ( BamgOpts bamgopts,
int  KeepVertices = 1 
)
private

Definition at line 4426 of file Mesh.cpp.

4426  { /*{{{*/
4427  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles1)*/
4428 
4429  /*Intermediaries*/
4430  int verbose=0;
4431 
4432  /*Get options*/
4433  if(bamgopts) verbose=bamgopts->verbose;
4434 
4435  Gh.NbRef++;// add a ref to Gh
4436 
4437  /*************************************************************************
4438  * method in 2 steps
4439  * 1 - compute the number of new edges to allocate
4440  * 2 - construct the edges
4441  * remark:
4442  * in this part we suppose to have a background mesh with the same geometry
4443  *
4444  * To construct the discretization of the new mesh we have to
4445  * rediscretize the boundary of background Mesh
4446  * because we have only the pointeur from the background mesh to the geometry.
4447  * We need the abcisse of the background mesh vertices on geometry
4448  * so a vertex is
4449  * 0 on GeomVertex ;
4450  * 1 on GeomEdge + abcisse
4451  * 2 internal
4452  *************************************************************************/
4453 
4454  //Check that background mesh and current mesh do have the same geometry
4455  _assert_(&BTh.Gh==&Gh);
4456  BTh.NbRef++; // add a ref to BackGround Mesh
4457 
4458  //Initialize new mesh
4460  int* bcurve = new int[Gh.nbcurves]; //
4461 
4462  /* There are 2 ways to make the loop
4463  * 1) on the geometry
4464  * 2) on the background mesh
4465  * if you do the loop on geometry, we don't have the pointeur on background,
4466  * and if you do the loop in background we have the pointeur on geometry
4467  * so do the walk on background */
4468 
4471 
4472  /*STEP 1 copy of Required vertices*/
4473 
4474  int i;
4475  for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
4476  printf("\n");
4478  delete [] bcurve;
4479  _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
4480  }
4481 
4482  VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex];
4483  VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
4484 
4485  //At this point there is NO vertex but vertices should have been allocated by Init
4486  _assert_(vertices);
4487  for (i=0;i<Gh.nbv;i++){
4488  if (Gh[i].Required()) {//Gh vertices Required
4489  vertices[nbv] =Gh[i];
4490  vertices[nbv].i=I2(0,0);
4491  Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
4492  VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
4493  nbv++;
4494  }
4495  else Gh[i].MeshVertexHook=0;
4496  }
4497  for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){
4498  VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
4499  if (vog.IsRequiredVertex()){
4500  GeomVertex* gv=vog;
4501  BamgVertex *bv = vog;
4502  _assert_(gv->MeshVertexHook); // use of Geom -> Th
4503  VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
4504  gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
4505  }
4506  }
4507  _assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/
4508 
4509  /*STEP 2: reseed boundary edges*/
4510 
4511  // find the begining of the curve in BTh
4512  Gh.UnMarkEdges();
4513  int bfind=0;
4514  for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1;
4515 
4516  /*Loop over the backgrounf mesh BTh edges*/
4517  for (int iedge=0;iedge<BTh.nbe;iedge++){
4518  Edge &ei = BTh.edges[iedge];
4519 
4520  /*Loop over the 2 vertices of the current edge*/
4521  for(int je=0;je<2;je++){
4522 
4523  /* If one of the vertex is required we are in a new curve*/
4524  if (ei[je].GeomEdgeHook->IsRequiredVertex()){
4525 
4526  /*Get curve number*/
4527  int nc=ei.GeomEdgeHook->CurveNumber;
4528 
4529  //_printf_("Dealing with curve number " << nc << "\n");
4530  //_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
4531  //if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
4532  // _printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n");
4533  // _printf_("Do we have the right extremity? curve last vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n");
4534  //}
4535  //BUG FIX from original bamg
4536  /*Check that we are on the same edge and right vertex (0 or 1) */
4537  if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
4538  bcurve[nc]=iedge*2+je;
4539  bfind++;
4540  }
4541  else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
4542  bcurve[nc]=iedge*2+je;
4543  bfind++;
4544  }
4545  }
4546  }
4547  }
4548  if (bfind!=Gh.nbcurves){
4549  delete [] bcurve;
4550  _error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
4551  }
4552 
4553  // method in 2 + 1 step
4554  // 0.0) compute the length and the number of vertex to do allocation
4555  // 1.0) recompute the length
4556  // 1.1) compute the vertex
4557 
4558  long nbex=0,NbVerticesOnGeomEdgex=0;
4559  for (int step=0; step <2;step++){
4560 
4561  long NbOfNewPoints=0;
4562  long NbOfNewEdge=0;
4563  long iedge;
4564  Gh.UnMarkEdges();
4565  double L=0;
4566 
4567  /*Go through all geometrical curve*/
4568  for (int icurve=0;icurve<Gh.nbcurves;icurve++){
4569 
4570  /*Get edge and vertex (index) of background mesh on this curve*/
4571  iedge=bcurve[icurve]/2;
4572  int jedge=bcurve[icurve]%2;
4573 
4574  /*Get edge of Bth with index iedge*/
4575  Edge &ei = BTh.edges[iedge];
4576 
4577  /*Initialize variables*/
4578  double Lstep=0; // step between two points (phase==1)
4579  long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
4580 
4581  /*Do phase 0 to step*/
4582  for(int phase=0;phase<=step;phase++){
4583 
4584  /*Current curve pointer*/
4585  Curve *curve= Gh.curves+icurve;
4586 
4587  /*Get index of current curve*/
4588  int icurveequi= Gh.GetId(curve);
4589 
4590  /*For phase 0, check that we are at the begining of the curve only*/
4591  if(phase==0 && icurveequi!=icurve) continue;
4592 
4593  int k0=jedge,k1;
4594  Edge* pe= BTh.edges+iedge;
4595  int iedgeequi=bcurve[icurveequi]/2;
4596  int jedgeequi=bcurve[icurveequi]%2;
4597 
4598  int k0equi=jedgeequi,k1equi;
4599  Edge * peequi= BTh.edges+iedgeequi;
4600  GeomEdge *ongequi = peequi->GeomEdgeHook;
4601 
4602  double sNew=Lstep;// abscisse of the new points (phase==1)
4603  L=0;// length of the curve
4604  long i=0;// index of new points on the curve
4605  GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
4606  BamgVertex *A0;
4607  A0 = GA0->MeshVertexHook; // the vertex in new mesh
4608  BamgVertex *A1;
4609  VertexOnGeom *GA1;
4610  Edge* PreviousNewEdge = 0;
4611 
4612  // New Curve phase
4613  _assert_(A0-vertices>=0 && A0-vertices<nbv);
4614  if(ongequi->Required()){
4615  GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
4616  A1 = GA1->MeshVertexHook; //
4617  }
4618  else {
4619  for(;;){
4620  Edge &ee=*pe;
4621  Edge &eeequi=*peequi;
4622  k1 = 1-k0; // next vertex of the edge
4623  k1equi= 1 - k0equi;
4624  _assert_(pe && ee.GeomEdgeHook);
4625  ee.GeomEdgeHook->SetMark();
4626  BamgVertex & v0=ee[0], & v1=ee[1];
4627  R2 AB=(R2)v1-(R2)v0;
4628  double L0=L,LAB;
4629  LAB=LengthInterpole(v0.m,v1.m,AB);
4630  L+= LAB;
4631 
4632  if (phase){
4633  // computation of the new points for the given curve
4634  while ((i!=NbCreatePointOnCurve) && sNew<=L) {
4635 
4636  //some checks
4637  _assert_(sNew>=L0);
4638  _assert_(LAB);
4640  _assert_(edges && nbe<nbex);
4641  _assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
4642 
4643  // new vertex on edge
4644  A1=vertices+nbv++;
4646  Edge* e = edges + nbe++;
4647  double se= (sNew-L0)/LAB;
4648  if (se<0 || se>=1.000000001){
4649  _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
4650  }
4651  se = abscisseInterpole(v0.m,v1.m,AB,se,1);
4652  if (se<0 || se>1){
4653  _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
4654  }
4655  se = k1 ? se : 1. - se;
4656  se = k1==k1equi ? se : 1. - se;
4657  VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
4658  ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
4659  A1->ReferenceNumber = eeequi.ReferenceNumber;
4660  e->GeomEdgeHook = ongequi;
4661  e->v[0]=A0;
4662  e->v[1]=A1;
4663  e->ReferenceNumber = eeequi.ReferenceNumber;
4664  e->adj[0]=PreviousNewEdge;
4665 
4666  if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
4667  PreviousNewEdge=e;
4668  A0=A1;
4669  sNew += Lstep;
4670  if (++i== NbCreatePointOnCurve) break;
4671  }
4672  }
4673 
4674  //some checks
4675  _assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
4676  if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
4677  _assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
4678  GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
4679  A1=GA1->MeshVertexHook;// the vertex in new mesh
4680  _assert_(A1-vertices>=0 && A1-vertices<nbv);
4681  break;
4682  }
4683  if (!ee.adj[k1]) {
4684  _error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
4685  }
4686  pe = ee.adj[k1]; // next edge
4687  k0 = pe->Intersection(ee);
4688  peequi= eeequi.adj[k1equi]; // next edge
4689  k0equi=peequi->Intersection(eeequi);
4690  }// for(;;) end of the curve
4691  }
4692 
4693  if (phase){ // construction of the last edge
4694  Edge* e=edges + nbe++;
4695  e->GeomEdgeHook = ongequi;
4696  e->v[0]=A0;
4697  e->v[1]=A1;
4698  e->ReferenceNumber = peequi->ReferenceNumber;
4699  e->adj[0]=PreviousNewEdge;
4700  e->adj[1]=0;
4701  if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
4702  PreviousNewEdge = e;
4703 
4704  _assert_(i==NbCreatePointOnCurve);
4705  }
4706 
4707  if (!phase) { //
4708  long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
4709  Lstep = L/NbSegOnCurve;
4710  NbCreatePointOnCurve = NbSegOnCurve-1;
4711  NbOfNewEdge += NbSegOnCurve;
4712  NbOfNewPoints += NbCreatePointOnCurve;
4713  }
4714  }
4715  }// end of curve loop
4716 
4717  //Allocate memory
4718  if(step==0){
4719  if(nbv+NbOfNewPoints > maxnbv) {
4720  _error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
4721  }
4722  edges = new Edge[NbOfNewEdge];
4723  nbex = NbOfNewEdge;
4724  if(NbOfNewPoints) {
4725  VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
4726  NbVertexOnBThEdge = NbOfNewPoints;
4727  VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];
4728  NbVerticesOnGeomEdgex = NbOfNewPoints;
4729  }
4730  NbOfNewPoints =0;
4731  NbOfNewEdge = 0;
4732  }
4733  }
4734  _assert_(nbe!=0);
4735  delete [] bcurve;
4736 
4737  //Insert points inside existing triangles
4738  if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4739  if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
4740  if (verbose>3) _printf_(" Inserting boundary points\n");
4741  Insert(bamgopts);
4742 
4743  //Force the boundary
4744  if (verbose>3) _printf_(" Forcing boundaries\n");
4745  ForceBoundary(bamgopts);
4746 
4747  //Extract SubDomains
4748  if (verbose>3) _printf_(" Extracting subdomains\n");
4749  FindSubDomain(bamgopts);
4750 
4751  if (verbose>3) _printf_(" Inserting internal points\n");
4752  NewPoints(BTh,bamgopts,KeepVertices) ;
4753  if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4754  }

◆ TriangulateFromGeom0()

void bamg::Mesh::TriangulateFromGeom0 ( BamgOpts bamgopts)
private

Definition at line 4112 of file Mesh.cpp.

4112  {/*{{{*/
4113  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/
4114  /*Generate mesh from geometry*/
4115 
4116  /*Intermediaries*/
4117  int i,k;
4118  int verbose=0;
4119  int nbcurves = 0;
4120  int NbNewPoints,NbEdgeCurve;
4121  double lcurve,lstep,s;
4122  const int MaxSubEdge = 10;
4123 
4124  R2 AB;
4125  GeomVertex *a, *b;
4126  BamgVertex *va,*vb;
4127  GeomEdge *e;
4128 
4129  // add a ref to GH to make sure that it is not destroyed by mistake
4130  Gh.NbRef++;
4131 
4132  /*Get options*/
4133  if(bamgopts) verbose=bamgopts->verbose;
4134 
4135  //build background mesh flag (1 if background, else 0)
4136  bool background=(&BTh != this);
4137 
4138  /*Build VerticesOnGeomVertex*/
4139 
4140  //Compute the number of geometrical vertices that we are going to use to mesh
4141  for (i=0;i<Gh.nbv;i++){
4142  if (Gh[i].Required()) NbVerticesOnGeomVertex++;
4143  }
4144  //allocate
4145  VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
4146  if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
4147  _assert_(nbv==0);
4148  //Build VerticesOnGeomVertex
4149  for (i=0;i<Gh.nbv;i++){
4150  /* Add vertex only if required*/
4151  if (Gh[i].Required()) {//Gh vertices Required
4152 
4153  //Add the vertex
4154  _assert_(nbv<maxnbv);
4155  vertices[nbv]=Gh[i];
4156 
4157  //Add pointer from geometry (Gh) to vertex from mesh (Th)
4158  Gh[i].MeshVertexHook=vertices+nbv;
4159 
4160  //Build VerticesOnGeomVertex for current point
4161  VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
4162 
4163  //nbv increment
4164  nbv++;
4165  }
4166  }
4167 
4168  /*Build VerticesOnGeomEdge*/
4169 
4170  //check that edges is still empty (Init)
4171  _assert_(!edges);
4172 
4173  /* Now we are going to create the first edges corresponding
4174  * to the one present in the geometry provided.
4175  * We proceed in 2 steps
4176  * -step 0: we count all the edges
4177  * we allocate the number of edges at the end of step 0
4178  * -step 1: the edges are created */
4179  for (int step=0;step<2;step++){
4180 
4181  //initialize number of edges and number of edges max
4182  long nbex=0;
4183  nbe=0;
4184  long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
4185  Gh.UnMarkEdges();
4186  nbcurves=0;
4187 
4188  //go through the edges of the geometry
4189  for (i=0;i<Gh.nbe;i++){
4190 
4191  //ei = current Geometrical edge
4192  GeomEdge &ei=Gh.edges[i];
4193 
4194  //loop over the two vertices of the edge ei
4195  for(int j=0;j<2;j++) {
4196 
4197  /*Take only required vertices (corner->beginning of a new curve)*/
4198  if (!ei.Mark() && ei[j].Required()){
4199 
4200  long nbvend=0;
4201  Edge* PreviousNewEdge=NULL;
4202  lstep = -1;
4203 
4204  /*If Edge is required (do that only once for the 2 vertices)*/
4205  if(ei.Required()){
4206  if (j==0){
4207  //do not create internal points if required (take it as is)
4208  if(step==0) nbe++;
4209  else{
4210  e=&ei;
4211  a=ei(0);
4212  b=ei(1);
4213 
4214  //check that edges has been allocated
4215  _assert_(edges);
4216  edges[nbe].v[0]=a->MeshVertexHook;
4217  edges[nbe].v[1]=b->MeshVertexHook;;
4218  edges[nbe].ReferenceNumber = e->ReferenceNumber;
4219  edges[nbe].GeomEdgeHook = e;
4220  edges[nbe].adj[0] = 0;
4221  edges[nbe].adj[1] = 0;
4222  nbe++;
4223  }
4224  }
4225  }
4226 
4227  /*If Edge is not required: we are on a curve*/
4228  else {
4229  for (int kstep=0;kstep<=step;kstep++){
4230  //kstep=0, compute number of edges (discretize curve)
4231  //kstep=1 create the points and edge
4232  PreviousNewEdge=0;
4233  NbNewPoints=0;
4234  NbEdgeCurve=0;
4235  if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
4236  lcurve =0;
4237  s = lstep; //-1 initially, then length of each sub edge
4238 
4239  /*reminder: i = edge number, j=[0;1] vertex index in edge*/
4240  k=j; // k = vertex index in edge (0 or 1)
4241  e=&ei; // e = reference of current edge
4242  a=ei(k); // a = pointer toward the kth vertex of the current edge
4243  va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
4244  e->SetMark(); // Mark edge
4245 
4246  /*Loop until we reach the end of the curve*/
4247  for(;;){
4248  k = 1-k; // other vertx index of the curve
4249  b = (*e)(k); // b = pointer toward the other vertex of the current edge
4250  AB= b->r - a->r; // AB = vector of the current edge
4251  Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A
4252  Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
4253  double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2; //Get edge length in metric
4254 
4255  /* We are now creating the mesh edges from the geometrical edge selected above.
4256  * The edge will be divided according to the metric previously computed and cannot
4257  * be divided more than 10 times (MaxSubEdge). */
4258 
4259  //By default, there is only one subedge that is the geometrical edge itself
4260  int NbSubEdge = 1;
4261 
4262  //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
4263  double lSubEdge[MaxSubEdge];
4264 
4265  //Build Subedges according to the edge length
4266  if (ledge < 1.5){
4267  //if ledge < 1.5 (between one and 2), take the edge as is
4268  lSubEdge[0] = ledge;
4269  }
4270  //else, divide the edge
4271  else {
4272  //compute number of subedges (division of the edge), Maximum is 10
4273  NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
4274  /*Now, we are going to divide the edge according to the metric.
4275  * Get segment by sement along the edge.
4276  * Build lSubEdge, which holds the distance between the first vertex
4277  * of the edge and the next point on the edge according to the
4278  * discretization (each SubEdge is AB)*/
4279  R2 A,B;
4280  A=a->r;
4281  Metric MAs=MA,MBs;
4282  ledge=0;
4283  double x =0, xstep= 1./NbSubEdge;
4284  for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
4285  x += xstep;
4286  B = e->F(k ? x : 1-x);
4287  MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
4288  AB = A-B;
4289  lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2);
4290  }
4291  }
4292 
4293  double lcurveb = lcurve+ledge;
4294 
4295  /*Now, create corresponding points*/
4296  while(s>=lcurve && s<=lcurveb && nbv<nbvend){
4297 
4298  /*Schematic of current curve
4299  *
4300  * a vb b // vertex
4301  * 0 ll0 ll1 ledge // length from a
4302  * + --- + - ... - + --S-- + --- + - ... - + // where is S
4303  * 0 kk0 kk1 NbSubEdge // Sub edge index
4304  *
4305  */
4306 
4307  double ss = s-lcurve;
4308 
4309  /*Find the SubEdge containing ss using Dichotomy*/
4310  int kk0=-1,kk1=NbSubEdge-1,kkk;
4311  double ll0=0,ll1=ledge,llk;
4312  while (kk1-kk0>1){
4313  if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
4314  kk1=kkk,ll1=llk;
4315  else
4316  kk0=kkk,ll0=llk;
4317  }
4318  _assert_(kk1!=kk0);
4319 
4320  /*Curvilinear coordinate in [0 1] of ss in current edge*/
4321  // WARNING: This is what we would do
4322  // ssa = (ss-ll0)/(ll1-ll0);
4323  // aa = (kk0+ssa)/NbSubEdge
4324  // This is what Bamg does:
4325  double sbb = (ss-ll0)/(ll1-ll0);
4326  /*Curvilinear coordinate in [0 1] of ss in current curve*/
4327  double bb = (kk1+sbb)/NbSubEdge;
4328  double aa = 1-bb;
4329 
4330  // new vertex on edge
4331  vb = &vertices[nbv++];
4332  vb->m = Metric(aa,a->m,bb,b->m);
4333  vb->ReferenceNumber = e->ReferenceNumber;
4334  double abcisse = k ? bb : aa;
4335  vb->r = e->F(abcisse);
4336  VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);
4337 
4338  // to take into account the direction of the edge
4339  s += lstep;
4340  edges[nbe].v[0]=va;
4341  edges[nbe].v[1]=vb;
4342  edges[nbe].ReferenceNumber =e->ReferenceNumber;
4343  edges[nbe].GeomEdgeHook = e;
4344  edges[nbe].adj[0] = PreviousNewEdge;
4345  if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
4346  PreviousNewEdge=edges+nbe;
4347  nbe++;
4348  va = vb;
4349  }
4350 
4351  /*We just added one edge to the curve: Go to the next one*/
4352  lcurve = lcurveb;
4353  e->SetMark();
4354  a=b;
4355 
4356  /*If b is required, we are on a new curve->break*/
4357  if (b->Required()) break;
4358  int kprev=k;
4359  k = e->AdjVertexIndex[kprev];// next vertices
4360  e = e->Adj[kprev];
4361  _assert_(e);
4362  }// for(;;)
4363  vb = b->MeshVertexHook;
4364 
4365  /*Number of edges in the last disretized curve*/
4366  NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
4367  /*Number of internal vertices in the last disretized curve*/
4368  NbNewPoints = NbEdgeCurve-1;
4369  if(!kstep){
4370  NbVerticesOnGeomEdge0 += NbNewPoints;
4371  nbcurves++;
4372  }
4373  nbvend=nbv+NbNewPoints;
4374  lstep = lcurve / NbEdgeCurve; //approximately one
4375  }// end of curve --
4376  if (edges) { // last edges of the curves
4377  edges[nbe].v[0]=va;
4378  edges[nbe].v[1]=vb;
4379  edges[nbe].ReferenceNumber = e->ReferenceNumber;
4380  edges[nbe].GeomEdgeHook = e;
4381  edges[nbe].adj[0] = PreviousNewEdge;
4382  edges[nbe].adj[1] = 0;
4383  if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
4384  nbe++;
4385  }
4386  else nbe += NbEdgeCurve;
4387  } // end on curve ---
4388  }
4389  }
4390  } // for (i=0;i<nbe;i++)
4391  if(!step) {
4392  _assert_(!edges);
4394 
4395  edges = new Edge[nbex=nbe];
4396  if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
4397 
4398  // do the vertex on a geometrical vertex
4399  _assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
4400  NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
4401  }
4402  else{
4403  _assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
4404  }
4405  }
4406 
4407  //Insert points inside existing triangles
4408  if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4409  if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
4410  if (verbose>3) _printf_(" Inserting boundary points\n");
4411  Insert(bamgopts);
4412 
4413  //Force the boundary
4414  if (verbose>3) _printf_(" Forcing boundaries\n");
4415  ForceBoundary(bamgopts);
4416 
4417  //Extract SubDomains
4418  if (verbose>3) _printf_(" Extracting subdomains\n");
4419  FindSubDomain(bamgopts);
4420 
4421  if (verbose>3) _printf_(" Inserting internal points\n");
4422  NewPoints(*this,bamgopts,0) ;
4423  if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
4424  }

◆ Triangulate()

void bamg::Mesh::Triangulate ( double *  x,
double *  y,
int  nods,
BamgOpts bamgopts 
)
private

Definition at line 4084 of file Mesh.cpp.

4084  {/*{{{*/
4085 
4086  int verbose=0;
4087  int i;
4088  Metric M1(1);
4089 
4090  /*Get options*/
4091  if(bamgopts) verbose=bamgopts->verbose;
4092 
4093  /*Initialize mesh*/
4094  Init(nods);//this resets nbv to 0
4095  nbv=nods;
4096 
4097  //Vertices
4098  if(verbose) _printf_("Reading vertices (" << nbv << ")\n");
4099  for(i=0;i<nbv;i++){
4100  vertices[i].r.x=x[i];
4101  vertices[i].r.y=y[i];
4103  vertices[i].m=M1;
4104  vertices[i].color=0;
4105  }
4106  maxnbt=2*maxnbv-2; // for filling The Holes and quadrilaterals
4107 
4108  /*Insert Vertices*/
4109  Insert(bamgopts);
4110  }

◆ Init()

void bamg::Mesh::Init ( long  maxnbv_in)
private

Definition at line 2632 of file Mesh.cpp.

2632  {/*{{{*/
2633  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/
2634 
2635  /*Initialize fields*/
2636  this->NbRef = 0;
2637  this->quadtree = NULL;
2638  this->nbv = 0;
2639  this->nbt = 0;
2640  this->nbe = 0;
2641  this->edges = NULL;
2642  this->nbsubdomains = 0;
2643  this->subdomains = NULL;
2644  this->maxnbv = maxnbv_in;
2645  this->maxnbt = 2 *maxnbv_in-2;
2646  this->NbVertexOnBThVertex = 0;
2647  this->VertexOnBThVertex = NULL;
2648  this->NbVertexOnBThEdge = 0;
2649  this->VertexOnBThEdge = NULL;
2650  this->NbCrackedVertices = 0;
2651  this->CrackedVertices = NULL;
2652  this->NbCrackedEdges = 0;
2653  this->CrackedEdges = NULL;
2654  this->NbVerticesOnGeomVertex = 0;
2655  this->VerticesOnGeomVertex = NULL;
2656  this->NbVerticesOnGeomEdge = 0;
2657  this->VerticesOnGeomEdge = NULL;
2658 
2659  /*Initialize random seed*/
2660  this->randomseed = 1;
2661 
2662  /*Allocate if maxnbv_in>0*/
2663  if(maxnbv_in){
2664  this->vertices=new BamgVertex[this->maxnbv];
2665  this->orderedvertices=new BamgVertex* [this->maxnbv];
2666  this->triangles=new Triangle[this->maxnbt];
2667  _assert_(this->vertices);
2668  _assert_(this->orderedvertices);
2669  _assert_(this->triangles);
2670  }
2671  else{
2672  this->vertices = NULL;
2673  this->orderedvertices = NULL;
2674  this->triangles = NULL;
2675  this->maxnbt = 0;
2676  }
2677  }

◆ ForceEdge()

int bamg::Mesh::ForceEdge ( BamgVertex a,
BamgVertex b,
AdjacentTriangle taret 
)
private

Definition at line 4763 of file Mesh.cpp.

4763  { /*{{{*/
4764  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/
4765 
4766  int NbSwap =0;
4767  if (!a.t || !b.t){ // the 2 vertex is in a mesh
4768  _error_("!a.t || !b.t");
4769  }
4770  int k=0;
4771  taret=AdjacentTriangle(0,0); // erreur
4772 
4773  AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);
4774  BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
4775  // we turn around a in the direct direction
4776 
4777  long long det2 = v2 ? det(*v2,a,b): -1 , det1;
4778  if(v2) // normal case
4779  det2 = det(*v2,a,b);
4780  else { // no chance infini vertex try the next
4781  tta= Previous(Adj(tta));
4782  v2 = tta.EdgeVertex(0);
4783  vbegin =v2;
4784  if (!v2){
4785  _error_("!v2");
4786  }
4787  det2 = det(*v2,a,b);
4788  }
4789 
4790  while (v2 != &b) {
4791  AdjacentTriangle tc = Previous(Adj(tta));
4792  v1 = v2;
4793  v2 = tc.EdgeVertex(0);
4794  det1 = det2;
4795  det2 = v2 ? det(*v2,a,b): det2;
4796 
4797  if((det1 < 0) && (det2 >0)) {
4798  // try to force the edge
4799  BamgVertex * va = &a, *vb = &b;
4800  tc = Previous(tc);
4801  if (!v1 || !v2){
4802  _error_("!v1 || !v2");
4803  }
4804  long long detss = 0,l=0;
4805  while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))
4806  if(l++ > 10000000) {
4807  _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
4808  }
4809  BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
4810  if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){
4811  tc.SetLock();
4812  a.Optim(1,0);
4813  b.Optim(1,0);
4814  taret = tc;
4815  return NbSwap;
4816  }
4817  else
4818  {
4819  taret = tc;
4820  return -2; // error boundary is crossing
4821  }
4822  }
4823  tta = tc;
4824  k++;
4825  if (k>=2000){
4826  _error_("k>=2000");
4827  }
4828  if ( vbegin == v2 ) return -1;// error
4829  }
4830 
4831  tta.SetLock();
4832  taret=tta;
4833  a.Optim(1,0);
4834  b.Optim(1,0);
4835  return NbSwap;
4836  }

◆ SwapForForcingEdge()

int bamg::Mesh::SwapForForcingEdge ( BamgVertex *&  pva,
BamgVertex *&  pvb,
AdjacentTriangle tt1,
long long &  dets1,
long long &  detsa,
long long &  detsb,
int &  nbswap 
)
private

Definition at line 4838 of file Mesh.cpp.

4838  {/*{{{*/
4839  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/
4840  // l'arete ta coupe l'arete pva pvb
4841  // de cas apres le swap sa coupe toujours
4842  // on cherche l'arete suivante
4843  // on suppose que detsa >0 et detsb <0
4844  // attention la routine echange pva et pvb
4845 
4846  if(tt1.Locked()) return 0; // frontiere croise
4847 
4848  AdjacentTriangle tt2 = Adj(tt1);
4849  Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
4850  short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
4851  if ( a1<0 || a1>=3 ){
4852  _error_("a1<0 || a1>=3");
4853  }
4854 
4855  BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
4856  BamgVertex & s1= (*t1)[OppositeVertex[a1]];
4857  BamgVertex & s2= (*t2)[OppositeVertex[a2]];
4858 
4859  long long dets2 = det(*pva,*pvb,s2);
4860  long long det1=t1->det , det2=t2->det ;
4861  long long detT = det1+det2;
4862  if ((det1<=0 ) || (det2<=0)){
4863  _error_("(det1<=0 ) || (det2<=0)");
4864  }
4865  if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb
4866  _error_("(detsa>=0) || (detsb<=0)");
4867  }
4868  long long ndet1 = bamg::det(s1,sa,s2);
4869  long long ndet2 = detT - ndet1;
4870 
4871  int ToSwap =0; //pas de swap
4872  if ((ndet1 >0) && (ndet2 >0))
4873  { // on peut swaper
4874  if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
4875  ToSwap =1;
4876  else // swap alleatoire
4877  if (this->RandomNumber(1))
4878  ToSwap =2;
4879  }
4880  if (ToSwap) NbSwap++,
4881  bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
4882 
4883  int ret=1;
4884 
4885  if (dets2 < 0) {// haut
4886  dets1 = ToSwap ? dets1 : detsa ;
4887  detsa = dets2;
4888  tt1 = Previous(tt2) ;}
4889  else if (dets2 > 0){// bas
4890  dets1 = ToSwap ? dets1 : detsb ;
4891  detsb = dets2;
4892  //xxxx tt1 = ToSwap ? tt1 : Next(tt2);
4893  if(!ToSwap) tt1 = Next(tt2);
4894  }
4895  else { // changement de direction
4896  ret = -1;
4897  Exchange(pva,pvb);
4898  Exchange(detsa,detsb);
4899  Exchange(dets1,dets2);
4900  Exchange(tt1,tt2);
4901  dets1=-dets1;
4902  dets2=-dets2;
4903  detsa=-detsa;
4904  detsb=-detsb;
4905 
4906  if(ToSwap){
4907  if (dets2 < 0) {// haut
4908  dets1 = (ToSwap ? dets1 : detsa) ;
4909  detsa = dets2;
4910  tt1 = Previous(tt2) ;}
4911  else if(dets2 > 0){// bas
4912  dets1 = (ToSwap ? dets1 : detsb) ;
4913  detsb = dets2;
4914  if(!ToSwap) tt1 = Next(tt2);
4915  }
4916  else {// on a fin ???
4917  tt1 = Next(tt2);
4918  ret =0;}
4919  }
4920 
4921  }
4922  return ret;
4923  }

Field Documentation

◆ Gh

Geometry& bamg::Mesh::Gh

Definition at line 25 of file Mesh.h.

◆ BTh

Mesh& bamg::Mesh::BTh

Definition at line 26 of file Mesh.h.

◆ vertices

BamgVertex* bamg::Mesh::vertices

Definition at line 27 of file Mesh.h.

◆ triangles

Triangle* bamg::Mesh::triangles

Definition at line 28 of file Mesh.h.

◆ edges

Edge* bamg::Mesh::edges

Definition at line 29 of file Mesh.h.

◆ quadtree

BamgQuadtree* bamg::Mesh::quadtree

Definition at line 30 of file Mesh.h.

◆ orderedvertices

BamgVertex** bamg::Mesh::orderedvertices

Definition at line 31 of file Mesh.h.

◆ subdomains

SubDomain* bamg::Mesh::subdomains

Definition at line 32 of file Mesh.h.

◆ NbRef

long bamg::Mesh::NbRef

Definition at line 33 of file Mesh.h.

◆ maxnbv

long bamg::Mesh::maxnbv

Definition at line 34 of file Mesh.h.

◆ maxnbt

long bamg::Mesh::maxnbt

Definition at line 34 of file Mesh.h.

◆ nbv

long bamg::Mesh::nbv

Definition at line 35 of file Mesh.h.

◆ nbt

long bamg::Mesh::nbt

Definition at line 35 of file Mesh.h.

◆ nbe

long bamg::Mesh::nbe

Definition at line 35 of file Mesh.h.

◆ nbsubdomains

long bamg::Mesh::nbsubdomains

Definition at line 36 of file Mesh.h.

◆ nbtout

long bamg::Mesh::nbtout

Definition at line 37 of file Mesh.h.

◆ pmin

R2 bamg::Mesh::pmin

Definition at line 39 of file Mesh.h.

◆ pmax

R2 bamg::Mesh::pmax

Definition at line 39 of file Mesh.h.

◆ coefIcoor

double bamg::Mesh::coefIcoor

Definition at line 40 of file Mesh.h.

◆ lIntTria

ListofIntersectionTriangles bamg::Mesh::lIntTria

Definition at line 41 of file Mesh.h.

◆ randomseed

long bamg::Mesh::randomseed

Definition at line 42 of file Mesh.h.

◆ NbVerticesOnGeomVertex

long bamg::Mesh::NbVerticesOnGeomVertex

Definition at line 44 of file Mesh.h.

◆ VerticesOnGeomVertex

VertexOnGeom* bamg::Mesh::VerticesOnGeomVertex

Definition at line 45 of file Mesh.h.

◆ NbVerticesOnGeomEdge

long bamg::Mesh::NbVerticesOnGeomEdge

Definition at line 46 of file Mesh.h.

◆ VerticesOnGeomEdge

VertexOnGeom* bamg::Mesh::VerticesOnGeomEdge

Definition at line 47 of file Mesh.h.

◆ NbVertexOnBThVertex

long bamg::Mesh::NbVertexOnBThVertex

Definition at line 48 of file Mesh.h.

◆ VertexOnBThVertex

VertexOnVertex* bamg::Mesh::VertexOnBThVertex

Definition at line 49 of file Mesh.h.

◆ NbVertexOnBThEdge

long bamg::Mesh::NbVertexOnBThEdge

Definition at line 50 of file Mesh.h.

◆ VertexOnBThEdge

VertexOnEdge* bamg::Mesh::VertexOnBThEdge

Definition at line 51 of file Mesh.h.

◆ NbCrackedVertices

long bamg::Mesh::NbCrackedVertices

Definition at line 52 of file Mesh.h.

◆ CrackedVertices

long* bamg::Mesh::CrackedVertices

Definition at line 53 of file Mesh.h.

◆ NbCrackedEdges

long bamg::Mesh::NbCrackedEdges

Definition at line 54 of file Mesh.h.

◆ CrackedEdges

CrackedEdge* bamg::Mesh::CrackedEdges

Definition at line 55 of file Mesh.h.


The documentation for this class was generated from the following files:
bamg::Mesh::coefIcoor
double coefIcoor
Definition: Mesh.h:40
bamg::Mesh::subdomains
SubDomain * subdomains
Definition: Mesh.h:32
BamgGeom::Edges
double * Edges
Definition: BamgGeom.h:13
bamg::Mesh::GetId
long GetId(const Triangle &t) const
Definition: Mesh.cpp:2608
bamg::Mesh::MakeGeomEdgeToEdge
Edge ** MakeGeomEdgeToEdge()
Definition: Mesh.cpp:2899
bamg::CrackedEdge::E
GeomEdge * E
Definition: CrackedEdge.h:17
bamg::Mesh::MinimalHmin
double MinimalHmin()
Definition: Mesh.cpp:3031
bamg::BamgVertex::MetricFromHessian
void MetricFromHessian(const double Hxx, const double Hyx, const double Hyy, const double smin, const double smax, const double s, const double err, BamgOpts *bamgopts)
Definition: BamgVertex.cpp:34
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
bamg::GeomEdge::type
int type
Definition: GeomEdge.h:20
bamg::BamgVertex::i
I2 i
Definition: BamgVertex.h:20
IssmDouble
double IssmDouble
Definition: types.h:37
bamg::Mesh::NbVerticesOnGeomVertex
long NbVerticesOnGeomVertex
Definition: Mesh.h:44
bamg::BamgVertex::MeshVertexHook
BamgVertex * MeshVertexHook
Definition: BamgVertex.h:30
bamg::GeomSubDomain::edge
GeomEdge * edge
Definition: GeomSubDomain.h:13
BamgMesh::VerticesOnGeomVertexSize
int VerticesOnGeomVertexSize[2]
Definition: BamgMesh.h:19
bamg::Triangle::SetSingleVertexToTriangleConnectivity
void SetSingleVertexToTriangleConnectivity()
Definition: Triangle.cpp:213
bamg::BamgVertex::r
R2 r
Definition: BamgVertex.h:21
bamg::Curve::FirstVertexIndex
int FirstVertexIndex
Definition: Curve.h:16
bamg::Triangle::Hidden
int Hidden(int a) const
Definition: Triangle.cpp:105
BamgOpts::err
double * err
Definition: BamgOpts.h:47
bamg::PreviousEdge
static const short PreviousEdge[3]
Definition: macros.h:18
bamg::Geometry::pmin
R2 pmin
Definition: Geometry.h:32
BamgMesh::EdgesOnGeomEdge
double * EdgesOnGeomEdge
Definition: BamgMesh.h:24
bamg::Mesh::Triangulate
void Triangulate(double *x, double *y, int nods, BamgOpts *bamgopts)
Definition: Mesh.cpp:4084
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
bamg::Geometry::nbsubdomains
long nbsubdomains
Definition: Geometry.h:25
bamg::Edge::adj
Edge * adj[2]
Definition: Edge.h:18
bamg::BamgVertex::GetReferenceNumber
int GetReferenceNumber() const
Definition: BamgVertex.cpp:30
BamgMesh::NodalConnectivitySize
int NodalConnectivitySize[2]
Definition: BamgMesh.h:42
BamgMesh::NodalElementConnectivity
double * NodalElementConnectivity
Definition: BamgMesh.h:45
bamg::Mesh::SetVertexFieldOn
void SetVertexFieldOn()
Definition: Mesh.h:130
bamg::abscisseInterpole
double abscisseInterpole(const Metric &Ma, const Metric &Mb, R2 AB, double s, int optim)
Definition: Metric.cpp:327
BamgMesh::VerticesOnGeomEdgeSize
int VerticesOnGeomEdgeSize[2]
Definition: BamgMesh.h:21
bamg::Geometry::edges
GeomEdge * edges
Definition: Geometry.h:28
bamg::SubDomain::edge
Edge * edge
Definition: SubDomain.h:19
bamg::Mesh::RandomNumber
long RandomNumber(long max)
Definition: Mesh.cpp:4756
BamgOpts::verbose
int verbose
Definition: BamgOpts.h:26
bamg::BamgVertex::PreviousNumber
long PreviousNumber
Definition: BamgVertex.h:24
BamgMesh::NodalConnectivity
double * NodalConnectivity
Definition: BamgMesh.h:43
bamg::AdjacentTriangle::Adj
AdjacentTriangle Adj() const
Definition: AdjacentTriangle.cpp:28
bamg::Mesh::ForceEdge
int ForceEdge(BamgVertex &a, BamgVertex &b, AdjacentTriangle &taret)
Definition: Mesh.cpp:4763
BamgMesh::Vertices
double * Vertices
Definition: BamgMesh.h:12
bamg::Geometry::pmax
R2 pmax
Definition: Geometry.h:32
bamg::Mesh::ReadMesh
void ReadMesh(int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts)
Definition: Mesh.cpp:249
bamg::Mesh::NbRef
long NbRef
Definition: Mesh.h:33
bamg::SubDomain::head
Triangle * head
Definition: SubDomain.h:16
bamg::Mesh::NbCrackedVertices
long NbCrackedVertices
Definition: Mesh.h:52
BamgMesh::SubDomainsSize
int SubDomainsSize[2]
Definition: BamgMesh.h:26
bamg::Adj
AdjacentTriangle Adj(const AdjacentTriangle &a)
Definition: Mesh.h:164
bamg::PreviousVertex
static const short PreviousVertex[3]
Definition: macros.h:20
BamgOpts::field
double * field
Definition: BamgOpts.h:45
bamg::Geometry::vertices
GeomVertex * vertices
Definition: Geometry.h:27
BamgOpts::Hessiantype
int Hessiantype
Definition: BamgOpts.h:18
bamg::Mesh::CrackedEdges
CrackedEdge * CrackedEdges
Definition: Mesh.h:55
bamg::Mesh::pmax
R2 pmax
Definition: Mesh.h:39
bamg::Mesh::VertexOnBThEdge
VertexOnEdge * VertexOnBThEdge
Definition: Mesh.h:51
bamg::Mesh::MetricAt
Metric MetricAt(const R2 &)
Definition: Mesh.cpp:3011
bamg::Mesh::AddVertex
void AddVertex(BamgVertex &s, Triangle *t, long long *=0)
Definition: Mesh.cpp:1020
BamgOpts::coeff
double coeff
Definition: BamgOpts.h:15
bamg::Mesh::ForceBoundary
void ForceBoundary(BamgOpts *bamgopts)
Definition: Mesh.cpp:2318
BamgOpts::nbjacobi
int nbjacobi
Definition: BamgOpts.h:22
bamg::Mesh::VerticesOnGeomVertex
VertexOnGeom * VerticesOnGeomVertex
Definition: Mesh.h:45
bamg::GeomEdge::v
GeomVertex * v[2]
Definition: GeomEdge.h:14
bamg::Mesh::maxnbv
long maxnbv
Definition: Mesh.h:34
bamg::NextVertex
static const short NextVertex[3]
Definition: macros.h:19
bamg::Mesh::nbtout
long nbtout
Definition: Mesh.h:37
bamg::Mesh::lIntTria
ListofIntersectionTriangles lIntTria
Definition: Mesh.h:41
bamg::Curve::FirstEdge
GeomEdge * FirstEdge
Definition: Curve.h:14
BamgMesh::CrackedVertices
double * CrackedVertices
Definition: BamgMesh.h:31
BamgMesh::ElementConnectivity
double * ElementConnectivity
Definition: BamgMesh.h:41
Exchange
void Exchange(T &a, T &b)
Definition: Exchange.h:4
bamg::BamgQuadtree::TooClose
BamgVertex * TooClose(BamgVertex *, double, int, int)
Definition: BamgQuadtree.cpp:424
bamg::BigPrimeNumber
long BigPrimeNumber(long n)
Definition: BigPrimeNumber.cpp:6
bamg::Max
T Max(const T &a, const T &b)
Definition: extrema.h:7
BamgMesh::SubDomains
double * SubDomains
Definition: BamgMesh.h:27
BamgMesh::EdgesOnGeomEdgeSize
int EdgesOnGeomEdgeSize[2]
Definition: BamgMesh.h:23
bamg::Geometry::subdomains
GeomSubDomain * subdomains
Definition: Geometry.h:30
bamg::Mesh::SetIntCoor
void SetIntCoor(const char *from=0)
Definition: Mesh.cpp:3678
BamgMesh::VerticesSize
int VerticesSize[2]
Definition: BamgMesh.h:11
bamg::Mesh::TriangulateFromGeom0
void TriangulateFromGeom0(BamgOpts *bamgopts)
Definition: Mesh.cpp:4112
bamg::Triangle::link
Triangle * link
Definition: Triangle.h:25
bamg::CrackedEdge::length
double length
Definition: CrackedEdge.h:20
bamg::Mesh::triangles
Triangle * triangles
Definition: Mesh.h:28
bamg::Mesh::BuildMetric1
void BuildMetric1(BamgOpts *bamgopts)
Definition: Mesh.cpp:1845
BamgMesh::CrackedVerticesSize
int CrackedVerticesSize[2]
Definition: BamgMesh.h:30
bamg::OppositeEdge
static const short OppositeEdge[3]
Definition: macros.h:16
bamg::Mesh::BuildGeometryFromMesh
void BuildGeometryFromMesh(BamgOpts *bamgopts=NULL)
Definition: Mesh.cpp:1197
BamgMesh::VerticesOnGeomEdge
double * VerticesOnGeomEdge
Definition: BamgMesh.h:22
bamg::Geometry::ProjectOnCurve
GeomEdge * ProjectOnCurve(const Edge &, double, BamgVertex &, VertexOnGeom &) const
Definition: Geometry.cpp:781
bamg::Mesh::NbVertexOnBThVertex
long NbVertexOnBThVertex
Definition: Mesh.h:48
bamg::Triangle::Adj
AdjacentTriangle Adj(int i) const
Definition: Triangle.cpp:48
BamgOpts::fieldSize
int fieldSize[2]
Definition: BamgOpts.h:44
BamgMesh::EdgesSize
int EdgesSize[2]
Definition: BamgMesh.h:14
bamg::Mesh::ReconstructExistingMesh
void ReconstructExistingMesh(BamgOpts *bamgopts)
Definition: Mesh.cpp:3302
bamg::Mesh::CrackedVertices
long * CrackedVertices
Definition: Mesh.h:53
bamg::Area2
RR Area2(const P2< R, RR > a, const P2< R, RR > b, const P2< R, RR > c)
Definition: R2.h:81
bamg::Triangle::SetLocked
void SetLocked(int a)
Definition: Triangle.cpp:202
bamg::GeomEdge::SetRequired
void SetRequired()
Definition: GeomEdge.cpp:91
bamg::Geometry::PostRead
void PostRead(bool checkcurve=false)
Definition: Geometry.cpp:440
bamg::Mesh::NbVertexOnBThEdge
long NbVertexOnBThEdge
Definition: Mesh.h:50
bamg::Next
AdjacentTriangle Next(const AdjacentTriangle &ta)
Definition: Mesh.h:161
bamg::Triangle::color
long color
Definition: Triangle.h:26
BamgMesh::SubDomainsFromGeom
double * SubDomainsFromGeom
Definition: BamgMesh.h:29
bamg::GeomEdge::tg
R2 tg[2]
Definition: GeomEdge.h:17
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
bamg::Curve::LastEdge
GeomEdge * LastEdge
Definition: Curve.h:15
bamg::GeomEdge::ReferenceNumber
long ReferenceNumber
Definition: GeomEdge.h:15
bamg::Geometry::ReadGeometry
void ReadGeometry(BamgGeom *bamggeom, BamgOpts *bamgopts)
Definition: Geometry.cpp:54
bamg::Mesh::CreateSingleVertexToTriangleConnectivity
void CreateSingleVertexToTriangleConnectivity()
Definition: Mesh.h:121
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
bamg::Mesh::MaximalHmax
double MaximalHmax()
Definition: Mesh.cpp:2949
bamg::Norme2
RR Norme2(const P2< R, RR > x)
Definition: R2.h:93
bamg::I2
P2< int, long long > I2
Definition: typedefs.h:11
bamg::Geometry::curves
Curve * curves
Definition: Geometry.h:31
BamgOpts::hmax
double hmax
Definition: BamgOpts.h:35
bamg::Geometry::GetId
long GetId(const GeomVertex &t) const
Definition: Geometry.cpp:425
bamg::Mesh::BuildMetric0
void BuildMetric0(BamgOpts *bamgopts)
Definition: Mesh.cpp:1638
bamg::Mesh::nbe
long nbe
Definition: Mesh.h:35
bamg::CrackedEdge::normal
R2 normal
Definition: CrackedEdge.h:21
bamg::Mesh::maxnbt
long maxnbt
Definition: Mesh.h:34
BamgMesh::PreviousNumbering
double * PreviousNumbering
Definition: BamgMesh.h:13
bamg::ListofIntersectionTriangles::NewPoints
long NewPoints(BamgVertex *, long &nbv, long maxnbv)
Definition: ListofIntersectionTriangles.cpp:119
bamg::GeomEdge::CurveNumber
long CurveNumber
Definition: GeomEdge.h:16
bamg::Geometry::UnMarkEdges
void UnMarkEdges()
Definition: Geometry.cpp:926
BamgMesh::Edges
double * Edges
Definition: BamgMesh.h:15
bamg::LengthInterpole
double LengthInterpole(const Metric &Ma, const Metric &Mb, R2 AB)
Definition: Metric.cpp:158
bamg::Mesh::Init
void Init(long)
Definition: Mesh.cpp:2632
bamg::Mesh::nbt
long nbt
Definition: Mesh.h:35
BamgMesh::SubDomainsFromGeomSize
int SubDomainsFromGeomSize[2]
Definition: BamgMesh.h:28
bamg::IsVertexOnVertex
const int IsVertexOnVertex
Definition: macros.h:11
bamg::SubDomain::direction
int direction
Definition: SubDomain.h:18
BamgMesh::IssmSegmentsSize
int IssmSegmentsSize[2]
Definition: BamgMesh.h:38
bamg::Mesh::NewPoints
void NewPoints(Mesh &, BamgOpts *bamgopts, int KeepVertices=1)
Definition: Mesh.cpp:3040
bamg::BamgQuadtree::NearestVertex
BamgVertex * NearestVertex(int i, int j)
Definition: BamgQuadtree.cpp:227
bamg::BamgVertex::ReferenceNumber
long ReferenceNumber
Definition: BamgVertex.h:23
bamg::Mesh::TriangleIntNumbering
void TriangleIntNumbering(long *renumbering)
Definition: Mesh.cpp:4036
bamg::Mesh::BTh
Mesh & BTh
Definition: Mesh.h:26
BamgOpts::metric
double * metric
Definition: BamgOpts.h:43
bamg::Metric::a11
double a11
Definition: Metric.h:22
bamg::Mesh::edges
Edge * edges
Definition: Mesh.h:29
bamg::Edge::v
BamgVertex * v[2]
Definition: Edge.h:15
bamg::SubDomain::ReferenceNumber
long ReferenceNumber
Definition: SubDomain.h:17
bamg::Mesh::pmin
R2 pmin
Definition: Mesh.h:39
BamgMesh::TrianglesSize
int TrianglesSize[2]
Definition: BamgMesh.h:16
R
const double R
Definition: Gembx.cpp:30
bamg::EdgesVertexTriangle
static const short EdgesVertexTriangle[3][2]
Definition: macros.h:14
BamgMesh::VerticesOnGeomVertex
double * VerticesOnGeomVertex
Definition: BamgMesh.h:20
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
BamgOpts::hmin
double hmin
Definition: BamgOpts.h:34
bamg::Mesh::I2ToR2
R2 I2ToR2(const I2 &P) const
Definition: Mesh.cpp:3941
bamg::Triangle::det
long long det
Definition: Triangle.h:23
bamg::CloseBoundaryEdge
AdjacentTriangle CloseBoundaryEdge(I2 A, Triangle *t, double &a, double &b)
Definition: Mesh.cpp:4927
bamg::VerticesOfTriangularEdge
static const short VerticesOfTriangularEdge[3][2]
Definition: macros.h:13
bamg::GeomSubDomain::direction
int direction
Definition: GeomSubDomain.h:14
BamgMesh::IssmEdgesSize
int IssmEdgesSize[2]
Definition: BamgMesh.h:36
bamg::Mesh::R2ToI2
I2 R2ToI2(const R2 &P) const
Definition: Mesh.cpp:3937
bamg::Geometry::nbe
long nbe
Definition: Geometry.h:24
bamg::Mesh::NbVerticesOnGeomEdge
long NbVerticesOnGeomEdge
Definition: Mesh.h:46
bamg::R2
P2< double, double > R2
Definition: typedefs.h:12
bamg::Mesh::SwapForForcingEdge
int SwapForForcingEdge(BamgVertex *&pva, BamgVertex *&pvb, AdjacentTriangle &tt1, long long &dets1, long long &detsa, long long &detsb, int &nbswap)
Definition: Mesh.cpp:4838
bamg::Mesh::FindSubDomain
void FindSubDomain(BamgOpts *bamgopts, int OutSide=0)
Definition: Mesh.cpp:2363
bamg::P2::x
R x
Definition: R2.h:15
bamg::CrackedEdge::a
Triangle * a
Definition: CrackedEdge.h:15
bamg::Geometry::nbcurves
long nbcurves
Definition: Geometry.h:26
bamg::Mesh::TriangleReferenceList
long TriangleReferenceList(long *) const
Definition: Mesh.cpp:4046
bamg::Mesh::Gh
Geometry & Gh
Definition: Mesh.h:25
bamg::BamgVertex::color
long color
Definition: BamgVertex.h:29
bamg::GeomSubDomain::ReferenceNumber
long ReferenceNumber
Definition: GeomSubDomain.h:15
BamgMesh::NodalElementConnectivitySize
int NodalElementConnectivitySize[2]
Definition: BamgMesh.h:44
bamg::Edge::ReferenceNumber
long ReferenceNumber
Definition: Edge.h:16
bamg::AdjacentTriangle::EdgeVertex
BamgVertex * EdgeVertex(const int &i) const
Definition: AdjacentTriangle.cpp:32
bamg::Mesh::nbsubdomains
long nbsubdomains
Definition: Mesh.h:36
bamg::Orthogonal
P2< R, RR > Orthogonal(const P2< R, RR > x)
Definition: R2.h:97
bamg::ListofIntersectionTriangles::SplitEdge
void SplitEdge(Mesh &, const R2 &, const R2 &, int nbegin=0)
Definition: ListofIntersectionTriangles.cpp:190
bamg::P2::y
R y
Definition: R2.h:15
bamg::Mesh::InsertNewPoints
long InsertNewPoints(long nbvold, long &NbTSwap, BamgOpts *bamgopts)
Definition: Mesh.cpp:2815
bamg::Mesh::nbv
long nbv
Definition: Mesh.h:35
bamg::BamgQuadtree::Add
void Add(BamgVertex &w)
Definition: BamgQuadtree.cpp:104
BamgMesh::IssmSegments
double * IssmSegments
Definition: BamgMesh.h:39
bamg::Geometry::nbv
long nbv
Definition: Geometry.h:23
bamg::Triangle::TriangleAdj
Triangle * TriangleAdj(int i) const
Definition: Triangle.cpp:223
bamg::Abs
T Abs(const T &a)
Definition: Abs.h:5
bamg::Mesh::TriangleFindFromCoord
Triangle * TriangleFindFromCoord(const I2 &, long long[3], Triangle *tstart=0)
Definition: Mesh.cpp:3945
bamg::swap
void swap(Triangle *t1, short a1, Triangle *t2, short a2, BamgVertex *s1, BamgVertex *s2, long long det1, long long det2)
Definition: Mesh.cpp:4966
bamg::Metric::a22
double a22
Definition: Metric.h:22
bamg::D2xD2
P2xP2< double, double > D2xD2
Definition: Metric.h:12
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
bamg::Mesh::TriangulateFromGeom1
void TriangulateFromGeom1(BamgOpts *bamgopts, int KeepVertices=1)
Definition: Mesh.cpp:4426
bamg::Previous
AdjacentTriangle Previous(const AdjacentTriangle &ta)
Definition: Mesh.h:158
bamg::Mesh::orderedvertices
BamgVertex ** orderedvertices
Definition: Mesh.h:31
bamg::Metric::a21
double a21
Definition: Metric.h:22
bamg::Min
T Min(const T &a, const T &b)
Definition: extrema.h:6
bamg::GeomVertex::SetRequired
void SetRequired()
Definition: GeomVertex.cpp:28
BamgMesh::ElementConnectivitySize
int ElementConnectivitySize[2]
Definition: BamgMesh.h:40
bamg::Mesh::vertices
BamgVertex * vertices
Definition: Mesh.h:27
bamg::IsVertexOnEdge
const int IsVertexOnEdge
Definition: macros.h:12
bamg::Mesh::quadtree
BamgQuadtree * quadtree
Definition: Mesh.h:30
bamg::BamgVertex::m
Metric m
Definition: BamgVertex.h:22
bamg::Edge::GeomEdgeHook
GeomEdge * GeomEdgeHook
Definition: Edge.h:17
bamg::Mesh::NbCrackedEdges
long NbCrackedEdges
Definition: Mesh.h:54
bamg::OppositeVertex
static const short OppositeVertex[3]
Definition: macros.h:15
bamg::CrackedEdge::e1
Edge * e1
Definition: CrackedEdge.h:18
BamgMesh::CrackedEdges
double * CrackedEdges
Definition: BamgMesh.h:33
bamg::CrackedEdge::b
Triangle * b
Definition: CrackedEdge.h:16
bamg::Triangle::SetAdj2
void SetAdj2(short a, Triangle *t, short aat)
Definition: Triangle.cpp:182
bamg::Geometry::NbRef
long NbRef
Definition: Geometry.h:22
bamg::Curve::LastVertexIndex
int LastVertexIndex
Definition: Curve.h:17
bamg::Triangle::SetHidden
void SetHidden(int a)
Definition: Triangle.cpp:194
bamg::Mesh::Insert
void Insert(BamgOpts *bamgopts)
Definition: Mesh.cpp:2679
BamgMesh::Triangles
double * Triangles
Definition: BamgMesh.h:17
bamg::Mesh::VertexOnBThVertex
VertexOnVertex * VertexOnBThVertex
Definition: Mesh.h:49
bamg::Geometry::coefIcoor
double coefIcoor
Definition: Geometry.h:33
BamgMesh::IssmEdges
double * IssmEdges
Definition: BamgMesh.h:37
BamgMesh::CrackedEdgesSize
int CrackedEdgesSize[2]
Definition: BamgMesh.h:32
bamg::Mesh::randomseed
long randomseed
Definition: Mesh.h:42
bamg::BamgVertex::IndexInTriangle
short IndexInTriangle
Definition: BamgVertex.h:25
bamg::NextEdge
static const short NextEdge[3]
Definition: macros.h:17
bamg::BamgVertex::t
Triangle * t
Definition: BamgVertex.h:28
bamg::Mesh::VerticesOnGeomEdge
VertexOnGeom * VerticesOnGeomEdge
Definition: Mesh.h:47
bamg::BamgVertex::Optim
long Optim(int=1, int=0)
Definition: BamgVertex.cpp:112