From 975652a97a01f5b6b00479953bab5c3b25fd5e25 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 16 Sep 2016 09:32:19 +0000
Subject: [PATCH] improvement of hybrid meshes
---
Geo/discreteDiskFace.cpp | 4 +--
Geo/discreteDiskFace.h | 16 ++++------
Geo/discreteFace.cpp | 13 +++++---
Mesh/CMakeLists.txt | 1 +
Mesh/Generator.cpp | 2 +-
Mesh/Generator.h | 4 +++
Mesh/meshGFaceDelaunayInsertion.cpp | 13 +++++---
Mesh/meshGRegion.cpp | 49 ++++++++++++++++++++++++++---
Mesh/meshGRegionRelocateVertex.cpp | 33 +++++++++++--------
Mesh/meshGRegionRelocateVertex.h | 7 +++--
Solver/dofManager.h | 2 +-
11 files changed, 100 insertions(+), 44 deletions(-)
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index b52d22a38c..5bb7361e62 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -1122,8 +1122,8 @@ static double computeDistanceLinePoint (MVertex *v1, MVertex *v2, MVertex *v){
}
*/
-static double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, MVertex *v)
-{
+inline double computeDistance (MVertex *v1, double d1, MVertex *v2, double d2, MVertex *v){
+
// o------------a
//
//
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 1719fe3d39..cad358b209 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -84,6 +84,7 @@ class triangulation {
{
double L = bord.rbegin()->first;
if (L == 0.0)return 1.e22;
+ if (maxD < 0) maxD = geodesicDistance();
// printf("%12.5E %12.5E\n",L,maxD);
return maxD / L;
}
@@ -136,7 +137,7 @@ class triangulation {
std::map<MVertex*,std::vector<MVertex*> > firstNode2Edge;
for (std::set<MEdge>::iterator ie = borderEdg.begin(); ie != borderEdg.end() ; ++ie) {
MEdge ed = *ie;
- std::vector<int> nT = ed2tri[ed];
+ const std::vector<int> &nT = ed2tri[ed];
MElement* t = tri[nT[0]];
std::vector<MVertex*> vecver;
@@ -163,7 +164,7 @@ class triangulation {
std::map<MVertex*,std::vector<MVertex*> >::iterator in = firstNode2Edge.begin();
MVertex* previous = in->first;
while(in != firstNode2Edge.end()) { // it didn't find it
- std::vector<MVertex*> myV = in->second;
+ const std::vector<MVertex*> &myV = in->second;
for(unsigned int i=0; i<myV.size(); i++){
loop.push_back(previous);
MVertex* current = myV[i];
@@ -189,7 +190,7 @@ class triangulation {
assignVert();
assignEd2tri();
assignBord();
- maxD = geodesicDistance();
+ maxD = -1.e22;
}
void print (char *name, int elementary) const {
@@ -203,8 +204,8 @@ class triangulation {
}
triangulation() : gf(0) {}
- triangulation(int id, std::vector<MElement*> input, GFace* gface)
- : idNum(id), tri(input), gf(gface) { assign(); }
+ triangulation(int id, const std::vector<MElement*> &input, GFace* gface)
+ : idNum(id), tri(input), gf(gface){ assign(); }
};
// triangles in the physical space xyz, with their parametric coordinates
@@ -247,12 +248,7 @@ class discreteDiskFace : public GFace {
void printAtlasMesh () ;
void printParamMesh () ;
-
-
std::vector<MElement*> discrete_triangles;
-
-
-
protected:
// a copy of the mesh that should not be destroyed
triangulation* initialTriangulation;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index ea1c672c1e..09a792deff 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -111,16 +111,16 @@ void discreteFace::secondDer(const SPoint2 ¶m,
void discreteFace::createGeometry()
{
-
checkAndFixOrientation();
#if defined(HAVE_SOLVER) && defined(HAVE_ANN)
int order = 1;
int nPart = 2;
double eta = 5/(2.*3.14);
-
if (!_atlas.empty())return;
+ double dtSplit = 0.0;
+
int id=1;
std::stack<triangulation*> toSplit;
std::vector<triangulation*> toParam;
@@ -134,7 +134,11 @@ void discreteFace::createGeometry()
std::vector<triangulation*> part;
triangulation* tosplit = toSplit.top();
toSplit.pop();
+
+ clock_t ts0 = clock();
split(tosplit,part,nPart);
+ clock_t ts1 = clock();
+ dtSplit += (double)(ts1-ts0)/CLOCKS_PER_SEC;
delete tosplit;
for(unsigned int i=0; i<part.size(); i++){
@@ -165,7 +169,6 @@ void discreteFace::createGeometry()
}
for(unsigned int i=0; i<toParam.size(); i++){
- std::vector<MElement*> mytri = toParam[i]->tri;
discreteDiskFace *df = new discreteDiskFace (this,toParam[i], order,(_CAD.empty() ? NULL : &_CAD));
df->printAtlasMesh();
df->replaceEdges(toParam[i]->my_GEdges);
@@ -228,7 +231,7 @@ void discreteFace::mesh(bool verbose)
_atlas[i]->mesh(verbose);/*
const char *name = "atlas%d";
char filename[256];
- sprintf(filename,name,i);
+ sprintf(filename,name,i);t0
_atlas[i]->printPhysicalMesh(filename);*/
}
@@ -461,7 +464,7 @@ void discreteFace::split(triangulation* trian,std::vector<triangulation*> &parti
for(int j=0; j<3; j++){// adjacency from edges
MEdge ed = current->getEdge(j);
if(trian->ed2tri[ed].size()>1){
- std::vector<int> oldnums = trian->ed2tri[ed];
+ const std::vector<int> &oldnums = trian->ed2tri[ed];
int on = trian->tri[oldnums[0]] == current ? oldnums[1] : oldnums[0];
if(check_todo.find(trian->tri[on])==check_todo.end() && el2part[trian->tri[on]]==p){
my_todo.push(trian->tri[on]);
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 7cf93c40fc..d870c95fbc 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -17,6 +17,7 @@ set(SRC
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
meshGFaceQuadrilateralize.cpp
meshGRegion.cpp
+ meshDiscreteRegion.cpp
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
meshGRegionLocalMeshMod.cpp meshGRegionMMG3D.cpp
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index b81bf497ce..4f14d5c3d2 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -802,7 +802,7 @@ static void Mesh3D(GModel *m)
// CTX::instance()->mesh.recombine3DConformity);
// 0: no pyramid, 1: single-step, 2: two-steps (conforming),
// true: fill non-conformities with trihedra
- RelocateVertices(gr);
+ RelocateVertices(gr, CTX::instance()->mesh.nbSmoothing);
// while(LaplaceSmoothing (gr)){
// }
nb_elements_recombination += post.get_nb_elements();
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index 52ec48a371..6d59607831 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -7,6 +7,8 @@
#define _GENERATOR_H_
class GModel;
+class GRegion;
+#include "fullMatrix.h"
void GetStatistics(double stat[50], double quality[4][100]=0);
void AdaptMesh(GModel *m);
@@ -17,5 +19,7 @@ void SmoothMesh(GModel *m);
void RefineMesh(GModel *m, bool linear, bool splitIntoQuads=false,
bool splitIntoHexas=false);
void RecombineMesh(GModel *m);
+GRegion * createTetrahedralMesh ( GModel *gm, fullMatrix<double> & pts, fullMatrix<int> &triangles ) ;
+ //GRegion * createTetrahedralMesh ( GModel *gm, unsigned int nbPts , double *pts, unsigned int nbTriangles, int *triangles );
#endif
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 259c4281b9..1e8b2e85fb 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1281,7 +1281,7 @@ double optimalPointFrontal (GFace *gf,
and the edge length
*/
-void optimalPointFrontalB (GFace *gf,
+bool optimalPointFrontalB (GFace *gf,
MTri3* worst,
int active_edge,
bidimMeshData & data,
@@ -1320,8 +1320,9 @@ void optimalPointFrontalB (GFace *gf,
if (gp.succeeded()){
newPoint[0] = gp.u();
newPoint[1] = gp.v();
- return ;
+ return true;
}
+ return false;
}
}
#endif
@@ -1333,7 +1334,7 @@ void optimalPointFrontalB (GFace *gf,
if (gp.succeeded()){
newPoint[0] = gp.u();
newPoint[1] = gp.v();
- return ;
+ return true;
}
}
}
@@ -1348,8 +1349,10 @@ void optimalPointFrontalB (GFace *gf,
}
else {
Msg::Debug("--- Non optimal point found -----------");
+ return false;
// Msg::Info("--- Non optimal point found -----------");
}
+ return true;
}
void bowyerWatsonFrontal(GFace *gf,
@@ -1406,8 +1409,8 @@ void bowyerWatsonFrontal(GFace *gf,
gf->mesh_vertices.size(), worst->getRadius());
double newPoint[2], metric[3];
//optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
- optimalPointFrontalB (gf,worst,active_edge,DATA,newPoint,metric);
- insertAPoint(gf, AllTris.end(), newPoint, metric, DATA, AllTris, &ActiveTris, worst);
+ if (optimalPointFrontalB (gf,worst,active_edge,DATA,newPoint,metric))
+ insertAPoint(gf, AllTris.end(), newPoint, metric, DATA, AllTris, &ActiveTris, worst);
}
/* if(ITER % 1== 0){
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 0eb02f7594..d04b788cfa 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -57,6 +57,41 @@ class splitQuadRecovery {
{
_data.insert(std::make_pair(ge, std::make_pair(v,f)));
}
+ void relocateVertices(GRegion *region, int niter) {
+ if(empty()) return ;
+ v2t_cont adj;
+ buildVertexToElement(region->tetrahedra, adj);
+ buildVertexToElement(region->pyramids, adj);
+ buildVertexToElement(region->prisms, adj);
+ buildVertexToElement(region->hexahedra, adj);
+
+ double minQual = 1;
+ double minQualOpti = 1;
+
+ std::list<GFace*> faces = region->faces();
+
+ for (int iter=0; iter < niter+2;iter++){
+ for (std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); ++it){
+ for (std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 =
+ _data.lower_bound(*it); it2 != _data.upper_bound(*it) ; ++it2){
+ const MFace &f = it2->second.second;
+ MVertex *v = it2->second.first;
+ MPyramid p (f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v);
+ minQual = std::min(minQual, fabs(p.minSICNShapeMeasure()));
+ std::vector<MElement*> e = adj[v];
+ e.push_back(&p);
+ v->setEntity (region);
+ double relax = std::min((double)(iter+1)/niter, 1.0);
+ // printf("%g (%d) --> ",e.size(),p.minSICNShapeMeasure());
+ _relocateVertexGolden( v, e, relax);
+ minQualOpti = std::min(minQualOpti, fabs(p.minSICNShapeMeasure()));
+ // printf("%g \n",p.minSICNShapeMeasure());
+ v->setEntity (*it);
+ }
+ }
+ }
+ printf("relocation improves %g --> %g\n", minQual, minQualOpti);
+ }
int buildPyramids(GModel *gm)
{
if(empty()) return 0;
@@ -226,9 +261,9 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
//// TEST
{
- std::vector<MVertex*>ALL;
- std::vector<MTetrahedron*> MESH;
- ALL.insert(ALL.begin(),allBoundingVertices.begin(),allBoundingVertices.end());
+ // std::vector<MVertex*>ALL;
+ // std::vector<MTetrahedron*> MESH;
+ // ALL.insert(ALL.begin(),allBoundingVertices.begin(),allBoundingVertices.end());
// delaunayMeshIn3D (ALL,MESH);
// exit(1);
}
@@ -561,6 +596,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
return;
}
TransferTetgenMesh(gr, in, out, numberedV);
+ sqr.relocateVertices(gr,3);
}
@@ -575,6 +611,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
+
// now do insertion of points
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL_DEL)
bowyerWatsonFrontalLayers(gr, false);
@@ -590,10 +627,12 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
insertVerticesInRegion(gr,2000000000,true);
}
}
-
+ // crete an initial mesh
+ printf("coucou1\n");
if (sqr.buildPyramids (gr->model())){
- RelocateVertices (regions);
+ RelocateVertices (regions, 3);
}
+ printf("coucou2\n");
}
#else
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index 3a4ac88484..cabedbfd06 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -25,7 +25,8 @@ static double objective_function (double xi, MVertex *ver,
if (lt[i]->getNumVertices () == 4)
minQual = std::min((lt[i]->minSICNShapeMeasure()), minQual);
else
- minQual = std::min(fabs(lt[i]->minSICNShapeMeasure()), minQual);
+ // minQual = std::min((lt[i]->specialQuality()), minQual);
+ minQual = std::min(fabs(lt[i]->minSICNShapeMeasure())*.2, minQual);
// minQual = std::min(lt[i]->minAnisotropyMeasure(), minQual);
}
ver->x() = x;
@@ -70,7 +71,7 @@ static int Stopping_Rule(double x0, double x1, double tol)
double Maximize_Quality_Golden_Section( MVertex *ver,
double xTarget, double yTarget, double zTarget,
- const std::vector<MElement*> < , double tol)
+ const std::vector<MElement*> < , double tol, double &q)
{
static const double lambda = 0.5 * (sqrt5 - 1.0);
@@ -104,6 +105,7 @@ double Maximize_Quality_Golden_Section( MVertex *ver,
}
}
// printf("finally : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
+ q = std::min(fx1,fx2);
return a;
}
@@ -155,8 +157,8 @@ double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
return a;
}
-static void _relocateVertex(MVertex *ver,
- const std::vector<MElement*> <, double tol)
+void _relocateVertexGolden(MVertex *ver,
+ const std::vector<MElement*> <, double relax , double tol)
{
if(ver->onWhat()->dim() != 3) return;
double x = 0, y=0, z=0;
@@ -174,14 +176,16 @@ static void _relocateVertex(MVertex *ver,
N += lt[i]->getNumVertices();
}
- double xi = Maximize_Quality_Golden_Section( ver, x/N, y/N, z/N, lt , tol);
+ double q;
+ double xi = relax * Maximize_Quality_Golden_Section( ver, x/N, y/N, z/N, lt , tol, q);
ver->x() = (1.-xi) * ver->x() + xi * x/N;
ver->y() = (1.-xi) * ver->y() + xi * y/N;
ver->z() = (1.-xi) * ver->z() + xi * z/N;
}
static double _relocateVertex(GFace* gf, MVertex *ver,
- const std::vector<MElement*> <, double tol) {
+ const std::vector<MElement*> <,
+ double tol) {
if(ver->onWhat()->dim() != 2) return 2.0;
SPoint2 p1(0,0);
@@ -239,22 +243,25 @@ void RelocateVertices (GFace* gf, int niter, double tol) {
}
-void RelocateVertices (GRegion* region, double tol) {
+void RelocateVertices (GRegion* region, int niter, double tol) {
v2t_cont adj;
buildVertexToElement(region->tetrahedra, adj);
buildVertexToElement(region->pyramids, adj);
buildVertexToElement(region->prisms, adj);
buildVertexToElement(region->hexahedra, adj);
- v2t_cont::iterator it = adj.begin();
- while (it != adj.end()){
- _relocateVertex( it->first, it->second, tol);
- ++it;
+ for (int i=0;i<niter+2;i++){
+ v2t_cont::iterator it = adj.begin();
+ double relax = std::min((double)(i+1)/niter, 1.0);
+ while (it != adj.end()){
+ _relocateVertexGolden( it->first, it->second, relax, tol);
+ ++it;
+ }
}
}
-void RelocateVertices (std::vector<GRegion*> ®ions, double tol) {
+void RelocateVertices (std::vector<GRegion*> ®ions, int niter, double tol) {
for(unsigned int k = 0; k < regions.size(); k++){
- RelocateVertices (regions[k], tol);
+ RelocateVertices (regions[k], niter, tol);
}
}
diff --git a/Mesh/meshGRegionRelocateVertex.h b/Mesh/meshGRegionRelocateVertex.h
index 56c8ac45f7..f267ef73d8 100644
--- a/Mesh/meshGRegionRelocateVertex.h
+++ b/Mesh/meshGRegionRelocateVertex.h
@@ -3,7 +3,10 @@
#include <vector>
class GRegion;
class GFace;
-void RelocateVertices (GRegion* region, double tol = 1.e-2);
-void RelocateVertices (std::vector<GRegion*> ®ions, double tol = 1.e-2);
+class MElement;
+void RelocateVertices (GRegion* region, int niter, double tol = 1.e-2);
+void RelocateVertices (std::vector<GRegion*> ®ions, int niter, double tol = 1.e-2);
void RelocateVertices (GFace*, int niter, double tol = 1.e-3);
+void _relocateVertexGolden(MVertex *ver, const std::vector<MElement*> <, double relax, double tol= 1.e-2);
+
#endif
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 1205cc3e59..6a10ab223f 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -357,7 +357,7 @@ class dofManager : public dofManagerBase{
insertInSparsityPatternLinConst(R, C);
}
}
-
+
virtual inline void sparsityDof(const std::vector<Dof> &keys){
for (unsigned int itR=0; itR< keys.size(); itR++){
for (unsigned int itC=0; itC<keys.size(); itC++){
--
GitLab