From dbee018ed13aadc55f60b16e1edfa5fe4de69b1e Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 1 Mar 2008 10:52:05 +0000
Subject: [PATCH] *** empty log message ***
---
Mesh/meshGFaceOptimize.cpp | 410 ++++++++++++++++++-------------------
doc/TODO | 12 +-
2 files changed, 197 insertions(+), 225 deletions(-)
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 3160b32f2b..34ab2d0976 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceOptimize.cpp,v 1.11 2008-02-21 13:34:40 geuzaine Exp $
+// $Id: meshGFaceOptimize.cpp,v 1.12 2008-03-01 10:52:05 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -65,10 +65,13 @@ void buidMeshGenerationDataStructures(GFace *gf, std::set<MTri3*, compareTri3Ptr
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);
+ 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){
+ 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);
@@ -199,9 +202,11 @@ 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]);
+ 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,
@@ -224,7 +229,8 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
if (configs.find(sq) != configs.end()) return false;
configs.insert(sq);
- 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);
MTriangle *t1b = new MTriangle(v2, v3, v4);
MTriangle *t2b = new MTriangle(v4, v3, v1);
@@ -271,13 +277,19 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
break;
case SWCR_CLOSE:
{
- double avg1[3] = {(v1->x() + v2->x()) *.5,(v1->y() + v2->y()) *.5,(v1->z() + v2->z()) *.5};
- double avg2[3] = {(v3->x() + v4->x()) *.5,(v3->y() + v4->y()) *.5,(v3->z() + v4->z()) *.5};
+ double avg1[3] = {(v1->x() + v2->x()) *.5,(v1->y() + v2->y()) *.5,
+ (v1->z() + v2->z()) *.5};
+ double avg2[3] = {(v3->x() + v4->x()) *.5,(v3->y() + v4->y()) *.5,
+ (v3->z() + v4->z()) *.5};
- GPoint gp1 = gf->point(SPoint2((Us[v1->getNum()]+Us[v2->getNum()])*.5,(Vs[v1->getNum()]+Vs[v2->getNum()])*.5));
- GPoint gp2 = gf->point(SPoint2((Us[v3->getNum()]+Us[v4->getNum()])*.5,(Vs[v3->getNum()]+Vs[v4->getNum()])*.5));
- double d1 = (avg1[0]-gp1.x()) *(avg1[0]-gp1.x()) +(avg1[1]-gp1.y()) *(avg1[1]-gp1.y()) +(avg1[2]-gp1.z()) *(avg1[2]-gp1.z());
- double d2 = (avg2[0]-gp2.x()) *(avg2[0]-gp2.x()) +(avg2[1]-gp2.y()) *(avg2[1]-gp2.y()) +(avg2[2]-gp2.z()) *(avg2[2]-gp2.z());
+ GPoint gp1 = gf->point(SPoint2((Us[v1->getNum()] + Us[v2->getNum()]) * .5,
+ (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5));
+ GPoint gp2 = gf->point(SPoint2((Us[v3->getNum()] + Us[v4->getNum()]) * .5,
+ (Vs[v3->getNum()] + Vs[v4->getNum()]) * .5));
+ double d1 = (avg1[0] - gp1.x()) * (avg1[0] - gp1.x()) + (avg1[1] - gp1.y()) *
+ (avg1[1]-gp1.y()) + (avg1[2] - gp1.z()) * (avg1[2] - gp1.z());
+ double d2 = (avg2[0] - gp2.x()) * (avg2[0] - gp2.x()) + (avg2[1] - gp2.y()) *
+ (avg2[1] - gp2.y()) + (avg2[2] - gp2.z()) * (avg2[2] - gp2.z());
if (d1 < d2){
delete t1b;
delete t2b;
@@ -290,44 +302,46 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
}
std::list<MTri3*> cavity;
- for(int i=0;i<3;i++){
+ for(int i = 0; i < 3; i++){
if (t1->getNeigh(i) && t1->getNeigh(i) != t2){
bool found = false;
- for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
- if (*it == t1->getNeigh(i))found = true;
+ for (std::list<MTri3*>::iterator it = cavity.begin(); it != cavity.end(); it++){
+ if (*it == t1->getNeigh(i)) found = true;
}
if (!found)cavity.push_back(t1->getNeigh(i));
}
}
- for(int i=0;i<3;i++){
+ for(int i = 0; i < 3; i++){
if (t2->getNeigh(i) && t2->getNeigh(i) != t1){
bool found = false;
- for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
- if (*it == t2->getNeigh(i))found = true;
+ for (std::list<MTri3*>::iterator it = cavity.begin(); it != cavity.end(); it++){
+ if (*it == t2->getNeigh(i)) found = true;
}
if (!found)cavity.push_back(t2->getNeigh(i));
}
}
- double lc1 = 0.3333333333*(vSizes [t1b->getVertex(0)->getNum()]+
- vSizes [t1b->getVertex(1)->getNum()]+
- vSizes [t1b->getVertex(2)->getNum()]);
- double lcBGM1 = 0.3333333333*(vSizesBGM [t1b->getVertex(0)->getNum()]+
- vSizesBGM [t1b->getVertex(1)->getNum()]+
- vSizesBGM [t1b->getVertex(2)->getNum()]);
- double lc2 = 0.3333333333*(vSizes [t2b->getVertex(0)->getNum()]+
- vSizes [t2b->getVertex(1)->getNum()]+
- vSizes [t2b->getVertex(2)->getNum()]);
- double lcBGM2 = 0.3333333333*(vSizesBGM [t2b->getVertex(0)->getNum()]+
- vSizesBGM [t2b->getVertex(1)->getNum()]+
- vSizesBGM [t2b->getVertex(2)->getNum()]);
- MTri3 *t1b3 = new MTri3 ( t1b , Extend1dMeshIn2dSurfaces() ? std::min(lc1,lcBGM1) : lcBGM1 );
- MTri3 *t2b3 = new MTri3 ( t2b , Extend1dMeshIn2dSurfaces() ? std::min(lc2,lcBGM2) : lcBGM2 );
+ double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getNum()] +
+ vSizes[t1b->getVertex(1)->getNum()] +
+ vSizes[t1b->getVertex(2)->getNum()]);
+ double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getNum()] +
+ vSizesBGM[t1b->getVertex(1)->getNum()] +
+ vSizesBGM[t1b->getVertex(2)->getNum()]);
+ double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getNum()] +
+ vSizes[t2b->getVertex(1)->getNum()] +
+ vSizes[t2b->getVertex(2)->getNum()]);
+ double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getNum()] +
+ vSizesBGM[t2b->getVertex(1)->getNum()] +
+ vSizesBGM[t2b->getVertex(2)->getNum()]);
+ MTri3 *t1b3 = new MTri3(t1b, Extend1dMeshIn2dSurfaces() ?
+ std::min(lc1, lcBGM1) : lcBGM1);
+ MTri3 *t2b3 = new MTri3(t2b, Extend1dMeshIn2dSurfaces() ?
+ std::min(lc2, lcBGM2) : lcBGM2);
cavity.push_back(t1b3);
cavity.push_back(t2b3);
t1->setDeleted(true);
t2->setDeleted(true);
- connectTriangles ( cavity );
+ connectTriangles(cavity);
newTris.push_back(t1b3);
newTris.push_back(t2b3);
@@ -335,71 +349,70 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
}
inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
- const std::vector<double> & Us,
- const std::vector<double> & Vs,
- const std::vector<double> & vSizes ,
- const std::vector<double> & vSizesBGM ){
- const double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
- (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
- GPoint GP = f->point (edgeCenter[0],edgeCenter[1]);
+ const std::vector<double> &Us,
+ const std::vector<double> &Vs,
+ const std::vector<double> &vSizes ,
+ const std::vector<double> &vSizesBGM)
+{
+ const double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()]) * .5,
+ (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5};
+ GPoint GP = f->point (edgeCenter[0], edgeCenter[1]);
- const double dx1 = v1->x()-GP.x();
- const double dy1 = v1->y()-GP.y();
- const double dz1 = v1->z()-GP.z();
- const double l1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
- const double dx2 = v2->x()-GP.x();
- const double dy2 = v2->y()-GP.y();
- const double dz2 = v2->z()-GP.z();
- const double l2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2);
- if (Extend1dMeshIn2dSurfaces())return 2*(l1+l2)/(std::min(vSizes[v1->getNum()],vSizesBGM[v1->getNum()]) +
- std::min(vSizes[v2->getNum()],vSizesBGM[v2->getNum()]));
- return 2*(l1+l2)/(vSizesBGM[v1->getNum()] + vSizesBGM[v2->getNum()]);
+ const double dx1 = v1->x() - GP.x();
+ const double dy1 = v1->y() - GP.y();
+ const double dz1 = v1->z() - GP.z();
+ const double l1 = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ const double dx2 = v2->x() - GP.x();
+ const double dy2 = v2->y() - GP.y();
+ const double dz2 = v2->z() - GP.z();
+ const double l2 = sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
+ if (Extend1dMeshIn2dSurfaces())
+ return 2 * (l1 + l2) / (std::min(vSizes[v1->getNum()], vSizesBGM[v1->getNum()]) +
+ std::min(vSizes[v2->getNum()], vSizesBGM[v2->getNum()]));
+ return 2 * (l1 + l2) / (vSizesBGM[v1->getNum()] + vSizesBGM[v2->getNum()]);
}
-bool gmshEdgeSplit(const double lMax,
- MTri3 *t1,
- GFace *gf,
- int iLocalEdge,
- std::vector<MTri3*> &newTris,
- const gmshSplitCriterion &cr,
- std::vector<double> & Us,
- std::vector<double> & Vs,
- std::vector<double> & vSizes ,
- std::vector<double> & vSizesBGM ){
-
+bool gmshEdgeSplit(const double lMax, MTri3 *t1, GFace *gf, int iLocalEdge,
+ std::vector<MTri3*> &newTris, const gmshSplitCriterion &cr,
+ std::vector<double> &Us, std::vector<double> &Vs,
+ std::vector<double> &vSizes, 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);
- double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
- (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
+ MVertex *v1 = t1->tri()->getVertex(iLocalEdge == 0 ? 2 : iLocalEdge - 1);
+ MVertex *v2 = t1->tri()->getVertex(iLocalEdge % 3);
+ double edgeCenter[2] ={(Us[v1->getNum()] + Us[v2->getNum()]) * .5,
+ (Vs[v1->getNum()] + Vs[v2->getNum()]) * .5};
- double al = computeEdgeAdimLength(v1, v2, gf,Us,Vs,vSizes,vSizesBGM);
- if (al <= lMax )return false;
- // printf("al,lMax %12.5E %12.5E \n",al,lMax);
+ double al = computeEdgeAdimLength(v1, v2, gf,Us, Vs, vSizes, vSizesBGM);
+ if (al <= lMax) return false;
- MVertex *v3 = t1->tri()->getVertex((iLocalEdge+1)%3);
+ MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
MVertex *v4 = 0 ;
- for (int i=0;i<3;i++)
+ for (int i = 0; i < 3; i++)
if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
v4 = t2->tri()->getVertex(i);
- GPoint p = gf->point (edgeCenter[0],edgeCenter[1]);
- MVertex *vnew = new MFaceVertex (p.x(),p.y(),p.z(),gf,edgeCenter[0],edgeCenter[1]);
+ GPoint p = gf->point(edgeCenter[0], edgeCenter[1]);
+ MVertex *vnew = new MFaceVertex(p.x(), p.y(), p.z(), gf,
+ edgeCenter[0], edgeCenter[1]);
- MTriangle *t1b = new MTriangle (v1,vnew,v3);
- MTriangle *t2b = new MTriangle (vnew,v2,v3);
- MTriangle *t3b = new MTriangle (v1,v4,vnew);
- MTriangle *t4b = new MTriangle (v4,v2,vnew);
+ MTriangle *t1b = new MTriangle(v1, vnew, v3);
+ MTriangle *t2b = new MTriangle(vnew, v2, v3);
+ MTriangle *t3b = new MTriangle(v1, v4, vnew);
+ MTriangle *t4b = new MTriangle(v4, v2, vnew);
switch(cr){
case SPCR_QUAL:
{
- const double triQualityRef = std::min(qmTriangle(t1->tri(),QMTRI_RHO),qmTriangle(t2->tri(),QMTRI_RHO));
+ const double triQualityRef = std::min(qmTriangle(t1->tri(), QMTRI_RHO),
+ qmTriangle(t2->tri(), QMTRI_RHO));
const double triQuality =
- std::min(std::min(std::min(qmTriangle(t1b,QMTRI_RHO),qmTriangle(t2b,QMTRI_RHO)),
- qmTriangle(t3b,QMTRI_RHO)),qmTriangle(t4b,QMTRI_RHO));
+ std::min(std::min(std::min(qmTriangle(t1b, QMTRI_RHO),
+ qmTriangle(t2b, QMTRI_RHO)),
+ qmTriangle(t3b, QMTRI_RHO)),
+ qmTriangle(t4b, QMTRI_RHO));
if (triQuality < triQualityRef){
delete t1b;
delete t2b;
@@ -418,61 +431,61 @@ bool gmshEdgeSplit(const double lMax,
gf->mesh_vertices.push_back(vnew);
vnew->setNum(Us.size());
- double lcN = 0.5 * (vSizes [v1->getNum()] + vSizes [v2->getNum()] );
- double lcBGM = BGM_MeshSize(gf,edgeCenter[0],edgeCenter[1],p.x(),p.y(),p.z());
+ double lcN = 0.5 * (vSizes [v1->getNum()] + vSizes [v2->getNum()] );
+ double lcBGM = BGM_MeshSize(gf, edgeCenter[0], edgeCenter[1], p.x(), p.y(), p.z());
- vSizesBGM.push_back( lcBGM );
- vSizes.push_back ( lcN);
- Us.push_back( edgeCenter[0] );
- Vs.push_back( edgeCenter[1] );
+ vSizesBGM.push_back(lcBGM);
+ vSizes.push_back(lcN);
+ Us.push_back(edgeCenter[0]);
+ Vs.push_back(edgeCenter[1]);
std::list<MTri3*> cavity;
- for(int i=0;i<3;i++){
+ for(int i = 0; i < 3; i++){
if (t1->getNeigh(i) && t1->getNeigh(i) != t2){
bool found = false;
- for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
- if (*it == t1->getNeigh(i))found = true;
+ for (std::list<MTri3*>::iterator it = cavity.begin();it != cavity.end(); it++){
+ if (*it == t1->getNeigh(i)) found = true;
}
- if (!found)cavity.push_back(t1->getNeigh(i));
+ if (!found) cavity.push_back(t1->getNeigh(i));
}
}
- for(int i=0;i<3;i++){
+ for(int i = 0; i < 3; i++){
if (t2->getNeigh(i) && t2->getNeigh(i) != t1){
bool found = false;
- for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
+ for (std::list<MTri3*>::iterator it = cavity.begin(); it != cavity.end(); it++){
if (*it == t2->getNeigh(i))found = true;
}
- if (!found)cavity.push_back(t2->getNeigh(i));
+ if (!found) cavity.push_back(t2->getNeigh(i));
}
}
- double lc1 = 0.3333333333*(vSizes [t1b->getVertex(0)->getNum()]+
- vSizes [t1b->getVertex(1)->getNum()]+
- vSizes [t1b->getVertex(2)->getNum()]);
- double lcBGM1 = 0.3333333333*(vSizesBGM [t1b->getVertex(0)->getNum()]+
- vSizesBGM [t1b->getVertex(1)->getNum()]+
- vSizesBGM [t1b->getVertex(2)->getNum()]);
- double lc2 = 0.3333333333*(vSizes [t2b->getVertex(0)->getNum()]+
- vSizes [t2b->getVertex(1)->getNum()]+
- vSizes [t2b->getVertex(2)->getNum()]);
- double lcBGM2 = 0.3333333333*(vSizesBGM [t2b->getVertex(0)->getNum()]+
- vSizesBGM [t2b->getVertex(1)->getNum()]+
- vSizesBGM [t2b->getVertex(2)->getNum()]);
- double lc3 = 0.3333333333*(vSizes [t3b->getVertex(0)->getNum()]+
- vSizes [t3b->getVertex(1)->getNum()]+
- vSizes [t3b->getVertex(2)->getNum()]);
- double lcBGM3 = 0.3333333333*(vSizesBGM [t3b->getVertex(0)->getNum()]+
- vSizesBGM [t3b->getVertex(1)->getNum()]+
- vSizesBGM [t3b->getVertex(2)->getNum()]);
- double lc4 = 0.3333333333*(vSizes [t4b->getVertex(0)->getNum()]+
- vSizes [t4b->getVertex(1)->getNum()]+
- vSizes [t4b->getVertex(2)->getNum()]);
- double lcBGM4 = 0.3333333333*(vSizesBGM [t4b->getVertex(0)->getNum()]+
- vSizesBGM [t4b->getVertex(1)->getNum()]+
- vSizesBGM [t4b->getVertex(2)->getNum()]);
- MTri3 *t1b3 = new MTri3 ( t1b , Extend1dMeshIn2dSurfaces() ? std::min(lc1,lcBGM1) : lcBGM1 );
- MTri3 *t2b3 = new MTri3 ( t2b , Extend1dMeshIn2dSurfaces() ? std::min(lc2,lcBGM2) : lcBGM2 );
- MTri3 *t3b3 = new MTri3 ( t3b , Extend1dMeshIn2dSurfaces() ? std::min(lc3,lcBGM3) : lcBGM3 );
- MTri3 *t4b3 = new MTri3 ( t4b , Extend1dMeshIn2dSurfaces() ? std::min(lc4,lcBGM4) : lcBGM4 );
+ double lc1 = 0.3333333333 * (vSizes[t1b->getVertex(0)->getNum()] +
+ vSizes[t1b->getVertex(1)->getNum()] +
+ vSizes[t1b->getVertex(2)->getNum()]);
+ double lcBGM1 = 0.3333333333 * (vSizesBGM[t1b->getVertex(0)->getNum()] +
+ vSizesBGM[t1b->getVertex(1)->getNum()] +
+ vSizesBGM[t1b->getVertex(2)->getNum()]);
+ double lc2 = 0.3333333333 * (vSizes[t2b->getVertex(0)->getNum()] +
+ vSizes[t2b->getVertex(1)->getNum()] +
+ vSizes[t2b->getVertex(2)->getNum()]);
+ double lcBGM2 = 0.3333333333 * (vSizesBGM[t2b->getVertex(0)->getNum()] +
+ vSizesBGM[t2b->getVertex(1)->getNum()] +
+ vSizesBGM[t2b->getVertex(2)->getNum()]);
+ double lc3 = 0.3333333333 * (vSizes[t3b->getVertex(0)->getNum()] +
+ vSizes[t3b->getVertex(1)->getNum()] +
+ vSizes[t3b->getVertex(2)->getNum()]);
+ double lcBGM3 = 0.3333333333 * (vSizesBGM[t3b->getVertex(0)->getNum()] +
+ vSizesBGM[t3b->getVertex(1)->getNum()] +
+ vSizesBGM[t3b->getVertex(2)->getNum()]);
+ double lc4 = 0.3333333333 * (vSizes[t4b->getVertex(0)->getNum()] +
+ vSizes[t4b->getVertex(1)->getNum()] +
+ vSizes[t4b->getVertex(2)->getNum()]);
+ double lcBGM4 = 0.3333333333 * (vSizesBGM[t4b->getVertex(0)->getNum()] +
+ vSizesBGM[t4b->getVertex(1)->getNum()] +
+ vSizesBGM[t4b->getVertex(2)->getNum()]);
+ MTri3 *t1b3 = new MTri3(t1b, Extend1dMeshIn2dSurfaces() ? std::min(lc1, lcBGM1) : lcBGM1);
+ MTri3 *t2b3 = new MTri3(t2b, Extend1dMeshIn2dSurfaces() ? std::min(lc2, lcBGM2) : lcBGM2);
+ MTri3 *t3b3 = new MTri3(t3b, Extend1dMeshIn2dSurfaces() ? std::min(lc3, lcBGM3) : lcBGM3);
+ MTri3 *t4b3 = new MTri3(t4b, Extend1dMeshIn2dSurfaces() ? std::min(lc4, lcBGM4) : lcBGM4);
cavity.push_back(t1b3);
cavity.push_back(t2b3);
@@ -480,7 +493,7 @@ bool gmshEdgeSplit(const double lMax,
cavity.push_back(t4b3);
t1->setDeleted(true);
t2->setDeleted(true);
- connectTriangles ( cavity );
+ connectTriangles(cavity);
newTris.push_back(t1b3);
newTris.push_back(t2b3);
newTris.push_back(t3b3);
@@ -511,17 +524,14 @@ void computeNeighboringTrisOfACavity(const std::vector<MTri3*> &cavity,
}
}
}
- if(!found)outside.push_back(neigh);
+ if(!found) outside.push_back(neigh);
}
}
}
}
-bool gmshBuildVertexCavity(MTri3 *t,
- int iLocalVertex,
- MVertex **v1,
- std::vector<MTri3*> &cavity,
- std::vector<MTri3*> &outside,
+bool gmshBuildVertexCavity(MTri3 *t, int iLocalVertex, MVertex **v1,
+ std::vector<MTri3*> &cavity, std::vector<MTri3*> &outside,
std::vector<MVertex*> &ring)
{
cavity.clear();
@@ -529,75 +539,51 @@ bool gmshBuildVertexCavity(MTri3 *t,
*v1 = t->tri()->getVertex(iLocalVertex);
-// printf("VERTEX %d\n",
-// t->tri()->getVertex(iLocalVertex)->getNum());
-
- MVertex *lastinring = t->tri()->getVertex((iLocalVertex+1)%3);
+ MVertex *lastinring = t->tri()->getVertex((iLocalVertex + 1) % 3);
ring.push_back(lastinring);
cavity.push_back(t);
-// printf("START triangle %d %d %d, vertex %d\n",
-// t->tri()->getVertex(0)->getNum(),
-// t->tri()->getVertex(1)->getNum(),
-// t->tri()->getVertex(2)->getNum(),
-// lastinring->getNum());
-
-
while (1){
int iEdge = -1;
- // printf("look for %d %d\n",(*v1)->getNum(),lastinring->getNum());
- for (int i=0;i<3;i++){
- MVertex *v2 = t->tri()->getVertex((i+2)%3);
+ for (int i = 0; i < 3; i++){
+ MVertex *v2 = t->tri()->getVertex((i + 2) % 3);
MVertex *v3 = t->tri()->getVertex(i);
- // printf("--> %d %d\n",v2->getNum(),v3->getNum());
- if ( (v2 == *v1 && v3 == lastinring ) ||
- (v2 == lastinring && v3 == *v1 )){
+ if ((v2 == *v1 && v3 == lastinring) ||
+ (v2 == lastinring && v3 == *v1)){
iEdge = i;
t = t->getNeigh(i);
- if (t==cavity[0]) {
- computeNeighboringTrisOfACavity (cavity,outside);
+ if (t == cavity[0]) {
+ computeNeighboringTrisOfACavity(cavity, outside);
return true;
}
- if (!t)return false;
- if (t->isDeleted()){printf("weird\n");throw;}
+ if (!t) return false;
+ if (t->isDeleted()){ printf("weird\n"); throw; }
cavity.push_back(t);
- for (int j=0;j<3;j++){
+ for (int j = 0; j < 3; j++){
if (t->tri()->getVertex(j) !=lastinring && t->tri()->getVertex(j) != *v1){
lastinring = t->tri()->getVertex(j);
ring.push_back(lastinring);
j = 100;
}
}
-// printf("CONTINUE (%d) triangle %p = %d %d %d, vertex %d size %d\n",i,
-// t,
-// t->tri()->getVertex(0)->getNum(),
-// t->tri()->getVertex(1)->getNum(),
-// t->tri()->getVertex(2)->getNum(),
-// lastinring->getNum(),ring.size());
break;
}
}
- if (iEdge == -1) {printf("not found\n"); throw;}
+ if (iEdge == -1) { printf("not found\n"); throw; }
}
}
-
-
-bool gmshVertexCollapse(const double lMin,
- MTri3 *t1,
- GFace *gf,
- int iLocalVertex,
- std::vector<MTri3*> &newTris,
- std::vector<double> & Us,
- std::vector<double> & Vs,
- std::vector<double> & vSizes ,
- std::vector<double> & vSizesBGM ){
+bool gmshVertexCollapse(const double lMin, MTri3 *t1, GFace *gf,
+ int iLocalVertex, std::vector<MTri3*> &newTris,
+ std::vector<double> &Us, std::vector<double> &Vs,
+ std::vector<double> &vSizes, std::vector<double> &vSizesBGM)
+{
MVertex *v;
std::vector<MTri3*> cavity;
std::vector<MTri3*> outside;
std::vector<MVertex*> ring ;
- // printf("%p \n",t1);
- if (!gmshBuildVertexCavity (t1,iLocalVertex, &v,cavity,outside,ring) )return false;
+
+ if (!gmshBuildVertexCavity(t1, iLocalVertex, &v, cavity, outside, ring)) return false;
double l_min = lMin;
int iMin = -1;
@@ -650,24 +636,25 @@ bool gmshVertexCollapse(const double lMin,
return true;
}
-int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
- const gmshSwapCriterion &cr,
- const std::vector<double> & Us ,
- const std::vector<double> & Vs,
- const std::vector<double> & vSizes ,
- const std::vector<double> & vSizesBGM)
+int edgeSwapPass (GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
+ const gmshSwapCriterion &cr,
+ const std::vector<double> &Us, const std::vector<double> &Vs,
+ const std::vector<double> &vSizes, const std::vector<double> &vSizesBGM)
{
- typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ typedef std::set<MTri3*, compareTri3Ptr> CONTAINER ;
- int nbSwapTot=0;
+ int nbSwapTot = 0;
std::set<swapquad> configs;
- for (int iter=0;iter<1200;iter++){
+ for (int iter = 0; iter < 1200; iter++){
int nbSwap = 0;
std::vector<MTri3*> newTris;
- for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
+ for (CONTAINER::iterator it = allTris.begin(); it != allTris.end(); ++it){
if (!(*it)->isDeleted()){
- for (int i=0;i<3;i++){
- if (gmshEdgeSwap(configs,*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSwap++;break;}
+ for (int i = 0; i < 3; i++){
+ if (gmshEdgeSwap(configs, *it, gf, i, newTris, cr, Us, Vs, vSizes, vSizesBGM)){
+ nbSwap++;
+ break;
+ }
}
}
else{
@@ -675,37 +662,33 @@ int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
CONTAINER::iterator itb = it;
++it;
allTris.erase(itb);
- if (it == allTris.end())break;
+ if (it == allTris.end()) break;
}
}
- allTris.insert(newTris.begin(),newTris.end());
- // printf("iter %d nbswam %d\n",iter,nbSwap);
- nbSwapTot+=nbSwap;
- if (nbSwap == 0)break;
+ allTris.insert(newTris.begin(), newTris.end());
+ nbSwapTot += nbSwap;
+ if (nbSwap == 0) break;
}
-
- // printf("B %d %d tris ",allTris.size(),newTris.size());
- // printf("A %d %d tris\n",allTris.size(),newTris.size());
return nbSwapTot;
}
-int edgeSplitPass (double maxLC,
- GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
- const gmshSplitCriterion &cr,
- std::vector<double> & Us ,
- std::vector<double> & Vs,
- std::vector<double> & vSizes ,
- std::vector<double> & vSizesBGM)
+int edgeSplitPass(double maxLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ const gmshSplitCriterion &cr,
+ std::vector<double> &Us, std::vector<double> &Vs,
+ std::vector<double> &vSizes, std::vector<double> &vSizesBGM)
{
- typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ typedef std::set<MTri3*, compareTri3Ptr> CONTAINER ;
std::vector<MTri3*> newTris;
int nbSplit = 0;
- for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
+ for (CONTAINER::iterator it = allTris.begin(); it != allTris.end(); ++it){
if (!(*it)->isDeleted()){
- for (int i=0;i<3;i++){
- if (gmshEdgeSplit (maxLC,*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSplit++;break;}
+ for (int i = 0; i < 3; i++){
+ if (gmshEdgeSplit(maxLC, *it, gf, i, newTris, cr, Us, Vs, vSizes, vSizesBGM)) {
+ nbSplit++;
+ break;
+ }
}
}
else{
@@ -713,7 +696,7 @@ int edgeSplitPass (double maxLC,
delete *it;
++it;
allTris.erase(itb);
- if (it == allTris.end())break;
+ if (it == allTris.end()) break;
}
}
printf("B %d %d tris ", (int)allTris.size(), (int)newTris.size());
@@ -722,22 +705,22 @@ int edgeSplitPass (double maxLC,
return nbSplit;
}
-int edgeCollapsePass (double minLC,
- GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
- std::vector<double> & Us ,
- std::vector<double> & Vs,
- std::vector<double> & vSizes ,
- std::vector<double> & vSizesBGM)
+int edgeCollapsePass(double minLC, GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ std::vector<double> &Us, std::vector<double> &Vs,
+ std::vector<double> &vSizes, std::vector<double> &vSizesBGM)
{
typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
std::vector<MTri3*> newTris;
int nbCollapse = 0;
- for (CONTAINER::reverse_iterator it = allTris.rbegin();it!=allTris.rend();++it){
+ for (CONTAINER::reverse_iterator it = allTris.rbegin(); it != allTris.rend(); ++it){
if (!(*it)->isDeleted()){
- for (int i=0;i<3;i++){
- if (gmshVertexCollapse (minLC,*it,gf,i,newTris,Us,Vs,vSizes,vSizesBGM)) {nbCollapse++; break;}
+ for (int i = 0; i < 3; i++){
+ if (gmshVertexCollapse(minLC, *it, gf, i, newTris, Us, Vs, vSizes, vSizesBGM)) {
+ nbCollapse++;
+ break;
+ }
}
}
// else{
@@ -754,4 +737,3 @@ int edgeCollapsePass (double minLC,
printf("A %d %d tris\n", (int)allTris.size(), (int)newTris.size());
return nbCollapse;
}
-
diff --git a/doc/TODO b/doc/TODO
index ed77c7ec93..39bca7d328 100644
--- a/doc/TODO
+++ b/doc/TODO
@@ -1,4 +1,4 @@
-$Id: TODO,v 1.68 2008-02-17 10:24:50 geuzaine Exp $
+$Id: TODO,v 1.69 2008-03-01 10:52:05 geuzaine Exp $
********************************************************************
@@ -215,13 +215,3 @@ create "Volume visualization" range type? (interpolate on regular grid
+ create cut planes // to viewpoint with transparency; can be done in
a straightforward way or using 3D textures)
-********************************************************************
-
-Yves Krahenbuhl wrote:
-
-> Lors de la creation des elements du 2eme ordre, et selon la courbure
-> du contour exterieur, le jacobien de l'element peut devenir negatif.
-> Cependant le programme genere le maillage sans protester. Il
-> faudrait peut-etre faire apparaitre un message d'erreurs / ou se
-> restreindre (automatiquement ou non) a une interpolation lineaire
-> entre les points en question.
--
GitLab