diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index fe7d2bf6b6b1e7221fcb5ee5c0a966e64981827d..8762efc0d71c8e99b2396c4d6fad72f25c1794cd 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -875,9 +875,9 @@ void GetOptions(int argc, char *argv[])
i++;
CTX::instance()->mesh.recombinationTestNoGreedyStrat = 1;
}
- else if(!strcmp(argv[i] + 1, "newstrat")) {
+ else if(!strcmp(argv[i] + 1, "meshdiscrete")) {
i++;
- CTX::instance()->mesh.recombinationTestNewStrat = 1;
+ CTX::instance()->meshDiscrete = 1;
}
else if(!strcmp(argv[i] + 1, "format") || !strcmp(argv[i] + 1, "f")) {
i++;
diff --git a/Common/Context.cpp b/Common/Context.cpp
index 8da9ca1132b09c6442f5339a9beca92f754f8d6f..c51fe053896cddec207a0f2dcef231661e17c5e6 100644
--- a/Common/Context.cpp
+++ b/Common/Context.cpp
@@ -41,6 +41,7 @@ CTX::CTX() : gamepad(0)
homeDir += "/";
batch = batchAfterMesh = 0;
+ meshDiscrete=0;
outputFileName = "";
bgmFileName = "";
createAppendMeshStatReport = 0;
diff --git a/Common/Context.h b/Common/Context.h
index 0be231c191aba8f171498228f838f4a744deb96c..a0f6c10930f300fe0e181db8f58c84dddca5b495 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -138,6 +138,8 @@ class CTX {
int batch;
// batch operations to apply after meshing (1: partition mesh)
int batchAfterMesh;
+ // mesh discrete faces / edges
+ int meshDiscrete;
// initial menu (0: automatic, 1: geom, 2: mesh, 3: solver, 4: post)
int initialContext;
// never popup dialogs in scripts (use default values instead)?
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index bd1b1b1424f206d90c0461903105547591eff3eb..4d3736589b84c6736cb0f41bd1fb0a0f9968d650 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1208,8 +1208,6 @@ StringXNumber MeshOptions_Number[] = {
"Depth start" },
{ F|O, "RecombinationTestNoGreedyStrat" , opt_mesh_recombination_no_greedy_strat , 0 ,
"No greedy (global) strategies" },
- { F|O, "RecombinationTestNewStrat" , opt_mesh_recombination_new_strat , 0 ,
- "New strategies" },
{ F|O, "RemeshAlgorithm" , opt_mesh_remesh_algo , 0 ,
"Remeshing algorithm (0=no split, 1=automatic, 2=automatic only with metis)" },
{ F|O, "RemeshParametrization" , opt_mesh_remesh_param , 4 ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index d03c46c379d19fc0c64e84dab1549fc14f0f258d..2b85b0ef08015558e324fbd472f73967664c0650 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -5885,14 +5885,6 @@ double opt_mesh_recombination_no_greedy_strat(OPT_ARGS_NUM)
return CTX::instance()->mesh.recombinationTestNoGreedyStrat;
}
-double opt_mesh_recombination_new_strat(OPT_ARGS_NUM)
-{
- if(action & GMSH_SET){
- CTX::instance()->mesh.recombinationTestNewStrat = (int)val;
- }
- return CTX::instance()->mesh.recombinationTestNewStrat;
-}
-
double opt_mesh_remesh_algo(OPT_ARGS_NUM)
{
if(action & GMSH_SET){
diff --git a/Common/Options.h b/Common/Options.h
index 1757b6b53d05fa296b130ee5a0b64f2b069b456d..e2ad04ad598f277ba95dab40fd0aeae6d7f9883a 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -480,7 +480,6 @@ double opt_mesh_flexible_transfinite(OPT_ARGS_NUM);
double opt_mesh_do_recombination_test(OPT_ARGS_NUM);
double opt_mesh_recombination_test_start(OPT_ARGS_NUM);
double opt_mesh_recombination_no_greedy_strat(OPT_ARGS_NUM);
-double opt_mesh_recombination_new_strat(OPT_ARGS_NUM);
double opt_mesh_remesh_algo(OPT_ARGS_NUM);
double opt_mesh_remesh_param(OPT_ARGS_NUM);
double opt_mesh_algo_subdivide(OPT_ARGS_NUM);
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index c3bbcf902223147ed2603fa2e5ed768893761950..7e8347dc917ff708c4251b46aca77375a89e346f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1271,20 +1271,29 @@ static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &el
for(int j = 0; j < elements[i]->getNumVertices(); j++){
if (force || !elements[i]->getVertex(j)->onWhat() ||
elements[i]->getVertex(j)->onWhat()->dim() > ge->dim())
- elements[i]->getVertex(j)->setEntity(ge);
+ elements[i]->getVertex(j)->setEntity(ge);
}
}
}
void GModel:: _createGeometryOfDiscreteEntities() {
- for(eiter it = firstEdge(); it != lastEdge(); ++it){
- if((*it)->geomType() == GEntity::DiscreteCurve) {
- discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
- if (!de)Msg::Fatal("no fun");
- de->createGeometry();
+ if (CTX::instance()->meshDiscrete){
+ Msg::Info("creating the geometry of discrete curves");
+ for(eiter it = firstEdge(); it != lastEdge(); ++it){
+ if((*it)->geomType() == GEntity::DiscreteCurve) {
+ discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
+ if (!de)Msg::Fatal("no fun");
+ de->createGeometry();
+ }
+ }
+ Msg::Info("creating the geometry of discrete surfaces");
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ if((*it)->geomType() == GEntity::DiscreteSurface) {
+ discreteFace *df = dynamic_cast<discreteFace*> (*it);
+ if (!df)Msg::Fatal("no fun");
+ df->createGeometry();
+ }
}
- }
- for(fiter it = firstFace(); it != lastFace(); ++it){
}
}
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 1115b54deaef1d71d0cefe1c9a32e8852dc4c0c1..4c986a60a48c28d16c47a79612f1e42f9e799b72 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -159,6 +159,7 @@ int GModel::readMSH(const std::string &name)
char str[256] = "";
+
// detect prehistoric MSH files
if(!fgets(str, sizeof(str), fp)){ fclose(fp); return 0; }
if(!strncmp(&str[1], "NOD", 3) || !strncmp(&str[1], "NOE", 3)){
@@ -353,7 +354,6 @@ int GModel::readMSH(const std::string &name)
_vertexMapCache.clear();
}
}
-
// $Elements section
else if(!strncmp(&str[1], "Elements", 8)) {
if(!fgets(str, sizeof(str), fp)){ fclose(fp); return 0; }
@@ -439,7 +439,7 @@ int GModel::readMSH(const std::string &name)
// entity does not exist, create a new (discrete) one.
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
-
+
// associate the correct geometrical entity with each mesh vertex
_associateEntityWithMeshVertices();
@@ -450,6 +450,7 @@ int GModel::readMSH(const std::string &name)
_storeVerticesInEntities(_vertexMapCache);
_createGeometryOfDiscreteEntities() ;
+
fclose(fp);
return postpro ? 2 : 1;
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 355e81eb9bba0c696388265095e1141c1b32df0b..810db951197b9700f93bf99452e618ee4898f02b 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -11,6 +11,7 @@
#include "GmshMessage.h"
#include "Octree.h"
#include "discreteDiskFace.h"
+#include "discreteEdge.h"
#include "MTriangle.h"
#include "MEdge.h"
#include "Geo.h"
@@ -116,6 +117,13 @@ discreteDiskFace::discreteDiskFace(GFace *gf, std::vector<MTriangle*> &mesh) :
}
else vs[j] = it->second;
}
+ else if (v->onWhat()->dim() == 1){
+ if (v->onWhat()->geomType() == DiscreteCurve){
+ discreteEdge *de = dynamic_cast<discreteEdge*> (v->onWhat());
+ vs[j] = de->getGeometricalVertex(v);
+ }
+ else Msg::Fatal("fatality");
+ }
else vs[j] = v;
}
discrete_triangles.push_back(new MTriangle (vs[0], vs[1], vs[2]));
@@ -449,10 +457,21 @@ SPoint2 discreteDiskFace::parFromVertex(MVertex *v) const
if(it != coordinates.end()) return SPoint2(it->second.x(),it->second.y());
// The 1D mesh has been re-done
if (v->onWhat()->dim()==1){
+ if (v->onWhat()->geomType() == DiscreteCurve){
+ discreteEdge *de = dynamic_cast<discreteEdge*> (v->onWhat());
+ if (de){
+ MVertex *v1,*v2;
+ double xi;
+ de->interpolateInGeometry (v,&v1,&v2,xi);
+ std::map<MVertex*,SPoint3>::iterator it1 = coordinates.find(v1);
+ std::map<MVertex*,SPoint3>::iterator it2 = coordinates.find(v2);
+ if(it1 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
+ if(it2 == coordinates.end()) Msg::Fatal("FIXME TO DO %d %s",__LINE__,__FILE__);
+ return SPoint2(it1->second.x(),it1->second.y()) * (1.-xi) +
+ SPoint2(it2->second.x(),it2->second.y()) * xi;
+ }
+ }
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");
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 5ccf30023d59d9ba6e58672a9a996f11478f19aa..3b0b037410f47b15ea234241e5bbd6f842158bd2 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -23,6 +23,7 @@
#include "GEdgeCompound.h"
#if defined(HAVE_MESH)
#include "meshGEdge.h"
+#include "Context.h"
#endif
discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
: GEdge(model, num, _v0, _v1)
@@ -42,6 +43,7 @@ void discreteEdge::createTopo()
}
}
+// FULL OF BUGS !!!!!!
void discreteEdge::orderMLines()
{
//printf("ordering line %d\n", tag());
@@ -143,8 +145,6 @@ void discreteEdge::orderMLines()
mesh_vertices.end()) mesh_vertices.push_back(v1);
if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v2) ==
mesh_vertices.end()) mesh_vertices.push_back(v2);
- v1->setEntity(this);
- v2->setEntity(this);
}
//special case reverse orientation
@@ -413,7 +413,7 @@ GPoint discreteEdge::point(double par) const
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);
+ return GPoint(x,y,z,this,par);
}
SVector3 discreteEdge::firstDer(double par) const
@@ -475,12 +475,11 @@ Range<double> discreteEdge::parBounds(int i) const
void discreteEdge::createGeometry(){
if (discrete_lines.empty()){
- Msg::Info("creating the geometry of discrete curve %d",tag());
createTopo();
// copy the mesh
- std::map<MVertex*,MVertex*> v2v;
for (unsigned int i = 0; i < mesh_vertices.size(); i++){
- MVertex *v = new MEdgeVertex(mesh_vertices[i]->x(),mesh_vertices[i]->y(),mesh_vertices[i]->z(),this,(double)(i+1));
+ MEdgeVertex *v = new MEdgeVertex(mesh_vertices[i]->x(),mesh_vertices[i]->y(),mesh_vertices[i]->z(),this,(double)(i+1));
+ // printf("%3d %p\n",tag(),mesh_vertices[i]);
v2v [mesh_vertices[i]] = v;
}
@@ -497,8 +496,6 @@ void discreteEdge::createGeometry(){
discrete_lines.push_back(new MLine(v01,v00));
}
}
-
-
// compute parameters and recompute the vertices
_pars.push_back(0.0);
for (unsigned int i = 1; i < discrete_lines.size(); i++){
@@ -510,13 +507,31 @@ void discreteEdge::createGeometry(){
}
}
+MVertex * discreteEdge::getGeometricalVertex (MVertex *v){
+ std::map<MVertex*,MVertex*>::const_iterator it = v2v.find(v);
+ if (it == v2v.end()){
+ printf("----> %p %d %d %g %g %g %d %d\n",v,v->getNum(),tag(),v->x(),v->y(),v->z(),v->onWhat()->tag(),v->onWhat()->dim());
+ Msg::Fatal("fatality %ld %ld %ld",v2v.size(),mesh_vertices.size(),discrete_vertices.size());
+ }
+ return it->second;
+}
+
+void discreteEdge::interpolateInGeometry (MVertex *v, MVertex **v1, MVertex **v2, double &xi) const {
+ double t;
+ if (v->onWhat() != this)Msg::Fatal("%s %d",__FILE__,__LINE__);
+ v->getParameter (0,t);
+ int i = (int) t;
+ MLine *l = discrete_lines [i];
+ *v1 = l->getVertex(0);
+ *v2 = l->getVertex(1);
+ xi = t - i;
+ // Msg::Fatal("fatality...");
+}
+
void discreteEdge::mesh(bool verbose){
#if defined(HAVE_MESH)
- // copy the mesh into geometrical entities
- // FIXME
- // return;
- // createGeometry();
+ if (!CTX::instance()->meshDiscrete) return;
meshGEdge mesher;
mesher(this);
#endif
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 7ece5ad825eda60dbd0a3c9877642798e8c7138b..a44c18559046c0383cb6700518bcc5e82a551f8b 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -12,6 +12,7 @@
class discreteEdge : public GEdge {
protected:
+ std::map<MVertex*,MVertex*> v2v;
std::vector<double> _pars;
std::vector<int> _orientation;
std::map<MVertex*, MLine*> boundv;
@@ -31,6 +32,7 @@ class discreteEdge : public GEdge {
virtual Range<double> parBounds(int) const;
bool getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
+ void interpolateInGeometry (MVertex *v, MVertex **v1, MVertex **v2, double &xi) const;
void parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
std::less<MVertex*> > > &face2Verts,
std::map<GRegion*, std::map<MVertex*, MVertex*,
@@ -43,6 +45,7 @@ class discreteEdge : public GEdge {
virtual void mesh(bool) ;
void writeGEO(FILE *fp);
int minimumDrawSegments() const {return 2*_pars.size();}
+ MVertex * getGeometricalVertex (MVertex *v);
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index d02ebff1e167cc44bc7ac5f04ba53e20517077e2..2e70b3caf90a122b7d5395790f2a0efb2abb2cb4 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -165,8 +165,7 @@ void discreteFace::gatherMeshes() {
void discreteFace::mesh(bool verbose) {
#if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
- return;
- createGeometry();
+ if (!CTX::instance()->meshDiscrete) return;
for (unsigned int i=0;i<_atlas.size();i++)
_atlas[i]->mesh(verbose);
gatherMeshes();
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 11a23f5a106261ee44fb7f2aa857ac5c9c1be186..6b1d0aa9dd244a5f5d422338cd0e03edb744a4e0 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -46,11 +46,12 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
pts[0]->iD, pts[1]->iD, pts[2]->iD, pts[3]->iD);
}
else{
- fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- pts[0]->X, pts[0]->Y, pts[0]->Z,
- pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
- pts[0]->iD, pts[1]->iD, pts[2]->iD);
+ if (pts[0]->iD >= 0 && pts[1]->iD >= 0 && pts[2]->iD >= 0)
+ fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z,
+ pts[1]->X, pts[1]->Y, pts[1]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[0]->iD, pts[1]->iD, pts[2]->iD);
}
}
else{
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 8f04b245277c5fc6ce2583268aeb01be450eb3f2..49c1d120f810197f53fb0a39ff030c0c64f6654c 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -305,11 +305,7 @@ void copyMesh(GEdge *from, GEdge *to, int direction)
void deMeshGEdge::operator() (GEdge *ge)
{
- if(ge->geomType() == GEntity::DiscreteCurve) {
- // FIXME : NOTHING SPECIAL TO DO
- return;
- }
-
+ if(ge->geomType() == GEntity::DiscreteCurve && !CTX::instance()->meshDiscrete)return;
ge->deleteMesh();
ge->meshStatistics.status = GEdge::PENDING;
ge->correspondingVertices.clear();
@@ -560,6 +556,7 @@ void meshGEdge::operator() (GEdge *ge)
const double d = norm(der);
double lc = d/(P1.lc + dlc / dp * (d - P1.p));
GPoint V = ge->point(t);
+ // printf("%d %g\n",NUMP-1,t);
mesh_vertices[NUMP - 1] = new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, lc);
NUMP++;
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 37cf03b57cd193b21470148f65655597960cf8a7..0aceb8e58a472cb24de9a0a3c307fc52f37957f9 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2368,14 +2368,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
void deMeshGFace::operator() (GFace *gf)
{
- if(gf->geomType() == GEntity::DiscreteSurface) {
- discreteFace *df = dynamic_cast<discreteFace*> (gf);
- // copy the mesh and parametrize it to have a "geometry"
- // of the discrete edge
- // FIXME
- return;
- df->createGeometry();
- }
+ if(gf->geomType() == GEntity::DiscreteSurface && !CTX::instance()->meshDiscrete)return;
gf->deleteMesh();
gf->meshStatistics.status = GFace::PENDING;
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
@@ -2383,7 +2376,7 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1; //-1;
+int debugSurface = -100; //-1;
void meshGFace::operator() (GFace *gf, bool print)
{
diff --git a/benchmarks/2d/stator_parametric.geo b/benchmarks/2d/stator_parametric.geo
index c96ba7d3604a1aa00691047b8023de865cc47cc9..2e66d821ee882dbf1ae41c6cc43b2d712351dfb7 100644
--- a/benchmarks/2d/stator_parametric.geo
+++ b/benchmarks/2d/stator_parametric.geo
@@ -1,4 +1,5 @@
-Mesh.Algorithm = 6; // test frontal algorithm
+Mesh.RecombinationAlgorithm=2;
+Mesh.Algorithm = 8; // test frontal algorithm
// Einheiten
mm = 1.e-3;