source: issm/trunk/src/c/Bamgx/Mesh2.h@ 3232

Last change on this file since 3232 was 3232, checked in by Mathieu Morlighem, 15 years ago

Separated all objects (still some work to do)
File size: 14.8 KB
Line 
1#ifndef _MESH2_H_
2#define _MESH2_H_
3
4#include "../objects/objects.h"
5#include "../shared/shared.h"
6#include "../include/macros.h"
7#include "../toolkits/toolkits.h"
8
9#include "meshtype.h"
10
11#include "objects/Metric.h"
12#include "objects/DoubleAndInt4.h"
13#include "objects/Direction.h"
14#include "objects/Vertex.h"
15#include "objects/TriangleAdjacent.h"
16#include "objects/Edge.h"
17#include "objects/GeometricalVertex.h"
18#include "objects/GeometricalEdge.h"
19#include "objects/Curve.h"
20#include "objects/Triangle.h"
21#include "objects/ListofIntersectionTriangles.h"
22#include "objects/GeometricalSubDomain.h"
23#include "objects/SubDomain.h"
24#include "objects/VertexOnGeom.h"
25#include "objects/VertexOnVertex.h"
26#include "objects/VertexOnEdge.h"
27#include "objects/CrackedEdge.h"
28#include "objects/Triangles.h"
29#include "objects/Geometry.h"
30
31namespace bamg {
32
33 /*INLINE functions{{{1*/
34 // to sort in descending order
35 template<class T> inline void HeapSort(T *c,long n){
36 /*Intermediary*/
37 int l,j,r,i;
38 T crit;
39 c--; //the array must starts at 1 and not 0
40 if(n<=1) return; //return if size <=1
41 l=n/2+1; //initialize l and r
42 r=n;
43 for(;;){
44 if(l<=1){
45 crit =c[r];
46 c[r--]=c[1];
47 if (r==1){c[1]=crit; return;}
48 }
49 else crit = c[--l];
50 j=l;
51 for(;;){
52 i=j;
53 j=2*j;
54 if (j>r) {c[i]=crit;break;}
55 if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)
56 if (crit < c[j]) c[i]=c[j]; //c[j] > crit -> stock this large value in i(<j)
57 else{c[i]=crit;break;} //c[j] < crit -> stock crit in i (<j)
58 }
59 }
60 }
61 // to sort in descending order and return ordering
62 template<class T> inline void HeapSort(int** porder,T* c,int n){
63 /*Intermediary*/
64 int l,j,r,i;
65 T crit;
66 int pos;
67 int* order=NULL;
68 order=(int*)xmalloc(n*sizeof(int));
69 for(i=0;i<n;i++) order[i]=i+1;
70 c--; //the array must starts at 1 and not 0
71 order--;
72 if(n<=1) return; //return if size <=1
73 l=n/2+1; //initialize l and r
74 r=n;
75 for(;;){
76 if(l<=1){
77 crit =c[r]; pos=order[r];
78 c[r--]=c[1]; order[r+1]=order[1];
79 if (r==1){
80 c[1]=crit; order[1]=pos;
81 order++;
82 *porder=order;
83 return;
84 }
85 }
86 else {crit=c[--l]; pos=order[l];}
87 j=l;
88 for(;;){
89 i=j;
90 j=2*j;
91 if (j>r) {c[i]=crit;order[i]=pos;break;}
92 if ((j<r) && (c[j] < c[j+1]))j++;
93 if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}
94 else{c[i]=crit;order[i]=pos;break;}
95 }
96 }
97 }
98 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){
99 return A[0]*( B[1]*C[2]-B[2]*C[1])
100 - A[1]*( B[0]*C[2]-B[2]*C[0])
101 + A[2]*( B[0]*C[1]-B[1]*C[0]);
102 }
103 inline TriangleAdjacent Adj(const TriangleAdjacent & a)
104 { return a.Adj();}
105
106 inline TriangleAdjacent Next(const TriangleAdjacent & ta)
107 { return TriangleAdjacent(ta.t,NextEdge[ta.a]);}
108
109 inline TriangleAdjacent Previous(const TriangleAdjacent & ta)
110 { return TriangleAdjacent(ta.t,PreviousEdge[ta.a]);}
111 inline void Adj(GeometricalEdge * & on,int &i)
112 {int j=i;i=on->DirAdj[i];on=on->Adj[j];}
113 inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc)
114 {
115 Real4 ret;
116 I2 ia=va,ib=vb,ic=vc;
117 I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;
118 Icoor2 deta=Det(ab,ac);
119 if (deta <=0) ret = -1;
120 else {
121 Real8 a = sqrt((Real8) (ac,ac)),
122 b = sqrt((Real8) (bc,bc)),
123 c = sqrt((Real8) (ab,ab)),
124 p = a+b+c;
125 Real8 h= Max(Max(a,b),c),ro=deta/p;
126 ret = ro/h;}
127 return ret;
128 }
129 Icoor2 inline det(const Vertex & a,const Vertex & b,const Vertex & c){
130 register Icoor2 bax = b.i.x - a.i.x ,bay = b.i.y - a.i.y;
131 register Icoor2 cax = c.i.x - a.i.x ,cay = c.i.y - a.i.y;
132 return bax*cay - bay*cax;
133 }
134 inline TriangleAdjacent FindTriangleAdjacent(Edge &E){
135 Vertex * a = E.v[0];
136 Vertex * b = E.v[1];
137
138 Triangle * t = a->t;
139 int i = a->vint;
140 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
141 if (!t || i<0 || i>=3){
142 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
143 }
144 if ( a!=(*t)(i)){
145 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
146 }
147 int k=0;
148 do { // turn around vertex in direct sens (trigo)
149 k++;
150 if (k>=20000){
151 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
152 }
153 // in no crack => ta.EdgeVertex(1) == a otherwise ???
154 if (ta.EdgeVertex(1) == a && ta.EdgeVertex(0) == b) return ta; // find
155 ta = ta.Adj();
156 if (ta.EdgeVertex(0) == a && ta.EdgeVertex(1) == b) return ta; // find
157 --ta;
158 } while (t != (Triangle *)ta);
159 throw ErrorException(__FUNCT__,exprintf("FindTriangleAdjacent: triangle not found"));
160 return TriangleAdjacent(0,0);//for compiler
161 }
162 inline Vertex* TheVertex(Vertex * a){
163 // give a unique vertex with smallest number
164 // in case on crack in mesh
165 Vertex * r(a), *rr;
166 Triangle * t = a->t;
167 int i = a->vint;
168 TriangleAdjacent ta(t,EdgesVertexTriangle[i][0]); // Previous edge
169 if (!t || i<0 || i>=3){
170 throw ErrorException(__FUNCT__,exprintf("!t || i<0 !! i>=3"));
171 }
172 if ( a!=(*t)(i)){
173 throw ErrorException(__FUNCT__,exprintf("a!=(*t)(i)"));
174 }
175 int k=0;
176 do { // turn around vertex in direct sens (trigo)
177 k++;
178 if (k>=20000){
179 throw ErrorException(__FUNCT__,exprintf("k>=20000"));
180 }
181 // in no crack => ta.EdgeVertex(1) == a
182 if ((rr=ta.EdgeVertex(0)) < r) r = rr;
183 ta = ta.Adj();
184 if ((rr=ta.EdgeVertex(1)) < r) r =rr;
185 --ta;
186 } while (t != (Triangle*) ta);
187 return r;
188 }
189
190 /*}}}1*/
191
192 /*INLINE functions of CLASS GeometricalVertex{{{1*/
193 inline void GeometricalVertex::Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){
194 *this = rec;
195 }
196 /*}}}1*/
197 /*INLINE functions of CLASS VertexOnVertex{{{1*/
198 inline void VertexOnVertex::Set(const Triangles & Th ,Int4 i,Triangles & ThNew) {
199 *this = Th.VertexOnBThVertex[i];
200 v = ThNew.vertices + Th.Number(v);
201
202 }
203/*}}}1*/
204 /*INLINE functions of CLASS Triangles{{{1*/
205 inline void Triangles::ReMakeTriangleContainingTheVertex(){
206 register Int4 i;
207 for ( i=0;i<nbv;i++){
208 vertices[i].vint=0;
209 vertices[i].t=0;
210 }
211 for ( i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
212 }
213
214 inline void Triangles::UnMarkUnSwapTriangle()
215 {
216 register Int4 i;
217 for ( i=0;i<nbt;i++)
218 for(int j=0;j<3;j++)
219 triangles[i].SetUnMarkUnSwap(j);
220 }
221
222 inline void Triangles::SetVertexFieldOn(){
223 for (Int4 i=0;i<nbv;i++)
224 vertices[i].onGeometry=0;
225 for (Int4 j=0;j<NbVerticesOnGeomVertex;j++ )
226 VerticesOnGeomVertex[j].SetOn();
227 for (Int4 k=0;k<NbVerticesOnGeomEdge;k++ )
228 VerticesOnGeomEdge[k].SetOn();
229 }
230 inline void Triangles::SetVertexFieldOnBTh(){
231 for (Int4 i=0;i<nbv;i++)
232 vertices[i].onGeometry=0;
233 for (Int4 j=0;j<NbVertexOnBThVertex;j++ )
234 VertexOnBThVertex[j].SetOnBTh();
235 for (Int4 k=0;k<NbVertexOnBThEdge;k++ )
236 VertexOnBThEdge[k].SetOnBTh();
237
238 }
239 /*}}}1*/
240 /*INLINE functions of CLASS Triangle{{{1*/
241 inline Triangle* Triangle::Quadrangle(Vertex * & v0,Vertex * & v1,Vertex * & v2,Vertex * & v3) const
242 {
243 // return the other triangle of the quad if a quad or 0 if not a quat
244 Triangle * t =0;
245 if (link) {
246 int a=-1;
247 if (TriaAdjSharedEdge[0] & 16 ) a=0;
248 if (TriaAdjSharedEdge[1] & 16 ) a=1;
249 if (TriaAdjSharedEdge[2] & 16 ) a=2;
250 if (a>=0) {
251 t = TriaAdjTriangles[a];
252 // if (t-this<0) return 0;
253 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
254 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
255 v1 = TriaVertices[OppositeEdge[a]];
256 v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
257 }
258 }
259 return t;
260 }
261
262 inline double Triangle::QualityQuad(int a,int option) const
263 { // first do the logique part
264 double q;
265 if (!link || TriaAdjSharedEdge[a] &4)
266 q= -1;
267 else {
268 Triangle * t = TriaAdjTriangles[a];
269 if (t-this<0) q= -1;// because we do 2 times
270 else if (!t->link ) q= -1;
271 else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16 || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )
272 q= -1;
273 else if(option)
274 {
275 const Vertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
276 const Vertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
277 const Vertex & v1 = *TriaVertices[OppositeEdge[a]];
278 const Vertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
279 q = QuadQuality(v0,v1,v2,v3); // do the float part
280 }
281 else q= 1;
282 }
283 return q;
284 }
285 inline void Triangle::Set(const Triangle & rec,const Triangles & Th ,Triangles & ThNew)
286 {
287 *this = rec;
288 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices + Th.Number(TriaVertices[0]);
289 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices + Th.Number(TriaVertices[1]);
290 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices + Th.Number(TriaVertices[2]);
291 if(TriaAdjTriangles[0]) TriaAdjTriangles[0] = ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
292 if(TriaAdjTriangles[1]) TriaAdjTriangles[1] = ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
293 if(TriaAdjTriangles[2]) TriaAdjTriangles[2] = ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
294 if (link >= Th.triangles && link < Th.triangles + Th.nbt)
295 link = ThNew.triangles + Th.Number(link);
296 }
297 inline Triangle::Triangle(Triangles *Th,Int4 i,Int4 j,Int4 k) {
298 Vertex *v=Th->vertices;
299 Int4 nbv = Th->nbv;
300 if (i<0 || j<0 || k<0){
301 throw ErrorException(__FUNCT__,exprintf("i<0 || j<0 || k<0"));
302 }
303 if (i>=nbv || j>=nbv || k>=nbv){
304 throw ErrorException(__FUNCT__,exprintf("i>=nbv || j>=nbv || k>=nbv"));
305 }
306 TriaVertices[0]=v+i;
307 TriaVertices[1]=v+j;
308 TriaVertices[2]=v+k;
309 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
310 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
311 det=0;
312 }
313 inline Triangle::Triangle(Vertex *v0,Vertex *v1,Vertex *v2){
314 TriaVertices[0]=v0;
315 TriaVertices[1]=v1;
316 TriaVertices[2]=v2;
317 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
318 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
319 if (v0) det=0;
320 else {
321 det=-1;
322 link=NULL;};
323 }
324 inline Real4 Triangle::qualite()
325 {
326 return det < 0 ? -1 : bamg::qualite(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
327 }
328
329 /*}}}1*/
330 /*INLINE functions of CLASS Edge{{{1*/
331 inline void Edge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
332 {
333 *this = Th.edges[i];
334 v[0] = ThNew.vertices + Th.Number(v[0]);
335 v[1] = ThNew.vertices + Th.Number(v[1]);
336 if (onGeometry)
337 onGeometry = ThNew.Gh.edges+Th.Gh.Number(onGeometry);
338 if (adj[0]) adj[0] = ThNew.edges + Th.Number(adj[0]);
339 if (adj[1]) adj[1] = ThNew.edges + Th.Number(adj[1]);
340
341 }
342
343 /*}}}1*/
344 /*INLINE functions of CLASS GeometricalEdge{{{1*/
345 inline void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew)
346 {
347 *this = rec;
348 v[0] = GhNew.vertices + Gh.Number(v[0]);
349 v[1] = GhNew.vertices + Gh.Number(v[1]);
350 if (Adj[0]) Adj[0] = GhNew.edges + Gh.Number(Adj[0]);
351 if (Adj[1]) Adj[1] = GhNew.edges + Gh.Number(Adj[1]);
352 }
353 /*}}}1*/
354 /*INLINE functions of CLASS Curve{{{1*/
355 inline void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew)
356 {
357 *this = rec;
358 be = GhNew.edges + Gh.Number(be);
359 ee = GhNew.edges + Gh.Number(ee);
360 if(next) next= GhNew.curves + Gh.Number(next);
361 }
362 /*}}}1*/
363 /*INLINE functions of CLASS SubDomain{{{1*/
364 inline void SubDomain::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
365 {
366 *this = Th.subdomains[i];
367 if ( head-Th.triangles<0 || head-Th.triangles>=Th.nbt){
368 throw ErrorException(__FUNCT__,exprintf("head-Th.triangles<0 || head-Th.triangles>=Th.nbt"));
369 }
370 head = ThNew.triangles + Th.Number(head) ;
371 if (edge-Th.edges<0 || edge-Th.edges>=Th.nbe);{
372 throw ErrorException(__FUNCT__,exprintf("edge-Th.edges<0 || edge-Th.edges>=Th.nbe"));
373 }
374 edge = ThNew.edges+ Th.Number(edge);
375 }
376 /*}}}1*/
377 /*INLINE functions of CLASS GeometricalSubDomain{{{1*/
378 inline void GeometricalSubDomain::Set(const GeometricalSubDomain & rec,const Geometry & Gh ,const Geometry & GhNew)
379 {
380 *this = rec;
381 edge = Gh.Number(edge) + GhNew.edges;
382 }
383 /*}}}1*/
384 /*INLINE functions of CLASS Vertex{{{1*/
385 Int4 inline Vertex::Optim(int i,int koption){
386 Int4 ret=0;
387 if ( t && (vint >= 0 ) && (vint <3) ){
388 ret = t->Optim(vint,koption);
389 if(!i){
390 t =0; // for no future optime
391 vint= 0;
392 }
393 }
394 return ret;
395 }
396 inline void Vertex::Set(const Vertex & rec,const Triangles & ,Triangles & ){
397 *this = rec;
398 }
399 /*}}}1*/
400 /*INLINE functions of CLASS VertexOnEdge{{{1*/
401 inline void VertexOnEdge::Set(const Triangles & Th ,Int4 i,Triangles & ThNew)
402 {
403 *this = Th.VertexOnBThEdge[i];
404 v = ThNew.vertices + Th.Number(v);
405 }
406 /*}}}1*/
407 /*INLINE functions of CLASS VertexOnGeom{{{1*/
408 inline void VertexOnGeom::Set(const VertexOnGeom & rec,const Triangles & Th ,Triangles & ThNew)
409 {
410 *this = rec;
411 mv = ThNew.vertices + Th.Number(mv);
412 if (gv)
413 if (abscisse < 0 )
414 gv = ThNew.Gh.vertices + Th.Gh.Number(gv);
415 else
416 ge = ThNew.Gh.edges + Th.Gh.Number(ge);
417
418 }
419 /*}}}1*/
420 /*INLINE functions of CLASS TriangleAdjacent{{{1*/
421 inline void TriangleAdjacent::SetAdj2(const TriangleAdjacent & ta, int l ){
422 //set Adjacent Triangle of a triangle
423 if(t) {
424 t->TriaAdjTriangles[a]=ta.t;
425 t->TriaAdjSharedEdge[a]=ta.a|l;
426 }
427 if(ta.t) {
428 ta.t->TriaAdjTriangles[ta.a] = t ;
429 ta.t->TriaAdjSharedEdge[ta.a] = a| l ;
430 }
431 }
432 inline int TriangleAdjacent::Locked() const
433 { return t->TriaAdjSharedEdge[a] &4;}
434 inline int TriangleAdjacent::Cracked() const
435 { return t->TriaAdjSharedEdge[a] &32;}
436 inline int TriangleAdjacent::GetAllFlag_UnSwap() const
437 { return t->TriaAdjSharedEdge[a] & 1012;} // take all flag except MarkUnSwap
438 inline int TriangleAdjacent::MarkUnSwap() const
439 { return t->TriaAdjSharedEdge[a] &8;}
440 inline void TriangleAdjacent::SetLock(){ t->SetLocked(a);}
441 inline void TriangleAdjacent::SetCracked() { t->SetCracked(a);}
442 inline TriangleAdjacent TriangleAdjacent::Adj() const
443 { return t->Adj(a);}
444 inline Vertex * TriangleAdjacent::EdgeVertex(const int & i) const
445 {return t->TriaVertices[VerticesOfTriangularEdge[a][i]]; }
446 inline Vertex * TriangleAdjacent::OppositeVertex() const
447 {return t->TriaVertices[bamg::OppositeVertex[a]]; }
448 inline Icoor2 & TriangleAdjacent::det() const
449 { return t->det;}
450 int inline TriangleAdjacent::swap()
451 { return t->swap(a);}
452 /*}}}1*/
453
454 /*Other prototypes {{{1*/
455 TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ;
456 TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b);
457 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);
458 void swap(Triangle *t1,Int1 a1,
459 Triangle *t2,Int1 a2,
460 Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2);
461 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
462 TriangleAdjacent & tt1,Icoor2 & dets1,
463 Icoor2 & detsa,Icoor2 & detsb, int & nbswap);
464 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ;
465 /*}}}1*/
466
467}
468#endif
Note: See TracBrowser for help on using the repository browser.