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

make sure to re-orient 2D mesh before meshing in 3D (fixes e.g. bug when 2D... 

make sure to re-orient 2D mesh before meshing in 3D (fixes e.g. bug when 2D Frontal leaves badly oriented elements, that 3D frontal cannot handle)
parent 80486220
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 33 additions and 29 deletions
Mesh/Generator.cpp
+ 33
− 29
View file @ 11ae0eb1
...@@ -48,7 +48,7 @@ ...@@ -48,7 +48,7 @@
#include "PViewData.h" #include "PViewData.h"
#endif #endif
template<class T> template<class T>
static void GetQualityMeasure(std::vector<T*> &ele, static void GetQualityMeasure(std::vector<T*> &ele,
double &gamma, double &gammaMin, double &gammaMax, double &gamma, double &gammaMin, double &gammaMax,
double &minSICN, double &minSICNMin, double &minSICNMax, double &minSICN, double &minSICNMin, double &minSICNMax,
...@@ -255,6 +255,9 @@ static void Mesh1D(GModel *m) ...@@ -255,6 +255,9 @@ static void Mesh1D(GModel *m)
if(nIter++ > 10) break; if(nIter++ > 10) break;
} }
// re-orient according to geometrical/user constraints
std::for_each(m->firstEdge(), m->lastEdge(), orientMeshGEdge());
double t2 = Cpu(); double t2 = Cpu();
CTX::instance()->meshTimer[0] = t2 - t1; CTX::instance()->meshTimer[0] = t2 - t1;
Msg::StatusBar(true, "Done meshing 1D (%g s)", CTX::instance()->meshTimer[0]); Msg::StatusBar(true, "Done meshing 1D (%g s)", CTX::instance()->meshTimer[0]);
...@@ -282,8 +285,8 @@ static void PrintMesh2dStatistics(GModel *m) ...@@ -282,8 +285,8 @@ static void PrintMesh2dStatistics(GModel *m)
if(CTX::instance()->createAppendMeshStatReport == 1){ if(CTX::instance()->createAppendMeshStatReport == 1){
fprintf(statreport, "2D stats\tname\t\t#faces\t\t#fail\t\t" fprintf(statreport, "2D stats\tname\t\t#faces\t\t#fail\t\t"
"#t\t\tQavg\t\tQbest\t\tQworst\t\t#Q>90\t\t#Q>90/#t\t" "#t\t\tQavg\t\tQbest\t\tQworst\t\t#Q>90\t\t#Q>90/#t\t"
"#e\t\ttau\t\t#Egood\t\t#Egood/#e\tCPU\n"); "#e\t\ttau\t\t#Egood\t\t#Egood/#e\tCPU\n");
if(m->empty()){ if(m->empty()){
fclose(statreport); fclose(statreport);
return; return;
...@@ -312,11 +315,11 @@ static void PrintMesh2dStatistics(GModel *m) ...@@ -312,11 +315,11 @@ static void PrintMesh2dStatistics(GModel *m)
fprintf(statreport,"\t%16s\t%d\t\t%d\t\t", m->getName().c_str(), numFaces, nUnmeshed); fprintf(statreport,"\t%16s\t%d\t\t%d\t\t", m->getName().c_str(), numFaces, nUnmeshed);
fprintf(statreport,"%d\t\t%8.7f\t%8.7f\t%8.7f\t%d\t\t%8.7f\t", fprintf(statreport,"%d\t\t%8.7f\t%8.7f\t%8.7f\t%d\t\t%8.7f\t",
nTotT, avg / (double)nTotT, best, worst, nTotGoodQuality, nTotT, avg / (double)nTotT, best, worst, nTotGoodQuality,
(double)nTotGoodQuality / nTotT); (double)nTotGoodQuality / nTotT);
fprintf(statreport,"%d\t\t%8.7f\t%d\t\t%8.7f\t%8.1f\n", fprintf(statreport,"%d\t\t%8.7f\t%d\t\t%8.7f\t%8.1f\n",
nTotE, exp(e_avg / (double)nTotE), nTotGoodLength, nTotE, exp(e_avg / (double)nTotE), nTotGoodLength,
(double)nTotGoodLength / nTotE, CTX::instance()->meshTimer[1]); (double)nTotGoodLength / nTotE, CTX::instance()->meshTimer[1]);
fclose(statreport); fclose(statreport);
} }
...@@ -355,19 +358,19 @@ static void Mesh2D(GModel *m) ...@@ -355,19 +358,19 @@ static void Mesh2D(GModel *m)
backgroundMesh::current()->unset(); backgroundMesh::current()->unset();
meshGFace mesher(true); meshGFace mesher(true);
mesher(temp[K]); mesher(temp[K]);
#if defined(HAVE_BFGS) #if defined(HAVE_BFGS)
if(CTX::instance()->mesh.optimizeLloyd){ if(CTX::instance()->mesh.optimizeLloyd){
if (temp[K]->geomType()==GEntity::CompoundSurface || if (temp[K]->geomType()==GEntity::CompoundSurface ||
temp[K]->geomType()==GEntity::Plane || temp[K]->geomType()==GEntity::RuledSurface) { temp[K]->geomType()==GEntity::Plane ||
temp[K]->geomType()==GEntity::RuledSurface) {
if (temp[K]->meshAttributes.method != MESH_TRANSFINITE && if (temp[K]->meshAttributes.method != MESH_TRANSFINITE &&
!temp[K]->meshAttributes.extrude) { !temp[K]->meshAttributes.extrude) {
smoothing smm(CTX::instance()->mesh.optimizeLloyd,6); smoothing smm(CTX::instance()->mesh.optimizeLloyd,6);
//m->writeMSH("beforeLLoyd.msh"); //m->writeMSH("beforeLLoyd.msh");
smm.optimize_face(temp[K]); smm.optimize_face(temp[K]);
int rec = ((CTX::instance()->mesh.recombineAll || int rec = ((CTX::instance()->mesh.recombineAll ||
temp[K]->meshAttributes.recombine) && temp[K]->meshAttributes.recombine) &&
!CTX::instance()->mesh.recombine3DAll); !CTX::instance()->mesh.recombine3DAll);
//m->writeMSH("afterLLoyd.msh"); //m->writeMSH("afterLLoyd.msh");
if (rec) recombineIntoQuads(temp[K]); if (rec) recombineIntoQuads(temp[K]);
//m->writeMSH("afterRecombine.msh"); //m->writeMSH("afterRecombine.msh");
...@@ -375,7 +378,6 @@ static void Mesh2D(GModel *m) ...@@ -375,7 +378,6 @@ static void Mesh2D(GModel *m)
} }
} }
#endif #endif
#if defined(_OPENMP) #if defined(_OPENMP)
#pragma omp critical #pragma omp critical
#endif #endif
...@@ -394,19 +396,19 @@ static void Mesh2D(GModel *m) ...@@ -394,19 +396,19 @@ static void Mesh2D(GModel *m)
backgroundMesh::current()->unset(); backgroundMesh::current()->unset();
meshGFace mesher(true); meshGFace mesher(true);
mesher(*it); mesher(*it);
#if defined(HAVE_BFGS) #if defined(HAVE_BFGS)
if(CTX::instance()->mesh.optimizeLloyd){ if(CTX::instance()->mesh.optimizeLloyd){
if ((*it)->geomType()==GEntity::CompoundSurface || if ((*it)->geomType()==GEntity::CompoundSurface ||
(*it)->geomType()==GEntity::Plane || (*it)->geomType()==GEntity::RuledSurface) { (*it)->geomType()==GEntity::Plane ||
(*it)->geomType()==GEntity::RuledSurface) {
if ((*it)->meshAttributes.method != MESH_TRANSFINITE && if ((*it)->meshAttributes.method != MESH_TRANSFINITE &&
!(*it)->meshAttributes.extrude) { !(*it)->meshAttributes.extrude) {
smoothing smm(CTX::instance()->mesh.optimizeLloyd,6); smoothing smm(CTX::instance()->mesh.optimizeLloyd,6);
//m->writeMSH("beforeLLoyd.msh"); //m->writeMSH("beforeLLoyd.msh");
smm.optimize_face(*it); smm.optimize_face(*it);
int rec = ((CTX::instance()->mesh.recombineAll || int rec = ((CTX::instance()->mesh.recombineAll ||
(*it)->meshAttributes.recombine) && (*it)->meshAttributes.recombine) &&
!CTX::instance()->mesh.recombine3DAll); !CTX::instance()->mesh.recombine3DAll);
//m->writeMSH("afterLLoyd.msh"); //m->writeMSH("afterLLoyd.msh");
if (rec) recombineIntoQuads(*it); if (rec) recombineIntoQuads(*it);
//m->writeMSH("afterRecombine.msh"); //m->writeMSH("afterRecombine.msh");
...@@ -414,7 +416,6 @@ static void Mesh2D(GModel *m) ...@@ -414,7 +416,6 @@ static void Mesh2D(GModel *m)
} }
} }
#endif #endif
nPending++; nPending++;
} }
if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D..."); if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
...@@ -426,6 +427,10 @@ static void Mesh2D(GModel *m) ...@@ -426,6 +427,10 @@ static void Mesh2D(GModel *m)
// collapseSmallEdges(*m); // collapseSmallEdges(*m);
// re-orient according to geometrical/user constraints
std::for_each(m->firstEdge(), m->lastEdge(), orientMeshGEdge());
std::for_each(m->firstFace(), m->lastFace(), orientMeshGFace());
double t2 = GetTimeInSeconds(); double t2 = GetTimeInSeconds();
CTX::instance()->meshTimer[1] = t2 - t1; CTX::instance()->meshTimer[1] = t2 - t1;
Msg::StatusBar(true, "Done meshing 2D (%g s)", CTX::instance()->meshTimer[1]); Msg::StatusBar(true, "Done meshing 2D (%g s)", CTX::instance()->meshTimer[1]);
...@@ -434,7 +439,7 @@ static void Mesh2D(GModel *m) ...@@ -434,7 +439,7 @@ static void Mesh2D(GModel *m)
} }
static void FindConnectedRegions(std::vector<GRegion*> &delaunay, static void FindConnectedRegions(std::vector<GRegion*> &delaunay,
std::vector<std::vector<GRegion*> > &connected) std::vector<std::vector<GRegion*> > &connected)
{ {
const unsigned int nbVolumes = delaunay.size(); const unsigned int nbVolumes = delaunay.size();
if (!nbVolumes)return; if (!nbVolumes)return;
...@@ -465,7 +470,7 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay, ...@@ -465,7 +470,7 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay,
delaunay=temp1; delaunay=temp1;
} }
Msg::Info("Delaunay Meshing %d volumes with %d connected components", Msg::Info("Delaunay Meshing %d volumes with %d connected components",
nbVolumes,connected.size()); nbVolumes,connected.size());
} }
static void Mesh3D(GModel *m) static void Mesh3D(GModel *m)
...@@ -567,6 +572,12 @@ static void Mesh3D(GModel *m) ...@@ -567,6 +572,12 @@ static void Mesh3D(GModel *m)
// Ensure that all volume Jacobians are positive // Ensure that all volume Jacobians are positive
m->setAllVolumesPositive(); m->setAllVolumesPositive();
// Ensure that all edge/surface meshes that could have been changed by the 3D
// algo (e.g. the Frontal), are re-oriented according to the geometrical/user
// constraints
std::for_each(m->firstEdge(), m->lastEdge(), orientMeshGEdge());
std::for_each(m->firstFace(), m->lastFace(), orientMeshGFace());
double t2 = Cpu(); double t2 = Cpu();
CTX::instance()->meshTimer[2] = t2 - t1; CTX::instance()->meshTimer[2] = t2 - t1;
Msg::StatusBar(true, "Done meshing 3D (%g s)", CTX::instance()->meshTimer[2]); Msg::StatusBar(true, "Done meshing 3D (%g s)", CTX::instance()->meshTimer[2]);
...@@ -668,17 +679,10 @@ void GenerateMesh(GModel *m, int ask) ...@@ -668,17 +679,10 @@ void GenerateMesh(GModel *m, int ask)
Mesh3D(m); Mesh3D(m);
} }
// Orient the line and surface meshes so that they match the orientation of
// the geometrical entities and/or the user orientation constraints
if(m->getMeshStatus() >= 1)
std::for_each(m->firstEdge(), m->lastEdge(), orientMeshGEdge());
if(m->getMeshStatus() >= 2)
std::for_each(m->firstFace(), m->lastFace(), orientMeshGFace());
// Optimize quality of 3D tet mesh // Optimize quality of 3D tet mesh
if(m->getMeshStatus() == 3){ if(m->getMeshStatus() == 3){
for(int i = 0; i < std::max(CTX::instance()->mesh.optimize, for(int i = 0; i < std::max(CTX::instance()->mesh.optimize,
CTX::instance()->mesh.optimizeNetgen); i++){ CTX::instance()->mesh.optimizeNetgen); i++){
if(CTX::instance()->mesh.optimize > i) OptimizeMesh(m); if(CTX::instance()->mesh.optimize > i) OptimizeMesh(m);
if(CTX::instance()->mesh.optimizeNetgen > i) OptimizeMeshNetgen(m); if(CTX::instance()->mesh.optimizeNetgen > i) OptimizeMeshNetgen(m);
} }
...@@ -696,7 +700,7 @@ void GenerateMesh(GModel *m, int ask) ...@@ -696,7 +700,7 @@ void GenerateMesh(GModel *m, int ask)
// Create high order elements // Create high order elements
if(m->getMeshStatus() && CTX::instance()->mesh.order > 1) if(m->getMeshStatus() && CTX::instance()->mesh.order > 1)
SetOrderN(m, CTX::instance()->mesh.order, CTX::instance()->mesh.secondOrderLinear, SetOrderN(m, CTX::instance()->mesh.order, CTX::instance()->mesh.secondOrderLinear,
CTX::instance()->mesh.secondOrderIncomplete); CTX::instance()->mesh.secondOrderIncomplete);
// Optimize high order elements // Optimize high order elements
if(CTX::instance()->mesh.hoOptimize){ if(CTX::instance()->mesh.hoOptimize){
...@@ -719,7 +723,7 @@ void GenerateMesh(GModel *m, int ask) ...@@ -719,7 +723,7 @@ void GenerateMesh(GModel *m, int ask)
} }
Msg::Info("%d vertices %d elements", Msg::Info("%d vertices %d elements",
m->getNumMeshVertices(), m->getNumMeshElements()); m->getNumMeshVertices(), m->getNumMeshElements());
Msg::PrintErrorCounter("Mesh generation error summary"); Msg::PrintErrorCounter("Mesh generation error summary");
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment