diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 3f3b477c5d843c8f2ff409767005d26588bbe68b..d0037053ed7d8b5e3296fbab8663c30ab4688840 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -3,6 +3,7 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+#include <stdlib.h>
#include "GModelFactory.h"
#include "ListUtils.h"
#include "Context.h"
@@ -62,7 +63,7 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > edges)
{
-
+
//create line loops
int nLoops = edges.size();
std::vector<EdgeLoop *> vecLoops;
@@ -77,8 +78,8 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
if(!c){
GVertex *gvb = ge->getBeginVertex();
GVertex *gve = ge->getEndVertex();
- Vertex *vertb = FindPoint((int)fabs(gvb->tag()));
- Vertex *verte = FindPoint((int)fabs(gve->tag()));
+ Vertex *vertb = FindPoint(abs(gvb->tag()));
+ Vertex *verte = FindPoint(abs(gve->tag()));
if (!vertb){
vertb = Create_Vertex(gvb->tag(), gvb->x(), gvb->y(), gvb->z(),
gvb->prescribedMeshSizeAtVertex(), 1.0);
@@ -121,7 +122,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
List_Delete(temp);
}
List_Add(temp, &numEdge);
- }
+ }
int num = gm->getMaxElementaryNumber(2) + 1+i;
while (FindSurfaceLoop(num)){
@@ -132,9 +133,9 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
EdgeLoop *l = Create_EdgeLoop(num, temp);
vecLoops.push_back(l);
Tree_Add(gm->getGEOInternals()->EdgeLoops, &l);
- List_Delete(temp);
- }
-
+ List_Delete(temp);
+ }
+
//create surface
int numf = gm->getMaxElementaryNumber(2)+1;
Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index aea0e12caec9407e5987c23d161a53a8927c1524..67f524934745798842800d8f25cb631b5c60fbb2 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -48,7 +48,7 @@ double computeLength(std::vector<MLine*> lines){
double length= 0.0;
for (int i = 0; i< lines.size(); i++){
- length += lines[i]->getLength();
+ length += lines[i]->getLength();
}
return length;
}
@@ -77,7 +77,7 @@ bool isClosed(std::set<MEdge, Less_Edge> &theCut){
// }
// // find a lonely MEdge
- // for (std::list<MEdge>::iterator it = segments.begin();
+ // for (std::list<MEdge>::iterator it = segments.begin();
// it != segments.end(); ++it){
// MVertex *vL = it->getVertex(0);
// MVertex *vR = it->getVertex(1);
@@ -89,7 +89,7 @@ bool isClosed(std::set<MEdge, Less_Edge> &theCut){
// else boundv.erase(it2);
// }
- // if (boundv.size() == 0) return true;
+ // if (boundv.size() == 0) return true;
// else return false;
}
@@ -107,7 +107,7 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
}
// find a lonely MLine
- for (std::list<MLine*>::iterator it = segments.begin();
+ for (std::list<MLine*>::iterator it = segments.begin();
it != segments.end(); ++it){
MVertex *vL = (*it)->getVertex(0);
MVertex *vR = (*it)->getVertex(1);
@@ -129,7 +129,7 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
if (v == vB) firstLine = (boundv.rbegin())->second;
else{ Msg::Error("begin vertex not found for branch"); exit(1);}
}
- for (std::list<MLine*>::iterator it = segments.begin();
+ for (std::list<MLine*>::iterator it = segments.begin();
it != segments.end(); ++it){
if (*it == firstLine){
segments.erase(it);
@@ -149,7 +149,7 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
while (first != last){
if (segments.empty())break;
bool found = false;
- for (std::list<MLine*>::iterator it = segments.begin();
+ for (std::list<MLine*>::iterator it = segments.begin();
it != segments.end(); ++it){
MLine *e = *it;
if (e->getVertex(0) == last){
@@ -185,12 +185,12 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
//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++)
+ for (unsigned int i = 0; i < lines.size(); i++)
_orientation[i] = !_orientation[i];
}
@@ -198,7 +198,7 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
// else if (junctions.find(lines[0]->getVertex(1)) != junctions.end()) vB = lines[0]->getVertex(1);
// else printf("no vB junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
// if (junctions.find(lines[lines.size()-1]->getVertex(0)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(0);
- // else if (junctions.find(lines[lines.size()-1]->getVertex(1)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(1);
+ // else if (junctions.find(lines[lines.size()-1]->getVertex(1)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(1);
// else printf("no vE junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
// printf("in order vB= %d =%d \n", vB->getNum(), vE->getNum());
}
@@ -206,7 +206,7 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
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> &touched,
std::set<MEdge, Less_Edge> &theCut)
{
if (touched.find(e) != touched.end()) return;
@@ -216,8 +216,8 @@ static void recurConnectByMEdge(const MEdge &e,
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;
+ if (theCut.find(me) != theCut.end()){
+ touched.insert(me); //break;
}
else recurConnectByMEdge(me, e2e, group, touched, theCut);
}
@@ -225,19 +225,19 @@ static void recurConnectByMEdge(const MEdge &e,
}
-void cutTriangle(MTriangle *tri,
- std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
+void cutTriangle(MTriangle *tri,
+ std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
std::set<MVertex*> &cutVertices,
- std::vector<MTriangle*> &newTris,
+ 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);
+ 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);
@@ -316,7 +316,7 @@ void cutTriangle(MTriangle *tri,
}
Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
-
+
recombine = CTX::instance()->mesh.recombineAll;
index = new ANNidx[1];
@@ -378,14 +378,14 @@ void Centerline::importFile(std::string fileName){
}
if(triangles.empty()){
- Msg::Error("Current GModel has no triangles ...");
+ Msg::Error("Current GModel has no triangles ...");
exit(1);
}
-
+
mod = new GModel();
mod->load(fileName);
mod->removeDuplicateMeshVertices(1.e-8);
- current->setAsCurrent();
+ current->setAsCurrent();
int maxN = 0.0;
std::vector<GEdge*> modEdges = mod->bindingsGetEdges();
@@ -408,7 +408,7 @@ void Centerline::importFile(std::string fileName){
}
createBranches(maxN);
-
+
}
@@ -435,7 +435,7 @@ void Centerline::createBranches(int maxN){
junctions.insert(*it);
}
}
-
+
//split edges
int tag = 0;
for(unsigned int i = 0; i < color_edges.size(); ++i){
@@ -446,7 +446,7 @@ void Centerline::createBranches(int maxN){
MVertex *vB = (*itl)->getVertex(0);
MVertex *vE = (*itl)->getVertex(1);
myLines.push_back(*itl);
- erase(lines, *itl);
+ erase(lines, *itl);
itl = lines.begin();
while ( !( junctions.find(vE) != junctions.end() &&
junctions.find(vB) != junctions.end()) ) {
@@ -481,14 +481,14 @@ void Centerline::createBranches(int maxN){
else itl++;
}
if (vB == vE) { Msg::Error("Begin and end points branch are the same \n");break;}
- orderMLines(myLines, vB, vE);
+ 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) ;
}
}
printf("*** Centerline in/outlets =%d branches =%d \n", (int)color_edges.size()+1, (int)edges.size());
-
+
//create children
for(unsigned int i = 0; i < edges.size(); ++i) {
MVertex *vE = edges[i].vE;
@@ -513,7 +513,7 @@ void Centerline::distanceToSurface(){
Msg::Info("Centerline: computing distance to surface mesh ");
//COMPUTE WITH REVERSE ANN TREE (SURFACE POINTS IN TREE)
- std::set<MVertex*> allVS;
+ std::set<MVertex*> allVS;
for(int j = 0; j < triangles.size(); j++)
for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
int nbSNodes = allVS.size();
@@ -522,13 +522,13 @@ void Centerline::distanceToSurface(){
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();
+ nodesR[ind][0] = v->x();
+ nodesR[ind][1] = v->y();
+ nodesR[ind][2] = v->z();
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);
@@ -536,9 +536,9 @@ void Centerline::distanceToSurface(){
double midp[3] = {0.5*(v1->x()+v2->x()), 0.5*(v1->y()+v1->y()),0.5*(v1->z()+v2->z())};
kdtreeR->annkSearch(midp, 1, index, dist);
double minRad = sqrt(dist[0]);
- radiusl.insert(std::make_pair(lines[i], minRad));
+ radiusl.insert(std::make_pair(lines[i], minRad));
}
-
+
}
void Centerline::computeRadii(){
@@ -552,7 +552,7 @@ void Centerline::computeRadii(){
edges[i].maxRad = std::max(itr->second, edges[i].maxRad);
}
else printf("ARGG line not found \n");
- }
+ }
}
}
@@ -570,12 +570,12 @@ void Centerline::buildKdTree(){
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();
+ 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){
+ 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());
@@ -583,9 +583,9 @@ void Centerline::buildKdTree(){
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();
+ nodes[ind][0] = Pj.x();
+ nodes[ind][1] = Pj.y();
+ nodes[ind][2] = Pj.z();
ind++;
}
}
@@ -619,23 +619,23 @@ void Centerline::createSplitCompounds(){
GEdge *pe = current->getEdgeByTag(i+1);//current edge
e_compound.push_back(pe);
int num_gec = NE+i+1;
- Msg::Info("Parametrize Compound Line (%d) = %d discrete edge",
+ Msg::Info("Parametrize Compound Line (%d) = %d discrete edge",
num_gec, pe->tag());
GEdgeCompound *gec = new GEdgeCompound(current, num_gec, e_compound);
current->add(gec);
}
-
+
// Parametrize Compound surfaces
std::list<GEdge*> U0;
for (int i=0; i < NF; i++){
std::list<GFace*> f_compound;
- GFace *pf = current->getFaceByTag(i+1);//current face
- f_compound.push_back(pf);
- int num_gfc = NF+i+1;
+ GFace *pf = current->getFaceByTag(i+1);//current face
+ f_compound.push_back(pf);
+ int num_gfc = NF+i+1;
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONICPLANE;
- //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
+ GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONICPLANE;
+ //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
@@ -663,8 +663,8 @@ void Centerline::cleanMesh(){
Msg::Info("Writing new splitted mesh mySPLITMESH.msh");
current->writeMSH("mySPLITMESH.msh", 2.2, false, false);
- std::set<MVertex*> allNod;
- discreteFace * mySplitMesh;
+ std::set<MVertex*> allNod;
+ discreteFace * mySplitMesh;
std::vector<std::set<MVertex*> > inOutNod;
std::vector<discreteFace* > inOutMesh;
@@ -693,7 +693,7 @@ void Centerline::cleanMesh(){
mySplitMesh->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3]));
}
}
-
+
//Removing discrete Vertices - Edges - Faces
for (int i=0; i < NV; i++){
GVertex *gv = current->getVertexByTag(i+1);
@@ -702,7 +702,7 @@ void Centerline::cleanMesh(){
for (int i=0; i < NE; i++){
GEdge *ge = current->getEdgeByTag(i+1);
GEdge *gec = current->getEdgeByTag(NE+i+1);
- current->remove(ge);
+ current->remove(ge);
current->remove(gec);
}
for (int i=0; i < NF; i++){
@@ -711,19 +711,19 @@ void Centerline::cleanMesh(){
}
for (int i=0; i < discFaces.size(); i++){
GFace *gfc = current->getFaceByTag(NF+i+1);
- current->remove(gfc);
+ current->remove(gfc);
}
-
+
//Put new mesh in a new discreteFace
for(std::set<MVertex*>::iterator it = allNod.begin(); it != allNod.end(); ++it)
mySplitMesh->mesh_vertices.push_back(*it);
- mySplitMesh->meshStatistics.status = GFace::DONE;
+ mySplitMesh->meshStatistics.status = GFace::DONE;
current->createTopologyFromMesh();
-
+
}
void Centerline::createFaces(){
-
+
std::vector<std::vector<MTriangle*> > faces;
std::multimap<MEdge, MTriangle*, Less_Edge> e2e;
@@ -799,14 +799,14 @@ void Centerline::createClosedVolume(){
GEdge * gec = current->getEdgeByTag(NE+boundEdges[i]->tag());
myEdges.push_back(gec);
myEdgeLoops.push_back(myEdges);
- GFace *newFace = current->addPlanarFace(myEdgeLoops);
+ GFace *newFace = current->addPlanarFace(myEdgeLoops);
newFace->addPhysicalEntity(200);
current->setPhysicalName("in/out", 2, 200);
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 = current->getFaceByTag(NF+i+1);
myFaces.push_back(gf);
@@ -829,7 +829,7 @@ void Centerline::closeVolume(){
}
void Centerline::cutMesh(){
-
+
is_cut = true;
if (update_needed){
@@ -854,16 +854,16 @@ void Centerline::cutMesh(){
// }
Msg::Info("Splitting surface mesh (%d tris) with centerline %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 = (edges[i].minRad+edges[i].maxRad);
double AR = L/D;
- printf("*** Centerline branch %d (AR=%d) \n", i, (int)round(AR));
+ printf("*** Centerline branch %d (AR=%d) \n", i, (int)floor(AR + 0.5));
if( AR > 4.0){
- int nbSplit = (int)round(AR/3.);
+ int nbSplit = (int)floor(AR / 3. + 0.5);
double li = L/nbSplit;
double lc = 0.0;
for (int j= 0; j < lines.size(); j++){
@@ -897,7 +897,7 @@ void Centerline::cutMesh(){
createFaces();
current->createTopologyFromFaces(discFaces);
current->exportDiscreteGEOInternals();
-
+
//write
Msg::Info("Writing splitted mesh 'myPARTS.msh'");
current->writeMSH("myPARTS.msh", 2.2, false, false);
@@ -911,25 +911,25 @@ void Centerline::cutMesh(){
}
void 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();
//printf("cut disk (R=%g)= %g %g %g %g \n", maxRad, a, b, c, d);
-
+
const double EPS = 0.007;
std::set<MEdge,Less_Edge> allEdges;
- for(unsigned int i = 0; i < triangles.size(); i++)
+ for(unsigned int i = 0; i < triangles.size(); i++)
for ( unsigned int j= 0; j < 3; j++)
- allEdges.insert(triangles[i]->getEdge(j));
+ allEdges.insert(triangles[i]->getEdge(j));
bool closedCut = false;
int step = 0;
while (!closedCut && step < 10){
double rad = 1.2*maxRad+0.1*step*maxRad;
std::map<MEdge,MVertex*,Less_Edge> cutEdges;
std::set<MVertex*> cutVertices;
- std::vector<MTriangle*> newTris;
+ std::vector<MTriangle*> newTris;
std::set<MEdge,Less_Edge> newCut;
cutEdges.clear();
cutVertices.clear();
@@ -943,7 +943,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
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;
+ 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);
@@ -955,7 +955,7 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
else if (inters && rdist >= 1.-EPS && inDisk)
cutVertices.insert(me.getVertex(1));
}
- for(unsigned int i = 0; i < triangles.size(); i++){
+ for(unsigned int i = 0; i < triangles.size(); i++){
cutTriangle(triangles[i], cutEdges,cutVertices, newTris, newCut);
}
if (isClosed(newCut)) {
@@ -981,14 +981,14 @@ void Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
// 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++;
+ // itp++;
// }
// fprintf(f2,"};\n");
// fclose(f2);
}
}
-
+
return;
@@ -1004,19 +1004,19 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
buildKdTree();
update_needed = false;
}
-
+
double xyz[3] = {x,y,z};
bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
- std::list<GFace*> cFaces;
+ std::list<GFace*> cFaces;
if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
-
+
//take xyz = closest point on boundary in case we are on the planar in/out faces
if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
(isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
int num_neighbours = 1;
kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
- xyz[0] = nodesR[index[0]][0];
+ xyz[0] = nodesR[index[0]][0];
xyz[1] = nodesR[index[0]][1];
xyz[2] = nodesR[index[0]][2];
}
@@ -1024,7 +1024,7 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
int num_neighbours = 1;
kdtree->annkSearch(xyz, num_neighbours, index, dist);
double d = sqrt(dist[0]);
- double lc = 2*M_PI*d/nbPoints;
+ double lc = 2*M_PI*d/nbPoints;
return lc;
}
@@ -1040,7 +1040,7 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
buildKdTree();
update_needed = false;
}
-
+
//double lc = operator()(x,y,z,ge);
//metr = SMetric3(1./(lc*lc));
@@ -1050,20 +1050,20 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
int num_neighbours = 2;
kdtree->annkSearch(xyz, num_neighbours, index2, dist2);
double d = sqrt(dist2[0]);
- double lc = 2*M_PI*d/nbPoints;
+ double lc = 2*M_PI*d/nbPoints;
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]);
SVector3 dir0 = p1-p0;
SVector3 dir1, dir2;
buildOrthoBasis(dir0,dir1,dir2);
-
+
double lcA = 4.*lc;
double lam1 = 1./(lcA*lcA);
double lam2 = 1./(lc*lc);
metr = SMetric3(lam1,lam2,lam2, dir0, dir1, dir2);
delete[]index2;
- delete[]dist2;
+ delete[]dist2;
return;
}
@@ -1093,7 +1093,7 @@ void Centerline::printSplit() const{
// fprintf(f3, "SP(%g,%g,%g){%g};\n",
// v->x(), v->y(), v->z(),
// (double)v->getNum());
- // itj++;
+ // itj++;
// }
// fprintf(f3,"};\n");
// fclose(f3);
@@ -1110,7 +1110,7 @@ void Centerline::printSplit() const{
}
fprintf(f4,"};\n");
fclose(f4);
-
+
}
#endif
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index c136d27b121b972be06409da0478adc40153744e..54089337372aa7f1019d1f12e174cb3a329ab14f 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -54,7 +54,7 @@ static bool isActive(MTri3 *t, double limit_, int &i, std::set<MEdge,Less_Edge>
MTri3 *neigh = t->getNeigh(i);
if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
int ip1 = i - 1 < 0 ? 2 : i - 1;
- int ip2 = i;
+ int ip2 = i;
MEdge me (t->tri()->getVertex(ip1),t->tri()->getVertex(ip2));
if(front->find(me) != front->end()){
return true;
@@ -72,7 +72,7 @@ static void updateActiveEdges(MTri3 *t, double limit_, std::set<MEdge,Less_Edge>
MTri3 *neigh = t->getNeigh(active);
if (!neigh || (neigh->getRadius() < limit_ && neigh->getRadius() > 0)) {
int ip1 = active - 1 < 0 ? 2 : active - 1;
- int ip2 = active;
+ int ip2 = active;
MEdge me (t->tri()->getVertex(ip1),t->tri()->getVertex(ip2));
front.insert(me);
}
@@ -90,9 +90,9 @@ bool circumCenterMetricInTriangle(MTriangle *base, const double *metric,
}
void circumCenterMetric(double *pa, double *pb, double *pc,
- const double *metric, double *x, double &Radius2)
+ const double *metric, double *x, double &Radius2)
{
- // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
+ // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
double sys[2][2];
double rhs[2];
@@ -106,25 +106,25 @@ void circumCenterMetric(double *pa, double *pb, double *pc,
sys[1][1] = 2. * d * (pa[1] - pc[1]) + 2. * b * (pa[0] - pc[0]);
rhs[0] =
- a * (pa[0] * pa[0] - pb[0] * pb[0]) +
- d * (pa[1] * pa[1] - pb[1] * pb[1]) +
+ a * (pa[0] * pa[0] - pb[0] * pb[0]) +
+ d * (pa[1] * pa[1] - pb[1] * pb[1]) +
2. * b * (pa[0] * pa[1] - pb[0] * pb[1]);
rhs[1] =
- a * (pa[0] * pa[0] - pc[0] * pc[0]) +
- d * (pa[1] * pa[1] - pc[1] * pc[1]) +
+ a * (pa[0] * pa[0] - pc[0] * pc[0]) +
+ d * (pa[1] * pa[1] - pc[1] * pc[1]) +
2. * b * (pa[0] * pa[1] - pc[0] * pc[1]);
if (!sys2x2(sys, rhs, x)){
// printf("%g %g %g\n",a,b,d);
}
- Radius2 =
+ Radius2 =
(x[0] - pa[0]) * (x[0] - pa[0]) * a +
(x[1] - pa[1]) * (x[1] - pa[1]) * d +
2. * (x[0] - pa[0]) * (x[1] - pa[1]) * b;
}
-void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric,
+void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric,
double *res, double *uv, double &radius)
{
double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
@@ -146,14 +146,14 @@ void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric
double mm[3] = {tra(0,0),tra(0,1),tra(1,1)};
circumCenterMetric(p1P, p2P, p3P, mm, resP, radius);
-
+
if(uv){
double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
{p2P[1] - p1P[1], p3P[1] - p1P[1]}};
double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
sys2x2(mat, rhs, uv);
}
-
+
res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
@@ -162,9 +162,9 @@ void circumCenterMetricXYZ(double *p1, double *p2, double *p3, SMetric3 & metric
void circumCenterMetric(MTriangle *base, const double *metric,
const std::vector<double> &Us,
const std::vector<double> &Vs,
- double *x, double &Radius2)
+ double *x, double &Radius2)
{
- // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
+ // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
double pa[2] = {Us[base->getVertex(0)->getIndex()],
Vs[base->getVertex(0)->getIndex()]};
double pb[2] = {Us[base->getVertex(1)->getIndex()],
@@ -177,15 +177,15 @@ void circumCenterMetric(MTriangle *base, const double *metric,
void buildMetric(GFace *gf, double *uv, double *metric)
{
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
-
+
metric[0] = dot(der.first(), der.first());
metric[1] = dot(der.second(), der.first());
metric[2] = dot(der.second(), der.second());
-}
+}
-// m 3x3
+// m 3x3
// d 3x2
-// M = d^T m d
+// M = d^T m d
void buildMetric(GFace *gf, double *uv, SMetric3 & m, double *metric)
{
@@ -214,10 +214,10 @@ void buildMetric(GFace *gf, double *uv, SMetric3 & m, double *metric)
metric[0] = dot(x1, der.first());
metric[1] = dot(x2, der.first());
metric[2] = dot(x2, der.second());
-}
+}
-int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
- double *uv, double *metric)
+int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
+ double *uv, double *metric)
{
double x[2], Radius2;
circumCenterMetric(p1, p2, p3, metric, x, Radius2);
@@ -227,22 +227,22 @@ int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3,
double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
+ (x[1] - uv[1]) * (x[1] - uv[1]) * d
+ 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
- return d2 < Radius2;
+ return d2 < Radius2;
}
-int inCircumCircleAniso(GFace *gf, MTriangle *base,
+int inCircumCircleAniso(GFace *gf, MTriangle *base,
const double *uv, const double *metricb,
const std::vector<double> &Us,
- const std::vector<double> &Vs)
+ const std::vector<double> &Vs)
{
double x[2], Radius2, metric[3];
- double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
Us[base->getVertex(1)->getIndex()] +
Us[base->getVertex(2)->getIndex()]) / 3.,
(Vs[base->getVertex(0)->getIndex()] +
- Vs[base->getVertex(1)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
Vs[base->getVertex(2)->getIndex()]) / 3.};
-
+
buildMetric(gf, pa, metric);
circumCenterMetric(base, metric, Us, Vs, x, Radius2);
@@ -253,11 +253,11 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
+ (x[1] - uv[1]) * (x[1] - uv[1]) * d
+ 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
-
- return d2 < Radius2;
+
+ return d2 < Radius2;
}
-MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double> *Us, const std::vector<double> *Vs, GFace *gf)
+MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double> *Us, const std::vector<double> *Vs, GFace *gf)
: deleted(false), base(t)
{
neigh[0] = neigh[1] = neigh[2] = 0;
@@ -270,13 +270,13 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double
{base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
// compute the circumradius of the triangle
-
+
if (!metric){
if (radiusNorm == 2){
circumCenterXYZ(pa, pb, pc, center);
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
- const double dz = base->getVertex(0)->z() - center[2];
+ const double dz = base->getVertex(0)->z() - center[2];
circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
circum_radius /= lc;
}
@@ -287,33 +287,33 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double
(*Vs)[base->getVertex(1)->getIndex()]};
double p3[2] = {(*Us)[base->getVertex(2)->getIndex()],
(*Vs)[base->getVertex(2)->getIndex()]};
-
+
double midpoint[2] = {(p1[0]+p2[0]+p3[0])/3.0,(p1[1]+p2[1]+p3[1])/3.0};
- double quadAngle = backgroundMesh::current() ? backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) : 0.0;
+ double quadAngle = backgroundMesh::current() ? backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) : 0.0;
double x0 = p1[0] * cos(quadAngle) + p1[1] * sin(quadAngle);
- double y0 = -p1[0] * sin(quadAngle) + p1[1] * cos(quadAngle);
+ double y0 = -p1[0] * sin(quadAngle) + p1[1] * cos(quadAngle);
double x1 = p2[0] * cos(quadAngle) + p2[1] * sin(quadAngle);
- double y1 = -p2[0] * sin(quadAngle) + p2[1] * cos(quadAngle);
+ double y1 = -p2[0] * sin(quadAngle) + p2[1] * cos(quadAngle);
double x2 = p3[0] * cos(quadAngle) + p3[1] * sin(quadAngle);
- double y2 = -p3[0] * sin(quadAngle) + p3[1] * cos(quadAngle);
+ double y2 = -p3[0] * sin(quadAngle) + p3[1] * cos(quadAngle);
double xmax = std::max(std::max(x0,x1),x2);
double ymax = std::max(std::max(y0,y1),y2);
double xmin = std::min(std::min(x0,x1),x2);
double ymin = std::min(std::min(y0,y1),y2);
-
+
double metric[3];
buildMetric(gf, midpoint, metric);
double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
circum_radius = std::max(xmax-xmin,ymax-ymin) / RATIO;
- circum_radius /= lc;
+ circum_radius /= lc;
}
}
else{
double uv[2];
- circumCenterMetricXYZ(pa, pb, pc, *metric, center, uv, circum_radius);
+ circumCenterMetricXYZ(pa, pb, pc, *metric, center, uv, circum_radius);
}
}
@@ -329,14 +329,14 @@ int MTri3::inCircumCircle(const double *p) const
double fourth[3];
fourthPoint(pa, pb, pc, fourth);
- double result = robustPredicates::insphere(pa, pb, pc, fourth, (double*)p) *
- robustPredicates::orient3d(pa, pb, pc,fourth);
- return (result > 0) ? 1 : 0;
+ double result = robustPredicates::insphere(pa, pb, pc, fourth, (double*)p) *
+ robustPredicates::orient3d(pa, pb, pc,fourth);
+ return (result > 0) ? 1 : 0;
}
-int inCircumCircle(MTriangle *base,
+int inCircumCircle(MTriangle *base,
const double *p, const double *param ,
- std::vector<double> &Us, std::vector<double> &Vs)
+ std::vector<double> &Us, std::vector<double> &Vs)
{
double pa[2] = {Us[base->getVertex(0)->getIndex()],
Vs[base->getVertex(0)->getIndex()]};
@@ -345,9 +345,9 @@ int inCircumCircle(MTriangle *base,
double pc[2] = {Us[base->getVertex(2)->getIndex()],
Vs[base->getVertex(2)->getIndex()]};
- double result = robustPredicates::incircle(pa, pb, pc, (double*)param) *
- robustPredicates::orient2d(pa, pb, pc);
- return (result > 0) ? 1 : 0;
+ double result = robustPredicates::incircle(pa, pb, pc, (double*)param) *
+ robustPredicates::orient2d(pa, pb, pc);
+ return (result > 0) ? 1 : 0;
}
template <class ITER>
@@ -410,7 +410,7 @@ void connectTriangles(std::set<MTri3*, compareTri3Ptr> &l)
connectTris(l.begin(), l.end());
}
-void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
+void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
double *v, double *param, MTri3 *t,
std::vector<double> &Us, std::vector<double> &Vs)
{
@@ -434,7 +434,7 @@ void recurFindCavity(std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
}
void recurFindCavityAniso(GFace *gf,
- std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
+ std::list<edgeXface> &shell, std::list<MTri3*> &cavity,
double *metric, double *param, MTri3 *t,
std::vector<double> &Us, std::vector<double> &Vs)
{
@@ -496,14 +496,14 @@ bool invMapUV(MTriangle *t, double *p,
b[1] = p[1] - v0;
sys2x2(mat, b, uv);
- if(uv[0] >= -tol &&
- uv[1] >= -tol &&
- uv[0] <= 1. + tol &&
- uv[1] <= 1. + tol &&
+ if(uv[0] >= -tol &&
+ uv[1] >= -tol &&
+ uv[0] <= 1. + tol &&
+ uv[1] <= 1. + tol &&
1. - uv[0] - uv[1] > -tol) {
return true;
}
- return false;
+ return false;
}
double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
@@ -516,22 +516,22 @@ double getSurfUV(MTriangle *t, std::vector<double> &Us, std::vector<double> &Vs)
double v3 = Vs[t->getVertex(2)->getIndex()];
const double vv1[2] = {u2 - u1, v2 - v1};
const double vv2[2] = {u3 - u1, v3 - v1};
- double s = vv1[0] * vv2[1] - vv1[1] * vv2[0];
+ double s = vv1[0] * vv2[1] - vv1[1] * vv2[0];
return s * 0.5;
}
bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
std::set<MTri3*, compareTri3Ptr> &allTets,
std::set<MTri3*, compareTri3Ptr> *activeTets,
- std::vector<double> &vSizes,
+ std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
std::vector<SMetric3> &vMetricsBGM,
- std::vector<double> &Us,
+ std::vector<double> &Us,
std::vector<double> &Vs,
double *metric = 0)
{
std::list<edgeXface> shell;
- std::list<MTri3*> cavity;
+ std::list<MTri3*> cavity;
std::list<MTri3*> new_cavity;
if (!metric){
@@ -539,10 +539,10 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
recurFindCavity(shell, cavity, p, param, t, Us, Vs);
}
else{
- recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
+ recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
// printf("RECURSIVELY FIND A CAVITY %d %d \n",shell.size(),cavity.size());
}
-
+
// check that volume is conserved
double newVolume = 0;
double oldVolume = 0;
@@ -553,7 +553,7 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
oldVolume += fabs(getSurfUV((*ittet)->tri(),Us,Vs));
++ittet;
}
-
+
MTri3 **newTris = new MTri3*[shell.size()];
int k = 0;
std::list<edgeXface>::iterator it = shell.begin();
@@ -571,15 +571,15 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
double LL = Extend1dMeshIn2dSurfaces() ? std::min(lc, lcBGM) : lcBGM;
MTri3 *t4;
- t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
-
+ t4 = new MTri3(t, LL,0,&Us,&Vs,gf);
+
double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
(it->v[0]->y() - v->y()) * (it->v[0]->y() - v->y()) +
(it->v[0]->z() - v->z()) * (it->v[0]->z() - v->z()));
double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
(it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
(it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
- const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),0.5*(it->v[0]->y()+it->v[1]->y()),0.5*(it->v[0]->z()+it->v[1]->z())};
+ const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),0.5*(it->v[0]->y()+it->v[1]->y()),0.5*(it->v[0]->z()+it->v[1]->z())};
double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) + (MID[1] - v->y()) * (MID[1] - v->y()) + (MID[2] - v->z()) * (MID[2] - v->z()));
if ((d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) && !force) {
// printf("TOO CLOSE %g %g %g %g %g %g\n",d1,d2,d3,LL,lc,lcBGM);
@@ -600,9 +600,9 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
newVolume += ss;
++it;
}
- if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() >= 3 &&
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume && shell.size() >= 3 &&
!onePointIsTooClose){
- connectTris(new_cavity.begin(), new_cavity.end());
+ connectTris(new_cavity.begin(), new_cavity.end());
allTets.insert(newTris, newTris + shell.size());
// printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
if (activeTets){
@@ -618,13 +618,13 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
return true;
}
// The cavity is NOT star shaped
- else{
+ else{
// printf("cavity not star shaped %22.15E vs %22.15E\n",oldVolume,newVolume);
// printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
- delete [] newTris;
+ delete [] newTris;
ittet = cavity.begin();
- ittete = cavity.end();
+ ittete = cavity.end();
while(ittet != ittete){
(*ittet)->setDeleted(false);
++ittet;
@@ -678,11 +678,11 @@ void _printTris(char *name, std::set<MTri3*, compareTri3Ptr> &AllTris,
fclose (ff);
}
-static MTri3* search4Triangle (MTri3 *t, double pt[2],
+static MTri3* search4Triangle (MTri3 *t, double pt[2],
std::vector<double> &Us, std::vector<double> &Vs,
std::set<MTri3*,compareTri3Ptr> &AllTris, double uv[2]) {
-
- bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
+
+ bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
if (inside) return t;
SPoint3 q1(pt[0],pt[1],0);
int ITER = 0;
@@ -699,16 +699,16 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
if (intersection_segments(p1,p2,q1,q2,xcc))break;
}
if (i>=3)break;
- t = t->getNeigh(i);
+ t = t->getNeigh(i);
if (!t)break;
- bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
+ bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
if (inside) {return t;}
if (ITER++ > AllTris.size())break;
}
return 0;
for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin(); itx != AllTris.end();++itx){
if (!(*itx)->isDeleted()){
- inside = invMapUV((*itx)->tri(), pt, Us, Vs, uv, 1.e-8);
+ inside = invMapUV((*itx)->tri(), pt, Us, Vs, uv, 1.e-8);
if (inside){
return *itx;
}
@@ -718,9 +718,9 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
}
static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it,
- double center[2], double metric[3],
+ double center[2], double metric[3],
std::vector<double> &Us, std::vector<double> &Vs,
- std::vector<double> &vSizes,
+ std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
std::vector<SMetric3> &vMetricsBGM,
std::set<MTri3*,compareTri3Ptr> &AllTris,
@@ -741,7 +741,7 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
double uv[2];
// FIXME !!! ----> FIXED (JFR)
if (MTri3::radiusNorm == 2){
- inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
+ inside = invMapUV(worst->tri(), center, Us, Vs, uv, 1.e-8);
if (inside)ptin = worst;
if (!inside && worst->getNeigh(0)){
inside |= invMapUV(worst->getNeigh(0)->tri(), center, Us, Vs, uv, 1.e-8);
@@ -779,13 +779,13 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
v->setIndex(Us.size());
double lc1,lc;
if (backgroundMesh::current()){
- lc1 = lc =
+ lc1 = lc =
backgroundMesh::current()->operator()(center[0], center[1], 0.0);
}
else {
- lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
- uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
- uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
+ lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getIndex()] +
+ uv[0] * vSizes[ptin->tri()->getVertex(1)->getIndex()] +
+ uv[1] * vSizes[ptin->tri()->getVertex(2)->getIndex()]);
lc = BGM_MeshSize(gf, center[0], center[1], p.x(), p.y(), p.z());
}
//SMetric3 metr = BGM_MeshMetric(gf, center[0], center[1], p.x(), p.y(), p.z());
@@ -794,18 +794,18 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
vSizes.push_back(lc1);
Us.push_back(center[0]);
Vs.push_back(center[1]);
-
+
if(!p.succeeded() || !insertVertex
- (false, gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
+ (false, gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
-
+
MTriangle *base = worst->tri();
-
+
// Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
-
+
AllTris.erase(it);
worst->forceRadius(-1);
- AllTris.insert(worst);
+ AllTris.insert(worst);
delete v;
return false;
}
@@ -816,15 +816,15 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
else {
- MTriangle *base = worst->tri();
- /*
+ MTriangle *base = worst->tri();
+ /*
Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
center[0], center[1],
- Us[base->getVertex(0)->getIndex()],
- Vs[base->getVertex(0)->getIndex()],
- Us[base->getVertex(1)->getIndex()],
- Vs[base->getVertex(1)->getIndex()],
- Us[base->getVertex(2)->getIndex()],
+ Us[base->getVertex(0)->getIndex()],
+ Vs[base->getVertex(0)->getIndex()],
+ Us[base->getVertex(1)->getIndex()],
+ Vs[base->getVertex(1)->getIndex()],
+ Us[base->getVertex(2)->getIndex()],
Vs[base->getVertex(2)->getIndex()]);
*/
AllTris.erase(it);
@@ -875,18 +875,18 @@ void bowyerWatson(GFace *gf)
if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0)) break;
circUV(worst->tri(), Us, Vs, center, gf);
MTriangle *base = worst->tri();
- double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
- Us[base->getVertex(1)->getIndex()] +
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
Us[base->getVertex(2)->getIndex()]) / 3.,
- (Vs[base->getVertex(0)->getIndex()] +
- Vs[base->getVertex(1)->getIndex()] +
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
Vs[base->getVertex(2)->getIndex()]) / 3.};
buildMetric(gf, pa, metric);
- circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
- insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes,
+ circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+ insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris);
}
- }
+ }
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
@@ -896,7 +896,7 @@ void bowyerWatson(GFace *gf)
if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
}
}
- transferDataStructure(gf, AllTris, Us, Vs);
+ transferDataStructure(gf, AllTris, Us, Vs);
}
/*
@@ -907,7 +907,7 @@ void bowyerWatson(GFace *gf)
The point isertion is done on the Voronoi Edges
*/
-double lengthInfniteNorm(const double p[2], const double q[2],
+double lengthInfniteNorm(const double p[2], const double q[2],
const double quadAngle){
double xp = p[0] * cos(quadAngle) + p[1] * sin(quadAngle);
double yp = -p[0] * sin(quadAngle) + p[1] * cos(quadAngle);
@@ -920,7 +920,7 @@ double lengthInfniteNorm(const double p[2], const double q[2],
return std::max(xmax-xmin,ymax-ymin);
}
-void circumCenterInfinite (MTriangle *base, double quadAngle,
+void circumCenterInfinite (MTriangle *base, double quadAngle,
const std::vector<double> &Us,
const std::vector<double> &Vs, double *x) {
double pa[2] = {Us[base->getVertex(0)->getIndex()],
@@ -944,7 +944,7 @@ void circumCenterInfinite (MTriangle *base, double quadAngle,
}
-static double lengthMetric(const double p[2], const double q[2],
+static double lengthMetric(const double p[2], const double q[2],
const double metric[3])
{
return sqrt ( (p[0] - q[0]) * metric[0] * (p[0] - q[0]) +
@@ -959,118 +959,118 @@ static double lengthMetric(const double p[2], const double q[2],
/ |
lc / | r
/ |
- / |
+ / |
/ x
- / |
+ / |
/ | r/2
-/ |
+/ |
-----------+
lc/2
- (3 r/2)^2 = lc^2 - lc^2/4
+ (3 r/2)^2 = lc^2 - lc^2/4
-> lc^2 3/4 = 9r^2/4
-> lc^2 = 3 r^2
r^2 /4 + lc^2/4 = r^2
-> lc^2 = 3 r^2
-
+
*/
-double optimalPointFrontal (GFace *gf,
- MTri3* worst,
+double optimalPointFrontal (GFace *gf,
+ MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM,
+ std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
double center[2],r2;
MTriangle *base = worst->tri();
circUV(base, Us, Vs, center, gf);
- double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
- Us[base->getVertex(1)->getIndex()] +
+ double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
+ Us[base->getVertex(1)->getIndex()] +
Us[base->getVertex(2)->getIndex()]) / 3.,
- (Vs[base->getVertex(0)->getIndex()] +
- Vs[base->getVertex(1)->getIndex()] +
+ (Vs[base->getVertex(0)->getIndex()] +
+ Vs[base->getVertex(1)->getIndex()] +
Vs[base->getVertex(2)->getIndex()]) / 3.};
buildMetric(gf, pa, metric);
- circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
+ circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
// compute the middle point of the edge
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
int ip3 = (active_edge+1)%3;
- double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
Vs[base->getVertex(ip1)->getIndex()]};
- double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
Vs[base->getVertex(ip2)->getIndex()]};
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
-
- // now we have the edge center and the center of the circumcircle,
+
+ // now we have the edge center and the center of the circumcircle,
// we try to find a point that would produce a perfect triangle while
// connecting the 2 points of the active edge
-
+
double dir[2] = {center[0] - midpoint[0], center[1] - midpoint[1]};
double norm = sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
dir[0] /= norm;
dir[1] /= norm;
const double RATIO = sqrt(dir[0] * dir[0] * metric[0] +
2 * dir[1] * dir[0] * metric[1] +
- dir[1] * dir[1] * metric[2]);
-
+ dir[1] * dir[1] * metric[2]);
+
// const double p = 0.5 * lengthMetric(P, Q, metric); // / RATIO;
// const double q = lengthMetric(center, midpoint, metric);
/*
- const double rhoM1 = 0.5 *
- (vSizes[base->getVertex(ip1)->getIndex()] +
+ const double rhoM1 = 0.5 *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
- const double rhoM2 = 0.5 *
- (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ const double rhoM2 = 0.5 *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
*/
-
+
// const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
// double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
-
+
//const double d = 100./RATIO;
// printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
// printf("size %12.5E\n",vSizesBGM[base->getVertex(ip1)->getIndex()]);
const double rhoM1 = 0.5*
- (vSizes[base->getVertex(ip1)->getIndex()] +
+ (vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
- const double rhoM2 = 0.5*
- (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ const double rhoM2 = 0.5*
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1,rhoM2) : rhoM2;
const double rhoM_hat = rhoM;
- const double d = rhoM_hat * sqrt(3)*0.5/RATIO;
-
+ const double d = rhoM_hat * sqrt(3.)*0.5/RATIO;
+
// printf("%12.5E %12.5E\n",d,RATIO);
- newPoint[0] = midpoint[0] + d * dir[0];
+ newPoint[0] = midpoint[0] + d * dir[0];
newPoint[1] = midpoint[1] + d * dir[1];
return d * RATIO;
}
/*
x
- |
- |
- | d = 3^{1/2}/2 h
- |
- |
+ |
+ |
+ | d = 3^{1/2}/2 h
+ |
+ |
------p-------> n
h
-
+
x point of the plane
h being some kind of average between the size field
and the edge length
*/
-struct intersectCurveSurfaceData
+struct intersectCurveSurfaceData
{
GFace *gf;
SVector3 &n1,&n2;
@@ -1081,7 +1081,7 @@ struct intersectCurveSurfaceData
{ }
};
-static void intersectCircleSurface(fullVector<double> &uvt,
+static void intersectCircleSurface(fullVector<double> &uvt,
fullVector<double> &res, void *_data){
intersectCurveSurfaceData *data = (intersectCurveSurfaceData*)_data;
GPoint gp = data->gf->point(uvt(0),uvt(1));
@@ -1094,13 +1094,13 @@ static void intersectCircleSurface(fullVector<double> &uvt,
}
-void optimalPointFrontalB (GFace *gf,
- MTri3* worst,
+void optimalPointFrontalB (GFace *gf,
+ MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM,
+ std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
// as a starting point, let us use the "fast algo"
@@ -1122,7 +1122,7 @@ void optimalPointFrontalB (GFace *gf,
SVector3 n2 = crossprod(v1v2,n1);
n1.normalize();
n2.normalize();
- // we look for a point that is
+ // we look for a point that is
// P = d * (n1 cos(t) + n2 sin(t)) that is on the surface
// so we have to find t, starting with t = 0
fullVector<double> uvt(3);
@@ -1153,7 +1153,7 @@ void bowyerWatsonFrontal(GFace *gf)
// delaunise the initial mesh
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
-
+
int ITER = 0, active_edge;
// compute active triangle
std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
@@ -1163,7 +1163,7 @@ void bowyerWatsonFrontal(GFace *gf)
else if ((*it)->getRadius() < LIMIT_)break;
}
-
+
// insert points
while (1){
/*
@@ -1180,7 +1180,7 @@ void bowyerWatsonFrontal(GFace *gf)
ActiveTris.erase(ActiveTris.begin());
// printf("active_tris.size = %d\n",ActiveTris.size());
- if (!worst->isDeleted() && isActive(worst, LIMIT_, active_edge) &&
+ if (!worst->isDeleted() && isActive(worst, LIMIT_, active_edge) &&
worst->getRadius() > LIMIT_){
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
@@ -1189,8 +1189,8 @@ void bowyerWatsonFrontal(GFace *gf)
optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
insertAPoint(gf, AllTris.end(), newPoint, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
- }
-
+ }
+
/* if(ITER % 1== 0){
char name[245];
sprintf(name,"frontal%d-ITER%d.pos",gf->tag(),ITER);
@@ -1204,8 +1204,8 @@ void bowyerWatsonFrontal(GFace *gf)
// _printTris (name, AllTris, Us, Vs,false);
// sprintf(name,"frontal%d-param.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,true);
- transferDataStructure(gf, AllTris, Us, Vs);
- // in case of boundary layer meshing
+ transferDataStructure(gf, AllTris, Us, Vs);
+ // in case of boundary layer meshing
{
FieldManager *fields = gf->model()->getFields();
BoundaryLayerField *blf = 0;
@@ -1215,40 +1215,40 @@ void bowyerWatsonFrontal(GFace *gf)
if (blf && !blf->iRecombine)quadsToTriangles(gf,10000);
}
}
-}
+}
-void optimalPointFrontalQuad (GFace *gf,
- MTri3* worst,
+void optimalPointFrontalQuad (GFace *gf,
+ MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM,
+ std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
MTriangle *base = worst->tri();
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
int ip3 = (active_edge+1)%3;
-
- double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
Vs[base->getVertex(ip1)->getIndex()]};
- double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
Vs[base->getVertex(ip2)->getIndex()]};
- double O[2] = {Us[base->getVertex(ip3)->getIndex()],
+ double O[2] = {Us[base->getVertex(ip3)->getIndex()],
Vs[base->getVertex(ip3)->getIndex()]};
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
-
+
// compute background mesh data
double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
double center[2];
- circumCenterInfinite (base, quadAngle,Us,Vs,center);
-
+ circumCenterInfinite (base, quadAngle,Us,Vs,center);
+
// rotate the points with respect to the angle
double XP1 = 0.5*(Q[0] - P[0]);
double YP1 = 0.5*(Q[1] - P[1]);
- double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
- double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
+ double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
+ double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
// ensure xp > yp
bool exchange = false;
if (fabs(xp) < fabs(yp)){
@@ -1257,23 +1257,23 @@ void optimalPointFrontalQuad (GFace *gf,
yp = temp;
exchange = true;
}
-
+
buildMetric(gf, midpoint, metric);
double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
-
- const double p = 0.5 * lengthInfniteNorm(P, Q, quadAngle);
+
+ const double p = 0.5 * lengthInfniteNorm(P, Q, quadAngle);
const double q = lengthInfniteNorm(center, midpoint, quadAngle);
- const double rhoM1 = 0.5 * RATIO *
- (vSizes[base->getVertex(ip1)->getIndex()] +
+ const double rhoM1 = 0.5 * RATIO *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
- const double rhoM2 = 0.5 * RATIO *
- (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ const double rhoM2 = 0.5 * RATIO *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] ) / sqrt(3.);// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
-
+
const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
const double factor = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) /(sqrt(3.)*p);
-
+
double npx,npy;
if (xp*yp > 0){
npx = - fabs(xp)*factor;
@@ -1287,9 +1287,9 @@ void optimalPointFrontalQuad (GFace *gf,
double temp = npx;
npx = npy;
npy = temp;
- }
-
-
+ }
+
+
newPoint[0] = midpoint[0] + cos(quadAngle) * npx - sin(quadAngle) * npy;
newPoint[1] = midpoint[1] + sin(quadAngle) * npx + cos(quadAngle) * npy;
@@ -1297,36 +1297,36 @@ void optimalPointFrontalQuad (GFace *gf,
(midpoint[1] - newPoint[1])*(midpoint[1] - O[1]) < 0){
newPoint[0] = midpoint[0] - cos(quadAngle) * npx + sin(quadAngle) * npy;
newPoint[1] = midpoint[1] - sin(quadAngle) * npx - cos(quadAngle) * npy;
- }
-}
+ }
+}
-void optimalPointFrontalQuadAniso (GFace *gf,
- MTri3* worst,
+void optimalPointFrontalQuadAniso (GFace *gf,
+ MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM,
+ std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
MTriangle *base = worst->tri();
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
int ip3 = (active_edge+1)%3;
-
- double P[2] = {Us[base->getVertex(ip1)->getIndex()],
+
+ double P[2] = {Us[base->getVertex(ip1)->getIndex()],
Vs[base->getVertex(ip1)->getIndex()]};
- double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
+ double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
Vs[base->getVertex(ip2)->getIndex()]};
- double O[2] = {Us[base->getVertex(ip3)->getIndex()],
+ double O[2] = {Us[base->getVertex(ip3)->getIndex()],
Vs[base->getVertex(ip3)->getIndex()]};
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
-
+
// compute background mesh data
double quadAngle = backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0);
double center[2];
- circumCenterInfinite (base, quadAngle,Us,Vs,center);
-
+ circumCenterInfinite (base, quadAngle,Us,Vs,center);
+
// rotate the points with respect to the angle
double XP1 = 0.5*(Q[0] - P[0]);
double YP1 = 0.5*(Q[1] - P[1]);
@@ -1334,14 +1334,14 @@ void optimalPointFrontalQuadAniso (GFace *gf,
double DX = 0.5*(Q[0] + P[0]);
double DY = 0.5*(Q[0] + P[0]);
- double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
- double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
+ double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
+ double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
double X4 = O[0] -DX;
double Y4 = O[1] -DY;
- double x4 = X4 * cos(quadAngle) + Y4 * sin(quadAngle);
- double y4 = -X4 * sin(quadAngle) + Y4 * cos(quadAngle);
+ double x4 = X4 * cos(quadAngle) + Y4 * sin(quadAngle);
+ double y4 = -X4 * sin(quadAngle) + Y4 * cos(quadAngle);
// ensure xp > yp
bool exchange = false;
@@ -1354,18 +1354,18 @@ void optimalPointFrontalQuadAniso (GFace *gf,
y4 = temp;
exchange = true;
}
-
+
buildMetric(gf, midpoint, metric);
- double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
- const double rhoM1 = 0.5 * RATIO *
- (vSizes[base->getVertex(ip1)->getIndex()] +
+ double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
+ const double rhoM1 = 0.5 * RATIO *
+ (vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] );
- const double rhoM2 = 0.5 * RATIO *
- (vSizesBGM[base->getVertex(ip1)->getIndex()] +
+ const double rhoM2 = 0.5 * RATIO *
+ (vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] );// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
-
-
+
+
double npx,npy;
const double L_edge = lengthInfniteNorm(P,Q,quadAngle);
@@ -1377,10 +1377,10 @@ void optimalPointFrontalQuadAniso (GFace *gf,
double xl[2] = {rhoM, rhoM+fabs(xp)-fabs(yp)};
const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
- const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
+ const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
double t = ( L_edge - rhoM)/(L_edge - R_Active);
t = std::max(1.,t);
- t = std::min(-L_el/L_ec,t);
+ t = std::min(-L_el/L_ec,t);
npx = xe[0] + t*(xc[0]-xe[0]);
npy = xe[1] + t*(xc[1]-xe[1]);
}
@@ -1390,11 +1390,11 @@ void optimalPointFrontalQuadAniso (GFace *gf,
double xl[2] = {rhoM, rhoM+xp-yp};
const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
- const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
+ const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
const double XP[2]={xp,yp};
double t = ( L_edge - rhoM)/(L_edge - R_Active);
t = std::max(1.,t);
- t = std::min(-L_el/L_ec,t);
+ t = std::min(-L_el/L_ec,t);
npx = xe[0] + t*(xc[0]-xe[0]);
npy = xe[1] + t*(xc[1]-xe[1]);
}
@@ -1413,9 +1413,9 @@ void optimalPointFrontalQuadAniso (GFace *gf,
double temp = npx;
npx = npy;
npy = temp;
- }
-
-
+ }
+
+
newPoint[0] = midpoint[0] + cos(quadAngle) * npx - sin(quadAngle) * npy;
newPoint[1] = midpoint[1] + sin(quadAngle) * npx + cos(quadAngle) * npy;
@@ -1423,8 +1423,8 @@ void optimalPointFrontalQuadAniso (GFace *gf,
(midpoint[1] - newPoint[1])*(midpoint[1] - O[1]) < 0){
newPoint[0] = midpoint[0] - cos(quadAngle) * npx + sin(quadAngle) * npy;
newPoint[1] = midpoint[1] - sin(quadAngle) * npx - cos(quadAngle) * npy;
- }
-}
+ }
+}
void buildBackGroundMesh (GFace *gf) {
@@ -1448,15 +1448,15 @@ void buildBackGroundMesh (GFace *gf) {
gf->quadrangles[i]->getVertex(3)));
}
*/
- // avoid computing curvatures on the fly : only on the
+ // avoid computing curvatures on the fly : only on the
// BGM computes once curvatures at each node
// Disable curvature control
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
- CTX::instance()->mesh.lcFromCurvature = 0;
+ CTX::instance()->mesh.lcFromCurvature = 0;
// Do a background mesh
bowyerWatson(gf);
// Re-enable curv control if asked
- CTX::instance()->mesh.lcFromCurvature = CurvControl;
+ CTX::instance()->mesh.lcFromCurvature = CurvControl;
// apply this to the BGM
// printf("1 end build bak mesh\n");
backgroundMesh::set(gf);
@@ -1470,7 +1470,7 @@ void buildBackGroundMesh (GFace *gf) {
}
gf->triangles = TR;
// gf->quadrangles = QR;
- }
+ }
// printf("end build bak mesh\n");
}
@@ -1488,14 +1488,14 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
//LIMIT_ = 4./3.;//sqrt(2.) * .99;
MTri3::radiusNorm =-1;
}
-
+
buildMeshGenerationDataStructures
(gf, AllTris, vSizes, vSizesBGM, vMetricsBGM,Us, Vs);
// delaunise the initial mesh
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, Us, Vs, vSizes, vSizesBGM);
Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
-
+
int ITER = 0, active_edge;
// compute active triangle
std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
@@ -1519,42 +1519,42 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
_printTris (name, AllTris, Us,Vs,true);
sprintf(name,"denInfinite_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
_printTris (name, ActiveTris, Us,Vs,true);
- }
-
+ }
+
std::set<MTri3*,compareTri3Ptr> ActiveTrisNotInFront;
// printf("%d active triangles\n",ActiveTris.size());
while (1){
-
+
if (!ActiveTris.size())break;
-
+
// if (gf->tag() == 1900){ char name[245];
// sprintf(name,"x_GFace_%d_Layer_%d.pos",gf->tag(),ITER);
// _printTris (name, AllTris, Us,Vs,true);
// }
std::set<MTri3*,compareTri3Ptr>::iterator WORST_ITER = ActiveTris.begin();
-
+
MTri3 *worst = (*WORST_ITER);
ActiveTris.erase(WORST_ITER);
- if (!worst->isDeleted() && (ITERATION > max_layers ? isActive(worst, LIMIT_, active_edge) : isActive(worst, LIMIT_, active_edge,&_front) ) &&
+ if (!worst->isDeleted() && (ITERATION > max_layers ? isActive(worst, LIMIT_, active_edge) : isActive(worst, LIMIT_, active_edge,&_front) ) &&
worst->getRadius() > LIMIT_){
// for (active_edge = 0 ; active_edge < 0 ; active_edge ++){
- // if (active_edges[active_edge])break;
+ // if (active_edges[active_edge])break;
// }
// Msg::Info("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
// vSizes.size(), worst->getRadius(),_front.size());
if(ITER++ % 5000 == 0)
Msg::Debug("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
vSizes.size(), worst->getRadius(),_front.size());
-
+
// compute the middle point of the edge
double newPoint[2],metric[3]={1,0,1};
if (quad)optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
else optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
-
-
+
+
// printf("start INSERT A POINT %g %g \n",newPoint[0],newPoint[1]);
insertAPoint(gf, AllTris.end(), newPoint, 0, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
@@ -1562,7 +1562,7 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
// ActiveTrisNotInFront.insert(worst);
// }
// printf("-----------------> size %d\n",AllTris.size());
-
+
/*
if(ITER % 1== 0){
char name[245];
@@ -1588,16 +1588,16 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
// Msg::Info("%d active tris %d front edges %d not in front",ActiveTris.size(),_front.size(),ActiveTrisNotInFront.size());
if (!ActiveTris.size())break;
}
-
+
char name[245];
// sprintf(name,"frontal%d-real.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,false);
// sprintf(name,"frontal%d-param.pos", gf->tag());
// _printTris (name, AllTris, Us, Vs,true);
- transferDataStructure(gf, AllTris, Us, Vs);
+ transferDataStructure(gf, AllTris, Us, Vs);
MTri3::radiusNorm = 2;
LIMIT_ = 0.5 * sqrt(2.0) * 1;
backgroundMesh::unset();
-}
+}
diff --git a/Post/PView.cpp b/Post/PView.cpp
index f89816ca6996300eb1511c8042961d68dacd9511..fa34cea63d71d2dab9ee9d0c537af79e16e422ba 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -90,7 +90,7 @@ PView::PView(const std::string &xname, const std::string &yname,
_data->setFileName(yname + ".pos");
_options = new PViewOptions(PViewOptions::reference);
_options->type = PViewOptions::Plot2D;
- _options->axes = 2;
+ _options->axes = 3;
//_options->lineWidth = 2.;
//_options->pointSize = 4.;
_options->axesLabel[0] = xname;