diff --git a/contrib/bamg/bamglib/Mesh2.cpp b/contrib/bamg/bamglib/Mesh2.cpp
index ebfac325dfb11afac7f8bd672393b2403d313d92..101df21b0445fdb463e9de50a8140500491801c6 100644
--- a/contrib/bamg/bamglib/Mesh2.cpp
+++ b/contrib/bamg/bamglib/Mesh2.cpp
@@ -1,26 +1,26 @@
// -*- Mode : c++ -*-
//
-// SUMMARY :
-// USAGE :
-// ORG :
+// SUMMARY :
+// USAGE :
+// ORG :
// AUTHOR : Frederic Hecht
// E-MAIL : hecht@ann.jussieu.fr
//
/*
-
+
This file is part of Freefem++
-
+
Freefem++ is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.
-
+
Freefem++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
-
+
You should have received a copy of the GNU Lesser General Public License
along with Freefem++; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
@@ -36,10 +36,10 @@
extern bool withrgraphique;
#include <stdio.h>
#include <string.h>
-#include <math.h>
+#include <math.h>
#include <time.h>
#include <iostream>
-using namespace std;
+using namespace std;
#include "Mesh2.h"
#include "QuadTree.h"
@@ -49,8 +49,8 @@ namespace bamg {
#ifdef DEBUG1
-extern int SHOW ; // for debugging
-int SHOW = 0; // for debugging
+extern int SHOW ; // for debugging
+int SHOW = 0; // for debugging
#endif
@@ -64,13 +64,13 @@ int hinterpole=1;
long NbUnSwap =0;
int ForDebugging = 0;
const Direction NoDirOfSearch = Direction();
-#ifndef NDEBUG
-inline void MyAssert(int i,char*ex,char * file,long line)
+#ifndef NDEBUG
+inline void MyAssert(int i,char*ex,char * file,long line)
{
if( i) {
cerr << "Error Assert:" << ex << " in " << file << " line: " << line << endl;
#ifdef NOTFREEFEM
- exit(1);
+ exit(1);
#else
throw(ErrorExec("exit",1000));
#endif
@@ -84,23 +84,23 @@ Int4 AGoodNumberPrimeWith(Int4 n)
567890431L, 567890437L, 567890461L, 567890471L,
567890483L, 567890489L, 567890497L, 567890507L,
567890591L, 567890599L, 567890621L, 567890629L , 0};
-
+
Int4 o = 0;
Int4 pi = BigPrimeNumber[1];
for (int i=0; BigPrimeNumber[i]; i++) {
Int4 r = BigPrimeNumber[i] % n;
Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));
- if ( o < oo)
+ if ( o < oo)
o=oo,pi=BigPrimeNumber[i];}
// cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl;
- return pi;
+ return pi;
}
class Triangles;
-void MeshError(int Err,Triangles *Th){
+void MeshError(int Err,Triangles *Th){
cerr << " Fatal error in the meshgenerator " << Err << endl ;
#ifdef NOTFREEFEM
- exit(1);
+ exit(1);
#else
throw(ErrorMesh("Bamg",Err,Th));
#endif
@@ -109,42 +109,42 @@ void MeshError(int Err,Triangles *Th){
ostream& operator <<(ostream& f, const Triangle & ta)
{
if(CurrentTh)
- f << "[" << CurrentTh->Number(ta) << "::"
- << CurrentTh->Number(ta.ns[0]) << ","
- << CurrentTh->Number(ta.ns[1]) << ","
- << CurrentTh->Number(ta.ns[2]) << ","
- << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} "
- << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} "
- << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} "
+ f << "[" << CurrentTh->Number(ta) << "::"
+ << CurrentTh->Number(ta.ns[0]) << ","
+ << CurrentTh->Number(ta.ns[1]) << ","
+ << CurrentTh->Number(ta.ns[2]) << ","
+ << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} "
+ << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} "
+ << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} "
<< "]" ;
else
- f << "["
- << ta.ns[0] << ","
- << ta.ns[1] << ","
- << ta.ns[2] << ","
- << "{" << ta.at[0] << " " << ta.aa[0] << "} "
- << "{" << ta.at[1] << " " << ta.aa[1] << "} "
- << "{" << ta.at[2] << " " << ta.aa[2] << "} "
+ f << "["
+ << ta.ns[0] << ","
+ << ta.ns[1] << ","
+ << ta.ns[2] << ","
+ << "{" << ta.at[0] << " " << ta.aa[0] << "} "
+ << "{" << ta.at[1] << " " << ta.aa[1] << "} "
+ << "{" << ta.at[2] << " " << ta.aa[2] << "} "
<< "]" ;
return f;}
void swap(Triangle *t1,Int1 a1,
Triangle *t2,Int1 a2,
Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2)
-{ // swap
+{ // swap
// --------------------------------------------------------------
// Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles
- //
- // sb sb
+ //
+ // sb sb
// / | \ / \ !
// as1/ | \ /a2 \ !
// / | \ / t2 \ !
// s1 /t1 | t2 \s2 --> s1 /___as2___\s2 !
- // \ a1|a2 / \ as1 /
- // \ | / \ t1 /
- // \ | / as2 \ a1/
- // \ | / \ /
- // sa sa
+ // \ a1|a2 / \ as1 /
+ // \ | / \ t1 /
+ // \ | / as2 \ a1/
+ // \ | / \ /
+ // sa sa
// -------------------------------------------------------------
int as1 = NextEdge[a1];
int as2 = NextEdge[a2];
@@ -161,7 +161,7 @@ void swap(Triangle *t1,Int1 a1,
#endif
(*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb
(*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa
- // mise a jour des 2 adjacences externes
+ // mise a jour des 2 adjacences externes
TriangleAdjacent taas1 = t1->Adj(as1),
taas2 = t2->Adj(as2),
tas1(t1,as1), tas2(t2,as2),
@@ -174,10 +174,10 @@ void swap(Triangle *t1,Int1 a1,
taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());
// externe bas droite
taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());
- // remove the Mark UnMarkSwap
+ // remove the Mark UnMarkSwap
t1->SetUnMarkUnSwap(ap1);
t2->SetUnMarkUnSwap(ap2);
- // interne
+ // interne
tas1.SetAdj2(tas2);
t1->det = det1;
@@ -189,17 +189,17 @@ void swap(Triangle *t1,Int1 a1,
t1->check();
t2->check();
#endif
-#ifdef DRAWING1
+#ifdef DRAWING1
couleur(1);
t1->Draw();
t2->Draw();
-#endif
+#endif
#ifdef DRAWING1
if( CurrentTh)
CurrentTh->inquire();
#endif
-} // end swap
+} // end swap
@@ -209,26 +209,26 @@ Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside)
{
CurrentTh=&Th;
assert(&Th);
- I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y)));
+ I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y)));
Icoor2 dete[3];
Triangle & tb = *Th.FindTriangleContening(I,dete);
-
- if (tb.link)
+
+ if (tb.link)
{ // internal point in a true triangles
a[0]= (Real8) dete[0]/ tb.det;
a[1]= (Real8) dete[1] / tb.det;
a[2] = (Real8) dete[2] / tb.det;
- inside = 1;
+ inside = 1;
return Th.Number(tb);
- }
- else
+ }
+ else
{
- inside = 0;
+ inside = 0;
double aa,bb;
- TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb);
+ TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb);
int k = ta;
Triangle * tc = ta;
- if (!tc->link)
+ if (!tc->link)
{ ta = ta.Adj();
tc=ta;
k = ta;
@@ -244,18 +244,18 @@ Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside)
TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {
-//
- // cout << " - ";
+//
+ // cout << " - ";
int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
int dir=0;
assert(k>=0);
- int kkk=0;
+ int kkk=0;
Icoor2 IJ_IA,IJ_AJ;
- TriangleAdjacent edge(t,OppositeEdge[k]);
- for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge))))
- {
-
- assert(kkk++<1000);
+ TriangleAdjacent edge(t,OppositeEdge[k]);
+ for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge))))
+ {
+
+ assert(kkk++<1000);
Vertex &vI = *edge.EdgeVertex(0);
Vertex &vJ = *edge.EdgeVertex(1);
I2 I=vI, J=vJ, IJ= J-I;
@@ -264,10 +264,10 @@ TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {
if (IJ_IA<0) {
if (dir>0) {a=1;b=0;return edge;}// change of signe => I
else {dir=-1;
- continue;}};// go in direction i
+ continue;}};// go in direction i
IJ_AJ = (IJ ,(J-A));
if (IJ_AJ<0) {
- if(dir<0) {a=0;b=1;return edge;}
+ if(dir<0) {a=0;b=1;return edge;}
else {dir = 1;
continue;}}// go in direction j
double IJ2 = IJ_IA + IJ_AJ;
@@ -276,7 +276,7 @@ TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {
b= IJ_IA/IJ2;
// cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl;
return edge;
- }
+ }
}
TriangleAdjacent Triangle::FindBoundaryEdge(int i) const
@@ -288,53 +288,53 @@ TriangleAdjacent Triangle::FindBoundaryEdge(int i) const
Triangle *t = (Triangle *) this , *ttc;
int k=0,j = EdgesVertexTriangle[i][0],jc;
int exterieur = !link ;
-
- do
+
+ do
{
int exterieurp = exterieur;
- k++;
+ k++;
#ifdef DEBUG
assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] );
#endif
ttc = t->at[j];
exterieur = !ttc->link;
- if (exterieur+exterieurp == 1)
+ if (exterieur+exterieurp == 1)
return TriangleAdjacent(t,j);
jc = NextEdge[t->aa[j]&3];
t = ttc;
j = NextEdge[jc];
assert(k<2000);
- } while ( (this!= t));
+ } while ( (this!= t));
return TriangleAdjacent(0,0);
-
+
}
-TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
-{
- // walk around the vertex
+TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
+{
+ // walk around the vertex
// version 2 for remove the probleme if we fill the hole
//int bug=1;
// Triangle *torigine = t;
// restart:
// int dir=0;
assert(t->link == 0);
- // to have a starting edges
- // try the 3 edge bourna-- in case of internal hole
- // and choice the best
- //
- //
+ // to have a starting edges
+ // try the 3 edge bourna-- in case of internal hole
+ // and choice the best
+ //
+ //
// the probleme is in case of the fine and long internal hole
// for exemple neart the training edge of a wing
- //
+ //
Vertex * s=0,*s1=0, *s0=0;
Icoor2 imax = MaxICoor22;
Icoor2 l0 = imax,l1 = imax;
double dd2 = imax;// infinity
- TriangleAdjacent er;
+ TriangleAdjacent er;
int cas=-2;
for (int j=0;j<3;j++)
- {
+ {
TriangleAdjacent ta=t->FindBoundaryEdge(j);
if (! (Triangle *) ta) continue;
s0 = ta.EdgeVertex(0);
@@ -345,11 +345,11 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
Icoor2 ACAC = (AC,AC), BCBC = (BC,BC);
Icoor2 AB2 = Norme2_2(AB); // ||AB||^2
Icoor2 ABAC = (AB,AC); // AB.AC|
-
+
double d2;
if ( ABAC < 0 ) // DIST A
{
- if ( (d2=(double) ACAC) < dd2)
+ if ( (d2=(double) ACAC) < dd2)
{
// cout << " A " << d2 << " " << dd2;
er = ta;
@@ -361,7 +361,7 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
}
else if (ABAC > AB2) // DIST B
{
- if ( (d2=(double) BCBC) < dd2)
+ if ( (d2=(double) BCBC) < dd2)
{
// cout << " B " << d2 << " " << dd2;
dd2 = d2;
@@ -373,16 +373,16 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
}
}
else // DIST AB
- {
+ {
- double det_2 = (double) Det(AB,AC);
+ double det_2 = (double) Det(AB,AC);
det_2 *= det_2; // square of area*2 of triangle ABC
- d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC
- // cout << " AB " << d2 << " " << dd2
- // << " " << CurrentTh->Number(ta.EdgeVertex(0))
+ d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC
+ // cout << " AB " << d2 << " " << dd2
+ // << " " << CurrentTh->Number(ta.EdgeVertex(0))
// << " " << CurrentTh->Number(ta.EdgeVertex(1)) << " " ;
- if (d2 < dd2)
+ if (d2 < dd2)
{
dd2 = d2;
er = ta;
@@ -401,33 +401,33 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
// l1 = ||C s1|| , l0 = ||C s0||
// where s0,s1 are the vertex of the edge er
- if ( s)
- {
+ if ( s)
+ {
t=er;
- TriangleAdjacent edge(er);
-
- int kkk=0;
+ TriangleAdjacent edge(er);
+
+ int kkk=0;
int linkp = t->link == 0;
-
+
Triangle * tt=t=edge=Adj(Previous(edge));
// cout << CurrentTh->Number(t) << " " << linkp << endl;
do { // loop around vertex s
-
+
assert(edge.EdgeVertex(0)==s && kkk++<10000);
-
+
int link = tt->link == 0;
- // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s)
- // << " " << CurrentTh->Number(er.EdgeVertex(0))
- // << " " << CurrentTh->Number(er.EdgeVertex(1))
- // << " " << CurrentTh->Number(edge.EdgeVertex(0))
- // << " " << CurrentTh->Number(edge.EdgeVertex(1))
+ // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s)
+ // << " " << CurrentTh->Number(er.EdgeVertex(0))
+ // << " " << CurrentTh->Number(er.EdgeVertex(1))
+ // << " " << CurrentTh->Number(edge.EdgeVertex(0))
+ // << " " << CurrentTh->Number(edge.EdgeVertex(1))
// << endl;
- if ((link + linkp) == 1)
- { // a boundary edge
+ if ((link + linkp) == 1)
+ { // a boundary edge
Vertex * st = edge.EdgeVertex(1);
I2 I=*st;
Icoor2 ll = Norme2_2 (C-I);
- if (ll < l1) { // the other vertex is neart
+ if (ll < l1) { // the other vertex is neart
s1=st;
l1=ll;
er = edge;
@@ -443,7 +443,7 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
}
}
}
-
+
linkp=link;
edge=Adj(Previous(edge));
tt = edge;
@@ -455,23 +455,23 @@ TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)
I2 AB=B-A,AC=C-A,CB=B-C;
double aa = (double) (AB,AC);
double bb = (double) (AB,CB);
- // cout << " " << aa << " " << bb
- // << " " << CurrentTh->Number(er.EdgeVertex(0))
+ // cout << " " << aa << " " << bb
+ // << " " << CurrentTh->Number(er.EdgeVertex(0))
// << " " << CurrentTh->Number(er.EdgeVertex(1)) ;
if (aa<0) a=1,b=0;
else if(bb<0) a=0,b=1;
- else
+ else
{
a = bb/(aa+bb);
b = aa/(aa+bb);
}
}
-
- // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " "
- // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a
+
+ // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " "
+ // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a
// << " " << b <<" " << l0 << " " <<l1 <<endl;
return er;
-}
+}
@@ -509,22 +509,22 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
// int OnAVertices =0;
Icoor2 dt[3];
I2 a = Bh.toI2(A) ,b= Bh.toI2(B);// compute the Icoor a,b
- I2 vi,vj;
+ I2 vi,vj;
int iedge =-1;// not a edge
- if(nbegin) {// optimisation
+ if(nbegin) {// optimisation
// we suppose knowing the starting triangle
t=tbegin=lIntTria[ilast=(Size-1)].t;
- if (tbegin->det>=0)
- ifirst = ilast;}
- else {// not optimisation
+ if (tbegin->det>=0)
+ ifirst = ilast;}
+ else {// not optimisation
init();
t=tbegin = Bh.FindTriangleContening(a,deta);
// if(SHOW) cout <<t << " " << Real8(deta[0])/t->det<< " " << Real8(deta[1])/t->det
// << " " << Real8(deta[2])/t->det << endl;
if( t->det>=0)
ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det);
- else
+ else
{// find the nearest boundary edge of the vertex A
// find a edge or such normal projection a the edge IJ is on the edge
// <=> IJ.IA >=0 && IJ.AJ >=0
@@ -543,7 +543,7 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
return;
}
} // find the nearest boundary edge of the vertex A
- } // end not optimisation
+ } // end not optimisation
if (t->det<0) { // outside departure
while (t->det <0) { // intersection boundary edge and a,b,
k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
@@ -557,20 +557,20 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
detj = bamg::det(a,b,vj);
// if(SHOW) { penthickness(3);
// Move(vi);Line(vj);CurrentTh->inquire();penthickness(1);
-// cout << Bh.Number(tbegin) << " " << Bh.Number(t) << " i= " << i <<" j= " << j << " k=" << k
-// << " deti= " << deti << " detj= " << detj
+// cout << Bh.Number(tbegin) << " " << Bh.Number(t) << " i= " << i <<" j= " << j << " k=" << k
+// << " deti= " << deti << " detj= " << detj
// << " v = " << Bh.Number((*t)[i]) << (*t)[i].r << " " << Bh.Number((*t)[j]) << (*t)[j].r << endl;}
if (deti>0) // go to i direction on gamma
- ocut = PreviousEdge[ocut];
+ ocut = PreviousEdge[ocut];
else if (detj<=0) // go to j direction on gamma
- ocut = NextEdge[ocut];
+ ocut = NextEdge[ocut];
TriangleAdjacent tadj =t->Adj(ocut);
t = tadj;
- iedge= tadj;
- if (t == tbegin) { //
+ iedge= tadj;
+ if (t == tbegin) { //
double ba,bb;
- if (verbosity>7)
- cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b)
+ if (verbosity>7)
+ cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b)
<< " deti= " << deti << " detj=" <<detj << endl;
TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);
Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
@@ -594,12 +594,12 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
Line(B);
penthickness(1);
Bh.inquire();
-#endif
+#endif
MeshError(997);*/
}
} // end while (t->det <0)
// theoriticaly we have: deti =<0 and detj>0
-
+
// computation of barycentric coor
// test if the point b is on size on t
// we revert vi,vj because vi,vj is def in Adj triangle
@@ -626,14 +626,14 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
ilast=NewItem(t,ba[0],ba[1],ba[2]); }
} // outside departure
-
-
+
+
// recherche the intersection of [a,b] with Bh Mesh.
// we know a triangle ta contening the vertex a
// we have 2 case for intersection [a,b] with a edge [A,B] of Bh
// 1) the intersection point is in ]A,B[
// 2) is A or B
- // first version ---
+ // first version ---
for (;;) {
// t->Draw();
if (iedge < 0) {
@@ -644,11 +644,11 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
else {
i2 = iedge;
i0 = NextEdge[i2];
- i1 = NextEdge[i0];
+ i1 = NextEdge[i0];
dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because
dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side
dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);}
-
+
// so we have just to see the transition from - to + of the det0..2 on edge of t
// because we are going from a to b
if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
@@ -657,7 +657,7 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
(dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
ocut =i1;
- else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
+ else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
(dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
ocut =i2;
else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) &&
@@ -666,7 +666,7 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) &&
(dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
ocut =i1;
- else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&
+ else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&
(dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
ocut =i2;
else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
@@ -675,13 +675,13 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
(dt[j=VerticesOfTriangularEdge[i1][1]] == 0))
ocut =i1;
- else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
+ else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
(dt[j=VerticesOfTriangularEdge[i2][1]] == 0))
ocut =i2;
else { // On a edge (2 zero)
k =0;
- if (dt[0]) ocut=0,k++;
- if (dt[1]) ocut=1,k++;
+ if (dt[0]) ocut=0,k++;
+ if (dt[1]) ocut=1,k++;
if (dt[2]) ocut=2,k++;
if(k == 1) {
if (dt[ocut] >0) // triangle upper AB
@@ -691,8 +691,8 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
}
else {
cerr << " Bug Split Edge " << endl;
- cerr << " dt[0]= " << dt[0]
- << " dt[1]= " << dt[1]
+ cerr << " dt[0]= " << dt[0]
+ << " dt[1]= " << dt[1]
<< " dt[2]= "<< dt[2] << endl;
cerr << i0 << " " << i1 << " " << i2 << endl;
cerr << " A = " << A << " B= " << B << endl;
@@ -700,36 +700,36 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl;
cerr << " nbt = " << nbt << endl;
MeshError(100);}}
-
+
k = OppositeVertex[ocut];
Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);
-
+
if (detbij >= 0) { //we find the triangle contening b
dt[0]=bamg::det((*t)[1],(*t)[2],b);
dt[1]=bamg::det((*t)[2],(*t)[0],b);
dt[2]=bamg::det((*t)[0],(*t)[1],b);
-#ifdef DEBUG
+#ifdef DEBUG
assert(dt[0] >= 0);
assert(dt[1] >= 0);
assert(dt[2] >= 0);
#endif
Real8 dd = t->det;
- NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);
+ NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);
return ;}
- else { // next triangle by adjacent by edge ocut
+ else { // next triangle by adjacent by edge ocut
deti = dt[i];
detj = dt[j];
Real4 dij = detj-deti;
ba[i] = detj/dij;
ba[j] = -deti/dij;
ba[3-i-j ] = 0;
- ilast=NewItem(t, ba[0],ba[1],ba[2]);
-
+ ilast=NewItem(t, ba[0],ba[1],ba[2]);
+
TriangleAdjacent ta =t->Adj(ocut);
t = ta;
- iedge= ta;
+ iedge= ta;
if (t->det <= 0) {
double ba,bb;
TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);
@@ -738,20 +738,20 @@ void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,
// ajoute le 03 frev 1997 par F. hecht
return;
}
- }// we go outside of omega
+ }// we go outside of omega
} // for(;;)
-
-
+
+
} // routine SplitEdge
-int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) {
+int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) {
register int n;
R2 x(0,0);
if ( d0) x = (*tt)[0].r * d0;
if ( d1) x = x + (*tt)[1].r * d1;
if ( d2) x = x + (*tt)[2].r * d2;
- // newer add same point
+ // newer add same point
if(!Size || Norme2_2(lIntTria[Size-1].x-x)) {
if (Size==MaxSize) ReShape();
lIntTria[Size].t=tt;
@@ -773,7 +773,7 @@ int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8
else n=Size-1;
return n;
}
-int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
+int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
register int n;
if(!Size || Norme2_2(lIntTria[Size-1].x-A)) {
if (Size==MaxSize) ReShape();
@@ -786,32 +786,32 @@ int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {
n=Size++;
}
else n=Size-1;
- return n;
+ return n;
}
Real8 ListofIntersectionTriangles::Length()
{
// cout << " n= " << Size << ":" ;
assert(Size>0);
- // computation of the length
+ // computation of the length
R2 C;
Metric Mx,My;
int ii,jj;
R2 x,y,xy;
-
+
SegInterpolation *SegI=lSegsI;
SegI=lSegsI;
- lSegsI[NbSeg].last=Size;// improvement
-
+ lSegsI[NbSeg].last=Size;// improvement
+
int EndSeg=Size;
-
+
y = lIntTria[0].x;
Real8 sxy, s = 0;
lIntTria[0].s =0;
SegI->lBegin=s;
- for (jj=0,ii=1;ii<Size;jj=ii++)
- {
+ for (jj=0,ii=1;ii<Size;jj=ii++)
+ {
// seg jj,ii
x=y;
y = lIntTria[ii].x;
@@ -826,12 +826,12 @@ Real8 ListofIntersectionTriangles::Length()
sxy = LengthInterpole(Mx,My,xy);
s += sxy;
lIntTria[ii].s = s;
- if (ii == EndSeg)
+ if (ii == EndSeg)
SegI->lEnd=s,
SegI++,
EndSeg=SegI->last,
SegI->lBegin=s;
-
+
// cout << ii << " " << jj << x<< y <<xy << s << lIntTria[ii].m ;
}
len = s;
@@ -847,7 +847,7 @@ Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 n
const Int4 nbvold = nbv;
Real8 s = Length();
if (s < 1.5 ) return 0;
- //////////////////////
+ //////////////////////
int ii = 1 ;
R2 y,x;
Metric My,Mx ;
@@ -858,35 +858,35 @@ Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 n
int EndSeg=Size;
SegInterpolation *SegI=0;
- if (NbSeg)
+ if (NbSeg)
SegI=lSegsI,EndSeg=SegI->last;
-
+
for (int k=1;k<nbi;k++)
{
- while ((ii < Size) && ( lIntTria[ii].s <= si ))
- if (ii++ == EndSeg)
+ while ((ii < Size) && ( lIntTria[ii].s <= si ))
+ if (ii++ == EndSeg)
SegI++,EndSeg=SegI->last;
int ii1=ii-1;
x =lIntTria[ii1].x;
sx =lIntTria[ii1].s;
Metric Mx=lIntTria[ii1].m;
-#ifdef DEBUG
+#ifdef DEBUG
double lx = lIntTria[ii-1].sn;
#endif
y =lIntTria[ii].x;
sy =lIntTria[ii].s;
Metric My=lIntTria[ii].m;
-#ifdef DEBUG
- double ly =lIntTria[ii].sp;
+#ifdef DEBUG
+ double ly =lIntTria[ii].sp;
assert( sx <= si);
assert( si <= sy);
assert( sy != sx);
-#endif
+#endif
Real8 lxy = sy-sx;
Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
-
+
R2 C;
Real8 cx = 1-cy;
C = SegI ? SegI->F(si): x * cx + y *cy;
@@ -917,36 +917,36 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap)
{ // l'arete ta coupe l'arete pva pvb
// de cas apres le swap sa coupe toujours
- // on cherche l'arete suivante
+ // on cherche l'arete suivante
// on suppose que detsa >0 et detsb <0
- // attention la routine echange pva et pvb
+ // attention la routine echange pva et pvb
- if(tt1.Locked()) return 0; // frontiere croise
+ if(tt1.Locked()) return 0; // frontiere croise
TriangleAdjacent tt2 = Adj(tt1);
Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent
Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles
assert ( a1 >= 0 && a1 < 3 );
-
+
Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];
Vertex & s1= (*t1)[OppositeVertex[a1]];
Vertex & s2= (*t2)[OppositeVertex[a2]];
-
+
Icoor2 dets2 = det(*pva,*pvb,s2);
#ifdef DEBUG
Vertex & sb= (*t1)[VerticesOfTriangularEdge[a1][1]];
- Icoor2 wdets1 = det(*pva,*pvb,s1);
+ Icoor2 wdets1 = det(*pva,*pvb,s1);
Icoor2 wdetsa = det(*pva,*pvb,sa);
Icoor2 wdetsb = det(*pva,*pvb,sb);
assert(wdets1 == dets1);
assert(wdetsa == detsa);
assert(wdetsb == detsb);
#endif
-
+
Icoor2 det1=t1->det , det2=t2->det ;
-#ifdef DEBUG
+#ifdef DEBUG
assert(det1>0 && det2 >0);
Icoor2 ddet1 = det((*t1)[0],(*t1)[1],(*t1)[2]);
Icoor2 ddet2 = det((*t2)[0],(*t2)[1],(*t2)[2]);
@@ -957,9 +957,9 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
}
Icoor2 detvasasb = det(*pva,sa,sb);
Icoor2 detvbsasb = det(*pvb,sa,sb);
- if ( CurrentTh && ! ( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0))))
+ if ( CurrentTh && ! ( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0))))
{
- cout << " detvasasb =" << detvasasb << "detvbsasb = " << detvbsasb
+ cout << " detvasasb =" << detvasasb << "detvbsasb = " << detvbsasb
<< " " << pva << " " << pvb << " " <<CurrentTh <<endl;
#ifdef DRAWING1
reffecran();
@@ -968,7 +968,7 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
pva->MoveTo();pvb->LineTo();
penthickness(1);
CurrentTh->inquire();
-#endif
+#endif
}
assert( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0)));
#endif
@@ -980,18 +980,18 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
Icoor2 ndet2 = detT - ndet1;
int ToSwap =0; //pas de swap
- if ((ndet1 >0) && (ndet2 >0))
- { // on peut swaper
+ if ((ndet1 >0) && (ndet2 >0))
+ { // on peut swaper
if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
- ToSwap =1;
- else // swap alleatoire
- if (BinaryRand())
- ToSwap =2;
+ ToSwap =1;
+ else // swap alleatoire
+ if (BinaryRand())
+ ToSwap =2;
}
#ifdef DEBUG
if (ForDebugging) {
- cerr << "swap = " << ToSwap << " ndet1 " << ndet1 << ", ndet2 " << ndet2 << "det1 " << det1 << " det2 " << det2
- << " if1 = " << ((ndet1 >0) && (ndet2 >0))
+ cerr << "swap = " << ToSwap << " ndet1 " << ndet1 << ", ndet2 " << ndet2 << "det1 " << det1 << " det2 " << det2
+ << " if1 = " << ((ndet1 >0) && (ndet2 >0))
<< " if2 = " << ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) << endl;
#ifdef DRAWING
couleur(0);
@@ -1002,7 +1002,7 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
#endif
if (ToSwap) NbSwap++,
bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);
-
+
#ifdef DEBUG
if (ForDebugging) {
#ifdef DRAWING
@@ -1017,15 +1017,15 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
if (dets2 < 0) {// haut
dets1 = ToSwap ? dets1 : detsa ;
- detsa = dets2;
+ detsa = dets2;
tt1 = Previous(tt2) ;}
- else if (dets2 > 0){// bas
+ else if (dets2 > 0){// bas
dets1 = ToSwap ? dets1 : detsb ;
detsb = dets2;
//xxxx tt1 = ToSwap ? tt1 : Next(tt2);
if(!ToSwap) tt1 = Next(tt2);
}
- else { // changement de sens
+ else { // changement de sens
if (ForDebugging) cout << "changement de sens" << endl;
ret = -1;
Exchange(pva,pvb);
@@ -1037,12 +1037,12 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
detsa=-detsa;
detsb=-detsb;
- if (ToSwap)
+ if (ToSwap)
if (dets2 < 0) {// haut
dets1 = (ToSwap ? dets1 : detsa) ;
- detsa = dets2;
+ detsa = dets2;
tt1 = Previous(tt2) ;}
- else if (dets2 > 0){// bas
+ else if (dets2 > 0){// bas
dets1 = (ToSwap ? dets1 : detsb) ;
detsb = dets2;
if(!ToSwap) tt1 = Next(tt2);
@@ -1055,22 +1055,22 @@ int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,
return ret;
}
-int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
-{
+int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
+{
#ifdef DEBUG
- restart: // for debug
+ restart: // for debug
#endif
int NbSwap =0;
- assert(a.t && b.t); // the 2 vertex is in a mesh
+ assert(a.t && b.t); // the 2 vertex is in a mesh
int k=0;
- taret=TriangleAdjacent(0,0); // erreur
-
+ taret=TriangleAdjacent(0,0); // erreur
+
TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);
Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
- // we turn around a in the direct sens
-
+ // we turn around a in the direct sens
+
Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;
- if(v2) // normal case
+ if(v2) // normal case
det2 = det(*v2,a,b);
else { // no chance infini vertex try the next
tta= Previous(Adj(tta));
@@ -1086,18 +1086,18 @@ int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
#endif
while (v2 != &b) {
- TriangleAdjacent tc = Previous(Adj(tta));
- v1 = v2;
+ TriangleAdjacent tc = Previous(Adj(tta));
+ v1 = v2;
v2 = tc.EdgeVertex(0);
det1 = det2;
#ifdef DEBUG
assert( v1 == tta.EdgeVertex(0));
assert( &a == tc.EdgeVertex(1) );
#endif
- det2 = v2 ? det(*v2,a,b): det2;
-
- if((det1 < 0) && (det2 >0)) {
- // try to force the edge
+ det2 = v2 ? det(*v2,a,b): det2;
+
+ if((det1 < 0) && (det2 >0)) {
+ // try to force the edge
Vertex * va = &a, *vb = &b;
tc = Previous(tc);
assert ( v1 && v2);
@@ -1113,17 +1113,17 @@ int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
ForDebugging=1;
goto restart;
}
-
-#endif
+
+#endif
while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))
if(l++ > 10000000) {
- cerr << " Loop in forcing Egde AB"
- <<"\n vertex A " << a
- <<"\n vertex B " << b
- <<"\n nb de swap " << NbSwap
+ cerr << " Loop in forcing Egde AB"
+ <<"\n vertex A " << a
+ <<"\n vertex B " << b
+ <<"\n nb de swap " << NbSwap
<<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl;
-
- if ( CurrentTh )
+
+ if ( CurrentTh )
cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl;
#ifdef DEBUG
ForDebugging = 1;
@@ -1144,7 +1144,7 @@ int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
while (ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap) && (l++ < 1000))
cerr << " " << CurrentTh->Number(tc.EdgeVertex(0))<<" " <<CurrentTh->Number(tc.EdgeVertex(1)) << " ";
}
-#endif
+#endif
MeshError(990);
}
Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
@@ -1155,7 +1155,7 @@ int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
taret = tc;
return NbSwap;
}
- else
+ else
{
taret = tc;
return -2; // error boundary is crossing
@@ -1166,25 +1166,25 @@ int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)
cerr << CurrentTh->Number(aa) << " " << CurrentTh->Number(bb) << " ] " << endl;
}
MeshError(991);
-*/
+*/
}
}
tta = tc;
assert(k++<2000);
- if ( vbegin == v2 ) return -1;// error
+ if ( vbegin == v2 ) return -1;// error
}
tta.SetLock();
taret=tta;
a.Optim(1,0);
b.Optim(1,0);
- return NbSwap;
+ return NbSwap;
}
int Triangle::swap(Int2 a,int koption){
#ifdef DEBUG
- if(a &4 ) return 0;// arete lock
+ if(a &4 ) return 0;// arete lock
int munswap1 = a/4;
a &=3;
#else
@@ -1200,22 +1200,22 @@ int Triangle::swap(Int2 a,int koption){
#else
if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
#endif
-
+
register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]];
register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]];
register Vertex *s1=t1->ns[OppositeVertex[a1]];
register Vertex *s2=t2->ns[OppositeVertex[a2]];
#ifdef DEBUG
- assert ( a >= 0 && a < 3 );
+ assert ( a >= 0 && a < 3 );
#endif
-
+
Icoor2 det1=t1->det , det2=t2->det ;
Icoor2 detT = det1+det2;
Icoor2 detA = Abs(det1) + Abs(det2);
Icoor2 detMin = Min(det1,det2);
- int OnSwap = 0;
+ int OnSwap = 0;
// si 2 triangle infini (bord) => detT = -2;
if (sa == 0) {// les deux triangles sont frontieres
det2=bamg::det(s2->i,sb->i,s1->i);
@@ -1227,9 +1227,9 @@ int Triangle::swap(Int2 a,int koption){
det1 = bamg::det(s1->i,sa->i,s2->i);
det2 = detT - det1;
OnSwap = (Abs(det1) + Abs(det2)) < detA;
-
+
Icoor2 detMinNew=Min(det1,det2);
- // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test
+ // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test
if (! OnSwap &&(detMinNew>0)) {
OnSwap = detMin ==0;
if (! OnSwap) {
@@ -1241,7 +1241,7 @@ int Triangle::swap(Int2 a,int koption){
x21 = s2->i.x - s1->i.x,
yb1 = sb->i.y - s1->i.y,
y21 = s2->i.y - s1->i.y,
- xba = sb->i.x - sa->i.x,
+ xba = sb->i.x - sa->i.x,
x2a = s2->i.x - sa->i.x,
yba = sb->i.y - sa->i.y,
y2a = s2->i.y - sa->i.y;
@@ -1251,18 +1251,18 @@ int Triangle::swap(Int2 a,int koption){
sinb12 = double(det2),
sinba2 = double(t2->det);
-
+
// angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)
OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);
-// if(CurrentTh)
+// if(CurrentTh)
// cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ;
-// cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12
+// cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb12
// << " Onswap = " << OnSwap << endl;
break;
}
- else
- {
- // critere de Delaunay anisotrope
+ else
+ {
+ // critere de Delaunay anisotrope
Real8 som;
I2 AB=(I2) *sb - (I2) *sa;
I2 MAB2=((I2) *sb + (I2) *sa);
@@ -1275,81 +1275,81 @@ int Triangle::swap(Int2 a,int koption){
Metric M=s1->m;
R2 ABo = M.Orthogonal(AB);
R2 A1o = M.Orthogonal(A1);
- // (A+B)+ x ABo = (S1+B)/2+ y A1
+ // (A+B)+ x ABo = (S1+B)/2+ y A1
// ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
if (Abs(d) > dd*1.e-3) {
R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
som = M(C - S2)/M(C - S1);
- } else
+ } else
{kopt=1;continue;}
-
+
}
{
Metric M=s2->m;
R2 ABo = M.Orthogonal(AB);
R2 A1o = M.Orthogonal(A1);
- // (A+B)+ x ABo = (S1+B)/2+ y A1
- // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
+ // (A+B)+ x ABo = (S1+B)/2+ y A1
+ // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
if(Abs(d) > dd*1.e-3) {
R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
som += M(C - S2)/M(C - S1);
- } else
+ } else
{kopt=1;continue;}
}
OnSwap = som < 2;
break;
}
-
- } // OnSwap
+
+ } // OnSwap
} // (! OnSwap &&(det1 > 0) && (det2 > 0) )
}
-#ifdef DEBUG1
+#ifdef DEBUG1
if (OnSwap && ( munswap1 || munswap2)) {
cout << " erreur Mark unswap T " << CurrentTh->Number(t1) << " " << CurrentTh->Number(t2) << endl
<< *t1 << endl
<< *t2 << endl;
return 0;
}
-#endif
- if( OnSwap )
+#endif
+ if( OnSwap )
bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);
else {
NbUnSwap ++;
- t1->SetMarkUnSwap(a1);
+ t1->SetMarkUnSwap(a1);
}
return OnSwap;
}
Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega)
{
-#ifdef DEBUG
+#ifdef DEBUG
register Int4 NbSwap =0;
#endif
register Vertex * s = this;
Vertex &vP = *s,vPsave=vP;
- // if (vP.on) return 0;// Don't move boundary vertex
-
+ // if (vP.on) return 0;// Don't move boundary vertex
+
register Triangle * tbegin= t , *tria = t , *ttc;
-
+
register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc;
R2 P(s->r),PNew(0,0);
// cout << BTh.quadtree << " " << BTh.quadtree->root << endl;
// assert(BTh.quadtree && BTh.quadtree->root);
do {
- k++;
-
+ k++;
+
#ifdef DEBUG
assert( s == & (*tria)[VerticesOfTriangularEdge[j][1]] );
assert( tria->det >0);
#endif
if (!tria->Hidden(j))
{
- Vertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
-
+ Vertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
+
R2 Q = vQ,QP(P-Q);
Real8 lQP = LengthInterpole(vP,vQ,QP);
PNew += Q+QP/Max(lQP,1e-20);
@@ -1360,28 +1360,28 @@ Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tsta
tria = ttc;
j = NextEdge[jc];
assert(k<2000);
- } while ( tbegin != tria);
+ } while ( tbegin != tria);
if (kk<4) return 0;
PNew = PNew/(Real8)kk;
R2 Xmove((PNew-P)*omega);
PNew = P+Xmove;
- Real8 delta=Norme2_2(Xmove);
-
-
- //
+ Real8 delta=Norme2_2(Xmove);
+
+
+ //
Icoor2 deta[3];
I2 IBTh = BTh.toI2(PNew);
-
- tstart=BTh.FindTriangleContening(IBTh,deta,tstart);
-
- if (tstart->det <0)
+
+ tstart=BTh.FindTriangleContening(IBTh,deta,tstart);
+
+ if (tstart->det <0)
{ // outside
double ba,bb;
TriangleAdjacent edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
tstart = edge;
vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));
}
- else
+ else
{ // inside
Real8 aa[3];
Real8 s = deta[0]+deta[1]+deta[2];
@@ -1390,28 +1390,28 @@ Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tsta
aa[2]=deta[2]/s;
vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]);
}
-
+
// recompute the det of the triangle
vP.r = PNew;
-
+
vP.i = Th.toI2(PNew);
Vertex vPnew = vP;
-
+
int ok=1;
int loop=1;
k=0;
- while (ok)
+ while (ok)
{
ok =0;
do {
- k++;
+ k++;
double detold = tria->det;
tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]);
- if (loop)
+ if (loop)
{
Vertex *v0,*v1,*v2,*v3;
- if (tria->det<0) ok =1;
+ if (tria->det<0) ok =1;
else if (tria->Quadrangle(v0,v1,v2,v3))
{
vP = vPsave;
@@ -1421,7 +1421,7 @@ Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tsta
if (qnew<qold) ok = 1;
}
else if ( (double)tria->det < detold/2 ) ok=1;
-
+
}
tria->SetUnMarkUnSwap(0);
tria->SetUnMarkUnSwap(1);
@@ -1431,8 +1431,8 @@ Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tsta
tria = ttc;
j = NextEdge[jc];
assert(k<2000);
- } while ( tbegin != tria);
- if (ok && loop) vP=vPsave; // no move
+ } while ( tbegin != tria);
+ if (ok && loop) vP=vPsave; // no move
loop=0;
}
return delta;
@@ -1441,7 +1441,7 @@ Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tsta
-void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
+void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
{
// -------------------------------------------
// s2
@@ -1456,8 +1456,8 @@ void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
// / . ` \ !
// ---------------- !
// s0 tt2 s1
- //--------------------------------------------
-
+ //--------------------------------------------
+
Triangle * tt[3]; // the 3 new Triangles
Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2];
Icoor2 det3local[3];
@@ -1466,52 +1466,52 @@ void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
register int nbd0 =0; // number of zero det3
register int izerodet=-1,iedge; // izerodet = egde contening the vertex s
Icoor2 detOld = t->det;
-
- if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) )
+
+ if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) )
{
cerr << " infv " << infv << " det = " << detOld << endl;
- cerr << Number(s) << " "<< Number(s0) << " "
+ cerr << Number(s) << " "<< Number(s0) << " "
<< Number(s1) << " " << Number(s2) << endl;
MeshError(3);
}
-
+
// if det3 do not exist then constuct det3
- if (!det3) {
- det3 = det3local; // alloc
+ if (!det3) {
+ det3 = det3local; // alloc
if ( infv<0 ) {
det3[0]=bamg::det(s ,s1,s2);
det3[1]=bamg::det(s0,s ,s2);
det3[2]=bamg::det(s0,s1,s );}
- else {
+ else {
// one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk
det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ;
det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ;
det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}}
-
+
if (!det3[0]) izerodet=0,nbd0++;
if (!det3[1]) izerodet=1,nbd0++;
if (!det3[2]) izerodet=2,nbd0++;
-
- if (nbd0 >0 ) // point s on a egde or on a vertex
+
+ if (nbd0 >0 ) // point s on a egde or on a vertex
if (nbd0 == 1) {
iedge = OppositeEdge[izerodet];
TriangleAdjacent ta = t->Adj(iedge);
-#ifdef DEBUG1
- cout << " the point " << Number(s) << " is the edge " << izerodet
- << " of " << Number(t) << " det3 = "
+#ifdef DEBUG1
+ cout << " the point " << Number(s) << " is the edge " << izerodet
+ << " of " << Number(t) << " det3 = "
<< det3[0] << " " << det3[1] << " " << det3[2] << " " << endl;
- cout << " ta = " << ta << "ta->det =" << ((Triangle*) ta)->det
+ cout << " ta = " << ta << "ta->det =" << ((Triangle*) ta)->det
<< " "<< t->det<< endl;
#endif
- // the point is on the edge
- // if the point is one the boundary
- // add the point in outside part
+ // the point is on the edge
+ // if the point is one the boundary
+ // add the point in outside part
if ( t->det >=0) { // inside triangle
if ((( Triangle *) ta)->det < 0 ) {
- // add in outside triangle
+ // add in outside triangle
Add(s,( Triangle *) ta);
return;}
}}
@@ -1525,14 +1525,14 @@ void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
MeshError(5,this);}
// remove de MarkUnSwap edge
- t->SetUnMarkUnSwap(0);
- t->SetUnMarkUnSwap(1);
+ t->SetUnMarkUnSwap(0);
+ t->SetUnMarkUnSwap(1);
t->SetUnMarkUnSwap(2);
tt[0]= t;
tt[1]= &triangles[nbt++];
tt[2]= &triangles[nbt++];
-
+
if (nbt>nbtx) {
cerr << " No enougth triangles " << endl;
MeshError(999,this);
@@ -1541,17 +1541,17 @@ void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
*tt[1]= *tt[2]= *t;
// gestion of the link
tt[0]->link=tt[1];
- tt[1]->link=tt[2];
-
+ tt[1]->link=tt[2];
+
(* tt[0])(OppositeVertex[0])=&s;
(* tt[1])(OppositeVertex[1])=&s;
(* tt[2])(OppositeVertex[2])=&s;
tt[0]->det=det3[0];
tt[1]->det=det3[1];
- tt[2]->det=det3[2];
+ tt[2]->det=det3[2];
- // update adj des triangles externe
+ // update adj des triangles externe
tt[0]->SetAdjAdj(0);
tt[1]->SetAdjAdj(1);
tt[2]->SetAdjAdj(2);
@@ -1563,53 +1563,53 @@ void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3)
tt[i0]->SetAdj2(i2,tt[i2],i0);
tt[i1]->SetAdj2(i0,tt[i0],i1);
tt[i2]->SetAdj2(i1,tt[i1],i2);
-
+
tt[0]->SetTriangleContainingTheVertex();
tt[1]->SetTriangleContainingTheVertex();
tt[2]->SetTriangleContainingTheVertex();
-
-
+
+
// swap if the point s is on a edge
if(izerodet>=0) {
// cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;
int rswap =tt[izerodet]->swap(iedge);
-
- if (!rswap)
+
+ if (!rswap)
{
cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;
#ifdef DRAWING
- if( CurrentTh && withrgraphique)
+ if( CurrentTh && withrgraphique)
{
reffecran();
-
+
DrawMark(s.r);
CurrentTh->inquire();
DrawMark(s.r);
- rattente(1);
- }
+ rattente(1);
+ }
#endif
}
assert(rswap);
}
-
-#ifdef DEBUG
+
+#ifdef DEBUG
tt[0]->check();
tt[1]->check();
tt[2]->check();
#endif
-#ifdef DRAWING1
+#ifdef DRAWING1
tt[0]->Draw();
tt[1]->Draw();
tt[2]->Draw();
#endif
-
+
}
Int4 Triangles::SplitInternalEdgeWithBorderVertices()
{
Int4 NbSplitEdge=0;
- SetVertexFieldOn();
+ SetVertexFieldOn();
Int4 it;
Int4 nbvold=nbv;
for (it=0;it<nbt;it++)
@@ -1619,8 +1619,8 @@ Int4 Triangles::SplitInternalEdgeWithBorderVertices()
for (int j=0;j<3;j++)
if(!t.Locked(j) && !t.Hidden(j)){
Triangle &tt = *t.TriangleAdj(j);
- if ( &tt && tt.link && it < Number(tt))
- { // an internal edge
+ if ( &tt && tt.link && it < Number(tt))
+ { // an internal edge
Vertex &v0 = t[VerticesOfTriangularEdge[j][0]];
Vertex &v1 = t[VerticesOfTriangularEdge[j][1]];
if (v0.on && v1.on)
@@ -1634,32 +1634,32 @@ Int4 Triangles::SplitInternalEdgeWithBorderVertices()
}
NbSplitEdge++;
if (verbosity>7)
- cout <<" Internal edge with two vertices on boundary"
+ cout <<" Internal edge with two vertices on boundary"
<< Number(v0) << " " << Number(v1) << " by " << endl;
}
}
}
}
- ReMakeTriangleContainingTheVertex();
- if (nbvold!=nbv)
+ ReMakeTriangleContainingTheVertex();
+ if (nbvold!=nbv)
{
Int4 iv = nbvold;
Int4 NbSwap = 0;
- Icoor2 dete[3];
- for (Int4 i=nbvold;i<nbv;i++)
+ Icoor2 dete[3];
+ for (Int4 i=nbvold;i<nbv;i++)
{// for all the new point
Vertex & vi = vertices[i];
vi.i = toI2(vi.r);
vi.r = toR2(vi.i);
// if (!quadtree->ToClose(vi,seuil,hi,hj)) {
- // a good new point
- vi.ReferenceNumber=0;
+ // a good new point
+ vi.ReferenceNumber=0;
vi.DirOfSearch =NoDirOfSearch;
- // cout << " Add " << Number(vi) << " " << vi
+ // cout << " Add " << Number(vi) << " " << vi
// << " " << Number(vi) << " <--> " << Number(vi) <<endl;
Triangle *tcvi = FindTriangleContening(vi.i,dete);
if (tcvi && !tcvi->link) {
- cout << i << " PB insert point " << Number(vi) << vi << Number(vi)
+ cout << i << " PB insert point " << Number(vi) << vi << Number(vi)
<< " tcvi = " << tcvi << " " << tcvi->link << endl;
cout << (*tcvi)[1] << (*tcvi)[2] << endl;
tcvi = FindTriangleContening(vi.i,dete);
@@ -1671,26 +1671,26 @@ Int4 Triangles::SplitInternalEdgeWithBorderVertices()
penthickness(1);
inquire();
#endif
-
+
MeshError(1001,this);
}
-
-
+
+
quadtree->Add(vi);
#ifdef DRAWING1
DrawMark(vi.r);
#endif
- assert (tcvi && tcvi->det >= 0) ;// internal
+ assert (tcvi && tcvi->det >= 0) ;// internal
Add(vi,tcvi,dete);
- NbSwap += vi.Optim(1);
+ NbSwap += vi.Optim(1);
iv++;
// }
}
- if (verbosity>3)
+ if (verbosity>3)
{
cout << " Nb Of New Point " << iv ;
cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;}
-
+
nbv = iv;
}
if (NbSplitEdge > nbv-nbvold)
@@ -1702,31 +1702,31 @@ Int4 Triangles::SplitInternalEdgeWithBorderVertices()
}
Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
{
- Real8 seuil= 1.414/2 ;// for two close point
+ Real8 seuil= 1.414/2 ;// for two close point
Int4 i;
- // insertion part ---
+ // insertion part ---
const Int4 nbvnew = nbv-nbvold;
- if (verbosity>5)
- cout << " Try to Insert the " <<nbvnew<< " new points " << endl;
+ if (verbosity>5)
+ cout << " Try to Insert the " <<nbvnew<< " new points " << endl;
Int4 NbSwap=0;
Icoor2 dete[3];
-
- // construction d'un ordre aleatoire
- if (! nbvnew)
- return 0;
+
+ // construction d'un ordre aleatoire
+ if (! nbvnew)
+ return 0;
if (nbvnew) {
const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
- Int4 k3 = rand()%nbvnew ;
+ Int4 k3 = rand()%nbvnew ;
for (Int4 is3=0; is3<nbvnew; is3++) {
register Int4 j = nbvold +(k3 = (k3 + PrimeNumber)% nbvnew);
- register Int4 i = nbvold+is3;
+ register Int4 i = nbvold+is3;
ordre[i]= vertices + j;
ordre[i]->ReferenceNumber=i;
}
- // be carefull
+ // be carefull
Int4 iv = nbvold;
- for (i=nbvold;i<nbv;i++)
+ for (i=nbvold;i<nbv;i++)
{// for all the new point
Vertex & vi = *ordre[i];
vi.i = toI2(vi.r);
@@ -1734,24 +1734,24 @@ Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
Real4 hx,hy;
vi.m.Box(hx,hy);
Icoor1 hi=(Icoor1) (hx*coefIcoor),hj=(Icoor1) (hy*coefIcoor);
- if (!quadtree->ToClose(vi,seuil,hi,hj))
+ if (!quadtree->ToClose(vi,seuil,hi,hj))
{
- // a good new point
+ // a good new point
Vertex & vj = vertices[iv];
- Int4 j = vj.ReferenceNumber;
+ Int4 j = vj.ReferenceNumber;
assert( &vj== ordre[j]);
if(i!=j)
- { // for valgring
+ { // for valgring
Exchange(vi,vj);
Exchange(ordre[j],ordre[i]);
}
- vj.ReferenceNumber=0;
- // cout << " Add " << Number(vj) << " " << vj
+ vj.ReferenceNumber=0;
+ // cout << " Add " << Number(vj) << " " << vj
// << " " << Number(vi) << " <--> " << Number(vj) <<endl;
Triangle *tcvj = FindTriangleContening(vj.i,dete);
- if (tcvj && !tcvj->link)
+ if (tcvj && !tcvj->link)
{
- cerr << i << " PB insert point " << Number(vj) << vj << Number(vi)
+ cerr << i << " PB insert point " << Number(vj) << vj << Number(vi)
<< " tcvj = " << tcvj << " " << tcvj->link << endl;
cerr << (*tcvj)[1] << (*tcvj)[2] << endl;
tcvj = FindTriangleContening(vj.i,dete);
@@ -1761,24 +1761,24 @@ Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
penthickness(5);
DrawMark(vj.r);
penthickness(1);
-
+
inquire();
#endif
-
+
MeshError(1001,this);
}
-
-
+
+
quadtree->Add(vj);
#ifdef DRAWING1
DrawMark(vj.r);
#endif
- assert (tcvj && tcvj->det >= 0) ;// internal
+ assert (tcvj && tcvj->det >= 0) ;// internal
Add(vj,tcvj,dete);
- NbSwap += vj.Optim(1);
+ NbSwap += vj.Optim(1);
iv++;
}
- }
+ }
if (verbosity>3) {
cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ;
cout << " Nb swap = " << NbSwap << " after " ;}
@@ -1790,37 +1790,37 @@ Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
inquire();
#endif
- for (i=nbvold;i<nbv;i++)
- NbSwap += vertices[i].Optim(1);
- if (verbosity>3)
+ for (i=nbvold;i<nbv;i++)
+ NbSwap += vertices[i].Optim(1);
+ if (verbosity>3)
cout << " NbSwap = " << NbSwap << endl;
NbTSwap += NbSwap ;
#ifdef DEBUG
-{
+{
Int4 NbErr=0;
Int4 i;
for (i=0;i<nbt;i++)
- if (triangles[i].link)
+ if (triangles[i].link)
{
double dd =Det(triangles[i][1].r-triangles[i][0].r,triangles[i][2].r-triangles[i][0].r);
- if(dd <=0)
+ if(dd <=0)
{
NbErr++;
cerr << " det triangle i " << i << " = " << triangles[i].det ;
cerr << " det triangle " << dd ;
cerr << " Les trois sommets " ;
- cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
+ cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
<< Number(triangles[i][2]) << endl;
cerr << "I2 " <<triangles[i][0].r << triangles[i][1].r << triangles[i][2].r << endl;
cerr << "R2 " << triangles[i][0].i << triangles[i][1].i << triangles[i][2].i << endl;
- cerr << "I2-R2 =" <<toR2(triangles[i][0].i)-triangles[i][0].r
+ cerr << "I2-R2 =" <<toR2(triangles[i][0].i)-triangles[i][0].r
<< toR2(triangles[i][1].i)-triangles[i][1].r
<< toR2(triangles[i][2].i)-triangles[i][2].r << endl;
}
}
- if(NbErr) {
+ if(NbErr) {
#ifdef DRAWING
Int4 kkk=0;
// UnMarkUnSwapTriangle();
@@ -1830,16 +1830,16 @@ Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
cout << " Nb of swap louche " << kkk << endl;
if(kkk) {
for (i=0;i<nbt;i++)
- if (triangles[i].link)
+ if (triangles[i].link)
{
double dd =Det(triangles[i][1].r-triangles[i][0].r,triangles[i][2].r-triangles[i][0].r);
- if(dd <=0)
+ if(dd <=0)
{
NbErr++;
cerr << " xxxdet triangle i " << i << " = " << triangles[i].det ;
cerr << " xxxdet triangle " << dd ;
cerr << " xxxLes trois sommets " ;
- cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
+ cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "
<< Number(triangles[i][2]) << endl;
cerr << triangles[i][0].r << triangles[i][1].r << triangles[i][2].r << endl;
cerr << triangles[i][0].i << triangles[i][1].i << triangles[i][2].i << endl;
@@ -1849,16 +1849,16 @@ Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap)
#endif
// MeshError(11);
}
-
+
}
#endif
return nbv-nbvold;
-
+
}
void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
{ // Triangles::NewPoints
Int4 nbtold(nbt),nbvold(nbv);
- if (verbosity>2)
+ if (verbosity>2)
cout << " -- Triangles::NewPoints ";
if (verbosity>3)cout << " nbv (in) on Boundary = " << nbv <<endl;
Int4 i,k;
@@ -1867,86 +1867,86 @@ void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
Int4 NbTSwap =0;
// insert old point
nbtold = nbt;
-
- if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx))
+
+ if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx))
{
// Bh.SetVertexFieldOn();
for (i=0;i<Bh.nbv;i++)
- {
+ {
Vertex & bv = Bh[i];
if (!bv.on) {
vertices[nbv].r = bv.r;
vertices[nbv++].m = bv.m;}
}
int nbv1=nbv;
- Bh.ReMakeTriangleContainingTheVertex();
- InsertNewPoints(nbvold,NbTSwap) ;
+ Bh.ReMakeTriangleContainingTheVertex();
+ InsertNewPoints(nbvold,NbTSwap) ;
if (verbosity>2)
- cout << " (Nb of Points from background mesh = "
+ cout << " (Nb of Points from background mesh = "
<< nbv-nbvold << " / " << nbv1-nbvold << ")" << endl;
- }
- else
- Bh.ReMakeTriangleContainingTheVertex();
+ }
+ else
+ Bh.ReMakeTriangleContainingTheVertex();
Triangle *t;
- // generation of the list of next Triangle
+ // generation of the list of next Triangle
// at 1 time we test all the triangles
Int4 Headt =0,next_t;
for(i=0;i<nbt;i++)
first_np_or_next_t[i]=-(i+1);
- // end list i >= nbt
- // the list of test triangle is
+ // end list i >= nbt
+ // the list of test triangle is
// the next traingle on i is -first_np_or_next_t[i]
int iter=0;
- // Big loop
+ // Big loop
do {
iter++;
nbtold = nbt;
nbvold = nbv;
-#ifdef DRAWING1
+#ifdef DRAWING1
inquire();
-#endif
+#endif
// default size of IntersectionTriangle
i=Headt;
next_t=-first_np_or_next_t[i];
- for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
+ for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
{ // for each triangle t
// we can change first_np_or_next_t[i]
// cout << " Do the triangle " << i << " Next_t=" << next_t << endl;
assert(i>=0 && i < nbt );
- first_np_or_next_t[i] = iter;
+ first_np_or_next_t[i] = iter;
for(j=0;j<3;j++)
- { // for each edge
+ { // for each edge
TriangleAdjacent tj(t,j);
Vertex & vA = * tj.EdgeVertex(0);
Vertex & vB = * tj.EdgeVertex(1);
-
+
if (!t->link) continue;// boundary
if (t->det <0) continue;
if (t->Locked(j)) continue;
- TriangleAdjacent tadjj = t->Adj(j);
+ TriangleAdjacent tadjj = t->Adj(j);
Triangle * ta= tadjj;
- if (ta->det <0) continue;
-
+ if (ta->det <0) continue;
+
R2 A = vA;
R2 B = vB;
-
+
k=Number(ta);
-
- if(first_np_or_next_t[k]==iter) // this edge is done before
- continue; // next edge of the triangle
-
+
+ if(first_np_or_next_t[k]==iter) // this edge is done before
+ continue; // next edge of the triangle
+
//const Int4 NbvOld = nbv;
lIntTria.SplitEdge(Bh,A,B);
lIntTria.NewPoints(vertices,nbv,nbvx);
- } // end loop for each edge
-
- }// for triangle
-
+ } // end loop for each edge
+
+ }// for triangle
+
#ifdef DRAWING1
cout << " -------------------------------------------- " << endl;
inquire();
@@ -1955,37 +1955,37 @@ void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
penthickness(2);
#endif
- if (!InsertNewPoints(nbvold,NbTSwap))
+ if (!InsertNewPoints(nbvold,NbTSwap))
break;
-
+
for (i=nbtold;i<nbt;i++)
first_np_or_next_t[i]=iter;
- Headt = nbt; // empty list
- for (i=nbvold;i<nbv;i++)
+ Headt = nbt; // empty list
+ for (i=nbvold;i<nbv;i++)
{ // for all the triangle contening the vertex i
Vertex * s = vertices + i;
TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
Triangle * tbegin= (Triangle*) ta;
Int4 kt;
- do {
+ do {
kt = Number((Triangle*) ta);
- if (first_np_or_next_t[kt]>0)
+ if (first_np_or_next_t[kt]>0)
first_np_or_next_t[kt]=-Headt,Headt=kt;
assert( ta.EdgeVertex(0) == s);
ta = Next(Adj(ta));
- } while ( (tbegin != (Triangle*) ta));
- }
-
+ } while ( (tbegin != (Triangle*) ta));
+ }
+
} while (nbv!=nbvold);
-
+
delete [] first_np_or_next_t;
Int4 NbSwapf =0,NbSwp;
-
- // bofbof
-
-
+
+ // bofbof
+
+
NbSwp = NbSwapf;
for (i=0;i<nbv;i++)
NbSwapf += vertices[i].Optim(0);
@@ -2001,8 +2001,8 @@ void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
*/
NbTSwap += NbSwapf ;
if (verbosity>3) cout << " " ;
- if (verbosity>2)
- cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT
+ if (verbosity>2)
+ cout << " Nb Of Vertices =" << nbv << " Nb of triangles = " << nbt-NbOutT
<< " NbSwap final = " << NbSwapf << " Nb Total Of Swap = " << NbTSwap << endl;
@@ -2010,8 +2010,8 @@ void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex)
void Triangles::NewPointsOld(Triangles & Bh)
{ // Triangles::NewPointsOld
- Real8 seuil= 0.7 ;// for two neart point
- if (verbosity>1)
+ Real8 seuil= 0.7 ;// for two neart point
+ if (verbosity>1)
cout << " begin : Triangles::NewPointsOld " << endl;
Int4 i,k;
int j;
@@ -2025,29 +2025,29 @@ void Triangles::NewPointsOld(Triangles & Bh)
Int4 ColorEdge[3];
Int4 color=-1;
Triangle *t;
- // generation of the list of next Triangle
+ // generation of the list of next Triangle
// at 1 time we test all the triangles
Int4 Headt =0,next_t;
for(i=0;i<nbt;i++)
first_np_or_next_t[i]=-(i+1);
- // end list i >= nbt
- // the list of test triangle is
+ // end list i >= nbt
+ // the list of test triangle is
// the next Triangle on i is -first_np_or_next_t[i]
Int4 nbtold,nbvold;
- // Big loop
+ // Big loop
do {
nbtold = nbt;
nbvold = nbv;
-#ifdef DRAWING1
+#ifdef DRAWING1
inquire();
-#endif
+#endif
// default size of IntersectionTriangle
i=Headt;
next_t=-first_np_or_next_t[i];
- for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
+ for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i])
{ // for each triangle t
// we can change first_np_or_next_t[i]
#ifdef TRACETRIANGLE
@@ -2057,20 +2057,20 @@ void Triangles::NewPointsOld(Triangles & Bh)
assert(i>=0 && i < nbt );
first_np_or_next_t[i] = nbv; // to save the fist new point of triangle
for(j=0;j<3;j++)
- { // for each edge
+ { // for each edge
TriangleAdjacent tj(t,j);
// color++;// the color is 3i+j
color = 3*i + j ;;
ColorEdge[j]=color;
BeginNewPoint[j]=nbv;
EndNewPoint[j]=nbv-1;
- step[j]=1;// right sens
+ step[j]=1;// right sens
Vertex & vA = * tj.EdgeVertex(0);
Vertex & vB = * tj.EdgeVertex(1);
-
+
#ifdef TRACETRIANGLE
if(trace) {
- cout << " i " << Number(vA) <<" j "<< Number(vB)
+ cout << " i " << Number(vA) <<" j "<< Number(vB)
<< " " << t->Locked(j) ;
}
#endif
@@ -2078,32 +2078,32 @@ void Triangles::NewPointsOld(Triangles & Bh)
if (t->det <0) continue;
if (t->Locked(j)) continue;
- TriangleAdjacent tadjj = t->Adj(j);
+ TriangleAdjacent tadjj = t->Adj(j);
Triangle * ta= tadjj;
- if (ta->det <0) continue;
-
+ if (ta->det <0) continue;
+
R2 A = vA;
R2 B = vB;
-
+
k=Number(ta);
- // the 2 opposite vertices
+ // the 2 opposite vertices
const Vertex & vC1 = *tj.OppositeVertex();
const Vertex & vC2 = *tadjj.OppositeVertex();
-
+
#ifdef TRACETRIANGLE
trace = trace || k == TRACETRIANGLE;
if(trace) {
- cout << "Test Arete " << i << " AB = " << A << B
- << "i " <<Number(vA)<< "j" <<Number(vB);
- cout << " link" <<(int)t->link << " ta=" << Number( ta)
+ cout << "Test Arete " << i << " AB = " << A << B
+ << "i " <<Number(vA)<< "j" <<Number(vB);
+ cout << " link" <<(int)t->link << " ta=" << Number( ta)
<< " det " <<ta->det ;
cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ;
cout << " loc " << ta->Locked(j) << endl;
}
#endif
-
- if(first_np_or_next_t[k]>0) { // this edge is done before
+
+ if(first_np_or_next_t[k]>0) { // this edge is done before
// find the color of the edge and begin , end of newpoint
register int kk = t->NuEdgeTriangleAdj(j);
assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1]));
@@ -2111,20 +2111,20 @@ void Triangles::NewPointsOld(Triangles & Bh)
register Int4 kolor =3*k + kk;
ColorEdge[j]=kolor;
register Int4 kkk= 1;
- step[j]=-1;// other sens
+ step[j]=-1;// other sens
BeginNewPoint[j]=0;
- EndNewPoint[j]=-1; // empty list
- for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++)
- if (vertices[iv].color > kolor)
- break; // the color is passed
+ EndNewPoint[j]=-1; // empty list
+ for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++)
+ if (vertices[iv].color > kolor)
+ break; // the color is passed
else if (vertices[iv].color == kolor) {
- EndNewPoint[j]=iv;
+ EndNewPoint[j]=iv;
if (kkk) // one time test
kkk=0,BeginNewPoint[j]=iv;}
- continue; // next edge of the triangle
- } // end if( k < i)
-
-
+ continue; // next edge of the triangle
+ } // end if( k < i)
+
+
#ifdef DRAWING1
penthickness(2); Move(A);Line(B); penthickness(1);
#endif
@@ -2146,51 +2146,51 @@ void Triangles::NewPointsOld(Triangles & Bh)
&& ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) )
nbv++;
}
-
+
EndNewPoint[j] = nbv-1;
- } // end loop for each edge
+ } // end loop for each edge
#ifdef TRACETRIANGLE
if(trace) {
- // verification des point cree
- cout << "\n ------------ " << t->link << " " << t->det
+ // verification des point cree
+ cout << "\n ------------ " << t->link << " " << t->det
<< " b " << BeginNewPoint[0] << " " << BeginNewPoint[1]
- << " " << BeginNewPoint[2] << " "
- << " e " << EndNewPoint[0] << " " << EndNewPoint[1]
- << " " << EndNewPoint[2] << " "
- << " s " << step[0] << " " << step[1] << " " << step[2] << " "
+ << " " << BeginNewPoint[2] << " "
+ << " e " << EndNewPoint[0] << " " << EndNewPoint[1]
+ << " " << EndNewPoint[2] << " "
+ << " s " << step[0] << " " << step[1] << " " << step[2] << " "
<< endl;
}
#endif
if (!t->link) continue;// boundary
- if (t->det<=0) continue;// outside
+ if (t->det<=0) continue;// outside
// continue;
for(int j0=0;j0<3;j0++)
for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++)
{
- // find the neart point one the opposite edge
+ // find the neart point one the opposite edge
// to compute i1
Vertex & vi0 = vertices[i0];
int kstack = 0;
Int4 stack[10];
// Int4 savRef[10];
- int j1 = j0;
+ int j1 = j0;
while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge
// computation of the intersection of edge j1 and DOrto
- // take the good sens
-
+ // take the good sens
+
if (BeginNewPoint[j1]> EndNewPoint[j1])
- continue; //
+ continue; //
else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) {
for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++)
stack[kstack++] = ii1;
continue;}
-
-
+
+
int k0,k1;
if (step[j1]<0) k0=1,k1=0; // reverse
- else k0=0,k1=1;
+ else k0=0,k1=1;
R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]];
R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]];
R2 D = V11-V10;
@@ -2202,9 +2202,9 @@ void Triangles::NewPointsOld(Triangles & Bh)
Real8 s;
//cout << " --i0 = " << i0 << D << V10 << V11 << endl ;
//cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl;
- if (( c10 < c0 ) && (c0 < c11))
+ if (( c10 < c0 ) && (c0 < c11))
s = (c11-c0)/(c11-c10);
- else if (( c11 < c0 ) && (c0 < c10))
+ else if (( c11 < c0 ) && (c0 < c10))
s = (c11-c0) /(c11-c10);
else break;
R2 VP = V10*s + V11*(1-s);
@@ -2218,11 +2218,11 @@ void Triangles::NewPointsOld(Triangles & Bh)
// find the 2 point by dichotomie
//cout << " t =" << Number(t) << " c0 " << c0 ;
Int4 ii0 = BeginNewPoint[j1];
- Int4 ii1 = EndNewPoint[j1];
+ Int4 ii1 = EndNewPoint[j1];
Real8 ciii=-1,cii0=-1,cii1=-1 ;
- if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )
+ if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )
stack[kstack++] = ii0;//cout << " add+0 " << ii0;
- else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )
+ else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )
stack[kstack++] = ii1;//cout << " add+1 " << ii1;
else {
while ((ii1-ii0)> 1) {
@@ -2230,31 +2230,31 @@ void Triangles::NewPointsOld(Triangles & Bh)
ciii = vi0.m( D ,(R2) vertices[iii]);
//cout << " (iii = " << iii << " " << ciii << ") ";
if ( sss * (ciii - c0) <0 ) ii0 = iii;
- else ii1 = iii;}
+ else ii1 = iii;}
stack[kstack++] = ii0;// cout << " add0 " << ii0;
if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1;
}
#ifdef DEBUG2
- cout << "ii1 = " << ii1
+ cout << "ii1 = " << ii1
<< " ii0 = " << ii0 << endl;
- cout << " cccc = " << cii0 << " " << ciii
+ cout << " cccc = " << cii0 << " " << ciii
<< " " << cii1 << " sss=" << sss << endl;
#endif
if (kstack >5) // bug ?
cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl;
}
-
+
stack[kstack++] = -1; // to stop
Int4 i1;
- kstack =0;
- while( (i1=stack[kstack++]) >= 0)
- { // the two parameter is i0 and i1
+ kstack =0;
+ while( (i1=stack[kstack++]) >= 0)
+ { // the two parameter is i0 and i1
assert(i1 < nbv && i1 >= 0);
assert(i0 < nbv && i0 >= 0);
assert(i1 != i0);
R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];
Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));
-
+
#ifdef DRAWING1
Move(vertices[i0].r);
@@ -2262,15 +2262,15 @@ void Triangles::NewPointsOld(Triangles & Bh)
#endif
#ifdef TRACETRIANGLE
if(trace) {
- cout << "\n test j" << j <<" " << i0
+ cout << "\n test j" << j <<" " << i0
<< " " << i1 << " d01=" << d01 <<endl;}
#endif
assert (i0 >= nbvold);
assert (i1 >= nbvold);
assert(i0 != i1);
- if (d01 == 0)
- break;
- if ( d01 < seuil)
+ if (d01 == 0)
+ break;
+ if ( d01 < seuil)
if (i1<nbvold) {
// remove all the points i0;
register Int4 ip,ipp;
@@ -2285,13 +2285,13 @@ void Triangles::NewPointsOld(Triangles & Bh)
// count the number of common points to compute weight w0,w1
for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++;
for (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++;
-
+
register Real8 w0 = ((Real8) kk0)/(kk0+kk1);
register Real8 w1 = ((Real8) kk1)/(kk0+kk1);
// make a circular link
Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber);
- // the new coordinate
+ // the new coordinate
R2 C //= vertices[i0] ;
= vertices[i0].r *w0 + vertices[i1].r* w1;
@@ -2299,20 +2299,20 @@ void Triangles::NewPointsOld(Triangles & Bh)
if(trace) {
cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl;
}
-#endif
-#ifdef DRAWING1
+#endif
+#ifdef DRAWING1
Move(vertices[i0].r);
Line(vertices[i1].r);
DrawMark(C);
-#endif
- // update the new point points of the list
+#endif
+ // update the new point points of the list
for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)
- vertices[ip0].r = C;
+ vertices[ip0].r = C;
vertices[ip0].r = C;
}
}
} // for (i0= ....
- }// for triangle
+ }// for triangle
#ifdef DRAWING1
cout << " -------------------------------------------- " << endl;
@@ -2321,91 +2321,91 @@ void Triangles::NewPointsOld(Triangles & Bh)
Draw();
penthickness(2);
#endif
-
+
// remove of all the double points
Int4 ip,ipp,kkk=nbvold;
- for (i=nbvold;i<nbv;i++)
+ for (i=nbvold;i<nbv;i++)
if (vertices[i].ReferenceNumber>=0) {// good points
// cout <<" i = " << i ;
for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp)
vertices[ip].ReferenceNumber = -1;// mark remove
vertices[ip].ReferenceNumber = -1;// mark remove
- // cout << i << " ---> " << kkk << endl;
+ // cout << i << " ---> " << kkk << endl;
vertices[kkk] = vertices[i];
vertices[kkk].i = toI2(vertices[kkk].r);
#ifdef DRAWING1
DrawMark(vertices[kkk]);
#endif
vertices[kkk++].ReferenceNumber = 0;
-
- }
-#ifdef DRAWING1
+
+ }
+#ifdef DRAWING1
penthickness(1);
#endif
-
- // insertion part ---
+
+ // insertion part ---
const Int4 nbvnew = kkk-nbvold;
cout << " Remove " << nbv - kkk << " to close vertex " ;
- cout << " and Insert the " <<nbvnew<< " new points " << endl;
+ cout << " and Insert the " <<nbvnew<< " new points " << endl;
nbv = kkk;
Int4 NbSwap=0;
Icoor2 dete[3];
-
- // construction d'un ordre aleatoire
- if (! nbvnew)
- break;
+
+ // construction d'un ordre aleatoire
+ if (! nbvnew)
+ break;
if (nbvnew) {
const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
- Int4 k3 = rand()%nbvnew ;
- for (Int4 is3=0; is3<nbvnew; is3++)
+ Int4 k3 = rand()%nbvnew ;
+ for (Int4 is3=0; is3<nbvnew; is3++)
ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)];
- for (i=nbvold;i<nbv;i++)
+ for (i=nbvold;i<nbv;i++)
{ Vertex * vi = ordre[i];
Triangle *tcvi = FindTriangleContening(vi->i,dete);
// Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y);
- // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is "
+ // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is "
// << Number(nv) << nv->i << endl;
// Int4 kt = Number(tcvi);
- //
-
+ //
+
quadtree->Add(*vi); //
#ifdef DRAWING1
DrawMark(vi->r);
#endif
- assert (tcvi->det >= 0) ;// internal
- Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);
- }
+ assert (tcvi->det >= 0) ;// internal
+ Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);
+ }
}
cout << " Nb swap = " << NbSwap << " after " ;
#ifdef DRAWING1
inquire();
#endif
- for (i=nbvold;i<nbv;i++)
- NbSwap += vertices[i].Optim(1);
+ for (i=nbvold;i<nbv;i++)
+ NbSwap += vertices[i].Optim(1);
cout << NbSwap << endl;
for (i=nbtold;i<nbt;i++)
first_np_or_next_t[i]=1;
- Headt = nbt; // empty list
- for (i=nbvold;i<nbv;i++)
+ Headt = nbt; // empty list
+ for (i=nbvold;i<nbv;i++)
{ // for all the triangle contening the vertex i
Vertex * s = vertices + i;
TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);
Triangle * tbegin= (Triangle*) ta;
Int4 kt;
- do {
+ do {
kt = Number((Triangle*) ta);
- if (first_np_or_next_t[kt]>0)
+ if (first_np_or_next_t[kt]>0)
first_np_or_next_t[kt]=-Headt,Headt=kt;
assert( ta.EdgeVertex(0) == s);
ta = Next(Adj(ta));
- } while ( (tbegin != (Triangle*) ta));
+ } while ( (tbegin != (Triangle*) ta));
}
@@ -2414,7 +2414,7 @@ void Triangles::NewPointsOld(Triangles & Bh)
#ifdef DEBUG
int nberr=0;
for (int it=0;it<nbt;it++)
- if(triangles[it].det==0)
+ if(triangles[it].det==0)
if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;
if(nberr) MeshError(992,this);
#endif
@@ -2422,7 +2422,7 @@ void Triangles::NewPointsOld(Triangles & Bh)
}
-void Triangles::Insert()
+void Triangles::Insert()
{
if (verbosity>2) cout << " -- Insert initial " << nbv << " vertices " << endl ;
Triangles * OldCurrentTh =CurrentTh;
@@ -2431,42 +2431,42 @@ void Triangles::Insert()
double time0=CPUtime(),time1,time2,time3;
SetIntCoor();
Int4 i;
- for (i=0;i<nbv;i++)
+ for (i=0;i<nbv;i++)
ordre[i]= &vertices[i] ;
- // construction d'un ordre aleatoire
+ // construction d'un ordre aleatoire
const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;
- Int4 k3 = rand()%nbv ;
- for (int is3=0; is3<nbv; is3++)
+ Int4 k3 = rand()%nbv ;
+ for (int is3=0; is3<nbv; is3++)
ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
- for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
+ for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
if ( ++i >= nbv) {
cerr << " All the vertices are aline " << endl;
MeshError(998,this); }
-
- // echange i et 2 dans ordre afin
+
+ // echange i et 2 dans ordre afin
// que les 3 premiers ne soit pas aligne
Exchange( ordre[2], ordre[i]);
-
+
// on ajoute un point a l'infini pour construire le maillage
// afin d'avoir une definition simple des aretes frontieres
nbt = 2;
-
-
+
+
// on construit un maillage trivale forme
// d'une arete et de 2 triangles
- // construit avec le 2 aretes orientes et
+ // construit avec le 2 aretes orientes et
Vertex * v0=ordre[0], *v1=ordre[1];
-
- triangles[0](0) = 0; // sommet pour infini
+
+ triangles[0](0) = 0; // sommet pour infini
triangles[0](1) = v0;
triangles[0](2) = v1;
-
- triangles[1](0) = 0;// sommet pour infini
+
+ triangles[1](0) = 0;// sommet pour infini
triangles[1](2) = v0;
triangles[1](1) = v1;
const int e0 = OppositeEdge[0];
@@ -2475,26 +2475,26 @@ void Triangles::Insert()
triangles[0].SetAdj2(e0, &triangles[1] ,e0);
triangles[0].SetAdj2(e1, &triangles[1] ,e2);
triangles[0].SetAdj2(e2, &triangles[1] ,e1);
-
+
triangles[0].det = -1; // faux triangles
triangles[1].det = -1; // faux triangles
-
+
triangles[0].SetTriangleContainingTheVertex();
triangles[1].SetTriangleContainingTheVertex();
-
+
triangles[0].link=&triangles[1];
triangles[1].link=&triangles[0];
-
-#ifdef DEBUG
+
+#ifdef DEBUG
triangles[0].check();
triangles[1].check();
-#endif
+#endif
// nbtf = 2;
if ( !quadtree ) quadtree = new QuadTree(this,0);
quadtree->Add(*v0);
quadtree->Add(*v1);
-
- // on ajoute les sommets un � un
+
+ // on ajoute les sommets un � un
Int4 NbSwap=0;
time1=CPUtime();
@@ -2506,67 +2506,67 @@ void Triangles::Insert()
// cout << " Insert " << Number(vi) << endl;
Icoor2 dete[3];
Triangle *tcvi = FindTriangleContening(vi->i,dete);
- quadtree->Add(*vi);
+ quadtree->Add(*vi);
Add(*vi,tcvi,dete);
NbSwap += vi->Optim(1,0);
#ifdef DRAWING1
inquire();
#endif
-
+
}// fin de boucle en icount
time2=CPUtime();
- if (verbosity>3)
- cout << " NbSwap of insertion " << NbSwap
- << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
- << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv "
- << (float)NbUnSwap /(float) nbv
+ if (verbosity>3)
+ cout << " NbSwap of insertion " << NbSwap
+ << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
+ << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv "
+ << (float)NbUnSwap /(float) nbv
<<endl;
NbUnSwap = 0;
- // construction d'un ordre aleatoire
+ // construction d'un ordre aleatoire
// const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ;
#ifdef NBLOOPOPTIM
- k3 = rand()%nbv ;
- for (int is4=0; is4<nbv; is4++)
+ k3 = rand()%nbv ;
+ for (int is4=0; is4<nbv; is4++)
ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv];
-
+
double timeloop = time2 ;
- for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++)
+ for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++)
{
double time000 = timeloop;
Int4 NbSwap = 0;
- for (int is1=0; is1<nbv; is1++)
+ for (int is1=0; is1<nbv; is1++)
NbSwap += ordre[is1]->Optim(0,0);
timeloop = CPUtime();
- if (verbosity>3)
- cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap
- << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
- << " CPU=" << timeloop - time000 << " s, "
- << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv
+ if (verbosity>3)
+ cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap
+ << " NbSwap/Nbv " << (float) NbSwap / (float) nbv
+ << " CPU=" << timeloop - time000 << " s, "
+ << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv
<< endl;
NbUnSwap = 0;
if(!NbSwap) break;
}
- ReMakeTriangleContainingTheVertex();
+ ReMakeTriangleContainingTheVertex();
// because we break the TriangleContainingTheVertex
#endif
#ifdef DEBUG
int nberr=0;
for (int it=0;it<nbt;it++)
- if(triangles[it].det==0)
+ if(triangles[it].det==0)
if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;
if(nberr) MeshError(991,this);
#endif
time3=CPUtime();
- if (verbosity>4)
- cout << " init " << time1 - time0 << " initialisation, "
- << time2 - time1 << "s, insert point "
+ if (verbosity>4)
+ cout << " init " << time1 - time0 << " initialisation, "
+ << time2 - time1 << "s, insert point "
<< time3 -time2 << "s, optim " << endl
<< " Init Total Cpu Time = " << time3 - time0 << "s " << endl;
-
+
#ifdef DRAWING1
inquire();
-#endif
+#endif
CurrentTh=OldCurrentTh;
}
@@ -2577,18 +2577,18 @@ void Triangles::ForceBoundary()
cout << " -- ForceBoundary nb of edge " << nbe << endl;
int k=0;
Int4 nbfe=0,nbswp=0,Nbswap=0;
- for (Int4 t = 0; t < nbt; t++)
+ for (Int4 t = 0; t < nbt; t++)
if (!triangles[t].det)
k++,cerr << " det T" << t << " = " << 0 << endl;
if (k!=0) {
cerr << " ther is " << k << " triangles of mes = 0 " << endl;
MeshError(11,this);}
-
+
TriangleAdjacent ta(0,0);
- for (Int4 i = 0; i < nbe; i++)
+ for (Int4 i = 0; i < nbe; i++)
{
nbswp = ForceEdge(edges[i][0],edges[i][1],ta);
-
+
if ( nbswp < 0) k++;
else Nbswap += nbswp;
if (nbswp) nbfe++;
@@ -2598,17 +2598,17 @@ void Triangles::ForceBoundary()
<<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ;
if(ta.t)
{
- Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1);
+ Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1);
cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n";
}
else cerr << endl;
-
+
}
if ( nbswp >=0 && edges[i].on->Cracked())
ta.SetCracked();
}
-
-
+
+
if (k!=0) {
cerr << " they is " << k << " lost edges " << endl;
cerr << " The boundary is crossing may be!" << endl;
@@ -2624,106 +2624,106 @@ void Triangles::ForceBoundary()
void Triangles::FindSubDomain(int OutSide=0)
{
//#define DRAWING1
-
+
if (verbosity >2)
{
if (OutSide)
- cout << " -- Find all external sub-domain ";
+ cout << " -- Find all external sub-domain ";
else
cout << " -- Find all internal sub-domain ";
if(verbosity>99)
{
-
+
for(int i=0;i<nbt;++i)
cout << i<< " " << triangles[i] << endl;
}
-
+
}
// if (verbosity > 4) cout << " OutSide=" << OutSide << endl;
short * HeapArete = new short[nbt];
Triangle ** HeapTriangle = new Triangle* [nbt];
Triangle *t,*t1;
Int4 k,it;
-
- for (Int4 itt=0;itt<nbt;itt++)
+
+ for (Int4 itt=0;itt<nbt;itt++)
triangles[itt].link=0; // par defaut pas de couleur
#ifdef DRAWING1
reffecran();
#endif
-
+
Int4 NbSubDomTot =0;
- for ( it=0;it<nbt;it++) {
+ for ( it=0;it<nbt;it++) {
if ( ! triangles[it].link ) {
t = triangles + it;
NbSubDomTot++;; // new composante connexe
- Int4 i = 0; // niveau de la pile
+ Int4 i = 0; // niveau de la pile
t->link = t ; // sd forme d'un triangle cicular link
#ifdef DRAWING1
t->Draw(NbSubDomTot-1);
#endif
-
-
- HeapTriangle[i] =t ;
+
+
+ HeapTriangle[i] =t ;
HeapArete[i] = 3;
-
+
while (i >= 0) // boucle sur la pile
- { while ( HeapArete[i]--) // boucle sur les 3 aretes
- {
+ { while ( HeapArete[i]--) // boucle sur les 3 aretes
+ {
int na = HeapArete[i];
Triangle * tc = HeapTriangle[i]; // triangle courant
if( ! tc->Locked(na)) // arete non frontiere
{
Triangle * ta = tc->TriangleAdj(na) ; // n� triangle adjacent
if (ta->link == 0 ) // non deja chainer => on enpile
- {
+ {
i++;
#ifdef DRAWING1
ta->Draw(NbSubDomTot-1);
#endif
ta->link = t->link ; // on chaine les triangles
- t->link = ta ; // d'un meme sous domaine
+ t->link = ta ; // d'un meme sous domaine
HeapArete[i] = 3; // pour les 3 triangles adjacents
HeapTriangle[i] = ta;
}}
} // deplie fin de boucle sur les 3 adjacences
i--;
- }
- }
+ }
+ }
}
-
+
// supression de tous les sous domaine infini <=> contient le sommet NULL
it =0;
NbOutT = 0;
while (it<nbt) {
- if (triangles[it].link)
- {
- if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))
+ if (triangles[it].link)
+ {
+ if (!( triangles[it](0) && triangles[it](1) && triangles[it](2) ))
{
- // infini triangle
+ // infini triangle
NbSubDomTot --;
// cout << " triangle infini " << it << triangles[it] << endl;
t=&triangles[it];
- NbOutT--; // on fait un coup de trop.
+ NbOutT--; // on fait un coup de trop.
while (t){ // cout << Number(t) << " " << endl;
NbOutT++;
t1=t;
t=t->link;
t1->link=0;}//while (t)
}
- }
+ }
it++;} // end while (it<nbt)
- if (nbt == NbOutT || !NbSubDomTot)
+ if (nbt == NbOutT || !NbSubDomTot)
{
cout << "\n error : " << NbOutT << " " << NbSubDomTot <<" " << nbt << endl;
cerr << "Error: The boundary is not close => All triangles are outside " << endl;
MeshError(888,this);
}
-
+
delete [] HeapArete;
delete [] HeapTriangle;
-
-
- if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains )
+
+
+ if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains )
{ // No geom sub domain
Int4 i;
if (subdomains) delete [] subdomains;
@@ -2736,7 +2736,7 @@ void Triangles::FindSubDomain(int OutSide=0)
Int4 * mark = new Int4[nbt];
for (it=0;it<nbt;it++)
mark[it]=triangles[it].link ? -1 : -2;
-
+
it =0;
k = 0;
while (it<nbt) {
@@ -2744,26 +2744,26 @@ void Triangles::FindSubDomain(int OutSide=0)
t1 = & triangles[it];
t = t1->link;
mark[it]=k;
-#ifdef DRAWING1
+#ifdef DRAWING1
t1->Draw(k);
#endif
subdomains[k].head = t1;
// cout << " New -- " << Number(t1) << " " << it << endl;
do {// cout << " k " << k << " " << Number(t) << endl;
mark[Number(t)]=k;
-#ifdef DRAWING1
+#ifdef DRAWING1
t->Draw(k);
#endif
t=t->link;
} while (t!=t1);
-#ifdef DRAWING1
+#ifdef DRAWING1
t1->Draw(k);
#endif
mark[it]=k++;}
// else if(mark[it] == -2 ) triangles[it].Draw(999);
it++;} // end white (it<nbt)
assert(k== NbSubDomains);
- if(OutSide)
+ if(OutSide)
{
// to remove all the sub domain by parity adjacents
// because in this case we have only the true boundary edge
@@ -2787,18 +2787,18 @@ void Triangles::FindSubDomain(int OutSide=0)
nbk--,subdomains[kl].ref=-1;
if(kl<0 && kr >=0 && subdomains[kr].ref>=0)
nbk--,subdomains[kr].ref=-1;
- // cout << " after \t "
- // << kl << subdomains[kl].ref << " rr " << kr
+ // cout << " after \t "
+ // << kl << subdomains[kl].ref << " rr " << kr
// << subdomains[kr].ref << endl;
}
}
// cout << subdomains[0].ref << subdomains[1].ref << endl;
Int4 j=0;
for ( i=0;i<NbSubDomains;i++)
- if((-subdomains[i].ref) %2) { // good
+ if((-subdomains[i].ref) %2) { // good
//cout << " sudom ok = " << i << " " << subdomains[i].ref
// << " " << (-subdomains[i].ref) %2 << endl;
- if(i != j)
+ if(i != j)
Exchange(subdomains[i],subdomains[j]);
j++;}
else
@@ -2810,20 +2810,20 @@ void Triangles::FindSubDomain(int OutSide=0)
t=t->link;
t1->link=0;}//while (t)
}
-
+
if(verbosity>4)
cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl;
NbSubDomains=j;
}
-
- delete [] mark;
-
+
+ delete [] mark;
+
}
else
{ // find the head for all sub domaine
if (Gh.NbSubDomains != NbSubDomains && subdomains)
delete [] subdomains, subdomains=0;
- if (! subdomains )
+ if (! subdomains )
subdomains = new SubDomain[ Gh.NbSubDomains];
NbSubDomains =Gh.NbSubDomains;
if(verbosity>4)
@@ -2832,24 +2832,24 @@ void Triangles::FindSubDomain(int OutSide=0)
ReMakeTriangleContainingTheVertex();
Int4 * mark = new Int4[nbt];
Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge();
-
+
for (it=0;it<nbt;it++)
mark[it]=triangles[it].link ? -1 : -2;
Int4 inew =0;
- for (Int4 i=0;i<NbSubDomains;i++)
+ for (Int4 i=0;i<NbSubDomains;i++)
{
GeometricalEdge &eg = *Gh.subdomains[i].edge;
subdomains[i].ref = Gh.subdomains[i].ref;
int ssdlab = subdomains[i].ref;
- // by carefull is not easy to find a edge create from a GeometricalEdge
+ // by carefull is not easy to find a edge create from a GeometricalEdge
// see routine MakeGeometricalEdgeToEdge
Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)];
assert(&e);
Vertex * v0 = e(0),*v1 = e(1);
Triangle *t = v0->t;
int sens = Gh.subdomains[i].sens;
- // test if ge and e is in the same sens
- // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref
+ // test if ge and e is in the same sens
+ // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref
// << " ref edge =" << eg.ref << " sens " << sens ;
if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0)
sens = -sens ;
@@ -2858,40 +2858,40 @@ void Triangles::FindSubDomain(int OutSide=0)
// cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl;
// cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl;
assert(t && sens);
-
+
TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges
-
- while (1)
+
+ while (1)
{
assert( v0 == ta.EdgeVertex(1) );
// cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;
if (ta.EdgeVertex(0) == v1) { // ok we find the edge
- if (sens>0)
+ if (sens>0)
subdomains[i].head=t=Adj(ta);
- else
+ else
subdomains[i].head=t=ta;
//cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl;
- if(t<triangles || t >= triangles+nbt || t->det < 0
- || t->link == 0) // Ajoute aout 200
+ if(t<triangles || t >= triangles+nbt || t->det < 0
+ || t->link == 0) // Ajoute aout 200
{
cerr << " Error in the def of sub domain "<<i
<< " form border " << NbSubDomains - i << "/" << NbSubDomains
- << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl;
+ << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl;
err++;
break;}
Int4 it = Number(t);
if (mark[it] >=0) {
if(verbosity>10)
- cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref
+ cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref
<< " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref
<< " skip this def " << endl;
break;}
- if(i != inew)
+ if(i != inew)
Exchange(subdomains[i],subdomains[inew]);
inew++;
Triangle *tt=t;
Int4 kkk=0;
- do
+ do
{
kkk++;
assert(mark[Number(tt)]<0);
@@ -2904,38 +2904,38 @@ void Triangles::FindSubDomain(int OutSide=0)
if(verbosity>7)
cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl;
break;}
- ta = Previous(Adj(ta));
+ ta = Previous(Adj(ta));
if(t == (Triangle *) ta) {
err++;
- cerr << " Error in the def of sub domain " << i
+ cerr << " Error in the def of sub domain " << i
<< " edge=" << Gh.Number(eg) << " " << sens << endl;
break;}
// cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl;
-
+
}
-
+
}
if (err) MeshError(777,this);
-
+
if (inew < NbSubDomains) {
- if (verbosity>5)
+ if (verbosity>5)
cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl;
NbSubDomains=inew;}
-
-
+
+
for (it=0;it<nbt;it++)
- if ( mark[it] ==-1 )
+ if ( mark[it] ==-1 )
NbOutT++,triangles[it].link =0;
delete [] GeometricalEdgetoEdge;
delete [] mark;
-
+
}
-
+
#ifdef DRAWING1
inquire();
#endif
NbOutT=0;
- for (it=0;it<nbt;it++)
+ for (it=0;it<nbt;it++)
if(!triangles[it].link) NbOutT++;
if (verbosity> 4)
cout << " " ;
@@ -2944,32 +2944,32 @@ void Triangles::FindSubDomain(int OutSide=0)
#ifdef DRAWING1
inquire();
#endif
-
+
//#undef DRAWING1
-
-
+
+
}
void Triangles::ReNumberingVertex(Int4 * renu)
{
// warning be carfull because pointeur
- // from on mesh to over mesh
+ // from on mesh to over mesh
// -- so do ReNumbering a the beginning
Vertex * ve = vertices+nbv;
Int4 it,ie,i;
-
- for ( it=0;it<nbt;it++)
+
+ for ( it=0;it<nbt;it++)
triangles[it].ReNumbering(vertices,ve,renu);
-
- for ( ie=0;ie<nbe;ie++)
+
+ for ( ie=0;ie<nbe;ie++)
edges[ie].ReNumbering(vertices,ve,renu);
-
+
for (i=0;i< NbVerticesOnGeomVertex;i++)
{
Vertex *v = VerticesOnGeomVertex[i].mv;
if (v>=vertices && v < ve)
VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)];
}
-
+
for (i=0;i< NbVerticesOnGeomEdge;i++)
{
Vertex *v =VerticesOnGeomEdge[i].mv;
@@ -2983,7 +2983,7 @@ void Triangles::ReNumberingVertex(Int4 * renu)
if (v>=vertices && v < ve)
VertexOnBThVertex[i].v=vertices+renu[Number(v)];
}
-
+
for (i=0;i< NbVertexOnBThEdge;i++)
{
Vertex *v=VertexOnBThEdge[i].v;
@@ -2991,44 +2991,44 @@ void Triangles::ReNumberingVertex(Int4 * renu)
VertexOnBThEdge[i].v=vertices+renu[Number(v)];
}
- // move the Vertices without a copy of the array
- // be carefull not trivial code
+ // move the Vertices without a copy of the array
+ // be carefull not trivial code
Int4 j;
for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
if (renu[it] >= 0) // a new sub cycle
- {
+ {
i=it;
Vertex ti=vertices[i],tj;
- while ( (j=renu[i]) >= 0)
- { // i is old, and j is new
- renu[i] = -1-renu[i]; // mark
+ while ( (j=renu[i]) >= 0)
+ { // i is old, and j is new
+ renu[i] = -1-renu[i]; // mark
tj = vertices[j]; // save new
vertices[j]= ti; // new <- old
- i=j; // next
+ i=j; // next
ti = tj;
- }
+ }
}
- if (quadtree)
+ if (quadtree)
{ delete quadtree;
quadtree = new QuadTree(this);
}
for ( it=0;it<nbv;it++)
renu[i]= -renu[i]-1;
-
+
}
void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress)
{
Int4 *renu= new Int4[nbt];
register Triangle *t0,*t,*te=triangles+nbt;
register Int4 k=0,it,i,j;
-
- for ( it=0;it<nbt;it++)
- renu[it]=-1; // outside triangle
+
+ for ( it=0;it<nbt;it++)
+ renu[it]=-1; // outside triangle
for ( i=0;i<NbSubDomains;i++)
- {
+ {
t=t0=subdomains[i].head;
assert(t0); // no empty sub domain
- do {
+ do {
Int4 kt = Number(t);
assert(kt>=0 && kt < nbt );
assert(renu[kt]==-1);
@@ -3038,59 +3038,59 @@ void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress)
}
if (verbosity>9)
cout << " number of inside triangles " << k << " nbt = " << nbt << endl;
- // take is same numbering if possible
+ // take is same numbering if possible
if(justcompress)
- for ( k=0,it=0;it<nbt;it++)
- if(renu[it] >=0 )
+ for ( k=0,it=0;it<nbt;it++)
+ if(renu[it] >=0 )
renu[it]=k++;
-
+
// put the outside triangles at the end
- for ( it=0;it<nbt;it++)
- if (renu[it]==-1)
+ for ( it=0;it<nbt;it++)
+ if (renu[it]==-1)
renu[it]=k++;
-
+
assert(k == nbt);
- // do the change on all the pointeur
+ // do the change on all the pointeur
for ( it=0;it<nbt;it++)
triangles[it].ReNumbering(triangles,te,renu);
-
+
for ( i=0;i<NbSubDomains;i++)
subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
-
- // move the Triangles without a copy of the array
- // be carefull not trivial code
+
+ // move the Triangles without a copy of the array
+ // be carefull not trivial code
for ( it=0;it<nbt;it++) // for all sub cycles of the permutation renu
if (renu[it] >= 0) // a new sub cycle
- {
+ {
i=it;
Triangle ti=triangles[i],tj;
- while ( (j=renu[i]) >= 0)
- { // i is old, and j is new
- renu[i] = -1; // mark
+ while ( (j=renu[i]) >= 0)
+ { // i is old, and j is new
+ renu[i] = -1; // mark
tj = triangles[j]; // save new
triangles[j]= ti; // new <- old
- i=j; // next
+ i=j; // next
ti = tj;
- }
+ }
}
delete [] renu;
nt = nbt - NbOutT;
#ifdef DEBUG
-// verif
+// verif
for ( it=0;it<nbt;it++)
triangles[it].check();
-#endif
+#endif
}
Int4 Triangles::ConsRefTriangle(Int4 *reft) const
{
assert(reft);
register Triangle *t0,*t;
- register Int4 k=0, num;
- for (Int4 it=0;it<nbt;it++)
- reft[it]=-1; // outside triangle
+ register Int4 k=0, num;
+ for (Int4 it=0;it<nbt;it++)
+ reft[it]=-1; // outside triangle
for (Int4 i=0;i<NbSubDomains;i++)
- {
+ {
t=t0=subdomains[i].head;
assert(t0); // no empty sub domain
// register Int4 color=i+1;// because the color 0 is outside triangle
@@ -3103,11 +3103,11 @@ Int4 Triangles::ConsRefTriangle(Int4 *reft) const
while (t0 != (t=t->link));
}
// NbOutT = nbt - k;
- if (verbosity>5)
- cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k
+ if (verbosity>5)
+ cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k
<< " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl;
-
- return k;
+
+ return k;
}
@@ -3118,8 +3118,8 @@ void Triangles::ConsLinkTriangle()
subdomains[i].head=0;
register Triangle * t=triangles,*tend = triangles+nbt,*hst;
for (;t<tend;t++)
- {
- if (t->link)
+ {
+ if (t->link)
{
register Int4 color = t->color-1;
assert(color<NbSubDomains && color>=0);
@@ -3128,18 +3128,18 @@ void Triangles::ConsLinkTriangle()
hst->link=t; }
else {
subdomains[color].head = t;
- t->link=t;}// circular link
+ t->link=t;}// circular link
}
}
{
for (Int4 i=0;i<NbSubDomains;i++)
assert(subdomains[i].head);
}
-
+
}
*/
-/* void Triangles::RandomInit()
-{
+/* void Triangles::RandomInit()
+{
// Meshbegin = vertices;
// Meshend = vertices + nbvx;
nbv = nbvx;
@@ -3150,8 +3150,8 @@ void Triangles::ConsLinkTriangle()
vertices[i].ref = 0;
}
}
-void Triangles::CubeInit(int n,int m)
-{
+void Triangles::CubeInit(int n,int m)
+{
// nbvx = n*m;
// Meshbegin = vertices;
// Meshend = vertices + nbvx;
@@ -3167,8 +3167,8 @@ void Triangles::CubeInit(int n,int m)
}
}
*/
-Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j)
- { return quadtree->NearestVertex(i,j); }
+Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j)
+ { return quadtree->NearestVertex(i,j); }
void Triangles::PreInit(Int4 inbvx,char *fname)
@@ -3191,17 +3191,17 @@ void Triangles::PreInit(Int4 inbvx,char *fname)
NbVertexOnBThEdge=0;
VertexOnBThVertex=0;
VertexOnBThEdge=0;
-
+
NbCrackedVertices=0;
NbCrackedEdges =0;
- CrackedEdges =0;
- nbe = 0;
+ CrackedEdges =0;
+ nbe = 0;
name = fname ;
if (inbvx) {
vertices=new Vertex[nbvx];
assert(vertices);
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
assert(ordre);
triangles=new Triangle[nbtx];
assert(triangles);}
@@ -3214,10 +3214,10 @@ void Triangles::PreInit(Int4 inbvx,char *fname)
if ( name || inbvx)
{
time_t timer =time(0);
- char buf[70];
+ char buf[70];
strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer));
- counter++;
- char countbuf[30];
+ counter++;
+ char countbuf[30];
sprintf(countbuf,"%d",counter);
int lg =0 ;
if (&BTh != this && BTh.name)
@@ -3232,8 +3232,8 @@ void Triangles::PreInit(Int4 inbvx,char *fname)
strcat(strcat(identity,";"),countbuf);
strcat(identity,buf);
// cout << "New MAILLAGE "<< identity << endl;
- }
-
+ }
+
quadtree=0;
// edgescomponante=0;
edges=0;
@@ -3247,36 +3247,36 @@ void Triangles::PreInit(Int4 inbvx,char *fname)
NbSubDomains=0;
// Meshbegin = vertices;
// Meshend = vertices + nbvx;
- if (verbosity>98)
- cout << "Triangles::PreInit() " << nbvx << " " << nbtx
- << " " << vertices
+ if (verbosity>98)
+ cout << "Triangles::PreInit() " << nbvx << " " << nbtx
+ << " " << vertices
<< " " << ordre << " " << triangles <<endl;
}
-void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
-{
+void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
+{
//#define DRAWING1
Gh.NbRef++;// add a ref to Gh
-
-/*************************************************************************
+
+/*************************************************************************
// methode in 2 step
// 1 - compute the number of new edge to allocate
// 2 - construct the edge
- remark:
+ remark:
in this part we suppose to have a background mesh with the same
- geometry
-
- To construct the discretisation of the new mesh we have to
- rediscretize the boundary of background Mesh
+ geometry
+
+ To construct the discretisation of the new mesh we have to
+ rediscretize the boundary of background Mesh
because we have only the pointeur from the background mesh to the geometry.
We need the abcisse of the background mesh vertices on geometry
- so a vertex is
+ so a vertex is
0 on GeometricalVertex ;
1 on GeometricalEdge + abcisse
- 2 internal
-
+ 2 internal
+
*************************************************************************/
assert(&BTh.Gh == &Gh);
@@ -3287,36 +3287,36 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
#ifdef DRAWING
if (withrgraphique)
{ BTh.InitDraw();
- reffecran();
+ reffecran();
CurrentTh = this;}
#endif
- int * bcurve = new int[Gh.NbOfCurves]; //
-
- // we have 2 ways to make the loop
- // 1) on the geometry
+ int * bcurve = new int[Gh.NbOfCurves]; //
+
+ // we have 2 ways to make the loop
+ // 1) on the geometry
// 2) on the background mesh
// if you do the loop on geometry, we don't have the pointeur on background,
// and if you do the loop in background we have the pointeur on geometry
// so do the walk on background
// Int4 NbVerticesOnGeomVertex;
- // VertexOnGeom * VerticesOnGeomVertex;
+ // VertexOnGeom * VerticesOnGeomVertex;
// Int4 NbVerticesOnGeomEdge;
// VertexOnGeom * VerticesOnGeomEdge;
-
-
+
+
NbVerticesOnGeomVertex=0;
NbVerticesOnGeomEdge=0;
//1 copy of the Required vertex
- int i;
+ int i;
for ( i=0;i<Gh.nbv;i++)
if (Gh[i].Required()) NbVerticesOnGeomVertex++;
-
+
VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
//
- if( NbVerticesOnGeomVertex >= nbvx)
+ if( NbVerticesOnGeomVertex >= nbvx)
{
- cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;
+ cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;
MeshError(1,this);
}
assert(vertices);
@@ -3329,9 +3329,9 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
// cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;
nbv++;}
else Gh[i].to=0;
- //
+ //
for (i=0;i<BTh.NbVerticesOnGeomVertex;i++)
- {
+ {
VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i];
if (vog.IsRequiredVertex()) {
GeometricalVertex * gv = vog;
@@ -3345,88 +3345,88 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
// new stuff FH with curve
// find the begin of the curve in BTh
{
- Gh.UnMarkEdges();
+ Gh.UnMarkEdges();
int bfind=0;
/*
cout << " nb curves = " << Gh.NbOfCurves << endl;
for(int i=0;i<Gh.NbOfCurves ;i++)
{
- cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb
+ cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb
<< " end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl;
}*/
for (int i=0;i<Gh.NbOfCurves;i++)
{
- bcurve[i]=-1;
+ bcurve[i]=-1;
}
-
- for (int iedge=0;iedge<BTh.nbe;iedge++)
- {
+
+ for (int iedge=0;iedge<BTh.nbe;iedge++)
+ {
Edge & ei = BTh.edges[iedge];
for(int je=0;je<2;je++) // for the 2 extremites
if (!ei.on->Mark() && ei[je].on->IsRequiredVertex() )
{
- // a begin of curve
+ // a begin of curve
int nc = ei.on->CurveNumber;
-
+
//cout << "curve " << nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " "
// << " e " << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl;
if(
- ei.on==Gh.curves[nc].be &&
- (GeometricalVertex *) *ei[je].on == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity
- )
- {
+ ei.on==Gh.curves[nc].be &&
+ (GeometricalVertex *) *ei[je].on == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity
+ )
+ {
// cout << " find " << endl;
bcurve[nc]=iedge*2+je;
- bfind++;
- }
+ bfind++;
+ }
}
- }
+ }
assert( bfind==Gh.NbOfCurves);
-}
-// method in 2 + 1 step
+}
+// method in 2 + 1 step
// 0.0) compute the length and the number of vertex to do allocation
// 1.0) recompute the length
-// 1.1) compute the vertex
+// 1.1) compute the vertex
Int4 nbex=0,NbVerticesOnGeomEdgex=0;
for (int step=0; step <2;step++)
{
Int4 NbOfNewPoints=0;
Int4 NbOfNewEdge=0;
Int4 iedge;
- Gh.UnMarkEdges();
-/* add Curve loop FH
- // find a starting back groud edges to walk
- for (iedge=0;iedge<BTh.nbe;iedge++) {
+ Gh.UnMarkEdges();
+/* add Curve loop FH
+ // find a starting back groud edges to walk
+ for (iedge=0;iedge<BTh.nbe;iedge++) {
Edge & ei = BTh.edges[iedge];
for(int jedge=0;jedge<2;jedge++) // for the 2 extremites
if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() )
{
-*/
-// new code FH 2004
+*/
+// new code FH 2004
Real8 L=0;
for (int icurve=0;icurve<Gh.NbOfCurves;icurve++)
- {
+ {
iedge=bcurve[icurve]/2;
int jedge=bcurve[icurve]%2;
if( ! Gh.curves[icurve].master) continue; // we skip all equi curve
Edge & ei = BTh.edges[iedge];
// warning: ei.on->Mark() can be change in
- // loop for(jedge=0;jedge<2;jedge++)
- // new curve
- // good the find a starting edge
- Real8 Lstep=0,Lcurve=0;// step between two points (phase==1)
- Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
+ // loop for(jedge=0;jedge<2;jedge++)
+ // new curve
+ // good the find a starting edge
+ Real8 Lstep=0,Lcurve=0;// step between two points (phase==1)
+ Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
+
-
// cout.precision(16);
- for(int phase=0;phase<=step;phase++)
+ for(int phase=0;phase<=step;phase++)
{
for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
{
int icurveequi= Gh.Number(curve);
-
+
if( phase == 0 && icurveequi != icurve) continue;
int k0=jedge,k1;
@@ -3435,11 +3435,11 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
int iedgeequi=bcurve[icurveequi]/2;
int jedgeequi=bcurve[icurveequi]%2;
- int k0equi=jedgeequi,k1equi;
+ int k0equi=jedgeequi,k1equi;
Edge * peequi= BTh.edges+iedgeequi;
GeometricalEdge *ongequi = peequi->on;
-
- Real8 sNew=Lstep;// abcisse of the new points (phase==1)
+
+ Real8 sNew=Lstep;// abcisse of the new points (phase==1)
L=0;// length of the curve
Int4 i=0;// index of new points on the curve
register GeometricalVertex * GA0 = *(*peequi)[k0equi].on;
@@ -3448,38 +3448,38 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
Vertex *A1;
VertexOnGeom *GA1;
Edge * PreviousNewEdge = 0;
- // cout << " --------------New Curve phase " << phase
+ // cout << " --------------New Curve phase " << phase
// << "---------- A0=" << *A0 << ei[k0] <<endl;
assert (A0-vertices>=0 && A0-vertices <nbv);
- if(ongequi->Required() )
+ if(ongequi->Required() )
{
GeometricalVertex *GA1 = *(*peequi)[1-k0equi].on;
A1 = GA1->to; //
- }
- else
- for(;;)
+ }
+ else
+ for(;;)
{
// assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe);
- Edge &ee=*pe;
- Edge &eeequi=*peequi;
- k1 = 1-k0; // next vertex of the edge
+ Edge &ee=*pe;
+ Edge &eeequi=*peequi;
+ k1 = 1-k0; // next vertex of the edge
k1equi= 1 - k0equi;
-
+
assert(pe && ee.on);
ee.on->SetMark();
Vertex & v0=ee[0], & v1=ee[1];
R2 AB= (R2) v1 - (R2) v0;
Real8 L0=L,LAB;
LAB = LengthInterpole(v0.m,v1.m,AB);
- L+= LAB;
+ L+= LAB;
if (phase) {// computation of the new points
- while ((i!=NbCreatePointOnCurve) && sNew <= L) {
- // cout << " L0= " << L0 << " L " << L << " sN="
+ while ((i!=NbCreatePointOnCurve) && sNew <= L) {
+ // cout << " L0= " << L0 << " L " << L << " sN="
// << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i << endl;
assert (sNew >= L0);
assert(LAB);
-
-
+
+
assert(vertices && nbv<nbvx);
assert(edges && nbe < nbex);
assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex);
@@ -3498,8 +3498,8 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
//((k1==1) != (k1==k1equi))
se = k1 ? se : 1. - se;
se = k1==k1equi ? se : 1. - se;
- VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
- ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
+ VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
+ ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
A1->ReferenceNumber = eeequi.ref;
A1->DirOfSearch =NoDirOfSearch;
//cout << icurveequi << " " << i << " " << *A1 << endl;
@@ -3507,10 +3507,10 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
e->v[0]= A0;
e->v[1]= A1;
if(verbosity>99)
- cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L0
- << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "
+ cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L0
+ << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "
<< se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl;
-
+
#ifdef DEBUG
// code \label(xxx)
R2 A1A0 = A1->r - A0->r;
@@ -3518,15 +3518,15 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
if (dp > 1.4) {
cerr << " PB new Points "<< nbv-1 ;
cerr << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= " ;
- cerr << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<endl;
- cerr << " PB calcul new on cuver points trop loin l=" << dp
+ cerr << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<endl;
+ cerr << " PB calcul new on cuver points trop loin l=" << dp
<< " v=" << nbv-1 << " " << nbv-2 << " Lcurve = " << Lcurve << AB <<v0.m<< v1.m <<endl;
}
#endif
e->ref = eeequi.ref;
e->adj[0]=PreviousNewEdge;
-
+
if (PreviousNewEdge)
PreviousNewEdge->adj[1] = e;
PreviousNewEdge = e;
@@ -3538,12 +3538,12 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
#endif
A0=A1;
sNew += Lstep;
- // cout << " sNew = " << sNew << " L = " << L
+ // cout << " sNew = " << sNew << " L = " << L
// << " ------" <<NbCreatePointOnCurve << " " << i << endl;
if (++i== NbCreatePointOnCurve) break;
}
-
- }
+
+ }
assert(ee.on->CurveNumber==ei.on->CurveNumber);
if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl;
if ( ee[k1].on->IsRequiredVertex()) {
@@ -3553,22 +3553,22 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
assert (A1-vertices>=0 && A1-vertices <nbv);
break;}
if (!ee.adj[k1])
- {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe
- << " Gh.vertices " << Gh.vertices
+ {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe
+ << " Gh.vertices " << Gh.vertices
<< " k1 = " << k1 << " on=" << *ee[k1].on << endl;
cerr << ee[k1].on->gv-Gh.vertices << endl;
}
pe = ee.adj[k1]; // next edge
- k0 = pe->Intersection(ee);
+ k0 = pe->Intersection(ee);
peequi= eeequi.adj[k1equi]; // next edge
- k0equi=peequi->Intersection(eeequi);
+ k0equi=peequi->Intersection(eeequi);
}// for(;;) end of the curve
-
+
if (phase) // construction of the last edge
{
Edge *e = edges + nbe++;
- if (verbosity>10)
+ if (verbosity>10)
cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" <<
NbCreatePointOnCurve << " == " <<i<<endl;
e->on = ongequi;
@@ -3580,50 +3580,50 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
if (PreviousNewEdge)
PreviousNewEdge->adj[1] = e;
PreviousNewEdge = e;
- // cout << "Last new edge " << nbe << " " << " on " << Gh.Number(pe->on)
+ // cout << "Last new edge " << nbe << " " << " on " << Gh.Number(pe->on)
// << " of curve =" <<pe->on->CurveNumber <<endl;
-
-#ifdef DRAWING1
+
+#ifdef DRAWING1
e->Draw();
A1->Draw();
A0->Draw();
// inquire();
#endif
assert(i==NbCreatePointOnCurve);
-
+
}
- } // end loop on equi curve
-
- if (!phase) { //
+ } // end loop on equi curve
+
+ if (!phase) { //
Int4 NbSegOnCurve = Max((Int4)(L+0.5),(Int4) 1);// nb of seg
- Lstep = L/NbSegOnCurve;
+ Lstep = L/NbSegOnCurve;
Lcurve = L;
NbCreatePointOnCurve = NbSegOnCurve-1;
-
+
for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
{
NbOfNewEdge += NbSegOnCurve;
NbOfNewPoints += NbCreatePointOnCurve;
}
if(verbosity>5)
- cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep="
+ cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep="
<< Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl;
- // do'nt
+ // do'nt
// if(NbCreatePointOnCurve<1) break;
}
}//for(phase;;)
-/* modif FH add Curve class
- }}//for (iedge=0;iedge<BTh.nbe;iedge++)
+/* modif FH add Curve class
+ }}//for (iedge=0;iedge<BTh.nbe;iedge++)
*/
-// new code Curve class
- } // end of curve loop
-// end new code
+// new code Curve class
+ } // end of curve loop
+// end new code
// do the allocation
- if(step==0)
+ if(step==0)
{
- //if(!NbOfNewPoints) break;// nothing ????? bug
- if(nbv+NbOfNewPoints > nbvx)
+ //if(!NbOfNewPoints) break;// nothing ????? bug
+ if(nbv+NbOfNewPoints > nbvx)
{
cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl;
MeshError(3,this);
@@ -3631,7 +3631,7 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
//cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl;
edges = new Edge[NbOfNewEdge];
nbex = NbOfNewEdge;
- if(NbOfNewPoints) { //
+ if(NbOfNewPoints) { //
VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
NbVertexOnBThEdge =NbOfNewPoints;
VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];
@@ -3643,19 +3643,19 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
assert(nbe);
delete [] bcurve;
-
-
+
+
#ifdef DRAWING1
reffecran();
InitDraw();
Draw();
inquire();
#endif
-
+
Insert();
ForceBoundary();
FindSubDomain();
-
+
#ifdef DRAWING1
reffecran();
Draw();
@@ -3667,10 +3667,10 @@ void Triangles::GeomToTriangles1(Int4 inbvx,int KeepBackVertices)
// BTh[iv].i = toI2(BTh[iv].r);
NewPoints(BTh,KeepBackVertices) ;
CurrentTh = 0;
-//#undef DRAWING1
+//#undef DRAWING1
}
-void Triangles::GeomToTriangles0(Int4 inbvx)
+void Triangles::GeomToTriangles0(Int4 inbvx)
{
Gh.NbRef++;// add a ref to GH
@@ -3683,7 +3683,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
Gh.InitDraw() ;
CurrentTh = this; }
#endif
-
+
R2 AB;
GeometricalVertex *a,*b;
Vertex *va,*vb;
@@ -3696,9 +3696,9 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
NbVerticesOnGeomEdge=0;
for (i=0;i<Gh.nbv;i++)
if (Gh[i].Required() && Gh[i].IsThe() ) NbVerticesOnGeomVertex++;
- VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
+ VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
//
- if( NbVerticesOnGeomVertex >= nbvx)
+ if( NbVerticesOnGeomVertex >= nbvx)
{
cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl;
MeshError(1,this);
@@ -3713,35 +3713,35 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
nbv++;
}
// assert( Gh.nbv < nbvx);
-
- // Method in 2 step: 0 and 1
- // 1) compute de nb of edge
- // 2) construct the edge
+
+ // Method in 2 step: 0 and 1
+ // 1) compute de nb of edge
+ // 2) construct the edge
// generation of the curves
assert(! edges);
#ifdef DRAWING1
reffecran();
#endif
-// 2 step
+// 2 step
// --step=0 to compute the number of edges + alloc at end
// --step=1 to construct the edges
- for (int step=0;step<2;step++)
- {// for (int step=0;step<2;step++)
+ for (int step=0;step<2;step++)
+ {// for (int step=0;step<2;step++)
Int4 nbex = 0;
nbe = 0;
Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
// cout << " -------------- step =" << step << endl;
- Gh.UnMarkEdges();
+ Gh.UnMarkEdges();
NbOfCurves = 0;
for (i=0;i<Gh.nbe;i++) {
- GeometricalEdge & ei = Gh.edges[i];
+ GeometricalEdge & ei = Gh.edges[i];
if (!ei.Dup()) // a good curve (not dup )
- for(int j=0;j<2;j++)
- if (!ei.Mark() && ei[j].Required()) {
- // warning ei.Mark() can be change in loop for(j=0;j<2;j++)
+ for(int j=0;j<2;j++)
+ if (!ei.Mark() && ei[j].Required()) {
+ // warning ei.Mark() can be change in loop for(j=0;j<2;j++)
// cout << " New curve = " << NbOfCurves << endl;
Int4 nbvend =0;
-
+
Edge * PreviousNewEdge=0;
lstep = -1;//to do not create points
@@ -3751,7 +3751,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
if(step==0)
nbe++;
else
- {
+ {
e = & ei;
a=ei(0)->The();
b=ei(1)->The();
@@ -3767,7 +3767,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
#endif
nbe++;}
}
- else
+ else
{ // on curve ------
for ( int kstep=0;kstep<= step;kstep++)
{ // begin for ( int kstep=0;kstep<= step;kstep++)
@@ -3777,7 +3777,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
PreviousNewEdge=0;
NbNewPoints=0;
NbEdgeCurve=0;
- assert(nbvend < nbvx);
+ assert(nbvend < nbvx);
lcurve =0;
s = lstep;
int k=j;
@@ -3786,32 +3786,32 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
va = a->to;
e->SetMark();
// cout << " New curve " ;
-
- // if SameGeo We have go in the background geometry
+
+ // if SameGeo We have go in the background geometry
// to find the discretisation of the curve
-
- for(;;)
- {
+
+ for(;;)
+ {
k = 1-k;
b= (*e)(k)->The();
AB = b->r - a->r;
Metric MA = background ? BTh.MetricAt(a->r) :a->m ;
Metric MB = background ? BTh.MetricAt(b->r) :b->m ;
Real8 ledge = (MA(AB) + MB(AB))/2;
- //
+ //
const int MaxSubEdge = 10;
int NbSubEdge = 1;
Real8 lSubEdge[MaxSubEdge];
R2 A,B;
- if (ledge < 1.5)
+ if (ledge < 1.5)
lSubEdge[0] = ledge;
else {
NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
A= a->r;
Metric MAs =MA,MBs;
- // cout << " lSubEdge old=" << ledge
+ // cout << " lSubEdge old=" << ledge
// << " new " << A << MA << endl;
- ledge = 0;
+ ledge = 0;
Real8 x =0, xstep= 1. / NbSubEdge;
for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
x += xstep;
@@ -3819,20 +3819,20 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB);
AB = A-B;
lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
- // cout << " " << lSubEdge[kk] << " x " << x
+ // cout << " " << lSubEdge[kk] << " x " << x
// << " " << A << B << MA << MB<< endl ;
}
// cout << endl;
}
-
+
Real8 lcurveb = lcurve+ ledge ;
while (lcurve<=s && s <= lcurveb && nbv < nbvend)
{
// New points
-
+
// Real8 aa=(lcurveb-s)/ledge;
// Real8 bb=(s-lcurve)/ledge;
-
+
Real8 ss = s-lcurve;
// 1) find the SubEdge containing ss by dichotomie
int kk0=-1,kk1=NbSubEdge-1,kkk;
@@ -3844,10 +3844,10 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
else
kk0=kkk,ll0=llk;}
assert(kk1 != kk0);
-
+
Real8 sbb = (ss-ll0 )/(ll1-ll0);
Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;
-
+
// new vertex on edge
vb = &vertices[nbv++];
vb->m = Metric(aa,a->m,bb,b->m);
@@ -3855,10 +3855,10 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
vb->DirOfSearch =NoDirOfSearch;
Real8 abcisse = k ? bb : aa;
vb->r = e->F( abcisse );
- VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);
-
+ VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);
+
// to take in account the sens of the edge
-
+
s += lstep;
edges[nbe].v[0]=va;
edges[nbe].v[1]=vb;
@@ -3873,25 +3873,25 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
#endif
PreviousNewEdge = edges + nbe;
nbe++;
-#ifdef DEBUG1
+#ifdef DEBUG1
cout << " new points " << nbv-1 << " " << vb->r ;
cout << " new edge " << nbe-1 << " " ;
- cout << va << vb << " kk0 = " << kk0
+ cout << va << vb << " kk0 = " << kk0
<< " " << kk1 << " ss=" << ss ;
cout << " " << sbb << endl;
- cout << " " << aa << va->r << bb << vb->r
+ cout << " " << aa << va->r << bb << vb->r
<<" length=" << Norme(va->r-vb->r) << endl;
- cout << " s " << s << " lstep= " << lstep
- << " ledge= " << ledge
+ cout << " s " << s << " lstep= " << lstep
+ << " ledge= " << ledge
<< " lcurve= " << lcurve << endl;
#endif
va = vb;
}
lcurve = lcurveb;
e->SetMark();
- // cout << e-Gh.edges << ", " << k << " "
- // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" "
- // << lcurve<< ";; " ;
+ // cout << e-Gh.edges << ", " << k << " "
+ // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" "
+ // << lcurve<< ";; " ;
a=b;
if (b->Required() ) break;
int kprev=k;
@@ -3906,15 +3906,15 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
if(!kstep)
{ NbVerticesOnGeomEdge0 += NbNewPoints;
NbOfCurves++;}
-
- nbvend=nbv+NbNewPoints;
-
+
+ nbvend=nbv+NbNewPoints;
+
lstep = lcurve / NbEdgeCurve;
- // cout <<"lstep " << lstep << " lcurve "
+ // cout <<"lstep " << lstep << " lcurve "
// << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<< endl;
- }
+ }
// end of curve --
- if (edges) { // last edges of the curves
+ if (edges) { // last edges of the curves
edges[nbe].v[0]=va;
edges[nbe].v[1]=vb;
edges[nbe].ref = e->ref;
@@ -3923,8 +3923,8 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
edges[nbe].adj[1] = 0;
if(PreviousNewEdge)
PreviousNewEdge->adj[1] = & edges[nbe];
-
-
+
+
#ifdef DRAWING1
edges[nbe].Draw();
#endif
@@ -3932,7 +3932,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
else
nbe += NbEdgeCurve;
} // end on curve ---
- } // if (edges[i][j].Corner())
+ } // if (edges[i][j].Corner())
} // for (i=0;i<nbe;i++)
if(!step) {
// cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl;
@@ -3944,11 +3944,11 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
assert(edges);
assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0);
// do the vertex on a geometrical vertex
- NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
+ NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
}
- else
+ else
assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0);
- // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe
+ // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe
// << "== " << nbex << " nbv = " << nbv << endl;
assert(nbex=nbe);
} // for (step=0;step<2;step++)
@@ -3959,7 +3959,7 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
Draw();
inquire();
#endif
-
+
Insert();
ForceBoundary();
FindSubDomain();
@@ -3977,42 +3977,42 @@ void Triangles::GeomToTriangles0(Int4 inbvx)
Edge** Triangles::MakeGeometricalEdgeToEdge()
{
assert(Gh.nbe);
- Edge **e= new Edge* [Gh.nbe];
-
+ Edge **e= new Edge* [this->Gh.nbe];
+
Int4 i;
for ( i=0;i<Gh.nbe ; i++)
e[i]=NULL;
- for ( i=0;i<nbe ; i++)
- {
+ for ( i=0;i<nbe ; i++)
+ {
Edge * ei = edges+i;
- GeometricalEdge *on = ei->on;
- e[Gh.Number(on)] = ei;
+ GeometricalEdge *on = ei->on;
+ e[Gh.Number(on)] = ei;
}
- for ( i=0;i<nbe ; i++)
- for (int ii=0;ii<2;ii++) {
+ for ( i=0;i<nbe ; i++)
+ for (int ii=0;ii<2;ii++) {
Edge * ei = edges+i;
GeometricalEdge *on = ei->on;
int j= ii;
- while (!(*on)[j].Required()) {
- // cout << i << " " << ii << " j= " << j << " curve = "
- // << on->CurveNumber << " " << Gh.Number(on) << " on " << j
- // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1])
+ while (!(*on)[j].Required()) {
+ // cout << i << " " << ii << " j= " << j << " curve = "
+ // << on->CurveNumber << " " << Gh.Number(on) << " on " << j
+ // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1])
// << " -> " ;
Adj(on,j); // next geom edge
j=1-j;
- // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)]
- // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl;
- if (e[Gh.Number(on)]) break; // optimisation
- e[Gh.Number(on)] = ei;
+ // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)]
+ // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl;
+ if (e[Gh.Number(on)]) break; // optimisation
+ e[Gh.Number(on)] = ei;
}
}
int kk=0;
for ( i=0;i<Gh.nbe ; i++)
- if (!e[i])
+ if (!e[i])
if(kk++<10) {
- cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber
- << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl;
+ cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber
+ << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl;
// assert( e[i]);
}
if(kk) MeshError(997,this);
@@ -4020,7 +4020,7 @@ Edge** Triangles::MakeGeometricalEdgeToEdge()
return e;
}
-Triangles::~Triangles()
+Triangles::~Triangles()
{
assert(NbRef<=0);
if (CurrentTh == this) CurrentTh=0;
@@ -4038,8 +4038,8 @@ Triangles::~Triangles()
if (identity) delete [] identity;
if (VertexOnBThVertex) delete [] VertexOnBThVertex;
if (VertexOnBThEdge) delete [] VertexOnBThEdge;
-
- if (&Gh)
+
+ if (&Gh)
{
if (Gh.NbRef>0) Gh.NbRef--;
else if (Gh.NbRef==0) delete &Gh;
@@ -4049,8 +4049,8 @@ Triangles::~Triangles()
if (BTh.NbRef>0) BTh.NbRef--;
else if (BTh.NbRef==0) delete &BTh;
}
- PreInit(0); // set all to zero
-
+ PreInit(0); // set all to zero
+
}
void Triangles::SetIntCoor(const char * strfrom)
@@ -4068,13 +4068,13 @@ void Triangles::SetIntCoor(const char * strfrom)
}
R2 DD = (pmax-pmin)*0.05;
pmin = pmin-DD;
- pmax = pmax+DD;
+ pmax = pmax+DD;
coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
assert(coefIcoor >0);
- // generation of integer coord
+ // generation of integer coord
for (i=0;i<nbv;i++) {
- vertices[i].i = toI2(vertices[i].r);
+ vertices[i].i = toI2(vertices[i].r);
}
#ifdef DRAWING
xGrafCoef = coefIcoor;
@@ -4086,7 +4086,7 @@ void Triangles::SetIntCoor(const char * strfrom)
#endif
#endif
- // computation of the det
+ // computation of the det
int Nberr=0;
for (i=0;i<nbt;i++)
{
@@ -4100,47 +4100,47 @@ void Triangles::SetIntCoor(const char * strfrom)
{
if(Nberr==1)
if (strfrom)
- cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl;
- else
+ cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl;
+ else
cerr << "+++ Fatal Error Triangle (in SetInCoor) area of Triangle < 0" << endl;
cerr << " Triangle " << i << " det (I2) = " << triangles[i].det ;
cerr << " (R2) " << Det(v1.r-v0.r,v2.r-v0.r);
cerr << "; The 3 vertices " << endl;
- cerr << Number(v0) << " " << Number(v1) << " "
+ cerr << Number(v0) << " " << Number(v1) << " "
<< Number(v2) << " : " ;
cerr << v0.r << v1.r << v2.r << " ; ";
cerr << v0.i << v1.i << v2.i << endl;
}
}
else
- triangles[i].det= -1; // Null triangle;
+ triangles[i].det= -1; // Null triangle;
}
if (Nberr) MeshError(899,this);
-
+
}
-void Triangles::FillHoleInMesh()
+void Triangles::FillHoleInMesh()
{
Triangles * OldCurrentTh =CurrentTh;
CurrentTh=this;
// Int4 NbTold = nbt;
// generation of the integer coor
{
-
+
// double coef = coefIcoor;
// recherche des extrema des vertices pmin,pmax
Int4 i;
if(verbosity>2)
- cout << " -- FillHoleInMesh: Nb of vertices =" << nbv
+ cout << " -- FillHoleInMesh: Nb of vertices =" << nbv
<< " Pmin = "<< pmin << " Pmax = "<< pmax << endl;
-
+
assert(ordre);
- for (i=0;i<nbv;i++)
+ for (i=0;i<nbv;i++)
ordre[i]= 0 ;
-
+
NbSubDomains =0;
-
+
// generation of the adjacence of the triangles
SetOfEdges4 * edge4= new SetOfEdges4(nbt*3,nbv);
Int4 * st = new Int4[nbt*3];
@@ -4150,7 +4150,7 @@ void Triangles::FillHoleInMesh()
for (i=0;i<nbe;i++)
kk += (i == edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])));
if (kk != nbe)
- {
+ {
cerr << " Some Double edge in the mesh, the number is " << kk-nbe << endl;
MeshError(1002,this);
}
@@ -4165,7 +4165,7 @@ void Triangles::FillHoleInMesh()
st[k]=3*i+j;
else if(st[k]>=0) {
assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)));
-
+
triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3));
if (invisible) triangles[i].SetHidden(j);
if (k<nbe) {
@@ -4173,16 +4173,16 @@ void Triangles::FillHoleInMesh()
}
st[k]=-2-st[k]; }
else {
- cerr << " The edge ("
+ cerr << " The edge ("
<< Number(triangles[i][VerticesOfTriangularEdge[j][0]])
- << " , "
+ << " , "
<< Number(triangles[i][VerticesOfTriangularEdge[j][1]])
<< " ) is in more than 2 triangles " <<k <<endl;
cerr << " Edge " << j << " Of Triangle " << i << endl;
cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3 << endl;
cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3))
<< " Of Triangle " << Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) << endl;
- MeshError(9999,this);}
+ MeshError(9999,this);}
}
@@ -4193,45 +4193,45 @@ void Triangles::FillHoleInMesh()
cout << " - The number of given edge = " << nbe << endl;
cout << " - The number of all edges = " << edge4->nb() << endl;
cout << " - The Euler number = 1-Nb Of Hole = " << nbt-edge4->nb()+nbv << endl; }
-
-
+
+
// check the consistant of edge[].adj and the geometrical required vertex
Int4 k=0;
for (i=0;i<edge4->nb();i++)
- if (st[i] >=0) // edge alone
- if (i < nbe)
+ if (st[i] >=0) // edge alone
+ if (i < nbe)
{
Int4 i0=edge4->i(i);ordre[i0] = vertices+i0;
Int4 i1=edge4->j(i);ordre[i1] = vertices+i1;
}
else {
k++;
- if (verbosity>20 && k <20)
+ if (verbosity>20 && k <20)
{
Int4 i0=edge4->i(i);
Int4 i1=edge4->j(i);
cerr << " Lose boundary edges " << i << " : " << i0 << " " << i1 << endl;
}
}
-
+
if(k != 0) {
if (verbosity>20)
{
cout << " The given edge are " << endl;
for (int i=0;i< nbe;i++)
- cout << " Edge " << i << " : " << Number(edges[i][0]) << " " << Number(edges[i][1])
- << " " << edges[i].ref << endl;
+ cout << " Edge " << i << " : " << Number(edges[i][0]) << " " << Number(edges[i][1])
+ << " " << edges[i].ref << endl;
}
cerr << k << " boundary edges are not defined as edges " << endl;
MeshError(9998,this);
}
- // generation of the mesh with boundary points
+ // generation of the mesh with boundary points
Int4 nbvb = 0;
for (i=0;i<nbv;i++)
- {
+ {
vertices[i].t=0;
vertices[i].vint=0;
- if (ordre[i])
+ if (ordre[i])
ordre[nbvb++] = ordre[i];
}
@@ -4244,13 +4244,13 @@ void Triangles::FillHoleInMesh()
Int4 Nbtriafillhole = 2*nbvb;
Triangle * triafillhole =new Triangle[Nbtriafillhole];
if (verbosity>9)
- cout << " Nbtriafillhole triafillhole*" << triafillhole << endl;
+ cout << " Nbtriafillhole triafillhole*" << triafillhole << endl;
triangles = triafillhole;
nbt=2;
nbtx= Nbtriafillhole;
-
- for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
+
+ for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;)
if ( ++i >= nbvb) {
cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;
MeshError(998,this); }
@@ -4259,11 +4259,11 @@ void Triangles::FillHoleInMesh()
Vertex * v0=ordre[0], *v1=ordre[1];
- triangles[0](0) = 0; // sommet pour infini
+ triangles[0](0) = 0; // sommet pour infini
triangles[0](1) = v0;
triangles[0](2) = v1;
-
- triangles[1](0) = 0;// sommet pour infini
+
+ triangles[1](0) = 0;// sommet pour infini
triangles[1](2) = v0;
triangles[1](1) = v1;
const int e0 = OppositeEdge[0];
@@ -4272,40 +4272,40 @@ void Triangles::FillHoleInMesh()
triangles[0].SetAdj2(e0, &triangles[1] ,e0);
triangles[0].SetAdj2(e1, &triangles[1] ,e2);
triangles[0].SetAdj2(e2, &triangles[1] ,e1);
-
+
triangles[0].det = -1; // faux triangles
triangles[1].det = -1; // faux triangles
-
+
triangles[0].SetTriangleContainingTheVertex();
triangles[1].SetTriangleContainingTheVertex();
-
+
triangles[0].link=&triangles[1];
triangles[1].link=&triangles[0];
-
-#ifdef DEBUG
+
+#ifdef DEBUG
triangles[0].check();
triangles[1].check();
-#endif
+#endif
// nbtf = 2;
- if ( !quadtree )
+ if ( !quadtree )
delete quadtree; // ->ReInitialise();
-
+
quadtree = new QuadTree(this,0);
quadtree->Add(*v0);
quadtree->Add(*v1);
-
- // on ajoute les sommets un a un
+
+ // on ajoute les sommets un a un
Int4 NbSwap=0;
for (Int4 icount=2; icount<nbvb; icount++) {
-
+
Vertex *vi = ordre[icount];
// cout << " Add vertex " << Number(vi) << endl;
Icoor2 dete[3];
Triangle *tcvi = FindTriangleContening(vi->i,dete);
- quadtree->Add(*vi);
+ quadtree->Add(*vi);
Add(*vi,tcvi,dete);
NbSwap += vi->Optim(1,1);
-
+
#ifdef DRAWING2
cout << Number(vi) << " " << NbSwap << endl;
reffecran();
@@ -4313,16 +4313,16 @@ void Triangles::FillHoleInMesh()
vi->Draw();
inquire();
#endif
- }// end loop on icount
+ }// end loop on icount
#ifdef DRAWING1
inquire();
#endif
-
+
//Int4 nbtfillhole = nbt;
- // inforce the boundary
+ // inforce the boundary
TriangleAdjacent ta(0,0);
Int4 nbloss = 0,knbe=0;
- for ( i = 0; i < nbe; i++)
+ for ( i = 0; i < nbe; i++)
if (st[i] >=0) // edge alone => on border ... FH oct 2009
{
Vertex & a=edges[i][0], & b = edges[i][1];
@@ -4339,7 +4339,7 @@ void Triangles::FillHoleInMesh()
MeshError(1100,this);
}
FindSubDomain(1);
- // remove all the hole
+ // remove all the hole
// remove all the good sub domain
Int4 krm =0;
for (i=0;i<nbt;i++)
@@ -4350,16 +4350,16 @@ void Triangles::FillHoleInMesh()
{
TriangleAdjacent ta = triangles[i].Adj(j);
Triangle & tta = * (Triangle *) ta;
- if(! tta.link) // edge between remove and not remove
+ if(! tta.link) // edge between remove and not remove
{ // change the link of ta;
int ja = ta;
Vertex *v0= ta.EdgeVertex(0);
Vertex *v1= ta.EdgeVertex(1);
Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv);
- assert(st[k] >=0);
+ assert(st[k] >=0);
tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));
ta.SetLock();
- st[k]=-2-st[k];
+ st[k]=-2-st[k];
}
}
}
@@ -4372,15 +4372,15 @@ void Triangles::FillHoleInMesh()
else
{
triangles[i].color= savenbt+ NbTfillHoll++;
-#ifdef DEBUG
+#ifdef DEBUG
triangles[i].check();
#endif
}
// cout << savenbt+NbTfillHoll << " " << savenbtx << endl;
assert(savenbt+NbTfillHoll <= savenbtx );
- // copy of the outside triangles in saveTriangles
+ // copy of the outside triangles in saveTriangles
for (i=0;i<nbt;i++)
- if(triangles[i].color>=0)
+ if(triangles[i].color>=0)
{
savetriangles[savenbt]=triangles[i];
savetriangles[savenbt].link=0;
@@ -4389,16 +4389,16 @@ void Triangles::FillHoleInMesh()
// gestion of the adj
k =0;
Triangle * tmax = triangles + nbt;
- for (i=0;i<savenbt;i++)
- {
+ for (i=0;i<savenbt;i++)
+ {
Triangle & ti = savetriangles[i];
for (int j=0;j<3;j++)
{
Triangle * ta = ti.TriangleAdj(j);
int aa = ti.NuEdgeTriangleAdj(j);
int lck = ti.Locked(j);
- if (!ta) k++; // bug
- else if ( ta >= triangles && ta < tmax)
+ if (!ta) k++; // bug
+ else if ( ta >= triangles && ta < tmax)
{
ta= savetriangles + ta->color;
ti.SetAdj2(j,ta,aa);
@@ -4408,7 +4408,7 @@ void Triangles::FillHoleInMesh()
}
// OutSidesTriangles = triangles;
// Int4 NbOutSidesTriangles = nbt;
-
+
// restore triangles;
nbt=savenbt;
nbtx=savenbtx;
@@ -4416,7 +4416,7 @@ void Triangles::FillHoleInMesh()
delete [] subdomains;
triangles = savetriangles;
subdomains = savesubdomains;
- // cout << triangles << " <> " << OutSidesTriangles << endl;
+ // cout << triangles << " <> " << OutSidesTriangles << endl;
/* k=0;
for (i=0;i<nbt;i++)
for (int j=0;j<3;j++)
@@ -4429,18 +4429,18 @@ void Triangles::FillHoleInMesh()
}
FindSubDomain();
// cout << " NbTOld = " << NbTold << " == " << nbt - NbOutT << " " << nbt << endl;
-
- //
+
+ //
delete edge4;
delete [] st;
for (i=0;i<nbv;i++)
quadtree->Add(vertices[i]);
-
+
SetVertexFieldOn();
for (i=0;i<nbe;i++)
- if(edges[i].on)
+ if(edges[i].on)
for(int j=0;j<2;j++)
if (!edges[i].adj[j])
if(!edges[i][j].on->IsRequiredVertex()) {
@@ -4454,11 +4454,11 @@ void Triangles::FillHoleInMesh()
cerr << " = " << edges[i][j].on ;
cerr << endl;
}
-
+
#ifdef DRAWING1
InitDraw();
#endif
-
+
}
CurrentTh=OldCurrentTh;
}
@@ -4467,17 +4467,17 @@ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)
: Gh(*(pGh?pGh:&Th.Gh)), BTh(*(pBth?pBth:this))
{
Gh.NbRef++;
- nbvxx = Max(nbvxx,Th.nbv);
+ nbvxx = Max(nbvxx,Th.nbv);
Int4 i;
// do all the allocation to be sure all the pointer existe
-
+
char * cname = 0;
- if (Th.name)
+ if (Th.name)
{
cname = new char[strlen(Th.name)+1];
strcpy(cname,Th.name);
}
- PreInit(nbvxx,cname);// to make the allocation
+ PreInit(nbvxx,cname);// to make the allocation
// copy of triangles
nt=Th.nt;
nbv = Th.nbv;
@@ -4494,7 +4494,7 @@ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)
NbVerticesOnGeomEdge = Th.NbVerticesOnGeomEdge;
if (NbVerticesOnGeomEdge)
VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge] ;
- if (& BTh == & Th.BTh) // same back ground
+ if (& BTh == & Th.BTh) // same back ground
{
BTh.NbRef++;
NbVertexOnBThVertex = Th.NbVertexOnBThVertex;
@@ -4504,15 +4504,15 @@ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)
if(NbVertexOnBThEdge)
VertexOnBThEdge = new VertexOnEdge[NbVertexOnBThEdge];
}
- else
- { // no add on back ground mesh
+ else
+ { // no add on back ground mesh
BTh.NbRef++;
NbVertexOnBThVertex=0;
VertexOnBThVertex=0;
NbVertexOnBThEdge=0;
VertexOnBThEdge=0;
- // assert (& BTh == this); // --- a voir
-
+ // assert (& BTh == this); // --- a voir
+
}
@@ -4529,7 +4529,7 @@ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)
edges[i].Set(Th,i,*this);
for(i=0;i<nbv;i++)
vertices[i].Set(Th.vertices[i],Th,*this);
- for(i=0;i<NbSubDomains;i++)
+ for(i=0;i<NbSubDomains;i++)
subdomains[i].Set(Th,i,*this);
for (i=0;i<NbVerticesOnGeomVertex;i++)
VerticesOnGeomVertex[i].Set(Th.VerticesOnGeomVertex[i],Th,*this);
@@ -4545,18 +4545,18 @@ Triangles::Triangles(Triangles & Th,Geometry * pGh,Triangles * pBth,Int4 nbvxx)
Int4 Triangle::Optim(Int2 i,int koption)
{
- // turn in the positif sens around vertex s
+ // turn in the positif sens around vertex s
register Int4 NbSwap =0;
register Vertex * s = ns[i];
register Triangle * tbegin=0 , *t = this , *ttc;
register int k=0,j = EdgesVertexTriangle[i][0],jc;
tbegin=t;
do {
- k++;
+ k++;
#ifdef DEBUG
assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] );
#endif
-#ifdef DRAWING1
+#ifdef DRAWING1
t->Draw();
DrawMark( s->r);
#endif
@@ -4573,8 +4573,8 @@ Int4 Triangle::Optim(Int2 i,int koption)
t = ttc;
j = NextEdge[jc];
assert(k<20000);
- } while ( (tbegin != t));
-
+ } while ( (tbegin != t));
+
return NbSwap;
}
*/
@@ -4611,7 +4611,7 @@ Int4 Triangle::Optim(Int2 i,int koption)
assert( s == & (*t)[OppositeVertex[j]] );
#endif
}
- // end on this Triangle
+ // end on this Triangle
tp = t;
jp = NextEdge[j];
@@ -4623,7 +4623,7 @@ Int4 Triangle::Optim(Int2 i,int koption)
}
void Triangles::SmoothingVertex(int nbiter,Real8 omega )
- {
+ {
// if quatree exist remove it end reconstruct
if (quadtree) delete quadtree;
quadtree=0;
@@ -4631,32 +4631,32 @@ Int4 Triangle::Optim(Int2 i,int koption)
Triangle vide; // a triangle to mark the boundary vertex
Triangle ** tstart= new Triangle* [nbv];
Int4 i,j,k;
- // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
+ // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide
if ( this == & BTh)
for ( i=0;i<nbv;i++)
- tstart[i]=vertices[i].t;
- else
+ tstart[i]=vertices[i].t;
+ else
for ( i=0;i<nbv;i++)
tstart[i]=0;
- for ( j=0;j<NbVerticesOnGeomVertex;j++ )
+ for ( j=0;j<NbVerticesOnGeomVertex;j++ )
tstart[ Number(VerticesOnGeomVertex[j].mv)]=&vide;
- for ( k=0;k<NbVerticesOnGeomEdge;k++ )
+ for ( k=0;k<NbVerticesOnGeomEdge;k++ )
tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide;
- if(verbosity>2)
+ if(verbosity>2)
cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl;
for (k=0;k<nbiter;k++)
{
Int4 i,NbSwap =0;
Real8 delta =0;
for ( i=0;i<nbv;i++)
- if (tstart[i] != &vide) // not a boundary vertex
+ if (tstart[i] != &vide) // not a boundary vertex
delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
if (!NbOfQuad)
for ( i=0;i<nbv;i++)
- if (tstart[i] != &vide) // not a boundary vertex
+ if (tstart[i] != &vide) // not a boundary vertex
NbSwap += vertices[i].Optim(1);
if (verbosity>3)
- cout << " Move max = " << sqrt(delta) << " iteration = "
+ cout << " Move max = " << sqrt(delta) << " iteration = "
<< k << " Nb of Swap = " << NbSwap << endl;
}
@@ -4664,12 +4664,12 @@ Int4 Triangle::Optim(Int2 i,int koption)
if (quadtree) quadtree= new QuadTree(this);
}
void Triangles::MakeQuadTree()
-{
+{
if(verbosity>8)
cout << " MakeQuadTree" << endl;
if ( !quadtree ) quadtree = new QuadTree(this);
-
+
#ifdef DRAWING1
quadtree->Draw();
rattente(1);
@@ -4677,20 +4677,20 @@ void Triangles::MakeQuadTree()
quadtree->Draw();
rattente(1);
#endif
-
+
}
-void Triangles::ShowRegulaty() const// Add FH avril 2007
+void Triangles::ShowRegulaty() const// Add FH avril 2007
{
- const Real8 sqrt32=sqrt(3.)*0.5;
+ const Real8 sqrt32=sqrt(3.)*0.5;
const Real8 aireKh=sqrt32*0.5;
D2 Beq(1,0),Heq(0.5,sqrt32);
D2xD2 Br(D2xD2(Beq,Heq).t());
D2xD2 B1r(Br.inv());
/* D2xD2 BB = Br.t()*Br;
- cout << " BB = " << BB << " " << Br*B1r << endl;
+ cout << " BB = " << BB << " " << Br*B1r << endl;
MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y);
MatVVP2x2 VMM(MMM);
- cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl;
+ cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl;
*/
double gammamn=1e100,hmin=1e100;
double gammamx=0,hmax=0;
@@ -4702,7 +4702,7 @@ void Triangles::ShowRegulaty() const// Add FH avril 2007
// Real8 cf2= 6.*cf*cf;
int nt=0;
for (int it=0;it<nbt;it++)
- if ( triangles[it].link)
+ if ( triangles[it].link)
{
nt++;
Triangle &K=triangles[it];
@@ -4718,7 +4718,7 @@ void Triangles::ShowRegulaty() const// Add FH avril 2007
alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
// cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl;
Real8 betaK=0;
-
+
for (int j=0;j<3;j++)
{
Real8 he= Norme2(R2(K[j],K[(j+1)%3]));
@@ -4727,27 +4727,27 @@ void Triangles::ShowRegulaty() const// Add FH avril 2007
Vertex & v=K[j];
D2xD2 M((MetricAnIso)v);
betaK += sqrt(M.det());
-
+
D2xD2 BMB = BK.t()*M*BK;
MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y);
MatVVP2x2 VM1(M1);
- //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl;
+ //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl;
gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
- gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);
+ gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);
}
betaK *= area3;// 1/2 (somme sqrt(det))* area(K)
Marea+= betaK;
// cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl;
beta=min(beta,betaK);
beta0=max(beta0,betaK);
- }
- area*=3;
+ }
+ area*=3;
gammamn=sqrt(gammamn);
- gammamx=sqrt(gammamx);
- cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl;
- cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl;
- cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl;
- cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh)
+ gammamx=sqrt(gammamx);
+ cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl;
+ cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl;
+ cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl;
+ cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh)
<< " min "<< 1./sqrt(beta0/aireKh) << endl;
}
void Triangles::ShowHistogram() const
@@ -4760,40 +4760,40 @@ void Triangles::ShowHistogram() const
Int4 i,it,k, nbedges =0;
for (i=0;i<=kmax;i++) histo[i]=0;
for (it=0;it<nbt;it++)
- if ( triangles[it].link)
+ if ( triangles[it].link)
{
-
+
for (int j=0;j<3;j++)
{
Triangle *ta = triangles[it].TriangleAdj(j);
- if ( !ta || !ta->link || Number(ta) >= it)
- {
+ if ( !ta || !ta->link || Number(ta) >= it)
+ {
Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
if ( !& vP || !&vQ) continue;
R2 PQ = vQ.r - vP.r;
Real8 l = log(LengthInterpole(vP,vQ,PQ));
-#ifdef DRAWING2
+#ifdef DRAWING2
if (l>1.4) {
penthickness(3);
vP.MoveTo(),vQ.LineTo();
penthickness(1);
cout << " l = " << l << Number(vP) << " edge = " << Number(vQ) << endl;
}
-#endif
+#endif
nbedges++;
k = (int) ((l - lmin)*delta);
k = Min(Max(k,0L),kmax);
histo[k]++;
}
}
- }
+ }
cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <<endl;
-
+
cout << " length of edge in | % of edge | Nb of edges " << endl;
cout << " ------------------- | ---------- | ----------- " << endl;
for (i=0;i<=kmax;i++)
- {
+ {
cout << " " ;
cout.width(10);
if (i==0) cout << " 0 " ;
@@ -4803,7 +4803,7 @@ void Triangles::ShowHistogram() const
if (i==kmax) cout << " +infty " ;
else cout << exp(lmin+(i+1)/delta) ;
cout.width();cout << " | " ;
-
+
cout.precision(4);
cout.width(6);
cout << ((long) ((10000.0 * histo[i])/ nbedges))/100.0 ;
@@ -4812,20 +4812,20 @@ void Triangles::ShowHistogram() const
cout << " | " << histo[i] <<endl;
}
cout << " ------------------- | ---------- | ----------- " << endl <<endl;
-
+
}
int Triangles::Crack()
- {
- assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
+ {
+ assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
for (int i=0;i<NbCrackedEdges;i++)
CrackedEdges[i].Crack();
return NbCrackedEdges;
}
-
-int Triangles::UnCrack()
-{
- assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
+
+int Triangles::UnCrack()
+{
+ assert(NbCrackedEdges ==0 || NbCrackedVertices >0);
for (int i=0;i<NbCrackedEdges;i++)
CrackedEdges[i].UnCrack();
return NbCrackedEdges;
@@ -4847,114 +4847,114 @@ int Triangles::CrackMesh()
Edge * e = new Edge[ nbe + k];
// copy
- for (i=0;i<nbe;i++)
+ for (i=0;i<nbe;i++)
e[i] = edges[i];
delete edges;
edges = e;
const int nbe0 = nbe;
- for (k=i=0;i<nbe0;i++) // on double les arete cracked
+ for (k=i=0;i<nbe0;i++) // on double les arete cracked
if(edges[i].on->Cracked())
{
e[nbe] = e[i];
- // return the edge
+ // return the edge
e[nbe].v[0] = e[i].v[1];
e[nbe].v[1] = e[i].v[0];
- e[nbe].on = e[i].on->link ; // fqux
+ e[nbe].on = e[i].on->link ; // fqux
CrackedEdges[k++]=CrackedEdge(edges,i,nbe);
nbe++;
}
- ReMakeTriangleContainingTheVertex() ;
- //
+ ReMakeTriangleContainingTheVertex() ;
+ //
int nbcrakev =0;
Vertex *vlast = vertices + nbv;
Vertex *vend = vertices + nbvx; // end of array
- for (int iv=0;iv<nbv;iv++) // vertex
+ for (int iv=0;iv<nbv;iv++) // vertex
{
Vertex & v= vertices[iv];
- Vertex * vv = & v;
+ Vertex * vv = & v;
int kk=0; // nb cracked
- int kc=0;
- int kkk =0; // nb triangle with same number
+ int kc=0;
+ int kkk =0; // nb triangle with same number
Triangle * tbegin = v.t;
- int i = v.vint;
+ int i = v.vint;
assert(tbegin && (i >= 0 ) && (i <3));
// turn around the vertex v
TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge
int k=0;
do {
int kv = VerticesOfTriangularEdge[ta][1];
- k++;
+ k++;
Triangle * tt (ta);
- if ( ta.Cracked() )
- {
+ if ( ta.Cracked() )
+ {
TriangleAdjacent tta=(ta.Adj());
assert(tta.Cracked());
- if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart
- if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
- kk++;// number of cracked edge view
+ if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart
+ if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
+ kk++;// number of cracked edge view
}
- if ( tt->link ) { // if good triangles store the value
+ if ( tt->link ) { // if good triangles store the value
int it = Number(tt);
assert(it < nt);
- (*tt)(kv)= vv; // Change the vertex of triangle
- if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref
+ (*tt)(kv)= vv; // Change the vertex of triangle
+ if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref
// tt->SetTriangleContainingTheVertex();
kkk++;
- } else if (kk) { // crack + boundary
- if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
+ } else if (kk) { // crack + boundary
+ if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use
}
-
- ta = Next(ta).Adj();
- } while ( (tbegin != ta));
+
+ ta = Next(ta).Adj();
+ } while ( (tbegin != ta));
assert(k);
if (kc) nbcrakev++;
}
-
- if ( nbcrakev )
+
+ if ( nbcrakev )
for (int iec =0;iec < NbCrackedEdges; iec ++)
CrackedEdges[iec].Set();
-
- // set the ref
+
+ // set the ref
cout << " set the ref " << endl ;
NbCrackedVertices = nbcrakev;
// int nbvo = nbv;
nbv = vlast - vertices;
- int nbnewv = nbv - nbv; // nb of new vrtices
+ int nbnewv = nbv - nbv; // nb of new vrtices
if (nbcrakev && verbosity > 1 )
cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl;
- // all the new vertices are on geometry
+ // all the new vertices are on geometry
// BOFBO-- A VOIR
if (nbnewv)
- { //
+ { //
Int4 n = nbnewv+NbVerticesOnGeomVertex;
Int4 i,j,k;
VertexOnGeom * vog = new VertexOnGeom[n];
- for ( i =0; i<NbVerticesOnGeomVertex;i++)
+ for ( i =0; i<NbVerticesOnGeomVertex;i++)
vog[i]=VerticesOnGeomVertex[i];
delete [] VerticesOnGeomVertex;
VerticesOnGeomVertex = vog;
- // loop on cracked edge
+ // loop on cracked edge
Vertex * LastOld = vertices + nbv - nbnewv;
for (int iec =0;iec < NbCrackedEdges; iec ++)
for (k=0;k<2;k++)
{
Edge & e = *( k ? CrackedEdges[iec].a.edge : CrackedEdges[iec].b.edge);
- for (j=0;j<2;j++)
+ for (j=0;j<2;j++)
{
Vertex * v = e(j);
if ( v >= LastOld)
- { // a new vertex
- Int4 old = v->ReferenceNumber ; // the old same vertex
+ { // a new vertex
+ Int4 old = v->ReferenceNumber ; // the old same vertex
Int4 i = ( v - LastOld);
// if the old is on vertex => warning
- // else the old is on edge => ok
+ // else the old is on edge => ok
vog[i] = vog[old];
// vog[i].mv = v;
//g[i].ge = ;
//og[i].abcisse = ;
}
-
+
}
}
@@ -4962,9 +4962,9 @@ int Triangles::CrackMesh()
}
SetVertexFieldOn();
-
+
if (vlast >= vend)
- {
+ {
cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl;
MeshError(555,this);
}
@@ -4972,12 +4972,12 @@ int Triangles::CrackMesh()
CurrentTh = CurrentThOld;
return NbCrackedVertices;
}
-
+
Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
: Gh(*(new Geometry())), BTh(*this)
{ // truncature
- //
-
+ //
+
char cname[] = "trunc";
int i,k,itadj;
@@ -4992,9 +4992,9 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
for (i=0;i<Tho.nbv;i++)
refv[i]=0;
int nbNewBedge =0;
- // int nbOldBedge =0;
+ // int nbOldBedge =0;
for (i=0;i<Tho.nbt;i++)
- if( reft[i] >=0 && flag[i])
+ if( reft[i] >=0 && flag[i])
{
const Triangle & t = Tho.triangles[i];
kt++;
@@ -5020,7 +5020,7 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
}
k=0;
for (i=0;i<Tho.nbv;i++)
- if (kk[i]>=0)
+ if (kk[i]>=0)
kk[i]=k++;
cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl;
cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl;
@@ -5028,7 +5028,7 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
Int4 inbvx =k;
PreInit(inbvx,cname);
for (i=0;i<Tho.nbv;i++)
- if (kk[i]>=0)
+ if (kk[i]>=0)
{
vertices[nbv] = Tho.vertices[i];
if (!vertices[nbv].ref())
@@ -5037,7 +5037,7 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
}
assert(inbvx == nbv);
for (i=0;i<Tho.nbt;i++)
- if( reft[i] >=0 && flag[i])
+ if( reft[i] >=0 && flag[i])
{
const Triangle & t = Tho.triangles[i];
int i0 = Tho.Number(t[0]);
@@ -5048,8 +5048,8 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
// cout <<i<< " F" << flag[i] << " T " << nbt << " = " << kk[i0] << " " << kk[i1] << " " << kk[i2] ;
// cout << " OT " << i0 << " " << i1 << " " << i2 << " " << reft[i] << endl;
triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
- triangles[nbt].color = Tho.subdomains[reft[i]].ref;
- nbt++;
+ triangles[nbt].color = Tho.subdomains[reft[i]].ref;
+ nbt++;
}
assert(kt==nbt);
if (nbt ==0 && nbv ==0) {
@@ -5061,80 +5061,80 @@ Triangles::Triangles(const Triangles & Tho,const int *flag ,const int *bb)
delete [] refv;
double cutoffradian = 10.0/180.0*Pi;
ConsGeometry(cutoffradian);
- Gh.AfterRead();
+ Gh.AfterRead();
SetIntCoor();
FillHoleInMesh();
-
+
assert(NbSubDomains);
assert(subdomains[0].head && subdomains[0].head->link);
-
+
}
-
+
Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const
-{ // in: B
+{ // in: B
// out: t
// out : dete[3]
// t the triangle and s0,s1,s2 the 3 vertices of t
// in dete[3] = { det(B,s1,s2) , det(s0,B,s2), det(s0,s1,B)}
- // with det(a,b,c ) = -1 if one of 3 vertices a,b,c is NULL
- Triangle * t=0;
+ // with det(a,b,c ) = -1 if one of 3 vertices a,b,c is NULL
+ Triangle * t=0;
int j,jp,jn,jj;
- if (tstart)
+ if (tstart)
t=tstart;
- else
+ else
{
assert(quadtree);
Vertex *a = quadtree->NearestVertex(B.x,B.y) ;
-
+
if (! a || !a->t ) {
- if (a)
+ if (a)
{cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl;
cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;}
cerr << " Pb with " << B << toR2(B) << endl;
MeshError(7777);
}
assert(a>= vertices && a < vertices+nbv);
-#ifdef DRAWING1
+#ifdef DRAWING1
a->Draw();
-#endif
+#endif
// int k=0;
t = a->t;
assert(t>= triangles && t < triangles+nbt);
-
+
}
Icoor2 detop ;
- int kkkk =0; // number of test triangle
+ int kkkk =0; // number of test triangle
- while ( t->det < 0)
- { // the initial triangles is outside
+ while ( t->det < 0)
+ { // the initial triangles is outside
int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
- assert(k0>=0); // k0 the NULL vertex
+ assert(k0>=0); // k0 the NULL vertex
int k1=NextVertex[k0],k2=PreviousVertex[k0];
dete[k0]=det(B,(*t)[k1],(*t)[k2]);
- dete[k1]=dete[k2]=-1;
- if (dete[k0] > 0) // outside B
- return t;
+ dete[k1]=dete[k2]=-1;
+ if (dete[k0] > 0) // outside B
+ return t;
t = t->TriangleAdj(OppositeEdge[k0]);
assert(kkkk++ < 2);
}
jj=0;
detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
-
- while(t->det > 0 )
- {
- assert( kkkk++ < 2000 );
+
+ while(t->det > 0 )
+ {
+ assert( kkkk++ < 2000 );
j= OppositeVertex[jj];
-
+
#ifdef DRAWING1
t->Draw();
-#endif
+#endif
dete[j] = detop; //det(*b,*s1,*s2);
jn = NextVertex[j];
jp = PreviousVertex[j];
dete[jp]= det(*(*t)(j),*(*t)(jn),B);
dete[jn] = t->det-dete[j] -dete[jp];
-
+
#ifdef DEBUG
const Vertex * s0 = (*t)(0);
const Vertex * s1 = (*t)(1);
@@ -5146,14 +5146,14 @@ Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangl
#endif
// count the number k of dete <0
int k=0,ii[3];
- if (dete[0] < 0 ) ii[k++]=0;
+ if (dete[0] < 0 ) ii[k++]=0;
if (dete[1] < 0 ) ii[k++]=1;
if (dete[2] < 0 ) ii[k++]=2;
// 0 => ok
// 1 => go in way 1
// 2 => two way go in way 1 or 2 randomly
-
- if (k==0)
+
+ if (k==0)
break;
if (k == 2 && BinaryRand())
Exchange(ii[0],ii[1]);
@@ -5166,10 +5166,10 @@ Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangl
detop = -dete[OppositeVertex[jj]];
jj = j;
}
-
- if (t->det<0) // outside triangle
+
+ if (t->det<0) // outside triangle
dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
- // NbOfTriangleSearchFind += kkkk;
+ // NbOfTriangleSearchFind += kkkk;
return t;
}
diff --git a/contrib/bamg/bamglib/MeshRead.cpp b/contrib/bamg/bamglib/MeshRead.cpp
index e23f9144f9bd02355e0e3d7468c34a2baa6bcf92..aa313d224bd8483a7b3ab7f95ed73ce4dfccedd2 100644
--- a/contrib/bamg/bamglib/MeshRead.cpp
+++ b/contrib/bamg/bamglib/MeshRead.cpp
@@ -1,26 +1,26 @@
// -*- Mode : c++ -*-
//
-// SUMMARY :
-// USAGE :
-// ORG :
+// SUMMARY :
+// USAGE :
+// ORG :
// AUTHOR : Frederic Hecht
// E-MAIL : hecht@ann.jussieu.fr
//
/*
-
+
This file is part of Freefem++
-
+
Freefem++ is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.
-
+
Freefem++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
-
+
You should have received a copy of the GNU Lesser General Public License
along with Freefem++; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
@@ -49,8 +49,8 @@ static const Direction NoDirOfSearch=Direction();
void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
{
- Real8 hmin = HUGE_VAL;// the infinie value
- // Real8 MaximalAngleOfCorner = 10*Pi/180;//
+ Real8 hmin = HUGE_VAL;// the infinie value
+ // Real8 MaximalAngleOfCorner = 10*Pi/180;//
Int4 i;
Int4 dim=0;
Int4 hvertices =0;
@@ -60,8 +60,8 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
cout << " -- ReadMesh " << f_in.CurrentFile << " Version = " << Version << endl;
int field=0;
int showfield=0;
- while (f_in.cm())
- {
+ while (f_in.cm())
+ {
field=0;
char fieldname[256];
if(f_in.eof()) break;
@@ -80,16 +80,16 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
assert(dim ==2);
}
else if (!strcmp(fieldname,"Geometry"))
- {
+ {
char * fgeom ;
f_in >> fgeom;
- // if (cutoffradian>=0) => bug if change edit the file geometry
+ // if (cutoffradian>=0) => bug if change edit the file geometry
// Gh.MaximalAngleOfCorner = cutoffradian;
if (strlen(fgeom))
Gh.ReadGeometry(fgeom);
- else
- {
- // include geometry
+ else
+ {
+ // include geometry
f_in.cm();
Gh.ReadGeometry(f_in,fgeom);
}
@@ -111,7 +111,7 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
vertices[i].m = Metric(h);}
}
else if (!strcmp(fieldname,"Vertices"))
- {
+ {
assert(dim ==2);
f_in >> nbv ;
if(verbosity>3)
@@ -119,25 +119,25 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
nbvx=nbv;
vertices=new Vertex[nbvx];
assert(vertices);
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
assert(ordre);
-
+
nbiv = nbv;
for (i=0;i<nbv;i++) {
- f_in >>
- vertices[i].r.x >>
- vertices[i].r.y >>
+ f_in >>
+ vertices[i].r.x >>
+ vertices[i].r.y >>
vertices[i].ReferenceNumber ;
vertices[i].DirOfSearch =NoDirOfSearch;
vertices[i].m=M1;
vertices[i].color =0;}
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
nbt =0;
}
else if (!strcmp(fieldname,"Triangles"))
- {
+ {
if (dim !=2)
cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
if(! vertices || !triangles || !nbv )
@@ -150,9 +150,9 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
if (nbt+NbOfTria >= nbtx)
cerr << "ReadMesh: We must have 2*NbOfQuad + NbOfTria = "
<< nbt + NbOfTria<<" < 2*nbv-2 =" << nbtx << endl,
- f_in.ShowIoErr(0);
+ f_in.ShowIoErr(0);
// begintria = nbt;
- for (i=0;i<NbOfTria;i++)
+ for (i=0;i<NbOfTria;i++)
{
Int4 i1,i2,i3,iref;
Triangle & t = triangles[nbt++];
@@ -163,7 +163,7 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
// endtria=nbt;
}
else if (!strcmp(fieldname,"Quadrilaterals"))
- {
+ {
if (dim !=2)
cerr << "ReadMesh: Dimension <> 2" <<endl , f_in.ShowIoErr(0);
@@ -178,8 +178,8 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
<< nbt + 2*NbOfQuad <<" < 2*nbv-2 =" << nbtx << endl,
f_in.ShowIoErr(0);
// beginquad=nbt;
- for (i=0;i<NbOfQuad;i++)
- {
+ for (i=0;i<NbOfQuad;i++)
+ {
Int4 i1,i2,i3,i4,iref;
Triangle & t1 = triangles[nbt++];
Triangle & t2 = triangles[nbt++];
@@ -188,26 +188,26 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
t1.color=iref;
t2 = Triangle(this,i3-1,i4-1,i1-1);
t2.color=iref;
- t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
- t2.SetHidden(OppositeEdge[1]);
+ t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
+ t2.SetHidden(OppositeEdge[1]);
}
// endquad=nbt;
}
else if (!strcmp(fieldname,"VertexOnGeometricVertex"))
- {
+ {
f_in >> NbVerticesOnGeomVertex ;
if(verbosity>5)
cout << " NbVerticesOnGeomVertex = " << NbVerticesOnGeomVertex << endl
<< " Gh.vertices " << Gh.vertices << endl;
- if( NbVerticesOnGeomVertex)
+ if( NbVerticesOnGeomVertex)
{
VerticesOnGeomVertex= new VertexOnGeom[NbVerticesOnGeomVertex] ;
if(verbosity>7)
cout << " VerticesOnGeomVertex = " << VerticesOnGeomVertex << endl
<< " Gh.vertices " << Gh.vertices << endl;
assert(VerticesOnGeomVertex);
- for (Int4 i0=0;i0<NbVerticesOnGeomVertex;i0++)
- {
+ for (Int4 i0=0;i0<NbVerticesOnGeomVertex;i0++)
+ {
Int4 i1,i2;
//VertexOnGeom & v =VerticesOnGeomVertex[i0];
f_in >> i1 >> i2 ;
@@ -216,16 +216,16 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
}
}
else if (!strcmp(fieldname,"VertexOnGeometricEdge"))
- {
+ {
f_in >> NbVerticesOnGeomEdge ;
if(verbosity>3)
cout << " NbVerticesOnGeomEdge = " << NbVerticesOnGeomEdge << endl;
- if(NbVerticesOnGeomEdge)
+ if(NbVerticesOnGeomEdge)
{
VerticesOnGeomEdge= new VertexOnGeom[NbVerticesOnGeomEdge] ;
assert(VerticesOnGeomEdge);
- for (Int4 i0=0;i0<NbVerticesOnGeomEdge;i0++)
- {
+ for (Int4 i0=0;i0<NbVerticesOnGeomEdge;i0++)
+ {
Int4 i1,i2;
Real8 s;
//VertexOnGeom & v =VerticesOnGeomVertex[i0];
@@ -235,7 +235,7 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
}
}
else if (!strcmp(fieldname,"Edges"))
- {
+ {
Int4 i,j, i1,i2;
f_in >> nbe ;
edges = new Edge[nbe];
@@ -243,16 +243,16 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
cout << " Record Edges: Nb of Edge " << nbe << " edges " << edges << endl;
assert(edges);
Real4 *len =0;
- if (!hvertices)
+ if (!hvertices)
{
len = new Real4[nbv];
for(i=0;i<nbv;i++)
len[i]=0;
}
- for (i=0;i<nbe;i++)
+ for (i=0;i<nbe;i++)
{
f_in >> i1 >> i2 >> edges[i].ref ;
-
+
assert(i1 >0 && i2 >0);
assert(i1 <= nbv && i2 <= nbv);
i1--;
@@ -263,7 +263,7 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
edges[i].adj[1]=0;
R2 x12 = vertices[i2].r-vertices[i1].r;
- Real8 l12=sqrt( (x12,x12));
+ Real8 l12=sqrt( (x12,x12));
if (!hvertices) {
vertices[i1].color++;
vertices[i2].color++;
@@ -271,24 +271,24 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
len[i2] += l12;}
hmin = Min(hmin,l12);
}
- // definition the default of the given mesh size
- if (!hvertices)
+ // definition the default of the given mesh size
+ if (!hvertices)
{
- for (i=0;i<nbv;i++)
- if (vertices[i].color > 0)
+ for (i=0;i<nbv;i++)
+ if (vertices[i].color > 0)
vertices[i].m= Metric(len[i] /(Real4) vertices[i].color);
- else
+ else
vertices[i].m= Metric(hmin);
delete [] len;
}
if(verbosity>5)
cout << " hmin " << hmin << endl;
- // construction of edges[].adj
- for (i=0;i<nbv;i++)
+ // construction of edges[].adj
+ for (i=0;i<nbv;i++)
vertices[i].color = (vertices[i].color ==2) ? -1 : -2;
for (i=0;i<nbe;i++)
- for (j=0;j<2;j++)
- {
+ for (j=0;j<2;j++)
+ {
Vertex *v=edges[i].v[j];
Int4 i0=v->color,j0;
if(i0==-1)
@@ -306,23 +306,23 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
}
}
-
+
/* ne peut pas marche car il n'y a pas de flag dans les aretes du maillages
else if (!strcmp(fieldname,"RequiredEdges"))
- {
+ {
int i,j,n;
f_in >> n ;
- for (i=0;i<n;i++) {
+ for (i=0;i<n;i++) {
f_in >> j ;
assert( j <= nbe );
assert( j > 0 );
j--;
edges[j].SetRequired(); }
}
-*/
+*/
else if (!strcmp(fieldname,"EdgeOnGeometricEdge"))
- {
+ {
assert(edges);
int i1,i2,i,j;
f_in >> i2 ;
@@ -339,15 +339,15 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
cerr << " Fatal error in file " << name << " line " << f_in.LineNumber << endl;
MeshError(1);
}
-
-
+
+
edges[i-1].on = Gh.edges + j-1;
}
// cout << "end EdgeOnGeometricEdge" << endl;
}
else if (!strcmp(fieldname,"SubDomain") || !strcmp(fieldname,"SubDomainFromMesh") )
- {
-
+ {
+
f_in >> NbSubDomains ;
subdomains = new SubDomain [ NbSubDomains ];
for (i=0;i< NbSubDomains;i++)
@@ -356,27 +356,27 @@ void Triangles::Read(MeshIstream & f_in,int Version,Real8 cutoffradian)
assert (i3==3);
head --;
assert(head < nbt && head >=0);
- subdomains[i].head = triangles+head;
+ subdomains[i].head = triangles+head;
}
}
else
{ // unkown field
field = ++showfield;
- if(showfield==1) // just to show one time
+ if(showfield==1) // just to show one time
if (verbosity>5)
cout << " Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;
}
- showfield=field; // just to show one time
+ showfield=field; // just to show one time
}
-
+
if (ifgeom==0)
{
if (verbosity)
- cout << " ## Warning: Mesh without geometry we construct a geometry (theta ="
+ cout << " ## Warning: Mesh without geometry we construct a geometry (theta ="
<< cutoffradian*180/Pi << " degres )" << endl;
- ConsGeometry(cutoffradian);
+ ConsGeometry(cutoffradian);
Gh.AfterRead();
- }
+ }
}
@@ -388,41 +388,41 @@ void Triangles::Read_am_fmt(MeshIstream & f_in)
if (verbosity>1)
cout << " -- ReadMesh .am_fmt file " << f_in.CurrentFile << endl;
-
+
Int4 i;
f_in.cm() >> nbv >> nbt ;
if (verbosity>3)
cout << " nbv = " << nbv << " nbt = " << nbt << endl;
- f_in.eol() ;//
+ f_in.eol() ;//
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
-
+ ordre=new Vertex* [this->nbvx];
+
for ( i=0;i<nbt;i++)
{
Int4 i1,i2,i3;
f_in >> i1 >> i2 >> i3 ;
- triangles[i] = Triangle(this,i1-1,i2-1,i3-1);
+ triangles[i] = Triangle(this,i1-1,i2-1,i3-1);
}
- f_in.eol() ;//
-
+ f_in.eol() ;//
+
for ( i=0;i<nbv;i++)
f_in >> vertices[i].r.x >> vertices[i].r.y,
vertices[i].m = M1,vertices[i].DirOfSearch =NoDirOfSearch;
- f_in.eol() ;//
-
- for ( i=0;i<nbt;i++)
+ f_in.eol() ;//
+
+ for ( i=0;i<nbt;i++)
f_in >> triangles[i].color;
- f_in.eol() ;//
-
- for ( i=0;i<nbv;i++)
+ f_in.eol() ;//
+
+ for ( i=0;i<nbv;i++)
f_in >> vertices[i].ReferenceNumber;
-
-
+
+
}
////////////////////////
@@ -431,8 +431,8 @@ void Triangles::Read_am(MeshIstream &ff)
{
if (verbosity>1)
cout << " -- ReadMesh .am_fmt file " << ff.CurrentFile << endl;
- Metric M1(1);
-
+ Metric M1(1);
+
IFortranUnFormattedFile f_in(ff);
Int4 l=f_in.Record();
@@ -442,33 +442,33 @@ void Triangles::Read_am(MeshIstream &ff)
assert((size_t) l==nbt*sizeof(long)*4 + nbv*(2*sizeof(float)+sizeof(long)));
if (verbosity>3)
cout << " nbv = " << nbv << " nbt = " << nbt << endl;
-
+
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
-
+ ordre=new Vertex* [this->nbvx];
+
Int4 i;
for ( i=0;i<nbt;i++) {
long i1,i2,i3;
f_in >> i1 >> i2 >> i3 ;
triangles[i] = Triangle(this,i1-1,i2-1,i3-1); }
-
+
for ( i=0;i<nbv;i++) {
float x,y;
f_in >> x >> y;
vertices[i].r.x =x;
vertices[i].r.y=y;
vertices[i].m=M1;}
-
+
for ( i=0;i<nbt;i++) {
long i;
f_in >> i;
triangles[i].color=i;}
-
+
for ( i=0;i<nbv;i++) {
long i;
f_in >> i;
@@ -483,10 +483,10 @@ void Triangles::Read_nopo(MeshIstream & ff)
if (verbosity>1)
cout << " -- ReadMesh .nopo file " << ff.CurrentFile << endl;
IFortranUnFormattedFile f_in(ff);
-
-
+
+
Int4 l,i,j;
- l=f_in.Record();
+ l=f_in.Record();
l=f_in.Record();
f_in >> i;
assert(i==32);
@@ -512,7 +512,7 @@ void Triangles::Read_nopo(MeshIstream & ff)
}
if(verbosity>2)
cout << " nb de tableau associe : " << ntacm << " niveau =" << niveau << endl;
-
+
for (i=0;i<ntacm;i++)
f_in.Record();
@@ -530,8 +530,8 @@ void Triangles::Read_nopo(MeshIstream & ff)
Int4 np = nop2[21];
// Int4 nef = nop2[13];
Metric M1(1);
- if(verbosity>2)
- cout << " ndim = " << ndim << " ncopnp= " << ncopnp << " ne = " << ne
+ if(verbosity>2)
+ cout << " ndim = " << ndim << " ncopnp= " << ncopnp << " ne = " << ne
<< " ntri = " << ntria << " nquad = " << nquad << " np = " << np << endl;
nbv = np;
nbt = 2*nquad + ntria;
@@ -543,18 +543,18 @@ void Triangles::Read_nopo(MeshIstream & ff)
if( nop2[24]>=0) f_in.Record();
NbOfQuad = nquad;
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
f_in >> l;
if(verbosity>9)
cout << " Read cnop4 nb of float " << l << endl;
-
+
assert(l==2*np);
for (i=0;i<np;i++)
{ float x,y;
@@ -573,7 +573,7 @@ void Triangles::Read_nopo(MeshIstream & ff)
Int4 k=0;
Int4 nbe4 = 3*ntria + 4*nquad;
// cout << " nbv = " << nbv << " nbe4 " << nbe4 << endl;
- SetOfEdges4 * edge4= new SetOfEdges4(nbe4,nbv);
+ SetOfEdges4 * edge4= new SetOfEdges4(nbe4,nbv);
Int4 * refe = new Int4[nbe4];
Int4 kr =0;
for (i=0;i<ne;i++)
@@ -582,42 +582,42 @@ void Triangles::Read_nopo(MeshIstream & ff)
Int4 np[4],rv[4],re[4];
Int4 ncge,nmae,ndsde,npo;
f_in >> ncge >> nmae >> ndsde >> npo ;
- //cout << " element " << i << " " << ncge << " "
+ //cout << " element " << i << " " << ncge << " "
// << nmae <<" " << npo << endl;
if (ncge != 3 && ncge != 4)
{
- cerr << " read nopo type element[" << i << "] ="
+ cerr << " read nopo type element[" << i << "] ="
<< ncge << " not 3 or 4 " << endl;
MeshError(115);
}
if (npo != 3 && npo != 4)
{
- cerr << " read nopo element[" << i << "] npo = "
+ cerr << " read nopo element[" << i << "] npo = "
<< npo << " not 3 or 4 " << endl;
MeshError(115);
}
-
+
for( j=0;j<npo;j++)
{f_in >>np[j];np[j]--;}
-
- if (ncopnp !=1)
+
+ if (ncopnp !=1)
{
f_in >> npo;
if (npo != 3 || npo != 4)
{
- cerr << " read nopo type element[" << i << "]= "
+ cerr << " read nopo type element[" << i << "]= "
<< ncge << " not 3 or 4 " << endl;
MeshError(115);
}
-
+
for(j=0;j<npo;j++)
{f_in >>np[j];np[j]--;}
-
+
}
- if (nmae>0)
+ if (nmae>0)
{
- Int4 ining; // no ref
-
+ Int4 ining; // no ref
+
f_in>>ining;
if (ining==1)
MeshError(116);
@@ -626,12 +626,12 @@ void Triangles::Read_nopo(MeshIstream & ff)
f_in >>re[j];
for (j=0;j<npo;j++)
f_in >>rv[j];
-
-
+
+
// set the ref on vertex and the shift np of -1 to start at 0
for (j=0;j<npo;j++)
vertices[np[j]].ReferenceNumber = rv[j];
-
+
if (ining==2)
for (j=0;j<npo;j++)
if (re[j])
@@ -643,21 +643,21 @@ void Triangles::Read_nopo(MeshIstream & ff)
refe[edge4->addtrie(i0,i1)]=re[j];
}
}
-
- if (npo==3)
- { // triangles
- triangles[k] = Triangle(this,np[0],np[1],np[2]);
+
+ if (npo==3)
+ { // triangles
+ triangles[k] = Triangle(this,np[0],np[1],np[2]);
triangles[k].color = ndsde;
k++;
}
else if (npo==4)
- { // quad
+ { // quad
Triangle & t1 = triangles[k++];
Triangle & t2 = triangles[k++];
t1 = Triangle(this,np[0],np[1],np[2]);
t2 = Triangle(this,np[2],np[3],np[0]);
- t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
- t2.SetHidden(OppositeEdge[1]);
+ t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
+ t2.SetHidden(OppositeEdge[1]);
t1.color = ndsde;
t2.color = ndsde;
}
@@ -666,8 +666,8 @@ void Triangles::Read_nopo(MeshIstream & ff)
cerr << " read nopo type element =" << npo << " not 3 or 4 " << endl;
MeshError(114);
}
-
-
+
+
}
// cout << k << " == " << nbt << endl;
assert(k==nbt);
@@ -695,37 +695,37 @@ void Triangles::Read_nopo(MeshIstream & ff)
}
void Triangles::Read_ftq(MeshIstream & f_in)
{
- //
+ //
if (verbosity>1)
cout << " -- ReadMesh .ftq file " << f_in.CurrentFile << endl;
-
+
Int4 i,ne,nt,nq;
f_in.cm() >> nbv >> ne >> nt >> nq ;
if (verbosity>3)
cout << " nbv = " << nbv << " nbtra = " << nt << " nbquad = " << nq << endl;
nbt = nt+2*nq;
-
+
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
Int4 k=0;
-
- for ( i=0;i<ne;i++)
+
+ for ( i=0;i<ne;i++)
{
long ii,i1,i2,i3,i4,ref;
f_in >> ii;
- if (ii==3)
- { // triangles
+ if (ii==3)
+ { // triangles
f_in >> i1>> i2 >> i3 >> ref ;
- triangles[k] = Triangle(this,i1-1,i2-1,i3-1);
+ triangles[k] = Triangle(this,i1-1,i2-1,i3-1);
triangles[k++].color = ref;
}
else if (ii==4)
- { // quad
+ { // quad
f_in >> i1>> i2 >> i3 >> i4 >> ref ;
Triangle & t1 = triangles[k++];
Triangle & t2 = triangles[k++];
@@ -733,8 +733,8 @@ void Triangles::Read_nopo(MeshIstream & ff)
t1.color=ref;
t2 = Triangle(this,i3-1,i4-1,i1-1);
t2.color=ref;
- t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
- t2.SetHidden(OppositeEdge[1]);
+ t1.SetHidden(OppositeEdge[1]); // two time because the adj was not created
+ t2.SetHidden(OppositeEdge[1]);
}
else
{
@@ -760,35 +760,35 @@ void Triangles::Read_msh(MeshIstream &f_in)
Metric M1(1.);
if (verbosity>1)
cout << " -- ReadMesh .msh file " << f_in.CurrentFile << endl;
-
+
Int4 i;
f_in.cm() >> nbv >> nbt ;
while (f_in.in.peek()==' ')
f_in.in.get();
- if(isdigit(f_in.in.peek()))
+ if(isdigit(f_in.in.peek()))
f_in >> nbe;
if (verbosity>3)
cout << " nbv = " << nbv << " nbt = " << nbt << " nbe = " << nbe << endl;
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
edges = new Edge[nbe];
for ( i=0;i<nbv;i++)
{
f_in >> vertices[i].r.x >> vertices[i].r.y >> vertices[i].ReferenceNumber;
vertices[i].on=0;
vertices[i].m=M1;
- //if(vertices[i].ReferenceNumber>NbRef) NbRef=vertices[i].ReferenceNumber;
+ //if(vertices[i].ReferenceNumber>NbRef) NbRef=vertices[i].ReferenceNumber;
}
for ( i=0;i<nbt;i++)
{
Int4 i1,i2,i3,r;
f_in >> i1 >> i2 >> i3 >> r;
triangles[i] = Triangle(this,i1-1,i2-1,i3-1);
- triangles[i].color = r;
+ triangles[i].color = r;
}
for (i=0;i<nbe;i++)
{
@@ -801,7 +801,7 @@ void Triangles::Read_msh(MeshIstream &f_in)
edges[i].ref = r;
edges[i].on=0;
}
-
+
}
//////////////////////////////////////////////////
@@ -811,18 +811,18 @@ void Triangles::Read_amdba(MeshIstream &f_in )
Metric M1(1);
if (verbosity>1)
cout << " -- ReadMesh .amdba file " << f_in.CurrentFile << endl;
-
+
Int4 i;
f_in.cm() >> nbv >> nbt ;
// if (verbosity>3)
cout << " nbv = " << nbv << " nbt = " << nbt << endl;
- f_in.eol() ;//
+ f_in.eol() ;//
nbvx = nbv;
- nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
+ nbtx = 2*nbv-2; // for filling The Holes and quadrilaterals
triangles =new Triangle[nbtx];
assert(triangles);
vertices=new Vertex[nbvx];
- ordre=new Vertex* [nbvx];
+ ordre=new Vertex* [this->nbvx];
Int4 j;
for ( i=0;i<nbv;i++)
{
@@ -833,7 +833,7 @@ void Triangles::Read_amdba(MeshIstream &f_in )
vertices[j].m=M1;
vertices[j].DirOfSearch=NoDirOfSearch;
}
-
+
for ( i=0;i<nbt;i++)
{
Int4 i1,i2,i3,ref;
@@ -844,18 +844,18 @@ void Triangles::Read_amdba(MeshIstream &f_in )
triangles[j] = Triangle(this,i1-1,i2-1,i3-1);
triangles[j].color =ref;
}
- f_in.eol() ;//
-
+ f_in.eol() ;//
+
// cerr << " a faire " << endl;
//MeshError(888);
}
-Triangles::Triangles(const char * filename,Real8 cutoffradian)
+Triangles::Triangles(const char * filename,Real8 cutoffradian)
: Gh(*(new Geometry())), BTh(*this)
-{
+{
#ifdef DRAWING1
- if (!withrgraphique) {initgraphique();withrgraphique=true;}
+ if (!withrgraphique) {initgraphique();withrgraphique=true;}
#endif
// Int4 beginquad=0,begintria=0;
@@ -877,36 +877,36 @@ Triangles::Triangles(const char * filename,Real8 cutoffradian)
Int4 inbvx =0;
PreInit(inbvx,cname);
OnDisk = 1;
- // allocGeometry = &Gh; // after Preinit ;
+ // allocGeometry = &Gh; // after Preinit ;
MeshIstream f_in (filename);
-
+
if (f_in.IsString("MeshVersionFormatted"))
{
int version ;
f_in >> version ;
Read(f_in,version,cutoffradian);
}
- else {
+ else {
if (am_fmt) Read_am_fmt(f_in);
else if (am) Read_am(f_in);
else if (amdba) Read_amdba(f_in);
else if (msh) Read_msh(f_in);
else if (nopo) Read_nopo(f_in);
else if (ftq) Read_ftq(f_in);
- else
- {
+ else
+ {
cerr << " Unknown type mesh " << filename << endl;
MeshError(2);
}
ConsGeometry(cutoffradian);
- Gh.AfterRead();
+ Gh.AfterRead();
}
SetIntCoor();
FillHoleInMesh();
// Make the quad ---
-
+
}
void Geometry::ReadGeometry(const char * filename)
@@ -920,7 +920,7 @@ void Geometry::ReadGeometry(const char * filename)
-void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
+void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
{
if(verbosity>1)
cout << " -- ReadGeometry " << filename << endl;
@@ -931,8 +931,8 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
// BeginOfCurve=0;
name=new char [strlen(filename)+1];
strcpy(name,filename);
- Real8 Hmin = HUGE_VAL;// the infinie value
-// Real8 MaximalAngleOfCorner = 10*Pi/180; ;
+ Real8 Hmin = HUGE_VAL;// the infinie value
+// Real8 MaximalAngleOfCorner = 10*Pi/180; ;
Int4 hvertices =0;
Int4 i;
Int4 Version,dim=0;
@@ -940,10 +940,10 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
int showfield=0;
int NbErr=0;
- while (f_in.cm())
- {
+ while (f_in.cm())
+ {
field=0;
- // warning ruse for on allocate fiedname at each time
+ // warning ruse for on allocate fiedname at each time
char fieldname[256];
f_in.cm() >> fieldname ;
f_in.err();
@@ -958,22 +958,22 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
else if (!strcmp(fieldname,"Dimension"))
{
f_in >> dim ;
- if(verbosity>5)
- cout << " Geom Record Dimension dim = " << dim << endl;
+ if(verbosity>5)
+ cout << " Geom Record Dimension dim = " << dim << endl;
assert(dim ==2);
}
else if (!strcmp(fieldname,"hVertices"))
- {
+ {
if (nbv <=0) {
cerr<<"Error: the field Vertex is not found before hVertices " << filename<<endl;
- NbErr++;}
- if(verbosity>5)
+ NbErr++;}
+ if(verbosity>5)
cout << " Geom Record hVertices nbv=" << nbv << endl;
hvertices =1;
- for (i=0;i< nbv;i++)
+ for (i=0;i< nbv;i++)
{
Real4 h;
- f_in >> h ;
+ f_in >> h ;
vertices[i].m = Metric(h);
}
}
@@ -981,13 +981,13 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
{ hvertices =1;
if (nbv <=0) {
cerr<<"Error: the field Vertex is not found before MetricVertices " << filename<<endl;
- NbErr++;}
- if(verbosity>5)
+ NbErr++;}
+ if(verbosity>5)
cout << " Geom Record MetricVertices nbv =" << nbv << endl;
- for (i=0;i< nbv;i++)
+ for (i=0;i< nbv;i++)
{
Real4 a11,a21,a22;
- f_in >> a11 >> a21 >> a22 ;
+ f_in >> a11 >> a21 >> a22 ;
vertices[i].m = Metric(a11,a21,a22);
}
}
@@ -995,26 +995,26 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
{ hvertices =1;
if (nbv <=0) {
cerr<<"Error: the field Vertex is not found before h1h2VpVertices " << filename<<endl;
- NbErr++;}
- if(verbosity>5)
+ NbErr++;}
+ if(verbosity>5)
cout << " Geom Record h1h2VpVertices nbv=" << nbv << endl;
- for (i=0;i< nbv;i++)
+ for (i=0;i< nbv;i++)
{
Real4 h1,h2,v1,v2;
- f_in >> h1 >> h2 >>v1 >>v2 ;
+ f_in >> h1 >> h2 >>v1 >>v2 ;
vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));
}
}
else if (!strcmp(fieldname,"Vertices"))
- {
+ {
assert(dim ==2);
f_in >> nbv ;
// if(LineError) break;
nbvx = nbv;
-
+
vertices = new GeometricalVertex[nbvx];
- if(verbosity>5)
+ if(verbosity>5)
cout << " Geom Record Vertices nbv = " << nbv << "vertices = " << vertices<<endl;
assert(nbvx >= nbv);
nbiv = nbv;
@@ -1038,21 +1038,21 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
pmax.x = Max(pmax.x,vertices[i].r.x);
pmax.y = Max(pmax.y,vertices[i].r.y);
}
-
+
R2 DD05 = (pmax-pmin)*0.05;
pmin -= DD05;
pmax += DD05;
-
+
coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
assert(coefIcoor >0);
if (verbosity>5) {
cout << " Geom: min="<< pmin << "max ="<< pmax << " hmin = " << MinimalHmin() << endl;}
}
else if(!strcmp(fieldname,"MaximalAngleOfCorner")||!strcmp(fieldname,"AngleOfCornerBound"))
- {
+ {
f_in >> MaximalAngleOfCorner;
-
- if(verbosity>5)
+
+ if(verbosity>5)
cout << " Geom Record MaximalAngleOfCorner " << MaximalAngleOfCorner <<endl;
MaximalAngleOfCorner *= Pi/180;
}
@@ -1060,30 +1060,30 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
{
if (nbv <=0) {
cerr<<"Error: the field edges is not found before MetricVertices " << filename<<endl;
- NbErr++;}
- else
+ NbErr++;}
+ else
{
int i1,i2;
R2 zero2(0,0);
f_in >> nbe ;
-
+
edges = new GeometricalEdge[nbe];
- if(verbosity>5)
+ if(verbosity>5)
cout << " Record Edges: Nb of Edge " << nbe <<endl;
assert(edges);
- assert (nbv >0);
+ assert (nbv >0);
Real4 *len =0;
- if (!hvertices)
+ if (!hvertices)
{
len = new Real4[nbv];
for(i=0;i<nbv;i++)
len[i]=0;
}
-
- for (i=0;i<nbe;i++)
+
+ for (i=0;i<nbe;i++)
{
f_in >> i1 >> i2 >> edges[i].ref ;
-
+
i1--;i2--; // for C index
edges[i].v[0]= vertices + i1;
edges[i].v[1]= vertices + i2;
@@ -1094,41 +1094,41 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;
edges[i].Adj[0] = edges[i].Adj[1] = 0;
edges[i].flag = 0;
- if (!hvertices)
+ if (!hvertices)
{
vertices[i1].color++;
vertices[i2].color++;
len[i1] += l12;
len[i2] += l12;
}
-
+
Hmin = Min(Hmin,l12);
}
- // definition the default of the given mesh size
- if (!hvertices)
+ // definition the default of the given mesh size
+ if (!hvertices)
{
- for (i=0;i<nbv;i++)
- if (vertices[i].color > 0)
+ for (i=0;i<nbv;i++)
+ if (vertices[i].color > 0)
vertices[i].m= Metric(len[i] /(Real4) vertices[i].color);
- else
+ else
vertices[i].m= Metric(Hmin);
delete [] len;
-
- if(verbosity>3)
+
+ if(verbosity>3)
cout << " Geom Hmin " << Hmin << endl;
}
-
+
}
}
else if (!strcmp(fieldname,"EdgesTangence") ||!strcmp(fieldname,"TangentAtEdges") )
- {
+ {
int n,i,j,k;
R2 tg;
f_in >> n ;
-
- if(verbosity>5)
+
+ if(verbosity>5)
cout << " Record TangentAtEdges: Nb of Edge " << n <<endl;
-
+
for (k=0;k<n;k++)
{
f_in >> i >> j ;
@@ -1141,13 +1141,13 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
}
}
else if (!strcmp(fieldname,"Corners"))
- {
+ {
int i,j,n;
f_in >> n ;
- if(verbosity>5)
+ if(verbosity>5)
cout << " Record Corner: Nb of Corner " << n <<endl;
-
- for (i=0;i<n;i++) {
+
+ for (i=0;i<n;i++) {
f_in >> j ;
assert( j <= nbv );
assert( j > 0 );
@@ -1156,11 +1156,11 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
vertices[j].SetRequired(); }
}
else if (!strcmp(fieldname,"RequiredVertices"))
- {
+ {
int i,j,n;
f_in >> n ;
- for (i=0;i<n;i++) {
+ for (i=0;i<n;i++) {
f_in >> j ;
assert( j <= nbv );
assert( j > 0 );
@@ -1168,11 +1168,11 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
vertices[j].SetRequired(); }
}
else if (!strcmp(fieldname,"RequiredEdges"))
- {
+ {
int i,j,n;
f_in >> n ;
- for (i=0;i<n;i++) {
+ for (i=0;i<n;i++) {
f_in >> j ;
assert( j <= nbe );
assert( j > 0 );
@@ -1180,17 +1180,17 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
edges[j].SetRequired(); }
}
else if (!strcmp(fieldname,"SubDomain") || !strcmp(fieldname,"SubDomainFromGeom"))
- {
+ {
f_in >> NbSubDomains ;
- if (NbSubDomains>0)
+ if (NbSubDomains>0)
{
subdomains = new GeometricalSubDomain[ NbSubDomains];
Int4 i0,i1,i2,i3;
- for (i=0;i<NbSubDomains;i++)
+ for (i=0;i<NbSubDomains;i++)
{
- f_in >> i0 >>i1
- >> i2 >>i3 ;
-
+ f_in >> i0 >>i1
+ >> i2 >>i3 ;
+
assert(i0 == 2);
assert(i1<=nbe && i1>0);
subdomains[i].edge=edges + (i1-1);
@@ -1202,16 +1202,16 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
else
{ // unkown field
field = ++showfield;
- if(showfield==1) // just to show one time
+ if(showfield==1) // just to show one time
if (verbosity>3)
cout << " Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;
}
- showfield=field; // just to show one time
+ showfield=field; // just to show one time
} // while !eof()
- // generation de la geometrie
- // 1 construction des aretes
- // construire des aretes en chaque sommets
-
+ // generation de la geometrie
+ // 1 construction des aretes
+ // construire des aretes en chaque sommets
+
if (nbv <=0) {
cerr<<"Error: the field Vertex is not found in " << filename<<endl;
NbErr++;}
@@ -1220,8 +1220,8 @@ void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)
;NbErr++;}
if(NbErr) MeshError(1);
-
+
}
-} // end of namespace bamg
+} // end of namespace bamg