diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 68295e6da7e23aab51aec931325452a8c26c34bc..032e814b5c0e0890ac837711890a895406ae416d 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1392,8 +1392,9 @@ static void meshCompound (GFace* gf, bool verbose) {
if (found == ec.end())ec.insert((*it)->tag());
else ec.erase(found);
}
- df->triangles.insert(df->triangles.begin(), c->triangles.begin(),c->triangles.end());
- df->mesh_vertices.insert(df->mesh_vertices.begin(), c->mesh_vertices.begin(),c->mesh_vertices.end());
+ df->triangles.insert(df->triangles.end(), c->triangles.begin(),c->triangles.end());
+ df->mesh_vertices.insert(df->mesh_vertices.end(), c->mesh_vertices.begin(),c->mesh_vertices.end());
+ for (unsigned int j=0;j<c->triangles.size();j++)df->_CAD.push_back(c);
c->triangles.clear();
c->mesh_vertices.clear();
}
@@ -1402,15 +1403,6 @@ static void meshCompound (GFace* gf, bool verbose) {
df->setBoundEdges(gf->model(), cedges);
df->createGeometry();
df->mesh(verbose);
- gf->triangles = df->triangles;
- for (unsigned int i=0;i<df->mesh_vertices.size();i++){
- MVertex *v = df->mesh_vertices[i];
- v->setEntity(gf);
- gf->mesh_vertices.push_back(v);
- }
- df->triangles.clear();
- df->mesh_vertices.clear();
-
delete df;
}
#endif
@@ -1489,68 +1481,6 @@ void GFace::moveToValidRange(SPoint2 &pt) const
}
}
-void GFace::addLayersOfQuads(int nLayers, GVertex *gv, double hmin, double ratio)
-{
- SVector3 ez (0, 0, 1);
- std::list<GEdgeLoop>::iterator it = edgeLoops.begin();
- FILE *f = Fopen ("coucou.pos","w");
- if(f) fprintf(f,"View \"\"{\n");
- for (; it != edgeLoops.end(); ++it){
- bool found = false;
- // look if this edge loop has the GVertex as an endpoint
- for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
- if (it2->ge->getBeginVertex() == gv || it2->ge->getEndVertex() == gv)
- found = true;
- }
- // we found an edge loop with the GVertex that was specified
- if (found){
- // first build a list of edges in the parametric space
- std::vector<std::pair<MVertex*,SPoint2> > contour;
- for (GEdgeLoop::iter it2 = it->begin(); it2 != it->end(); ++it2){
- GEdge *ge = it2->ge;
- SPoint2 p[2];
- if (it2->_sign == 1){
- for (unsigned int i = 0; i < ge->lines.size(); i++){
- reparamMeshEdgeOnFace(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),
- this, p[0], p[1]);
- contour.push_back(std::make_pair(ge->lines[i]->getVertex(0), p[0]));
- }
- }
- else {
- for(int i = ge->lines.size() - 1; i >= 0; i--){
- reparamMeshEdgeOnFace(ge->lines[i]->getVertex(0), ge->lines[i]->getVertex(1),
- this,p[0],p[1]);
- contour.push_back(std::make_pair(ge->lines[i]->getVertex(1),p[1]));
- }
- }
- }
- double hlayer = hmin;
- for (int j = 0; j < nLayers; j++){
- for (unsigned int i = 0; i < contour.size(); i++){
- SPoint2 p0 = contour[(i+0) % contour.size()].second;
- SPoint2 p1 = contour[(i+1) % contour.size()].second;
- SPoint2 p2 = contour[(i+2) % contour.size()].second;
- SVector3 p0p1 (p1.x()-p0.x(),p1.y()-p0.y(),0.0);
- SVector3 p1p2 (p2.x()-p1.x(),p2.y()-p1.y(),0.0);
- SVector3 n01 = crossprod(ez,p0p1);
- SVector3 n12 = crossprod(ez,p1p2);
- SVector3 n = (n01+n12)*-0.5;
- n.normalize();
- double u = p1.x() + n.x() * hmin;
- double v = p1.y() + n.y() * hmin;
- GPoint gp = point(SPoint2(u,v));
- additionalVertices.push_back(new MFaceVertex(gp.x(),gp.y(),gp.z(),this,u,v));
- if(f) fprintf(f,"SP(%g, %g, 0){1};\n",gp.x(),gp.y());
- }
- hlayer *= ratio;
- hmin += hlayer;
- }
- if(f) fprintf(f,"};\n");
- }
- }
- if(f) fclose(f);
-}
-
void GFace::relocateMeshVertices()
{
for(unsigned int i = 0; i < mesh_vertices.size(); i++){
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 0e06f8c8f4ea10b4084ea88f1366f0725439a0e2..2a3c020442e6d5719725a7aa13bda6207e379f35 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -289,9 +289,6 @@ class GFace : public GEntity{
// new interface for meshing
virtual void mesh(bool verbose);
- // add layers of quads
- void addLayersOfQuads(int nLayers, GVertex *start, double hmin, double factor);
-
struct {
// do we recombine the triangles of the mesh?
int recombine;
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 11693baf264644ded75773dbc0672d3901628ada..7262b86987e1004b5000968d2ba58e6ffc13e179 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -220,7 +220,7 @@ static bool orderVertices(const double &tot_length, const std::vector<MVertex*>
}
/*BUILDER*/
-discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int p) :
+discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int p, std::vector<GFace*> *CAD) :
GFace(gf->model(),123), _parent (gf),_ddft(NULL), oct(NULL)
{
_order = p;
@@ -237,7 +237,14 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int
if (v->onWhat()->dim() == 2) {
std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v);
if (it == v2v.end()){
- MFaceVertex *vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), 0, 0);
+ MFaceVertex *vv;
+ if (!CAD) vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), 0, 0);
+ else{
+ GFace *cad = (*CAD)[i];
+ if(cad != v->onWhat())Msg::Fatal("Line %d FILE %s : erroneous cad list",__LINE__,__FILE__);
+ double pu,pv; v->getParameter(0,pu);v->getParameter(1,pv);
+ vv = new MFaceVertex ( v->x(), v->y(), v->z(), v->onWhat(), pu, pv);
+ }
v2v[v] = vv;
discrete_vertices.push_back(vv);
vs.push_back(vv);
@@ -280,11 +287,11 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh, int
orderVertices(_totLength, _U0, _coords);
parametrize(false);
- buildOct();
+ buildOct(CAD);
if (!checkOrientationUV()){
parametrize(true);
- buildOct();
+ buildOct(CAD);
}
putOnView();
@@ -453,7 +460,7 @@ void discreteDiskFace::getBoundingEdges()
}
-void discreteDiskFace::buildOct() const
+void discreteDiskFace::buildOct(std::vector<GFace*> *CAD) const
{
if (oct)Octree_Delete(oct);
@@ -464,6 +471,9 @@ void discreteDiskFace::buildOct() const
discreteDiskFaceCentroid, discreteDiskFaceInEle);
_ddft = new discreteDiskFaceTriangle[discrete_triangles.size()];
+
+ // if (CAD) printf("-------------> %ld %ld\n",CAD->size(),discrete_triangles.size());
+
for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
MTriangle *t = discrete_triangles[i];
discreteDiskFaceTriangle* my_ddft = &_ddft[i];
@@ -471,7 +481,7 @@ void discreteDiskFace::buildOct() const
for(int io=0; io<_N; io++)
my_ddft->p[io] = coordinates[t->getVertex(io)];
- my_ddft->gf = _parent;
+ my_ddft->gf = CAD ? (*CAD)[i] : _parent;
my_ddft->tri = t;
Octree_Insert(my_ddft, oct);
@@ -589,6 +599,7 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,
_eta = Xi[1];
}
+
bool discreteDiskFace::checkOrientationUV(){
double initial, current; // initial and current orientation
@@ -814,11 +825,11 @@ void discreteDiskFace::putOnView()
for (int j=0; j<_N-1; j++){
fprintf(view_u,"%g,",my_ddft->p[j].x());
fprintf(view_v,"%g,",my_ddft->p[j].y());
- fprintf(UVxyz,"%g,",sqrt(my_ddft->tri->getVertex(j)->x()*my_ddft->tri->getVertex(j)->x()+my_ddft->tri->getVertex(j)->z()*my_ddft->tri->getVertex(j)->z()+my_ddft->tri->getVertex(j)->y()*my_ddft->tri->getVertex(j)->y()));
+ fprintf(UVxyz,"%d,",my_ddft->gf->tag());
}
fprintf(view_u,"%g};\n",my_ddft->p[_N-1].x());
fprintf(view_v,"%g};\n",my_ddft->p[_N-1].y());
- fprintf(UVxyz,"%g};\n",sqrt(my_ddft->tri->getVertex(_N-1)->x()*my_ddft->tri->getVertex(_N-1)->x()+my_ddft->tri->getVertex(_N-1)->z()*my_ddft->tri->getVertex(_N-1)->z()+my_ddft->tri->getVertex(_N-1)->y()*my_ddft->tri->getVertex(_N-1)->y()));
+ fprintf(UVxyz,"%d};\n",my_ddft->gf->tag());
}
fprintf(view_u,"};\n");
fprintf(view_v,"};\n");
@@ -827,28 +838,6 @@ void discreteDiskFace::putOnView()
fclose(view_v);
fclose(UVxyz);
}
-
- /*
-#ifdef HAVE_POST
- std::map<int, std::vector<double> > u;
- std::map<int, std::vector<double> > v;
- for(std::set<MVertex*>::iterator iv = allNodes.begin(); iv!=allNodes.end(); ++iv){
- MVertex *p = *iv;
- u[p->getNum()].push_back(coordinates[p].x());
- v[p->getNum()].push_back(coordinates[p].y());
- }
- PView* view_u = new PView("u", "NodeData", GModel::current(), u);
- PView* view_v = new PView("v", "NodeData", GModel::current(), v);
- view_u->setChanged(true);
- view_v->setChanged(true);
- view_u->write("param_u.msh", 5);
- view_v->write("param_v.msh", 5);
- delete view_u;
- delete view_v;
-#else
- Msg::Error("discreteDiskFace: cannot export without post module")
-#endif
- */
}
// useful for mesh generators ----------------------------------------
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index c4ea3f4366d4c7f2fd50b5a8b9d26583329252e7..1fb99a2df7302f82179af44bdf502a417c7283ea 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -41,7 +41,7 @@ class discreteDiskFaceTriangle {
class discreteDiskFace : public GFace {
GFace *_parent;
- void buildOct() const;
+ void buildOct(std::vector<GFace*> *CAD = NULL) const;
bool parametrize(bool one2one = false) const;
void buildAllNodes() ;
void getBoundingEdges();
@@ -49,7 +49,7 @@ class discreteDiskFace : public GFace {
bool checkOrientationUV();
public:
- discreteDiskFace(GFace *parent, std::vector<MTriangle*> &mesh, int p=1);// MTriangle -> MTriangle 6
+ discreteDiskFace(GFace *parent, std::vector<MTriangle*> &mesh, int p=1, std::vector<GFace*> *CAD = NULL);// MTriangle -> MTriangle 6
virtual ~discreteDiskFace();
void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v) const;
GPoint point(double par1, double par2) const;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 67f170dae42f5094fb1cd27c391f8ebc7222392e..039681c756f427f2d29ccd680daf05f01a3050e8 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -109,7 +109,8 @@ void discreteFace::createGeometry()
#if defined(HAVE_ANN) && defined(HAVE_SOLVER)
if (!_atlas.empty())return;
// parametrization is done here !!!
- discreteDiskFace *df = new discreteDiskFace (this, triangles,2);
+ // if (_CAD != empty) --> _triangles[i] is classified on _CAD[i]
+ discreteDiskFace *df = new discreteDiskFace (this, triangles,1,(_CAD.empty() ? NULL : &_CAD));
df->replaceEdges(l_edges);
_atlas.push_back(df);
#endif
@@ -123,10 +124,35 @@ void discreteFace::gatherMeshes()
triangles.clear();
mesh_vertices.clear();
for (unsigned int i=0;i<_atlas.size();i++){
- triangles.insert(triangles.begin(), _atlas[i]->triangles.begin(),
- _atlas[i]->triangles.end());
- mesh_vertices.insert(mesh_vertices.begin(), _atlas[i]->mesh_vertices.begin(),
- _atlas[i]->mesh_vertices.end());
+ for (unsigned int j=0;j<_atlas[i]->triangles.size(); j++){
+ MTriangle *t = _atlas[i]->triangles[j];
+ SPoint2 p0,p1,p2;
+ reparamMeshVertexOnFace(t->getVertex(0),_atlas[i], p0);
+ reparamMeshVertexOnFace(t->getVertex(1),_atlas[i], p1);
+ reparamMeshVertexOnFace(t->getVertex(2),_atlas[i], p2);
+ SPoint2 pc = (p0+p1+p2)*(1./3.0);
+ discreteDiskFaceTriangle *mt=NULL;
+ double xi, eta;
+ _atlas[i]->getTriangleUV(pc.x(),pc.y(), &mt, xi,eta);
+ if (mt && mt->gf)mt->gf->triangles.push_back(t);
+ else Msg::Warning ("FILE %s LINE %d Triangle has no classification",__FILE__,__LINE__);
+ }
+ _atlas[i]->triangles.clear();
+
+ for (unsigned int j=0;j<_atlas[i]->mesh_vertices.size(); j++){
+ MVertex *v = _atlas[i]->mesh_vertices[j];
+ double pu,pv; v->getParameter(0,pu);v->getParameter(1,pv);
+ discreteDiskFaceTriangle *mt;
+ double xi, eta;
+ _atlas[i]->getTriangleUV(pu,pv, &mt, xi,eta);
+ if(mt && mt->gf){
+ v->setEntity(mt->gf);
+ // here we should recompute on the CAD if necessary
+ mt->gf->mesh_vertices.push_back(v);
+ }
+ else Msg::Warning ("FILE %s LINE %d Vertex has no classification",__FILE__,__LINE__);
+ }
+ _atlas[i]->mesh_vertices.clear();
}
#endif
}
@@ -135,8 +161,9 @@ void discreteFace::mesh(bool verbose)
{
#if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
if (!CTX::instance()->meshDiscrete) return;
- for (unsigned int i=0;i<_atlas.size();i++)
+ for (unsigned int i=0;i<_atlas.size();i++){
_atlas[i]->mesh(verbose);
+ }
gatherMeshes();
meshStatistics.status = GFace::DONE;
#endif
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 4faf8673758eb32e782d78e75f5e263aad998566..c489a934cf38cb98e51f3694007371b364b9de6d 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -39,6 +39,7 @@ class discreteFace : public GFace {
void gatherMeshes();
virtual void mesh (bool verbose);
std::vector<discreteDiskFace*> _atlas;
+ std::vector<GFace*> _CAD;
};
#endif
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 91088d1cff6e9b7f2759c43295d752a04eca084b..60dad18307d15d6d444e090f821bf3f4663270d6 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -817,16 +817,35 @@ void _relocateVertex(GFace *gf, MVertex *ver,
}
}
+void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs){
+ vs.clear();
+ BoundaryLayerColumns* _columns = gf->getColumns();
+ if (!_columns)return;
+ for ( std::map<MElement*,std::vector<MElement*> >::iterator it = _columns->_elemColumns.begin();
+ it != _columns->_elemColumns.end();it++){
+ std::vector<MElement *> e = it->second;
+ for (unsigned int i=0;i<e.size();i++){
+ for (int j=0;j<e[i]->getNumVertices();j++){
+ vs.insert(e[i]->getVertex(j));
+ }
+ }
+ }
+}
+
void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
{
if (!niter)return;
+ std::set<MVertex*> vs;
+ getAllBoundaryLayerVertices (gf, vs);
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
for(int i = 0; i < niter; i++){
v2t_cont::iterator it = adj.begin();
while (it != adj.end()){
- _relocateVertex(gf, it->first, it->second);
+ if (vs.find(it->first) == vs.end()){
+ _relocateVertex(gf, it->first, it->second);
+ }
++it;
}
}
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index 32522126bbd8537cd5cb8c0a361e360c42a463ce..ac2c899880d55e499d8e9aea9b34c52b201ba9ab 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -217,14 +217,21 @@ static double _relocateVertex(GFace* gf, MVertex *ver,
}
+void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs);
+
void RelocateVertices (GFace* gf, int niter, double tol) {
+ std::set<MVertex*> vs;
+ getAllBoundaryLayerVertices (gf, vs);
+
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
for (int i=0;i<niter;i++){
v2t_cont::iterator it = adj.begin();
while (it != adj.end()){
- _relocateVertex( gf, it->first, it->second, tol);
+ if (vs.find(it->first) == vs.end()){
+ _relocateVertex( gf, it->first, it->second, tol);
+ }
++it;
}
}
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index cc4192043dfcb59ca3f2dbe1e03e81f07fd39965..65bc9c890752094d2480312ed83bde18809734d2 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -22,7 +22,6 @@
#include "MVertex.h"
#include "MElement.h"
-//#include "directions3D.h"
#include "BackgroundMesh.h"
#include "intersectCurveSurface.h"
@@ -30,173 +29,6 @@
using namespace std;
-//static const double FACTOR = .71;
-//static const int NUMDIR = 3;
-//static const double DIRS [NUMDIR] = {0.0, M_PI/20.,-M_PI/20.};
-//PE MODIF
-//static const int NUMDIR = 1;
-//static const double DIRS [NUMDIR] = {0.0};
-// END PE MODIF
-
-
-
-/// a rectangle in the tangent plane is transformed
-/// into a parallelogram. We define an exclusion zone
-/// that is centered around a vertex and that is used
-/// in a r-tree structure for generating points with the
-/// right spacing in the tangent plane
-
-
-//struct surfacePointWithExclusionRegion {
-// MVertex *_v;
-// SPoint2 _center;
-// SPoint2 _p[4][NUMDIR];
-// SPoint2 _q[4];
-// SMetric3 _meshMetric;
-// double _distanceSummed;
-// /*
-// + p3
-// p4 |
-// +----c-----+ p2
-// |
-// + p1
-//
-//*/
-//
-// surfacePointWithExclusionRegion (MVertex *v, SPoint2 p[4][NUMDIR], SPoint2 &_mp, SMetric3 & meshMetric, surfacePointWithExclusionRegion *father = 0){
-// _v = v;
-// _meshMetric = meshMetric;
-// _center = _mp;
-// for (int i=0;i<4;i++)_q[i] = _center + (p[i][0]+p[(i+1)%4][0]-_center*2)*FACTOR;
-// for (int i=0;i<4;i++)for (int j=0;j<NUMDIR;j++)_p[i][j] = p[i][j];
-//
-// if (!father){
-// fullMatrix<double> V(3,3);
-// fullVector<double> S(3);
-// meshMetric.eig(V,S);
-// double l = std::max(std::max(S(0),S(1)),S(2));
-// _distanceSummed = sqrt(1/(l*l));
-// }
-// else {
-// _distanceSummed = father->_distanceSummed + distance (father->_v,_v);
-// }
-// }
-// bool inExclusionZone (const SPoint2 &p){
-// double mat[2][2];
-// double b[2] , uv[2];
-// mat[0][0]= _q[1].x()-_q[0].x();
-// mat[0][1]= _q[2].x()-_q[0].x();
-// mat[1][0]= _q[1].y()-_q[0].y();
-// mat[1][1]= _q[2].y()-_q[0].y();
-// b[0] = p.x() - _q[0].x();
-// b[1] = p.y() - _q[0].y();
-// sys2x2(mat, b, uv);
-// // printf("inversion 1 : %g %g \n",uv[0],uv[1]);
-// if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
-// mat[0][0]= _q[3].x()-_q[2].x();
-// mat[0][1]= _q[0].x()-_q[2].x();
-// mat[1][0]= _q[3].y()-_q[2].y();
-// mat[1][1]= _q[0].y()-_q[2].y();
-// b[0] = p.x() - _q[2].x();
-// b[1] = p.y() - _q[2].y();
-// sys2x2(mat, b, uv);
-// // printf("inversion 2 : %g %g \n",uv[0],uv[1]);
-// if (uv[0] >= 0 && uv[1] >= 0 && 1.-uv[0] - uv[1] >= 0)return true;
-// return false;
-// }
-// void minmax (double _min[2], double _max[2]) const{
-// _min[0] = std::min(std::min(std::min(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
-// _min[1] = std::min(std::min(std::min(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
-// _max[0] = std::max(std::max(std::max(_q[0].x(),_q[1].x()),_q[2].x()),_q[3].x());
-// _max[1] = std::max(std::max(std::max(_q[0].y(),_q[1].y()),_q[2].y()),_q[3].y());
-// }
-// void print (FILE *f, int i){
-// fprintf(f,"SP(%g,%g,%g){%d};\n",_center.x(),_center.y(),0.0,i);
-// fprintf(f,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
-// _q[0].x(),_q[0].y(),0.0,
-// _q[1].x(),_q[1].y(),0.0,
-// _q[2].x(),_q[2].y(),0.0,
-// _q[3].x(),_q[3].y(),0.0,i,i,i,i);
-//
-// }
-//};
-//
-//struct my_wrapper {
-// bool _tooclose;
-// SPoint2 _p;
-// my_wrapper (SPoint2 sp) : _tooclose (false), _p(sp) {}
-//};
-//
-//struct smoothness_point_pair{
-// double rank;
-// surfacePointWithExclusionRegion* ptr;
-//};
-//
-//class compareSurfacePointWithExclusionRegionPtr_Smoothness
-//{
-// public:
-// inline bool operator () (const smoothness_point_pair &a, const smoothness_point_pair &b) const
-// {
-// if (a.rank == b.rank){
-// if(a.ptr->_distanceSummed > b.ptr->_distanceSummed) return false;
-// if(a.ptr->_distanceSummed < b.ptr->_distanceSummed) return true;
-// return a.ptr<b.ptr;
-// }
-// // else
-// return (a.rank < b.rank);
-// }
-//};
-//
-//
-//class compareSurfacePointWithExclusionRegionPtr
-//{
-// public:
-// inline bool operator () (const surfacePointWithExclusionRegion *a, const surfacePointWithExclusionRegion *b) const
-// {
-// if(a->_distanceSummed > b->_distanceSummed) return false;
-// if(a->_distanceSummed < b->_distanceSummed) return true;
-// return a<b;
-// }
-//};
-//
-//
-//
-//
-//bool rtree_callback(surfacePointWithExclusionRegion *neighbour,void* point){
-// my_wrapper *w = static_cast<my_wrapper*>(point);
-//
-// if (neighbour->inExclusionZone(w->_p)){
-// w->_tooclose = true;
-// return false;
-// }
-//
-// return true;
-//}
-//
-//bool inExclusionZone (SPoint2 &p,
-// RTree<surfacePointWithExclusionRegion*,double,2,double> &rtree,
-// std::vector<surfacePointWithExclusionRegion*> & all ){
-// // should assert that the point is inside the domain
-// if (!backgroundMesh::current()->inDomain(p.x(),p.y(),0)) return true;
-//
-// my_wrapper w (p);
-// double _min[2] = {p.x()-1.e-1, p.y()-1.e-1},_max[2] = {p.x()+1.e-1,p.y()+1.e-1};
-// rtree.Search(_min,_max,rtree_callback,&w);
-//
-// return w._tooclose;
-//
-// for (unsigned int i=0;i<all.size();++i){
-// if (all[i]->inExclusionZone(p)){
-// // printf("%g %g is in exclusion zone of %g %g\n",p.x(),p.y(),all[i]._center.x(),all[i]._center.y());
-// return true;
-// }
-// }
-// return false;
-//}
-
-
-
-
// assume a point on the surface, compute the 4 possible neighbors.
//
// ^ t2
@@ -219,15 +51,11 @@ bool compute4neighbors (GFace *gf, // the surface
SPoint2 newP[4][NUMDIR], // look into other directions
SMetric3 &metricField, FILE *crossf = 0) // the mesh metric
{
-
- //Range<double> rangeU = gf->parBounds(0);
- //Range<double> rangeV = gf->parBounds(1);
-
// we assume that v is on surface gf
-
+
// get the parameter of the point on the surface
reparamMeshVertexOnFace(v_center, gf, midpoint);
-
+
double L = backgroundMesh::current()->operator()(midpoint[0],midpoint[1],0.0);
// printf("L = %12.5E\n",L);
metricField = SMetric3(1./(L*L));
@@ -322,23 +150,6 @@ bool compute4neighbors (GFace *gf, // the surface
size_param_1 = size_param_2 = std::min (size_param_1,size_param_2);
}
- // printf("%12.5E %12.5E\n", size_param_1, size_param_2);
-
-
- // if (v_center->onWhat() != gf && gf->tag() == 3)
- // printf("M = (%g %g %g) L = %g %g LP = %g %g\n",metricField(0,0),metricField(1,1),metricField(0,1),l1,l2,size_param_1,size_param_2);
- //if (l1 == 0.0 || l2 == 0.0) printf("bouuuuuuuuuuuuh %g %g %g %g --- %g %g %g %g %g %g\n",l1,l2,t1.norm(),t2.norm(),s1.x(),s1.y(),s1.z(),s2.x(),s2.y(),s2.z());
-
- /* printf("%12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E %g %g %g %g %g %g %g %g %g %g %g\n",
- M*covar1[0]*covar1[0]+
- 2*E*covar1[1]*covar1[0]+
- N*covar1[1]*covar1[1],
- M*covar2[0]*covar2[0]+
- 2*E*covar2[1]*covar2[0]+
- N*covar2[1]*covar2[1]
- ,covar1[0],covar1[1],covar2[0],covar2[1],l1,l2,size_1,size_2,size_param_1,size_param_2,M,N,E,s1.x(),s1.y(),s2.x(),s2.y());*/
-
- // this is the rectangle in the parameter plane.
const double EPS = 1.e-7;
double r1 = EPS*(double)rand() / RAND_MAX;
double r2 = EPS*(double)rand() / RAND_MAX;
@@ -360,19 +171,14 @@ bool compute4neighbors (GFace *gf, // the surface
// a nonlinear problem in order to find a better approximation in the real
// surface
double ERR[4];
- for (int i=0;i<4;i++){ //
- // if (newPoint[i][0] < rangeU.low())newPoint[i][0] = rangeU.low();
- // if (newPoint[i][0] > rangeU.high())newPoint[i][0] = rangeU.high();
- // if (newPoint[i][1] < rangeV.low())newPoint[i][1] = rangeV.low();
- // if (newPoint[i][1] > rangeV.high())newPoint[i][1] = rangeV.high();
+ for (int i=0;i<4;i++){
GPoint pp = gf->point(SPoint2(newPoint[i][0], newPoint[i][1]));
double D = sqrt ((pp.x() - v_center->x())*(pp.x() - v_center->x()) +
(pp.y() - v_center->y())*(pp.y() - v_center->y()) +
(pp.z() - v_center->z())*(pp.z() - v_center->z()) );
ERR[i] = 100*fabs(D-L)/(D+L);
- // printf("L = %12.5E D = %12.5E ERR = %12.5E\n",L,D,100*fabs(D-L)/(D+L));
}
-
+
if (1 && goNonLinear){
surfaceFunctorGFace ss (gf); //
SVector3 dirs[4] = {t1*(-1.0),t2*(-1.0),t1*(1.0),t2*(1.0)}; //
@@ -966,30 +772,21 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
// printf("initially : %d vertices in the domain\n",vertices.size());
-
+
while(!fifo.empty()){
- //surfacePointWithExclusionRegion & parent = fifo.top();
- // surfacePointWithExclusionRegion * parent = fifo.front();
surfacePointWithExclusionRegion * parent = *fifo.begin();
- // fifo.pop();
fifo.erase(fifo.begin());
for (int dir=0;dir<NUMDIR;dir++){
- // printf("dir = %d\n",dir);
int countOK = 0;
for (int i=0;i<4;i++){
- // printf("i = %d %12.5E %12.5E \n",i,parent._p[i][dir].x(),parent._p[i][dir].y());
-
- // if (!w._tooclose){
if (!inExclusionZone (parent->_p[i][dir], rtree, vertices) ){
countOK++;
GPoint gp = gf->point(parent->_p[i][dir]);
MFaceVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
- // printf(" %g %g %g %g\n",parent._center.x(),parent._center.y(),gp.u(),gp.v());
SPoint2 midpoint;
compute4neighbors (gf, v, midpoint, goNonLinear, newp, metricField,crossf);
surfacePointWithExclusionRegion *sp =
new surfacePointWithExclusionRegion (v, newp, midpoint, metricField, parent);
- // fifo.push(sp);
fifo.insert(sp);
vertices.push_back(sp);
double _min[2],_max[2];
@@ -999,41 +796,31 @@ void packingOfParallelograms(GFace* gf, std::vector<MVertex*> &packed, std::vec
}
if (countOK)break;
}
- // printf("%d\n",vertices.size());
- }
- if (crossf){
- fprintf(crossf,"};\n");
- fclose (crossf);
- }
- // printf("done\n");
-
+ }
+ if (crossf){
+ fprintf(crossf,"};\n");
+ fclose (crossf);
+ }
// add the vertices as additional vertices in the
// surface mesh
- char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
- FILE *f = Fopen(ccc,"w");
- if(f) fprintf(f,"View \"\"{\n");
- for (unsigned int i=0;i<vertices.size();i++){
- // if(vertices[i]->_v->onWhat() != gf)
- if(f) vertices[i]->print(f,i);
- if(vertices[i]->_v->onWhat() == gf) {
- packed.push_back(vertices[i]->_v);
- metrics.push_back(vertices[i]->_meshMetric);
- SPoint2 midpoint;
- reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
- // fprintf(f,"TP(%22.15E,%22.15E,%g){%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E};\n",vertices[i]->_v->x(),vertices[i]->_v->y(),vertices[i]->_v->z(),
- // vertices[i]->_meshMetric(0,0),vertices[i]->_meshMetric(0,1),vertices[i]->_meshMetric(0,2),
- // vertices[i]->_meshMetric(1,0),vertices[i]->_meshMetric(1,1),vertices[i]->_meshMetric(1,2),
- // vertices[i]->_meshMetric(2,0),vertices[i]->_meshMetric(2,1),vertices[i]->_meshMetric(2,2));
- //fprintf(f,"SP(%22.15E,%22.15E,%g){1};\n",midpoint.x(),midpoint.y(),0.0);
- }
- delete vertices[i];
- }
- if(f){
- fprintf(f,"};");
- fclose(f);
+ char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
+ FILE *f = Fopen(ccc,"w");
+ if(f) fprintf(f,"View \"\"{\n");
+ for (unsigned int i=0;i<vertices.size();i++){
+ // if(vertices[i]->_v->onWhat() != gf)
+ if(f) vertices[i]->print(f,i);
+ if(vertices[i]->_v->onWhat() == gf) {
+ packed.push_back(vertices[i]->_v);
+ metrics.push_back(vertices[i]->_meshMetric);
+ SPoint2 midpoint;
+ reparamMeshVertexOnFace(vertices[i]->_v, gf, midpoint);
}
- // printf("packed.size = %d\n",packed.size());
- // delete rtree;
+ delete vertices[i];
+ }
+ if(f){
+ fprintf(f,"};");
+ fclose(f);
+ }
}