Changeset 3243
- Timestamp:
- 03/10/10 09:47:13 (16 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 29 edited
-
BamgObjects.h (modified) (12 diffs)
-
Bamgx.cpp (modified) (1 diff)
-
meshtype.h (modified) (5 diffs)
-
objects/CrackedEdge.h (modified) (1 diff)
-
objects/DoubleAndInt4.h (modified) (1 diff)
-
objects/Edge.h (modified) (3 diffs)
-
objects/GeometricalEdge.cpp (modified) (7 diffs)
-
objects/GeometricalEdge.h (modified) (2 diffs)
-
objects/GeometricalSubDomain.h (modified) (1 diff)
-
objects/Geometry.cpp (modified) (21 diffs)
-
objects/Geometry.h (modified) (5 diffs)
-
objects/ListofIntersectionTriangles.cpp (modified) (10 diffs)
-
objects/ListofIntersectionTriangles.h (modified) (4 diffs)
-
objects/Metric.cpp (modified) (8 diffs)
-
objects/Metric.h (modified) (7 diffs)
-
objects/QuadTree.cpp (modified) (5 diffs)
-
objects/QuadTree.h (modified) (1 diff)
-
objects/SetOfE4.cpp (modified) (4 diffs)
-
objects/SetOfE4.h (modified) (2 diffs)
-
objects/SubDomain.h (modified) (1 diff)
-
objects/Triangle.cpp (modified) (2 diffs)
-
objects/Triangle.h (modified) (4 diffs)
-
objects/Triangles.cpp (modified) (156 diffs)
-
objects/Triangles.h (modified) (6 diffs)
-
objects/Vertex.cpp (modified) (5 diffs)
-
objects/Vertex.h (modified) (3 diffs)
-
objects/VertexOnEdge.h (modified) (2 diffs)
-
objects/VertexOnGeom.h (modified) (4 diffs)
-
objects/VertexOnVertex.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/BamgObjects.h
r3242 r3243 34 34 35 35 /*INLINE functions{{{1*/ 36 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8C[3]){36 inline double det3x3(double A[3] ,double B[3],double C[3]){ 37 37 return A[0]*( B[1]*C[2]-B[2]*C[1]) 38 38 - A[1]*( B[0]*C[2]-B[2]*C[0]) … … 57 57 if (deta <=0) ret = -1; 58 58 else { 59 Real8 a = sqrt((Real8) (ac,ac)),60 b = sqrt(( Real8) (bc,bc)),61 c = sqrt(( Real8) (ab,ab)),59 double a = sqrt((double) (ac,ac)), 60 b = sqrt((double) (bc,bc)), 61 c = sqrt((double) (ab,ab)), 62 62 p = a+b+c; 63 Real8h= Max(Max(a,b),c),ro=deta/p;63 double h= Max(Max(a,b),c),ro=deta/p; 64 64 ret = ro/h;} 65 65 return ret; … … 79 79 /*}}}1*/ 80 80 /*INLINE functions of CLASS VertexOnVertex{{{1*/ 81 inline void VertexOnVertex::Set(const Triangles & Th , Int4i,Triangles & ThNew) {81 inline void VertexOnVertex::Set(const Triangles & Th ,long i,Triangles & ThNew) { 82 82 *this = Th.VertexOnBThVertex[i]; 83 83 v = ThNew.vertices + Th.Number(v); … … 87 87 /*INLINE functions of CLASS Triangles{{{1*/ 88 88 inline void Triangles::ReMakeTriangleContainingTheVertex(){ 89 register Int4i;89 register long i; 90 90 for ( i=0;i<nbv;i++){ 91 91 vertices[i].vint=0; … … 97 97 inline void Triangles::UnMarkUnSwapTriangle() 98 98 { 99 register Int4i;99 register long i; 100 100 for ( i=0;i<nbt;i++) 101 101 for(int j=0;j<3;j++) … … 104 104 105 105 inline void Triangles::SetVertexFieldOn(){ 106 for ( Int4i=0;i<nbv;i++)106 for (long i=0;i<nbv;i++) 107 107 vertices[i].onGeometry=0; 108 for ( Int4j=0;j<NbVerticesOnGeomVertex;j++ )108 for (long j=0;j<NbVerticesOnGeomVertex;j++ ) 109 109 VerticesOnGeomVertex[j].SetOn(); 110 for ( Int4k=0;k<NbVerticesOnGeomEdge;k++ )110 for (long k=0;k<NbVerticesOnGeomEdge;k++ ) 111 111 VerticesOnGeomEdge[k].SetOn(); 112 112 } 113 113 inline void Triangles::SetVertexFieldOnBTh(){ 114 for ( Int4i=0;i<nbv;i++)114 for (long i=0;i<nbv;i++) 115 115 vertices[i].onGeometry=0; 116 for ( Int4j=0;j<NbVertexOnBThVertex;j++ )116 for (long j=0;j<NbVertexOnBThVertex;j++ ) 117 117 VertexOnBThVertex[j].SetOnBTh(); 118 for ( Int4k=0;k<NbVertexOnBThEdge;k++ )118 for (long k=0;k<NbVertexOnBThEdge;k++ ) 119 119 VertexOnBThEdge[k].SetOnBTh(); 120 120 … … 177 177 link = ThNew.triangles + Th.Number(link); 178 178 } 179 inline Triangle::Triangle(Triangles *Th, Int4 i,Int4 j,Int4k) {179 inline Triangle::Triangle(Triangles *Th,long i,long j,long k) { 180 180 Vertex *v=Th->vertices; 181 Int4nbv = Th->nbv;181 long nbv = Th->nbv; 182 182 if (i<0 || j<0 || k<0){ 183 183 throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0")); … … 211 211 /*}}}1*/ 212 212 /*INLINE functions of CLASS Edge{{{1*/ 213 inline void Edge::Set(const Triangles & Th , Int4i,Triangles & ThNew)213 inline void Edge::Set(const Triangles & Th ,long i,Triangles & ThNew) 214 214 { 215 215 *this = Th.edges[i]; … … 244 244 /*}}}1*/ 245 245 /*INLINE functions of CLASS SubDomain{{{1*/ 246 inline void SubDomain::Set(const Triangles & Th , Int4i,Triangles & ThNew)246 inline void SubDomain::Set(const Triangles & Th ,long i,Triangles & ThNew) 247 247 { 248 248 *this = Th.subdomains[i]; … … 265 265 /*}}}1*/ 266 266 /*INLINE functions of CLASS Vertex{{{1*/ 267 Int4inline Vertex::Optim(int i,int koption){268 Int4ret=0;267 long inline Vertex::Optim(int i,int koption){ 268 long ret=0; 269 269 if ( t && (vint >= 0 ) && (vint <3) ){ 270 270 ret = t->Optim(vint,koption); … … 281 281 /*}}}1*/ 282 282 /*INLINE functions of CLASS VertexOnEdge{{{1*/ 283 inline void VertexOnEdge::Set(const Triangles & Th , Int4i,Triangles & ThNew)283 inline void VertexOnEdge::Set(const Triangles & Th ,long i,Triangles & ThNew) 284 284 { 285 285 *this = Th.VertexOnBThEdge[i]; … … 335 335 336 336 /*Other prototypes IN TRIANGLES.CPP (TO BE REORGANIZED){{{1*/ 337 Int4 AGoodNumberPrimeWith(Int4n);337 long AGoodNumberPrimeWith(long n); 338 338 TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ; 339 339 TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b); -
issm/trunk/src/c/Bamgx/Bamgx.cpp
r3242 r3243 130 130 // init with hmax 131 131 Metric Mhmax(hmax); 132 for ( Int4iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;132 for (int iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax; 133 133 // change using hVertices if required 134 134 if (bamgmesh_in->hVertices){ -
issm/trunk/src/c/Bamgx/meshtype.h
r3242 r3243 17 17 18 18 //typedefs 19 typedef double Real8;20 typedef long Int4;21 19 typedef int Icoor1; 22 20 #if LONG_BIT > 63 //64 bits or more … … 25 23 typedef double Icoor2; 26 24 #endif 27 typedef P2<Icoor1,Icoor2> I2;28 typedef P2xP2<short, Int4> I2xI2;29 typedef P2< Real8,Real8>R2;30 typedef P2<double, Real8>R2xR2;25 typedef P2<Icoor1,Icoor2> I2; 26 typedef P2xP2<short,long> I2xI2; 27 typedef P2<double,double> R2; 28 typedef P2<double,double> R2xR2; 31 29 32 30 //Some parameters 33 const double Pi =3.141592653589793238462643383279502884197169399375105820974944592308;34 const float fPi=3.141592653589793238462643383279502884197169399375105820974944592308;35 const int IsVertexOnGeom = 8;36 const int IsVertexOnVertex = 16;37 const int IsVertexOnEdge = 32;38 static const short VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}};39 static const short EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}};40 static const short OppositeVertex[3] = {0,1,2};41 static const short OppositeEdge[3] = {0,1,2};42 static const short NextEdge[3] = {1,2,0};43 static const short PreviousEdge[3] = {2,0,1};44 static const short NextVertex[3] = {1,2,0};45 static const short PreviousVertex[3] = {2,0,1};31 const double Pi =3.141592653589793238462643383279502884197169399375105820974944592308; 32 const float fPi=3.141592653589793238462643383279502884197169399375105820974944592308; 33 const int IsVertexOnGeom = 8; 34 const int IsVertexOnVertex = 16; 35 const int IsVertexOnEdge = 32; 36 static const short VerticesOfTriangularEdge[3][2] = {{1,2},{2,0},{0,1}}; 37 static const short EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}}; 38 static const short OppositeVertex[3] = {0,1,2}; 39 static const short OppositeEdge[3] = {0,1,2}; 40 static const short NextEdge[3] = {1,2,0}; 41 static const short PreviousEdge[3] = {2,0,1}; 42 static const short NextVertex[3] = {1,2,0}; 43 static const short PreviousVertex[3] = {2,0,1}; 46 44 #if LONG_BIT > 63 47 45 const Icoor1 MaxICoor = 1073741823; // 2^30-1 =111...111 (29 times) … … 61 59 template<class T> inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} 62 60 template<class T> inline void HeapSort(T *c,long n){ 63 /*Intermediary*/64 61 int l,j,r,i; 65 62 T crit; … … 87 84 } 88 85 template<class T> inline void HeapSort(int** porder,T* c,int n){ 89 /*Intermediary*/90 86 int l,j,r,i; 91 87 T crit; … … 142 138 return bax*cay - bay*cax; 143 139 } 144 145 140 } 146 141 #endif -
issm/trunk/src/c/Bamgx/objects/CrackedEdge.h
r3237 r3243 73 73 //Constructors 74 74 CrackedEdge() :a(),b() {} 75 CrackedEdge(Edge * start, Int4 i,Int4j) : a(start+i),b(start+j) {};75 CrackedEdge(Edge * start, long i,long j) : a(start+i),b(start+j) {}; 76 76 CrackedEdge(Edge * e0, Edge * e1 ) : a(e0),b(e1) {}; 77 77 -
issm/trunk/src/c/Bamgx/objects/DoubleAndInt4.h
r3232 r3243 11 11 public: 12 12 double q; 13 Int4i3j;13 long i3j; 14 14 15 15 //Operators -
issm/trunk/src/c/Bamgx/objects/Edge.h
r3237 r3243 21 21 public: 22 22 Vertex* v[2]; 23 Int4ref;23 long ref; 24 24 GeometricalEdge* onGeometry; 25 25 Edge* adj[2]; // the 2 adj edges if on the same curve … … 32 32 33 33 //Methods 34 void ReNumbering(Vertex *vb,Vertex *ve, Int4*renu){34 void ReNumbering(Vertex *vb,Vertex *ve, long *renu){ 35 35 if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb]; 36 36 if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb]; … … 47 47 48 48 //Inline methods 49 inline void Set(const Triangles &, Int4,Triangles &);49 inline void Set(const Triangles &,long,Triangles &); 50 50 51 51 }; -
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
r3232 r3243 17 17 /*Methods*/ 18 18 /*FUNCTION GeometricalEdge::R1tg{{{1*/ 19 Real8 GeometricalEdge::R1tg(Real8theta,R2 & t) const // 1/R of radius of cuvature19 double GeometricalEdge::R1tg(double theta,R2 & t) const // 1/R of radius of cuvature 20 20 { R2 A=v[0]->r,B=v[1]->r; 21 Real8dca,dcb,dcta,dctb;22 Real8ddca,ddcb,ddcta,ddctb;23 // Real8t1 = 1 -theta;24 // Real8t1t1 = t1*t1;25 Real8tt = theta*theta;21 double dca,dcb,dcta,dctb; 22 double ddca,ddcb,ddcta,ddctb; 23 // double t1 = 1 -theta; 24 // double t1t1 = t1*t1; 25 double tt = theta*theta; 26 26 if ( theta<0){ 27 27 throw ErrorException(__FUNCT__,exprintf("theta<0")); … … 48 48 } 49 49 else { // 1-t*t, t-t*t, t*t 50 Real8t = theta;50 double t = theta; 51 51 // cb = t*t; 52 52 dcb = 2*t; … … 64 64 else 65 65 if (TgB()){ 66 Real8t = 1-theta;66 double t = 1-theta; 67 67 //ca = t*t; 68 68 dca = -2*t; … … 82 82 83 83 R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; 84 Real8d2=(d,d);85 Real8sd2 = sqrt(d2);84 double d2=(d,d); 85 double sd2 = sqrt(d2); 86 86 t=d; 87 87 if(d2>1.0e-20) {t/=sd2;return Abs(Det(d,dd))/(d2*sd2);} … … 90 90 /*}}}1*/ 91 91 /*FUNCTION GeometricalEdge::F{{{1*/ 92 R2 GeometricalEdge::F( Real8theta) const // parametrization of the curve edge92 R2 GeometricalEdge::F(double theta) const // parametrization of the curve edge 93 93 { R2 A=v[0]->r,B=v[1]->r; 94 Real8ca,cb,cta,ctb;94 double ca,cb,cta,ctb; 95 95 if ( theta<-1e-12){ 96 96 throw ErrorException(__FUNCT__,exprintf("theta<-1e-12")); … … 107 107 } 108 108 else { // 1-t*t, t-t*t, t*t 109 Real8t = theta;109 double t = theta; 110 110 cb = t*t; 111 111 ca = 1-cb; … … 115 115 else 116 116 if (TgB()){ 117 Real8t = 1-theta;117 double t = 1-theta; 118 118 ca = t*t; 119 119 cb = 1-ca; -
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h
r3232 r3243 21 21 public: 22 22 GeometricalVertex* v[2]; 23 Int4ref;24 Int4CurveNumber;23 long ref; 24 long CurveNumber; 25 25 R2 tg[2]; // the 2 tangentes (tg[0] =0 => no continuity) 26 26 GeometricalEdge* Adj[2]; … … 35 35 36 36 //Methods 37 R2 F( Real8theta) const ; // parametrization of the curve edge38 Real8 R1tg(Real8theta,R2 &t) const ; // 1/radius of curvature + tangente37 R2 F(double theta) const ; // parametrization of the curve edge 38 double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente 39 39 int Cracked() const {return flag & 1;} 40 40 int Dup() const { return flag & 32;} -
issm/trunk/src/c/Bamgx/objects/GeometricalSubDomain.h
r3232 r3243 19 19 GeometricalEdge *edge; 20 20 int sens; // -1 or 1 21 Int4ref;21 long ref; 22 22 23 23 //Inline methods -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r3242 r3243 16 16 /*FUNCTION Geometry::Geometry(const Geometry & Gh){{{1*/ 17 17 Geometry::Geometry(const Geometry & Gh) { 18 Int4i;18 long i; 19 19 *this = Gh; 20 20 NbRef =0; … … 223 223 NbOfCurves=0; 224 224 225 Real8Hmin = HUGE_VAL;// the infinie value226 Int4hvertices=0;225 double Hmin = HUGE_VAL;// the infinie value 226 long hvertices=0; 227 227 int i,j,k,n,i1,i2; 228 228 … … 250 250 vertices[i].r.x=(double)bamggeom->Vertices[i*3+0]; 251 251 vertices[i].r.y=(double)bamggeom->Vertices[i*3+1]; 252 vertices[i].ReferenceNumber=( Int4)bamggeom->Vertices[i*3+2];252 vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2]; 253 253 vertices[i].DirOfSearch=NoDirOfSearch; 254 254 vertices[i].color =0; … … 297 297 i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing 298 298 i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing 299 edges[i].ref=( Int4)bamggeom->Edges[i*3+2];299 edges[i].ref=(long)bamggeom->Edges[i*3+2]; 300 300 301 301 edges[i].v[0]= vertices + i1; … … 304 304 //get length of edge 305 305 R2 x12=vertices[i2].r-vertices[i1].r; 306 Real8l12=sqrt((x12,x12));306 double l12=sqrt((x12,x12)); 307 307 Hmin=Min(Hmin,l12); 308 308 … … 466 466 //SubDomain 467 467 if(bamggeom->SubDomains){ 468 Int4i0,i1,i2,i3;468 long i0,i1,i2,i3; 469 469 if(verbose>5) printf(" processing SubDomains\n"); 470 470 NbSubDomains=bamggeom->NumSubDomains; … … 493 493 long int verbosity=0; 494 494 495 Int4i,j,k;495 long i,j,k; 496 496 int jj; 497 Int4* head_v=new Int4[nbv];498 Int4* next_p=new Int4[2*nbe];497 long* head_v=new long[nbv]; 498 long* next_p=new long[2*nbe]; 499 499 float* eangle=new float[nbe]; 500 500 double eps=1e-20; … … 635 635 //compute vector of edge i that goes from vertex 0 to vertex 1 636 636 R2 v10=edges[i].v[1]->r - edges[i].v[0]->r; 637 Real8lv10=Norme2(v10);637 double lv10=Norme2(v10); 638 638 //check that its length is not 0 639 639 if(lv10==0) { … … 644 644 //build chains head_v and next_p 645 645 for (j=0;j<2;j++){ 646 Int4v=Number(edges[i].v[j]);646 long v=Number(edges[i].v[j]); 647 647 next_p[k]=head_v[v]; 648 648 head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1; … … 656 656 //exchange vertices position in head_v and next_p till tey are sorted 657 657 while (exch){ 658 Int4*p=head_v+i; // pointer toward head_v[vertex i]659 Int4*po=p; // copy of pointer p660 Int4n=*p; // next value of edge holding i658 long *p=head_v+i; // pointer toward head_v[vertex i] 659 long *po=p; // copy of pointer p 660 long n=*p; // next value of edge holding i 661 661 register float angleold=-1000 ; // angle = - infinity 662 662 ord=0; exch=0; … … 665 665 while (n >=0){ 666 666 ord++; 667 register Int4i1=n/2; // i1 = floor (n/2)668 register Int4j1=n%2; // j1 = 1 if n is odd669 register Int4* pn=next_p+n; // pointer to next_p[n]667 register long i1=n/2; // i1 = floor (n/2) 668 register long j1=n%2; // j1 = 1 if n is odd 669 register long* pn=next_p+n; // pointer to next_p[n] 670 670 671 671 // n = next_p[n] = position in edge of next vertex i … … 695 695 // angular test on current vertex to guess whether it is a corner (ord = number of edges horlding i) 696 696 if(ord == 2) { 697 Int4n1 = head_v[i];698 Int4n2 = next_p[n1];699 Int4i1 = n1/2, i2 = n2/2; // edge number700 Int4j1 = n1%2, j2 = n2%2; // vertex in the edge697 long n1 = head_v[i]; 698 long n2 = next_p[n1]; 699 long i1 = n1/2, i2 = n2/2; // edge number 700 long j1 = n1%2, j2 = n2%2; // vertex in the edge 701 701 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; 702 702 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; … … 718 718 719 719 // close the list around the vertex 720 Int4no=-1, ne = head_v[i];720 long no=-1, ne = head_v[i]; 721 721 while (ne >=0) ne = next_p[no=ne]; 722 722 if(no>=0) next_p[no] = head_v[i]; … … 727 727 for (i=0;i<nbe;i++){ 728 728 for (j=0;j<2;j++){ 729 Int4n1 = next_p[k++];730 Int4i1 = n1/2 ,j1=n1%2;729 long n1 = next_p[k++]; 730 long i1 = n1/2 ,j1=n1%2; 731 731 if( edges[i1].v[j1] != edges[i].v[j]) { 732 732 throw ErrorException(__FUNCT__,exprintf("Bug Adj edge")); … … 743 743 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r; 744 744 //Get length of AB 745 Real8lAB=Norme2(AB);745 double lAB=Norme2(AB); 746 746 //initialize tangent 747 Real8ltg2[2];747 double ltg2[2]; 748 748 ltg2[0]=0;ltg2[1]=0; 749 749 … … 751 751 for (jj=0;jj<2;jj++) { 752 752 R2 tg =edges[i].tg[jj]; 753 Real8ltg=Norme2(tg);753 double ltg=Norme2(tg); 754 754 755 755 //by default, tangent=[0 0] … … 777 777 for (i=0;i<nbe;i++) edges[i].SetUnMark(); 778 778 NbOfCurves = 0; 779 Int4nbgem=0;779 long nbgem=0; 780 780 for (int level=0;level < 2 && nbgem != nbe;level++) 781 781 for (i=0;i<nbe;i++) { … … 920 920 /*}}}1*/ 921 921 /*FUNCTION Geometry::ProjectOnCurve {{{1*/ 922 GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e, Real8s,Vertex &V,VertexOnGeom &GV ) const {923 Real8save_s=s;922 GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,double s,Vertex &V,VertexOnGeom &GV ) const { 923 double save_s=s; 924 924 int NbTry=0; 925 925 retry: … … 937 937 GeometricalEdge *ge[mxe+1]; 938 938 int sensge[mxe+1]; 939 Real8lge[mxe+1];939 double lge[mxe+1]; 940 940 int bge=mxe/2,tge=bge; 941 941 ge[bge] = e.onGeometry; … … 1000 1000 vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1); 1001 1001 1002 Real8sg;1002 double sg; 1003 1003 if (eg0 == eg1) { 1004 register Real8s0= vg0,s1=vg1;1004 register double s0= vg0,s1=vg1; 1005 1005 sg = s0 * (1.0-s) + s * s1; 1006 1006 on=eg0;} 1007 1007 else { 1008 1008 R2 AA=V0,BB; 1009 Real8s0,s1;1009 double s0,s1; 1010 1010 int i=bge; 1011 Real8ll=0;1011 double ll=0; 1012 1012 for(i=bge;i<tge;i++){ 1013 1013 if ( i<0 || i>mxe){ … … 1022 1022 throw ErrorException(__FUNCT__,exprintf("s>1.0")); 1023 1023 } 1024 Real8ls= s*ll;1024 double ls= s*ll; 1025 1025 on =0; 1026 1026 s0 = vg0; 1027 1027 s1= sensge[bge]; 1028 Real8l0=0,l1;1028 double l0=0,l1; 1029 1029 i=bge; 1030 1030 while ( (l1=lge[i]) < ls ) { -
issm/trunk/src/c/Bamgx/objects/Geometry.h
r3232 r3243 25 25 26 26 public: 27 Int4NbRef; // counter of ref on the this class if 0 we can delete28 Int4nbvx,nbtx; // maximum number of vertices29 Int4nbv,nbt,nbiv,nbe; // number of vertices30 Int4NbSubDomains; //31 Int4NbOfCurves;27 long NbRef; // counter of ref on the this class if 0 we can delete 28 long nbvx,nbtx; // maximum number of vertices 29 long nbv,nbt,nbiv,nbe; // number of vertices 30 long NbSubDomains; // 31 long NbOfCurves; 32 32 GeometricalVertex* vertices; 33 33 Triangle* triangles; … … 37 37 Curve* curves; 38 38 R2 pmin,pmax; // extrema 39 Real8coefIcoor; // coef to integer Icoor1;40 Real8MaximalAngleOfCorner;39 double coefIcoor; // coef to integer Icoor1; 40 double MaximalAngleOfCorner; 41 41 42 42 //Constructor/Destructors … … 46 46 47 47 //Operators 48 const GeometricalVertex & operator[] ( Int4i) const { return vertices[i];};49 GeometricalVertex & operator[]( Int4i) { return vertices[i];};50 const GeometricalEdge & operator() ( Int4i) const { return edges[i];};51 GeometricalEdge & operator()( Int4i) { return edges[i];};48 const GeometricalVertex & operator[] (long i) const { return vertices[i];}; 49 GeometricalVertex & operator[](long i) { return vertices[i];}; 50 const GeometricalEdge & operator() (long i) const { return edges[i];}; 51 GeometricalEdge & operator()(long i) { return edges[i];}; 52 52 53 53 //Methods … … 57 57 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) 58 58 ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );} 59 Real8MinimalHmin() {return 2.0/coefIcoor;}60 Real8MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}59 double MinimalHmin() {return 2.0/coefIcoor;} 60 double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 61 61 void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); 62 62 void EmptyGeometry(); … … 64 64 void AfterRead(); 65 65 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();} 66 Int4Number(const GeometricalVertex & t) const { return &t - vertices;}67 Int4Number(const GeometricalVertex * t) const { return t - vertices;}68 Int4Number(const GeometricalEdge & t) const { return &t - edges;}69 Int4Number(const GeometricalEdge * t) const { return t - edges;}70 Int4Number(const Curve * c) const { return c - curves;}71 void UnMarkEdges() {for ( Int4i=0;i<nbe;i++) edges[i].SetUnMark();}72 GeometricalEdge * ProjectOnCurve(const Edge & , Real8,Vertex &,VertexOnGeom &) const ;66 long Number(const GeometricalVertex & t) const { return &t - vertices;} 67 long Number(const GeometricalVertex * t) const { return t - vertices;} 68 long Number(const GeometricalEdge & t) const { return &t - edges;} 69 long Number(const GeometricalEdge * t) const { return t - edges;} 70 long Number(const Curve * c) const { return c - curves;} 71 void UnMarkEdges() {for (long i=0;i<nbe;i++) edges[i].SetUnMark();} 72 GeometricalEdge * ProjectOnCurve(const Edge & ,double,Vertex &,VertexOnGeom &) const ; 73 73 GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const; 74 74 void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); -
issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
r3242 r3243 19 19 long int verbosity=2; 20 20 Icoor2 deta[3], deti,detj; 21 Real8ba[3];21 double ba[3]; 22 22 int nbt =0,ifirst=-1,ilast; 23 23 int i0,i1,i2; … … 38 38 t=tbegin = Bh.FindTriangleContaining(a,deta); 39 39 if( t->det>=0) 40 ilast=NewItem(t, Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det);40 ilast=NewItem(t,double(deta[0])/t->det,double(deta[1])/t->det,double(deta[2])/t->det); 41 41 else 42 42 {// find the nearest boundary edge of the vertex A … … 94 94 if ( det(vi,vj,b)>=0) { 95 95 t=tbegin; 96 Real8ba,bb;96 double ba,bb; 97 97 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); 98 98 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); … … 104 104 i=VerticesOfTriangularEdge[iedge][0]; 105 105 j=VerticesOfTriangularEdge[iedge][1]; 106 Real8dij = detj-deti;106 double dij = detj-deti; 107 107 if (i+j+k != 0 + 1 +2){ 108 108 throw ErrorException(__FUNCT__,exprintf("i+j+k != 0 + 1 +2")); … … 191 191 dt[1]=bamg::det((*t)[2],(*t)[0],b); 192 192 dt[2]=bamg::det((*t)[0],(*t)[1],b); 193 Real8dd = t->det;193 double dd = t->det; 194 194 NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd); 195 195 return ;} … … 216 216 } 217 217 /*}}}1*/ 218 /*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt, Real8 d0,Real8 d1,Real8d2) {{{1*/219 int ListofIntersectionTriangles::NewItem(Triangle * tt, Real8 d0,Real8 d1,Real8d2) {218 /*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) {{{1*/ 219 int ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) { 220 220 register int n; 221 221 R2 x(0,0); … … 257 257 /*}}}1*/ 258 258 /*FUNCTION ListofIntersectionTriangles::Length{{{1*/ 259 Real8ListofIntersectionTriangles::Length(){259 double ListofIntersectionTriangles::Length(){ 260 260 // computation of the length 261 261 … … 275 275 276 276 y = lIntTria[0].x; 277 Real8sxy, s = 0;277 double sxy, s = 0; 278 278 lIntTria[0].s =0; 279 279 SegI->lBegin=s; … … 303 303 /*}}}1*/ 304 304 /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/ 305 Int4 ListofIntersectionTriangles::NewPoints(Vertex* vertices,Int4 &nbv,Int4nbvx){305 long ListofIntersectionTriangles::NewPoints(Vertex* vertices,long &nbv,long nbvx){ 306 306 307 307 //If length<1.5, do nothing 308 Real8s=Length();308 double s=Length(); 309 309 if (s<1.5) return 0; 310 310 311 const Int4nbvold=nbv;311 const long nbvold=nbv; 312 312 int ii = 1 ; 313 313 R2 y,x; 314 314 Metric My,Mx ; 315 Real8sx =0,sy;315 double sx =0,sy; 316 316 int nbi=Max(2,(int) (s+0.5)); 317 Real8sint=s/nbi;318 Real8si =sint;317 double sint=s/nbi; 318 double si =sint; 319 319 320 320 int EndSeg=Size; … … 337 337 sy =lIntTria[ii].s; 338 338 Metric My=lIntTria[ii].m; 339 Real8lxy = sy-sx;340 Real8cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);339 double lxy = sy-sx; 340 double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy); 341 341 342 342 R2 C; 343 Real8cx = 1-cy;343 double cx = 1-cy; 344 344 C = SegI ? SegI->F(si): x * cx + y *cy; 345 345 //C.Echo(); -
issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.h
r3232 r3243 23 23 public: 24 24 Triangle* t; 25 Real8bary[3]; // use if t != 025 double bary[3]; // use if t != 0 26 26 R2 x; 27 27 Metric m; 28 Real8s; // curvilinear coordinate29 Real8sp;// length of the previous segment in m30 Real8sn;// length of the next segment in m28 double s; // curvilinear coordinate 29 double sp;// length of the previous segment in m 30 double sn;// length of the next segment in m 31 31 }; 32 32 … … 34 34 public: 35 35 GeometricalEdge * e; 36 Real8sBegin,sEnd; // abscisse of the seg on edge parameter37 Real8lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length36 double sBegin,sEnd; // abscisse of the seg on edge parameter 37 double lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length 38 38 int last;// last index in ListofIntersectionTriangles for this Sub seg of edge 39 39 40 40 //Methods 41 R2 F( Real8s){42 Real8c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01;41 R2 F(double s){ 42 double c01=lEnd-lBegin, c0=(lEnd-s)/c01, c1=(s-lBegin)/c01; 43 43 if (lBegin>s || s>lEnd){ 44 44 throw ErrorException(__FUNCT__,exprintf("lBegin>s || s>lEnd")); … … 52 52 int MaxSize; 53 53 int Size; 54 Real8len;54 double len; 55 55 int state; 56 56 IntersectionTriangles * lIntTria; … … 77 77 //Methods 78 78 void init(){state=0;len=0;Size=0;} 79 int NewItem(Triangle * tt, Real8 d0,Real8 d1,Real8d2);79 int NewItem(Triangle * tt,double d0,double d1,double d2); 80 80 int NewItem(R2,const Metric & ); 81 81 void SplitEdge(const Triangles & ,const R2 &,const R2 &,int nbegin=0); 82 Real8Length();83 Int4 NewPoints(Vertex *,Int4 & nbv,Int4nbvx);84 void NewSubSeg(GeometricalEdge *e, Real8 s0,Real8s1){82 double Length(); 83 long NewPoints(Vertex *,long & nbv,long nbvx); 84 void NewSubSeg(GeometricalEdge *e,double s0,double s1){ 85 85 long int verbosity=0; 86 86 if (NbSeg>=MaxNbSeg) { -
issm/trunk/src/c/Bamgx/objects/Metric.cpp
r3232 r3243 16 16 17 17 /*Constructor/Destructor*/ 18 /*FUNCTION Metric::Metric(const Real8a[3],const Metric m0, const Metric m1,const Metric m2 ){{{1*/19 Metric::Metric(const Real8a[3],const Metric m0, const Metric m1,const Metric m2 ){18 /*FUNCTION Metric::Metric(const double a[3],const Metric m0, const Metric m1,const Metric m2 ){{{1*/ 19 Metric::Metric(const double a[3],const Metric m0, const Metric m1,const Metric m2 ){ 20 20 Metric mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11, 21 21 a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21, … … 27 27 R2 v2(-v1.y,v1.x); 28 28 29 Real8h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1);30 Real8h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2);29 double h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1); 30 double h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2); 31 31 32 32 vab.lambda1 = 1 / (h1*h1); … … 35 35 } 36 36 /*}}}1*/ 37 /*FUNCTION Metric::Metric( Real8 a,const Metric ma, Real8b,const Metric mb){{{1*/38 Metric::Metric( Real8 a,const Metric ma, Real8b,const Metric mb) {37 /*FUNCTION Metric::Metric( double a,const Metric ma, double b,const Metric mb){{{1*/ 38 Metric::Metric( double a,const Metric ma, double b,const Metric mb) { 39 39 Metric mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22); 40 40 MatVVP2x2 vab(mab); … … 44 44 45 45 46 Real8h1 = a / ma(v1) + b / mb(v1);47 Real8h2 = a / ma(v2) + b / mb(v2);46 double h1 = a / ma(v1) + b / mb(v1); 47 double h2 = a / ma(v2) + b / mb(v2); 48 48 vab.lambda1 = 1 / (h1*h1); 49 49 vab.lambda2 = 1 / (h2*h2); … … 117 117 /*Intermediary*/ 118 118 /*FUNCTION LengthInterpole{{{1*/ 119 Real8LengthInterpole(const Metric Ma,const Metric Mb, R2 AB) {120 Real8k=1./2.;119 double LengthInterpole(const Metric Ma,const Metric Mb, R2 AB) { 120 double k=1./2.; 121 121 int level=0; 122 122 static int kkk=0; 123 123 static Metric Ms1[32],Ms2[32]; 124 static Real8lMs1[32],lMs2[32];124 static double lMs1[32],lMs2[32]; 125 125 static double K[32]; 126 Real8l=0,sss=0;126 double l=0,sss=0; 127 127 Ms1[level]=Ma; 128 128 Ms2[level]=Mb; 129 Real8sa = Ma(AB);130 Real8sb = Mb(AB);129 double sa = Ma(AB); 130 double sb = Mb(AB); 131 131 lMs1[level]=sa; 132 132 lMs2[level]=sb; … … 134 134 level++; 135 135 int i=0; 136 Real8* L= LastMetricInterpole.L, *S = LastMetricInterpole.S;137 Real8sstop = 0.1; // Max(0.6,(sa+sb)/5000);136 double * L= LastMetricInterpole.L, *S = LastMetricInterpole.S; 137 double sstop = 0.1; // Max(0.6,(sa+sb)/5000); 138 138 while (level) { 139 139 level--; … … 141 141 Metric M2=Ms2[level]; 142 142 k=K[level]; 143 Real8s1= lMs1[level];144 Real8s2= lMs2[level];145 146 Real8s= (s1+s2)*k;143 double s1= lMs1[level]; 144 double s2= lMs2[level]; 145 146 double s= (s1+s2)*k; 147 147 if( s > sstop && level < 30 && i < 500-level ) { 148 148 Metric Mi(0.5,M1,0.5,M2); 149 Real8si = Mi(AB);149 double si = Mi(AB); 150 150 if( Abs((s1+s2)-(si+si)) > s1*0.001) 151 151 { … … 286 286 /*}}}1*/ 287 287 /*FUNCTION abscisseInterpole{{{1*/ 288 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8s,int optim) {288 double abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,double s,int optim) { 289 289 if(!optim) LengthInterpole(Ma,Mb,AB); 290 Real8l = s* LastMetricInterpole.lab,r;290 double l = s* LastMetricInterpole.lab,r; 291 291 int j=LastMetricInterpole.opt-1,i=0,k; 292 292 293 Real8* L= LastMetricInterpole.L, *S = LastMetricInterpole.S;293 double * L= LastMetricInterpole.L, *S = LastMetricInterpole.S; 294 294 // warning for optimisation S is the abcisse in [0:0.5] 295 295 // and L is le lenght -
issm/trunk/src/c/Bamgx/objects/Metric.h
r3242 r3243 23 23 public: 24 24 //fields 25 Real8a11,a21,a22;25 double a11,a21,a22; 26 26 //friends 27 27 friend class MatVVP2x2; … … 29 29 Metric(){}; 30 30 Metric(const MatVVP2x2); 31 Metric( Real8a): a11(1/(a*a)),a21(0),a22(1/(a*a)){}32 Metric( Real8 a,Real8 b,Real8c) :a11(a),a21(b),a22(c){}33 Metric( Real8 a,const Metric ma, Real8b,const Metric mb);34 Metric(const Real8a[3],const Metric m0,const Metric m1,const Metric m2 );31 Metric(double a): a11(1/(a*a)),a21(0),a22(1/(a*a)){} 32 Metric(double a,double b,double c) :a11(a),a21(b),a22(c){} 33 Metric( double a,const Metric ma, double b,const Metric mb); 34 Metric(const double a[3],const Metric m0,const Metric m1,const Metric m2 ); 35 35 void Echo(); 36 36 R2 mul(const R2 x)const {return R2(a11*x.x+a21*x.y,a21*x.x+a22*x.y);} 37 Real8det() const {return a11*a22-a21*a21;}37 double det() const {return a11*a22-a21*a21;} 38 38 R2 Orthogonal(const R2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} 39 39 R2 Orthogonal(const I2 x){return R2(-(a21*x.x+a22*x.y),a11*x.x+a21*x.y);} … … 41 41 inline void Box(double &hx,double &hy) const ; 42 42 //operators 43 Metric operator*( Real8 c) const {Real8c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);}44 Metric operator/( Real8 c) const {Real8c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);}43 Metric operator*(double c) const {double c2=c*c;return Metric(a11*c2,a21*c2,a22*c2);} 44 Metric operator/(double c) const {double c2=1/(c*c);return Metric(a11*c2,a21*c2,a22*c2);} 45 45 operator D2xD2(){ return D2xD2(a11,a21,a21,a22);} 46 Real8operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};47 Real8operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};46 double operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);}; 47 double operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;}; 48 48 49 49 }; … … 67 67 void Maxh(double h) {Max(1.0/(h*h));} 68 68 void Isotrope() {lambda1=lambda2=bamg::Max(lambda1,lambda2);} 69 Real8hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));}70 Real8hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));}71 Real8lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);}72 Real8lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);}73 Real8Aniso2() const { return lmax()/lmin();}74 Real8Aniso() const { return sqrt( Aniso2());}75 void BoundAniso(const Real8c){ BoundAniso2(1/(c*c));}76 inline void BoundAniso2(const Real8coef);69 double hmin() const {return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));} 70 double hmax() const {return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));} 71 double lmax() const {return bamg::Max3(lambda1,lambda2,1e-30);} 72 double lmin() const {return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);} 73 double Aniso2() const { return lmax()/lmin();} 74 double Aniso() const { return sqrt( Aniso2());} 75 void BoundAniso(const double c){ BoundAniso2(1/(c*c));} 76 inline void BoundAniso2(const double coef); 77 77 //operators 78 78 void operator *=(double coef){ lambda1*=coef;lambda2*=coef;} … … 80 80 81 81 class SaveMetricInterpole { 82 friend Real8LengthInterpole(const Metric ,const Metric , R2 );83 friend Real8 abscisseInterpole(const Metric ,const Metric , R2 ,Real8,int );82 friend double LengthInterpole(const Metric ,const Metric , R2 ); 83 friend double abscisseInterpole(const Metric ,const Metric , R2 ,double ,int ); 84 84 int opt; 85 Real8lab;86 Real8L[1024],S[1024];85 double lab; 86 double L[1024],S[1024]; 87 87 }; 88 88 … … 90 90 //Functions 91 91 void SimultaneousMatrixReduction( Metric M1, Metric M2,D2xD2 &V); 92 Real8LengthInterpole(const Metric Ma,const Metric Mb, R2 AB);93 Real8 abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,Real8s,int optim=0);92 double LengthInterpole(const Metric Ma,const Metric Mb, R2 AB); 93 double abscisseInterpole(const Metric Ma,const Metric Mb, R2 AB,double s,int optim=0); 94 94 95 95 //inlines 96 inline void MatVVP2x2::BoundAniso2(const Real8coef){96 inline void MatVVP2x2::BoundAniso2(const double coef){ 97 97 if (coef<=1.00000000001){ 98 98 if (lambda1 < lambda2) … … 117 117 } 118 118 inline void Metric::Box(double &hx,double &hy) const { 119 Real8d= a11*a22-a21*a21;119 double d= a11*a22-a21*a21; 120 120 hx = sqrt(a22/d); 121 121 hy = sqrt(a11/d); 122 122 } 123 inline Real8 LengthInterpole(Real8 la,Real8lb) {123 inline double LengthInterpole(double la,double lb) { 124 124 return ( Abs(la - lb) < 1.0e-6*Max3(la,lb,1.0e-20) ) ? (la+lb)/2 : la*lb*log(la/lb)/(la-lb); 125 125 } 126 inline Real8 abscisseInterpole(Real8 la,Real8 lb,Real8 lab,Real8s){126 inline double abscisseInterpole(double la,double lb,double lab,double s){ 127 127 return ( Abs(la - lb) <1.0e-6*Max3(la,lb,1.0e-20)) ? s : (exp(s*lab*(la-lb)/(la*lb))-1)*lb/(la-lb); 128 128 } -
issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
r3233 r3243 94 94 throw ErrorException(__FUNCT__,exprintf("MaxISize <= MaxICoor")); 95 95 } 96 for ( Int4i=0;i<nbv;i++)96 for (long i=0;i<nbv;i++) 97 97 Add(t->vertices[i]); 98 98 } … … 212 212 Icoor1 jj[MaxDeep]; 213 213 register int level=0; // levelevelevel 214 register Int4n0;214 register long n0; 215 215 register QuadTreeBox* b; 216 216 IntQuad h=MaxISize,h0; … … 358 358 // init for optimisation --- 359 359 b = root; 360 register Int4n0;360 register long n0; 361 361 if (!root->n) 362 362 return vn; // empty tree … … 466 466 /*}}}1*/ 467 467 /*FUNCTION QuadTree::ToClose {{{1*/ 468 Vertex * QuadTree::ToClose(Vertex & v, Real8seuil,Icoor1 hx,Icoor1 hy){468 Vertex * QuadTree::ToClose(Vertex & v,double seuil,Icoor1 hx,Icoor1 hy){ 469 469 const Icoor1 i=v.i.x; 470 470 const Icoor1 j=v.i.y; … … 505 505 { 506 506 R2 XY(X,b->v[k]->r); 507 Real8dd;507 double dd; 508 508 if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ){ 509 509 return b->v[k]; -
issm/trunk/src/c/Bamgx/objects/QuadTree.h
r3232 r3243 56 56 Vertex* NearestVertex(Icoor1 i,Icoor1 j); 57 57 Vertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j); 58 Vertex* ToClose(Vertex & , Real8,Icoor1,Icoor1);58 Vertex* ToClose(Vertex & ,double ,Icoor1,Icoor1); 59 59 long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();} 60 60 void Add( Vertex & w); -
issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp
r3233 r3243 9 9 10 10 /*Constructor*/ 11 /*FUNCTION SetOfEdges4::SetOfEdges4( Int4 mmx,Int4nnx){{{1*/12 SetOfEdges4::SetOfEdges4( Int4 mmx,Int4nnx){11 /*FUNCTION SetOfEdges4::SetOfEdges4(long mmx,long nnx){{{1*/ 12 SetOfEdges4::SetOfEdges4(long mmx,long nnx){ 13 13 14 14 /*Intermediary*/ … … 19 19 nbax =mmx; // 3 * number of triangles 20 20 NbOfEdges=0; 21 head = new Int4[nx];21 head = new long [nx]; 22 22 Edges= new Int4Edge[nbax]; 23 23 … … 30 30 /*Methods*/ 31 31 /*FUNCTION SetOfEdges4::find {{{1*/ 32 Int4 SetOfEdges4::find(Int4 ii,Int4jj) {32 long SetOfEdges4::find(long ii,long jj) { 33 33 34 34 /*Intermediary*/ … … 58 58 /*}}}1*/ 59 59 /*FUNCTION SetOfEdges4::add{{{1*/ 60 Int4 SetOfEdges4::add(Int4 ii,Int4jj) {60 long SetOfEdges4::add(long ii,long jj) { 61 61 62 62 /*Intermediary*/ -
issm/trunk/src/c/Bamgx/objects/SetOfE4.h
r3232 r3243 9 9 friend class SetOfEdges4; 10 10 public: 11 Int4i,j;12 Int4next;11 long i,j; 12 long next; 13 13 }; 14 14 … … 16 16 17 17 private: 18 Int4nx,nbax,NbOfEdges;19 Int4* head;18 long nx,nbax,NbOfEdges; 19 long* head; 20 20 Int4Edge* Edges; 21 21 22 22 public: 23 SetOfEdges4( Int4 ,Int4);// nb Edges mx , nb de sommet23 SetOfEdges4(long ,long);// nb Edges mx , nb de sommet 24 24 ~SetOfEdges4() {delete [] head; delete [] Edges;} 25 Int4 add (Int4 ii,Int4jj);26 Int4 SortAndAdd (Int4 ii,Int4jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;}27 Int4nb(){return NbOfEdges;}28 Int4 find (Int4 ii,Int4jj);29 Int4 SortAndFind (Int4 ii,Int4jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;}30 Int4 i(Int4k){return Edges[k].i;}31 Int4 j(Int4k){return Edges[k].j;}32 Int4 newarete(Int4k){return NbOfEdges == k+1;}25 long add (long ii,long jj); 26 long SortAndAdd (long ii,long jj) {return ii <=jj ? add (ii,jj) : add (jj,ii) ;} 27 long nb(){return NbOfEdges;} 28 long find (long ii,long jj); 29 long SortAndFind (long ii,long jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} 30 long i(long k){return Edges[k].i;} 31 long j(long k){return Edges[k].j;} 32 long newarete(long k){return NbOfEdges == k+1;} 33 33 34 34 //operators 35 Int4Edge & operator[]( Int4k){return Edges[k];}35 Int4Edge & operator[](long k){return Edges[k];} 36 36 }; 37 37 } -
issm/trunk/src/c/Bamgx/objects/SubDomain.h
r3232 r3243 19 19 public: 20 20 Triangle * head; 21 Int4ref;21 long ref; 22 22 int sens; // -1 or 1 23 23 Edge* edge; // to geometrical 24 24 25 25 //Inline methods 26 inline void Set(const Triangles &, Int4,Triangles &);26 inline void Set(const Triangles &,long,Triangles &); 27 27 }; 28 28 -
issm/trunk/src/c/Bamgx/objects/Triangle.cpp
r3242 r3243 92 92 /*}}}*/ 93 93 /*FUNCTION Triangle::Optim{{{1*/ 94 Int4Triangle::Optim(short i,int koption) {94 long Triangle::Optim(short i,int koption) { 95 95 // turn around (positive direction) 96 96 Triangle *t=this; 97 Int4NbSwap =0;97 long NbSwap =0; 98 98 int k = 0; 99 99 int j = OppositeEdge[i]; … … 183 183 else { 184 184 // critere de Delaunay anisotrope 185 Real8som;185 double som; 186 186 I2 AB=(I2) *sb - (I2) *sa; 187 187 I2 MAB2=((I2) *sb + (I2) *sa); -
issm/trunk/src/c/Bamgx/objects/Triangle.h
r3242 r3243 29 29 union { 30 30 Triangle * link ; 31 Int4color;31 long color; 32 32 }; 33 33 34 34 //Constructors/Destructors 35 35 Triangle() {} 36 Triangle(Triangles *Th, Int4 i,Int4 j,Int4k);36 Triangle(Triangles *Th,long i,long j,long k); 37 37 Triangle(Vertex *v0,Vertex *v1,Vertex *v2); 38 38 … … 46 46 void Echo(); 47 47 int swap(short a1,int=0); 48 Int4Optim(short a,int =0);48 long Optim(short a,int =0); 49 49 int Locked(int a)const { return TriaAdjSharedEdge[a]&4;} 50 50 int Hidden(int a)const { return TriaAdjSharedEdge[a]&16;} … … 58 58 Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];} 59 59 Triangle* Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const ; 60 void ReNumbering(Triangle *tb,Triangle *te, Int4*renu){60 void ReNumbering(Triangle *tb,Triangle *te, long *renu){ 61 61 if (link >=tb && link <te) link = tb + renu[link -tb]; 62 62 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb]; … … 64 64 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb]; 65 65 } 66 void ReNumbering(Vertex *vb,Vertex *ve, Int4*renu){66 void ReNumbering(Vertex *vb,Vertex *ve, long *renu){ 67 67 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb]; 68 68 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb]; -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3242 r3243 42 42 int kt=0; 43 43 int * kk = new int [Tho.nbv]; 44 Int4 * reft = new Int4[Tho.nbt];45 Int4nbInT = Tho.TriangleReferenceList(reft);46 Int4 * refv = new Int4[Tho.nbv];44 long * reft = new long[Tho.nbt]; 45 long nbInT = Tho.TriangleReferenceList(reft); 46 long * refv = new long[Tho.nbv]; 47 47 48 48 for (i=0;i<Tho.nbv;i++) … … 84 84 printf(" number of triangles %i, remove = \n",kt,nbInT-kt); 85 85 printf(" number of New boundary edge %i\n",nbNewBedge); 86 Int4inbvx =k;86 long inbvx =k; 87 87 PreInit(inbvx); 88 88 for (i=0;i<Tho.nbv;i++) … … 138 138 } 139 139 /*}}}1*/ 140 /*FUNCTION Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth, Int4nbvxx) COPY{{{1*/141 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth, Int4nbvxx)140 /*FUNCTION Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,long nbvxx) COPY{{{1*/ 141 Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,long nbvxx) 142 142 : Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this)) { 143 143 Gh.NbRef++; 144 144 nbvxx = Max(nbvxx,Th.nbv); 145 Int4i;145 long i; 146 146 // do all the allocation to be sure all the pointer existe 147 147 … … 234 234 void Triangles::ReadMesh(double* index,double* x,double* y,int nods,int nels){ 235 235 236 Real8Hmin = HUGE_VAL;// the infinie value237 Int4i1,i2,i3,iref;238 Int4i,j;239 Int4hvertices =0;236 double Hmin = HUGE_VAL;// the infinie value 237 long i1,i2,i3,iref; 238 long i,j; 239 long hvertices =0; 240 240 Metric M1(1); 241 241 int Version=1,dim=2; … … 266 266 for (i=0;i<nbt;i++){ 267 267 Triangle & t = triangles[i]; 268 i1=( Int4)index[i*3+0]-1; //for C indexing269 i2=( Int4)index[i*3+1]-1; //for C indexing270 i3=( Int4)index[i*3+2]-1; //for C indexing268 i1=(long)index[i*3+0]-1; //for C indexing 269 i2=(long)index[i*3+1]-1; //for C indexing 270 i3=(long)index[i*3+2]-1; //for C indexing 271 271 t=Triangle(this,i1,i2,i3); 272 272 t.color=i; … … 284 284 285 285 int verbose; 286 Real8Hmin = HUGE_VAL;// the infinie value287 Int4i1,i2,i3,iref;288 Int4i,j;289 Int4hvertices =0;290 Int4ifgeom=0;286 double Hmin = HUGE_VAL;// the infinie value 287 long i1,i2,i3,iref; 288 long i,j; 289 long hvertices =0; 290 long ifgeom=0; 291 291 Metric M1(1); 292 292 int Version=1,dim=2; … … 312 312 vertices[i].DirOfSearch =NoDirOfSearch; 313 313 vertices[i].m=M1; 314 vertices[i].color=( Int4)bamgmesh->Vertices[i*3+2];314 vertices[i].color=(long)bamgmesh->Vertices[i*3+2]; 315 315 } 316 316 nbtx=2*nbvx-2; // for filling The Holes and quadrilaterals … … 327 327 for (i=0;i<nbt;i++){ 328 328 Triangle & t = triangles[i]; 329 i1=( Int4)bamgmesh->Triangles[i*4+0]-1; //for C indexing330 i2=( Int4)bamgmesh->Triangles[i*4+1]-1; //for C indexing331 i3=( Int4)bamgmesh->Triangles[i*4+2]-1; //for C indexing329 i1=(long)bamgmesh->Triangles[i*4+0]-1; //for C indexing 330 i2=(long)bamgmesh->Triangles[i*4+1]-1; //for C indexing 331 i3=(long)bamgmesh->Triangles[i*4+2]-1; //for C indexing 332 332 t=Triangle(this,i1,i2,i3); 333 t.color=( Int4)bamgmesh->Triangles[i*4+3];333 t.color=(long)bamgmesh->Triangles[i*4+3]; 334 334 } 335 335 } … … 341 341 if(bamgmesh->Quadrilaterals){ 342 342 if(verbose>5) printf(" processing Quadrilaterals\n"); 343 Int4i1,i2,i3,i4,iref;343 long i1,i2,i3,i4,iref; 344 344 triangles =new Triangle[nbt]; 345 345 for (i=0;i<bamgmesh->NumQuadrilaterals;i++){ … … 347 347 Triangle & t1 = triangles[2*i]; 348 348 Triangle & t2 = triangles[2*i+1]; 349 i1=( Int4)bamgmesh->Quadrilaterals[i*4+0]-1; //for C indexing350 i2=( Int4)bamgmesh->Quadrilaterals[i*4+1]-1; //for C indexing351 i3=( Int4)bamgmesh->Quadrilaterals[i*4+2]-1; //for C indexing352 i4=( Int4)bamgmesh->Quadrilaterals[i*4+3]-1; //for C indexing349 i1=(long)bamgmesh->Quadrilaterals[i*4+0]-1; //for C indexing 350 i2=(long)bamgmesh->Quadrilaterals[i*4+1]-1; //for C indexing 351 i3=(long)bamgmesh->Quadrilaterals[i*4+2]-1; //for C indexing 352 i4=(long)bamgmesh->Quadrilaterals[i*4+3]-1; //for C indexing 353 353 t1=Triangle(this,i1,i2,i3); 354 354 t2=Triangle(this,i3,i4,i1); … … 383 383 VerticesOnGeomEdge= new VertexOnGeom[NbVerticesOnGeomEdge] ; 384 384 for (i=0;i<NbVerticesOnGeomEdge;i++){ 385 Int4i1,i2;386 Real8s;387 i1=( Int4)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing388 i2=( Int4)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing389 s =( Int4)bamgmesh->VerticesOnGeometricEdge[i*3+2];385 long i1,i2; 386 double s; 387 i1=(long)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing 388 i2=(long)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing 389 s =(long)bamgmesh->VerticesOnGeometricEdge[i*3+2]; 390 390 VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s); 391 391 } … … 419 419 edges[i].adj[1]=0; 420 420 R2 x12 = vertices[i2].r-vertices[i1].r; 421 Real8l12=sqrt( (x12,x12));421 double l12=sqrt( (x12,x12)); 422 422 423 423 if (!hvertices) { … … 447 447 for (j=0;j<2;j++) { 448 448 Vertex *v=edges[i].v[j]; 449 Int4i0=v->color,j0;449 long i0=v->color,j0; 450 450 if(i0==-1){ 451 451 v->color=i*2+j; … … 503 503 //SubDomain 504 504 if(bamgmesh->SubDomains){ 505 Int4i3,head,sens;505 long i3,head,sens; 506 506 if(verbose>5) printf(" processing SubDomains\n"); 507 507 NbSubDomains=bamgmesh->NumSubDomains; … … 548 548 549 549 //Build reft that holds the number the subdomain number of each triangle 550 Int4* reft = new Int4[nbt];551 Int4nbInT = TriangleReferenceList(reft);550 long* reft = new long[nbt]; 551 long nbInT = TriangleReferenceList(reft); 552 552 553 553 //Vertices … … 810 810 /*}}}1*/ 811 811 /*FUNCTION Triangles::ReadMetric{{{1*/ 812 void Triangles::ReadMetric(BamgOpts* bamgopts,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8coef=1) {812 void Triangles::ReadMetric(BamgOpts* bamgopts,const double hmin1=1.0e-30,const double hmax1=1.0e30,const double coef=1) { 813 813 int i,j; 814 814 815 815 if(bamgopts->verbose>3) printf(" processing metric\n"); 816 816 817 Real8hmin = Max(hmin1,MinimalHmin());818 Real8hmax = Min(hmax1,MaximalHmax());817 double hmin = Max(hmin1,MinimalHmin()); 818 double hmax = Min(hmax1,MaximalHmax()); 819 819 820 820 //for now we only use j==3 … … 822 822 823 823 for (i=0;i<nbv;i++){ 824 Real8h;824 double h; 825 825 if (j == 1){ 826 826 h=bamgopts->metric[i]; … … 830 830 //do not erase metric computed by hVertices 831 831 if (vertices[i].m.a11==1 && vertices[i].m.a21==0 && vertices[i].m.a22==1){ 832 Real8a,b,c;832 double a,b,c; 833 833 a=bamgopts->metric[i*3+0]; 834 834 b=bamgopts->metric[i*3+1]; … … 867 867 double errg =bamgopts->errg; 868 868 869 Real8ss[2]={0.00001,0.99999};870 Real8errC = 2*sqrt(2*errg);871 Real8hmax = Gh.MaximalHmax();872 Real8hmin = Gh.MinimalHmin();869 double ss[2]={0.00001,0.99999}; 870 double errC = 2*sqrt(2*errg); 871 double hmax = Gh.MaximalHmax(); 872 double hmin = Gh.MinimalHmin(); 873 873 874 874 //check that hmax is positive … … 884 884 885 885 //loop over all the vertices on edges 886 for ( Int4i=0;i<nbe;i++){886 for (long i=0;i<nbe;i++){ 887 887 for (int j=0;j<2;j++){ 888 888 … … 892 892 893 893 GeometricalEdge * eg = GV; 894 Real8s = GV;894 double s = GV; 895 895 R2 tg; 896 Real8R1= eg->R1tg(s,tg);897 Real8ht=hmax;896 double R1= eg->R1tg(s,tg); 897 double ht=hmax; 898 898 // err relative to the length of the edge 899 899 if (R1>1.0e-20) { 900 900 ht = Min(Max(errC/R1,hmin),hmax); 901 901 } 902 Real8hn=Min(hmax,ht*anisomax);902 double hn=Min(hmax,ht*anisomax); 903 903 if (ht<=0 || hn<=0){ 904 904 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); … … 1072 1072 /*}}}1*/ 1073 1073 /*FUNCTION Triangles::BoundAnisotropy{{{1*/ 1074 void Triangles::BoundAnisotropy( Real8 anisomax,Real8hminaniso) {1074 void Triangles::BoundAnisotropy(double anisomax,double hminaniso) { 1075 1075 1076 1076 long int verbosity=0; … … 1080 1080 if (verbosity > 1) printf(" BoundAnisotropy by %g\n",anisomax); 1081 1081 1082 Real8h1=1.e30,h2=1e-30;1083 Real8coef = 1./(anisomax*anisomax);1084 Real8hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0;1082 double h1=1.e30,h2=1e-30; 1083 double coef = 1./(anisomax*anisomax); 1084 double hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0; 1085 1085 1086 1086 //loop over all vertices 1087 for ( Int4i=0;i<nbv;i++){1087 for (long i=0;i<nbv;i++){ 1088 1088 MatVVP2x2 Vp(vertices[i]); 1089 1089 double lmax=Vp.lmax(); … … 1140 1140 //initialize st and edge4 1141 1141 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 1142 Int4* st = new Int4[nbt*3];1142 long* st = new long[nbt*3]; 1143 1143 1144 1144 //initialize st as -1 (chaining algorithm) … … 1154 1154 } 1155 1155 //keep nbe in nbeold 1156 Int4nbeold = nbe;1156 long nbeold = nbe; 1157 1157 1158 1158 //Go through the triangles and ass the edges in edge4 if they are not there yet … … 1161 1161 for (j=0;j<3;j++) { 1162 1162 //Add Edge to edge4 (k=numberofedges in edge4) 1163 Int4k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]]));1164 Int4invisible = triangles[i].Hidden(j);1163 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 1164 long invisible = triangles[i].Hidden(j); 1165 1165 1166 1166 //if st[k] has not been changed yet, add 3*i+j (= vertex position in the index) … … 1195 1195 1196 1196 //delete edge4 1197 Int4nbedges = edge4->nb(); // the total number of edges1197 long nbedges = edge4->nb(); // the total number of edges 1198 1198 delete edge4; edge4=NULL; 1199 1199 … … 1239 1239 1240 1240 for (i=0;i<nbedges;i++){ 1241 Int4add= -1;1241 long add= -1; 1242 1242 1243 1243 //internal edge (belongs to two triangles) … … 1300 1300 Vertex* v=edges[i].v[j]; 1301 1301 //get vertex color (i0) 1302 Int4i0=v->color,j0;1302 long i0=v->color,j0; 1303 1303 1304 1304 //if color<0 (first time), no adjacent edge … … 1338 1338 1339 1339 //color the subdomains 1340 Int4* colorT= new Int4[nbt];1340 long* colorT= new long[nbt]; 1341 1341 Triangle *tt,*t; 1342 1342 … … 1380 1380 1381 1381 //build subdomains 1382 Int4isd;1382 long isd; 1383 1383 subdomains = new SubDomain[NbSubDomains]; 1384 1384 … … 1410 1410 1411 1411 //build colorV -1 for all vertex and 0 for the vertices belonging to edges 1412 Int4* colorV = new Int4[nbv];1412 long* colorV = new long[nbv]; 1413 1413 for (i=0;i<nbv;i++) colorV[i]=-1; 1414 1414 for (i=0;i<nbe;i++){ … … 1472 1472 //initialize edge4 again 1473 1473 edge4= new SetOfEdges4(nbe,nbv); 1474 Real8hmin = HUGE_VAL;1474 double hmin = HUGE_VAL; 1475 1475 int kreq=0; 1476 1476 for (i=0;i<nbe;i++){ 1477 1477 1478 Int4i0 = Number(edges[i][0]);1479 Int4i1 = Number(edges[i][1]);1480 Int4j0 = colorV[i0];1481 Int4j1 = colorV[i1];1478 long i0 = Number(edges[i][0]); 1479 long i1 = Number(edges[i][1]); 1480 long j0 = colorV[i0]; 1481 long j1 = colorV[i1]; 1482 1482 1483 1483 Gh.edges[i].v[0] = Gh.vertices + j0; … … 1499 1499 1500 1500 R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; 1501 Real8l12=Norme2(x12);1501 double l12=Norme2(x12); 1502 1502 hmin = Min(hmin,l12); 1503 1503 … … 1530 1530 it = Number(subdomains[i].head); 1531 1531 j = subdomains[i].sens; 1532 Int4i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]);1533 Int4i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]);1532 long i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]); 1533 long i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]); 1534 1534 k = edge4->SortAndFind(i0,i1); 1535 1535 if(k>=0){ … … 1562 1562 const int dim = 2; 1563 1563 double* s=NULL; 1564 Int4nbsol;1564 long nbsol; 1565 1565 int verbosity; 1566 1566 … … 1578 1578 double* ss=(double*)s; 1579 1579 double sA,sB,sC; 1580 Real8* detT = new Real8[nbt];1581 Real8* sumareas = new Real8[nbv];1582 Real8* alpha= new Real8[nbt*3];1583 Real8* beta = new Real8[nbt*3];1584 Real8* dx_elem = new Real8[nbt];1585 Real8* dy_elem = new Real8[nbt];1586 Real8* dx_vertex = new Real8[nbv];1587 Real8* dy_vertex = new Real8[nbv];1588 Real8* dxdx_elem = new Real8[nbt];1589 Real8* dxdy_elem = new Real8[nbt];1590 Real8* dydy_elem = new Real8[nbt];1591 Real8* dxdx_vertex= new Real8[nbv];1592 Real8* dxdy_vertex= new Real8[nbv];1593 Real8* dydy_vertex= new Real8[nbv];1580 double* detT = new double[nbt]; 1581 double* sumareas = new double[nbv]; 1582 double* alpha= new double[nbt*3]; 1583 double* beta = new double[nbt*3]; 1584 double* dx_elem = new double[nbt]; 1585 double* dy_elem = new double[nbt]; 1586 double* dx_vertex = new double[nbv]; 1587 double* dy_vertex = new double[nbv]; 1588 double* dxdx_elem = new double[nbt]; 1589 double* dxdy_elem = new double[nbt]; 1590 double* dydy_elem = new double[nbt]; 1591 double* dxdx_vertex= new double[nbv]; 1592 double* dxdy_vertex= new double[nbv]; 1593 double* dydy_vertex= new double[nbv]; 1594 1594 1595 1595 //display infos … … 1623 1623 1624 1624 //compute triangle determinant (2*Area) 1625 Real8dett = bamg::Area2(A,B,C);1625 double dett = bamg::Area2(A,B,C); 1626 1626 detT[i]=dett; 1627 1627 … … 1655 1655 1656 1656 //for all Solutions 1657 for ( Int4nusol=0;nusol<nbsol;nusol++) {1658 Real8smin=ss[nusol],smax=ss[nusol];1657 for (long nusol=0;nusol<nbsol;nusol++) { 1658 double smin=ss[nusol],smax=ss[nusol]; 1659 1659 1660 1660 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) … … 1663 1663 smax=Max(smax,ss[iv*nbsol+nusol]); 1664 1664 } 1665 Real8sdelta=smax-smin;1666 Real8absmax=Max(Abs(smin),Abs(smax));1665 double sdelta=smax-smin; 1666 double absmax=Max(Abs(smin),Abs(smax)); 1667 1667 1668 1668 //display info … … 1701 1701 for(p=head_s[i];p!=-1;p=next_p[p]){ 1702 1702 //Get triangle number 1703 k=( Int4)(p/3);1703 k=(long)(p/3); 1704 1704 dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i]; 1705 1705 dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i]; … … 1726 1726 for(p=head_s[i];p!=-1;p=next_p[p]){ 1727 1727 //Get triangle number 1728 k=( Int4)(p/3);1728 k=(long)(p/3); 1729 1729 dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i]; 1730 1730 dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i]; … … 1765 1765 const int dim = 2; 1766 1766 double* s=NULL; 1767 Int4nbsol;1767 long nbsol; 1768 1768 int NbJacobi; 1769 1769 int verbosity; … … 1778 1778 1779 1779 //initialization of some variables 1780 Int4i,k,iA,iB,iC,iv;1780 long i,k,iA,iB,iC,iv; 1781 1781 R2 O(0,0); 1782 1782 double* ss=(double*)s; 1783 1783 double sA,sB,sC; 1784 Real8* detT = new Real8[nbt];1785 Real8* Mmass= new Real8[nbv];1786 Real8* Mmassxx= new Real8[nbv];1787 Real8* dxdx= new Real8[nbv];1788 Real8* dxdy= new Real8[nbv];1789 Real8* dydy= new Real8[nbv];1790 Real8* workT= new Real8[nbt];1791 Real8* workV= new Real8[nbv];1784 double* detT = new double[nbt]; 1785 double* Mmass= new double[nbv]; 1786 double* Mmassxx= new double[nbv]; 1787 double* dxdx= new double[nbv]; 1788 double* dxdy= new double[nbv]; 1789 double* dydy= new double[nbv]; 1790 double* workT= new double[nbt]; 1791 double* workV= new double[nbv]; 1792 1792 int* OnBoundary = new int[nbv]; 1793 1793 … … 1824 1824 1825 1825 //compute triangle determinant (2*Area) 1826 Real8dett = bamg::Area2(A,B,C);1826 double dett = bamg::Area2(A,B,C); 1827 1827 detT[i]=dett; 1828 1828 dett /= 6; … … 1861 1861 1862 1862 //for all Solution 1863 for ( Int4nusol=0;nusol<nbsol;nusol++) {1864 1865 Real8smin=ss[nusol],smax=ss[nusol];1866 Real8h1=1.e30,h2=1e-30,rx=0;1867 Real8hn1=1.e30,hn2=1e-30,rnx =1.e-30;1863 for (long nusol=0;nusol<nbsol;nusol++) { 1864 1865 double smin=ss[nusol],smax=ss[nusol]; 1866 double h1=1.e30,h2=1e-30,rx=0; 1867 double hn1=1.e30,hn2=1e-30,rnx =1.e-30; 1868 1868 1869 1869 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) … … 1873 1873 smax=Max(smax,ss[iv*nbsol+nusol]); 1874 1874 } 1875 Real8sdelta=smax-smin;1876 Real8absmax=Max(Abs(smin),Abs(smax));1875 double sdelta=smax-smin; 1876 double absmax=Max(Abs(smin),Abs(smax)); 1877 1877 1878 1878 //display info … … 1967 1967 } 1968 1968 1969 Int4kk=0;1969 long kk=0; 1970 1970 for ( iv=0,k=0 ; iv<nbv; iv++){ 1971 1971 if(Mmassxx[iv]>0){ … … 1988 1988 // correction of second derivative 1989 1989 // by a laplacien 1990 Real8* d2[3] = {dxdx, dxdy, dydy};1991 Real8* dd;1990 double* d2[3] = {dxdx, dxdy, dydy}; 1991 double* dd; 1992 1992 for (int xy = 0;xy<3;xy++) { 1993 1993 dd = d2[xy]; … … 2000 2000 iB = Number(triangles[i][1]); 2001 2001 iC = Number(triangles[i][2]); 2002 Real8cc=3;2002 double cc=3; 2003 2003 if(ijacobi==0) 2004 cc = Max(( Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);2004 cc = Max((double) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 2005 2005 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 2006 2006 } … … 2013 2013 iB = Number(triangles[i][1]); 2014 2014 iC = Number(triangles[i][2]); 2015 Real8cc = workT[i]*detT[i];2015 double cc = workT[i]*detT[i]; 2016 2016 workV[iA] += cc; 2017 2017 workV[iB] += cc; … … 2157 2157 if (nbnewv) 2158 2158 { // 2159 Int4n = nbnewv+NbVerticesOnGeomVertex;2160 Int4i,j,k;2159 long n = nbnewv+NbVerticesOnGeomVertex; 2160 long i,j,k; 2161 2161 VertexOnGeom * vog = new VertexOnGeom[n]; 2162 2162 for ( i =0; i<NbVerticesOnGeomVertex;i++) … … 2175 2175 if ( v >= LastOld) 2176 2176 { // a new vertex 2177 Int4old = v->ReferenceNumber ; // the old same vertex2178 Int4i = ( v - LastOld);2177 long old = v->ReferenceNumber ; // the old same vertex 2178 long i = ( v - LastOld); 2179 2179 // if the old is on vertex => warning 2180 2180 // else the old is on edge => ok … … 2208 2208 2209 2209 //check that there is no triangle with 0 determinant 2210 for ( Int4t = 0; t < nbt; t++){2210 for (long t = 0; t < nbt; t++){ 2211 2211 if (!triangles[t].det) k++; 2212 2212 } … … 2217 2217 //Force Edges 2218 2218 TriangleAdjacent ta(0,0); 2219 for ( Int4i = 0; i < nbe; i++){2219 for (long i = 0; i < nbe; i++){ 2220 2220 2221 2221 //Force edge i … … 2226 2226 if (nbswp) nbfe++; 2227 2227 if ( nbswp < 0 && k < 5){ 2228 for ( Int4j = 0; j < nbe; j++){2228 for (long j = 0; j < nbe; j++){ 2229 2229 printf("Edge %i: %i %i\n",j,Number(edges[j][0]),Number(edges[j][1])); 2230 2230 } … … 2238 2238 throw ErrorException(__FUNCT__,exprintf("There are %i lost edges, the boundary might be crossing",k)); 2239 2239 } 2240 for ( Int4j=0;j<nbv;j++){2240 for (long j=0;j<nbv;j++){ 2241 2241 Nbswap += vertices[j].Optim(1,0); 2242 2242 } … … 2255 2255 Triangle ** HeapTriangle = new Triangle* [nbt]; 2256 2256 Triangle *t,*t1; 2257 Int4k,it;2258 2259 for ( Int4itt=0;itt<nbt;itt++)2257 long k,it; 2258 2259 for (long itt=0;itt<nbt;itt++) 2260 2260 triangles[itt].link=0; // par defaut pas de couleur 2261 2261 2262 Int4NbSubDomTot =0;2262 long NbSubDomTot =0; 2263 2263 for ( it=0;it<nbt;it++) { 2264 2264 if ( ! triangles[it].link ) { 2265 2265 t = triangles + it; 2266 2266 NbSubDomTot++;; // new composante connexe 2267 Int4i = 0; // niveau de la pile2267 long i = 0; // niveau de la pile 2268 2268 t->link = t ; // sd forme d'un triangle cicular link 2269 2269 … … 2324 2324 if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains ) 2325 2325 { // No geom sub domain 2326 Int4i;2326 long i; 2327 2327 if (subdomains) delete [] subdomains; 2328 2328 subdomains = new SubDomain[ NbSubDomTot]; … … 2332 2332 subdomains[i].ref=i+1; 2333 2333 } 2334 Int4 * mark = new Int4[nbt];2334 long * mark = new long[nbt]; 2335 2335 for (it=0;it<nbt;it++) 2336 2336 mark[it]=triangles[it].link ? -1 : -2; … … 2359 2359 // because in this case we have only the true boundary edge 2360 2360 // so teh boundary is manifold 2361 Int4nbk = NbSubDomains;2361 long nbk = NbSubDomains; 2362 2362 while (nbk) 2363 2363 for (it=0;it<nbt && nbk ;it++) … … 2365 2365 { 2366 2366 Triangle *ta = triangles[it].TriangleAdj(na); 2367 Int4kl = ta ? mark[Number(ta)] : -2;2368 Int4kr = mark[it];2367 long kl = ta ? mark[Number(ta)] : -2; 2368 long kr = mark[it]; 2369 2369 if(kr !=kl) { 2370 2370 if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0) … … 2378 2378 } 2379 2379 } 2380 Int4j=0;2380 long j=0; 2381 2381 for ( i=0;i<NbSubDomains;i++) 2382 2382 if((-subdomains[i].ref) %2) { // good … … 2407 2407 subdomains = new SubDomain[ Gh.NbSubDomains]; 2408 2408 NbSubDomains =Gh.NbSubDomains; 2409 Int4err=0;2409 long err=0; 2410 2410 ReMakeTriangleContainingTheVertex(); 2411 Int4 * mark = new Int4[nbt];2411 long * mark = new long[nbt]; 2412 2412 Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge(); 2413 2413 2414 2414 for (it=0;it<nbt;it++) 2415 2415 mark[it]=triangles[it].link ? -1 : -2; 2416 Int4inew =0;2417 for ( Int4i=0;i<NbSubDomains;i++) {2416 long inew =0; 2417 for (long i=0;i<NbSubDomains;i++) { 2418 2418 GeometricalEdge &eg = *Gh.subdomains[i].edge; 2419 2419 subdomains[i].ref = Gh.subdomains[i].ref; … … 2450 2450 throw ErrorException(__FUNCT__,exprintf("bad definition of SubSomain %i",i)); 2451 2451 } 2452 Int4it = Number(t);2452 long it = Number(t); 2453 2453 if (mark[it] >=0) { 2454 2454 break; … … 2458 2458 inew++; 2459 2459 Triangle *tt=t; 2460 Int4kkk=0;2460 long kkk=0; 2461 2461 do 2462 2462 { … … 2595 2595 2596 2596 // find extrema coordinates of vertices pmin,pmax 2597 Int4i;2597 long i; 2598 2598 if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv); 2599 2599 … … 2611 2611 2612 2612 //initialize st 2613 Int4* st = new Int4[nbt*3];2613 long* st = new long[nbt*3]; 2614 2614 for (i=0;i<nbt*3;i++) st[i]=-1; 2615 2615 2616 2616 //check number of edges 2617 Int4kk =0;2617 long kk =0; 2618 2618 for (i=0;i<nbe;i++){ 2619 2619 kk=kk+(i == edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); … … 2626 2626 for (i=0;i<nbt;i++){ 2627 2627 for (int j=0;j<3;j++) { 2628 Int4k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),2628 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), 2629 2629 Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 2630 Int4invisible = triangles[i].Hidden(j);2630 long invisible = triangles[i].Hidden(j); 2631 2631 if(st[k]==-1){ 2632 2632 st[k]=3*i+j; … … 2660 2660 2661 2661 // check the consistant of edge[].adj and the geometrical required vertex 2662 Int4k=0;2662 long k=0; 2663 2663 for (i=0;i<edge4->nb();i++){ 2664 2664 if (st[i] >=0){ // edge alone 2665 2665 if (i < nbe) { 2666 Int4i0=edge4->i(i);2666 long i0=edge4->i(i); 2667 2667 ordre[i0] = vertices+i0; 2668 Int4i1=edge4->j(i);2668 long i1=edge4->j(i); 2669 2669 ordre[i1] = vertices+i1; 2670 2670 } … … 2682 2682 2683 2683 /* mesh generation with boundary points*/ 2684 Int4nbvb = 0;2684 long nbvb = 0; 2685 2685 for (i=0;i<nbv;i++){ 2686 2686 vertices[i].t=0; … … 2692 2692 2693 2693 Triangle* savetriangles= triangles; 2694 Int4savenbt=nbt;2695 Int4savenbtx=nbtx;2694 long savenbt=nbt; 2695 long savenbtx=nbtx; 2696 2696 SubDomain * savesubdomains = subdomains; 2697 2697 subdomains = 0; 2698 2698 2699 Int4Nbtriafillhole = 2*nbvb;2699 long Nbtriafillhole = 2*nbvb; 2700 2700 Triangle* triafillhole =new Triangle[Nbtriafillhole]; 2701 2701 triangles = triafillhole; … … 2744 2744 2745 2745 // We add the vertices one by one 2746 Int4NbSwap=0;2747 for ( Int4icount=2; icount<nbvb; icount++) {2746 long NbSwap=0; 2747 for (long icount=2; icount<nbvb; icount++) { 2748 2748 Vertex *vi = ordre[icount]; 2749 2749 Icoor2 dete[3]; … … 2756 2756 // inforce the boundary 2757 2757 TriangleAdjacent ta(0,0); 2758 Int4nbloss = 0,knbe=0;2758 long nbloss = 0,knbe=0; 2759 2759 for ( i = 0; i < nbe; i++){ 2760 2760 if (st[i] >=0){ // edge alone => on border ... FH oct 2009 … … 2775 2775 // remove all the hole 2776 2776 // remove all the good sub domain 2777 Int4krm =0;2777 long krm =0; 2778 2778 for (i=0;i<nbt;i++){ 2779 2779 if (triangles[i].link){ // remove triangles … … 2788 2788 Vertex *v0= ta.EdgeVertex(0); 2789 2789 Vertex *v1= ta.EdgeVertex(1); 2790 Int4k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);2790 long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); 2791 2791 if (st[k]<0){ 2792 2792 throw ErrorException(__FUNCT__,exprintf("st[k]<0")); … … 2799 2799 } 2800 2800 } 2801 Int4NbTfillHoll =0;2801 long NbTfillHoll =0; 2802 2802 for (i=0;i<nbt;i++){ 2803 2803 if (triangles[i].link) { … … 2842 2842 } 2843 2843 // OutSidesTriangles = triangles; 2844 // Int4NbOutSidesTriangles = nbt;2844 // long NbOutSidesTriangles = nbt; 2845 2845 2846 2846 // restore triangles; … … 2873 2873 /*}}}1*/ 2874 2874 /*FUNCTION Triangles::GeomToTriangles0{{{1*/ 2875 void Triangles::GeomToTriangles0( Int4inbvx,BamgOpts* bamgopts) {2875 void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts) { 2876 2876 /*Generate mesh from geometry*/ 2877 2877 … … 2880 2880 int i,j,k; 2881 2881 int NbOfCurves=0,NbNewPoints,NbEdgeCurve; 2882 Real8lcurve,lstep,s;2882 double lcurve,lstep,s; 2883 2883 const int MaxSubEdge = 10; 2884 2884 … … 2951 2951 2952 2952 //initialize number of edges and number of edges max 2953 Int4nbex=0;2953 long nbex=0; 2954 2954 nbe=0; 2955 Int4NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;2955 long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 2956 2956 Gh.UnMarkEdges(); 2957 2957 NbOfCurves = 0; … … 2970 2970 if (!ei.Mark() && ei[j].Required()) { 2971 2971 // warning ei.Mark() can be change in loop for(j=0;j<2;j++) 2972 Int4nbvend=0;2972 long nbvend=0; 2973 2973 Edge* PreviousNewEdge=NULL; 2974 2974 … … 3031 3031 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 3032 3032 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 3033 Real8ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric3033 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 3034 3034 3035 3035 /* We are now creating the edges of the mesh from the … … 3043 3043 3044 3044 //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10) 3045 Real8lSubEdge[MaxSubEdge];3045 double lSubEdge[MaxSubEdge]; 3046 3046 3047 3047 //Build Subedges according to the edge length … … 3057 3057 Metric MAs=MA,MBs; 3058 3058 ledge=0; 3059 Real8x =0, xstep= 1./NbSubEdge;3059 double x =0, xstep= 1./NbSubEdge; 3060 3060 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 3061 3061 x += xstep; … … 3067 3067 } 3068 3068 3069 Real8lcurveb = lcurve+ ledge ;3069 double lcurveb = lcurve+ ledge ; 3070 3070 3071 3071 //Now, create corresponding points … … 3073 3073 while (lcurve<=s && s <= lcurveb && nbv < nbvend){ 3074 3074 3075 Real8ss = s-lcurve;3075 double ss = s-lcurve; 3076 3076 3077 3077 /*Find the SubEdge containing ss using Dichotomy*/ 3078 3078 3079 3079 int kk0=-1,kk1=NbSubEdge-1,kkk; 3080 Real8ll0=0,ll1=ledge,llk;3080 double ll0=0,ll1=ledge,llk; 3081 3081 while (kk1-kk0>1){ 3082 3082 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) … … 3089 3089 } 3090 3090 3091 Real8sbb = (ss-ll0 )/(ll1-ll0);3092 Real8bb = (kk1+sbb)/NbSubEdge, aa=1-bb;3091 double sbb = (ss-ll0 )/(ll1-ll0); 3092 double bb = (kk1+sbb)/NbSubEdge, aa=1-bb; 3093 3093 3094 3094 // new vertex on edge … … 3097 3097 vb->ReferenceNumber = e->ref; 3098 3098 vb->DirOfSearch =NoDirOfSearch; 3099 Real8abcisse = k ? bb : aa;3099 double abcisse = k ? bb : aa; 3100 3100 vb->r = e->F( abcisse ); 3101 3101 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); … … 3125 3125 }// for(;;) 3126 3126 vb = b->to; 3127 NbEdgeCurve = Max(( Int4) (lcurve +0.5), (Int4) 1);3127 NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1); 3128 3128 NbNewPoints = NbEdgeCurve-1; 3129 3129 if(!kstep){ … … 3191 3191 /*}}}1*/ 3192 3192 /*FUNCTION Triangles::GeomToTriangles1{{{1*/ 3193 void Triangles::GeomToTriangles1( Int4inbvx,BamgOpts* bamgopts,int KeepVertices){3193 void Triangles::GeomToTriangles1(long inbvx,BamgOpts* bamgopts,int KeepVertices){ 3194 3194 3195 3195 /*Get options*/ … … 3231 3231 // and if you do the loop in background we have the pointeur on geometry 3232 3232 // so do the walk on background 3233 // Int4NbVerticesOnGeomVertex;3233 // long NbVerticesOnGeomVertex; 3234 3234 // VertexOnGeom * VerticesOnGeomVertex; 3235 // Int4NbVerticesOnGeomEdge;3235 // long NbVerticesOnGeomEdge; 3236 3236 // VertexOnGeom * VerticesOnGeomEdge; 3237 3237 … … 3314 3314 // 1.0) recompute the length 3315 3315 // 1.1) compute the vertex 3316 Int4nbex=0,NbVerticesOnGeomEdgex=0;3316 long nbex=0,NbVerticesOnGeomEdgex=0; 3317 3317 for (int step=0; step <2;step++) 3318 3318 { 3319 Int4NbOfNewPoints=0;3320 Int4NbOfNewEdge=0;3321 Int4iedge;3319 long NbOfNewPoints=0; 3320 long NbOfNewEdge=0; 3321 long iedge; 3322 3322 Gh.UnMarkEdges(); 3323 Real8L=0;3323 double L=0; 3324 3324 for (int icurve=0;icurve<Gh.NbOfCurves;icurve++) 3325 3325 { … … 3332 3332 // new curve 3333 3333 // good the find a starting edge 3334 Real8Lstep=0,Lcurve=0;// step between two points (phase==1)3335 Int4NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)3334 double Lstep=0,Lcurve=0;// step between two points (phase==1) 3335 long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 3336 3336 3337 3337 for(int phase=0;phase<=step;phase++) … … 3355 3355 GeometricalEdge *ongequi = peequi->onGeometry; 3356 3356 3357 Real8sNew=Lstep;// abcisse of the new points (phase==1)3357 double sNew=Lstep;// abcisse of the new points (phase==1) 3358 3358 L=0;// length of the curve 3359 Int4i=0;// index of new points on the curve3359 long i=0;// index of new points on the curve 3360 3360 register GeometricalVertex * GA0 = *(*peequi)[k0equi].onGeometry; 3361 3361 Vertex *A0; … … 3386 3386 Vertex & v0=ee[0], & v1=ee[1]; 3387 3387 R2 AB= (R2) v1 - (R2) v0; 3388 Real8L0=L,LAB;3388 double L0=L,LAB; 3389 3389 LAB = LengthInterpole(v0.m,v1.m,AB); 3390 3390 L+= LAB; … … 3410 3410 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; 3411 3411 Edge *e = edges + nbe++; 3412 Real8se= (sNew-L0)/LAB;3412 double se= (sNew-L0)/LAB; 3413 3413 if (se<0 || se>=1.000000001){ 3414 3414 throw ErrorException(__FUNCT__,exprintf("se<0 || se>=1.000000001")); … … 3484 3484 3485 3485 if (!phase) { // 3486 Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg3486 long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg 3487 3487 Lstep = L/NbSegOnCurve; 3488 3488 Lcurve = L; … … 3573 3573 * [0 nbv[ all distincts*/ 3574 3574 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ; 3575 const Int4PrimeNumber= AGoodNumberPrimeWith(nbv) ;3575 const long PrimeNumber= AGoodNumberPrimeWith(nbv) ; 3576 3576 int k0=rand()%nbv; 3577 3577 for (int is3=0; is3<nbv; is3++){ … … 3630 3630 /*Now, add the vertices One by One*/ 3631 3631 3632 Int4NbSwap=0;3632 long NbSwap=0; 3633 3633 if (verbosity>3) printf(" Begining of insertion process...\n"); 3634 3634 3635 for ( Int4icount=2; icount<nbv; icount++) {3635 for (long icount=2; icount<nbv; icount++) { 3636 3636 3637 3637 //Get new vertex … … 3665 3665 3666 3666 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){ 3667 Int4NbSwap = 0;3667 long NbSwap = 0; 3668 3668 for (int is1=0; is1<nbv; is1++) 3669 3669 NbSwap += ordre[is1]->Optim(0,0); … … 3680 3680 /*}}}1*/ 3681 3681 /*FUNCTION Triangles::InsertNewPoints{{{1*/ 3682 Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4& NbTSwap) {3682 long Triangles::InsertNewPoints(long nbvold,long & NbTSwap) { 3683 3683 3684 3684 long int verbosity=0; 3685 Real8seuil= 1.414/2 ;// for two close point3686 Int4i;3687 Int4NbSwap=0;3685 double seuil= 1.414/2 ;// for two close point 3686 long i; 3687 long NbSwap=0; 3688 3688 Icoor2 dete[3]; 3689 3689 3690 3690 //number of new points 3691 const Int4nbvnew=nbv-nbvold;3691 const long nbvnew=nbv-nbvold; 3692 3692 3693 3693 //display info if required … … 3698 3698 3699 3699 /*construction of a random order*/ 3700 const Int4PrimeNumber= AGoodNumberPrimeWith(nbv) ;3700 const long PrimeNumber= AGoodNumberPrimeWith(nbv) ; 3701 3701 //remainder of the division of rand() by nbvnew 3702 Int4k3 = rand()%nbvnew;3702 long k3 = rand()%nbvnew; 3703 3703 //loop over the new points 3704 for ( Int4is3=0; is3<nbvnew; is3++){3705 register Int4j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);3706 register Int4i=nbvold+is3;3704 for (long is3=0; is3<nbvnew; is3++){ 3705 register long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew); 3706 register long i=nbvold+is3; 3707 3707 ordre[i]= vertices + j; 3708 3708 ordre[i]->ReferenceNumber=i; … … 3710 3710 3711 3711 // for all the new point 3712 Int4iv=nbvold;3712 long iv=nbvold; 3713 3713 for (i=nbvold;i<nbv;i++){ 3714 3714 Vertex &vi=*ordre[i]; … … 3721 3721 // a good new point 3722 3722 Vertex &vj = vertices[iv]; 3723 Int4j=vj.ReferenceNumber;3723 long j=vj.ReferenceNumber; 3724 3724 if (&vj!=ordre[j]){ 3725 3725 throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]")); … … 3762 3762 Edge **e= new (Edge* [Gh.nbe]); 3763 3763 3764 Int4i;3764 long i; 3765 3765 for ( i=0;i<Gh.nbe ; i++) 3766 3766 e[i]=NULL; … … 3807 3807 3808 3808 3809 Int4nbqq = (nbt*3)/2;3810 DoubleAndInt4 *qq = new DoubleAndInt4[nbqq];3811 3812 Int4i,ij;3809 long nbqq = (nbt*3)/2; 3810 DoubleAndInt4 *qq = new DoubleAndInt4[nbqq]; 3811 3812 long i,ij; 3813 3813 int j; 3814 Int4k=0;3814 long k=0; 3815 3815 for (i=0;i<nbt;i++) 3816 3816 for (j=0;j<3;j++) … … 3820 3820 HeapSort(qq,k); 3821 3821 3822 Int4kk=0;3822 long kk=0; 3823 3823 for (ij=0;ij<k;ij++) { 3824 3824 i=qq[ij].i3j/3; … … 3845 3845 /*}}}1*/ 3846 3846 /*FUNCTION Triangles::MaxSubDivision{{{1*/ 3847 void Triangles::MaxSubDivision( Real8maxsubdiv) {3847 void Triangles::MaxSubDivision(double maxsubdiv) { 3848 3848 long int verbosity=0; 3849 3849 3850 const Real8maxsubdiv2 = maxsubdiv*maxsubdiv;3850 const double maxsubdiv2 = maxsubdiv*maxsubdiv; 3851 3851 if(verbosity>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv); 3852 3852 // for all the edges 3853 3853 // if the len of the edge is to long 3854 Int4it,nbchange=0;3855 Real8lmax=0;3854 long it,nbchange=0; 3855 double lmax=0; 3856 3856 for (it=0;it<nbt;it++) 3857 3857 { … … 3866 3866 R2 AB= (R2) v1-(R2) v0; 3867 3867 Metric M = v0; 3868 Real8l = M(AB,AB);3868 double l = M(AB,AB); 3869 3869 lmax = Max(lmax,l); 3870 3870 if(l> maxsubdiv2) 3871 3871 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3872 Real8lc = M(AC,AC);3872 double lc = M(AC,AC); 3873 3873 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 3874 3874 D2xD2 Rt1(Rt.inv()); … … 3883 3883 if(l> maxsubdiv2) 3884 3884 { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3885 Real8lc = M(AC,AC);3885 double lc = M(AC,AC); 3886 3886 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; 3887 3887 D2xD2 Rt1(Rt.inv()); … … 3911 3911 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));} 3912 3912 else { // inside 3913 Real8aa[3];3914 Real8s = deta[0]+deta[1]+deta[2];3913 double aa[3]; 3914 double s = deta[0]+deta[1]+deta[2]; 3915 3915 aa[0]=deta[0]/s; 3916 3916 aa[1]=deta[1]/s; … … 3929 3929 3930 3930 int i,j,k; 3931 Int4NbTSwap=0;3932 Int4nbtold=nbt;3933 Int4nbvold=nbv;3934 Int4Headt=0;3935 Int4next_t;3936 Int4* first_np_or_next_t=new Int4[nbtx];3931 long NbTSwap=0; 3932 long nbtold=nbt; 3933 long nbvold=nbv; 3934 long Headt=0; 3935 long next_t; 3936 long* first_np_or_next_t=new long[nbtx]; 3937 3937 Triangle* t=NULL; 3938 3938 … … 4018 4018 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]); 4019 4019 Triangle* tbegin= (Triangle*) ta; 4020 Int4kt;4020 long kt; 4021 4021 do { 4022 4022 kt = Number((Triangle*) ta); … … 4036 4036 delete [] first_np_or_next_t; 4037 4037 4038 Int4NbSwapf =0,NbSwp;4038 long NbSwapf =0,NbSwp; 4039 4039 4040 4040 NbSwp = NbSwapf; … … 4045 4045 /*}}}1*/ 4046 4046 /*FUNCTION Triangles::PreInit{{{1*/ 4047 void Triangles::PreInit( Int4inbvx) {4047 void Triangles::PreInit(long inbvx) { 4048 4048 long int verbosity=0; 4049 4049 … … 4103 4103 /*FUNCTION Triangles::ProjectOnCurve{{{1*/ 4104 4104 GeometricalEdge* Triangles::ProjectOnCurve( Edge & BhAB, Vertex & vA, Vertex & vB, 4105 Real8theta,Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR) {4105 double theta,Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR) { 4106 4106 void *pA=0,*pB=0; 4107 Real8tA=0,tB=0;4107 double tA=0,tB=0; 4108 4108 R2 A=vA,B=vB; 4109 4109 Vertex * pvA=&vA, * pvB=&vB; … … 4166 4166 R2 AB=B-A; 4167 4167 4168 Real8cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);4168 double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); 4169 4169 int kkk=0; 4170 4170 int sens = (cosE01AB>0) ? 1 : 0; 4171 4171 4172 // Real8l=0; // length of the edge AB4173 Real8abscisse = -1;4172 // double l=0; // length of the edge AB 4173 double abscisse = -1; 4174 4174 4175 4175 for (int cas=0;cas<2;cas++){ … … 4180 4180 Vertex *v0=pvA,*v1; 4181 4181 Edge *neee,*eee; 4182 Real8lg =0; // length of the curve4183 Real8te0;4182 double lg =0; // length of the curve 4183 double te0; 4184 4184 // we suppose take the curve's abcisse 4185 4185 for ( eee=e,iii=sens,te0=tA; … … 4194 4194 throw ErrorException(__FUNCT__,exprintf("!eee")); 4195 4195 } 4196 Real8lg0 = lg;4197 Real8dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);4196 double lg0 = lg; 4197 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); 4198 4198 lg += dp; 4199 4199 if (cas && abscisse <= lg) { // ok we find the geom edge 4200 Real8sss = (abscisse-lg0)/dp;4201 Real8thetab = te0*(1-sss)+ sss*iii;4200 double sss = (abscisse-lg0)/dp; 4201 double thetab = te0*(1-sss)+ sss*iii; 4202 4202 if (thetab<0 || thetab>1){ 4203 4203 throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); … … 4212 4212 tB = iii; 4213 4213 4214 Real8lg0 = lg;4214 double lg0 = lg; 4215 4215 if (!eee){ 4216 4216 throw ErrorException(__FUNCT__,exprintf("!eee")); 4217 4217 } 4218 4218 v1 = pvB; 4219 Real8dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);4219 double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); 4220 4220 lg += dp; 4221 4221 abscisse = lg*theta; 4222 4222 if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end 4223 4223 { // ok we find the geom edge 4224 Real8sss = (abscisse-lg0)/dp;4225 Real8thetab = te0*(1-sss)+ sss*tB;4224 double sss = (abscisse-lg0)/dp; 4225 double thetab = te0*(1-sss)+ sss*tB; 4226 4226 if (thetab<0 || thetab>1){ 4227 4227 throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); … … 4321 4321 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) { 4322 4322 long int verbosity=0; 4323 Int4 *renu= new Int4[nbt];4323 long *renu= new long[nbt]; 4324 4324 register Triangle *t0,*t,*te=triangles+nbt; 4325 register Int4k=0,it,i,j;4325 register long k=0,it,i,j; 4326 4326 4327 4327 for ( it=0;it<nbt;it++) … … 4334 4334 } 4335 4335 do { 4336 Int4kt = Number(t);4336 long kt = Number(t); 4337 4337 if (kt<0 || kt >= nbt ){ 4338 4338 throw ErrorException(__FUNCT__,exprintf("kt<0 || kt >= nbt")); … … 4387 4387 /*}}}1*/ 4388 4388 /*FUNCTION Triangles::ReNumberingVertex{{{1*/ 4389 void Triangles::ReNumberingVertex( Int4* renu) {4389 void Triangles::ReNumberingVertex(long * renu) { 4390 4390 // warning be carfull because pointeur 4391 4391 // from on mesh to over mesh 4392 4392 // -- so do ReNumbering a the beginning 4393 4393 Vertex * ve = vertices+nbv; 4394 Int4it,ie,i;4394 long it,ie,i; 4395 4395 4396 4396 for ( it=0;it<nbt;it++) … … 4430 4430 // move the Vertices without a copy of the array 4431 4431 // be carefull not trivial code 4432 Int4j;4432 long j; 4433 4433 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu 4434 4434 if (renu[it] >= 0) // a new sub cycle … … 4461 4461 pmin = vertices[0].r; 4462 4462 pmax = vertices[0].r; 4463 Int4i;4463 long i; 4464 4464 for (i=0;i<nbv;i++) { 4465 4465 pmin.x = Min(pmin.x,vertices[i].r.x); … … 4511 4511 /*FUNCTION Triangles::ShowRegulaty{{{1*/ 4512 4512 void Triangles::ShowRegulaty() const { 4513 const Real8sqrt32=sqrt(3.)*0.5;4514 const Real8aireKh=sqrt32*0.5;4513 const double sqrt32=sqrt(3.)*0.5; 4514 const double aireKh=sqrt32*0.5; 4515 4515 D2 Beq(1,0),Heq(0.5,sqrt32); 4516 4516 D2xD2 Br(D2xD2(Beq,Heq).t()); … … 4522 4522 double alpha2=0; 4523 4523 double area=0,Marea=0; 4524 // Real8 cf= Real8(coefIcoor);4525 // Real8cf2= 6.*cf*cf;4524 // double cf= double(coefIcoor); 4525 // double cf2= 6.*cf*cf; 4526 4526 int nt=0; 4527 4527 for (int it=0;it<nbt;it++) … … 4530 4530 nt++; 4531 4531 Triangle &K=triangles[it]; 4532 Real8area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;4532 double area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.; 4533 4533 area+= area3; 4534 4534 D2xD2 B_Kt(K[0],K[1],K[2]); … … 4540 4540 MatVVP2x2 VMK(MK); 4541 4541 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1)); 4542 Real8betaK=0;4542 double betaK=0; 4543 4543 4544 4544 for (int j=0;j<3;j++) 4545 4545 { 4546 Real8he= Norme2(R2(K[j],K[(j+1)%3]));4546 double he= Norme2(R2(K[j],K[(j+1)%3])); 4547 4547 hmin=Min(hmin,he); 4548 4548 hmax=Max(hmax,he); … … 4576 4576 void Triangles::ShowHistogram() const { 4577 4577 4578 const Int4kmax=10;4579 const Real8llmin = 0.5,llmax=2;4580 const Real8lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);4581 Int4histo[kmax+1];4582 Int4i,it,k, nbedges =0;4578 const long kmax=10; 4579 const double llmin = 0.5,llmax=2; 4580 const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin); 4581 long histo[kmax+1]; 4582 long i,it,k, nbedges =0; 4583 4583 for (i=0;i<=kmax;i++) histo[i]=0; 4584 4584 for (it=0;it<nbt;it++) … … 4595 4595 if ( !& vP || !&vQ) continue; 4596 4596 R2 PQ = vQ.r - vP.r; 4597 Real8l = log(LengthInterpole(vP,vQ,PQ));4597 double l = log(LengthInterpole(vP,vQ,PQ)); 4598 4598 nbedges++; 4599 4599 k = (int) ((l - lmin)*delta); … … 4618 4618 /*}}}1*/ 4619 4619 /*FUNCTION Triangles::SmoothingVertex{{{1*/ 4620 void Triangles::SmoothingVertex(int nbiter, Real8omega ) {4620 void Triangles::SmoothingVertex(int nbiter,double omega ) { 4621 4621 long int verbosity=0; 4622 4622 // if quatree exist remove it end reconstruct … … 4626 4626 Triangle vide; // a triangle to mark the boundary vertex 4627 4627 Triangle ** tstart= new Triangle* [nbv]; 4628 Int4i,j,k;4628 long i,j,k; 4629 4629 // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide 4630 4630 if ( this == & BTh) … … 4641 4641 for (k=0;k<nbiter;k++) 4642 4642 { 4643 Int4i,NbSwap =0;4644 Real8delta =0;4643 long i,NbSwap =0; 4644 double delta =0; 4645 4645 for ( i=0;i<nbv;i++) 4646 4646 if (tstart[i] != &vide) // not a boundary vertex … … 4658 4658 /*}}}1*/ 4659 4659 /*FUNCTION Triangles::SmoothMetric{{{1*/ 4660 void Triangles::SmoothMetric( Real8raisonmax) {4660 void Triangles::SmoothMetric(double raisonmax) { 4661 4661 long int verbosity=0; 4662 4662 … … 4664 4664 if(verbosity > 1) printf(" Triangles::SmoothMetric raisonmax = %g\n",raisonmax); 4665 4665 ReMakeTriangleContainingTheVertex(); 4666 Int4i,j,kch,kk,ip;4667 Int4 *first_np_or_next_t0 = new Int4[nbv];4668 Int4 *first_np_or_next_t1 = new Int4[nbv];4669 Int4Head0 =0,Head1=-1;4670 Real8logseuil= log(raisonmax);4666 long i,j,kch,kk,ip; 4667 long *first_np_or_next_t0 = new long[nbv]; 4668 long *first_np_or_next_t1 = new long[nbv]; 4669 long Head0 =0,Head1=-1; 4670 double logseuil= log(raisonmax); 4671 4671 4672 4672 for(i=0;i<nbv-1;i++) … … 4699 4699 } 4700 4700 R2 Aij = (R2) vj - (R2) vi; 4701 Real8ll = Norme2(Aij);4701 double ll = Norme2(Aij); 4702 4702 if (0) { 4703 Real8hi = ll/vi.m(Aij);4704 Real8hj = ll/vj.m(Aij);4703 double hi = ll/vi.m(Aij); 4704 double hj = ll/vj.m(Aij); 4705 4705 if(hi < hj) 4706 4706 { 4707 Real8dh=(hj-hi)/ll;4707 double dh=(hj-hi)/ll; 4708 4708 if (dh>logseuil) { 4709 4709 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); … … 4715 4715 else 4716 4716 { 4717 Real8li = vi.m(Aij);4717 double li = vi.m(Aij); 4718 4718 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) 4719 4719 if(first_np_or_next_t1[j]<0) // if the metrix change … … 4742 4742 ReNumberingTheTriangleBySubDomain(); 4743 4743 //int nswap =0; 4744 const Int4nfortria( choice ? 4 : 6);4744 const long nfortria( choice ? 4 : 6); 4745 4745 if(withBackground) 4746 4746 { … … 4751 4751 BTh.SetVertexFieldOn(); 4752 4752 4753 Int4newnbt=0,newnbv=0;4754 Int4* kedge = 0;4755 Int4newNbOfQuad=NbOfQuad;4756 Int4* ksplit= 0, * ksplitarray=0;4757 Int4kkk=0;4753 long newnbt=0,newnbv=0; 4754 long * kedge = 0; 4755 long newNbOfQuad=NbOfQuad; 4756 long * ksplit= 0, * ksplitarray=0; 4757 long kkk=0; 4758 4758 int ret =0; 4759 4759 if (nbvx<nbv+nbe) return 1;// 4760 4760 // 1) create the new points by spliting the internal edges 4761 4761 // set the 4762 Int4nbvold = nbv;4763 Int4nbtold = nbt;4764 Int4NbOutTold = NbOutT;4765 Int4NbEdgeOnGeom=0;4766 Int4i;4762 long nbvold = nbv; 4763 long nbtold = nbt; 4764 long NbOutTold = NbOutT; 4765 long NbEdgeOnGeom=0; 4766 long i; 4767 4767 4768 4768 nbt = nbt - NbOutT; // remove all the the ouside triangles 4769 Int4nbtsave = nbt;4769 long nbtsave = nbt; 4770 4770 Triangle * lastT = triangles + nbt; 4771 4771 for (i=0;i<nbe;i++) 4772 4772 if(edges[i].onGeometry) NbEdgeOnGeom++; 4773 Int4newnbe=nbe+nbe;4774 // Int4newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;4775 Int4newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;4776 // Int4newNbVertexOnBThVertex=NbVertexOnBThVertex;4777 Int4newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;4773 long newnbe=nbe+nbe; 4774 // long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex; 4775 long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom; 4776 // long newNbVertexOnBThVertex=NbVertexOnBThVertex; 4777 long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0; 4778 4778 4779 4779 // do allocation for pointeur to the geometry and background … … 4787 4787 // memcpy(newedges,edges,sizeof(Edge)*nbe); 4788 4788 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv); 4789 Int4k=nbv;4790 Int4kk=0;4791 Int4kvb = NbVertexOnBThEdge;4792 Int4kvg = NbVerticesOnGeomEdge;4793 Int4ie =0;4789 long k=nbv; 4790 long kk=0; 4791 long kvb = NbVertexOnBThEdge; 4792 long kvg = NbVerticesOnGeomEdge; 4793 long ie =0; 4794 4794 Edge ** edgesGtoB=0; 4795 4795 if (withBackground) 4796 4796 edgesGtoB= BTh.MakeGeometricalEdgeToEdge(); 4797 Int4ferr=0;4797 long ferr=0; 4798 4798 for (i=0;i<nbe;i++) 4799 4799 newedges[ie].onGeometry=0; … … 4828 4828 ; 4829 4829 // get the Info on background mesh 4830 Real8s = newVertexOnBThEdge[kvb];4830 double s = newVertexOnBThEdge[kvb]; 4831 4831 Vertex & bv0 = newVertexOnBThEdge[kvb][0]; 4832 4832 Vertex & bv1 = newVertexOnBThEdge[kvb][1]; … … 4878 4878 4879 4879 4880 kedge = new Int4[3*nbt+1];4881 ksplitarray = new Int4[nbt+1];4880 kedge = new long[3*nbt+1]; 4881 ksplitarray = new long[nbt+1]; 4882 4882 ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0] 4883 4883 … … 4900 4900 const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 4901 4901 const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]]; 4902 Int4ke =edge4->SortAndFind(Number(v0),Number(v1));4902 long ke =edge4->SortAndFind(Number(v0),Number(v1)); 4903 4903 if (ke>0) 4904 4904 { 4905 Int4ii = Number(tt);4905 long ii = Number(tt); 4906 4906 int jj = ta; 4907 Int4ks = ke + nbvold;4907 long ks = ke + nbvold; 4908 4908 kedge[3*i+j] = ks; 4909 4909 if (ii<nbt) // good triangle 4910 4910 kedge[3*ii+jj] = ks; 4911 4911 Vertex &A=vertices[ks]; 4912 Real8aa,bb,cc,dd;4912 double aa,bb,cc,dd; 4913 4913 if ((dd=Area2(v0.r,v1.r,A.r)) >=0){ 4914 4914 // warning PB roundoff error … … 4951 4951 if ( kedge[3*i+j] < 0) 4952 4952 { 4953 Int4ke =edge4->SortAndFind(Number(v0),Number(v1));4953 long ke =edge4->SortAndFind(Number(v0),Number(v1)); 4954 4954 if (ke<0) // new 4955 4955 { 4956 4956 if (&tt) // internal triangles all the boundary 4957 4957 { // new internal edges 4958 Int4ii = Number(tt);4958 long ii = Number(tt); 4959 4959 int jj = ta; 4960 4960 … … 5010 5010 for (i=0;i<nbtsave;i++){ 5011 5011 int nbmkadj=0; 5012 Int4mkadj [100];5012 long mkadj [100]; 5013 5013 mkadj[0]=i; 5014 Int4kk=ksplit[i]/10;5014 long kk=ksplit[i]/10; 5015 5015 int ke=(int) (ksplit[i]%10); 5016 5016 if (kk>=7 || kk<=0){ … … 5140 5140 vertices[nbv].DirOfSearch =NoDirOfSearch; 5141 5141 //vertices[nbv].i = toI2(vertices[nbv].r); 5142 Real8a3[]={1./3.,1./3.,1./3.};5142 double a3[]={1./3.,1./3.,1./3.}; 5143 5143 vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m); 5144 5144 Vertex * vc = vertices +nbv++; … … 5158 5158 // save all the new triangles 5159 5159 mkadj[nbmkadj++]=i; 5160 Int4jj;5160 long jj; 5161 5161 if (t0.link) 5162 5162 for (jj=nbt;jj<kkk;jj++) … … 5201 5201 for (i=0;i<NbSubDomains;i++) 5202 5202 { 5203 Int4k = subdomains[i].edge- edges;5203 long k = subdomains[i].edge- edges; 5204 5204 subdomains[i].edge = edges+2*k; // spilt all edge in 2 5205 5205 } … … 5242 5242 /*}}}1*/ 5243 5243 /*FUNCTION Triangles::SplitInternalEdgeWithBorderVertices{{{1*/ 5244 Int4Triangles::SplitInternalEdgeWithBorderVertices(){5245 Int4NbSplitEdge=0;5244 long Triangles::SplitInternalEdgeWithBorderVertices(){ 5245 long NbSplitEdge=0; 5246 5246 SetVertexFieldOn(); 5247 Int4it;5248 Int4nbvold=nbv;5247 long it; 5248 long nbvold=nbv; 5249 5249 long int verbosity=2; 5250 5250 for (it=0;it<nbt;it++){ … … 5273 5273 ReMakeTriangleContainingTheVertex(); 5274 5274 if (nbvold!=nbv){ 5275 Int4iv = nbvold;5276 Int4NbSwap = 0;5275 long iv = nbvold; 5276 long NbSwap = 0; 5277 5277 Icoor2 dete[3]; 5278 for ( Int4i=nbvold;i<nbv;i++) {// for all the new point5278 for (long i=nbvold;i<nbv;i++) {// for all the new point 5279 5279 Vertex & vi = vertices[i]; 5280 5280 vi.i = toI2(vi.r); … … 5311 5311 /*}}}1*/ 5312 5312 /*FUNCTION Triangles::TriangleReferenceList{{{1*/ 5313 Int4 Triangles::TriangleReferenceList(Int4* reft) const {5313 long Triangles::TriangleReferenceList(long* reft) const { 5314 5314 long int verbosity=0; 5315 5315 register Triangle *t0,*t; 5316 register Int4k=0, num;5316 register long k=0, num; 5317 5317 5318 5318 //initialize all triangles as -1 (outside) 5319 for ( Int4it=0;it<nbt;it++) reft[it]=-1;5319 for (long it=0;it<nbt;it++) reft[it]=-1; 5320 5320 5321 5321 //loop over all subdomains 5322 for ( Int4i=0;i<NbSubDomains;i++){5322 for (long i=0;i<NbSubDomains;i++){ 5323 5323 5324 5324 //first triangle of the subdomain i … … 5591 5591 /*}}}1*/ 5592 5592 /*FUNCTION AGoodNumberPrimeWith{{{1*/ 5593 Int4 AGoodNumberPrimeWith(Int4n){5593 long AGoodNumberPrimeWith(long n){ 5594 5594 5595 5595 //list of big prime numbers 5596 const Int4BigPrimeNumber[] ={ 567890359L,5596 const long BigPrimeNumber[] ={ 567890359L, 5597 5597 567890431L, 567890437L, 567890461L, 567890471L, 5598 5598 567890483L, 567890489L, 567890497L, 567890507L, … … 5600 5600 5601 5601 //initialize o and pi 5602 Int4o =0;5603 Int4pi=BigPrimeNumber[1];5602 long o =0; 5603 long pi=BigPrimeNumber[1]; 5604 5604 5605 5605 //loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber) … … 5607 5607 5608 5608 //compute r, rest of the remainder of the division of BigPrimeNumber[i] by n 5609 Int4r = BigPrimeNumber[i] % n;5609 long r = BigPrimeNumber[i] % n; 5610 5610 5611 5611 /*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/ 5612 Int4oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));5612 long oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); 5613 5613 if ( o < oo){ 5614 5614 o=oo; -
issm/trunk/src/c/Bamgx/objects/Triangles.h
r3232 r3243 38 38 Geometry & Gh; // Geometry 39 39 Triangles & BTh; // Background Mesh Bth==*this =>no background 40 Int4NbRef; // counter of ref on the this class if 0 we can delete41 Int4nbvx,nbtx; // nombre max de sommets , de triangles42 Int4nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference,43 Int4NbOfQuad; // nb of quadrangle44 Int4NbSubDomains;45 Int4NbOutT; // Nb of oudeside triangle46 Int4NbOfTriangleSearchFind;47 Int4NbOfSwapTriangle;40 long NbRef; // counter of ref on the this class if 0 we can delete 41 long nbvx,nbtx; // nombre max de sommets , de triangles 42 long nt,nbv,nbt,nbiv,nbe; // nb of legal triangles, nb of vertex, of triangles, of initial vertices, of edges with reference, 43 long NbOfQuad; // nb of quadrangle 44 long NbSubDomains; 45 long NbOutT; // Nb of oudeside triangle 46 long NbOfTriangleSearchFind; 47 long NbOfSwapTriangle; 48 48 Vertex* vertices; 49 Int4NbVerticesOnGeomVertex;49 long NbVerticesOnGeomVertex; 50 50 VertexOnGeom * VerticesOnGeomVertex; 51 Int4NbVerticesOnGeomEdge;51 long NbVerticesOnGeomEdge; 52 52 VertexOnGeom * VerticesOnGeomEdge; 53 Int4NbVertexOnBThVertex;53 long NbVertexOnBThVertex; 54 54 VertexOnVertex *VertexOnBThVertex; 55 Int4NbVertexOnBThEdge;55 long NbVertexOnBThEdge; 56 56 VertexOnEdge *VertexOnBThEdge; 57 Int4NbCrackedVertices;58 Int4NbCrackedEdges;57 long NbCrackedVertices; 58 long NbCrackedEdges; 59 59 CrackedEdge *CrackedEdges; 60 60 R2 pmin,pmax; // extrema 61 Real8coefIcoor; // coef to integer Icoor1;61 double coefIcoor; // coef to integer Icoor1; 62 62 Triangle* triangles; 63 63 Edge* edges; … … 70 70 Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); 71 71 Triangles(double* index,double* x,double* y,int nods,int nels); 72 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0, Int4nbvxx=0 ); //copy operator72 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,long nbvxx=0 ); //copy operator 73 73 Triangles(const Triangles &,const int *flag,const int *bb,BamgOpts* bamgopts); // truncature 74 Triangles( Int4nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {74 Triangles(long nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) { 75 75 try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);} 76 76 catch(...) { this->~Triangles(); throw; } 77 77 } 78 Triangles( Int4nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){78 Triangles(long nbvx,Geometry & G,BamgOpts* bamgopts) :Gh(G),BTh(*this){ 79 79 try { GeomToTriangles0(nbvx,bamgopts);} 80 80 catch(...) { this->~Triangles(); throw; } … … 83 83 84 84 //Operators 85 const Vertex & operator[] ( Int4i) const { return vertices[i];};86 Vertex & operator[]( Int4i) { return vertices[i];};87 const Triangle & operator() ( Int4i) const { return triangles[i];};88 Triangle & operator()( Int4i) { return triangles[i];};85 const Vertex & operator[] (long i) const { return vertices[i];}; 86 Vertex & operator[](long i) { return vertices[i];}; 87 const Triangle & operator() (long i) const { return triangles[i];}; 88 Triangle & operator()(long i) { return triangles[i];}; 89 89 90 90 //Methods 91 91 void SetIntCoor(const char * from =0); 92 Real8MinimalHmin() {return 2.0/coefIcoor;}93 Real8MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}92 double MinimalHmin() {return 2.0/coefIcoor;} 93 double MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 94 94 I2 toI2(const R2 & P) const { 95 95 return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) … … 103 103 void Renumber(BamgOpts* bamgopts); 104 104 void FindSubDomain(int OutSide=0); 105 Int4 TriangleReferenceList(Int4*) const;105 long TriangleReferenceList(long *) const; 106 106 void ShowHistogram() const; 107 107 void ShowRegulaty() const; 108 108 void ReMakeTriangleContainingTheVertex(); 109 109 void UnMarkUnSwapTriangle(); 110 void SmoothMetric( Real8raisonmax) ;111 void BoundAnisotropy( Real8anisomax,double hminaniso= 1e-100) ;112 void MaxSubDivision( Real8maxsubdiv);110 void SmoothMetric(double raisonmax) ; 111 void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ; 112 void MaxSubDivision(double maxsubdiv); 113 113 Edge** MakeGeometricalEdgeToEdge(); 114 114 void SetVertexFieldOn(); 115 115 void SetVertexFieldOnBTh(); 116 Int4SplitInternalEdgeWithBorderVertices();116 long SplitInternalEdgeWithBorderVertices(); 117 117 void MakeQuadrangles(double costheta); 118 118 int SplitElement(int choice); 119 119 void MakeQuadTree(); 120 120 void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1); 121 Int4 InsertNewPoints(Int4 nbvold,Int4& NbTSwap) ;121 long InsertNewPoints(long nbvold,long & NbTSwap) ; 122 122 void ReNumberingTheTriangleBySubDomain(bool justcompress=false); 123 void ReNumberingVertex( Int4* renu);124 void SmoothingVertex(int =3, Real8=0.3);123 void ReNumberingVertex(long * renu); 124 void SmoothingVertex(int =3,double=0.3); 125 125 Metric MetricAt (const R2 &) const; 126 GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B, Real8theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR);127 Int4Number(const Triangle & t) const { return &t - triangles;}128 Int4Number(const Triangle * t) const { return t - triangles;}129 Int4Number(const Vertex & t) const { return &t - vertices;}130 Int4Number(const Vertex * t) const { return t - vertices;}131 Int4Number(const Edge & t) const { return &t - edges;}132 Int4Number(const Edge * t) const { return t - edges;}133 Int4Number2(const Triangle * t) const {126 GeometricalEdge* ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,double theta, Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 127 long Number(const Triangle & t) const { return &t - triangles;} 128 long Number(const Triangle * t) const { return t - triangles;} 129 long Number(const Vertex & t) const { return &t - vertices;} 130 long Number(const Vertex * t) const { return t - vertices;} 131 long Number(const Edge & t) const { return &t - edges;} 132 long Number(const Edge * t) const { return t - edges;} 133 long Number2(const Triangle * t) const { 134 134 return t - triangles; 135 135 } … … 139 139 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); 140 140 void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts); 141 void ReadMetric(BamgOpts* bamgopts,const Real8 hmin,const Real8 hmax,const Real8coef);141 void ReadMetric(BamgOpts* bamgopts,const double hmin,const double hmax,const double coef); 142 142 void WriteMetric(BamgOpts* bamgopts); 143 143 void AddMetric(BamgOpts* bamgopts); … … 153 153 154 154 private: 155 void GeomToTriangles1( Int4nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption156 void GeomToTriangles0( Int4nbvx,BamgOpts* bamgopts);// the real constructor mesh generator157 void PreInit( Int4);155 void GeomToTriangles1(long nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption 156 void GeomToTriangles0(long nbvx,BamgOpts* bamgopts);// the real constructor mesh generator 157 void PreInit(long); 158 158 159 159 }; -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r3236 r3243 14 14 /*Methods*/ 15 15 /*FUNCTION Vertex::Smoothing{{{1*/ 16 Real8 Vertex::Smoothing(Triangles &Th,const Triangles &BTh,Triangle* &tstart ,Real8omega){16 double Vertex::Smoothing(Triangles &Th,const Triangles &BTh,Triangle* &tstart ,double omega){ 17 17 18 18 register Vertex* s=this; … … 30 30 31 31 R2 Q = vQ,QP(P-Q); 32 Real8lQP = LengthInterpole(vP,vQ,QP);32 double lQP = LengthInterpole(vP,vQ,QP); 33 33 PNew += Q+QP/Max(lQP,1e-20); 34 34 kk ++; … … 43 43 } while ( tbegin != tria); 44 44 if (kk<4) return 0; 45 PNew = PNew/( Real8)kk;45 PNew = PNew/(double)kk; 46 46 R2 Xmove((PNew-P)*omega); 47 47 PNew = P+Xmove; 48 Real8delta=Norme2_2(Xmove);48 double delta=Norme2_2(Xmove); 49 49 50 50 Icoor2 deta[3]; … … 60 60 } 61 61 else { // inside 62 Real8aa[3];63 Real8s = deta[0]+deta[1]+deta[2];62 double aa[3]; 63 double s = deta[0]+deta[1]+deta[2]; 64 64 aa[0]=deta[0]/s; 65 65 aa[1]=deta[1]/s; … … 90 90 { 91 91 vP = vPsave; 92 Real8qold =QuadQuality(*v0,*v1,*v2,*v3);92 double qold =QuadQuality(*v0,*v1,*v2,*v3); 93 93 vP = vPnew; 94 Real8qnew =QuadQuality(*v0,*v1,*v2,*v3);94 double qnew =QuadQuality(*v0,*v1,*v2,*v3); 95 95 if (qnew<qold) ok = 1; 96 96 } -
issm/trunk/src/c/Bamgx/objects/Vertex.h
r3242 r3243 25 25 R2 r; // real coordinates 26 26 Metric m; 27 Int4ReferenceNumber;27 long ReferenceNumber; 28 28 Direction DirOfSearch; 29 29 short vint; // the vertex number in triangle; varies between 0 and 2 in t 30 30 union { 31 31 Triangle* t; // one triangle which is containing the vertex 32 Int4color;32 long color; 33 33 Vertex* to; // use in geometry Vertex to now the Mesh Vertex associed 34 34 VertexOnGeom* onGeometry; // if vint == 8; // set with Triangles::SetVertexFieldOn() … … 41 41 operator const R2 & () const {return r;} // Cast operator 42 42 operator Metric () const {return m;} // Cast operator 43 Real8operator()(R2 x) const { return m(x);} // Get x in the metric m43 double operator()(R2 x) const { return m(x);} // Get x in the metric m 44 44 45 45 //methods (No constructor and no destructors...) 46 Real8 Smoothing(Triangles & ,const Triangles & ,Triangle * & ,Real8=1);46 double Smoothing(Triangles & ,const Triangles & ,Triangle * & ,double =1); 47 47 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); 48 48 void Echo(); … … 50 50 51 51 //inline functions 52 inline Int4Optim(int =1,int =0);52 inline long Optim(int =1,int =0); 53 53 inline void Set(const Vertex & rec,const Triangles &,Triangles &); 54 54 }; -
issm/trunk/src/c/Bamgx/objects/VertexOnEdge.h
r3232 r3243 21 21 Vertex* v; 22 22 Edge* be; 23 Real8abcisse;23 double abcisse; 24 24 25 25 //Constructors 26 VertexOnEdge( Vertex * w, Edge *bw, Real8s) :v(w),be(bw),abcisse(s) {}26 VertexOnEdge( Vertex * w, Edge *bw,double s) :v(w),be(bw),abcisse(s) {} 27 27 VertexOnEdge(){} 28 28 29 29 //Operators 30 operator Real8() const { return abcisse;}30 operator double () const { return abcisse;} 31 31 operator Vertex* () const { return v;} 32 32 operator Edge* () const { return be;} … … 37 37 38 38 //Inline methods 39 inline void Set(const Triangles &, Int4,Triangles &);39 inline void Set(const Triangles &,long,Triangles &); 40 40 }; 41 41 -
issm/trunk/src/c/Bamgx/objects/VertexOnGeom.h
r3232 r3243 22 22 23 23 Vertex* mv; 24 Real8abscisse;24 double abscisse; 25 25 union{ 26 26 GeometricalVertex * gv; // if abscisse <0; … … 31 31 VertexOnGeom(): mv(0),abscisse(0){gv=0;} 32 32 VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} 33 VertexOnGeom(Vertex & m,GeometricalEdge &g, Real8s) : mv(&m),abscisse(s){ge=&g;}33 VertexOnGeom(Vertex & m,GeometricalEdge &g,double s) : mv(&m),abscisse(s){ge=&g;} 34 34 35 35 //Operators … … 37 37 operator GeometricalVertex * () const {return gv;} 38 38 operator GeometricalEdge * () const {return ge;} 39 operator const Real8& () const {return abscisse;}39 operator const double & () const {return abscisse;} 40 40 41 41 //Methods … … 46 46 47 47 //Inline methods 48 inline void Set(const Triangles &, Int4,Triangles &);48 inline void Set(const Triangles &,long,Triangles &); 49 49 inline void Set(const VertexOnGeom&,const Triangles &,Triangles &); 50 50 -
issm/trunk/src/c/Bamgx/objects/VertexOnVertex.h
r3232 r3243 29 29 30 30 //Inline methods 31 inline void Set(const Triangles &, Int4,Triangles &);31 inline void Set(const Triangles &,long,Triangles &); 32 32 }; 33 33
Note:
See TracChangeset
for help on using the changeset viewer.
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)