diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a2fbfc862e261dd431ff9ce48f184e9cec58fff1..a9631630a9419b27588a58ac802d4616c16e9215 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -45,6 +45,9 @@
#include "MPoint.h"
#include "Numeric.h"
#include "meshGFace.h"
+#if defined(HAVE_ANN)
+#include <ANN/ANN.h>
+#endif
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
@@ -936,6 +939,7 @@ bool GFaceCompound::one2OneMap() const
bool GFaceCompound::parametrize() const
{
+
if (_compound.size() > 1) coherencePatches();
bool paramOK = true;
@@ -1236,9 +1240,11 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfCompound toc,
int allowPartition,
linearSystem<double>* lsys)
- : GFace(m, tag), _compound(compound), oct(0), _U0(U0),
+ : GFace(m, tag), _compound(compound), oct(0), octNew(0), _U0(U0),
_toc(toc), _allowPartition(allowPartition), _lsys(lsys)
{
+
+
ONE = new simpleFunction<double>(1.0);
MONE = new simpleFunction<double>(-1.0);
@@ -1272,6 +1278,12 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
nbSplit = 0;
fillTris.clear();
+
+#if defined(HAVE_ANN)
+ index = new ANNidx[2];
+ dist = new ANNdist[2];
+#endif
+
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
@@ -1279,7 +1291,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
typeOfCompound toc,
int allowPartition,
linearSystem<double>* lsys)
- : GFace(m, tag), _compound(compound), oct(0),
+ : GFace(m, tag), _compound(compound), oct(0), octNew(0),
_U0(U0), _V0(V0), _U1(U1), _V1(V1),
_toc(toc), _allowPartition(allowPartition), _lsys(lsys)
{
@@ -1318,6 +1330,12 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
nbSplit = 0;
fillTris.clear();
+
+#if defined(HAVE_ANN)
+ index = new ANNidx[2];
+ dist = new ANNdist[2];
+#endif
+
}
GFaceCompound::~GFaceCompound()
@@ -1327,9 +1345,17 @@ GFaceCompound::~GFaceCompound()
Octree_Delete(oct);
delete [] _gfct;
}
+ if(octNew) delete(octNew);
if (_lsys)delete _lsys;
delete ONE;
delete MONE;
+#if defined(HAVE_ANN)
+ if(uv_kdtree) delete uv_kdtree;
+ if(nodes) annDeallocPts(nodes);
+ delete[]index;
+ delete[]dist;
+#endif
+
}
SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
@@ -1346,71 +1372,79 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
return SPoint2(it->second.x(),it->second.y());
}
else{
- double tGlob, tLoc;
- double tL, tR;
- int iEdge;
-
- // getParameter Point
- v->getParameter(0,tGlob);
-
- // find compound Edge
- GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
-
- if(gec){
- // compute local parameter on Edge
- gec->getLocalParameter(tGlob,iEdge,tLoc);
- std::vector<GEdge*> gev = gec->getCompounds();
- GEdge *ge = gev[iEdge];
-
- // left and right vertex of the Edge
- MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
- MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
- std::map<MVertex*,SPoint3>::iterator itL = coordinates.find(v0);
- std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
-
- // for the Edge, find the left and right vertices of the initial
- // 1D mesh and interpolate to find (u,v)
- MVertex *vL = v0;
- MVertex *vR = v1;
- double tB = ge->parBounds(0).low();
- double tE = ge->parBounds(0).high();
- int j = 0;
- tL=tB;
- bool found = false;
- while(j < (int)ge->mesh_vertices.size()){
- vR = ge->mesh_vertices[j];
- vR->getParameter(0,tR);
- if(!vR->getParameter(0,tR)) {
- Msg::Error("vertex vr %p not MedgeVertex", vR);
- Msg::Exit(1);
- }
- if(tLoc > tL && tLoc < tR){
- found = true;
- itR = coordinates.find(vR);
- if(itR == coordinates.end()){
- Msg::Error("vertex %p (%g %g %g) not found", vR, vR->x(), vR->y(), vR->z());
- Msg::Exit(1);
- }
- break;
- }
- else{
- itL = coordinates.find(vR);
- vL = vR;
- tL = tR;
- }
- j++;
- }
- if(!found) {
- vR = v1;
- tR = tE;
+
+ //if(v->onWhat()->dim() == 1){
+ double tGlob, tLoc;
+ double tL, tR;
+ int iEdge;
+
+ // getParameter Point
+ v->getParameter(0,tGlob);
+
+ //find compound Edge
+ GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
+
+ if(gec){
+ // compute local parameter on Edge
+ gec->getLocalParameter(tGlob,iEdge,tLoc);
+ std::vector<GEdge*> gev = gec->getCompounds();
+ GEdge *ge = gev[iEdge];
+
+ // left and right vertex of the Edge
+ MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
+ MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
+ std::map<MVertex*,SPoint3>::iterator itL = coordinates.find(v0);
+ std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
+
+ // for the Edge, find the left and right vertices of the initial
+ // 1D mesh and interpolate to find (u,v)
+ MVertex *vL = v0;
+ MVertex *vR = v1;
+ double tB = ge->parBounds(0).low();
+ double tE = ge->parBounds(0).high();
+ int j = 0;
+ tL=tB;
+ bool found = false;
+ while(j < (int)ge->mesh_vertices.size()){
+ vR = ge->mesh_vertices[j];
+ vR->getParameter(0,tR);
+ if(!vR->getParameter(0,tR)) {
+ Msg::Error("vertex vr %p not MedgeVertex", vR);
+ Msg::Exit(1);
+ }
+ if(tLoc > tL && tLoc < tR){
+ found = true;
+ itR = coordinates.find(vR);
+ if(itR == coordinates.end()){
+ Msg::Error("vertex %p (%g %g %g) not found", vR, vR->x(), vR->y(), vR->z());
+ Msg::Exit(1);
+ }
+ break;
+ }
+ else{
+ itL = coordinates.find(vR);
+ vL = vR;
+ tL = tR;
+ }
+ j++;
+ }
+ if(!found) {
+ vR = v1;
+ tR = tE;
+ }
+
+ // linear interpolation between tL and tR
+ double uloc, vloc;
+ uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
+ vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
+ return SPoint2(uloc,vloc);
}
+
+ // }
+ // else if(v->onWhat()->dim() == 2){
+ // Msg::Debug("Get coordinates of Compound Face cannot find vertex");
+ // }
- // linear interpolation between tL and tR
- double uloc, vloc;
- uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
- vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
- return SPoint2(uloc,vloc);
- }
}
// never here
@@ -1921,6 +1955,27 @@ double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
return fabs(t->interpolateDiv(val, u, v, 3));
}
+SPoint2 GFaceCompound::parFromVertex(MVertex *v) const
+{
+ double U=0.0,V=0.0;
+ if(v->onWhat()->dim()==2 && meshStatistics.status == GFace::DONE) {
+ v->getParameter(0, U);
+ v->getParameter(1, V);
+ }
+ if ( v->onWhat()->dim()==1 ||
+ (v->onWhat()->dim()==2 && U == -1 && V==-1)){ //if MFaceVertex created on edge in bunin
+ SPoint2 sp = getCoordinates(v);
+ U = sp.x();
+ V = sp.y();
+ }
+ if (v->onWhat()->dim()==0){
+ SPoint2 sp = parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ U = sp.x();
+ V = sp.y();
+ }
+ return SPoint2(U,V);
+
+}
SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
{
if(!oct) parametrize();
@@ -1930,6 +1985,155 @@ SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
return SPoint2(sp.x(), sp.y());
}
+GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const {
+
+ //printf("in remeshed oct for par =%g %g\n", par1,par2);
+
+ //if not meshed yet
+ if (meshStatistics.status != GFace::DONE || triangles.size()+quadrangles.size() == 0) {
+ GPoint gp (3,3,0,this);
+ gp.setNoSuccess();
+ return gp;
+ }
+
+ //create new octree with new mesh elements
+ if(!octNew){
+ //printf("****************** creating new octree \n");
+
+ //FILE * of = fopen("myOCTREE.pos","w");
+ //fprintf(of, "View \"\"{\n");
+
+ std::vector<MElement *> myElems;
+ for (int i = 0; i < triangles.size(); i++) myElems.push_back(triangles[i]);
+ for (int i = 0; i < quadrangles.size(); i++) myElems.push_back(quadrangles[i]);
+
+ std::vector<MElement *> myParamElems;
+ std::set<SPoint2> myBCNodes;
+ for (int i = 0; i < myElems.size(); i++){
+ MElement *e = myElems[i];
+ std::vector<MVertex*> myParamVert;
+ for (unsigned int j = 0; j < e->getNumVertices(); j++){
+ MVertex *v = e->getVertex(j);
+ SPoint2 sp = parFromVertex(v);
+ myParamVert.push_back(new MVertex(sp.x(),sp.y(), 0.0, NULL,v->getNum()));
+ if(v->onWhat()->dim()<2) myBCNodes.insert(sp);
+ XYZoct.insert(std::make_pair(v->getNum(),SPoint3(v->x(), v->y(), v->z())));
+ }
+ if(e->getType() == TYPE_TRI)
+ myParamElems.push_back(new MTriangle(myParamVert[0],myParamVert[1],myParamVert[2],e->getNum()));
+ else if (e->getType() == TYPE_QUA) {
+ myParamElems.push_back(new MQuadrangle(myParamVert[0],myParamVert[1],myParamVert[2],myParamVert[3],e->getNum()));
+ // fprintf(of, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ // myParamVert[0]->x(), myParamVert[0]->y(), myParamVert[0]->z(),
+ // myParamVert[1]->x(), myParamVert[1]->y(), myParamVert[1]->z(),
+ // myParamVert[2]->x(), myParamVert[2]->y(), myParamVert[2]->z(),
+ // myParamVert[3]->x(), myParamVert[3]->y(), myParamVert[3]->z(),
+ // (double)e->getNum(), (double)e->getNum(), (double)e->getNum(), (double)e->getNum());
+ }
+ }
+
+ octNew = new MElementOctree(myParamElems);
+
+ //build kdtree boundary nodes in parametric space
+#if defined(HAVE_ANN)
+ //printf("creating uv kdtree %d \n", myBCNodes.size());
+ nodes = annAllocPts(myBCNodes.size(), 3);
+ std::set<SPoint2>::iterator itp = myBCNodes.begin();
+ int ind = 0;
+ while (itp != myBCNodes.end()){
+ SPoint2 pt = *itp;
+ //fprintf(of, "SP(%g,%g,%g){%g};\n", pt.x(), pt.y(), 0.0, 10000);
+ nodes[ind][0] = pt.x();
+ nodes[ind][1] = pt.y();
+ nodes[ind][2] = 0.0;
+ itp++; ind++;
+ }
+ uv_kdtree = new ANNkd_tree(nodes, myBCNodes.size(), 3);
+#endif
+
+ //fprintf(of,"};\n");
+ //fclose(of);
+
+ }
+
+ //now find point
+ double uvw[3]={par1,par2, 0.0};
+ double UV[3];
+ double initialTol = MElement::getTolerance();
+ MElement::setTolerance(1.e-2);
+ MElement *e = octNew->find(par1,par2, 0.0,-1,false);
+ MElement::setTolerance(initialTol);
+ if (e){
+ e->xyz2uvw(uvw,UV);
+ double valX[8], valY[8], valZ[8];
+ for (int i=0;i<e->getNumPrimaryVertices();i++){
+ SPoint3 pt = XYZoct[e->getVertex(i)->getNum()];
+ valX[i] = pt.x();
+ valY[i] = pt.y();
+ valZ[i] = pt.z();
+ }
+ GPoint gp;
+ gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
+ gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
+ gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
+ //printf("found point (UV=%g %g) %g %g %g (E=%D)\n",par1, par2, gp.x(), gp.y(), gp.z(), e->getNum());
+ return gp;
+ }
+ //if element not found find closest point
+ else{
+
+ //printf("point not found look in kdtree \n");
+ GPoint gp(0,0,0);
+
+#if defined(HAVE_ANN)
+ double pt[3] = {par1,par2,0.0};
+ uv_kdtree->annkSearch(pt, 1, index, dist);
+ SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
+ SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
+ SPoint3 pnew; double d;
+ signedDistancePointLine(p1,p2, SPoint3(par1,par2,0.0), d, pnew);
+ //printf("UV=%g %g UVnew =%g %g \n", par1,par2,pnew.x(), pnew.y());
+
+ //now find point
+ double uvw[3]={pnew.x(),pnew.y(), 0.0};
+ double UV[3];
+ MElement *e = octNew->find(pnew.x(), pnew.y(), 0.0, -1);
+ if (e){
+ e->xyz2uvw(uvw,UV);
+ double valX[8], valY[8], valZ[8];
+ for (int i=0;i<e->getNumPrimaryVertices();i++){
+ SPoint3 pt = XYZoct[e->getVertex(i)->getNum()];
+ valX[i] = pt.x();
+ valY[i] = pt.y();
+ valZ[i] = pt.z();
+ }
+ GPoint gp;
+ gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
+ gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
+ gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
+ // printf("found point after kdtree (UV=%g %g) %g %g %g (E=%D)\n",par1, par2, gp.x(), gp.y(), gp.z(), e->getNum());
+ }
+ else{
+ //printf("new point kdtree not found \n");
+ gp.setNoSuccess();
+ }
+#else
+ gp.setNoSuccess();
+#endif
+
+ if (gp.succeeded()){
+ //printf("found closest point with ANN \n");
+ return gp;
+ }
+ else{
+ //printf("NOT found point %g %g \n", par1, par2);
+ GPoint gp (3,3,0,this);
+ gp.setNoSuccess();
+ return gp;
+ }
+ }
+
+}
GPoint GFaceCompound::point(double par1, double par2) const
{
@@ -1942,11 +2146,10 @@ GPoint GFaceCompound::point(double par1, double par2) const
double U,V;
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
- getTriangle (par1, par2, <, U,V);
- SPoint3 p(3, 3, 0);
+ getTriangle(par1, par2, <, U,V);
if(!lt && _mapping != RBF){
- GPoint gp (p.x(),p.y(),p.z(),this);
- gp.setNoSuccess();
+ //printf("POINT no success %d tris %d quad \n", triangles.size(), quadrangles.size());
+ GPoint gp = pointInRemeshedOctree(par1,par2);
return gp;
}
else if (!lt && _mapping == RBF){
@@ -1970,7 +2173,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
const bool LINEARMESH = true; //false
if(LINEARMESH){
// linear Lagrange mesh
- p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
+ SPoint3 p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
return GPoint(p.x(),p.y(),p.z(),this,par);
}
else{
@@ -2181,12 +2384,19 @@ void GFaceCompound::buildOct() const
count++;
}
}
+
+ // make bounding box larger up to (absolute) geometrical tolerance
+ double eps = 0.0; //CTX::instance()->geom.tolerance;
+ SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
+ bbmin -= bbeps;
+ bbmax += bbeps;
+ double origin[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
+ double ssize[3] = {bbmax.x() - bbmin.x(),
+ bbmax.y() - bbmin.y(),
+ bbmax.z() - bbmin.z()};
+
_gfct = new GFaceCompoundTriangle[count];
- double origin[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
- double ssize[3] = {bb.max().x() - bb.min().x(),
- bb.max().y() - bb.min().y(),
- bb.max().z() - bb.min().z()};
- const int maxElePerBucket = 11;
+ const int maxElePerBucket = 15;
oct = Octree_Create(maxElePerBucket, origin, ssize, GFaceCompoundBB,
GFaceCompoundCentroid, GFaceCompoundInEle);
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index d9fffb7fa7cd0e71242fc91411c4b566ea480217..3d146edba9409a3ab8841c24b9e701e82050f282 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -20,6 +20,12 @@ template <class scalar> class simpleFunction;
#include "meshGFaceOptimize.h"
#include "linearSystem.h"
#include "GRbf.h"
+#include "MElementOctree.h"
+
+#if defined(HAVE_ANN)
+#include <ANN/ANN.h>
+class ANNkd_tree;
+#endif
#define AR_MAX 5 //maximal geometrical aspect ratio
@@ -70,6 +76,8 @@ class GFaceCompound : public GFace {
mutable int nbT;
mutable GFaceCompoundTriangle *_gfct;
mutable Octree *oct;
+ mutable MElementOctree *octNew;
+ mutable std::map<int,SPoint3> XYZoct;
mutable std::set<MVertex*> allNodes;
mutable v2t_cont adjv;
mutable bool mapv2Tri;
@@ -83,6 +91,12 @@ class GFaceCompound : public GFace {
mutable std::vector<double> _coords;
mutable std::map<MVertex*, int> _mapV;
linearSystem <double> *_lsys;
+#if defined(HAVE_ANN)
+ mutable ANNkd_tree *uv_kdtree;
+ mutable ANNpointArray nodes;
+ ANNidxArray index;
+ ANNdistArray dist;
+#endif
void buildOct() const ;
void buildAllNodes() const;
@@ -133,8 +147,10 @@ class GFaceCompound : public GFace {
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
- virtual GPoint point(double par1, double par2) const;
+ GPoint point(double par1, double par2) const;
+ GPoint pointInRemeshedOctree(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
+ SPoint2 parFromVertex(MVertex *v) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index c7e6a94ce84b69024f70665c2cb9cab2acbe35ab..578b4a15c890fc686e215e2cca62e44177d3c9b0 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -89,7 +89,7 @@ MElementOctree::MElementOctree(GModel *m) : _gm(m)
Octree_Arrange(_octree);
}
-MElementOctree::MElementOctree(std::vector<MElement*> &v) : _gm(0)
+MElementOctree::MElementOctree(std::vector<MElement*> &v) : _gm(0), _elems(v)
{
SBoundingBox3d bb;
for (unsigned int i = 0; i < v.size(); i++){
@@ -156,7 +156,7 @@ MElement *MElementOctree::find(double x, double y, double z, int dim, bool stric
if (!strict && _gm) {
double initialTol = MElement::getTolerance();
double tol = initialTol;
- while (tol < 1){
+ while (tol < 1.){
tol *= 10;
MElement::setTolerance(tol);
std::vector<GEntity*> entities;
@@ -174,7 +174,26 @@ MElement *MElementOctree::find(double x, double y, double z, int dim, bool stric
}
}
MElement::setTolerance(initialTol);
- Msg::Warning("Point %g %g %g not found",x,y,z);
+ //Msg::Warning("Point %g %g %g not found",x,y,z);
+ }
+ else if (!strict && !_gm){
+ double initialTol = MElement::getTolerance();
+ double tol = initialTol;
+ while (tol < 0.1){
+ tol *= 10.0;
+ MElement::setTolerance(tol);
+ for(unsigned int i = 0; i < _elems.size(); i++){
+ e = _elems[i];
+ if (dim == -1 || e->getDim() == dim){
+ if (MElementInEle(e, P)){
+ MElement::setTolerance(initialTol);
+ return e;
+ }
+ }
+ }
+ }
+ MElement::setTolerance(initialTol);
+ //Msg::Warning("Point %g %g %g not found",x,y,z);
}
return NULL;
}
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 3b46e887f2d6e689c4d72d5a2336c43be006050e..86d226fa98efb7a88604af7615ef83c7e7dbdd6e 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -16,6 +16,7 @@ class MElementOctree{
private:
Octree *_octree;
GModel *_gm;
+ std::vector<MElement*> _elems;
public:
MElementOctree(GModel *);
MElementOctree(std::vector<MElement*> &);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 5ff5db93e80cf1885480eb45dbba48f14ea6ab8d..1b61055ccc06df9b9ba502f32604766b5bd26acd 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -263,10 +263,9 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms
{
params.clear();
- if (gf->geomType() == GEntity::CompoundSurface &&
- v->onWhat()->dim() < 2){
+ if (gf->geomType() == GEntity::CompoundSurface ){
GFaceCompound *gfc = (GFaceCompound*) gf;
- params.push_back(gfc->getCoordinates(v));
+ params.push_back(gfc->parFromVertex(v));
return;
}
@@ -364,16 +363,23 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
-bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m,
+bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 ¶m,
bool onSurface)
{
- if (gf->geomType() == GEntity::CompoundSurface &&
- v->onWhat()->dim() < 2){
+
+ if (gf->geomType() == GEntity::CompoundSurface ){
GFaceCompound *gfc = (GFaceCompound*) gf;
- param = gfc->getCoordinates(const_cast<MVertex*>(v));
+ param = gfc->parFromVertex(v);
return true;
}
+ // if (gf->geomType() == GEntity::CompoundSurface &&
+ // v->onWhat()->dim() < 2){
+ // GFaceCompound *gfc = (GFaceCompound*) gf;
+ // param = gfc->getCoordinates(const_cast<MVertex*>(v));
+ // return true;
+ // }
+
if(v->onWhat()->geomType() == GEntity::DiscreteCurve ||
v->onWhat()->geomType() == GEntity::BoundaryLayerCurve){
param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()), onSurface);
@@ -416,7 +422,7 @@ bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m,
return true;
}
-bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double ¶m)
+bool reparamMeshVertexOnEdge( MVertex *v, const GEdge *ge, double ¶m)
{
param = 1.e6;
Range<double> bounds = ge->parBounds(0);
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index b34676b630620495c238123f51c0a3bf5a9fe3e3..6013a87bf410620bd35c165aff7b4bf4098a73c4 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -164,9 +164,9 @@ class MFaceVertex : public MVertex{
bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
SPoint2 ¶m1, SPoint2 ¶m2);
-bool reparamMeshVertexOnFace(const MVertex *v, const GFace *gf, SPoint2 ¶m,
+bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 ¶m,
bool onSurface=true);
-bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double ¶m);
+bool reparamMeshVertexOnEdge(MVertex *v, const GEdge *ge, double ¶m);
double angle3Vertices(MVertex *p1, MVertex *p2, MVertex *p3);
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 60fd0e49e2a4c8b48414fb08122ab6f4dec3156c..dd19a7459e5faeabb8dd710abb046d5cecf61476 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -259,28 +259,24 @@ void meshGFaceBamg(GFace *gf){
//If point not found because compound edges have been remeshed and boundary triangles have changed
//then we call our new octree
if ( !gp.succeeded() && hasCompounds){
- SPoint3 uvw(v[0],v[1], 0.0), UV;
+ double uvw[3] = {v[0],v[1], 0.0};
+ double UV[3];
double initialTol = MElement::getTolerance();
MElement::setTolerance(1.e-2);
MElement *e = _octree->find(v[0],v[1], 0.0, -1);
MElement::setTolerance(initialTol);
if (e){
e->xyz2uvw(uvw,UV);
- double *valX = new double[e->getNumVertices()];
- double *valY = new double[e->getNumVertices()];
- double *valZ = new double[e->getNumVertices()];
- for (int i=0;i<e->getNumVertices();i++){
+ double valX[8], valY[8], valZ[8];
+ for (int i=0;i<e->getNumPrimaryVertices();i++){
int numTri = e->getNum();
- valX[i] = gf->triangles[numTri]->getVertex(0)->x();
- valY[i] = gf->triangles[numTri]->getVertex(0)->y();
- valZ[i] = gf->triangles[numTri]->getVertex(0)->z();
+ valX[i] = gf->triangles[numTri]->getVertex(i)->x();
+ valY[i] = gf->triangles[numTri]->getVertex(i)->y();
+ valZ[i] = gf->triangles[numTri]->getVertex(i)->z();
}
gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
- delete [] valX;
- delete [] valY;
- delete [] valZ;
}
}
MFaceVertex *x = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, v[0], v[1]);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 1679b14c1468e53394fcb16fe472065921190705..f0a99fb19ece8a4b67d7ef4965415b154416cbb5 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -620,6 +620,7 @@ static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
double qualityNew[256];
GPoint gp = gf->point(after);
+ if (!gp.succeeded())return false;
SPoint3 pafter (gp.x(),gp.y(),gp.z());
SPoint3 pbefore (v1->x(),v1->y(),v1->z());
@@ -772,13 +773,19 @@ static std::vector<MVertex*> closestPathBetweenTwoDefects (v2t_cont &adj,
static MVertex * createNewVertex (GFace *gf, SPoint2 p){
GPoint gp = gf->point(p);
+ if (!gp.succeeded()) {
+ //printf("**** ARRG new vertex not created \n");
+ return NULL;
+ }
return new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
}
static std::vector<MVertex*> saturateEdge (GFace *gf, SPoint2 p1, SPoint2 p2, int n){
std::vector<MVertex*> pts;
for (int i=1;i<n;i++){
SPoint2 p = p1 + (p2-p1)*((double)i/((double)(n)));
- pts.push_back(createNewVertex(gf,p));
+ MVertex *v = createNewVertex(gf,p);
+ if (!v){ pts.clear(); pts.resize(0); return pts;}
+ pts.push_back(v);
}
return pts;
}
@@ -812,7 +819,7 @@ void createRegularMesh (GFace *gf,
std::vector<MQuadrangle*> &q){
int N = e12.size();
int M = e23.size();
-
+
//create points using transfinite interpolation
std::vector<std::vector<MVertex*> > tab(M+2);
@@ -828,7 +835,7 @@ void createRegularMesh (GFace *gf,
}
for (int i=0;i<M;i++){
tab[i+1][N+1] = sign23 > 0 ? e23 [i] : e23 [M-i-1];
- tab[i+1][0] = sign41 < 0 ? e41 [i] : e41 [M-i-1];
+ tab[i+1][0] = sign41 < 0 ? e41 [i] : e41 [M-i-1];
}
for (int i=0;i<N;i++){
@@ -847,6 +854,14 @@ void createRegularMesh (GFace *gf,
c1.x(),c2.x(),c3.x(),c4.x(),u,v);
double Vp = TRAN_QUA(p12.y(), p23.y(), p34.y(), p41.y(),
c1.y(),c2.y(),c3.y(),c4.y(),u,v);
+ // printf("v1=%d v2=%d v3=%d v4=%d \n", v1->getNum(), v2->getNum(), v3->getNum(), v4->getNum());
+ // printf("c1=%g %g, c2=%g %g, c3=%g %g, c4=%g,%g \n", c1.x(),c1.y(),c2.x(),c2.y(),c3.x(),c3.y(),c4.x(),c4.y());
+ // printf("p1=%g %g, p2=%g %g, p3=%g %g, p4=%g,%g \n", p12.x(),p12.x(),p23.x(),p23.y(),p34.x(),p34.y(),p41.x(),p41.y());
+ if (p12.x() && p12.y() == -1.0) {printf("wrong param -1 \n"); exit(1);}
+ if (p23.x() && p23.y() == -1.0) {printf("wrong param -1 \n"); exit(1);}
+ if (p34.x() && p34.y() == -1.0) {printf("wrong param -1 \n"); exit(1);}
+ if (p41.x() && p41.y() == -1.0) {printf("wrong param -1 \n"); exit(1);}
+
MVertex *vnew = createNewVertex (gf, SPoint2(Up,Vp));
tab[j+1][i+1] = vnew;
}
@@ -854,7 +869,6 @@ void createRegularMesh (GFace *gf,
// create quads
for (int i=0;i<=N;i++){
for (int j=0;j<=M;j++){
- //MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j+1][i],tab[j+1][i+1],tab[j][i+1]);
MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
q.push_back(qnew);
}
@@ -1163,7 +1177,6 @@ struct quadBlob {
else count[side]++;
}
-
if (ncorners == 3){
fullVector<double> rhs(6);
@@ -1188,6 +1201,9 @@ struct quadBlob {
std::vector<MVertex*> e012_01 = saturateEdge (gf,p012,p01,a2);
std::vector<MVertex*> e012_12 = saturateEdge (gf,p012,p12,a3);
std::vector<MVertex*> e012_20 = saturateEdge (gf,p012,p20,a1);
+ if (e012_01.size() == 0) return false;
+ if (e012_12.size() == 0) return false;
+ if (e012_20.size() == 0) return false;
std::vector<MVertex*> e0_01,e01_1,e1_12,e12_2,e2_20,e20_0;
for (int i=0;i<a1-1;i++)e0_01.push_back(bnodes[i+1]);
@@ -1337,6 +1353,11 @@ struct quadBlob {
std::vector<MVertex*> e01234_23 = saturateEdge (gf,p01234,p23,a4);
std::vector<MVertex*> e01234_34 = saturateEdge (gf,p01234,p34,a5);
std::vector<MVertex*> e01234_40 = saturateEdge (gf,p01234,p40,a1);
+ if (e01234_01.size() == 0) return false;
+ if (e01234_12.size() == 0) return false;
+ if (e01234_23.size() == 0) return false;
+ if (e01234_34.size() == 0) return false;
+ if (e01234_40.size() == 0) return false;
std::vector<MVertex*> e0_01,e01_1,e1_12,e12_2,e2_23,e23_3,e3_34,e34_4,e4_40,e40_0;
for (int i=0;i<a1-1;i++)e0_01.push_back(bnodes[i+1]);
@@ -1519,6 +1540,10 @@ static int _splitFlatQuads(GFace *gf, double minQual)
GPoint gb2 = gf->point(pb2);
GPoint gb3 = gf->point(pb3);
GPoint gb4 = gf->point(pb4);
+ if (!gb1.succeeded())return -1;
+ if (!gb2.succeeded())return -1;
+ if (!gb3.succeeded())return -1;
+ if (!gb4.succeeded())return -1;
MFaceVertex *b1 = new MFaceVertex (gb1.x(),gb1.y(),gb1.z(),gf,gb1.u(), gb1.v());
MFaceVertex *b2 = new MFaceVertex (gb2.x(),gb2.y(),gb2.z(),gf,gb2.u(), gb2.v());
MFaceVertex *b3 = new MFaceVertex (gb3.x(),gb3.y(),gb3.z(),gf,gb3.u(), gb3.v());
@@ -1748,6 +1773,7 @@ struct opti_data_vertex_relocation {
}
void set_(double U, double V){
GPoint gp = gf->point(U,V);
+ if (!gp.succeeded()) Msg::Error("point not OK \n");
v->x() = gp.x();
v->y() = gp.y();
v->z() = gp.z();
@@ -1807,7 +1833,7 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
// if (rep.terminationtype != 4){
// data.set_(U,V);
// }
- printf("heyhey -- end\n");
+
}
#endif
@@ -1873,7 +1899,7 @@ static void _relocateVertex(GFace *gf, MVertex *ver,
}
else{
#if defined(HAVE_BFGS)
- _relocateVertexOpti(gf, ver, lt);
+ _relocateVertexOpti(gf, ver, lt);
#endif
}
}
@@ -1938,15 +1964,15 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
ed = e2->getEdge(i);
if (ed.getVertex(0) == v1 && ed.getVertex(1) != v2) v21 = ed.getVertex(1);
if (ed.getVertex(1) == v1 && ed.getVertex(0) != v2) {rot2 = -1.0; v21 = ed.getVertex(0);}
- if (ed.getVertex(0) == v2 && ed.getVertex(1) != v1) {rot2 = -1.0;v22 = ed.getVertex(1);}
+ if (ed.getVertex(0) == v2 && ed.getVertex(1) != v1) {rot2 = -1.0; v22 = ed.getVertex(1);}
if (ed.getVertex(1) == v2 && ed.getVertex(0) != v1) v22 = ed.getVertex(0);
}
if (!v11 || !v12 || !v21 || !v22){
- Msg::Error("You found a BUG in the quad optimization routines of Gmsh Line %d of FILE %s",__LINE__,__FILE__);
+ Msg::Error("Quad vertices are wrong (in edgeSwapQuads opti)");
return 0;
}
if (rot1 != rot2){
- Msg::Error("Quads not oriented the same (in edgeSwapQuads opti)");
+ Msg::Error("Quads are not oriented the same (in edgeSwapQuads opti)");
//exit(1);
return 0;
}