From 33a33c12b306c1b4b8c04e7fc480491135e19a90 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 9 Feb 2012 10:33:46 +0000
Subject: [PATCH] improved GFaceCompound and centerlines
---
Geo/GFaceCompound.cpp | 87 ++++++++++++------------
Geo/GModel.cpp | 2 +-
Geo/GModel.h | 2 +-
Mesh/CenterlineField.cpp | 139 ++++++++++++++++++++-------------------
Mesh/CenterlineField.h | 1 +
Mesh/meshMetric.cpp | 2 +-
6 files changed, 118 insertions(+), 115 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index ca7e1ad31c..1568dcba48 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -234,10 +234,10 @@ static SPoint3 myNeighVert(std::map<MVertex*,SPoint3> &coordinates,
if (v->onWhat()->dim() < 2) vN.insert(v);
}
}
- if (vN.size()!=3) {
- printf("ARGG more than 3 points on line =%d \n", vN.size());
- exit(1);
- }
+ // if (vN.size()!=3) {
+ // printf("ARGG more than 3 points on line =%d \n", vN.size());
+ // exit(1);
+ // }
double ucg = 0.0;
double vcg = 0.0;
@@ -446,7 +446,7 @@ void GFaceCompound::fillNeumannBCS() const
if (fillTris.size() > 0){
char name[256];
std::list<GFace*>::const_iterator itf = _compound.begin();
- sprintf(name, "fillTris-%d.pos", (*itf)->tag());
+ sprintf(name, "fillTris-%d.pos", tag());
FILE * ftri = fopen(name,"w");
fprintf(ftri,"View \"\"{\n");
for (std::list<MTriangle*>::iterator it2 = fillTris.begin();
@@ -570,13 +570,13 @@ bool GFaceCompound::checkOrientation(int iter, bool moveBoundaries) const
int iterMax = 15;
if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Debug("--- Flipping : applying cavity checks.");
+ if (iter == 0) Msg::Info("--- Flipping : applying cavity checks.");
Msg::Debug("--- Cavity Check - iter %d -",iter);
one2OneMap(moveBoundaries);
return checkOrientation(iter+1, moveBoundaries);
}
- else if (oriented && iter < iterMax){
- Msg::Debug("Parametrization is bijective (no flips)");
+ else if (iter > 0 && iter < iterMax){
+ Msg::Info("--- Flipping : no more flips (%d iter)", iter);
}
return oriented;
@@ -662,19 +662,27 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- bool hasOverlap = parametrize_conformal_spectral();
+ std::vector<MVertex *> vert;
+ bool oriented, hasOverlap;
+ hasOverlap = parametrize_conformal_spectral();
printStuff(11);
- if (hasOverlap || !checkOrientation(0, true) ){
- Msg::Warning("!!! Overlap or Flipping: parametrization switched to 'FE conformal' map");
- printStuff(22);
- hasOverlap = parametrize_conformal(0, NULL, NULL);
- }
- if (hasOverlap || !checkOrientation(0, true) ){
+ if (hasOverlap) oriented = checkOrientation(0);
+ else oriented = checkOrientation(0, true);
+ printStuff(22);
+ hasOverlap = checkOverlap(vert);
+ if ( !oriented || hasOverlap ){
+ Msg::Warning("!!! parametrization switched to 'FE conformal' map");
+ parametrize_conformal(0, NULL, NULL);
printStuff(33);
- Msg::Warning("$$$ Overlap or Flipping: parametrization switched to 'harmonic' map");
+ checkOrientation(0, true);
+ printStuff(44);
+ }
+ if (!checkOrientation(0) || checkOverlap(vert)){
+ printStuff(55);
+ Msg::Warning("$$$ parametrization switched to 'harmonic' map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
- }
+ }
}
// Radial-Basis Function parametrization
else if (_mapping == RBF){
@@ -708,7 +716,7 @@ bool GFaceCompound::parametrize() const
if (_mapping != RBF){
if (!checkOrientation(0)){
- Msg::Info("### Flipping: parametrization switched to convex combination map");
+ Msg::Info("### parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
@@ -1185,7 +1193,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
laplaceTerm laplace2(model(), 2, ONE);
crossConfTerm cross12(model(), 1, 2, ONE);
crossConfTerm cross21(model(), 2, 1, MONE);
- std::list<GFace*>::const_iterator it = _compound.begin();
+ std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
SElement se((*it)->triangles[i]);
@@ -1204,7 +1212,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
cross21.addToMatrix(myAssembler, &se);
}
- double epsilon = 1.e-7;
+ double epsilon = 1.e-6;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
if (std::find(_ordered.begin(), _ordered.end(), v) == _ordered.end() ){
@@ -1241,7 +1249,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
// FIXME: force non-hermitian. For some reason (roundoff errors?)
// on some machines and for some meshes slepc complains about bad IP
// norm otherwise
- eigenSolver eig(&myAssembler, "B" , "A", false);
+ eigenSolver eig(&myAssembler, "B" , "A", true);
bool converged = eig.solve(1, "largest");
if(converged) {
@@ -1270,15 +1278,8 @@ bool GFaceCompound::parametrize_conformal_spectral() const
return parametrize_conformal(0,NULL,NULL);
}
- std::vector<MVertex *> vert;
- bool hasOverlap = checkOverlap(vert);
- if (hasOverlap){
- Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
- printStuff(3);
- return hasOverlap = parametrize_conformal(0, vert[0], vert[1]);
- }
-
- return hasOverlap;
+ std::vector<MVertex *> vert;
+ return checkOverlap(vert);;
#endif
}
@@ -1370,7 +1371,7 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
// check for overlap and compute new mapping with new pinned
// vertices
std::vector<MVertex *> vert;
- bool hasOverlap = checkOverlap(vert);;
+ bool hasOverlap = checkOverlap(vert);
if ( hasOverlap && iter < 3){
Msg::Info("Loop FE conformal iter (%d) v1=%d v2=%d", iter,
vert[0]->getNum(), vert[1]->getNum());
@@ -1873,11 +1874,11 @@ bool GFaceCompound::checkTopology() const
tag(), nbSplit);
}
else if (_allowPartition == 0){
- Msg::Warning("The geometrical aspect ratio of your geometry is quite high. "
- "You should enable partitioning of the mesh by activating the "
- "automatic remeshin algorithm. Add 'Mesh.RemeshAlgorithm=1;' "
- "in your geo file or through the Fltk window (Options > Mesh > "
- "General)");
+ Msg::Warning("The geometrical aspect ratio of your geometry is quite high.\n "
+ "You should enable partitioning of the mesh by activating the\n"
+ "automatic remeshing algorithm. Add 'Mesh.RemeshAlgorithm=1;'\n "
+ "in your geo file or through the Fltk window (Options > Mesh >\n "
+ "General) \n");
}
}
else{
@@ -2096,13 +2097,13 @@ void GFaceCompound::printStuff(int iNewton) const
char name0[256], name1[256], name2[256], name3[256];
char name4[256], name5[256], name6[256];
char name7[256];
- sprintf(name0, "UVAREA-%d.pos", (*it)->tag());
- sprintf(name1, "UVX-%d_%d.pos", (*it)->tag(), iNewton);
- sprintf(name2, "UVY-%d_%d.pos", (*it)->tag(), iNewton);
- sprintf(name3, "UVZ-%d_%d.pos", (*it)->tag(), iNewton);
- sprintf(name4, "XYZU-%d_%d.pos", (*it)->tag(), iNewton);
- sprintf(name5, "XYZV-%d_%d.pos", (*it)->tag(), iNewton);
- sprintf(name6, "XYZC-%d.pos", (*it)->tag());
+ sprintf(name0, "UVAREA-%d.pos", tag()); //(*it)->tag()
+ sprintf(name1, "UVX-%d_%d.pos", tag(), iNewton);
+ sprintf(name2, "UVY-%d_%d.pos", tag(), iNewton);
+ sprintf(name3, "UVZ-%d_%d.pos", tag(), iNewton);
+ sprintf(name4, "XYZU-%d_%d.pos", tag(), iNewton);
+ sprintf(name5, "XYZV-%d_%d.pos", tag(), iNewton);
+ sprintf(name6, "XYZC-%d.pos", tag());
sprintf(name7, "UVM-%d.pos", (*it)->tag());
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 68b06968aa..7e470f7ed1 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1739,7 +1739,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
for (std::vector<discreteFace*>::iterator it = discFaces.begin();
it != discFaces.end(); it++)
(*it)->findEdges(map_edges);
-
+
// return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 3afcc3a34a..cedde36c71 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -384,7 +384,7 @@ class GModel
// Available algorithms:
// 1) Assume that the function is a levelset -> adapt using Coupez technique (technique = 1)
// parameters[0] = thickness of the interface (mandatory)
- // 2) Assume that the function is a physical quantity -> adapt using the Hessain (technique = 2)
+ // 2) Assume that the function is a physical quantity -> adapt using the Hessian (technique = 2)
// parameters[0] = N, the final number of elements
// 3) A variant of 1) by P. Frey (= Coupez + takes curvature function into account)
// parameters[0] = thickness of the interface (mandatory)
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 51afb7ad26..ef21db9c10 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -16,6 +16,7 @@
#include <string>
#include <fstream>
#include <sstream>
+#include <iostream>
#include "OS.h"
#include "GModel.h"
#include "MElement.h"
@@ -324,6 +325,7 @@ Centerline::Centerline(std::string fileName): kdtree(0), nodes(0){
printf("centerline filename =%s \n", fileName.c_str());
importFile(fileName);
buildKdTree();
+ nbPoints = 25;
update_needed = false;
is_cut = false;
@@ -336,12 +338,13 @@ Centerline::Centerline(): kdtree(0), nodes(0){
recombine = CTX::instance()->mesh.recombineAll;
fileName = "centerlines.vtk";//default
+ nbPoints = 25;
options["FileName"] = new FieldOptionString (fileName, "File name for the centerlines", &update_needed);
+ options["nbPoints"] = new FieldOptionInt(nbPoints, "Number of mesh elements in a circle");
callbacks["cutMesh"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::cutMesh, "Cut the initial mesh in different mesh partitions using the centerlines \n");
-
is_cut = false;
-
+
}
Centerline::~Centerline(){
@@ -373,7 +376,7 @@ void Centerline::importFile(std::string fileName){
}
mod = new GModel();
- mod->readVTK(fileName.c_str());
+ mod->load(fileName);
mod->removeDuplicateMeshVertices(1.e-8);
current->setAsCurrent();
@@ -624,7 +627,8 @@ void Centerline::buildKdTree(){
fprintf(f, "View \"\"{\n");
int nbPL = 3; //10 points per line
- int nbNodes = (lines.size()+1) + (nbPL*lines.size());
+ //int nbNodes = (lines.size()+1) + (nbPL*lines.size());
+ int nbNodes = (colorp.size()) + (nbPL*lines.size());
nodes = annAllocPts(nbNodes, 3);
int ind = 0;
@@ -684,7 +688,7 @@ void Centerline::createSplitCompounds(){
current->add(gec);
//gec->parametrize();
}
-
+
// Parametrize Compound surfaces
std::list<GEdge*> U0;
for (int i=0; i < NF; i++){
@@ -698,70 +702,47 @@ void Centerline::createSplitCompounds(){
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
+ gfc->addPhysicalEntity(i+1);
current->add(gfc);
//gfc->parametrize();
+ // for(unsigned int j = 0; j < pf->triangles.size(); ++j){
+ // MTriangle *t = pf->triangles[j];
+ // for(int k = 0; k < 3; k++){
+ // MVertex *v = t->getVertex(k);
+ // v->setEntity(pf);
+ // }
+ //}
}
}
void Centerline::cleanMesh(){
+ for (int i=0; i < NF; i++){
+ std::ostringstream oss;
+ oss << "new" << NF+i+1 ;
+ std::string name = oss.str();
+ current->setPhysicalName(name, 2, i+1);
+ }
+ Msg::Info("Writing new splitted mesh mySPLITMESH.msh");
+ current->writeMSH("mySPLITMESH.msh", 2.2, false, false);
+
if (!is_cut) return;
std::set<MVertex*> allNod;
std::list<GEdge*> U0;
- discreteFace *mySplitMesh = new discreteFace(current, 2*NF+1);
- current->add(mySplitMesh);
-
- //set centerline field
- // FieldManager *fields = current->getFields();
- // fields->reset();
- // field->setBackgroundField(this);
- // Msg::Info("*** Starting meshing 1D edges ...:");
- // printf("NE=%d \n", NE);
- // for (int i = 0; i < NE; i++){
- // printf("getting %d \n", NE+i+1);
- // GEdge *gec = current->getEdgeByTag(NE+i+1);
- // printf("edge (%p) =%d lines =%d \n", gec, gec->tag(), gec->lines.size());
- // meshGEdge mge;
- // mge(gec);
- // }
- // double t2 = Cpu();
- // Msg::Info("*** Meshing 1D edges done (%gs)", t2-t1);
-
- // Msg::Info("*** Starting mesh of surface ");
- // for (int i=0; i < NF; i++){
- // GFace *gfc = current->getFaceByTag(NF+i+1);
- // meshGFace mgf;
- // mgf(gfc);
- // for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
- // MTriangle *t = gfc->triangles[j];
- // std::vector<MVertex *> v(3);
- // for(int k = 0; k < 3; k++){
- // v[k] = t->getVertex(k);
- // allNod.insert(v[k]);
- // }
- // mySplitMesh->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
- // }
- // for(unsigned int j = 0; j < gfc->quadrangles.size(); ++j){
- // MQuadrangle *q = gfc->quadrangles[j];
- // std::vector<MVertex *> v(4);
- // for(int k = 0; k < 4; k++){
- // v[k] = q->getVertex(k);
- // allNod.insert(v[k]);
- // }
- // mySplitMesh->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
- // }
- // }
+ discreteFace * mySplitMesh;
+ mySplitMesh= new discreteFace(current, 2*NF+1);
+ current->add(mySplitMesh);
for (int i=0; i < NF; i++){
GFace *gfc = current->getFaceByTag(NF+i+1);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
std::vector<MVertex *> v(3);
for(int k = 0; k < 3; k++){
- v[k] = t->getVertex(k);
- allNod.insert(v[k]);
+ v[k] = t->getVertex(k);
+ allNod.insert(v[k]);
}
mySplitMesh->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
}
@@ -769,13 +750,13 @@ void Centerline::cleanMesh(){
MQuadrangle *q = gfc->quadrangles[j];
std::vector<MVertex *> v(4);
for(int k = 0; k < 4; k++){
- v[k] = q->getVertex(k);
- allNod.insert(v[k]);
+ v[k] = q->getVertex(k);
+ allNod.insert(v[k]);
}
mySplitMesh->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
}
}
-
+
//Removing discrete Vertices - Edges - Faces
for (int i=0; i < NV; i++){
GVertex *gv = current->getVertexByTag(i+1);
@@ -793,15 +774,13 @@ void Centerline::cleanMesh(){
current->remove(gf);
current->remove(gfc);
}
-
+
//Put new mesh in a new discreteFace
- for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it){
+ for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it)
mySplitMesh->mesh_vertices.push_back(*it);
- }
- current->createTopologyFromMesh();
-
mySplitMesh->meshStatistics.status = GFace::DONE;
-
+ current->createTopologyFromMesh();
+
}
void Centerline::createFaces(){
@@ -831,7 +810,7 @@ void Centerline::createFaces(){
it != touched.end(); ++it)
e2e.erase(*it);
}
- printf("%d faces created \n", faces.size());
+ printf("%d faces created \n", (int)faces.size());
//create discFaces
int numBef = current->getMaxElementaryNumber(2) + 1;
@@ -848,6 +827,7 @@ void Centerline::createFaces(){
for (int k= 0; k< 3; k++){
MVertex *v = t->getVertex(k);
myVertices.insert(v);
+ v->setEntity(f);
}
}
f->mesh_vertices.insert(f->mesh_vertices.begin(),
@@ -887,11 +867,11 @@ void Centerline::cutMesh(){
for(unsigned int i = 0; i < edges.size(); i++){
std::vector<MLine*> lines = edges[i].lines;
double L = edges[i].length;
- double R = edges[i].maxRad;
+ double R = edges[i].minRad;
double AR = L/R;
printf("*** branch =%d \n", i, AR);
- if( AR > 8.){
- int nbSplit = (int)ceil(AR/3.5);
+ if( AR > 4.5){
+ int nbSplit = (int)ceil(AR/4.0);
double li = L/nbSplit;
double lc = 0.0;
for (int j= 0; j < lines.size(); j++){
@@ -951,8 +931,8 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
allEdges.insert(triangles[i]->getEdge(j));
bool closedCut = false;
int step = 0;
- while (!closedCut && step < 20){
- double rad = 1.3*maxRad+0.1*step*maxRad;
+ while (!closedCut && step < 10){
+ double rad = 1.2*maxRad+0.1*step*maxRad;
std::map<MEdge,MVertex*,Less_Edge> cutEdges;
std::set<MVertex*> cutVertices;
std::vector<MTriangle*> newTris;
@@ -1035,16 +1015,37 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
int num_neighbours = 1;
kdtree->annkSearch(xyz, num_neighbours, index, dist);
double d = sqrt(dist[0]);
- double lc = 2*M_PI*d/30.0; //30 mesh elements along circle
+ double lc = 2*M_PI*d/nbPoints;
return lc;
}
void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
- printf("in operator xyz centerline BADDD \n");
- Msg::Error("This anisotropic operator of CenterlineField is not implemnted yet ");
- return;
+ //printf("in operator metric xyz centerline \n");
+ if (update_needed){
+ std::ifstream input;
+ input.open(fileName.c_str());
+ if(StatFile(fileName)) Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
+ importFile(fileName);
+ buildKdTree();
+ update_needed = false;
+ }
+
+ //double lc = operator()(x,y,z,ge);
+ //metr = SMetric3(1./(lc*lc));
+
+ double xyz[3] = {x,y,z };
+ ANNidxArray index2 = = new ANNidx[1];
+ ANNdistArray dist2 = new ANNdist[1];
+ int num_neighbours = 2;
+ kdtree->annkSearch(xyz, num_neighbours, index2, dist2);
+ double d = sqrt(dist2[0]);
+
+ delete[]index2;
+ delete[]dist2;
+
+ return;
}
void Centerline::printSplit() const{
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index bb0f70fad5..c88b91d274 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -60,6 +60,7 @@ class Centerline : public Field{
ANNidxArray index;
ANNdistArray dist;
std::string fileName;
+ int nbPoints;
double recombine;
int NF, NV, NE;
bool is_cut;
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index ce465eb00d..464b325127 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -192,7 +192,7 @@ void meshMetric::computeHessian( v2t_cont adj){
while (it != adj.end()) {
std::vector<MElement*> lt = it->second;
MVertex *ver = it->first;
- while (lt.size() < 9) increaseStencil(ver,adj,lt); //<7
+ while (lt.size() < 11) increaseStencil(ver,adj,lt); //<7
// if ( ver->onWhat()->dim() < _dim ){
// while (lt.size() < 12){
// increaseStencil(ver,adj,lt);
--
GitLab