diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index e8bd05b0b061f6c825f89008d29b9da10b2e5156..5aafe141e37b02933acc4b60ae5350eec84eddd6 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -15,6 +15,9 @@
#include "GaussLegendre1D.h"
#include "Context.h"
#include "closestPoint.h"
+#if defined(HAVE_MESH)
+#include "meshGEdge.h"
+#endif
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
: GEntity(model, tag), _length(0.), _tooSmall(false), _cp(0),
@@ -627,3 +630,10 @@ void GEdge::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<doubl
ts.push_back(upts[p].t);
}
}
+
+void GEdge::mesh(bool verbose){
+#if defined(HAVE_MESH)
+ meshGEdge mesher;
+ mesher(this);
+#endif
+}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 009f4a488fd7c042aec95af9c3ff79bd71b81f99..9c96c4c66ad163d26d732edf9743b42ada6aa2f3 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -219,6 +219,7 @@ class GEdge : public GEntity{
virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
SPoint3 closestPoint (SPoint3 &p, double tolerance);
+ virtual void mesh(bool) ;
};
#endif
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 740cf21fa075dfbabe68169725bb6050c82c524a..ef26f57d22674fc622a6c119e5d0203f41c68e46 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -178,6 +178,10 @@ class GEntity {
virtual ~GEntity(){}
+
+ // mesh generation of the entity
+ virtual void mesh(bool verbose) {}
+
// delete the mesh data
virtual void deleteMesh(){}
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index f7230c902fc3fe1eee606e16b475458efb4ea3fe..739d988edc3a14edd895e19a4ae7f12d2dd57ba2 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -20,6 +20,7 @@
#include "OS.h"
#if defined(HAVE_MESH)
+#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "meshGFaceLloyd.h"
#include "BackgroundMeshTools.h"
@@ -1387,6 +1388,14 @@ bool GFace::fillPointCloud(double maxDist,
return true;
}
+void GFace::mesh(bool verbose) {
+#if defined(HAVE_MESH)
+ meshGFace mesher;
+ mesher (this, verbose);
+#endif
+}
+
+
void GFace::lloyd(int nbiter, int infn)
{
#if defined(HAVE_MESH) && defined(HAVE_BFGS)
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 237c5a8366cb7a305b409923fc635c44f2c0479a..f4fbcf28fb384cb99c2038e13e7e41ddb6fcc851 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -296,6 +296,9 @@ class GFace : public GEntity{
// tells if it's a sphere, and if it is, returns parameters
virtual bool isSphere(double &radius, SPoint3 ¢er) const { return false; }
+ // new interface for meshing
+ virtual void mesh(bool verbose);
+
// add layers of quads
void addLayersOfQuads(int nLayers, GVertex *start, double hmin, double factor);
@@ -326,7 +329,7 @@ class GFace : public GEntity{
int getCurvatureControlParameter () const;
void setCurvatureControlParameter(int);
virtual double getMeshSize() const { return meshAttributes.meshSize; }
-
+
struct {
mutable GEntity::MeshGenerationStatus status;
double worst_element_shape, best_element_shape, average_element_shape;
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index cc64b66fd395ed526ff9680a40091b0f85255c8e..178f05fb404112401b0e4259ef11fef26f8824c3 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -11,6 +11,7 @@
#include "GEdge.h"
#include "GFace.h"
#include "GFaceCompound.h"
+#include "discreteDiskFace.h"
#include "GmshMessage.h"
#include "StringUtils.h"
@@ -465,6 +466,13 @@ bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 ¶m,
return true;
}
+ if (gf->geomType() == GEntity::DiscreteDiskSurface ){
+ discreteDiskFace *gfc = (discreteDiskFace*) gf;
+ param = gfc->parFromVertex(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);
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 3bfbdacabf42e47be3e4b302a8bf234ab77d64b0..f67a4e2aab72663086c19827bc9db3bfb54264f6 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -19,7 +19,6 @@
#include "laplaceTerm.h"
#include "convexLaplaceTerm.h"
#include "convexCombinationTerm.h" // #
-#include "../Numeric/BasisFactory.h" // for derivatives
// The three things that are mandatory to manipulate octrees (octree in (u;v)).
static void discreteDiskFaceBB(void *a, double*mmin, double*mmax)
@@ -97,9 +96,29 @@ static bool orderVertices(const double &tot_length, const std::vector<MVertex*>
/*BUILDER*/
discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh) :
- GFace(GEntity::model(),0), _parent (gf)
+ GFace(gf->model(),123), _parent (gf)
{
- triangles = mesh;
+ std::map<MVertex*,MVertex*> v2v;
+ for (unsigned int i=0;i<mesh.size();i++){
+ MVertex *vs[3] = {NULL, NULL, NULL};
+ for (int j=0;j<3;j++){
+ MVertex *v = mesh[i]->getVertex(j);
+ if (v->onWhat() == gf) {
+ std::map<MVertex*,MVertex*> :: iterator it = v2v.find(v);
+ if (it == v2v.end()){
+ MFaceVertex *vv = new MFaceVertex ( v->x(), v->y(), v->z(), this, 0, 0);
+ v2v[v] = vv;
+ discrete_vertices.push_back(vv);
+ vs[j] = vv;
+ }
+ else vs[j] = it->second;
+ }
+ else vs[j] = v;
+ }
+ discrete_triangles.push_back(new MTriangle (vs[0], vs[1], vs[2]));
+ }
+
+ // triangles = mesh;
uv_kdtree = NULL;
kdtree = NULL;
// checkOrientation(); // mark, need to check after parametrization ? --> no ! (3D, cf. Seb)
@@ -108,16 +127,17 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh) :
orderVertices(_totLength, _U0, _coords);
parametrize();
buildOct();
- putOnView();
+ // putOnView();
}
+
void discreteDiskFace::getBoundingEdges()
{
// allEdges contains all edges of the boundary
std::set<MEdge,Less_Edge> allEdges;
- for(unsigned int i = 0; i < getNumMeshElements(); ++i){
- MElement *e = getMeshElement(i);
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ MElement *e = discrete_triangles[i];
for(int j = 0; j < e->getNumEdges() ; j++){
MEdge ed = e->getEdge(j);
std::set<MEdge,Less_Edge>::iterator it = allEdges.find(ed);
@@ -163,12 +183,12 @@ void discreteDiskFace::getBoundingEdges()
}
void discreteDiskFace::buildOct() const
-{ // mark ?
+{
SBoundingBox3d bb;
int count = 0;
- for(unsigned int i = 0; i < triangles.size(); ++i){
- MTriangle *t = triangles[i];
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ MTriangle *t = discrete_triangles[i];
//create bounding box
for(int j = 0; j < 3; j++){
std::map<MVertex*,SPoint3>::const_iterator itj = coordinates.find(t->getVertex(j));
@@ -196,8 +216,8 @@ void discreteDiskFace::buildOct() const
count = 0;
- for(unsigned int i = 0; i < triangles.size(); ++i){
- MTriangle *t = triangles[i];
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ MTriangle *t = discrete_triangles[i];
std::map<MVertex*,SPoint3>::const_iterator it0 =
coordinates.find(t->getVertex(0));
std::map<MVertex*,SPoint3>::const_iterator it1 =
@@ -242,21 +262,12 @@ void discreteDiskFace::buildOct() const
_ddft[count].tri = t;
- /*
- SVector3 dXdxi (_ddft[count].v2 - _ddft[count].v1); // constant
- SVector3 dXdeta(_ddft[count].v3 - _ddft[count].v1); // constant
- firstElemDerivatives[(MElement*)t] = Pair<SVector3,SVector3>(dXdxi,dXdeta);
- */
-
-
//compute first derivatives for every triangle , to disappear (erase)
double mat[2][2] = {{_ddft[count].p2.x() - _ddft[count].p1.x(),
_ddft[count].p3.x() - _ddft[count].p1.x()},
{_ddft[count].p2.y() - _ddft[count].p1.y(),
_ddft[count].p3.y() - _ddft[count].p1.y()}}; // modified higher
-
-
double inv[2][2];
inv2x2(mat, inv);
SVector3 dXdxi (_ddft[count].v2 - _ddft[count].v1); // constant
@@ -264,9 +275,7 @@ void discreteDiskFace::buildOct() const
SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]); // to be modified for higher order
SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]); // to be modified for higher order
firstElemDerivatives[(MElement*)t] = Pair<SVector3,SVector3>(dXdu,dXdv);
-
- // build ANN kdtree
nodes[count][0] = (it0->second.x() + it1->second.x() + it2->second.x())/3.0 ;
nodes[count][1] = (it0->second.y() + it1->second.y() + it2->second.y())/3.0 ;
nodes[count][2] = 0.0;
@@ -280,26 +289,24 @@ void discreteDiskFace::buildOct() const
// USELESS, laplacian
//smooth first derivatives at vertices
if(adjv.size() == 0){
- std::vector<MTriangle*> allTri;
- allTri.insert(allTri.end(), triangles.begin(), triangles.end() );
- buildVertexToTriangle(allTri, adjv);
+ std::vector<MTriangle*> allTri;
+ allTri.insert(allTri.end(), discrete_triangles.begin(), discrete_triangles.end() );
+ buildVertexToTriangle(allTri, adjv);
}
for(v2t_cont::iterator it = adjv.begin(); it!= adjv.end(); ++it){
- MVertex *v = it->first;
- std::vector<MElement*> vTri = it->second;
- SVector3 dXdu(0.0), dXdv(0.0);
- int nbTri = vTri.size();
- for (int j = 0; j < nbTri; j++){ // elements's contribution for a vertex
- dXdu += firstElemDerivatives[vTri[j]].first();
- dXdv += firstElemDerivatives[vTri[j]].second();
- }
- dXdu *= 1./nbTri;
- dXdv *= 1./nbTri;
- firstDerivatives[v] = Pair<SVector3, SVector3>(dXdu, dXdv);
+ MVertex *v = it->first;
+ std::vector<MElement*> vTri = it->second;
+ SVector3 dXdu(0.0), dXdv(0.0);
+ int nbTri = vTri.size();
+ for (int j = 0; j < nbTri; j++){ // elements's contribution for a vertex
+ dXdu += firstElemDerivatives[vTri[j]].first();
+ dXdv += firstElemDerivatives[vTri[j]].second();
+ }
+ dXdu *= 1./nbTri;
+ dXdv *= 1./nbTri;
+ firstDerivatives[v] = Pair<SVector3, SVector3>(dXdu, dXdv);
}
-
- //build ANN kdtree
kdtree = new ANNkd_tree(nodes, count, 3);
}
@@ -326,8 +333,8 @@ bool discreteDiskFace::parametrize() const
myAssemblerV.fixVertex(v, 0, 1, sin(theta));
}
- for(size_t i = 0; i < triangles.size(); ++i){
- MTriangle *t = triangles[i];
+ for(size_t i = 0; i < discrete_triangles.size(); ++i){
+ MTriangle *t = discrete_triangles[i];
myAssemblerU.numberVertex(t->getVertex(0), 0, 1);
myAssemblerU.numberVertex(t->getVertex(1), 0, 1);
@@ -348,12 +355,9 @@ bool discreteDiskFace::parametrize() const
convexLaplaceTerm mappingU(0, 1, &ONE);
convexLaplaceTerm mappingV(0, 1, &ONE);
- /*
- convexCombinationTerm mappingU(0,1,&ONE);
- convexCombinationTerm mappingV(0,1,&ONE);
- */
- for(unsigned int i = 0; i < triangles.size(); ++i){
- SElement se(triangles[i]);
+
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ SElement se(discrete_triangles[i]);
mappingU.addToMatrix(myAssemblerU, &se);
mappingV.addToMatrix(myAssemblerV, &se);
}
@@ -393,12 +397,12 @@ bool discreteDiskFace::parametrize() const
}
-// ------
-
-void discreteDiskFace::getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v)const{ // does it change with higher order ? no, I don't think
-
+void discreteDiskFace::getTriangleUV(const double u,const double v,
+ discreteDiskFaceTriangle **mt,
+ double &_u, double &_v)const{
double uv[3] = {u,v,0.};
*mt = (discreteDiskFaceTriangle*) Octree_Search(uv,oct);
+ if (!(*mt))return;
double M[2][2],X[2],R[2];
const SPoint3 p0 = (*mt)->p1;
@@ -416,7 +420,6 @@ void discreteDiskFace::getTriangleUV(const double u,const double v,discreteDiskF
}
-
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
@@ -425,47 +428,33 @@ GPoint discreteDiskFace::point(const double par1, const double par2) const
double par[2] = {par1,par2};
discreteDiskFaceTriangle* dt;
getTriangleUV(par1,par2,&dt,U,V);
+ if (!dt) {
+ GPoint gp = GPoint(1.e22,1.e22,1.e22,_parent,par);
+ gp.setNoSuccess();
+ return gp;
+ }
-
SPoint3 p = dt->v1*(1.-U-V) + dt->v2*U + dt->v3*V;
return GPoint(p.x(),p.y(),p.z(),_parent,par);
}
-// (x;y;z) |-> (u;v)
-SPoint2 discreteDiskFace::parFromPoint(const MVertex* v)
+SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
{
- double xyz[3] = {v->x(),v->y(),v->z()};
- double uvw[3] = {0.,0.,0.};
-
- unsigned int i;
- for(i = 0; i<triangles.size(); ++i){ // loop on every mesh triangles
-
- MTriangle* t = _ddft[i].tri;
-
- t->xyz2uvw(xyz,uvw); // (x;y;z) |-> (xsi,eta) [ref triangle]
-
- if (t->isInside(uvw[0],uvw[1],uvw[2])){ // the mesh vertex belongs to this triangle
-
- discreteDiskFaceTriangle* theTri = &_ddft[i];
-
- SPoint3 p = theTri->p1*(1-uvw[0]-uvw[1]) + theTri->p2*uvw[0] + theTri->p3*uvw[1];
-
- double U = p.x();
- double V = p.y();
-
- break;
- return SPoint2(U,V);
-
- }
-
+ if (v->onWhat()->dim()==2)Msg::Fatal("discreteDiskFace::parFromVertex should not be called with a vertex that is classified on a surface");
+ std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
+ if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
+ // The 1D mesh has been re-done
+ if (v->onWhat()->dim()==1){
+ Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
+ // get the discrete edge on which v is classified
+ // get the two ending vertices A and B on both sides of v A-------v----------B
+ // interpolate between the two uv coordinates of A and B
+ }
+ else if (v->onWhat()->dim()==0){
+ Msg::Fatal("discreteDiskFace::parFromVertex vertex classified on a model vertex that is not part of the face");
}
- if (i==triangles.size())
- Msg::Error("discreteDiskFace::parFromPoint << no triangle contains this vertex !");
-
-
- return SPoint2(0.,0.);
}
SVector3 discreteDiskFace::normal(const SPoint2 ¶m) const
@@ -485,7 +474,6 @@ double discreteDiskFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVec
return false;
}
-
Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
{
@@ -493,38 +481,32 @@ Pair<SVector3, SVector3> discreteDiskFace::firstDer(const SPoint2 ¶m) const
discreteDiskFaceTriangle* ddft;
getTriangleUV(param.x(),param.y(),&ddft,U,V);
-
MTriangle* tri = NULL;
if (ddft) tri = ddft->tri;
- else Msg::Error("discreteDiskFace::firstDer << triangle not found");
-
- /*
- const nodalBasis* base = BasisFactory::getNodalBasis(MSH_TRI_3);
- double grads[2][3];
- double mat[2][2];
- double inv[2][2];
- base->df(U,V,0.,grads);
-
- for (int i=0; i<2; i++)
- for (int j=0; j<2; j++)
- mat[i][j] = grads[i][j];
-
- inv2x2(mat,inv);
+ else Msg::Error("discreteDiskFace::firstDer << triangle not found");
+
+ std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it0 =
+ firstDerivatives.find(tri->getVertex(0));
+ if (it0 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (0) has no firstDerivatives",
+ tri->getVertex(0)->getNum());
+ std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it1 =
+ firstDerivatives.find(tri->getVertex(1));
+ if (it1 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (1) has no firstDerivatives",
+ tri->getVertex(1)->getNum());
+ std::map<MVertex*, Pair<SVector3,SVector3> >::iterator it2 =
+ firstDerivatives.find(tri->getVertex(2));
+ if (it2 == firstDerivatives.end()) Msg::Fatal ("Vertex %d (2) has no firstDerivatives",
+ tri->getVertex(2)->getNum());
+
- SVector3 dXdxi = firstElemDerivatives[tri].first();
- SVector3 dXdeta = firstElemDerivatives[tri].second();
+ SVector3 dXdu1 = it0->second.first();
+ SVector3 dXdu2 = it1->second.first();
+ SVector3 dXdu3 = it2->second.first();
- SVector3 dXdu = dXdxi*inv[0][0] + dXdeta*inv[1][0];
- SVector3 dXdv = dXdxi*inv[0][1] + dXdeta*inv[1][1];
- */
+ SVector3 dXdv1 = it0->second.second();
+ SVector3 dXdv2 = it1->second.second();
+ SVector3 dXdv3 = it2->second.second();
-
- SVector3 dXdu1 = firstDerivatives[tri->getVertex(0)].first();
- SVector3 dXdu2 = firstDerivatives[tri->getVertex(1)].first();
- SVector3 dXdu3 = firstDerivatives[tri->getVertex(2)].first();
- SVector3 dXdv1 = firstDerivatives[tri->getVertex(0)].second();
- SVector3 dXdv2 = firstDerivatives[tri->getVertex(1)].second();
- SVector3 dXdv3 = firstDerivatives[tri->getVertex(2)].second();
SVector3 dXdu = dXdu1*(1.-U-V) + dXdu2*U + dXdu3*V;
SVector3 dXdv = dXdv1*(1.-U-V) + dXdv2*U + dXdv3*V;
@@ -544,8 +526,8 @@ void discreteDiskFace::secondDer(const SPoint2 ¶m,
void discreteDiskFace::buildAllNodes()
{
// allNodes contains all mesh-nodes
- for(unsigned int i = 0; i < getNumMeshElements(); ++i){
- MElement *e = getMeshElement(i);
+ for(unsigned int i = 0; i < discrete_triangles.size(); ++i){
+ MElement *e = discrete_triangles[i];
for(int j = 0; j < e->getNumVertices() ; j++){
MVertex* ev = e->getVertex(j);
std::set<MVertex*>::iterator it = allNodes.find(ev); // several times the same nodes !
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index 2bd4ac92d40c10edbb25d43093714581c9f11d37..f82f6437097dae3e7aaeae3b1c29836d9e4b4947 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -49,7 +49,7 @@ class discreteDiskFace : public GFace {
virtual ~discreteDiskFace() {triangles.clear();}
void getTriangleUV(const double u,const double v,discreteDiskFaceTriangle **mt, double &_u, double &_v) const;
GPoint point(double par1, double par2) const;
- SPoint2 parFromPoint(const MVertex* v);
+ SPoint2 parFromVertex(MVertex *v) const;
SVector3 normal(const SPoint2&) const;
double curvatureMax(const SPoint2&) const;
double curvatures(const SPoint2&,SVector3*,SVector3*,double*,double*) const;
@@ -58,6 +58,11 @@ class discreteDiskFace : public GFace {
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
GEntity::GeomType geomType() const { return DiscreteDiskSurface; }
protected:
+ //------------------------------------------------
+ // a copy of the mesh that should not be destroyed
+ std::vector<MTriangle*> discrete_triangles;
+ std::vector<MVertex*> discrete_vertices;
+ //------------------------------------------------
double _totLength;
std::map<double,std::vector<MVertex*> > _loops;
std::vector<MVertex*> _U0; // dirichlet's bc's
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 4b5cc2652a70b97b2b2a3500b16bb1a8dc7a6539..793c30ece00b2f9431f1d97da5ff2c67e542043a 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -21,7 +21,9 @@
#include "OS.h"
#include "Curvature.h"
#include "GEdgeCompound.h"
-
+#if defined(HAVE_MESH)
+#include "meshGEdge.h"
+#endif
discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
: GEdge(model, num, _v0, _v1)
{
@@ -377,30 +379,15 @@ void discreteEdge::computeNormals () const
for ( ; itn != _normals.end(); ++itn){
itn->second.normalize();
}
-
-// //smooth normals
-// smooth_normals *normals = 0;
-// normals = new smooth_normals(180.);
-
-// for ( itn = _normals.begin(); itn != _normals.end(); ++itn){
-// MVertex *v = itn->first;
-// SVector3 n = itn->second;
-// normals->add(v->x(),v->y(),v->z(),n.x(),n.y(),n.z());
-// }
-// for ( itn = _normals.begin(); itn != _normals.end(); ++itn){
-// double nx,ny,nz;
-// MVertex *v = itn->first;
-// normals->get(v->x(),v->y(),v->z(),nx,ny,nz);
-// itn->second = SVector3(nx,ny,nz);
-// }
-// delete normals;
}
+/// WE REWRITE WHAT FOLLOWS : THIS WAS WRONG (PAB -JFR)
+
bool discreteEdge::getLocalParameter(const double &t, int &iLine,
double &tLoc) const
{
- if(_pars.empty()) return false;
- for (iLine = 0; iLine < (int)lines.size(); iLine++){
+
+ for (iLine = 0; iLine < (int)discrete_lines.size(); iLine++){
double tmin = _pars[iLine];
double tmax = _pars[iLine+1];
if (t >= tmin && t <= tmax){
@@ -422,55 +409,11 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const bool LINEARMESH = true; //false;
-
- if (LINEARMESH){
- //linear Lagrange mesh
- x = vB->x() + tLoc * (vE->x()- vB->x());
- y = vB->y() + tLoc * (vE->y()- vB->y());
- z = vB->z() + tLoc * (vE->z()- vB->z());
- return GPoint(x,y,z);
- }
- else{
- //curved PN triangle
-
- if (_normals.empty() ) computeNormals();
-
- const SVector3 n1 = _normals[vB];
- const SVector3 n2 = _normals[vE];
-
- SPoint3 v1(vB->x(),vB->y(),vB->z());
- SPoint3 v2(vE->x(),vE->y(),vE->z());
-
- SVector3 b300,b030,b003;
- SVector3 b210,b120;
- double w12,w21;
-
- b300 = v1;
- b030 = v2;
-
- w12 = dot(v2 - v1, n1);
- w21 = dot(v1 - v2, n2);
- b210 = (2*v1 + v2 -w12*n1)*0.333;
- b120 = (2*v2 + v1 -w21*n2)*0.333;
-
- // tagged PN trinagles (sigma=1)
- double theta = 0.0;
- SVector3 d1 = v1+.33*(1-theta)*(v2-v1);
- SVector3 d2 = v2+.33*(1-theta)*(v1-v2);
- SVector3 X1 = 1/norm(n1)*n1;
- SVector3 X2 = 1/norm(n2)*n2;
- b210 = d1 - dot(X1,d1-v1)*X1;
- b120 = d2 - dot(X2,d2-v2)*X2;
-
- double U = tLoc;
- double W = 1-U;
- SVector3 point = b300*W*W*W+b030*U*U*U+b210*3.*W*W*U+b120*3.*W*U*U;
-
- SPoint3 PP(point.x(), point.y(), point.z());
- return GPoint(PP.x(),PP.y(),PP.z());
- }
-
+ //linear Lagrange mesh
+ x = vB->x() + tLoc * (vE->x()- vB->x());
+ y = vB->y() + tLoc * (vE->y()- vB->y());
+ z = vB->z() + tLoc * (vE->z()- vB->z());
+ return GPoint(x,y,z);
}
SVector3 discreteEdge::firstDer(double par) const
@@ -527,9 +470,53 @@ double discreteEdge::curvatures(const double par, SVector3 *dirMax, SVector3 *di
Range<double> discreteEdge::parBounds(int i) const
{
- return Range<double>(0, lines.size());
+ return Range<double>(0, 1);
+}
+
+void discreteEdge::createGeometry(){
+ if (discrete_lines.empty()){
+ createTopo();
+ // copy the mesh
+ std::map<MVertex*,MVertex*> v2v;
+ for (unsigned int i = 0; i < mesh_vertices.size(); i++){
+ MVertex *v = new MVertex(mesh_vertices[i]->x(),mesh_vertices[i]->y(),mesh_vertices[i]->z());
+ v2v[mesh_vertices[i]] = v;
+ discrete_vertices.push_back(v);
+ }
+ for (unsigned int i = 0; i < lines.size(); i++){
+ MVertex *v0 = lines[i]->getVertex(0);
+ MVertex *v1 = lines[i]->getVertex(1);
+ discrete_lines.push_back(new MLine(v0,v1));
+ }
+ // compute parameters
+ double L = 0.0;
+ _pars.push_back(L);
+ for (unsigned int i = 0; i < discrete_lines.size(); i++){
+ MVertex *v0 = discrete_lines[i]->getVertex(0);
+ MVertex *v1 = discrete_lines[i]->getVertex(1);
+ L += distance(v0,v1);
+ _pars.push_back(L);
+ }
+ for (unsigned int i = 0; i < _pars.size(); i++){
+ _pars[i] /= L;
+ }
+ }
}
+
+void discreteEdge::mesh(bool verbose){
+#if defined(HAVE_MESH)
+ // copy the mesh into geometrical entities
+ // FIXME
+ return;
+ createGeometry();
+ meshGEdge mesher;
+ mesher(this);
+#endif
+}
+
+
+
void discreteEdge::writeGEO(FILE *fp)
{
if (getBeginVertex() && getEndVertex())
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 0d80563c3a014240bbc1853844161a5011fb28a9..6a09a39555eb8de914f19a31531d8cc526eb7a9b 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -16,6 +16,8 @@ class discreteEdge : public GEdge {
std::vector<int> _orientation;
std::map<MVertex*, MLine*> boundv;
bool createdTopo;
+ std::vector<MVertex*> discrete_vertices;
+ std::vector<MLine*> discrete_lines;
public:
discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
virtual ~discreteEdge() {}
@@ -30,13 +32,15 @@ class discreteEdge : public GEdge {
bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
void parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
- std::less<MVertex*> > > &face2Verts,
- std::map<GRegion*, std::map<MVertex*, MVertex*,
- std::less<MVertex*> > > ®ion2Vert);
+ std::less<MVertex*> > > &face2Verts,
+ std::map<GRegion*, std::map<MVertex*, MVertex*,
+ std::less<MVertex*> > > ®ion2Vert);
void orderMLines();
void setBoundVertices();
void createTopo();
+ void createGeometry();
void computeNormals () const;
+ virtual void mesh(bool) ;
void writeGEO(FILE *fp);
};
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index b730bacba4f342260657e5851aded544d55e25c2..46aa13a4105d573b8f9e0e3520d658662d941bbd 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -140,12 +140,34 @@ void discreteFace::secondDer(const SPoint2 ¶m,
// FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
void discreteFace::createAtlas() {
+ if (!_atlas.empty())return;
// parametrization is done here !!!
discreteDiskFace *df = new discreteDiskFace (this, triangles);
df->replaceEdges(l_edges);
_atlas.push_back(df);
}
+void discreteFace::gatherMeshes() {
+ for (unsigned int i=0;i<triangles.size();i++)delete triangles[i];
+ for (unsigned int i=0;i<mesh_vertices.size();i++)delete mesh_vertices[i];
+ 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());
+ }
+}
+
+void discreteFace::mesh(bool verbose) {
+#if defined(HAVE_MESH)
+ createAtlas();
+ for (unsigned int i=0;i<_atlas.size();i++)
+ _atlas[i]->mesh(verbose);
+ gatherMeshes();
+#endif
+ meshStatistics.status = GFace::DONE;
+}
+
// delete all discrete disk faces
//void discreteFace::deleteAtlas() {
//}
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index c68f55074507930859b1d9babc8d0fa10ac187d8..f8a91f1a91e8517547354969539568f171be301a 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -35,6 +35,8 @@ class discreteFace : public GFace {
void findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
void writeGEO(FILE *fp);
void createAtlas();
+ void gatherMeshes();
+ virtual void mesh (bool verbose);
std::vector<discreteDiskFace*> _atlas;
};
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 55adba595236d8e52d3c629c6ee23df6eba23617..54f8e06a976d093fa93e3f17e572686846a8dc7a 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -242,11 +242,11 @@ static void Mesh1D(GModel *m)
int nIter = 0, nTot = m->getNumEdges();
while(1){
- meshGEdge mesher;
+ // meshGEdge mesher;
int nPending = 0;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
if ((*it)->meshStatistics.status == GEdge::PENDING){
- mesher(*it);
+ (*it)->mesh(true);
nPending++;
}
if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 1D...");
@@ -354,8 +354,8 @@ static void Mesh2D(GModel *m)
for(size_t K = 0 ; K < temp.size() ; K++){
if (temp[K]->meshStatistics.status == GFace::PENDING){
backgroundMesh::current()->unset();
- meshGFace mesher(true);
- mesher(temp[K]);
+ // meshGFace mesher(true);
+ temp[K]->mesh(true);
#if defined(HAVE_BFGS)
if(CTX::instance()->mesh.optimizeLloyd){
if (temp[K]->geomType()==GEntity::CompoundSurface ||
@@ -392,8 +392,8 @@ static void Mesh2D(GModel *m)
it != cf.end(); ++it){
if ((*it)->meshStatistics.status == GFace::PENDING){
backgroundMesh::current()->unset();
- meshGFace mesher(true);
- mesher(*it);
+ // meshGFace mesher(true);
+ (*it)->mesh(true);
#if defined(HAVE_BFGS)
if(CTX::instance()->mesh.optimizeLloyd){
if ((*it)->geomType()==GEntity::CompoundSurface ||
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index b3d05859db0f66f2f3880e25d4703a6732b35674..59f2412b3bfe9431bc2c96b1730eea778634167d 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -327,6 +327,11 @@ void computeAdjacencies (Tet *t, int iFace, connSet &faceToTet){
}
}
+bool edgeSwaps(tetContainer &T, int myThread)
+{
+ // TODO
+}
+
void edgeBasedRefinement (const int numThreads,
const int nptsatonce,
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 91bc4bf2a1366b6075db2bee26671f35220055eb..652b7d351faa7df64a395ab261e4352df00b202b 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -340,7 +340,7 @@ void meshGEdge::operator() (GEdge *ge)
ge->model()->setCurrentMeshEntity(ge);
- if(ge->geomType() == GEntity::DiscreteCurve) return;
+ // if(ge->geomType() == GEntity::DiscreteCurve) return;
if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
if(ge->meshAttributes.method == MESH_NONE) return;
if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 1cd649da3449be17b65459d78cf5639f297ec188..d42f09ff280d32cbb75e42c8c722b50548bea67a 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -991,6 +991,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
++ite;
}
+
if(boundary.size()){
Msg::Error("The 1D mesh seems not to be forming a closed loop");
gf->meshStatistics.status = GFace::FAILED;
@@ -1073,8 +1074,6 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
count++;
}
- // here check if some boundary layer nodes should be added
-
bbox.makeCube();
// use a divide & conquer type algorithm to create a triangulation.
@@ -1118,6 +1117,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
doc.points[points.size() + ip].adjacent = 0;
doc.points[points.size() + ip].data = pp;
}
+
+
// Use "fast" inhouse recursive algo to generate the triangulation.
// At this stage the triangulation is not what we need
@@ -1133,6 +1134,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
}
Msg::Debug("Meshing of the convex hull (%d points) done", points.size());
+
for(int i = 0; i < doc.numTriangles; i++) {
int a = doc.triangles[i].a;
int b = doc.triangles[i].b;
@@ -1465,7 +1467,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// gf->deleteMesh() would also destroy e.g. the data in a compound face, which
// we should not do)
gf->GFace::deleteMesh();
-
+
Msg::Debug("Starting to add internal points");
// start mesh generation
if(!algoDelaunay2D(gf) && !onlyInitialMesh){
@@ -2386,16 +2388,18 @@ void meshGFace::operator() (GFace *gf, bool print)
//-----------------------------------------------------------
// FIXME PAB & JFR FIRST IMPLEMENTATION OF THE COMPOUND REPLACEMENT
- if(gf->geomType() == GEntity::DiscreteSurface) {
- discreteFace *df = dynamic_cast<discreteFace*> (gf);
- if (df) {
- df->createAtlas();
- for (unsigned int i=0;i<df->_atlas.size();i++){
- (*this)(df->_atlas[i],print);
- }
- }
- return;
- }
+ // if(gf->geomType() == GEntity::DiscreteSurface) {
+ // discreteFace *df = dynamic_cast<discreteFace*> (gf);
+ // if (df) {
+ // df->createAtlas();
+ // for (unsigned int i=0;i<df->_atlas.size();i++){
+ // (*this)(df->_atlas[i],print);
+ // }
+ // df->gatherMeshes();
+ // }
+ // gf->meshStatistics.status = GFace::DONE;
+ // return;
+ // }
//-----------------------------------------------------------
//-----------------------------------------------------------
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index fbb46a1a60121aaae729a3a7dee94ac3edca5bea..6542f9a54d9fb5415342d9d75693bf1b59528223 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -592,7 +592,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> ®ions)
}
// uncomment this to test the new code
-#define NEW_CODE
+//#define NEW_CODE
static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> ®ions) {
GRegion *gr = regions[0];