diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index 1facbedceaa576d26faf6b0b1b247aa31ad0af2f..4801a3d8520857cb72d3334e9c0386faeb96307f 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -187,11 +187,11 @@ class GSL_Matrix
private:
gsl_matrix_view view;
public:
- inline size_t size1() const { return data->size1; }
- inline size_t size2() const { return data->size2; }
+ inline int size1() const { return data->size1; }
+ inline int size2() const { return data->size2; }
gsl_matrix *data;
GSL_Matrix(gsl_matrix_view _data) : view(_data), data(&view.matrix) {}
- GSL_Matrix(size_t R, size_t C) { data = gsl_matrix_calloc(R, C); }
+ GSL_Matrix(int r, int c) { data = gsl_matrix_calloc(r, c); }
GSL_Matrix() : data(0) {}
GSL_Matrix(const GSL_Matrix &other) : data(0)
{
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index b71eb8622112c840962d82c57d348dfa48189bd7..34b31bfdc3187e1679a0191b395365613907da31 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.60 2008-02-21 07:48:49 geuzaine Exp $
+// $Id: GModel.cpp,v 1.61 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -438,7 +438,7 @@ MVertex *GModel::getMeshVertex(int num)
if(num < 0) return 0;
if(_vertexVectorCache.empty() && _vertexMapCache.empty())
buildMeshVertexCache();
- if(num < _vertexVectorCache.size())
+ if(num < (int)_vertexVectorCache.size())
return _vertexVectorCache[num];
else
return _vertexMapCache[num];
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index d40b8375cb67114c7d33a71c4a438d3e5c21bbc5..a9b20d7f8e521729767acccb051151f9142e99c2 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.35 2008-02-20 09:20:45 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.36 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -273,9 +273,9 @@ bool OCCFace::buildSTLTriangulation()
TopLoc_Location loc;
int p1, p2, p3;
Bnd_Box aBox;
- Standard_Boolean bWithShare;
+ Standard_Boolean bWithShare = Standard_False;
Standard_Real aDiscret, aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
- Standard_Real dX, dY, dZ, dMax, aCoeff, aAngle;
+ Standard_Real dX, dY, dZ, dMax, aCoeff;
BRepBndLib::Add(s, aBox);
aBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 115077a035fff900198d1b5d5b7afa16952dcc0f..df6c20cceabf2f5a483d11afd5db90ebd05826ea 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.117 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.118 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -877,9 +877,9 @@ inline double dist2 (const SPoint2 &p1,const SPoint2 &p2)
// }
-static void printMesh1d (int iEdge, int seam, std::vector<SPoint2> &m){
+void printMesh1d (int iEdge, int seam, std::vector<SPoint2> &m){
printf("Mesh1D for edge %d seam %d\n",iEdge,seam);
- for (int i=0;i<m.size();i++){
+ for (unsigned int i = 0; i < m.size(); i++){
printf("%12.5E %12.5E\n",m[i].x(),m[i].y());
}
}
@@ -957,7 +957,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
while (unordered.size())
{
- if (MYDEBUG)printf("unordered.size() = %d\n",unordered.size());
+ if (MYDEBUG)printf("unordered.size() = %d\n", (int)unordered.size());
std::list<GEdgeSigned>::iterator it = unordered.begin();
std::vector<SPoint2> coords;
@@ -1048,7 +1048,7 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
}
Finalize:
- if (MYDEBUG)printf("Finalize, found %d points\n",coords.size());
+ if (MYDEBUG) printf("Finalize, found %d points\n", (int)coords.size());
if (coords.size() == 0){
// It has not worked : either tolerance is wrong or the first seam edge
// has to be taken with the other parametric coordinates (because it is
@@ -1074,7 +1074,9 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
edgeLoop.push_back(found.ge->mesh_vertices[i]);
}
- if (MYDEBUG)printf("edge %d size %d size %d\n",found.ge->tag(), (int)edgeLoop.size(), (int)coords.size());
+ if (MYDEBUG)
+ printf("edge %d size %d size %d\n",
+ found.ge->tag(), (int)edgeLoop.size(), (int)coords.size());
std::vector<BDS_Point*> edgeLoop_BDS;
for (unsigned int i=0;i<edgeLoop.size();i++)
@@ -1109,7 +1111,10 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
m->add_geom (ge->tag(), ge->dim());
BDS_GeomEntity *g = m->get_geom(ge->tag(),ge->dim());
pp->g = g;
- if (MYDEBUG)printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n",count,pp->u,pp->v,param.x(),param.y(),pp->g->classif_tag,pp->g->classif_degree);
+ if (MYDEBUG)
+ printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n",
+ count,pp->u,pp->v,param.x(),param.y(),pp->g->classif_tag,
+ pp->g->classif_degree);
bbox += SPoint3(U,V,0);
edgeLoop_BDS.push_back(pp);
recover_map[pp] = here;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 3feb8b99def5c25a1e126dc87a7c90fb7344977e..a96a3c866062523303a163b0bd9cdc722fb90830 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceBDS.cpp,v 1.7 2008-02-21 12:11:12 geuzaine Exp $
+// $Id: meshGFaceBDS.cpp,v 1.8 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -217,9 +217,9 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
normal_triangle(p41, p42, p43, norm22);
double cosa; prosca(norm11, norm12, &cosa);
double cosb; prosca(norm21, norm22, &cosb);
- double smoothIndicator = cosb - cosa;
- bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
- bool smoothCouldSwap = !(cosb < cosa * .5);
+ // double smoothIndicator = cosb - cosa;
+ // bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
+ // bool smoothCouldSwap = !(cosb < cosa * .5);
double la = computeEdgeLinearLength(p11, p12);
double la_ = computeEdgeLinearLength(p11, p12, gf, m.scalingU, m.scalingV);
@@ -440,7 +440,7 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
std::sort(edges.begin(), edges.end());
- for (int i = 0; i < edges.size(); ++i){
+ for (unsigned int i = 0; i < edges.size(); ++i){
BDS_Edge *e = edges[i].second;
if (!e->deleted){
const double coord = 0.5;
@@ -477,7 +477,7 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
std::sort(edges.begin(), edges.end());
- for (int i = 0; i < edges.size(); i++){
+ for (unsigned int i = 0; i < edges.size(); i++){
BDS_Edge *e = edges[i].second;
if(!e->deleted){
bool res = false;
@@ -562,7 +562,6 @@ void gmshRefineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const double MINE_ = 0.67, MAXE_ = 1.4;
- double mesh_quality = 0;
while (1){
// we count the number of local mesh modifs.
int nb_split =0;
@@ -733,7 +732,6 @@ void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
gmshDelaunayizeBDS(gf, m, nb_swap);
for (int ITER = 0; ITER < 3; ITER++){
- double LIMIT = .1;
for (int KK = 0; KK < 4; KK++){
// swap edges that provide a better configuration
int NN1 = m.edges.size();
@@ -751,7 +749,6 @@ void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
}
}
- int nbSplit = 0;
if (recover_map){
while(gmshSolveInvalidPeriodic(gf, m, recover_map)){
}
@@ -762,12 +759,6 @@ void gmshOptimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
void delaunayPointInsertionBDS(GFace *gf, BDS_Mesh &m, BDS_Point *v, BDS_Face *f)
{
- const double p[2] = {v->u, v->v};
-
- BDS_Edge *e1 = f->e1;
- BDS_Edge *e2 = f->e2;
- BDS_Edge *e3 = f->e3;
- // face is splitted,
m.split_face(f, v);
int nb_swap = 0;
gmshDelaunayizeBDS(gf, m, nb_swap);
@@ -782,7 +773,7 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
GFace *gf = *it;
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g2 = m->get_geom(gf->tag(), 2);
- for (int i = 0; i < gf->triangles.size(); i++){
+ for (unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *e = gf->triangles[i];
BDS_Point *p[3];
for (int j = 0; j < 3; j++){
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 5327b53b24d3dce6570cc24c9222bca476264c84..657732f708247ed8dcbfdbe20e22b3278654c371 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.10 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.11 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -31,21 +31,19 @@
#include <map>
#include <algorithm>
-int III = 1;
-
-bool circumCenterMetricInTriangle ( MTriangle *base,
- const double *metric,
- const std::vector<double> & Us,
- const std::vector<double> & Vs)
+bool circumCenterMetricInTriangle(MTriangle *base,
+ const double *metric,
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs)
{
- double R,x[2],uv[2];
- circumCenterMetric(base,metric,Us,Vs,x,R);
- return invMapUV(base,x,Us,Vs,uv,1.e-8);
+ double R, x[2], uv[2];
+ circumCenterMetric(base, metric, Us, Vs, x, R);
+ return invMapUV(base, x, Us, Vs, uv, 1.e-8);
}
-void circumCenterMetric ( double *pa, double *pb, double *pc,
- const double *metric,
- double *x, double &Radius2)
+void circumCenterMetric(double *pa, double *pb, double *pc,
+ const double *metric,
+ double *x, double &Radius2)
{
// d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
double sys[2][2];
@@ -77,73 +75,57 @@ void circumCenterMetric ( double *pa, double *pb, double *pc,
2. * (x[0] - pa[0]) * (x[1] - pa[1]) * b;
}
-void circumCenterMetric ( MTriangle *base,
- const double *metric,
- const std::vector<double> & Us,
- const std::vector<double> & Vs,
- double *x, double &Radius2)
+void circumCenterMetric(MTriangle *base,
+ const double *metric,
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs,
+ double *x, double &Radius2)
{
// d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
-
-
- double sys[2][2];
- double rhs[2];
double pa[2] = {Us[base->getVertex(0)->getNum()],
Vs[base->getVertex(0)->getNum()]};
double pb[2] = {Us[base->getVertex(1)->getNum()],
Vs[base->getVertex(1)->getNum()]};
double pc[2] = {Us[base->getVertex(2)->getNum()],
Vs[base->getVertex(2)->getNum()]};
- circumCenterMetric ( pa,pb,pc,metric,x,Radius2);
+ circumCenterMetric(pa, pb, pc, metric, x, Radius2);
}
-
-void buildMetric ( GFace *gf , double *uv, double *metric){
- Pair<SVector3,SVector3> der = gf->firstDer(SPoint2(uv[0],uv[1]));
- metric[0] = dot(der.first(),der.first());
- metric[1] = dot(der.second(),der.first());
- metric[2] = dot(der.second(),der.second());
+void buildMetric(GFace *gf, double *uv, double *metric)
+{
+ Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
+ metric[0] = dot(der.first(), der.first());
+ metric[1] = dot(der.second(), der.first());
+ metric[2] = dot(der.second(), der.second());
}
-int inCircumCircleAniso ( GFace *gf, double *p1, double *p2, double *p3, double *uv, double *metric)
+int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
+ double *uv, double *metric)
{
- double x[2],Radius2;
- circumCenterMetric ( p1,p2,p3,metric, x,Radius2);
+ double x[2], Radius2;
+ circumCenterMetric(p1, p2, p3, metric, x, Radius2);
const double a = metric[0];
const double b = metric[1];
const double d = metric[2];
double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
+ (x[1] - uv[1]) * (x[1] - uv[1]) * d
+ 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
- // printf("%22.15E %22.15E %d\n",d2,Radius2,d2 < Radius2);
-
return d2 < Radius2;
}
-
-int inCircumCircleAniso ( GFace *gf,
- MTriangle *base,
- const double *uv ,
- const double *metricb,
- const std::vector<double> & Us,
- const std::vector<double> & Vs)
+int inCircumCircleAniso(GFace *gf, MTriangle *base,
+ const double *uv, const double *metricb,
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs)
{
-
- double x[2],Radius2,metric[3];
-
-// double pa[2] = {(Us[base->getVertex(0)->getNum()]+Us[base->getVertex(1)->getNum()]+Us[base->getVertex(2)->getNum()])/3.,
-// (Vs[base->getVertex(0)->getNum()]+Vs[base->getVertex(1)->getNum()]+Vs[base->getVertex(2)->getNum()])/3.};
- double pa[2] = {(Us[base->getVertex(0)->getNum()]+Us[base->getVertex(1)->getNum()]+Us[base->getVertex(2)->getNum()])/3.,
- (Vs[base->getVertex(0)->getNum()]+Vs[base->getVertex(1)->getNum()]+Vs[base->getVertex(2)->getNum()])/3.};
+ double x[2], Radius2, metric[3];
+ double pa[2] = {(Us[base->getVertex(0)->getNum()] + Us[base->getVertex(1)->getNum()] +
+ Us[base->getVertex(2)->getNum()]) / 3.,
+ (Vs[base->getVertex(0)->getNum()] + Vs[base->getVertex(1)->getNum()] +
+ Vs[base->getVertex(2)->getNum()]) / 3.};
- buildMetric (gf,pa,metric);
- circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
-// buildMetric (gf,pa,metric);
-// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
-// buildMetric (gf,pa,metric);
-// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
-// buildMetric (gf,pa,metric);
-// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
+ buildMetric(gf, pa, metric);
+ circumCenterMetric(base, metric, Us, Vs, x, Radius2);
const double a = metric[0];
const double b = metric[1];
@@ -153,46 +135,43 @@ int inCircumCircleAniso ( GFace *gf,
+ (x[1] - uv[1]) * (x[1] - uv[1]) * d
+ 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
- // printf("%22.15E %22.15E %d\n",d2,Radius2,d2 < Radius2);
return d2 < Radius2;
}
-MTri3::MTri3 ( MTriangle * t, double lc) : deleted(false), base (t)
+MTri3::MTri3(MTriangle *t, double lc) : deleted(false), base(t)
{
neigh[0] = neigh[1] = neigh[2] = 0;
// compute the circumradius of the triangle
- double pa[3] = {base->getVertex(0)->x(),base->getVertex(0)->y(),base->getVertex(0)->z()};
- double pb[3] = {base->getVertex(1)->x(),base->getVertex(1)->y(),base->getVertex(1)->z()};
- double pc[3] = {base->getVertex(2)->x(),base->getVertex(2)->y(),base->getVertex(2)->z()};
+ double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(), base->getVertex(0)->z()};
+ double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
+ double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
double center[3];
- base->circumcenterXYZ(pa,pb,pc,center);
+ base->circumcenterXYZ(pa, pb, pc, center);
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
- circum_radius = sqrt ( dx*dx + dy*dy + dz*dz);
+ circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
circum_radius /= lc;
}
-
-int MTri3::inCircumCircle ( const double *p ) const
+int MTri3::inCircumCircle(const double *p) const
{
- double pa[3] = {base->getVertex(0)->x(),base->getVertex(0)->y(),base->getVertex(0)->z()};
- double pb[3] = {base->getVertex(1)->x(),base->getVertex(1)->y(),base->getVertex(1)->z()};
- double pc[3] = {base->getVertex(2)->x(),base->getVertex(2)->y(),base->getVertex(2)->z()};
+ double pa[3] = {base->getVertex(0)->x(), base->getVertex(0)->y(), base->getVertex(0)->z()};
+ double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
+ double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
double fourth[3];
- fourthPoint ( pa,pb,pc,fourth);
+ fourthPoint(pa, pb, pc, fourth);
- double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) * gmsh::orient3d(pa, pb, pc,fourth);
+ double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) *
+ gmsh::orient3d(pa, pb, pc,fourth);
return (result > 0) ? 1 : 0;
}
-int inCircumCircle ( MTriangle *base,
- const double *p ,
- const double *param ,
- std::vector<double> & Us,
- std::vector<double> & Vs)
+int inCircumCircle(MTriangle *base,
+ const double *p, const double *param ,
+ std::vector<double> &Us, std::vector<double> &Vs)
{
double pa[2] = {Us[base->getVertex(0)->getNum()],
Vs[base->getVertex(0)->getNum()]};
@@ -203,147 +182,110 @@ int inCircumCircle ( MTriangle *base,
double result = gmsh::incircle(pa, pb, pc, (double*)param) * gmsh::orient2d(pa, pb, pc);
return (result > 0) ? 1 : 0;
-
-
-
-// double pa[3] = {base->getVertex(0)->x(),base->getVertex(0)->y(),base->getVertex(0)->z()};
-// double pb[3] = {base->getVertex(1)->x(),base->getVertex(1)->y(),base->getVertex(1)->z()};
-// double pc[3] = {base->getVertex(2)->x(),base->getVertex(2)->y(),base->getVertex(2)->z()};
-// double fourth[3];
-// fourthPoint ( pa,pb,pc,fourth);
-
-// double result = gmsh::insphere(pa, pb, pc, fourth, (double*)p) * gmsh::orient3d(pa, pb, pc,fourth);
-// return (result > 0) ? 1 : 0;
}
-
template <class ITER>
-void connectTris ( ITER beg, ITER end)
+void connectTris(ITER beg, ITER end)
{
std::set<edgeXface> conn;
- {
- while (beg != end)
- {
- if (!(*beg)->isDeleted())
- {
- for (int i=0;i<3;i++)
- {
- edgeXface fxt ( *beg , i );
- std::set<edgeXface>::iterator found = conn.find (fxt);
- if (found == conn.end())
- conn.insert ( fxt );
- else if (found->t1 != *beg)
- {
- found->t1->setNeigh(found->i1,*beg);
- (*beg)->setNeigh ( i, found->t1);
- }
- }
- }
- ++beg;
+ while (beg != end){
+ if (!(*beg)->isDeleted()){
+ for (int i = 0; i < 3; i++){
+ edgeXface fxt(*beg, i);
+ std::set<edgeXface>::iterator found = conn.find(fxt);
+ if (found == conn.end())
+ conn.insert(fxt);
+ else if (found->t1 != *beg){
+ found->t1->setNeigh(found->i1, *beg);
+ (*beg)->setNeigh(i, found->t1);
+ }
}
+ }
+ ++beg;
}
}
-void connectTriangles ( std::list<MTri3*> & l)
+void connectTriangles(std::list<MTri3*> &l)
{
- connectTris(l.begin(),l.end());
+ connectTris(l.begin(), l.end());
}
-void connectTriangles ( std::vector<MTri3*> & l)
+
+void connectTriangles(std::vector<MTri3*> &l)
{
- connectTris(l.begin(),l.end());
+ connectTris(l.begin(), l.end());
}
-void connectTriangles ( std::set<MTri3*,compareTri3Ptr> & l){
- connectTris(l.begin(),l.end());
+
+void connectTriangles(std::set<MTri3*, compareTri3Ptr> &l)
+{
+ connectTris(l.begin(), l.end());
}
-void recurFindCavity (std::list<edgeXface> & shell,
- std::list<MTri3*> & cavity,
- double *v,
- double *param,
- MTri3 *t,
- std::vector<double> & Us,
- std::vector<double> & Vs)
+void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
+ double *v, double *param, MTri3 *t,
+ std::vector<double> &Us, std::vector<double> &Vs)
{
t->setDeleted(true);
// the cavity that has to be removed
// because it violates delaunay criterion
cavity.push_back(t);
- for (int i=0;i<3;i++)
- {
- MTri3 *neigh = t->getNeigh(i) ;
- if (!neigh)
- shell.push_back ( edgeXface ( t, i ) );
- else if (!neigh->isDeleted())
- {
- int circ = inCircumCircle (neigh->tri(), v , param, Us, Vs);
- if (circ)
- recurFindCavity ( shell, cavity,v, param, neigh,Us,Vs);
- else
- shell.push_back ( edgeXface ( t, i ) );
- }
+ for (int i = 0; i < 3; i++){
+ MTri3 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back(edgeXface(t, i));
+ else if (!neigh->isDeleted()){
+ int circ = inCircumCircle(neigh->tri(), v , param, Us, Vs);
+ if (circ)
+ recurFindCavity(shell, cavity, v, param, neigh, Us, Vs);
+ else
+ shell.push_back(edgeXface(t, i));
}
+ }
}
void recurFindCavityAniso (GFace *gf,
- std::list<edgeXface> & shell,
- std::list<MTri3*> & cavity,
- double *metric,
- double *param,
- MTri3 *t,
- std::vector<double> & Us,
- std::vector<double> & Vs)
+ std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
+ double *metric, double *param, MTri3 *t,
+ std::vector<double> &Us, std::vector<double> &Vs)
{
t->setDeleted(true);
// the cavity that has to be removed
// because it violates delaunay criterion
cavity.push_back(t);
- for (int i=0;i<3;i++)
- {
- MTri3 *neigh = t->getNeigh(i) ;
- if (!neigh)
- shell.push_back ( edgeXface ( t, i ) );
- else if (!neigh->isDeleted())
- {
- int circ = inCircumCircleAniso (gf, neigh->tri(), param, metric, Us, Vs);
- if (circ)
- recurFindCavityAniso ( gf, shell, cavity,metric, param, neigh,Us,Vs);
- else
- shell.push_back ( edgeXface ( t, i ) );
- }
+ for (int i = 0; i < 3; i++){
+ MTri3 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back(edgeXface(t, i));
+ else if (!neigh->isDeleted()){
+ int circ = inCircumCircleAniso(gf, neigh->tri(), param, metric, Us, Vs);
+ if (circ)
+ recurFindCavityAniso(gf, shell, cavity,metric, param, neigh, Us, Vs);
+ else
+ shell.push_back(edgeXface(t, i));
}
+ }
}
-
-
-bool circUV ( MTriangle *t ,
- std::vector<double> & Us,
- std::vector<double> & Vs , double *res, GFace *gf)
+bool circUV(MTriangle *t, std::vector<double> & Us, std::vector<double> &Vs,
+ double *res, GFace *gf)
{
- double u1 [3] = {Us[t->getVertex(0)->getNum()],Vs[t->getVertex(0)->getNum()],0};
- double u2 [3] = {Us[t->getVertex(1)->getNum()],Vs[t->getVertex(1)->getNum()],0};
- double u3 [3] = {Us[t->getVertex(2)->getNum()],Vs[t->getVertex(2)->getNum()],0};
- t->circumcenterXY ( u1, u2 , u3, res);
+ double u1 [3] = {Us[t->getVertex(0)->getNum()], Vs[t->getVertex(0)->getNum()], 0};
+ double u2 [3] = {Us[t->getVertex(1)->getNum()], Vs[t->getVertex(1)->getNum()], 0};
+ double u3 [3] = {Us[t->getVertex(2)->getNum()], Vs[t->getVertex(2)->getNum()], 0};
+ t->circumcenterXY(u1, u2, u3, res);
return true;
- double p1 [3] = {t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z()};
- double p2 [3] = {t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z()};
- double p3 [3] = {t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z()};
+ double p1 [3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
+ double p2 [3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
+ double p3 [3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
double resxy[3], uv[2];
- t->circumcenterXYZ ( p1, p2 , p3, resxy,uv);
+ t->circumcenterXYZ(p1, p2, p3, resxy,uv);
return true;
- // printf("uv = %g %g\n",uv[0],uv[1]);
- // printf("%g %g vs ",res[0],res[1]);
- res[0] = u1[0] + uv[0] * ( u2[0] - u1[0] ) + uv[1] * ( u3[0] - u1[0] );
- res[1] = u1[1] + uv[0] * ( u2[1] - u1[1] ) + uv[1] * ( u3[1] - u1[1] );
- // printf("%g %g \n",res[0],res[1]);
-
}
-bool invMapUV ( MTriangle *t ,
- double *p,
- const std::vector<double> & Us,
- const std::vector<double> & Vs , double *uv, double tol)
+bool invMapUV(MTriangle *t, double *p,
+ const std::vector<double> &Us, const std::vector<double> &Vs,
+ double *uv, double tol)
{
double mat[2][2];
double b[2];
@@ -374,10 +316,7 @@ bool invMapUV ( MTriangle *t ,
return false;
}
-
-double getSurfUV ( MTriangle *t ,
- std::vector<double> & Us,
- std::vector<double> & Vs )
+double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
{
double u1 = Us[t->getVertex(0)->getNum()];
double v1 = Vs[t->getVertex(0)->getNum()];
@@ -389,279 +328,190 @@ double getSurfUV ( MTriangle *t ,
const double vv2 [2] = {u3 - u1, v3 - v1};
double s = vv1[0] * vv2[1] - vv1[1] * vv2[0];
return s * 0.5;
-
}
-bool insertVertex (GFace *gf,
- MVertex *v ,
- double *param ,
- MTri3 *t ,
- std::set<MTri3*,compareTri3Ptr> &allTets,
- std::vector<double> & vSizes,
- std::vector<double> & vSizesBGM,
- std::vector<double> & Us,
- std::vector<double> & Vs,
- double *metric = 0)
+bool insertVertex(GFace *gf, MVertex *v, double *param , MTri3 *t,
+ std::set<MTri3*, compareTri3Ptr> &allTets,
+ std::vector<double> &vSizes, std::vector<double> &vSizesBGM,
+ std::vector<double> &Us, std::vector<double> &Vs,
+ double *metric = 0)
{
- std::list<edgeXface> shell;
- std::list<MTri3*> cavity;
- std::list<MTri3*> new_cavity;
+ std::list<edgeXface> shell;
+ std::list<MTri3*> cavity;
+ std::list<MTri3*> new_cavity;
if (!metric){
- double p[3] = {v->x(),v->y(),v->z()};
- recurFindCavity ( shell, cavity, p , param, t, Us, Vs);
+ double p[3] = {v->x(), v->y(), v->z()};
+ recurFindCavity(shell, cavity, p, param, t, Us, Vs);
}
else{
- recurFindCavityAniso ( gf, shell, cavity, metric , param, t, Us, Vs);
+ recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
}
// check that volume is conserved
double newVolume = 0;
double oldVolume = 0;
-
std::list<MTri3*>::iterator ittet = cavity.begin();
std::list<MTri3*>::iterator ittete = cavity.end();
- while ( ittet != ittete )
- {
- oldVolume += fabs(getSurfUV((*ittet)->tri(),Us,Vs));
- ++ittet;
- }
-
-// fprintf(ff,"};\n");
-// fclose (ff);
-// char name[256];
-// sprintf(name,"test.pos");
-// FILE *ff = fopen (name,"w");
-// fprintf(ff,"View\"test\"{\n");
-
-
- MTri3** newTris = new MTri3*[shell.size()];;
+ while(ittet != ittete){
+ oldVolume += fabs(getSurfUV((*ittet)->tri(),Us,Vs));
+ ++ittet;
+ }
+
+ MTri3 **newTris = new MTri3*[shell.size()];
int k = 0;
- std::list<edgeXface>::iterator it = shell.begin();
-
- while (it != shell.end())
- {
- MTriangle *t = new MTriangle ( it->v[0],
- it->v[1],
- v);
-
- double lc = 0.3333333333*(vSizes [t->getVertex(0)->getNum()]+
- vSizes [t->getVertex(1)->getNum()]+
- vSizes [t->getVertex(2)->getNum()]);
- double lcBGM = 0.3333333333*(vSizesBGM [t->getVertex(0)->getNum()]+
+ std::list<edgeXface>::iterator it = shell.begin();
+
+ while (it != shell.end()){
+ MTriangle *t = new MTriangle(it->v[0], it->v[1], v);
+ double lc = 0.3333333333 * (vSizes [t->getVertex(0)->getNum()]+
+ vSizes [t->getVertex(1)->getNum()]+
+ vSizes [t->getVertex(2)->getNum()]);
+ double lcBGM = 0.3333333333 * (vSizesBGM [t->getVertex(0)->getNum()]+
vSizesBGM [t->getVertex(1)->getNum()]+
vSizesBGM [t->getVertex(2)->getNum()]);
-
- MTri3 *t4 = new MTri3 ( t , Extend1dMeshIn2dSurfaces() ? std::min(lc,lcBGM) : lcBGM );
-// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
-// Us [t->getVertex(0)->getNum()],
-// Vs [t->getVertex(0)->getNum()],
-// 0.0,
-// Us [t->getVertex(1)->getNum()],
-// Vs [t->getVertex(1)->getNum()],
-// 0.0,
-// Us [t->getVertex(2)->getNum()],
-// Vs [t->getVertex(2)->getNum()],
-// 0.0);
-// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
-// t->getVertex(0)->x(),
-// t->getVertex(0)->y(),
-// t->getVertex(0)->z(),
-// t->getVertex(1)->x(),
-// t->getVertex(1)->y(),
-// t->getVertex(1)->z(),
-// t->getVertex(2)->x(),
-// t->getVertex(2)->y(),
-// t->getVertex(2)->z());
- newTris[k++]=t4;
- // all new triangles are pushed front in order to
- // ba able to destroy them if the cavity is not
- // star shaped around the new vertex.
- new_cavity.push_back (t4);
- MTri3 *otherSide = it->t1->getNeigh(it->i1);
-
- if (otherSide)
- new_cavity.push_back (otherSide);
- // if (!it->t1->isDeleted())throw;
- double ss = fabs(getSurfUV(t4->tri(),Us,Vs));
- if (ss < 1.e-25)ss = 1256172121;
- newVolume += ss;
- ++it;
- }
-// fprintf(ff,"};\n");
-// fclose (ff);
-// printf("%d %d newVolume %g oldVolume %g\n",cavity.size(),new_cavity.size(),newVolume,oldVolume);
-// getchar();
-
-
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume )
- {
- connectTris ( new_cavity.begin(),new_cavity.end() );
- allTets.insert(newTris,newTris+shell.size());
- delete [] newTris;
- return true;
- }
+ MTri3 *t4 = new MTri3(t, Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM);
+ newTris[k++] = t4;
+ // all new triangles are pushed front in order to
+ // ba able to destroy them if the cavity is not
+ // star shaped around the new vertex.
+ new_cavity.push_back (t4);
+ MTri3 *otherSide = it->t1->getNeigh(it->i1);
+ if (otherSide)
+ new_cavity.push_back (otherSide);
+ double ss = fabs(getSurfUV(t4->tri(),Us,Vs));
+ if (ss < 1.e-25)ss = 1256172121;
+ newVolume += ss;
+ ++it;
+ }
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume){
+ connectTris(new_cavity.begin(),new_cavity.end());
+ allTets.insert(newTris,newTris+shell.size());
+ delete [] newTris;
+ return true;
+ }
// The cavity is NOT star shaped
- else
- {
- // printf("the cavity (%d members) is not star shaped %g vs %g\n",shell.size(),oldVolume,newVolume);
-
- for (unsigned int i=0;i<shell.size();i++)delete newTris[i];
- delete [] newTris;
- ittet = cavity.begin();
- ittete = cavity.end();
- while ( ittet != ittete )
- {
- (*ittet)->setDeleted(false);
- ++ittet;
- }
- return false;
+ else{
+ for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
+ delete [] newTris;
+ ittet = cavity.begin();
+ ittete = cavity.end();
+ while(ittet != ittete){
+ (*ittet)->setDeleted(false);
+ ++ittet;
}
+ return false;
+ }
}
-void _printTris (char *name, std::set<MTri3*,compareTri3Ptr>&AllTris, std::vector<double>&Us, std::vector<double>&Vs)
+void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
+ std::vector<double> &Us, std::vector<double> &Vs)
{
FILE *ff = fopen (name,"w");
fprintf(ff,"View\"test\"{\n");
std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
- while (it != AllTris.end() )
- {
- MTri3 *worst = *it;
- if (!worst->isDeleted())
- {
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
- Us [(worst)->tri()->getVertex(0)->getNum()],
- Vs [(worst)->tri()->getVertex(0)->getNum()],
- 0.0,
- Us [(worst)->tri()->getVertex(1)->getNum()],
- Vs [(worst)->tri()->getVertex(1)->getNum()],
- 0.0,
- Us [(worst)->tri()->getVertex(2)->getNum()],
- Vs [(worst)->tri()->getVertex(2)->getNum()],
- 0.0);
- }
- ++it;
+ while (it != AllTris.end() ){
+ MTri3 *worst = *it;
+ if (!worst->isDeleted()){
+ fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
+ Us [(worst)->tri()->getVertex(0)->getNum()],
+ Vs [(worst)->tri()->getVertex(0)->getNum()],
+ 0.0,
+ Us [(worst)->tri()->getVertex(1)->getNum()],
+ Vs [(worst)->tri()->getVertex(1)->getNum()],
+ 0.0,
+ Us [(worst)->tri()->getVertex(2)->getNum()],
+ Vs [(worst)->tri()->getVertex(2)->getNum()],
+ 0.0);
}
+ ++it;
+ }
fprintf(ff,"};\n");
fclose (ff);
}
-void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
+void insertVerticesInFace(GFace *gf, BDS_Mesh *bds)
{
-
std::set<MTri3*,compareTri3Ptr> AllTris;
- std::vector<double> vSizes, vSizesBGM,Us,Vs;
-
+ std::vector<double> vSizes, vSizesBGM, Us, Vs;
- buidMeshGenerationDataStructures (gf,AllTris,vSizes,vSizesBGM,Us,Vs);
+ buidMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
- // _printTris ("before.pos", AllTris, Us,Vs);
- int nbSwaps = edgeSwapPass(gf,AllTris,SWCR_DEL,Us,Vs,vSizes,vSizesBGM);
- // _printTris ("after2.pos", AllTris, Us,Vs);
- Msg(DEBUG,"Delaunization of the initial mesh done (%d swaps)",nbSwaps);
+ // _printTris ("before.pos", AllTris, Us,Vs);
+ int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
+ // _printTris ("after2.pos", AllTris, Us,Vs);
+ Msg(DEBUG,"Delaunization of the initial mesh done (%d swaps)", nbSwaps);
int ITER = 0;
-
- while (1)
- {
- MTri3 *worst = *AllTris.begin();
-
- if (worst->isDeleted()){
- delete worst->tri();
- delete worst;
- AllTris.erase(AllTris.begin());
- // Msg(INFO,"Worst tet is deleted");
- }
- else{
- if(ITER++%5000 ==0)
- Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
- double center[2],uv[2],metric[3],r2;
- if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
-
- circUV(worst->tri(),Us,Vs,center,gf);
- MTriangle *base = worst->tri();
- double pa[2] = {(Us[base->getVertex(0)->getNum()]+Us[base->getVertex(1)->getNum()]+Us[base->getVertex(2)->getNum()])/3.,
- (Vs[base->getVertex(0)->getNum()]+Vs[base->getVertex(1)->getNum()]+Vs[base->getVertex(2)->getNum()])/3.};
- buildMetric ( gf , pa, metric);
- circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
-// buildMetric ( gf , center, metric);
-// circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
-// buildMetric ( gf , center, metric);
-// circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
-
- // printf("%g %g %g\n",metric[0],metric[1],metric[2]);
-
- bool inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
- // if (!inside)
- // circUV(worst->tri(),Us,Vs,center,gf);
-
- // inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
+ while (1){
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
+ else{
+ if(ITER++ % 5000 == 0)
+ Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",
+ vSizes.size(), worst->getRadius());
+ double center[2],uv[2],metric[3],r2;
+ if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
+ circUV(worst->tri(), Us, Vs, center, gf);
+ MTriangle *base = worst->tri();
+ double pa[2] = {(Us[base->getVertex(0)->getNum()] +
+ Us[base->getVertex(1)->getNum()] +
+ Us[base->getVertex(2)->getNum()]) / 3.,
+ (Vs[base->getVertex(0)->getNum()] +
+ Vs[base->getVertex(1)->getNum()] +
+ Vs[base->getVertex(2)->getNum()]) / 3.};
+ buildMetric(gf, pa, metric);
+ circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+ bool inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
+ if (inside) {
+ // we use here local coordinates as real coordinates
+ // x,y and z will be computed hereafter
+ // Msg(INFO,"Point is inside");
+ GPoint p = gf->point(center[0], center[1]);
+ MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
+ v->setNum(Us.size());
+ double lc1 = ((1. - uv[0] - uv[1]) * vSizes[worst->tri()->getVertex(0)->getNum()] +
+ uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] +
+ uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()]);
+ // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
+ double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
+ vSizesBGM.push_back(lc);
+ vSizes.push_back(lc1);
+ Us.push_back(center[0]);
+ Vs.push_back(center[1]);
- if (inside) {
- // we use here local coordinates as real coordinates
- // x,y and z will be computed hereafter
- // Msg(INFO,"Point is inside");
- GPoint p = gf->point (center[0],center[1]);
- MVertex *v = new MFaceVertex (p.x(),p.y(),p.z(),gf,center[0],center[1]);
- v->setNum(Us.size());
- double lc1 = ((1.-uv[0]-uv[1]) * vSizes [worst->tri()->getVertex(0)->getNum()] +
- uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] +
- uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()] );
- // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
- double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
- // printf("lc1 %12.5E lc %12.5E\n",lc1,lc);
-
- vSizesBGM.push_back( lc );
- vSizes.push_back ( lc1);
- Us.push_back( center[0] );
- Vs.push_back( center[1] );
-
- if (!insertVertex ( gf, v , center, worst, AllTris,vSizes,vSizesBGM,Us,Vs,metric)) {
- Msg(DEBUG,"2D Delaunay : a cavity is not star shaped");
- AllTris.erase(AllTris.begin());
- worst->forceRadius(-1);
- AllTris.insert(worst);
- delete v;
- }
- else
- gf->mesh_vertices.push_back(v);
- }
- else {
- // printf("%g %g %g\n",metric[0],metric[1],metric[2]);
- // Msg(DEBUG,"Point %g %g is outside %g %g - %g %g - %g %g (%22.15E %22.15E)",
- // center[0],center[1],
- // Us [(worst)->tri()->getVertex(0)->getNum()],
- // Vs [(worst)->tri()->getVertex(0)->getNum()],
- // Us [(worst)->tri()->getVertex(1)->getNum()],
- // Vs [(worst)->tri()->getVertex(1)->getNum()],
- // Us [(worst)->tri()->getVertex(2)->getNum()],
- // Vs [(worst)->tri()->getVertex(2)->getNum()],
- // uv[0],uv[1]);
- // throw;
+ if (!insertVertex(gf, v, center, worst, AllTris, vSizes, vSizesBGM,
+ Us, Vs, metric)) {
+ Msg(DEBUG,"2D Delaunay : a cavity is not star shaped");
AllTris.erase(AllTris.begin());
- worst->forceRadius(0);
- AllTris.insert(worst);
- // break;
+ worst->forceRadius(-1);
+ AllTris.insert(worst);
+ delete v;
}
+ else
+ gf->mesh_vertices.push_back(v);
}
- }
+ else {
+ AllTris.erase(AllTris.begin());
+ worst->forceRadius(0);
+ AllTris.insert(worst);
+ }
+ }
+ }
for (int i=0;i<10;i++){
- // bool splits = edgeSplitPass(1.4,gf,AllTris,SPCR_ALLWAYS,Us,Vs,vSizes,vSizesBGM);
- // edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM);
- // bool collapses = edgeCollapsePass(0.67,gf,AllTris,Us,Vs,vSizes,vSizesBGM);
- // edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM);
- // if (!collapses && !splits)break;
+ // bool splits = edgeSplitPass(1.4, gf, AllTris, SPCR_ALLWAYS, Us, Vs, vSizes, vSizesBGM);
+ // edgeSwapPass(gf, AllTris, SWCR_QUAL, Us, Vs, vSizes, vSizesBGM);
+ // bool collapses = edgeCollapsePass(0.67, gf, AllTris, Us, Vs, vSizes, vSizesBGM);
+ // edgeSwapPass(gf, AllTris, SWCR_QUAL, Us, Vs, vSizes, vSizesBGM);
+ // if (!collapses && !splits) break;
}
- // for (int i=0;i<1;i++)if(!edgeCollapsePass(100,gf,AllTris,Us,Vs,vSizes,vSizesBGM))break;
- // for (int i=0;i<100;i++)if(!edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM))break;
-
- // char name[245];
- // sprintf(name,"paramFinal%d.pos",gf->tag());
- // _printTris (name, AllTris, Us,Vs);
- _printTris ("yo.pos", AllTris, Us,Vs);
+ //_printTris ("yo.pos", AllTris, Us, Vs);
- transferDataStructure (gf,AllTris);
- // fill new gmsh structures with triangles
+ transferDataStructure(gf, AllTris);
}
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index eae9ef8a427ea23a68e24279dc086ba16f354cda..3160b32f2bbe05ad01f80c6748b48e9e8750b1df 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceOptimize.cpp,v 1.10 2008-02-17 08:48:01 geuzaine Exp $
+// $Id: meshGFaceOptimize.cpp,v 1.11 2008-02-21 13:34:40 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -28,76 +28,67 @@
#include "MElement.h"
#include "BackgroundMesh.h"
-static void setLcsMax ( MTriangle *t, std::map<MVertex*,double> &vSizes)
+static void setLcsMax(MTriangle *t, std::map<MVertex*, double> &vSizes)
{
- for (int i=0;i<3;i++)
- {
- for (int j=i+1;j<3;j++)
- {
- MVertex *vi = t->getVertex(i);
- MVertex *vj = t->getVertex(j);
- vSizes[vi] = 1.e12;
- vSizes[vj] = 1.e12;
- }
+ for (int i = 0; i < 3; i++){
+ for (int j = i + 1; j < 3; j++){
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ vSizes[vi] = 1.e12;
+ vSizes[vj] = 1.e12;
}
+ }
}
-
-static void setLcs ( MTriangle *t, std::map<MVertex*,double> &vSizes)
+static void setLcs(MTriangle *t, std::map<MVertex*,double> &vSizes)
{
- for (int i=0;i<3;i++)
- {
- for (int j=i+1;j<3;j++)
- {
- MVertex *vi = t->getVertex(i);
- MVertex *vj = t->getVertex(j);
-
- double dx = vi->x()-vj->x();
- double dy = vi->y()-vj->y();
- double dz = vi->z()-vj->z();
- double l = sqrt(dx*dx+dy*dy+dz*dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
- if (iti->second > l)iti->second = l;
- if (itj->second > l)itj->second = l;
- }
+ for (int i = 0; i < 3; i++){
+ for (int j = i + 1; j < 3; j++){
+ MVertex *vi = t->getVertex(i);
+ MVertex *vj = t->getVertex(j);
+ double dx = vi->x()-vj->x();
+ double dy = vi->y()-vj->y();
+ double dz = vi->z()-vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if (iti->second > l) iti->second = l;
+ if (itj->second > l) itj->second = l;
}
+ }
}
-
-void buidMeshGenerationDataStructures (GFace *gf, std::set<MTri3*,compareTri3Ptr> &AllTris,
- std::vector<double> & vSizes,
- std::vector<double> & vSizesBGM,
- std::vector<double> & Us,
- std::vector<double> & Vs ){
-
- std::map<MVertex*,double> vSizesMap;
- for (unsigned int i=0;i<gf->triangles.size();i++)setLcsMax ( gf->triangles[i] , vSizesMap);
- for (unsigned int i=0;i<gf->triangles.size();i++)setLcs ( gf->triangles[i] , vSizesMap);
- int NUM=0;
- for (std::map<MVertex*,double>::iterator it = vSizesMap.begin();it!=vSizesMap.end();++it)
- {
- it->first->setNum(NUM++);
- vSizes.push_back(it->second);
- vSizesBGM.push_back(it->second);
- double u0,v0;
- parametricCoordinates ( it->first, gf, u0, v0);
- Us.push_back(u0);
- Vs.push_back(v0);
- }
- for (unsigned int i=0;i<gf->triangles.size();i++)
- {
- double lc = 0.3333333333*(vSizes [gf->triangles[i]->getVertex(0)->getNum()]+
- vSizes [gf->triangles[i]->getVertex(1)->getNum()]+
- vSizes [gf->triangles[i]->getVertex(2)->getNum()]);
- AllTris.insert ( new MTri3 ( gf->triangles[i] ,lc ) );
- }
+void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM,
+ std::vector<double> &Us,
+ std::vector<double> &Vs)
+{
+ std::map<MVertex*, double> vSizesMap;
+ for (unsigned int i = 0;i < gf->triangles.size(); i++) setLcsMax(gf->triangles[i], vSizesMap);
+ for (unsigned int i = 0;i < gf->triangles.size(); i++) setLcs(gf->triangles[i], vSizesMap);
+ int NUM = 0;
+ for (std::map<MVertex*, double>::iterator it = vSizesMap.begin(); it != vSizesMap.end(); ++it){
+ it->first->setNum(NUM++);
+ vSizes.push_back(it->second);
+ vSizesBGM.push_back(it->second);
+ double u0, v0;
+ parametricCoordinates(it->first, gf, u0, v0);
+ Us.push_back(u0);
+ Vs.push_back(v0);
+ }
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ double lc = 0.3333333333 * (vSizes [gf->triangles[i]->getVertex(0)->getNum()] +
+ vSizes [gf->triangles[i]->getVertex(1)->getNum()] +
+ vSizes [gf->triangles[i]->getVertex(2)->getNum()]);
+ AllTris.insert(new MTri3(gf->triangles[i], lc));
+ }
gf->triangles.clear();
- connectTriangles ( AllTris );
- // Msg(DEBUG,"All %d tris were connected",AllTris.size());
+ connectTriangles(AllTris );
}
-void transferDataStructure (GFace *gf,std::set<MTri3*,compareTri3Ptr> &AllTris){
+void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris)
+{
while (1) {
if (AllTris.begin() == AllTris.end() ) break;
MTri3 *worst = *AllTris.begin();
@@ -110,173 +101,139 @@ void transferDataStructure (GFace *gf,std::set<MTri3*,compareTri3Ptr> &AllTris){
}
}
-
-void buildVertexToTriangle ( std::vector<MTriangle*> &triangles, v2t_cont &adj )
+void buildVertexToTriangle(std::vector<MTriangle*> &triangles, v2t_cont &adj)
{
adj.clear();
- for (unsigned int i=0;i<triangles.size();i++)
- {
- MTriangle *t = triangles[i];
- for (unsigned int j=0;j<3;j++)
- {
- MVertex *v = t->getVertex(j);
- v2t_cont :: iterator it = adj.find ( v );
- if (it == adj.end())
- {
- std::vector<MTriangle*> one;
- one.push_back(t);
- adj[v] = one;
- }
- else
- {
- it->second.push_back(t);
- }
- }
+ for (unsigned int i = 0; i < triangles.size(); i++){
+ MTriangle *t = triangles[i];
+ for (unsigned int j = 0; j < 3; j++){
+ MVertex *v = t->getVertex(j);
+ v2t_cont :: iterator it = adj.find(v);
+ if (it == adj.end()){
+ std::vector<MTriangle*> one;
+ one.push_back(t);
+ adj[v] = one;
+ }
+ else{
+ it->second.push_back(t);
+ }
}
+ }
}
-
-void buildEdgeToTriangle ( GFace *gf , e2t_cont &adj )
+void buildEdgeToTriangle(GFace *gf, e2t_cont &adj)
{
- buildEdgeToTriangle(gf->triangles,adj);
+ buildEdgeToTriangle(gf->triangles, adj);
}
-void buildEdgeToTriangle ( std::vector<MTriangle*> &triangles , e2t_cont &adj )
+void buildEdgeToTriangle(std::vector<MTriangle*> &triangles, e2t_cont &adj)
{
adj.clear();
- for (unsigned int i=0;i<triangles.size();i++)
- {
- MTriangle *t = triangles[i];
- for (unsigned int j=0;j<3;j++)
+ for (unsigned int i = 0; i < triangles.size(); i++){
+ MTriangle *t = triangles[i];
+ for (unsigned int j = 0; j < 3; j++){
+ MVertex *v1 = t->getVertex(j);
+ MVertex *v2 = t->getVertex((j + 1) % 3);
+ MEdge e(v1, v2);
+ e2t_cont::iterator it = adj.find(e);
+ if (it == adj.end()){
+ std::pair<MTriangle*, MTriangle*> one = std::make_pair(t, (MTriangle*)0);
+ adj[e] = one;
+ }
+ else
{
- MVertex *v1 = t->getVertex(j);
- MVertex *v2 = t->getVertex((j+1)%3);
- MEdge e(v1,v2);
- e2t_cont :: iterator it = adj.find ( e );
- if (it == adj.end())
- {
- std::pair<MTriangle*,MTriangle*> one = std::make_pair (t,(MTriangle*)0);
- adj[e] = one;
- }
- else
- {
- it->second.second = t;
- }
+ it->second.second = t;
}
}
+ }
}
-
-void parametricCoordinates ( MTriangle *t , GFace *gf, double u[3], double v[3])
+void parametricCoordinates(MTriangle *t, GFace *gf, double u[3], double v[3])
{
- for (unsigned int j=0;j<3;j++)
- {
- MVertex *ver = t->getVertex(j);
- parametricCoordinates ( ver , gf, u[j], v[j]);
- }
+ for (unsigned int j = 0; j < 3; j++){
+ MVertex *ver = t->getVertex(j);
+ parametricCoordinates(ver, gf, u[j], v[j]);
+ }
}
-void laplaceSmoothing (GFace *gf)
+void laplaceSmoothing(GFace *gf)
{
v2t_cont adj;
- buildVertexToTriangle ( gf->triangles , adj );
-
- for (int i=0;i<5;i++)
- {
- v2t_cont :: iterator it = adj.begin();
-
- while (it != adj.end())
- {
- MVertex *ver= it->first;
- GEntity *ge = ver->onWhat();
- // this vertex in internal to the face
- if (ge->dim() == 2)
- {
- double initu,initv;
- ver->getParameter ( 0,initu);
- ver->getParameter ( 1,initv);
-
- const std::vector<MTriangle*> & lt = it->second;
-
- double fact = lt.size() ? 1./(3.*lt.size()):0;
-
- double cu=0,cv=0;
- double pu[3],pv[3];
- for (unsigned int i=0;i<lt.size();i++)
- {
- parametricCoordinates ( lt[i] , gf, pu, pv);
- cu += fact*(pu[0]+pu[1]+pu[2]);
- cv += fact*(pv[0]+pv[1]+pv[2]);
- // have to test validity !
- }
- ver->setParameter(0,cu);
- ver->setParameter(1,cv);
- GPoint pt = gf->point(SPoint2(cu,cv));
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
- }
- ++it;
- }
- }
+ buildVertexToTriangle(gf->triangles, adj);
+
+ for (int i = 0; i < 5; i++){
+ v2t_cont :: iterator it = adj.begin();
+ while (it != adj.end()){
+ MVertex *ver= it->first;
+ GEntity *ge = ver->onWhat();
+ // this vertex in internal to the face
+ if (ge->dim() == 2){
+ double initu,initv;
+ ver->getParameter(0, initu);
+ ver->getParameter(1, initv);
+ const std::vector<MTriangle*> < = it->second;
+ double fact = lt.size() ? 1. / (3. * lt.size()) : 0;
+ double cu = 0, cv = 0;
+ double pu[3], pv[3];
+ for (unsigned int i = 0; i < lt.size(); i++){
+ parametricCoordinates(lt[i], gf, pu, pv);
+ cu += fact * (pu[0] + pu[1] + pu[2]);
+ cv += fact * (pv[0] + pv[1] + pv[2]);
+ // have to test validity !
+ }
+ ver->setParameter(0, cu);
+ ver->setParameter(1, cv);
+ GPoint pt = gf->point(SPoint2(cu, cv));
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
+ }
+ ++it;
+ }
+ }
}
-extern void fourthPoint (double *p1, double *p2, double *p3, double *p4);
+extern void fourthPoint(double *p1, double *p2, double *p3, double *p4);
-double surfaceTriangleUV (MVertex *v1, MVertex *v2, MVertex *v3,
- const std::vector<double> & Us,
- const std::vector<double> & Vs){
- const double v12[2] = {Us[v2->getNum()]-Us[v1->getNum()],Vs[v2->getNum()]-Vs[v1->getNum()]};
- const double v13[2] = {Us[v3->getNum()]-Us[v1->getNum()],Vs[v3->getNum()]-Vs[v1->getNum()]};
+double surfaceTriangleUV(MVertex *v1, MVertex *v2, MVertex *v3,
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs)
+{
+ const double v12[2] = {Us[v2->getNum()] - Us[v1->getNum()],Vs[v2->getNum()] - Vs[v1->getNum()]};
+ const double v13[2] = {Us[v3->getNum()] - Us[v1->getNum()],Vs[v3->getNum()] - Vs[v1->getNum()]};
return 0.5*fabs (v12[0]*v13[1]-v12[1]*v13[0]);
}
-
-bool gmshEdgeSwap(std::set<swapquad> & configs,
- MTri3 *t1,
- GFace *gf,
- int iLocalEdge,
- std::vector<MTri3*> &newTris,
- const gmshSwapCriterion &cr,
- const std::vector<double> & Us,
- const std::vector<double> & Vs,
- const std::vector<double> & vSizes ,
- const std::vector<double> & vSizesBGM){
-
+bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
+ std::vector<MTri3*> &newTris, const gmshSwapCriterion &cr,
+ const std::vector<double> &Us, const std::vector<double> &Vs,
+ const std::vector<double> &vSizes, const std::vector<double> &vSizesBGM)
+{
MTri3 *t2 = t1->getNeigh(iLocalEdge);
if (!t2) return false;
- MVertex *v1 = t1->tri()->getVertex(iLocalEdge==0 ? 2 : iLocalEdge -1);
- MVertex *v2 = t1->tri()->getVertex((iLocalEdge)%3);
- MVertex *v3 = t1->tri()->getVertex((iLocalEdge+1)%3);
- MVertex *v4 = 0 ;
- for (int i=0;i<3;i++)
+ MVertex *v1 = t1->tri()->getVertex(iLocalEdge == 0 ? 2 : iLocalEdge - 1);
+ MVertex *v2 = t1->tri()->getVertex((iLocalEdge) % 3);
+ MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
+ MVertex *v4 = 0;
+ for (int i = 0; i < 3; i++)
if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
v4 = t2->tri()->getVertex(i);
- // printf("%d %d %d %d %d\n",Us.size(),v1->getNum(),v2->getNum(),v3->getNum(),v4->getNum());
-
- // printf("%d %d %d %d\n",tv1,tv2,tv3,tv4);
-
- swapquad sq (v1,v2,v3,v4);
- if (configs.find(sq) != configs.end())return false;
+ swapquad sq (v1, v2, v3, v4);
+ if (configs.find(sq) != configs.end()) return false;
configs.insert(sq);
- // if (tv1 != 0 && tv1 == tv2 && tv1 == tv3 && tv1 == tv4)return false;
-
- const double volumeRef = surfaceTriangleUV (v1,v2,v3,Us,Vs) + surfaceTriangleUV (v1,v2,v4,Us,Vs);
+ const double volumeRef = surfaceTriangleUV(v1, v2, v3, Us, Vs) + surfaceTriangleUV(v1, v2, v4, Us, Vs);
- /// printf("%p %p %p %p\n",v2,v3,v4);
- MTriangle *t1b = new MTriangle (v2,v3,v4);
- /// printf("%p %p %p %p\n",v2,v3,v4);
- MTriangle *t2b = new MTriangle (v4,v3,v1);
-
- const double v1b = surfaceTriangleUV (v2,v3,v4,Us,Vs);
- const double v2b = surfaceTriangleUV (v4,v3,v1,Us,Vs);
- const double volume = v1b+v2b;
- if (fabs(volume-volumeRef) > 1.e-10 * (volume+volumeRef) ||
- v1b < 1.e-8 * (volume+volumeRef) ||
- v2b < 1.e-8 * (volume+volumeRef)){
+ MTriangle *t1b = new MTriangle(v2, v3, v4);
+ MTriangle *t2b = new MTriangle(v4, v3, v1);
+ const double v1b = surfaceTriangleUV(v2, v3, v4, Us, Vs);
+ const double v2b = surfaceTriangleUV(v4, v3, v1, Us, Vs);
+ const double volume = v1b + v2b;
+ if (fabs(volume - volumeRef) > 1.e-10 * (volume + volumeRef) ||
+ v1b < 1.e-8 * (volume + volumeRef) ||
+ v2b < 1.e-8 * (volume + volumeRef)){
delete t1b;
delete t2b;
return false;
@@ -285,8 +242,10 @@ bool gmshEdgeSwap(std::set<swapquad> & configs,
switch(cr){
case SWCR_QUAL:
{
- const double triQualityRef = std::min(qmTriangle(t1->tri(),QMTRI_RHO),qmTriangle(t2->tri(),QMTRI_RHO));
- const double triQuality = std::min(qmTriangle(t1b,QMTRI_RHO),qmTriangle(t2b,QMTRI_RHO));
+ const double triQualityRef = std::min(qmTriangle(t1->tri(), QMTRI_RHO),
+ qmTriangle(t2->tri(), QMTRI_RHO));
+ const double triQuality = std::min(qmTriangle(t1b, QMTRI_RHO),
+ qmTriangle(t2b, QMTRI_RHO));
if (triQuality < triQualityRef){
delete t1b;
delete t2b;
@@ -296,12 +255,14 @@ bool gmshEdgeSwap(std::set<swapquad> & configs,
}
case SWCR_DEL:
{
- double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()]+Us[v3->getNum()]+Us[v4->getNum()])*.25,
- (Vs[v1->getNum()]+Vs[v2->getNum()]+Vs[v3->getNum()]+Vs[v4->getNum()])*.25};
- double uv4[2] ={Us[v4->getNum()],Vs[v4->getNum()]};
+ double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()] + Us[v3->getNum()] +
+ Us[v4->getNum()]) * .25,
+ (Vs[v1->getNum()] + Vs[v2->getNum()] + Vs[v3->getNum()] +
+ Vs[v4->getNum()]) * .25};
+ double uv4[2] ={Us[v4->getNum()], Vs[v4->getNum()]};
double metric[3];
- buildMetric ( gf , edgeCenter , metric);
- if (!inCircumCircleAniso (gf,t1->tri(),uv4,metric,Us,Vs)){
+ buildMetric(gf, edgeCenter, metric);
+ if (!inCircumCircleAniso(gf, t1->tri(), uv4, metric, Us, Vs)){
delete t1b;
delete t2b;
return false;