Select Git revision
-
Christophe Geuzaine authoredChristophe Geuzaine authored
CenterlineField.cpp 42.30 KiB
// Gmsh - Copyright (C) 1997-2014 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@geuz.org>.
//
// 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 "BackgroundMesh.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");
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);
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");
// 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");
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");
// 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");
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