From c0da83017c0218133d56ac7428bd3cb2de490cd9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 8 Jul 2013 09:58:42 +0000
Subject: [PATCH] trying to fix Delaunay point insertion when the surface
triangulation is anisotropic: instead of using the min edge length of tets in
the initial mesh, we - use the max edge length of boundary triangles - and
the min edge length of tets for the remaining vertices
---
Mesh/meshGRegionDelaunayInsertion.cpp | 49 ++++++++++++++++++++++++---
1 file changed, 44 insertions(+), 5 deletions(-)
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 7038278d3a..d1455b602a 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -513,8 +513,28 @@ bool insertVertex(MVertex *v,
}
+static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes,
+ std::set<MVertex*> &bndVertices)
+{
+ for(int i = 0; i < 3; i++){
+ bndVertices.insert(t->getVertex(i));
+ MEdge e = t->getEdge(i);
+ MVertex *vi = e.getVertex(0);
+ MVertex *vj = e.getVertex(1);
+ 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);
+ // use largest edge length
+ if (iti == vSizes.end() || iti->second < l) vSizes[vi] = l;
+ if (itj == vSizes.end() || itj->second < l) vSizes[vj] = l;
+ }
+}
-static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
+static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes,
+ std::set<MVertex*> &bndVertices)
{
for (int i = 0; i < 4; i++){
for (int j = i + 1; j < 4; j++){
@@ -526,8 +546,13 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
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 == vSizes.end() || iti->second > l) vSizes[vi] = l;
- if (itj == vSizes.end() || itj->second > l) vSizes[vj] = l;
+ std::set<MVertex*>::iterator itvi = bndVertices.find(vi);
+ std::set<MVertex*>::iterator itvj = bndVertices.find(vj);
+ // smallest tet edge
+ if (itvi == bndVertices.end() &&
+ (iti == vSizes.end() || iti->second > l)) vSizes[vi] = l;
+ if (itvj == bndVertices.end() &&
+ (itj == vSizes.end() || itj->second > l)) vSizes[vj] = l;
}
}
}
@@ -1130,8 +1155,15 @@ void insertVerticesInRegion (GRegion *gr)
{ // leave this in a block so the map gets deallocated directly
std::map<MVertex*, double> vSizesMap;
+ std::set<MVertex*> bndVertices;
+ for(GModel::fiter it = gr->model()->firstFace(); it != gr->model()->lastFace(); ++it){
+ GFace *gf = *it;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ setLcs(gf->triangles[i], vSizesMap, bndVertices);
+ }
+ }
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- setLcs(gr->tetrahedra[i], vSizesMap);
+ setLcs(gr->tetrahedra[i], vSizesMap, bndVertices);
for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
it->first->setIndex(NUM++);
@@ -1402,8 +1434,15 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
{ // leave this in a block so the map gets deallocated directly
std::map<MVertex*, double> vSizesMap;
+ std::set<MVertex*> bndVertices;
+ for(GModel::fiter it = gr->model()->firstFace(); it != gr->model()->lastFace(); ++it){
+ GFace *gf = *it;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ setLcs(gf->triangles[i], vSizesMap, bndVertices);
+ }
+ }
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
- setLcs(gr->tetrahedra[i], vSizesMap);
+ setLcs(gr->tetrahedra[i], vSizesMap, bndVertices);
for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
it->first->setIndex(NUM++);
--
GitLab