diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 7daeb0c63028a09c9dc495650b2c05a37fb9ef12..440b8dceb82452cd40231bb098f943b9952776bb 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -12,7 +12,6 @@
#include "GModel.h"
#include "GModelIO_GEO.h"
#include "GModelIO_OCC.h"
-#include "Geo.h"
#include "GModelFactory.h"
#include "GFaceCompound.h"
#include "GEdgeCompound.h"
@@ -3105,60 +3104,18 @@ int GModel::readGEO(const std::string &name)
GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num)
{
- if (num ==-1) num = getMaxElementaryNumber(1) + 1;
+ if (num < 0) num = getMaxElementaryNumber(1) + 1;
GEdgeCompound *gec = new GEdgeCompound(this, num, edges);
add(gec);
-
- //create old geo format for the compound edge; necessary for boundary layers
- if(FindCurve(num)){
- Msg::Error("Curve %d already exists", num);
- }
- else{
- Curve *c = Create_Curve(num, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
- for(unsigned int i= 0; i < edges.size(); i++)
- c->compound.push_back(edges[i]->tag());
- List_T *points = Tree2List(getGEOInternals()->Points);
- GVertex *gvb = gec->getBeginVertex();
- GVertex *gve = gec->getEndVertex();
- int nb = 2 ;
- c->Control_Points = List_Create(nb, 1, sizeof(Vertex *));
- for(int i = 0; i < List_Nbr(points); i++) {
- Vertex *v;
- List_Read(points, i, &v);
- if (v->Num == gvb->tag()) {
- List_Add(c->Control_Points, &v);
- c->beg = v;
- }
- if (v->Num == gve->tag()) {
- List_Add(c->Control_Points, &v);
- c->end = v;
- }
- }
-
- End_Curve(c);
- Tree_Add(getGEOInternals()->Curves, &c);
- CreateReversedCurve(c);
-
- // c->Method = gec->meshAttributes.method;
- // c->nbPointsTransfinite = gec->meshAttributes.nbPointsTransfinite;
- // c->typeTransfinite = gec->meshAttributes.typeTransfinite ;
- // c->coeffTransfinite = gec->meshAttributes.coeffTransfinite ;
- // c->Extrude = gec->meshAttributes.extrude ;
- // c->ReverseMesh = gec->meshAttributes.reverseMesh ;
-
- }
-
return gec;
}
GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split, int num)
{
#if defined(HAVE_SOLVER)
- if (num ==-1) num = getMaxElementaryNumber(2) + 1;
-
+ if (num < 0) num = getMaxElementaryNumber(2) + 1;
std::list<GFace*> faces_comp(faces.begin(), faces.end());
std::list<GEdge*> U0;
-
GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
if (param == 1) typ = GFaceCompound::CONFORMAL_SPECTRAL;
if (param == 2) typ = GFaceCompound::RADIAL_BASIS;
@@ -3167,51 +3124,7 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split,
if (param == 5) typ = GFaceCompound::CONVEX_PLANE;
if (param == 6) typ = GFaceCompound::HARMONIC_SQUARE;
if (param == 7) typ = GFaceCompound::CONFORMAL_FE;
-
GFaceCompound *gfc = new GFaceCompound(this, num, faces_comp, U0, typ, split);
-
- //create old geo format for the compound face
- //necessary for boundary layers
- if(FindSurface(num)){
- Msg::Error("Surface %d already exists", num);
- }
- else{
- Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
- for(unsigned int i= 0; i < faces.size(); i++)
- s->compound.push_back(faces[i]->tag());
-
- std::list<GEdge*> edges = gfc->edges();
-
- //Replace edges by their compounds
- std::set<GEdge*> mySet;
- std::list<GEdge*>::iterator it = edges.begin();
- while(it != edges.end()){
- if((*it)->getCompound()){
- mySet.insert((*it)->getCompound());
- }
- else{
- mySet.insert(*it);
- }
- ++it;
- }
- edges.clear();
- edges.insert(edges.begin(), mySet.begin(), mySet.end());
-
- s->Generatrices = List_Create(edges.size(), 1, sizeof(Curve *));
- List_T *curves = Tree2List(_geo_internals->Curves);
- Curve *c;
- for(std::list<GEdge*>::iterator ite = edges.begin(); ite != edges.end(); ite++){
- for(int i = 0; i < List_Nbr(curves); i++) {
- List_Read(curves, i, &c);
- if (c->Num == (*ite)->tag()) {
- List_Add(s->Generatrices, &c);
- }
- }
- }
-
- Tree_Add(_geo_internals->Surfaces, &s);
- }
-
add(gfc);
return gfc;
#else
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 66040e93f823c264fc0fcda2751bbdda28afc245..62a90298b27a74fe0a51bf1f5aac20a378fbe98b 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -11,7 +11,6 @@ set(SRC
meshGEdge.cpp
meshGEdgeExtruded.cpp
meshGFace.cpp
- meshGFaceElliptic.cpp
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
meshGFaceBamg.cpp meshGFaceBDS.cpp meshGFaceDelaunayInsertion.cpp
meshGFaceLloyd.cpp meshGFaceOptimize.cpp
@@ -21,7 +20,7 @@ set(SRC
meshGRegionDelaunayInsertion.cpp meshGRegionTransfinite.cpp
meshGRegionExtruded.cpp meshGRegionCarveHole.cpp
meshGRegionLocalMeshMod.cpp meshGRegionMMG3D.cpp
-meshGRegionRelocateVertex.cpp
+ meshGRegionRelocateVertex.cpp
meshMetric.cpp
BackgroundMeshTools.cpp
BackgroundMesh.cpp
@@ -49,7 +48,6 @@ meshGRegionRelocateVertex.cpp
filterElements.cpp
yamakawa.cpp
Field.cpp
- CenterlineField.cpp
surfaceFiller.cpp
)
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
deleted file mode 100644
index 923cb9d645833498e5b3b7a5a4d55425ce4622a9..0000000000000000000000000000000000000000
--- a/Mesh/CenterlineField.cpp
+++ /dev/null
@@ -1,1336 +0,0 @@
-// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
-// Contributor(s):
-// Emilie Marchandise
-//
-
-#include "CenterlineField.h"
-
-#include <vector>
-#include <map>
-#include <set>
-#include <list>
-#include <algorithm>
-#include <string>
-#include <fstream>
-#include <sstream>
-#include <iostream>
-#include <math.h>
-#include "OS.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "MTriangle.h"
-#include "MVertex.h"
-#include "MLine.h"
-#include "StringUtils.h"
-#include "GEntity.h"
-#include "Field.h"
-#include "GFace.h"
-#include "discreteEdge.h"
-#include "discreteFace.h"
-#include "GEdgeCompound.h"
-#include "GFaceCompound.h"
-#include "BackgroundMeshTools.h"
-#include "meshGFace.h"
-#include "meshGEdge.h"
-#include "MQuadrangle.h"
-#include "Curvature.h"
-#include "MElement.h"
-#include "Context.h"
-#include "directions3D.h"
-#include "meshGRegion.h"
-
-#if defined(HAVE_ANN)
-#include "ANN/ANN.h"
-
-static void erase(std::vector<MLine*>& lines, MLine* l)
-{
- std::vector<MLine*>::iterator it = std::remove(lines.begin(), lines.end(), l);
- lines.erase(it, lines.end());
-}
-
-static double computeLength(std::vector<MLine*> lines)
-{
- double length= 0.0;
- for (unsigned int i = 0; i< lines.size(); i++){
- length += lines[i]->getLength();
- }
- return length;
-}
-
-static bool isClosed(std::set<MEdge, Less_Edge> &theCut)
-{
- std::multiset<MVertex*> boundV;
- std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
- for (; it!= theCut.end(); it++){
- boundV.insert(it->getVertex(0));
- boundV.insert(it->getVertex(1));
- }
- std::multiset<MVertex*>::iterator itb = boundV.begin();
- for ( ; itb != boundV.end(); ++itb){
- if (boundV.count(*itb) != 2) {
- return false;
- }
- }
- return true;
-}
-
-static void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE)
-{
- std::vector<MLine*> _m;
- std::list<MLine*> segments;
- std::map<MVertex*, MLine*> boundv;
- std::vector<int> _orientation;
-
- // store all lines in a list : segments
- for (unsigned int i = 0; i < lines.size(); i++){
- segments.push_back(lines[i]);
- }
-
- // find a lonely MLine
- for (std::list<MLine*>::iterator it = segments.begin();
- it != segments.end(); ++it){
- MVertex *vL = (*it)->getVertex(0);
- MVertex *vR = (*it)->getVertex(1);
- std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
- if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
- else boundv.erase(it1);
- std::map<MVertex*,MLine*>::iterator it2 = boundv.find(vR);
- if (it2 == boundv.end()) boundv.insert(std::make_pair(vR,*it));
- else boundv.erase(it2);
- }
-
- // find the first MLine and erase it from the list segments
- MLine *firstLine;
- if (boundv.size() == 2){ // non periodic
- MVertex *v = (boundv.begin())->first;
- if ( v == vB) firstLine = (boundv.begin())->second;
- else{
- MVertex *v = (boundv.rbegin())->first;
- if (v == vB) firstLine = (boundv.rbegin())->second;
- else{
- Msg::Error("begin vertex not found for branch");
- return;
- }
- }
- for (std::list<MLine*>::iterator it = segments.begin();
- it != segments.end(); ++it){
- if (*it == firstLine){
- segments.erase(it);
- break;
- }
- }
- }
- else{
- Msg::Error("line is wrong (it has %d end points)", boundv.size());
- }
-
- // loop over all segments to order segments and store it in the list _m
- _m.push_back(firstLine);
- _orientation.push_back(1);
- MVertex *first = _m[0]->getVertex(0);
- MVertex *last = _m[0]->getVertex(1);
- while (first != last){
- if (segments.empty())break;
- bool found = false;
- for (std::list<MLine*>::iterator it = segments.begin();
- it != segments.end(); ++it){
- MLine *e = *it;
- if (e->getVertex(0) == last){
- _m.push_back(e);
- segments.erase(it);
- _orientation.push_back(1);
- last = e->getVertex(1);
- found = true;
- break;
- }
- else if (e->getVertex(1) == last){
- _m.push_back(e);
- segments.erase(it);
- _orientation.push_back(0);
- last = e->getVertex(0);
- found = true;
- break;
- }
- }
- if (!found && _orientation[0]==1){ //reverse orientation of first Line
- if (_m.size() == 1 ){
- MVertex *temp = first;
- first = last;
- last = temp;
- _orientation[0] = 0;
- }
- else {
- Msg::Error("lines is wrong");
- return;
- }
- }
- }
-
- //lines is now a list of ordered MLines
- lines = _m;
-
- //special case reverse orientation
- if (lines.size() < 2) return;
- if (_orientation[0] && lines[0]->getVertex(1) != lines[1]->getVertex(1)
- && lines[0]->getVertex(1) != lines[1]->getVertex(0)){
- for (unsigned int i = 0; i < lines.size(); i++)
- _orientation[i] = !_orientation[i];
- }
-}
-
-static void recurConnectByMEdge(const MEdge &e,
- std::multimap<MEdge, MTriangle*, Less_Edge> &e2e,
- std::set<MTriangle*> &group,
- std::set<MEdge, Less_Edge> &touched,
- std::set<MEdge, Less_Edge> &theCut)
-{
- if (touched.find(e) != touched.end()) return;
- touched.insert(e);
- for (std::multimap <MEdge, MTriangle*, Less_Edge>::iterator it = e2e.lower_bound(e);
- it != e2e.upper_bound(e); ++it){
- group.insert(it->second);
- for (int i = 0; i < it->second->getNumEdges(); ++i){
- MEdge me = it->second->getEdge(i);
- if (theCut.find(me) != theCut.end()){
- touched.insert(me); //break;
- }
- else recurConnectByMEdge(me, e2e, group, touched, theCut);
- }
- }
-}
-
-static void cutTriangle(MTriangle *tri,
- std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
- std::vector<MVertex*> &cutVertices,
- std::vector<MTriangle*> &newTris,
- std::set<MEdge,Less_Edge> &newCut)
-{
- MVertex *c[3] = {0,0,0};
- for (int j=0;j<3;j++){
- MEdge ed = tri->getEdge(j);
- std::map<MEdge,MVertex*,Less_Edge> :: iterator it = cutEdges.find(ed);
- if (it != cutEdges.end()){
- c[j] = it->second;
-
- }
- }
- MVertex *old_v0 = tri->getVertex(0);
- MVertex *old_v1 = tri->getVertex(1);
- MVertex *old_v2 = tri->getVertex(2);
- if (c[0] && c[1]){
- newTris.push_back(new MTriangle (c[0],old_v1,c[1]));
- newTris.push_back(new MTriangle (old_v0,c[0],old_v2));
- newTris.push_back(new MTriangle (old_v2,c[0],c[1]));
- newCut.insert(MEdge(c[0],c[1]));
- }
- else if (c[0] && c[2]){
- newTris.push_back(new MTriangle (old_v0,c[0],c[2]));
- newTris.push_back(new MTriangle (c[0],old_v1,old_v2));
- newTris.push_back(new MTriangle (old_v2,c[2],c[0]));
- newCut.insert(MEdge(c[0],c[2]));
- }
- else if (c[1] && c[2]){
- newTris.push_back(new MTriangle (old_v2,c[2],c[1]));
- newTris.push_back(new MTriangle (old_v0,old_v1,c[2]));
- newTris.push_back(new MTriangle (c[2],old_v1,c[1]));
- newCut.insert(MEdge(c[1],c[2]));
- }
- else if (c[0]){
- newTris.push_back(new MTriangle (old_v0,c[0],old_v2));
- newTris.push_back(new MTriangle (old_v2,c[0],old_v1));
- if (std::find(cutVertices.begin(), cutVertices.end(), old_v0) != cutVertices.end()){
- newCut.insert(MEdge(c[0],old_v0));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v1) != cutVertices.end()) {
- newCut.insert(MEdge(c[0],old_v1));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v2) != cutVertices.end()){
- newCut.insert(MEdge(c[0],old_v2));
- }
- }
- else if (c[1]){
- newTris.push_back(new MTriangle (old_v1,c[1],old_v0));
- newTris.push_back(new MTriangle (old_v0,c[1],old_v2));
- if (std::find(cutVertices.begin(), cutVertices.end(), old_v0) != cutVertices.end()){
- newCut.insert(MEdge(c[1],old_v0));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v1) != cutVertices.end()) {
- newCut.insert(MEdge(old_v1, c[1]));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v2) != cutVertices.end()){
- newCut.insert(MEdge(c[1],old_v2));
- }
- }
- else if (c[2]){
- newTris.push_back(new MTriangle (old_v0,old_v1, c[2]));
- newTris.push_back(new MTriangle (old_v1,old_v2, c[2]));
- if (std::find(cutVertices.begin(), cutVertices.end(), old_v0) != cutVertices.end()){
- newCut.insert(MEdge(c[2],old_v0));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v1) != cutVertices.end()) {
- newCut.insert(MEdge(c[2], old_v1));
- }
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v2) != cutVertices.end()){
- newCut.insert(MEdge(c[2], old_v2));
- }
- }
- else {
- newTris.push_back(tri);
- if (std::find(cutVertices.begin(), cutVertices.end(), old_v0) != cutVertices.end() &&
- std::find(cutVertices.begin(), cutVertices.end(), old_v1) != cutVertices.end())
- newCut.insert(MEdge(old_v0,old_v1));
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v1) != cutVertices.end() &&
- std::find(cutVertices.begin(), cutVertices.end(), old_v2) != cutVertices.end())
- newCut.insert(MEdge(old_v1,old_v2));
- else if (std::find(cutVertices.begin(), cutVertices.end(), old_v2) != cutVertices.end() &&
- std::find(cutVertices.begin(), cutVertices.end(), old_v0) != cutVertices.end())
- newCut.insert(MEdge(old_v2,old_v0));
- }
-}
-
-Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0)
-{
- recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll);
- nbPoints = 25;
- hLayer = 0.3;
- hSecondLayer = 0.3;
- nbElemLayer = 3;
- nbElemSecondLayer = 0;
- is_cut = 0;
- is_closed = 0;
- is_extruded = 0;
-
- importFile(fileName);
- buildKdTree();
- update_needed = false;
-}
-
-Centerline::Centerline(): kdtree(0), kdtreeR(0)
-{
-
- recombine = (CTX::instance()->mesh.recombineAll) || (CTX::instance()->mesh.recombine3DAll);
- fileName = "centerlines.vtk";//default
- nbPoints = 25;
- hLayer = 0.3;
- hSecondLayer = 0.3;
- nbElemLayer = 3;
- nbElemSecondLayer = 0;
- is_cut = 0;
- is_closed = 0;
- is_extruded = 0;
-
- options["closeVolume"] = new FieldOptionInt
- (is_closed, "Action: Create In/Outlet planar faces");
- options["extrudeWall"] = new FieldOptionInt
- (is_extruded, "Action: Extrude wall");
- options["reMesh"] = new FieldOptionInt
- (is_cut, "Action: Cut the initial mesh in different mesh partitions using the "
- "centerlines");
- callbacks["run"] = new FieldCallbackGeneric<Centerline>
- (this, &Centerline::run, "Run actions (closeVolume, extrudeWall, cutMesh) \n");
-
- options["FileName"] = new FieldOptionString
- (fileName, "File name for the centerlines", &update_needed);
- options["nbPoints"] = new FieldOptionInt
- (nbPoints, "Number of mesh elements in a circle");
- options["nbElemLayer"] = new FieldOptionInt
- (nbElemLayer, "Number of mesh elements the extruded layer");
- options["hLayer"] = new FieldOptionDouble
- (hLayer, "Thickness (% of radius) of the extruded layer");
- options["nbElemSecondLayer"] = new FieldOptionInt
- (nbElemSecondLayer, "Number of mesh elements the second extruded layer");
- options["hSecondLayer"] = new FieldOptionDouble
- (hSecondLayer, "Thickness (% of radius) of the second extruded layer");
-}
-
-Centerline::~Centerline()
-{
- if (mod) delete mod;
- if(kdtree){
- ANNpointArray nodes = kdtree->thePoints();
- if(nodes) annDeallocPts(nodes);
- delete kdtree;
- }
- if(kdtreeR){
- ANNpointArray nodesR = kdtreeR->thePoints();
- if(nodesR) annDeallocPts(nodesR);
- delete kdtreeR;
- }
-}
-
-void Centerline::importFile(std::string fileName)
-{
- current = GModel::current();
- std::vector<GFace*> currentFaces(current->firstFace(), current->lastFace());
- for (unsigned int i = 0; i < currentFaces.size(); i++){
- GFace *gf = currentFaces[i];
- if (gf->geomType() == GEntity::DiscreteSurface){
- for(unsigned int j = 0; j < gf->triangles.size(); j++)
- triangles.push_back(gf->triangles[j]);
- if (is_cut){
- gf->triangles.clear();
- gf->deleteVertexArrays();
- current->remove(gf);
- }
- }
- }
-
- if(triangles.empty()){
- Msg::Error("Current GModel has no triangles ...");
- return;
- }
-
- mod = new GModel();
- mod->load(fileName);
- mod->removeDuplicateMeshVertices(1.e-8);
- current->setAsCurrent();
- current->setVisibility(1);
-
- int maxN = 0.0;
- std::vector<GEdge*> modEdges(mod->firstEdge(), mod->lastEdge());
- MVertex *vin = modEdges[0]->lines[0]->getVertex(0);
- ptin = SPoint3(vin->x(), vin->y(), vin->z());
- for (unsigned int i = 0; i < modEdges.size(); i++){
- GEdge *ge = modEdges[i];
- for(unsigned int j = 0; j < ge->lines.size(); j++){
- MLine *l = ge->lines[j];
- MVertex *v0 = l->getVertex(0);
- MVertex *v1 = l->getVertex(1);
- std::map<MVertex*, int>::iterator it0 = colorp.find(v0);
- std::map<MVertex*, int>::iterator it1 = colorp.find(v1);
- if (it0 == colorp.end() || it1 == colorp.end()){
- lines.push_back(l);
- colorl.insert(std::make_pair(l, ge->tag()));
- maxN = std::max(maxN, ge->tag());
- }
- if (it0 == colorp.end()) colorp.insert(std::make_pair(v0, ge->tag()));
- if (it1 == colorp.end()) colorp.insert(std::make_pair(v1, ge->tag()));
- }
- }
-
- createBranches(maxN);
-}
-
-void Centerline::createBranches(int maxN)
-{
- //sort colored lines and create edges
- std::vector<std::vector<MLine*> > color_edges;
- color_edges.resize(maxN);
- std::multiset<MVertex*> allV;
- std::map<MLine*, int>::iterator itl = colorl.begin();
- while (itl != colorl.end()){
- int color = itl->second;
- MLine* l = itl->first;
- allV.insert(l->getVertex(0));
- allV.insert(l->getVertex(1));
- color_edges[color-1].push_back(l);
- itl++;
- }
-
- //detect junctions
- std::multiset<MVertex*>::iterator it = allV.begin();
- for ( ; it != allV.end(); ++it){
- if (allV.count(*it) != 2) {
- junctions.insert(*it);
- }
- }
-
- //split edges
- int tag = 0;
- for(unsigned int i = 0; i < color_edges.size(); ++i){
- std::vector<MLine*> lines = color_edges[i];
- while (!lines.empty()) {
- std::vector<MLine*> myLines;
- std::vector<MLine*>::iterator itl = lines.begin();
- MVertex *vB = (*itl)->getVertex(0);
- MVertex *vE = (*itl)->getVertex(1);
- myLines.push_back(*itl);
- erase(lines, *itl);
- itl = lines.begin();
- while ( !( junctions.find(vE) != junctions.end() &&
- junctions.find(vB) != junctions.end()) ) {
- MVertex *v1 = (*itl)->getVertex(0);
- MVertex *v2 = (*itl)->getVertex(1);
- bool goVE = (junctions.find(vE) == junctions.end()) ? true : false ;
- bool goVB = (junctions.find(vB) == junctions.end()) ? true : false;
- if (v1 == vE && goVE){
- myLines.push_back(*itl);
- erase(lines, *itl);
- itl = lines.begin();
- vE = v2;
- }
- else if ( v2 == vE && goVE){
- myLines.push_back(*itl);
- erase(lines, *itl);
- itl = lines.begin();
- vE = v1;
- }
- else if ( v1 == vB && goVB){
- myLines.push_back(*itl);
- erase(lines, *itl);
- itl = lines.begin();
- vB = v2;
- }
- else if ( v2 == vB && goVB){
- myLines.push_back(*itl);
- erase(lines, *itl);
- itl = lines.begin();
- vB = v1;
- }
- else itl++;
- }
- if (vB == vE) {
- Msg::Error("Begin and end points branch are the same \n");
- break;
- }
- orderMLines(myLines, vB, vE);
- std::vector<Branch> children;
- Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE,
- children, 1.e6, 0.0};
- edges.push_back(newBranch) ;
- }
- }
-
- Msg::Info("Centerline: in/outlets =%d branches =%d ",
- (int)color_edges.size()+1, (int)edges.size());
-
- //create children
- for(unsigned int i = 0; i < edges.size(); ++i) {
- MVertex *vE = edges[i].vE;
- std::vector<Branch> myChildren;
- for (std::vector<Branch>::iterator it = edges.begin(); it != edges.end(); ++it){
- Branch myBranch = *it;
- if (myBranch.vB == vE) myChildren.push_back(myBranch);
- }
- edges[i].children = myChildren;
- }
-
- //compute radius
- distanceToSurface();
- computeRadii();
-
- //print for debug
- printSplit();
-
-}
-
-void Centerline::distanceToSurface()
-{
- Msg::Info("Centerline: computing distance to surface mesh ");
-
- //COMPUTE WITH REVERSE ANN TREE (SURFACE POINTS IN TREE)
- std::set<MVertex*> allVS;
- for(unsigned int j = 0; j < triangles.size(); j++)
- for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
- int nbSNodes = allVS.size();
- ANNpointArray nodesR = annAllocPts(nbSNodes, 3);
- vertices.resize(nbSNodes);
- int ind = 0;
- std::set<MVertex*>::iterator itp = allVS.begin();
- while (itp != allVS.end()){
- MVertex *v = *itp;
- nodesR[ind][0] = v->x();
- nodesR[ind][1] = v->y();
- nodesR[ind][2] = v->z();
- vertices[ind] = v;
- itp++; ind++;
- }
- kdtreeR = new ANNkd_tree(nodesR, nbSNodes, 3);
-
- for(unsigned int i = 0; i < lines.size(); i++){
- MLine *l = lines[i];
- MVertex *v1 = l->getVertex(0);
- MVertex *v2 = l->getVertex(1);
- double midp[3] = {0.5*(v1->x()+v2->x()), 0.5*(v1->y()+v1->y()),0.5*(v1->z()+v2->z())};
- ANNidx index[1];
- ANNdist dist[1];
- kdtreeR->annkSearch(midp, 1, index, dist);
- double minRad = sqrt(dist[0]);
- radiusl.insert(std::make_pair(lines[i], minRad));
- }
-
-}
-
-void Centerline::computeRadii()
-{
- for(unsigned int i = 0; i < edges.size(); ++i) {
- std::vector<MLine*> lines = edges[i].lines;
- for(unsigned int j = 0; j < lines.size(); ++j) {
- MLine *l = lines[j];
- std::map<MLine*,double>::iterator itr = radiusl.find(l);
- if (itr != radiusl.end()){
- edges[i].minRad = std::min(itr->second, edges[i].minRad);
- edges[i].maxRad = std::max(itr->second, edges[i].maxRad);
- }
- else printf("ARGG line not found \n");
- }
- }
-
-}
-
-void Centerline::buildKdTree()
-{
- FILE *f = Fopen("myPOINTS.pos","w");
-
- if(f) fprintf(f, "View \"\"{\n");
-
- int nbPL = 3; //10 points per line
- //int nbNodes = (lines.size()+1) + (nbPL*lines.size());
- int nbNodes = (colorp.size()) + (nbPL*lines.size());
-
- ANNpointArray nodes = annAllocPts(nbNodes, 3);
- int ind = 0;
- std::map<MVertex*, int>::iterator itp = colorp.begin();
- while (itp != colorp.end()){
- MVertex *v = itp->first;
- nodes[ind][0] = v->x();
- nodes[ind][1] = v->y();
- nodes[ind][2] = v->z();
- itp++; ind++;
- }
- for(unsigned int k = 0; k < lines.size(); ++k){
- MVertex *v0 = lines[k]->getVertex(0);
- MVertex *v1 = lines[k]->getVertex(1);
- SVector3 P0(v0->x(),v0->y(), v0->z());
- SVector3 P1(v1->x(),v1->y(), v1->z());
- for (int j= 1; j < nbPL+1; j++){
- double inc = (double)j/(double)(nbPL+1);
- SVector3 Pj = P0+inc*(P1-P0);
- nodes[ind][0] = Pj.x();
- nodes[ind][1] = Pj.y();
- nodes[ind][2] = Pj.z();
- ind++;
- }
- }
-
- kdtree = new ANNkd_tree(nodes, nbNodes, 3);
-
- if(f){
- for(int i = 0; i < nbNodes; ++i){
- fprintf(f, "SP(%g,%g,%g){%g};\n",
- nodes[i][0], nodes[i][1],nodes[i][2],1.0);
- }
- fprintf(f,"};\n");
- fclose(f);
- }
-}
-
-void Centerline::createSplitCompounds()
-{
- //number of discrete vertices, edges, faces and regions for the mesh
- NV = current->getMaxElementaryNumber(0);
- NE = current->getMaxElementaryNumber(1);
- NF = current->getMaxElementaryNumber(2);
- NR = current->getMaxElementaryNumber(3);
-
- // Remesh new faces (Compound Lines and Compound Surfaces)
- Msg::Info("Centerline: creating split compounds ...");
-
- //Parametrize Compound Lines
- for (int i=0; i < NE; i++){
- std::vector<GEdge*>e_compound;
- GEdge *pe = current->getEdgeByTag(i+1);//current edge
- e_compound.push_back(pe);
- int num_gec = NE+i+1;
- Msg::Info("Create Compound Line (%d) = %d discrete edge",
- num_gec, pe->tag());
- GEdge *gec = current->addCompoundEdge(e_compound,num_gec);
- if (CTX::instance()->mesh.algo2d != ALGO_2D_BAMG){
- gec->meshAttributes.method = MESH_TRANSFINITE;
- gec->meshAttributes.nbPointsTransfinite = nbPoints+1;
- gec->meshAttributes.typeTransfinite = 0;
- gec->meshAttributes.coeffTransfinite = 1.0;
- }
- }
-
- // Parametrize Compound surfaces
- std::list<GEdge*> U0;
- for (int i=0; i < NF; i++){
- std::vector<GFace*> f_compound;
- GFace *pf = current->getFaceByTag(i+1);//current face
- f_compound.push_back(pf);
- int num_gfc = NF+i+1;
- Msg::Info("Create Compound Surface (%d) = %d discrete face",
- num_gfc, pf->tag());
-
- //1=conf_spectral 4=convex_circle, 7=conf_fe
- GFace *gfc = current->addCompoundFace(f_compound, 7, 0, num_gfc);
-
- gfc->meshAttributes.recombine = recombine;
- gfc->addPhysicalEntity(1);
- current->setPhysicalName("wall", 2, 1);//tag 1
-
- }
-}
-
-void Centerline::createFaces()
-{
- std::vector<std::vector<MTriangle*> > faces;
-
- std::multimap<MEdge, MTriangle*, Less_Edge> e2e;
- for(unsigned int i = 0; i < triangles.size(); ++i)
- for(int j = 0; j < 3; j++)
- e2e.insert(std::make_pair(triangles[i]->getEdge(j), triangles[i]));
- while(!e2e.empty()){
- std::set<MTriangle*> group;
- std::set<MEdge, Less_Edge> touched;
- group.clear();
- touched.clear();
- std::multimap<MEdge, MTriangle*, Less_Edge>::iterator ite = e2e.begin();
- MEdge me = ite->first;
- while (theCut.find(me) != theCut.end()){
- ite++;
- me = ite->first;
- }
- recurConnectByMEdge(me,e2e, group, touched, theCut);
- std::vector<MTriangle*> temp;
- temp.insert(temp.begin(), group.begin(), group.end());
- faces.push_back(temp);
- for(std::set<MEdge, Less_Edge>::iterator it = touched.begin();
- it != touched.end(); ++it)
- e2e.erase(*it);
- }
- Msg::Info("Centerline: action (cutMesh) has cut surface mesh in %d faces ",
- (int)faces.size());
-
- //create discFaces
- for(unsigned int i = 0; i < faces.size(); ++i){
- int numF = current->getMaxElementaryNumber(2) + 1;
- discreteFace *f = new discreteFace(current, numF);
- current->add(f);
- discFaces.push_back(f);
- std::set<MVertex*> myVertices;
- std::vector<MTriangle*> myFace = faces[i];
- for(unsigned int j= 0; j< myFace.size(); j++){
- MTriangle *t = myFace[j];
- f->triangles.push_back(t);
- 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(),
- myVertices.begin(), myVertices.end());
- }
-}
-
-void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges)
-{
- current->setFactory("Gmsh");
- std::vector<std::vector<GFace *> > myFaceLoops;
- std::vector<GFace *> myFaces;
- for (unsigned int i = 0; i< boundEdges.size(); i++){
- std::vector<std::vector<GEdge *> > myEdgeLoops;
- std::vector<GEdge *> myEdges;
- GEdge * gec;
- if(is_cut) gec = current->getEdgeByTag(NE+boundEdges[i]->tag());
- else gec = current->getEdgeByTag(boundEdges[i]->tag());
- myEdges.push_back(gec);
- myEdgeLoops.push_back(myEdges);
- GFace *newFace = current->addPlanarFace(myEdgeLoops);
- if (gin==boundEdges[i]) {
- newFace->addPhysicalEntity(2);
- current->setPhysicalName("inlet", 2, 2);//tag 2
- }
- else{
- newFace->addPhysicalEntity(3);
- current->setPhysicalName("outlets", 2, 3);//tag 3
- }
- myFaces.push_back(newFace);
- }
-
- Msg::Info("Centerline: action (closeVolume) has created %d in/out planar faces ",
- (int)boundEdges.size());
-
- for (int i = 0; i < NF; i++){
- GFace * gf;
- if(is_cut) gf = current->getFaceByTag(NF+i+1);
- else gf = current->getFaceByTag(i+1);
- myFaces.push_back(gf);
- }
- myFaceLoops.push_back(myFaces);
- GRegion *reg = current->addVolume(myFaceLoops);
- reg->addPhysicalEntity(reg->tag());
- current->setPhysicalName("lumenVolume", 3, reg->tag());
-
- Msg::Info("Centerline: action (closeVolume) has created volume %d ", reg->tag());
-
-}
-
-void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundEdges)
-{
- Msg::Info("Centerline: extrude boundary layer wall (%d, %g%%R) ", nbElemLayer, hLayer);
-
- //orient extrude direction outward
- int dir = 0;
- MElement *e = current->getFaceByTag(1)->getMeshElement(0);
- SVector3 ne = e->getFace(0).normal();
- SVector3 ps(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z());
- double xyz[3] = {ps.x(), ps.y(), ps.z()};
- ANNidx index[1];
- ANNdist dist[1];
- kdtree->annkSearch(xyz, 1, index, dist);
- ANNpointArray nodes = kdtree->thePoints();
- SVector3 pc(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
- SVector3 nc = ps-pc;
- if (dot(ne,nc) < 0) dir = 1;
- if (dir == 1 && hLayer > 0 ) hLayer *= -1.0;
-
- //int shift = 0;
- //if(is_cut) shift = NE;
- for (int i= 0; i< NF; i++){
- GFace *gfc ;
- if (is_cut) gfc = current->getFaceByTag(NF+i+1);
- else gfc = current->getFaceByTag(i+1);
- current->setFactory("Gmsh");
- //view -5 to scale hLayer by radius in BoundaryLayers.cpp
- std::vector<GEntity*> extrudedE = current->extrudeBoundaryLayer
- (gfc, nbElemLayer, hLayer, dir, -5);
- GFace *eFace = (GFace*) extrudedE[0];
- eFace->addPhysicalEntity(5);
- current->setPhysicalName("outerWall", 2, 5);//dim 2 tag 5
- GRegion *eRegion = (GRegion*) extrudedE[1];
- eRegion->addPhysicalEntity(6);
- current->setPhysicalName("wallVolume", 3, 6);//dim 3 tag 6
-
- //if double extruded layer
- if (nbElemSecondLayer > 0){
- std::vector<GEntity*> extrudedESec = current->extrudeBoundaryLayer
- (eFace, nbElemSecondLayer, hSecondLayer, dir, -5);
- GFace *eFaceSec = (GFace*) extrudedESec[0];
- eFaceSec->addPhysicalEntity(9); //tag 9
- current->setPhysicalName("outerSecondWall", 2, 9);//dim 2 tag 9
- GRegion *eRegionSec = (GRegion*) extrudedESec[1];
- eRegionSec->addPhysicalEntity(10); //tag 10
- current->setPhysicalName("wallVolume", 3, 10);//dim 3 tag 10
- }
- //end double extrusion
-
- for (unsigned int j = 2; j < extrudedE.size(); j++){
- GFace *elFace = (GFace*) extrudedE[j];
- std::list<GEdge*> l_edges = elFace->edges();
- for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
- GEdge *myEdge = *it;
- if (is_cut) myEdge = current->getEdgeByTag((*it)->tag()-NE);
- if( std::find(boundEdges.begin(), boundEdges.end(), myEdge) != boundEdges.end() ){
- if (myEdge==gin){
- elFace->addPhysicalEntity(7);
- current->setPhysicalName("inletRing", 2, 7);//tag 7
- }
- else{
- elFace->addPhysicalEntity(8);
- current->setPhysicalName("outletRings", 2, 8);//tag 8
- }
- }
- }
- }
- }
-
-}
-
-void Centerline::run()
-{
- double t1 = Cpu();
- if (update_needed){
- std::ifstream input;
- //std::string pattern = FixRelativePath(fileName, "./");
- //Msg::StatusBar(true, "Reading TEST '%s'...", pattern.c_str());
- //input.open(pattern.c_str());
- 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;
- }
-
- if (is_cut) cutMesh();
- else{
- GFace *gf = current->getFaceByTag(1);
- gf->addPhysicalEntity(1);
- current->setPhysicalName("wall", 2, 1);//tag 1
- current->createTopologyFromMesh();
- NV = current->getMaxElementaryNumber(0);
- NE = current->getMaxElementaryNumber(1);
- NF = current->getMaxElementaryNumber(2);
- NR = current->getMaxElementaryNumber(3);
- }
-
- //identify the boundary edges by looping over all discreteFaces
- std::vector<GEdge*> boundEdges;
- double dist_inlet = 1.e6;
- GEdge *gin = NULL;
- for (int i= 0; i< NF; i++){
- GFace *gf = current->getFaceByTag(i+1);
- std::list<GEdge*> l_edges = gf->edges();
- for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
- std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(),
- boundEdges.end(), *it);
- if (ite != boundEdges.end()) boundEdges.erase(ite);
- else boundEdges.push_back(*it);
- GVertex *gv = (*it)->getBeginVertex();
- SPoint3 pt(gv->x(), gv->y(), gv->z());
- double dist = pt.distance(ptin);
- if(dist < dist_inlet){
- dist_inlet = dist;
- gin = *it;
- }
- }
- }
-
- if (is_closed) createClosedVolume(gin, boundEdges);
- if (is_extruded) extrudeBoundaryLayerWall(gin, boundEdges);
-
- double t2 = Cpu();
- Msg::Info("Centerline operators computed in %g (s) ",t2-t1);
-}
-
-void Centerline::cutMesh()
-{
- Msg::Info("Centerline: action (cutMesh) splits surface mesh (%d tris) using %s ",
- triangles.size(), fileName.c_str());
-
- //splitMesh
- for(unsigned int i = 0; i < edges.size(); i++){
- std::vector<MLine*> lines = edges[i].lines;
- double L = edges[i].length;
- double D = 2.*edges[i].minRad; //(edges[i].minRad+edges[i].maxRad);
- double AR = L/D;
- // printf("*** Centerline branch %d (AR=%.1f) \n", edges[i].tag, AR);
-
- int nbSplit = (int)ceil(AR/2 + 1.1); //AR/2 + 0.9
- if( nbSplit > 1 ){
- //printf("->> cut branch in %d parts \n", nbSplit);
- double li = L/nbSplit;
- double lc = 0.0;
- for (unsigned int j= 0; j < lines.size(); j++){
- lc += lines[j]->getLength();
- if (lc > li && nbSplit > 1) {
- MVertex *v1 = lines[j]->getVertex(0);
- MVertex *v2 = lines[j]->getVertex(1);
- SVector3 pt(v1->x(), v1->y(), v1->z());
- SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
- std::map<MLine*,double>::iterator itr = radiusl.find(lines[j]);
- cutByDisk(pt, dir, itr->second);
- nbSplit--;
- lc = 0.0;
- }
- }
- }
- if(edges[i].children.size() > 0.0 && AR > 1.0){
- MVertex *v1 = lines[lines.size()-1]->getVertex(1);//end vertex
- MVertex *v2;
- if(AR < 1.5) v2 = lines[0]->getVertex(0);
- else if (lines.size() > 4) v2 = lines[lines.size()-4]->getVertex(0);
- else v2 = lines[lines.size()-1]->getVertex(0);
- SVector3 pt(v1->x(), v1->y(), v1->z());
- SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
- //printf("-->> cut branch at bifurcation \n");
- std::map<MLine*,double>::iterator itr = radiusl.find(lines[lines.size()-1]);
- //bool cutted =
- cutByDisk(pt, dir, itr->second);
- // if(!cutted){
- // int l = lines.size()-1-lines.size()/(4*nbSplit); //chech this!
- // v1 = lines[l]->getVertex(1);
- // v2 = lines[l]->getVertex(0);
- // pt = SVector3(v1->x(), v1->y(), v1->z());
- // dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
- // printf("-->> cut bifurcation NEW \n");
- // itr = radiusl.find(lines[l]);
- // cutted = cutByDisk(pt, dir, itr->second);
- // }
- }
- }
-
- //create discreteFaces
- createFaces();
- current->createTopologyFromFaces(discFaces);
- current->exportDiscreteGEOInternals();
-
- //write
- Msg::Info("Centerline: writing splitted mesh 'myPARTS.msh'");
- current->writeMSH("myPARTS.msh", 2.2, false, false);
-
- //create compounds
- createSplitCompounds();
-
- Msg::Info("Done splitting mesh by centerlines");
-}
-
-bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad)
-{
- double a = NORM.x();
- double b = NORM.y();
- double c = NORM.z();
- double d = -a * PT.x() - b * PT.y() - c * PT.z();
-
- int maxStep = 20;
- const double EPS = 0.007;
-
- //std::set<MEdge,Less_Edge> allEdges;
- std::vector<MEdge> allEdges;
- for(unsigned int i = 0; i < triangles.size(); i++){
- for ( unsigned int j= 0; j < 3; j++){
- allEdges.push_back(triangles[i]->getEdge(j));
- //allEdges.insert(triangles[i]->getEdge(j));
- }
- }
- std::unique(allEdges.begin(), allEdges.end());
-
- bool closedCut = false;
- int step = 0;
- while (!closedCut && step < maxStep){
- double rad = 1.1*maxRad+0.05*step*maxRad;
- std::map<MEdge,MVertex*,Less_Edge> cutEdges;
- std::vector<MVertex*> cutVertices;
- std::vector<MTriangle*> newTris;
- std::set<MEdge,Less_Edge> newCut;
- cutEdges.clear();
- cutVertices.clear();
- newTris.clear();
- newCut.clear();
- // for (std::set<MEdge,Less_Edge>::iterator it = allEdges.begin();
- // it != allEdges.end() ; ++it){
- // MEdge me = *it;
- for (unsigned int j = 0; j < allEdges.size(); j++){
- MEdge me = allEdges[j];
- SVector3 P1(me.getVertex(0)->x(),me.getVertex(0)->y(), me.getVertex(0)->z());
- SVector3 P2(me.getVertex(1)->x(),me.getVertex(1)->y(), me.getVertex(1)->z());
- double V1 = a * P1.x() + b * P1.y() + c * P1.z() + d;
- double V2 = a * P2.x() + b * P2.y() + c * P2.z() + d;
- bool inters = (V1*V2<=0.0) ? true: false;
- bool inDisk = ((norm(P1-PT) < rad ) || (norm(P2-PT) < rad)) ? true : false;
- double rdist = -V1/(V2-V1);
- if (inters && rdist > EPS && rdist < 1.-EPS){
- SVector3 PZ = P1+rdist*(P2-P1);
- if (inDisk){
- MVertex *newv = new MVertex (PZ.x(), PZ.y(), PZ.z());
- cutEdges.insert(std::make_pair(me,newv));
- }
- }
- else if (inters && rdist <= EPS && inDisk )
- cutVertices.push_back(me.getVertex(0));
- else if (inters && rdist >= 1.-EPS && inDisk)
- cutVertices.push_back(me.getVertex(1));
- }
- for(unsigned int i = 0; i < triangles.size(); i++){
- cutTriangle(triangles[i], cutEdges,cutVertices, newTris, newCut);
- }
- if (isClosed(newCut)) {
- triangles.clear();
- triangles = newTris;
- theCut.insert(newCut.begin(),newCut.end());
- break;
- }
- else {
- step++;
- //if (step == maxStep) {printf("no closed cut %d \n", (int)newCut.size()); };
- // // triangles = newTris;
- // // theCut.insert(newCut.begin(),newCut.end());
- // char name[256];
- // sprintf(name, "myCUT-%d.pos", step);
- // FILE * f2 = Fopen(name,"w");
- // if(f2){
- // fprintf(f2, "View \"\"{\n");
- // std::set<MEdge,Less_Edge>::iterator itp = newCut.begin();
- // while (itp != newCut.end()){
- // MEdge l = *itp;
- // fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- // l.getVertex(0)->x(), l.getVertex(0)->y(), l.getVertex(0)->z(),
- // l.getVertex(1)->x(), l.getVertex(1)->y(), l.getVertex(1)->z(),
- // 1.0,1.0);
- // itp++;
- // }
- // fprintf(f2,"};\n");
- // fclose(f2);
- // }
- }
- }
-
-
- if (step < maxStep){
- //printf("cutByDisk OK step =%d \n", step);
- return true;
- }
- else {
- //printf("cutByDisk not succeeded \n");
- return false;
- }
-
-}
-
-double Centerline::operator() (double x, double y, double z, GEntity *ge)
-{
-
- 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 xyz[3] = {x,y,z};
- //take xyz = closest point on boundary in case we are on the planar in/out faces
- //or in the volume
- bool isCompound = false;
- if(ge){
- if (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) isCompound = true;
- std::list<GFace*> cFaces;
- if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
- if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
- (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
- const int num_neighbours = 1;
- ANNidx index[num_neighbours];
- ANNdist dist[num_neighbours];
- kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
- ANNpointArray nodesR = kdtreeR->thePoints();
- xyz[0] = nodesR[index[0]][0];
- xyz[1] = nodesR[index[0]][1];
- xyz[2] = nodesR[index[0]][2];
- }
- }
-
- const int num_neighbours = 1;
- ANNidx index[num_neighbours];
- ANNdist dist[num_neighbours];
- kdtree->annkSearch(xyz, num_neighbours, index, dist);
- double rad = sqrt(dist[0]);
-
- //double cmax, cmin;
- //SVector3 dirMax,dirMin;
- //cmax = ge->curvatures(SPoint2(u, v),&dirMax, &dirMin, &cmax,&cmin);
- //cmax = ge->curvatureMax(SPoint2(u,v));
- //double radC = 1./cmax;
-
- double lc = 2*M_PI*rad/nbPoints;
-
- if(!ge) { return rad;}
- else return lc;
-
-}
-
-void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
-{
-
- 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;
- }
-
-
- //take xyz = closest point on boundary in case we are on
- //the planar IN/OUT FACES or in VOLUME
- double xyz[3] = {x,y,z};
- bool onTubularSurface = true;
- double ds = 0.0;
- bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ?
- true : false;
- bool onInOutlets = (ge->geomType() == GEntity::Plane) ? true: false;
- std::list<GFace*> cFaces;
- if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
- if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
- (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
- onTubularSurface = false;
- }
-
- ANNidx index[1];
- ANNdist dist[1];
- kdtreeR->annkSearch(xyz, 1, index, dist);
- if (! onTubularSurface){
- ANNpointArray nodesR = kdtreeR->thePoints();
- ds = sqrt(dist[0]);
- xyz[0] = nodesR[index[0]][0];
- xyz[1] = nodesR[index[0]][1];
- xyz[2] = nodesR[index[0]][2];
- }
-
- ANNidx index2[2];
- ANNdist dist2[2];
- kdtree->annkSearch(xyz, 2, index2, dist2);
- double radMax = sqrt(dist2[0]);
- ANNpointArray nodes = kdtree->thePoints();
- SVector3 p0(nodes[index2[0]][0], nodes[index2[0]][1], nodes[index2[0]][2]);
- SVector3 p1(nodes[index2[1]][0], nodes[index2[1]][1], nodes[index2[1]][2]);
-
- //dir_a = direction along the centerline
- //dir_n = normal direction of the disk
- //dir_t = tangential direction of the disk
- SVector3 dir_a = p1-p0; dir_a.normalize();
- SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); dir_n.normalize();
- SVector3 dir_cross = crossprod(dir_a,dir_n); dir_cross.normalize();
- // if (ge->dim()==3){
- // printf("coucou dim ==3 !!!!!!!!!!!!!!! \n");
- // SVector3 d1,d2,d3;
- // computeCrossField(x,y,z,d1,d2,d3);
- // exit(1);
- // }
-
- //find discrete curvature at closest vertex
- Curvature& curvature = Curvature::getInstance();
- if( !Curvature::valueAlreadyComputed() ) {
- Msg::Info("Need to compute discrete curvature");
- Curvature::typeOfCurvature type = Curvature::RUSIN;
- curvature.computeCurvature(current, type);
- }
- double curv, cMin, cMax;
- SVector3 dMin, dMax;
- int isAbs = 1.0;
- curvature.vertexNodalValuesAndDirections(vertices[index[0]],&dMax, &dMin, &cMax, &cMin, isAbs);
- curvature.vertexNodalValues(vertices[index[0]], curv, 1);
- if (cMin == 0) cMin =1.e-12;
- if (cMax == 0) cMax =1.e-12;
- double rhoMin = 1./cMin;
- double rhoMax = 1./cMax;
- //double signMin = (rhoMin > 0.0) ? -1.0: 1.0;
- //double signMax = (rhoMax > 0.0) ? -1.0: 1.0;
-
- //user-defined parameters
- //define h_n, h_t1, and h_t2
- double thickness = radMax/3.;
- double h_far = radMax/5.;
- double beta = (ds <= thickness) ? 1.2 : 2.1; //CTX::instance()->mesh.smoothRatio;
- double ddist = (ds <= thickness) ? ds: thickness;
-
- double h_n_0 = thickness/20.;
- double h_n = std::min( (h_n_0+ds*log(beta)), h_far);
-
- double betaMin = 10.0;
- double betaMax = 3.1;
- double oneOverD2_min = 1./(2.*rhoMin*rhoMin*(betaMin*betaMin-1)) *
- (sqrt(1+ (4.*rhoMin*rhoMin*(betaMin*betaMin-1))/(h_n*h_n))-1.);
- double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(betaMax*betaMax-1)) *
- (sqrt(1+ (4.*rhoMax*rhoMax*(betaMax*betaMax-1))/(h_n*h_n))-1.);
- double h_t1_0 = sqrt(1./oneOverD2_min);
- double h_t2_0 = sqrt(1./oneOverD2_max);
- //double h_t1 = h_t1_0*(rhoMin+signMin*ddist)/rhoMin ;
- //double h_t2 = h_t2_0*(rhoMax+signMax*ddist)/rhoMax ;
- double h_t1 = std::min( (h_t1_0+(ddist*log(beta))), radMax);
- double h_t2 = std::min( (h_t2_0+(ddist*log(beta))), h_far);
-
- double dCenter = radMax-ds;
- double h_a_0 = 0.5*radMax;
- double h_a = h_a_0 - (h_a_0-h_t1_0)/(radMax)*dCenter;
-
- //length between min and max
- double lcMin = ((2 * M_PI *radMax) /( 50*nbPoints )); //CTX::instance()->mesh.lcMin;
- double lcMax = lcMin*2000.; //CTX::instance()->mesh.lcMax;
- h_n = std::max(h_n, lcMin); h_n = std::min(h_n, lcMax);
- h_t1 = std::max(h_t1, lcMin); h_t1 = std::min(h_t1, lcMax);
- h_t2 = std::max(h_t2, lcMin); h_t2 = std::min(h_t2, lcMax);
-
- //curvature metric
- SMetric3 curvMetric, curvMetric1, curvMetric2;
- SMetric3 centMetric1, centMetric2, centMetric;
- if (onInOutlets){
- metr = buildMetricTangentToCurve(dir_n,h_n,h_t2);
- }
- else {
- //on surface and in volume boundary layer
- if ( ds < thickness ){
- metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2);
- }
- //in volume
- else {
- //curvMetric = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, h_n, h_t1, h_t2);
- metr = SMetric3( 1./(h_a*h_a), 1./(h_n*h_n), 1./(h_n*h_n), dir_a, dir_n, dir_cross);
-
- //metr = intersection_conserveM1_bis(metr, curvMetric);
- //metr = intersection_conserveM1(metr,curvMetric);
- //metr = intersection_conserve_mostaniso(metr, curvMetric);
- //metr = intersection(metr,curvMetric);
- }
- }
-
- return;
-
-}
-
-SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dirMax,
- double cmin, double cmax,
- double h_n, double h_t1, double h_t2)
-{
-
- SVector3 dirNorm = crossprod(dirMax,dirMin);
- SMetric3 curvMetric (1./(h_t1*h_t1),1./(h_t2*h_t2),1./(h_n*h_n), dirMin, dirMax, dirNorm);
-
- return curvMetric;
-}
-
-void Centerline::computeCrossField(double x,double y,double z,
- SVector3 &d1, SVector3 &d2, SVector3 &d3)
-{
-
- //Le code suivant permet d'obtenir les trois vecteurs orthogonaux en un point.
- // int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
- // Matrix m2;
- // if(NumSmooth)
- // m2 = Frame_field::findCross(x,y,z);
- // else
- // m2 = Frame_field::search(x,y,z);
-
-}
-
-void Centerline::printSplit() const
-{
- FILE * f = Fopen("mySPLIT.pos","w");
- if(f){
- fprintf(f, "View \"\"{\n");
- for(unsigned int i = 0; i < edges.size(); ++i){
- std::vector<MLine*> lines = edges[i].lines;
- for(unsigned int k = 0; k < lines.size(); ++k){
- MLine *l = lines[k];
- fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
- l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
- (double)edges[i].tag, (double)edges[i].tag);
- }
- }
- fprintf(f,"};\n");
- fclose(f);
- }
-
- // FILE * f3 = Fopen("myJUNCTIONS.pos","w");
- // if(f3){
- // fprintf(f3, "View \"\"{\n");
- // std::set<MVertex*>::const_iterator itj = junctions.begin();
- // while (itj != junctions.end()){
- // MVertex *v = *itj;
- // fprintf(f3, "SP(%g,%g,%g){%g};\n",
- // v->x(), v->y(), v->z(),
- // (double)v->getNum());
- // itj++;
- // }
- // fprintf(f3,"};\n");
- // fclose(f3);
- // }
-
- FILE * f4 = Fopen("myRADII.pos","w");
- if(f4){
- fprintf(f4, "View \"\"{\n");
- for(unsigned int i = 0; i < lines.size(); ++i){
- MLine *l = lines[i];
- std::map<MLine*,double>::const_iterator itc = radiusl.find(l);
- fprintf(f4, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
- l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
- itc->second,itc->second);
- }
- fprintf(f4,"};\n");
- fclose(f4);
- }
-}
-
-#endif
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
deleted file mode 100644
index 28c3a1976cbdb202941ebb600a2c3c947ee4a4d0..0000000000000000000000000000000000000000
--- a/Mesh/CenterlineField.h
+++ /dev/null
@@ -1,201 +0,0 @@
-// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-//
-// Contributor(s):
-// Emilie Marchandise
-
-#ifndef _CENTERLINEFIELD_H_
-#define _CENTERLINEFIELD_H_
-
-#include <vector>
-#include <map>
-#include <set>
-#include <string>
-#include "Field.h"
-#include "MEdge.h"
-
-#include "meshGFaceDelaunayInsertion.h"
-class GModel;
-class GFace;
-class MLine;
-class MVertex;
-class GEntity;
-class MTriangle;
-class discreteEdge;
-class discreteFace;
-class MElement;
-class SPoint3;
-
-// A branch of a 1D tree
-struct Branch{
- int tag;
- std::vector<MLine*> lines;
- double length;
- MVertex *vB;
- MVertex *vE;
- std::vector<Branch> children;
- double minRad;
- double maxRad;
-};
-
-#if defined(HAVE_ANN)
-class ANNkd_tree;
-
-// This class takes as input A 1D mesh which is the centerline
-// of a tubular 2D surface mesh
-// It computes a mesh size field function of the distance to the centerlines
-// and a cross field that is the direction of the centerline
-// It splits the tubular structure in many mesh partitions
-// along planes perpendicuar to the centerlines
-
-class Centerline : public Field{
-
- protected:
- GModel *current; //current GModel
- GModel *mod; //centerline GModel
- GModel *split; //split GModel
- ANNkd_tree *kdtree, *kdtreeR;
- std::string fileName;
- int nbPoints;
- double recombine;
- int NF, NV, NE, NR;
- int is_cut, is_closed, is_extruded;
- double hLayer;
- double hSecondLayer;
- int nbElemLayer;
- int nbElemSecondLayer;
-
- //inlet point
- SPoint3 ptin;
- //all (unique) lines of centerlines
- std::vector<MLine*> lines;
- //the stuctured tree of the centerlines
- std::vector<Branch> edges;
- //the radius of the surface mesh at a given line
- std::map<MLine*,double> radiusl;
- //the junctions of the tree
- std::set<MVertex*> junctions;
- //some colors (int) for all points and lines
- std::map<MVertex*,int> colorp;
- std::map<MLine*,int> colorl;
-
- //the tubular surface mesh
- std::vector<MTriangle*> triangles;
- std::vector<MVertex*> vertices;
-
- //the lines cut of the tubular mesh by planes
- std::set<MEdge,Less_Edge> theCut;
- std::set<MVertex*> theCutV;
-
- //discrete edes and faces created by the cut
- std::vector<discreteEdge*> discEdges;
- std::vector<discreteFace*> discFaces;
-
- public:
- Centerline(std::string fileName);
- Centerline();
- ~Centerline();
-
- virtual bool isotropic () const {return false;}
- virtual const char *getName()
- {
- return "centerline Field";
- }
- virtual std::string getDescription()
- {
- return "The value of this field is the distance to the centerline.\n\n"
-" You should specify a fileName that contains the centerline."
-" The centerline of a surface can be obtained with the open source software vmtk (http://www.vmtk.org/)"
-" using the following script:\n\n"
-"vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
- }
-
- //isotropic operator for mesh size field function of distance to centerline
- double operator() (double x, double y, double z, GEntity *ge=0);
- //anisotropic operator
- void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
-
- void computeCrossField(double x,double y,double z,
- SVector3 &d1, SVector3 &d2, SVector3 &d3);
-
- //import the 1D mesh of the centerlines (in vtk format)
- //and fill the vector of lines
- void importFile(std::string fileName);
-
- //refine the 1D mesh to have many points on the 1D centerlines
- //for annKDTree distance computations
- void buildKdTree();
-
- //Creates the branch structure (topology, connectivity) from the
- //vector of lines
- void createBranches(int maxN);
-
- //Computes for the Branches the min and maxRadius of the tubular structure
- //this function needs the current GModel
- void computeRadii();
-
- //Computes for each MLine the minRadius
- void distanceToSurface();
-
- //actions
- void run();
-
- // Cut the mesh in different parts of small aspect ratio
- void cutMesh();
- //Create In and Outlet Planar Faces
- void createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges);
- //extrude outer wall
- void extrudeBoundaryLayerWall(GEdge *gin, std::vector<GEdge*> boundEdges);
-
- // Cut the tubular structure with a disk
- // perpendicular to the tubular structure
- bool cutByDisk(SVector3 &pt, SVector3 &dir, double &maxRad);
-
- //create discrete faces
- void createFaces();
- void createSplitCompounds();
-
- //Print for debugging
- void printSplit() const;
-
- SMetric3 metricBasedOnSurfaceCurvature(SVector3 dMin, SVector3 dMax, double cMin, double cMax,
- double lc_n, double lc_t1, double lc_t2);
-
-};
-#else
-class Centerline : public Field{
-
- public:
- Centerline(std::string fileName){ Msg::Error("Gmsh has to be compiled with ANN support to use CenterlineFields");}
- Centerline(){ Msg::Error("Gmsh has to be compiled with ANN support to use CenterlineFields");}
- ~Centerline();
-
- virtual bool isotropic () const {return false;}
- virtual const char *getName()
- {
- return "centerline Field";
- }
- virtual std::string getDescription()
- {
- return "The value of this field is the distance to the centerline.\n\n"
-" You should specify a fileName that contains the centerline."
-" The centerline of a surface can be obtained with the open source software vmtk (http://www.vmtk.org/)"
-" using the following script:\n\n"
-"vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
- }
-
- //isotropic operator for mesh size field function of distance to centerline
- double operator() (double x, double y, double z, GEntity *ge=0);
- //anisotropic operator
- void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
-
- //temporary operator where v1, v2 and v3 are three orthonormal directions
- void operator()(double x,double y,double z,SVector3& v1,SVector3& v2,SVector3& v3,GEntity* ge=0);
-
-};
-
-#endif
-
-#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 6289d212687175ba46e74d87b646a9c164d30aad..79ee1ef00f10cf49b2c3c5c967887c07c946507a 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -26,7 +26,6 @@
#include "Numeric.h"
#include "mathEvaluator.h"
#include "BackgroundMeshTools.h"
-#include "CenterlineField.h"
#include "STensor3.h"
#include "meshMetric.h"
#if defined(HAVE_POST)
@@ -2702,7 +2701,6 @@ FieldManager::FieldManager()
map_type_name["Threshold"] = new FieldFactoryT<ThresholdField>();
#if defined(HAVE_ANN)
map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
- map_type_name["Centerline"] = new FieldFactoryT<Centerline>();
#endif
map_type_name["Box"] = new FieldFactoryT<BoxField>();
map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 8f4042452deb0364259f1ea37b567eb21fe5430e..19e3e05b367c21a01c9e856cfd380fc3147f012c 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -29,14 +29,12 @@
#include "HighOrder.h"
#include "Generator.h"
#include "meshGFaceLloyd.h"
-#include "CenterlineField.h"
#include "GFaceCompound.h"
#include "Field.h"
#include "Options.h"
#include "simple3D.h"
#include "yamakawa.h"
#include "meshGRegionRelocateVertex.h"
-
#include "pointInsertion.h"
#if defined(HAVE_OPTHOM)
@@ -70,7 +68,7 @@ public:
faces.insert(std::make_pair(MFace((*itf)->triangles[i]->getVertex(0),(*itf)->triangles[i]->getVertex(1),(*itf)->triangles[i]->getVertex(2)),*itf));
}
}
- Msg::Info ("Searching for %d embedded mesh edges and %d embedded mesh faces in region %d", edges.size(), faces.size(), gr->tag());
+ Msg::Info ("Searching for %d embedded mesh edges and %d embedded mesh faces in region %d", edges.size(), faces.size(), gr->tag());
for (unsigned int k=0;k<gr->getNumMeshElements();k++){
for (int j=0;j<gr->getMeshElement(k)->getNumEdges();j++){
edges.erase (gr->getMeshElement(k)->getEdge(j));
@@ -80,7 +78,7 @@ public:
}
}
if (edges.empty() && faces.empty()) {
- Msg::Info ("All embedded edges and faces are present in the final mesh");
+ Msg::Info ("All embedded edges and faces are present in the final mesh");
}
if (edges.size()) {
char name[256];
@@ -904,9 +902,9 @@ static void Mesh3D(GModel *m)
m->setAllVolumesPositive();
// std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
- if (Msg::GetVerbosity() > 98)
+ if (Msg::GetVerbosity() > 98)
std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
-
+
CTX::instance()->mesh.changed = ENT_ALL;
double t2 = Cpu();
CTX::instance()->meshTimer[2] = t2 - t1;
@@ -922,9 +920,9 @@ void OptimizeMeshNetgen(GModel *m)
// Ensure that all volume Jacobians are positive
m->setAllVolumesPositive();
- if (Msg::GetVerbosity() > 98)
+ if (Msg::GetVerbosity() > 98)
std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
-
+
double t2 = Cpu();
Msg::StatusBar(true, "Done optimizing 3D mesh with Netgen (%g s)", t2 - t1);
}
@@ -938,9 +936,9 @@ void OptimizeMesh(GModel *m)
// Ensure that all volume Jacobians are positive
m->setAllVolumesPositive();
- if (Msg::GetVerbosity() > 98)
+ if (Msg::GetVerbosity() > 98)
std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
-
+
CTX::instance()->mesh.changed = ENT_ALL;
double t2 = Cpu();
Msg::StatusBar(true, "Done optimizing 3D mesh (%g s)", t2 - t1);
diff --git a/Mesh/meshDiscreteRegion.cpp b/Mesh/meshDiscreteRegion.cpp
index 3765fec11e6319865e863aff2137b635e05b43fe..facda1ced753e5831cfe137020299f7a46eba81f 100644
--- a/Mesh/meshDiscreteRegion.cpp
+++ b/Mesh/meshDiscreteRegion.cpp
@@ -5,26 +5,18 @@
#include <stdlib.h>
#include <vector>
-// general
-#include "OS.h"
-#include "Geo.h"
-#include "Context.h"
-#include "GmshConfig.h"
#include "GmshMessage.h"
#include "fullMatrix.h"
-// mesh
#include "meshGFaceOptimize.h"
#include "meshGRegionBoundaryRecovery.h"
#include "meshGRegionDelaunayInsertion.h"
#include "meshGRegionRelocateVertex.h"
#include "delaunay3d.h"
-// model
#include "GModel.h"
#include "GRegion.h"
#include "GFace.h"
#include "discreteFace.h"
#include "discreteRegion.h"
-// elements
#include "MLine.h"
#include "MTriangle.h"
#include "MTetrahedron.h"
@@ -35,8 +27,8 @@
// Recursive function to generate all combinations of 4 indices between start
// and end indices included
// Jeanne - HEXTREME
-void combination_of_4( int combination[4], int start, int end, int index,
- const std::vector<MVertex*>& vertices, std::vector<MTetrahedron*>& tets)
+void combination_of_4( int combination[4], int start, int end, int index,
+ const std::vector<MVertex*>& vertices, std::vector<MTetrahedron*>& tets)
{
if (index == 4 ) {
// Create the tet and get out
@@ -57,13 +49,15 @@ void combination_of_4( int combination[4], int start, int end, int index,
// Fill a region with all possible tets genereated from the
// combination of points assigned to it
// Jeanne - HEXTREME
-void create_all_possible_tets(GRegion* region, const std::vector<MVertex*>& vertices) {
+void create_all_possible_tets(GRegion* region, const std::vector<MVertex*>& vertices)
+{
unsigned int nb_points = vertices.size();
int combinaison[4];
std::vector<MTetrahedron*> tets;
combination_of_4(combinaison, 0, nb_points - 1, 0, vertices, tets);
- std::cout << " Number of tets created - all possible combinations - " << tets.size() << std::endl;
+ std::cout << " Number of tets created - all possible combinations - "
+ << tets.size() << std::endl;
for (int i = 0; i < tets.size(); ++i) {
region->addTetrahedron(tets[i]);
@@ -113,7 +107,8 @@ GRegion * createDiscreteRegionFromRawData(GModel *gm, fullMatrix<double> &pts,
gf->triangles.push_back(t);
}
-// delaunayTriangulation(1, 1, vs, gr->tetrahedra); // tetrahedralization is done when recovering the boundary
+ // delaunayTriangulation(1, 1, vs, gr->tetrahedra);
+ // tetrahedralization is done when recovering the boundary
// create all tets
if (all_tets) {
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index b69abd3faea404161a84b4e6316476ceda4a445a..cc15f89f51ad1ad9139671f2a210104f5ed52e73 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -26,8 +26,6 @@
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MElementCut.h"
-#include "CenterlineField.h"
-#include "meshGFaceElliptic.h"
#include "Context.h"
#include "GPoint.h"
#include "GmshMessage.h"
@@ -1868,32 +1866,6 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
return true;
}
-static bool meshGeneratorElliptic(GFace *gf, bool debug = true)
-{
-#if defined(HAVE_ANN)
- Centerline *center = 0;
- FieldManager *fields = GModel::current()->getFields();
- if (fields->getBackgroundField() > 0 ){
- Field *myField = fields->get(fields->getBackgroundField());
- center = dynamic_cast<Centerline*> (myField);
- }
-
- bool recombine = (CTX::instance()->mesh.recombineAll) ;
- int nbBoundaries = gf->edges().size();
-
- if (center && recombine && nbBoundaries == 2) {
- printf("--> regular periodic grid generator (elliptic smooth) \n");
- //bool success = createRegularTwoCircleGrid(center, gf);
- bool success = createRegularTwoCircleGridPeriodic(center, gf);
- return success;
- }
- else return false;
-
-#else
- return false;
-#endif
-}
-
static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
{
std::map<BDS_Point*, MVertex*, PointLessThan> recoverMap;
@@ -2475,11 +2447,6 @@ void meshGFace::operator() (GFace *gf, bool print)
Msg::Debug("Generating the mesh");
- if(meshGeneratorElliptic(gf)){
- gf->meshStatistics.status = GFace::DONE;
- return;
- }
-
quadMeshRemoveHalfOfOneDMesh halfmesh (gf);
if ((gf->getNativeType() != GEntity::AcisModel ||
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
deleted file mode 100644
index 065cf3dfa98b4dc107ca1a230f0aaf5a8cde75d0..0000000000000000000000000000000000000000
--- a/Mesh/meshGFaceElliptic.cpp
+++ /dev/null
@@ -1,880 +0,0 @@
-// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-
-#include <stack>
-#include "GmshConfig.h"
-#include "meshGFaceElliptic.h"
-#include "qualityMeasures.h"
-#include "GFace.h"
-#include "GEdge.h"
-#include "GEdgeCompound.h"
-#include "GVertex.h"
-#include "GModel.h"
-#include "MVertex.h"
-#include "MTriangle.h"
-#include "MQuadrangle.h"
-#include "MLine.h"
-#include "BackgroundMeshTools.h"
-#include "Numeric.h"
-#include "GmshMessage.h"
-#include "Generator.h"
-#include "Context.h"
-#include "OS.h"
-#include "SVector3.h"
-#include "SPoint3.h"
-#include "fullMatrix.h"
-#include "CenterlineField.h"
-#if defined(HAVE_ANN)
-#include "ANN/ANN.h"
-#endif
-
-#define TRAN_QUAD(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
- (1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
-
-static void printQuads(GFace *gf, fullMatrix<SPoint2> uv,
- std::vector<MQuadrangle*> quads, int iter)
-{
- if(!CTX::instance()->mesh.saveAll) return;
-
- char name[234];
- sprintf(name, "quadUV_%d_%d.pos", gf->tag(), iter);
- FILE *f = Fopen(name, "w");
- if(!f){
- Msg::Error("Could not open file '%s'", name);
- return;
- }
-
- fprintf(f,"View \"%s\" {\n",name);
-
- for (int i = 1; i < uv.size1()-1; i++)
- for (int j = 0; j < uv.size2(); j++)
- fprintf(f,"SP(%g,%g,%g) {%d};\n", uv(i,j).x(), uv(i,j).y(), 0.0, i);
-
- fprintf(f,"};\n");
- fclose(f);
-
- char name3[234];
- sprintf(name3, "quadXYZ_%d_%d.pos", gf->tag(), iter);
- FILE *f3 = Fopen(name3,"w");
- if(!f3){
- Msg::Error("Could not open file '%s'", name3);
- return;
- }
-
- fprintf(f3,"View \"%s\" {\n",name3);
- for (unsigned int i = 0; i < quads.size(); i++){
- quads[i]->writePOS(f3,true,false,false,false,false,false);
- }
- fprintf(f3,"};\n");
- fclose(f3);
-}
-
-static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<MVertex*> vert2,
- std::vector<MVertex*> e01, std::vector<MVertex*> e10,
- std::vector<MVertex*> e23, std::vector<MVertex*> e32,
- std::vector<MVertex*> e02, std::vector<MVertex*> e13,
- std::vector<MQuadrangle*> quads)
-{
-
- if(!CTX::instance()->mesh.saveAll) return;
-
- std::vector<SPoint2> p1,p2;
- for (unsigned int i = 0; i< vert1.size(); i++){
- SPoint2 pi;
- reparamMeshVertexOnFace(vert1[i], gf, pi);
- p1.push_back(pi);
- }
- for (unsigned int j = 0; j< vert2.size(); j++){
- SPoint2 pj;
- reparamMeshVertexOnFace(vert2[j], gf, pj);
- p2.push_back(pj);
- }
-
- char name[234];
- sprintf(name,"paramGrid_%d.pos", gf->tag());
- FILE *f = Fopen(name,"w");
- if(!f){
- Msg::Error("Could not open file '%s'", name);
- return;
- }
-
- fprintf(f,"View \"%s\" {\n",name);
-
- // for (unsigned int i = 0; i < p1.size(); i++)
- // fprintf(f,"SP(%g,%g,%g) {%d};\n", p1[i].x(), p1[i].y(), 0.0, i);
- // for (unsigned int j = 0; j < p2.size(); j++)
- // fprintf(f,"SP(%g,%g,%g) {%d};\n", p2[j].x(), p2[j].y(), 0.0, 100+j);
-
- for (unsigned int i = 0; i < p1.size()-1; i++)
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[i].x(), p1[i].y(), 0.0, p1[i+1].x(), p1[i+1].y(), 0.0, 1, 1);
- for (unsigned int i = 0; i < p2.size()-1; i++)
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[i].x(), p2[i].y(), 0.0, p2[i+1].x(), p2[i+1].y(), 0.0, 1, 1);
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[p1.size()-1].x(), p1[ p1.size()-1].y(), 0.0, p1[0].x(), p1[0].y(), 0.0, 1, 1);
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[p2.size()-1].x(), p2[ p2.size()-1].y(), 0.0, p2[0].x(), p2[0].y(), 0.0, 1, 1);
-
- fprintf(f,"};\n");
- fclose(f);
-
- char name2[234];
- sprintf(name2,"paramEdges_%d.pos", gf->tag());
- FILE *f2 = Fopen(name2,"w");
- if(!f2){
- Msg::Error("Could not open file '%s'", name2);
- return;
- }
-
- fprintf(f2,"View \"%s\" {\n",name2);
- for (unsigned int i = 0; i < e01.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e01[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 1);
- }
-
- for (unsigned int i = 0; i < e10.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e10[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 10);
- }
-
- for (unsigned int i = 0; i < e23.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e23[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 23);
- }
-
- for (unsigned int i = 0; i < e32.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e32[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 32);
- }
-
- for (unsigned int i = 0; i < e02.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e02[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 2);
- }
-
- for (unsigned int i = 0; i < e13.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e13[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 13);
- }
- fprintf(f2,"};\n");
- fclose(f2);
-
- char name3[234];
- sprintf(name3,"quadXYZ_%d.pos", gf->tag());
- FILE *f3 = Fopen(name3,"w");
- if(!f3){
- Msg::Error("Could not open file '%s'", name2);
- return;
- }
-
- fprintf(f3,"View \"%s\" {\n",name3);
- for (unsigned int i = 0; i < quads.size(); i++){
- quads[i]->writePOS(f3,true,false,false,false,false,false);
- }
- fprintf(f3,"};\n");
- fclose(f3);
-}
-
-static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv,
- std::vector<std::vector<MVertex*> > &tab,
- std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv)
-{
- newq.clear();
- newv.clear();
-
- int M = uv.size1();
- int N = uv.size2();
-
- for (int i = 1; i < M-1; i++){
- for (int j = 0; j < N-1; j++){
- GPoint gp = gf->point(uv(i,j));
- if (!gp.succeeded()){
- printf("** QUADS FROM UV new vertex not created p=%g %g \n", uv(i,j).x(), uv(i,j).y()); //exit(1);
- }
- tab[i][j] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- }
- }
- for (int i = 1; i < M-1; i++) tab[i][N-1] = tab[i][0];
-
- //create vertices
- for (int i = 1; i < M-1; i++)
- for (int j = 0; j < N-1; j++)
- newv.push_back(tab[i][j]);
-
- // create quads
- for (int i=0;i<M-1;i++){
- for (int j=0;j<N-1;j++){
- MQuadrangle *qnew = new MQuadrangle (tab[i][j],tab[i][j+1],tab[i+1][j+1],tab[i+1][j]);
- newq.push_back(qnew);
- }
- }
-
- return;
-
-}
-static std::vector<MVertex*> saturateEdgeRegular (GFace *gf, SPoint2 p1, SPoint2 p2,
- int M, std::vector<SPoint2> &pe)
-{
- std::vector<MVertex*> pts;
- for (int i=1;i<M;i++){
- double s = ((double)i/((double)(M)));
- SPoint2 p = p1 + (p2-p1)*s;
- pe.push_back(p);
-
- GPoint gp = gf->point(p);
- if (!gp.succeeded()) printf("** SATURATE EDGE new vertex not created p=%g %g \n", p.x(), p.y());
- MVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
-
- if (!v){ pts.clear(); pts.resize(0); return pts;}
- pts.push_back(v);
- }
- return pts;
-}
-
-static std::vector<MVertex*> saturateEdgeHarmonic (GFace *gf, SPoint2 p1, SPoint2 p2,
- double H, double L,
- int M, std::vector<SPoint2> &pe)
-{
- std::vector<MVertex*> pts;
- for (int i=1;i<M;i++){
- double y = ((double)(i))*H/M;
- double Yy = cosh(M_PI*y/L) - tanh(M_PI*H/L)*sinh(M_PI*y/L);
- double YyH = cosh(M_PI*H/L) - tanh(M_PI*H/L)*sinh(M_PI*H/L);
- double s = (1 - Yy)/(1.-YyH);
-
- SPoint2 p = p1 + (p2-p1)*s;
- pe.push_back(p);
-
- GPoint gp = gf->point(p);
- if (!gp.succeeded()) printf("** SATURATE EDGE new vertex not created p=%g %g \n", p.x(), p.y());
- MVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
-
- if (!v){ pts.clear(); pts.resize(0); return pts;}
- pts.push_back(v);
- }
- return pts;
-}
-
-static void transfiniteSmoother(GFace* gf,
- fullMatrix<SPoint2> &uv,
- std::vector<std::vector<MVertex*> > &tab,
- std::vector<MQuadrangle*> &newq,
- std::vector<MVertex*> &newv,
- bool isPeriodic=false)
-{
- int M = uv.size1();
- int N = uv.size2();
- int jStart = isPeriodic ? 0 : 1;
-
- int numSmooth = 150;
- fullMatrix<SPoint2> uvold = uv;
- for(int k = 0; k < numSmooth; k++){
- double norm = 0.0;
- for (int i=1; i<M-1; i++){
- for (int j = jStart; j<N-1; j++){
- int jm1 = (j==0) ? N-2: j-1;
- int jp1 = (isPeriodic ) ? (j+1)%(N-1) : j+1;
- double alpha = 0.25 * (gmsh_SQU(uv(i,jp1).x() - uv(i,jm1).x()) +
- gmsh_SQU(uv(i,jp1).y() - uv(i,jm1).y())) ;
- double gamma = 0.25 * (gmsh_SQU(uv(i+1,j).x() - uv(i-1,j).x()) +
- gmsh_SQU(uv(i+1,j).y() - uv(i-1,j).y()));
- double beta = 0.0625 *
- ((uv(i+1,j).x() - uv(i-1,j).x()) * (uv(i,jp1).x() - uv(i,jm1).x()) +
- (uv(i+1,j).y() - uv(i-1,j).y()) * (uv(i,jp1).y() - uv(i,jm1).y()));
- double uk = 0.5 *
- (alpha * (uv(i+1,j).x() + uv(i-1,j).x()) +
- gamma * (uv(i,jp1).x() + uv(i,jm1).x()) -
- 2. * beta * (uv(i+1,jp1).x() - uv(i-1,jp1).x() -
- uv(i+1,jm1).x() + uv(i-1,jm1).x())) / (alpha + gamma);
- double vk = 0.5 *
- (alpha * (uv(i+1,j).y() + uv(i-1,j).y()) +
- gamma * (uv(i,jp1).y()+ uv(i,jm1).y()) -
- 2. * beta * (uv(i+1,jp1).y() - uv(i-1,jp1).y() -
- uv(i+1,jm1).y() + uv(i-1,jm1).y())) / (alpha + gamma);
- uvold(i,j) = uv(i,j);
- norm += sqrt((uv(i,j).x()-uk)*(uv(i,j).x()-uk)+(uv(i,j).y()-vk)*(uv(i,j).y()-vk));
- uv(i,j) = SPoint2(uk,vk);
- }
- }
- if (isPeriodic){
- for (int i = 1; i<M-1; i++) {
- uv(i,N-1) = uv(i,0);
- uvold(i,N-1) = uvold(i,0);
- }
- }
-
- if(k%20==0){
- printf("Transfinite smoother iter (%d) = %g \n", k, norm);
- createQuadsFromUV(gf, uv, tab, newq, newv);
- printQuads(gf, uv, newq, k+1);
- }
-
- }
-
- //final print
- createQuadsFromUV(gf, uv, tab, newq, newv);
- printQuads(gf, uv, newq, numSmooth);
-}
-
-//elliptic surface grid generator
-//see eq. 9.26 page 9-24 in handbook of grid generation
-static void ellipticSmoother( GFace* gf,
- fullMatrix<SPoint2> &uv,
- std::vector<std::vector<MVertex*> > &tab,
- std::vector<MQuadrangle*> &newq,
- std::vector<MVertex*> &newv,
- bool isPeriodic=false)
-{
- printQuads(gf, uv, newq, 0);
-
- int nbSmooth = 70;
- int M = uv.size1();
- int N = uv.size2();
- int jStart = isPeriodic ? 0 : 1;
- int jEnd = N-1;
-
- fullMatrix<SPoint2> uvold = uv;
- for (int k = 0; k< nbSmooth; k++){
- double norm = 0.0;
- for (int i=1; i<M-1; i++){
- for (int j = jStart; j<jEnd; j++){
- int jm1 = (j==0) ? N-2: j-1;
- int jp1 = (isPeriodic ) ? (j+1)%(N-1) : j+1;
- Pair<SVector3, SVector3> der = gf->firstDer(uv(i,j));
- SVector3 xu = der.first();
- SVector3 xv = der.second();
- // double g11 = dot(xu,xu);
- // double g12 = dot(xu,xv);
- // double g22 = dot(xv,xv);
- SVector3 xuu = SVector3();
- SVector3 xvv = SVector3();
- SVector3 xuv = SVector3();
- gf->secondDer(uv(i,j), &xuu, &xvv, &xuv);
- double g11_bar = dot(xu,xu);
- double g12_bar = dot(xu,xv);
- double g22_bar = dot(xv,xv);
- double dxsi = 1.;
- double deta = 1.;
- double u_xsi = (uv(i,jp1).x()-uv(i,j).x())/dxsi;
- double u_eta = (uv(i+1,j).x()-uv(i,j).x())/deta;
- double v_xsi = (uv(i,jp1).y()-uv(i,j).y())/dxsi;
- double v_eta = (uv(i+1,j).y()-uv(i,j).y())/deta;
- double g11 = g11_bar*u_xsi*u_xsi+2.*g12_bar*u_xsi*v_xsi+g22_bar*v_xsi*v_xsi;
- double g12 = g11_bar*u_xsi*u_eta+g12_bar*(u_xsi*v_eta+u_eta*v_xsi)+g22_bar*v_xsi*v_eta;
- double g22 = g11_bar*u_eta*u_eta+2.*g12_bar*u_eta*v_eta+g22_bar*v_eta*v_eta;
- double jac = u_xsi*v_eta-u_eta*v_xsi; //sqrt(g11*g22-g12*g12);
- double jac_bar = sqrt(g11_bar*g22_bar-g12_bar*g12_bar);
- double g11_bar_u = 2.*dot(xu,xuu);
- double g11_bar_v = 2.*dot(xu,xuv);
- double g22_bar_u = 2.*dot(xv,xuv);
- double g22_bar_v = 2.*dot(xv,xvv);
- double g12_bar_u = dot(xu,xuv)+dot(xv,xuu);
- double g12_bar_v = dot(xu,xvv)+dot(xv,xuv);
- double jac_bar_u = 1./(2.*jac_bar)*(g11_bar*g22_bar_u+g22_bar*g11_bar_u-2.*g12_bar*g12_bar_u);
- double jac_bar_v = 1./(2.*jac_bar)*(g11_bar*g22_bar_v+g22_bar*g11_bar_v-2.*g12_bar*g12_bar_v);
- double lapl_u = 1./jac_bar*(jac_bar*(g22_bar_u-g12_bar_v)-(g22_bar*jac_bar_u-g12_bar*jac_bar_v));
- double lapl_v = 1./jac_bar*(jac_bar*(g11_bar_v-g12_bar_u)-(g11_bar*jac_bar_v-g12_bar*jac_bar_u));
- double over = 1./(4.*(g11+g22));
- double uk = over*(2.*g22*(uv(i+1,j).x()+uv(i-1,j).x())+
- 2.*g11*(uv(i,jp1).x()+uv(i,jm1).x())-
- 2*g12*(uv(i+1,jp1).x()-uv(i-1,jp1).x()-uv(i+1,jm1).x()+uv(i-1,jm1).x())-
- 2.*jac*jac*lapl_u);
- double vk = over*(2.*g22*(uv(i+1,j).y()+uv(i-1,j).y())+
- 2.*g11*(uv(i,jp1).y()+uv(i,jm1).y())-
- 2.*g12*(uv(i+1,jp1).y()-uv(i-1,jp1).y()-uv(i+1,jm1).y()+uv(i-1,jm1).y())-
- 2.*jac*jac*lapl_v);
- uvold(i,j) = uv(i,j);
- norm += sqrt((uv(i,j).x()-uk)*(uv(i,j).x()-uk)+(uv(i,j).y()-vk)*(uv(i,j).y()-vk));
- uv(i,j) = SPoint2(uk,vk);
- }
- }
- if (isPeriodic){
- for (int i = 1; i<M-1; i++) {
- uv(i,N-1) = uv(i,0);
- uvold(i,N-1) = uvold(i,0);
- }
- }
-
- //under-relaxation
- // double omega = 0.8;
- // for (int i=0; i<M; i++){
- // for (int j = 0; j<N; j++){
- // uv(i,j) = SPoint2(omega*uv(i,j).x() + (1.-omega)*uvold(i,j).x(),
- // omega*uv(i,j).y() + (1.-omega)*uvold(i,j).y());
- // }
- // }
-
-
- if(k%10==0){
- printf("Elliptic smoother iter (%d) = %g \n", k, norm);
- createQuadsFromUV(gf, uv, tab, newq, newv);
- printQuads(gf, uv, newq, k+1);
- }
-
- }
-
- createQuadsFromUV(gf, uv, tab, newq, newv);
- printQuads(gf, uv, newq, nbSmooth);
-}
-
-//create initial grid points MxN using transfinite interpolation
-/*
- c4 N c3
- +--------------------------+
- | | | | |
- | | | | |
- +--------------------------+ M
- | | | | |
- | | | | |
- +--------------------------+
- c1 c2
- N
-*/
-static void createRegularGrid (GFace *gf,
- MVertex *v1, SPoint2 &c1,
- std::vector<MVertex*> &e12, std::vector<SPoint2> &pe12, int sign12,
- MVertex *v2, SPoint2 &c2,
- std::vector<MVertex*> &e23, std::vector<SPoint2> &pe23, int sign23,
- MVertex *v3, SPoint2 &c3,
- std::vector<MVertex*> &e34,std::vector<SPoint2> &pe34, int sign34,
- MVertex *v4, SPoint2 &c4,
- std::vector<MVertex*> &e41, std::vector<SPoint2> &pe41,int sign41,
- std::vector<MQuadrangle*> &q,
- std::vector<MVertex*> &newv,
- fullMatrix<SPoint2> &uv,
- std::vector<std::vector<MVertex*> > &tab)
-{
- int M = e23.size();
- int N = e12.size();
-
- uv.resize(M+2,N+2);
-
- //std::vector<std::vector<MVertex*> > tab(M+2);
- tab.resize(M+2);
- for(int i = 0; i <= M+1; i++) tab[i].resize(N+2);
-
- tab[0][0] = v1; uv(0,0) = c1;
- tab[0][N+1] = v2; uv(0,N+1) = c2;
- tab[M+1][N+1] = v3; uv(M+1,N+1) = c3;
- tab[M+1][0] = v4; uv(M+1,0) = c4;
-
- for (int i=0;i<N;i++){
- tab[0][i+1] = sign12 > 0 ? e12 [i] : e12 [N-i-1];
- tab[M+1][i+1] = sign34 < 0 ? e34 [i] : e34 [N-i-1];
- uv(0,i+1) = sign12 > 0 ? pe12 [i] : pe12 [N-i-1];
- uv(M+1,i+1) = sign34 < 0 ? pe34 [i] : pe34 [N-i-1];
- }
- for (int i=0;i<M;i++){
- tab[i+1][N+1] = sign23 > 0 ? e23 [i] : e23 [M-i-1];
- tab[i+1][0] = sign41 < 0 ? e41 [i] : e41 [M-i-1];
- uv(i+1,N+1) = sign23 > 0 ? pe23 [i] : pe23 [M-i-1];
- uv(i+1,0) = sign41 < 0 ? pe41 [i] : pe41 [M-i-1];
- }
-
- for (int i=0;i<N;i++){
- for (int j=0;j<M;j++){
- double u = (double) (i+1)/ ((double)(N+1));
- double v = (double) (j+1)/ ((double)(M+1));
- SPoint2 p12 = (sign12 >0) ? pe12[i] : pe12 [N-1-i];
- SPoint2 p23 = (sign23 >0) ? pe23[j] : pe23 [M-1-j];
- SPoint2 p34 = (sign34 <0) ? pe34[i] : pe34 [N-1-i];
- SPoint2 p41 = (sign41 <0) ? pe41[j] : pe41 [M-1-j];
- double Up = TRAN_QUAD(p12.x(), p23.x(), p34.x(), p41.x(),
- c1.x(),c2.x(),c3.x(),c4.x(),u,v);
- double Vp = TRAN_QUAD(p12.y(), p23.y(), p34.y(), p41.y(),
- c1.y(),c2.y(),c3.y(),c4.y(),u,v);
- if ((p12.x() && p12.y() == -1.0) || (p23.x() && p23.y() == -1.0) ||
- (p34.x() && p34.y() == -1.0) || (p41.x() && p41.y() == -1.0)) {
- Msg::Error("Wrong param -1");
- return;
- }
-
- SPoint2 pij(Up,Vp);
-
- GPoint gp = gf->point(pij);
- if (!gp.succeeded()) printf("** INITIAL GRID new vertex not created p=%g %g \n", Up, Vp);
- MVertex *vnew = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- newv.push_back(vnew);
-
- uv(j+1,i+1) = pij;
- tab[j+1][i+1] = vnew;
- }
- }
-
- // create quads
- for (int i=0;i<=N;i++){
- for (int j=0;j<=M;j++){
- MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
- q.push_back(qnew);
- }
- }
-
-}
-
-//create initial grid points MxN using transfinite interpolation
-/*
- c2 N c2
- +--------------------------+
- | | | | |
- | | | | |
- +--------------------------+ M
- | | | | |
- | | | | |
- +--------------------------+
- c0 c0
- N
-*/
-static void createRegularGridPeriodic (GFace *gf,int sign2,
- MVertex *v0, SPoint2 &c0,
- MVertex *v2, SPoint2 &c2,
- std::vector<MVertex*> &e02, std::vector<SPoint2> &pe02,
- std::vector<MVertex*> &e00, std::vector<SPoint2> &pe00,
- std::vector<MVertex*> &e22, std::vector<SPoint2> &pe22,
- std::vector<MQuadrangle*> &q,
- std::vector<MVertex*> &newv,
- fullMatrix<SPoint2> &uv,
- std::vector<std::vector<MVertex*> > &tab)
-{
-
- int M = e02.size();
- int N = e00.size();
-
- uv.resize(M+2,N+2);
-
- char name3[234];
- sprintf(name3,"quadParam_%d.pos", gf->tag());
- FILE *f3 = Fopen(name3,"w");
-
- if(f3) fprintf(f3,"View \"%s\" {\n",name3);
-
- tab.resize(M+2);
- for(int i = 0; i < M+2; i++) tab[i].resize(N+2);
-
- //periodic boundary mesh vertices
- tab[0][0] = v0; uv(0,0) = c0;
- tab[0][N+1] = v0; uv(0,N+1) = c0;
- tab[M+1][N+1] = v2; uv(M+1,N+1) = c2;
- tab[M+1][0] = v2; uv(M+1,0) = c2;
- for (int i=0;i<N;i++){
- tab[0][i+1] = e00 [i]; uv(0,i+1) = pe00 [i] ;
- tab[M+1][i+1] = (sign2 > 0) ? e22 [i] : e22 [N-i-1] ;
- uv(M+1,i+1) = (sign2 > 0) ? pe22 [i] : pe22 [N-i-1] ;
- }
- for (int i=0;i<M;i++){
- tab[i+1][N+1] = e02 [i]; uv(i+1,N+1) = pe02 [i] ;
- tab[i+1][0] = e02 [i]; uv(i+1,0) = pe02 [i];
- }
-
- //interior mesh_vertices
- for (int i=0;i<N;i++){
- SPoint2 p0 = pe00[i] ;
- SPoint2 p2 = (sign2>0) ? pe22[i] : pe22 [N-i-1] ;
- std::vector<SPoint2> pei;
- std::vector<MVertex*> ei = saturateEdgeRegular(gf,p0,p2,M+1,pei);//M+1
- for (int j=0;j<M;j++){
- SPoint2 pij = pei[j];
- GPoint gp = gf->point(pij);
- if (!gp.succeeded()) printf("** INITIAL GRID new vertex not created p=%g %g \n", pij.x(), pij.y());
- tab[j+1][i+1] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- uv(j+1,i+1) = pij;
- }
- }
-
- //create vertices
- for (int i=0;i<N+1;i++)
- for (int j=1;j<M+1;j++)
- newv.push_back(tab[j][i]);
-
- // create quads
- for (int i=0;i<=N;i++){
- for (int j=0;j<=M;j++){
- MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
- q.push_back(qnew);
- }
- }
-
- //print
- if(f3){
- for (int i=0;i<N+2;i++)
- for (int j=0;j<M+2;j++)
- fprintf(f3,"SP(%g,%g,%g) {%d};\n", uv(j,i).x(), uv(j,i).y(), 0.0, j);
- fprintf(f3,"};\n");
- fclose(f3);
- }
-
-}
-
-static void updateFaceQuads(GFace *gf, std::vector<MQuadrangle*> &quads, std::vector<MVertex*> &newv)
-{
- for (unsigned int i = 0; i < quads.size(); i++){
- gf->quadrangles.push_back(quads[i]);
- }
- for(unsigned int i = 0; i < newv.size(); i++){
- gf->mesh_vertices.push_back(newv[i]);
- }
-}
-
-static bool computeRingVertices(GFace *gf, Centerline *center,
- std::vector<MVertex*> &vert1, std::vector<MVertex*> &vert2,
- int &N, int &M, int &close_ind, int &sign2,
- double &arc, double &length)
-{
-#if defined(HAVE_ANN)
- std::list<GEdge*> bedges = gf->edges();
- std::list<GEdge*>::iterator itb = bedges.begin();
- std::list<GEdge*>::iterator ite = bedges.end(); ite--;
- GEdge *ge1 = (*itb)->getCompound();
- GEdge *ge2 = (*ite)->getCompound();
- N = ge1->mesh_vertices.size() + 1;
- int N2 = ge2->mesh_vertices.size() + 1;
- if (N != N2 || N%2 != 0){
- Msg::Error("You should an even number nbPoints in centerline =%d \n", N);
- exit(1);
- return false;
- }
-
- vert1.push_back(ge1->getBeginVertex()->mesh_vertices[0]);
- vert2.push_back(ge2->getBeginVertex()->mesh_vertices[0]);
- for(int i = 0; i < N-1; i++) {
- vert1.push_back(ge1->mesh_vertices[i]);
- vert2.push_back(ge2->mesh_vertices[i]);
- }
- SPoint2 pv1; reparamMeshVertexOnFace(vert1[0], gf, pv1);
- SPoint2 pv2; reparamMeshVertexOnFace(vert2[0], gf, pv2);
- SPoint2 pv1b; reparamMeshVertexOnFace(vert1[1], gf, pv1b);
- SPoint2 pv2b; reparamMeshVertexOnFace(vert2[1], gf, pv2b);
- SPoint2 pv1c; reparamMeshVertexOnFace(vert1[2], gf, pv1c);
- SPoint2 pv2c; reparamMeshVertexOnFace(vert2[2], gf, pv2c);
- SVector3 vec1(pv1b.x()-pv1.x(),pv1b.y()-pv1.y() , 0.0);
- SVector3 vec1b(pv1c.x()-pv1b.x(),pv1c.y()-pv1b.y() ,0.0);
- SVector3 vec2(pv2b.x()-pv2.x(),pv2b.y()-pv2.y() ,0.0);
- SVector3 vec2b(pv2c.x()-pv2b.x(),pv2c.y()-pv2b.y() , 0.0);
- sign2 = (dot(crossprod(vec1,vec1b),crossprod(vec2,vec2b)) < 0) ? -1 : +1 ;
- double n1 = sqrt(pv1.x()*pv1.x()+pv1.y()*pv1.y());
- double n2 = sqrt(pv2.x()*pv2.x()+pv2.y()*pv2.y());
- std::vector<MVertex*> vert_temp;
- if (n2 > n1) {
- vert_temp = vert1;
- vert1.clear();
- vert1 = vert2;
- vert2.clear();
- vert2 = vert_temp;
- }
- double rad = center->operator()(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
-
- ANNpointArray nodes = annAllocPts(N, 3);
- ANNidxArray index = new ANNidx[1];
- ANNdistArray dist = new ANNdist[1];
- for (int ind = 0; ind < N; ind++){
- SPoint2 pp2; reparamMeshVertexOnFace(vert2[ind], gf, pp2);
- nodes[ind][0] = pp2.x();
- nodes[ind][1] = pp2.y();
- nodes[ind][2] = 0.0;
- }
- ANNkd_tree *kdtree = new ANNkd_tree(nodes, N, 3);
- SPoint2 pp1; reparamMeshVertexOnFace(vert1[0], gf, pp1);
- double xyz[3] = {pp1.x(), pp1.y(),0.0};
- kdtree->annkSearch(xyz, 1, index, dist);
- close_ind = index[0];
-
- SPoint2 p1; reparamMeshVertexOnFace(vert1[0], gf, p1);
- SPoint2 p2; reparamMeshVertexOnFace(vert2[close_ind], gf, p2);
- GPoint gp_i(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
- SPoint2 p1b; reparamMeshVertexOnFace(vert1[N/2], gf, p1b);
- SPoint2 p2b; reparamMeshVertexOnFace(vert2[(close_ind+N/2)%N], gf, p2b);
- GPoint gp_ib(vert1[N/2]->x(), vert1[N/2]->y(), vert1[N/2]->z());
- length = 0.0;
- double lengthb = 0.0;
- int nb = 100;
- for (int i = 0; i< nb; i++){
- SPoint2 pii(p1.x() + (double)(i+1)/nb*(p2.x()-p1.x()), p1.y() + (double)(i+1)/nb*(p2.y()-p1.y()));
- GPoint gp_ii = gf->point(pii);
- SPoint2 piib(p1b.x() + (double)(i+1)/nb*(p2b.x()-p1b.x()), p1b.y() + (double)(i+1)/nb*(p2b.y()-p1b.y()));
- GPoint gp_iib = gf->point(piib);
- length += gp_i.distance(gp_ii);
- lengthb += gp_ib.distance(gp_iib);
- gp_i = gp_ii;
- gp_ib = gp_iib;
- }
-
- arc = 2*M_PI*rad;
- double lc = arc/N;
- M = (int)(0.5*(length+lengthb)/lc) ;
-
- delete kdtree;
- delete[]index;
- delete[]dist;
- annDeallocPts(nodes);
-
- return true;
-#else
- return false;
-#endif
-}
-
-//vert1 is the outer circle
-//vert2 is the inner circle
-// - - - -
-// - -
-// - - -
-// v0- v2- - -
-// - - -
-// - -
-// - - - -
-bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf)
-{
-#if defined(HAVE_ANN)
- std::vector<MVertex*> vert1, vert2;
- int N, M, close_ind, sign2;
- double arc, length;
- bool success = computeRingVertices(gf, center, vert1, vert2, N, M, close_ind, sign2, arc,length);
- if(!success) return false;
-
- MVertex *v0 = vert1[0]; SPoint2 p0; reparamMeshVertexOnFace(v0, gf, p0);
- MVertex *v2 = vert2[close_ind]; SPoint2 p2; reparamMeshVertexOnFace(v2, gf, p2);
-
- //printf("grid N = %d M = %d\n", N , M);
- std::vector<MVertex*> e00,e22;//edges without first and last vertex
- for (int i=1;i<N;i++) e00.push_back(vert1[i]);
- for (int i=1;i<N;i++) e22.push_back(vert2[(close_ind+i)%N]);
-
- std::vector<SPoint2> pe00, pe22;
- for (unsigned int i = 0; i < e00.size(); i++){
- SPoint2 p00; reparamMeshVertexOnFace(e00[i], gf, p00);
- SPoint2 p22; reparamMeshVertexOnFace(e22[i], gf, p22);
- pe00.push_back(p00);
- pe22.push_back(p22);
- }
-
- std::vector<SPoint2> pe02;
- std::vector<MVertex*> e02 = saturateEdgeRegular(gf,p0,p2,M+1,pe02);
-
- std::vector<MQuadrangle*> quads;
- std::vector<MVertex*> newv;
- fullMatrix<SPoint2> uv;
- std::vector<std::vector<MVertex*> > tab;
-
- createRegularGridPeriodic (gf,sign2,
- v0, p0,
- v2, p2,
- e02, pe02,
- e00, pe00,
- e22, pe22,
- quads, newv, uv, tab);
-
- printQuads(gf, uv, quads, 0);
-
- transfiniteSmoother(gf, uv, tab, quads, newv, true);
- //ellipticSmoother(gf, uv, tab, quads, newv, true);
-
- updateFaceQuads(gf, quads, newv);
-
- //exit(1);
- //printParamGrid(gf, vert1, vert2, e00,e22,e02,e02,e02,e02, quads);
-
- return true;
-#else
- return false;
-#endif
-}
-
-//vert1 is the outer circle
-//vert2 is the inner circle
-// - - - -
-// - -
-// - - -
-// v0- v2- -v3 -v1
-// - - -
-// - -
-// - - - -
-
-bool createRegularTwoCircleGrid (Centerline *center, GFace *gf)
-{
-#if defined(HAVE_ANN)
- std::vector<MVertex*> vert1, vert2;
- int N, M, close_ind, sign2;
- double arc, length;
- bool success = computeRingVertices(gf, center, vert1, vert2, N, M, close_ind, sign2, arc,length);
- if(!success) return false;
-
- MVertex *v0 = vert1[0]; SPoint2 p0; reparamMeshVertexOnFace(v0, gf, p0);
- MVertex *v1 = vert1[N/2]; SPoint2 p1; reparamMeshVertexOnFace(v1, gf, p1);
- MVertex *v2 = vert2[close_ind]; SPoint2 p2; reparamMeshVertexOnFace(v2, gf, p2);
- MVertex *v3 = vert2[(close_ind+N/2)%N]; SPoint2 p3; reparamMeshVertexOnFace(v3, gf, p3);
-
- printf("grid N = %d M = %d\n", N , M);
- std::vector<MVertex*> e01,e10,e23,e32;//edges without first and last vertex
- for (int i=1;i<N/2;i++) e01.push_back(vert1[i]);
- for (int i=N/2+1;i<N;i++) e10.push_back(vert1[i]);
- for (int i=1;i<N/2;i++) e23.push_back(vert2[(close_ind+i)%N]);
- for (int i=N/2+1;i<N;i++) e32.push_back(vert2[(close_ind+i)%N]);
-
- std::vector<SPoint2> pe01, pe10, pe23, pe32;
- for (unsigned int i = 0; i < e01.size(); i++){
- SPoint2 p01; reparamMeshVertexOnFace(e01[i], gf, p01);
- SPoint2 p10; reparamMeshVertexOnFace(e10[i], gf, p10);
- SPoint2 p23; reparamMeshVertexOnFace(e23[i], gf, p23);
- SPoint2 p32; reparamMeshVertexOnFace(e32[i], gf, p32);
- pe01.push_back(p01);
- pe10.push_back(p10);
- pe23.push_back(p23);
- pe32.push_back(p32);
- }
-
- std::vector<SPoint2> pe02, pe13;
- std::vector<MVertex*> e02 = saturateEdgeHarmonic (gf,p0,p2,length, arc, M+1, pe02);
- std::vector<MVertex*> e13 = saturateEdgeHarmonic (gf,p1,p3,length, arc, M+1, pe13);
-
- std::vector<MVertex*> e_inner1 = e23;
- std::vector<MVertex*> e_inner2 = e32;
- std::vector<SPoint2> pe_inner1 = pe23;
- std::vector<SPoint2> pe_inner2 = pe32;
- if (sign2 <0){
- e_inner1 = e32;
- e_inner2 = e23;
- pe_inner1 = pe32;
- pe_inner2 = pe23;
- }
-
- std::vector<MQuadrangle*> quads, quadS1, quadS2;
- std::vector<MVertex*> newv, newv1, newv2;
- fullMatrix<SPoint2> uv1, uv2;
- std::vector<std::vector<MVertex*> > tab1, tab2;
- createRegularGrid (gf,
- v0, p0,
- e01, pe01, +1,
- v1,p1,
- e13, pe13, +1,
- v3,p3,
- e_inner1, pe_inner1, -sign2,
- v2,p2,
- e02, pe02, -1,
- quads, newv, uv1, tab1);
-
- createRegularGrid (gf,
- v0,p0,
- e02, pe02, +1,
- v2,p2,
- e_inner2, pe_inner2, -sign2,
- v3,p3,
- e13, pe13, -1,
- v1,p1,
- e10, pe10, +1,
- quads, newv, uv2, tab2);
-
- printQuads(gf, uv1, quads, 0);
-
- ellipticSmoother(gf, uv1, tab1, quadS1, newv1);
- ellipticSmoother(gf, uv2, tab2, quadS2, newv2);
-
- quads.clear();
- for (unsigned int i= 0; i< quadS1.size(); i++) quads.push_back(quadS1[i]);
- for (unsigned int i= 0; i< quadS2.size(); i++) quads.push_back(quadS2[i]);
- //updateFaceQuads(gf, quads, newv);
-
- printParamGrid(gf, vert1, vert2, e01,e10,e23,e32,e02,e13, quads);
-
- return true;
-#else
- return false;
-#endif
-}
diff --git a/Mesh/meshGFaceElliptic.h b/Mesh/meshGFaceElliptic.h
deleted file mode 100644
index 42fd9e9fe7116b33703ca603562788b4cefe9b34..0000000000000000000000000000000000000000
--- a/Mesh/meshGFaceElliptic.h
+++ /dev/null
@@ -1,25 +0,0 @@
-// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to the public mailing list <gmsh@onelab.info>.
-
-#ifndef _MESH_GFACE_ELLIPTIC_H_
-#define _MESH_GFACE_ELLIPTIC_H_
-
-#include <map>
-#include <vector>
-#include "MElement.h"
-#include "MEdge.h"
-#include "STensor3.h"
-
-class GFace;
-class GVertex;
-class MVertex;
-class Centerline;
-
-
-bool createRegularTwoCircleGrid (Centerline *center, GFace *gf);
-bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf);
-
-
-#endif
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index dd988806ae780dac50114851a81001bc01ba034a..c53d4884245e5cf42262b7e915bfeb2b530661b0 100644
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -6,22 +6,21 @@
// Contributor(s):
// Tristan Carrier François Henrotte
-
+#include <queue>
+#include <fstream>
+#include <algorithm>
+#include <iostream>
+#include <string>
#include "simple3D.h"
#include "GModel.h"
#include "MElement.h"
#include "MElementOctree.h"
#include "meshGRegion.h"
-#include <queue>
-#include <fstream>
-#include "CenterlineField.h"
-#include <algorithm>
#include "directions3D.h"
#include "ThinLayer.h"
#include "Context.h"
-#include <iostream>
-#include <string>
#include "rtree.h"
+#include "Field.h"
#define k1 0.7 //k1*h is the minimal distance between two nodes
#define k2 0.5 //k2*h is the minimal distance to the boundary