Skip to content
Snippets Groups Projects
Commit dbee018e authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent 5564ebca
No related branches found
No related tags found
No related merge requests found
// $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,8 +202,10 @@ 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()]};
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]);
}
......@@ -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};
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());
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());
if (d1 < d2){
delete t1b;
delete t2b;
......@@ -320,8 +332,10 @@ bool gmshEdgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalE
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 );
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);
......@@ -338,7 +352,8 @@ 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 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]);
......@@ -351,33 +366,27 @@ inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
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()]) +
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);
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);
MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
MVertex *v4 = 0 ;
......@@ -386,7 +395,8 @@ bool gmshEdgeSplit(const double lMax,
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]);
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);
......@@ -396,10 +406,13 @@ bool gmshEdgeSplit(const double lMax,
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;
......@@ -517,11 +530,8 @@ void computeNeighboringTrisOfACavity(const std::vector<MTri3*> &cavity,
}
}
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,27 +539,15 @@ 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);
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);
MVertex *v3 = t->tri()->getVertex(i);
// printf("--> %d %d\n",v2->getNum(),v3->getNum());
if ((v2 == *v1 && v3 == lastinring) ||
(v2 == lastinring && v3 == *v1)){
iEdge = i;
......@@ -568,12 +566,6 @@ bool gmshBuildVertexCavity(MTri3 *t,
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;
}
}
......@@ -581,22 +573,16 @@ bool gmshBuildVertexCavity(MTri3 *t,
}
}
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;
double l_min = lMin;
......@@ -652,10 +638,8 @@ bool gmshVertexCollapse(const double lMin,
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)
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 ;
......@@ -667,7 +651,10 @@ int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
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;}
if (gmshEdgeSwap(configs, *it, gf, i, newTris, cr, Us, Vs, vSizes, vSizesBGM)){
nbSwap++;
break;
}
}
}
else{
......@@ -679,23 +666,16 @@ int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
}
}
allTris.insert(newTris.begin(), newTris.end());
// printf("iter %d nbswam %d\n",iter,nbSwap);
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,
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)
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;
......@@ -705,7 +685,10 @@ int edgeSplitPass (double maxLC,
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;}
if (gmshEdgeSplit(maxLC, *it, gf, i, newTris, cr, Us, Vs, vSizes, vSizesBGM)) {
nbSplit++;
break;
}
}
}
else{
......@@ -722,12 +705,9 @@ 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;
......@@ -737,7 +717,10 @@ int edgeCollapsePass (double minLC,
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;}
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;
}
$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.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment