diff --git a/Geo/GFace.h b/Geo/GFace.h
index 8ca6d9ab2f36c13b1b561cc8bf76ff8b6cf908d0..ae98641424a071d14d0b4d63e02b916cf0248e99 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -250,8 +250,6 @@ class GFace : public GEntity
// all diagonals of the triangulation are left (-1), right (1) or
// alternated (0)
int transfiniteArrangement;
- // type of mapping is cad (0), harmonic (1) or conformal (-1)
- int typeOfMapping;
// do we smooth (transfinite) mesh? (<0 to use default smoothing)
int transfiniteSmoothing;
// the extrusion parameters (if any)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index d0de4112fdfe2c67aa29b173af8f347b15d5750e..dc96ff257e3e9c102658416bef34e2aa8ba58818 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -475,7 +475,8 @@ bool GFaceCompound::parametrize() const
if(trivial()) return paramOK;
coordinates.clear();
-
+ computeNormals();
+
if(allNodes.empty()) buildAllNodes();
@@ -517,7 +518,7 @@ bool GFaceCompound::parametrize() const
printStuff();
if (!checkOrientation(0)){
- Msg::Info("*** Parametrization switched to convex combination map");
+ Msg::Info("--- Parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
@@ -529,7 +530,7 @@ bool GFaceCompound::parametrize() const
- computeNormals();
+
if (checkAspectRatio() > AR_MAX){
Msg::Warning("Geometrical aspect ratio too high");
@@ -1305,15 +1306,14 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
return lt->gf->curvatureMax(pv);
}
else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
- //printf("!!!! compute here curvatureMax \n");
double curv= 0.;
- curv = curvature(lt->tri,U,V);
+ curv = locCurvature(lt->tri,U,V);
return curv;
}
return 0.;
}
-double GFaceCompound::curvature(MTriangle *t, double u, double v) const
+double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
{
SVector3 n1 = _normals[t->getVertex(0)];
@@ -1323,9 +1323,8 @@ double GFaceCompound::curvature(MTriangle *t, double u, double v) const
n1.y(), n2.y(), n3.y(),
n1.z(), n2.z(), n3.z()};
- return fabs(t->interpolateDiv(val, u, v, 0.0));
+ return fabs(t->interpolateDiv(val, u, v, 3));
- return 0.;
}
SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p) const
@@ -1358,10 +1357,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
gp.setNoSuccess();
return gp;
}
- // if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
- // SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
- // return lt->gf->point(pv.x(),pv.y());
- // }
const bool LINEARMESH = true; //false
@@ -1785,7 +1780,7 @@ double GFaceCompound::checkAspectRatio() const
void GFaceCompound::coherenceNormals()
{
- Msg::Info("Coherence Normals ");
+ Msg::Info("Re-orient all triangles (face normals) coherently");
std::map<MEdge, std::set<MTriangle*>, Less_Edge > edge2tris;
for(unsigned int i = 0; i < triangles.size(); i++){
@@ -1876,6 +1871,9 @@ int GFaceCompound::genusGeom() const
void GFaceCompound::printStuff() const
{
+
+ if( !CTX::instance()->mesh.saveAll) return;
+
std::list<GFace*>::const_iterator it = _compound.begin();
char name0[256], name1[256], name2[256], name3[256];
@@ -1888,13 +1886,6 @@ void GFaceCompound::printStuff() const
sprintf(name5, "XYZV-%d.pos", (*it)->tag());
sprintf(name6, "XYZC-%d.pos", (*it)->tag());
- // sprintf(name1, "UVX.pos");
-// sprintf(name2, "UVY.pos");
-// sprintf(name3, "UVZ.pos");
-// sprintf(name4, "XYZU.pos");
-// sprintf(name5, "XYZV.pos");
-// sprintf(name6, "XYZC.pos");
-
FILE * uva = fopen(name0,"w");
FILE * uvx = fopen(name1,"w");
FILE * uvy = fopen(name2,"w");
@@ -1920,13 +1911,6 @@ void GFaceCompound::printStuff() const
coordinates.find(t->getVertex(1));
std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
-// fprintf(xyzu,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
-// t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
-// t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
-// t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
-// (35.*it0->second.x()-t->getVertex(0)->x()), -t->getVertex(0)->y(), (35.*it0->second.y()-t->getVertex(0)->z()),
-// (35.*it1->second.x()-t->getVertex(1)->x()), -t->getVertex(1)->y(), (35.*it1->second.y()-t->getVertex(1)->z()),
-// (35.*it2->second.x()-t->getVertex(2)->x()), -t->getVertex(2)->y(), (35.*it2->second.y()-t->getVertex(2)->z()));
fprintf(xyzv,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
@@ -1937,16 +1921,14 @@ void GFaceCompound::printStuff() const
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
it0->second.x(),it1->second.x(),it2->second.x());
- //double K1 = curvature(t,it0->second.x(),it0->second.y());
- //double K2 = curvature(t,it1->second.x(),it1->second.y());
- //double K3 = curvature(t,it2->second.x(),it2->second.y());
- // const double K = fabs(curvature (t));
+ double K1 = locCurvature(t,it0->second.x(),it0->second.y());
+ double K2 = locCurvature(t,it1->second.x(),it1->second.y());
+ double K3 = locCurvature(t,it2->second.x(),it2->second.y());
fprintf(xyzc,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- it0->second.z(),it1->second.z(),it2->second.z());
- //K1, K2, K3);
+ K1, K2, K3);
double p0[3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
double p1[3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 023db4e6dbd1111db178428f832b90f582d3b232..286c81dd15a569cdd125f2d122c46872e9a2cf07 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -85,7 +85,7 @@ class GFaceCompound : public GFace {
void computeALoop(std::set<GEdge*> & _unique, std::list<GEdge*> &);
void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
double &_u, double &_v) const;
- virtual double curvature(MTriangle *t, double u, double v) const;
+ virtual double locCurvature(MTriangle *t, double u, double v) const;
void printStuff() const;
bool trivial() const ;
linearSystem <double> *_lsys;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 619415efe615bdc25f0f9a396daec22200d94fe5..9ff2c2e4fbea20bb6d5e88404e4b1c007cf6f789 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -73,7 +73,6 @@ SPoint2 discreteFace::parFromPoint(const SPoint3 &p) const
{
if (getCompound()){
- printf("par from point in GFaceCompound %d \n", getCompound()->tag());
return getCompound()->parFromPoint(p);
}
else{
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 65f17842934c68261ac8c34d48d2b25f1eb8c475..0087e967a32fb32152b239b9443804bdb518b584 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -39,13 +39,10 @@ static double max_surf_curvature(const GEdge *ge, double u)
double val = 0;
std::list<GFace *> faces = ge->faces();
std::list<GFace *>::iterator it = faces.begin();
- //printf("for gedge=%d facesnb=%d \n", ge->tag(), faces.size());
while(it != faces.end()){
if ((*it)->geomType() != GEntity::CompoundSurface){
- printf("for gedge %d face=%d \n",ge->tag(), (*it)->tag());
SPoint2 par = ge->reparamOnFace((*it), u, 1);
double cc = (*it)->curvature(par);
- printf("for face compute curvature =%g \n", cc);
val = std::max(cc, val);
}
++it;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index c13d5071eed22864d795b4dc8657c9c31b9e8007..e7f86b426d4898c97d603fb07722b13faf9d9726 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1356,12 +1356,12 @@ void partitionAndRemesh(GFaceCompound *gf)
gf->model()->createTopologyFromFaces(pFaces);
- Msg::Info("*** Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF );
- CreateOutputFile("multiscalePARTS.msh", CTX::instance()->mesh.format);
+ Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF );
+ gf->model()->writeMSH("multiscalePARTS.msh", 2.0, false, true);
//Remesh new faces (Compound Lines and Compound Surfaces)
//-----------------------------------------------------
- Msg::Info("*** Parametrize Compounds:");
+ Msg::Info("--- Parametrize Compounds:");
//Parametrize Compound Lines
int NE = gf->model()->maxEdgeNum() - nume + 1;
@@ -1395,11 +1395,10 @@ void partitionAndRemesh(GFaceCompound *gf)
Msg::Info("*** Mesh Compounds:");
- //Mesh 1D and 2D
for (int i=0; i < NE; i++){
GEdge *gec = gf->model()->getEdgeByTag(nume + NE + i);
meshGEdge mge;
- mge(gec);//meshing 1D
+ mge(gec);
}
Msg::Info("*** Starting Mesh of surface %d ...", gf->tag());
@@ -1407,7 +1406,7 @@ void partitionAndRemesh(GFaceCompound *gf)
for (int i=0; i < NF; i++){
GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;
- mgf(gfc);//meshing 2D
+ mgf(gfc);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index dcc7a010e541616f25db671b39af450567cee468..6ae61efc40fc931dcdaaa4f7edfeaa80bc3f2572 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -343,7 +343,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
int genus, AR, NB;
getGenusAndRatio(regions[i], genus, AR, NB);
- printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
+ //printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
@@ -356,7 +356,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfP
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
}
- else if (genus == 0 && AR > 4 || genus == 0 && NB > 1){
+ else if (genus == 0 && AR > 3 || genus == 0 && NB > 1){
int nbParts = 2;
Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
nextLevel->recur,nextLevel->region, AR, NB, nbParts);
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 9e69c27f7eb85c84aea3f9b3780818352bae88d6..5852f1d9bff9f2d699b1ced4bfc612c84f890bb3 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -590,7 +590,8 @@ static void printLevel ( const char* fn,
std::map<MVertex*,SPoint2> *coordinates,
double version){
-
+ if( !CTX::instance()->mesh.saveAll) return;
+
std::set<MVertex*> vs;
for (int i=0;i<elements.size();i++)
for (int j=0;j<elements[i]->getNumVertices();j++)
@@ -1073,9 +1074,9 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- printLevel ("multiscale-all.msh",elements, 0,2.0);
- printLevel ("multiscale-left.msh",left,0,2.0);
- printLevel ("multiscale-right.msh",right,0,2.0);
+ //printLevel ("multiscale-all.msh",elements, 0,2.0);
+ //printLevel ("multiscale-left.msh",left,0,2.0);
+ //printLevel ("multiscale-right.msh",right,0,2.0);
}
diff --git a/benchmarks/extrude/sphere_boundary_layer.geo b/benchmarks/extrude/sphere_boundary_layer.geo
index 6e177750aab1fdf9efe312946f7aaf8bdea17740..02c20ae9ebc022f7912716da415f6ea38fc36035 100644
--- a/benchmarks/extrude/sphere_boundary_layer.geo
+++ b/benchmarks/extrude/sphere_boundary_layer.geo
@@ -38,7 +38,7 @@ Line Loop(27) = {-4,12,-6};
Ruled Surface(28) = {27};
Extrude {
- Surface{14:28:2}; Layers{10, 0.1}; // Recombine;
+ Surface{14:28:2}; Layers{10, 0.2}; // Recombine;
}
diff --git a/benchmarks/stl/artery.geo b/benchmarks/stl/artery.geo
index 92f345193d8bb55a7880c621b12f8cbc916eb5e9..1462ff23849e07caf1b41148e680a55df66681bf 100644
--- a/benchmarks/stl/artery.geo
+++ b/benchmarks/stl/artery.geo
@@ -1,10 +1,9 @@
-Mesh.CharacteristicLengthFactor=0.1;
-Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)
+Mesh.CharacteristicLengthFactor=0.05;
//merge the stl triangulation
Merge "artery.stl";
CreateTopology;
-Compound Surface(100)={1} Harmonic;
+Compound Surface(100)={1} Conformal;
Physical Surface(101)={100};
\ No newline at end of file
diff --git a/benchmarks/stl/implant.geo b/benchmarks/stl/implant.geo
new file mode 100644
index 0000000000000000000000000000000000000000..fe6d02377a5e135bbbc89f5b94e09891a7b7b1d7
--- /dev/null
+++ b/benchmarks/stl/implant.geo
@@ -0,0 +1,9 @@
+Mesh.CharacteristicLengthFactor=0.2;
+
+Merge "implant2.stl";
+CreateTopology;
+
+Compound Surface(100)={1};
+
+Physical Surface (501)={100};
+
diff --git a/benchmarks/stl/implant_CLASS_GEO.geo b/benchmarks/stl/implant_CLASS_GEO.geo
deleted file mode 100644
index fcfa5a5a1766c3fb188f9f4cfe1dc5a8331f8229..0000000000000000000000000000000000000000
--- a/benchmarks/stl/implant_CLASS_GEO.geo
+++ /dev/null
@@ -1,20 +0,0 @@
-Mesh.CharacteristicLengthFactor=0.1;
-
-Merge "implant_CLASS.msh";
-CreateTopology;
-
-Compound Line(20)={3};
-Compound Surface(100)={2} Boundary {{}};
-Compound Surface(200)={3} Boundary {{}};
-
-//Compound Line(20)={6};
-//Compound Line(30)={7};
-//Compound Line(40)={8};
-
-//Compound Surface(100)={2} Boundary {{}};
-//Compound Surface(200)={3} Boundary {{}};
-//Compound Surface(300)={4} Boundary {{}};
-//Compound Surface(400)={5} Boundary {{}};
-
-
-
diff --git a/benchmarks/stl/pelvis-BAD.stl b/benchmarks/stl/pelvis-BAD.stl
new file mode 100644
index 0000000000000000000000000000000000000000..4922ea70d83e6b4a07a6bc2bf383afc42eeab330
Binary files /dev/null and b/benchmarks/stl/pelvis-BAD.stl differ
diff --git a/benchmarks/stl/pelvis.geo b/benchmarks/stl/pelvis.geo
index ed4c625bc814c51f6be5ad10d911228acc36c422..0a3ac76140b1dd521923479b98b070705d406e30 100644
--- a/benchmarks/stl/pelvis.geo
+++ b/benchmarks/stl/pelvis.geo
@@ -1,10 +1,10 @@
-//Mesh.CharacteristicLengthFactor=0.133;
+
+Mesh.CharacteristicLengthFactor=0.1;
Mesh.Algorithm3D = 4; //Frontal (4) Delaunay(1)
-Merge "pelvisSMOOTH.stl";
-CreateTopology;
+Merge "pelvis.stl";
-Compound Surface(200)={1} Conformal;
+Compound Surface(200)={1};
Surface Loop(300)={200};
Volume(301)={300};
diff --git a/benchmarks/stl/pelvis.stl b/benchmarks/stl/pelvis.stl
index 4922ea70d83e6b4a07a6bc2bf383afc42eeab330..15e69a5510d7c509f58ce0ea57ca2a875a5d8e3e 100644
Binary files a/benchmarks/stl/pelvis.stl and b/benchmarks/stl/pelvis.stl differ
diff --git a/benchmarks/stl/pelvisSMOOTH.stl b/benchmarks/stl/pelvisSMOOTH.stl
index db4770d05ccb88d4db1b52ecda4fb56f5648aed6..15e69a5510d7c509f58ce0ea57ca2a875a5d8e3e 100644
Binary files a/benchmarks/stl/pelvisSMOOTH.stl and b/benchmarks/stl/pelvisSMOOTH.stl differ
diff --git a/benchmarks/stl/reparam.geo b/benchmarks/stl/reparam.geo
deleted file mode 100644
index 2319536c39ab262af32ce4b8a555c18019d600ef..0000000000000000000000000000000000000000
--- a/benchmarks/stl/reparam.geo
+++ /dev/null
@@ -1,12 +0,0 @@
-
-Mesh.CharacteristicLengthFactor=0.5;
-Geometry.Lines = 0;
-
-Merge "reparam_input.msh";
-CreateTopology;
-
-Compound Line(1000)={1};
-Compound Line(2000)={2};
-Compound Line(3000)={3};
-Compound Line(4000)={4};
-Compound Surface(-1000)={5};
diff --git a/contrib/Netgen/libsrc/meshing/meshclass.cpp b/contrib/Netgen/libsrc/meshing/meshclass.cpp
index 7dbec43e959bc56cb5fc3052aaa3129231b6d58a..c47ffa9d9c7002f50d7b3ba836d49150e68f02d4 100644
--- a/contrib/Netgen/libsrc/meshing/meshclass.cpp
+++ b/contrib/Netgen/libsrc/meshing/meshclass.cpp
@@ -3334,16 +3334,16 @@ namespace netgen
if (IntersectTriangleTriangle (&trip1[0], &trip2[0]))
{
overlap = 1;
- PrintWarning ("Intersecting elements "
+ PrintError ("Intersecting elements "
,i, " and ", inters.Get(j));
(*testout) << "Intersecting: " << endl;
(*testout) << "openelement " << i << " with open element " << inters.Get(j) << endl;
- cout << "el1 = " << tri << endl;
- cout << "el2 = " << tri2 << endl;
- cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl;
- cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl;
+// cout << "el1 = " << tri << endl;
+// cout << "el2 = " << tri2 << endl;
+// cout << "layer1 = " << (*this)[tri[0]].GetLayer() << endl;
+// cout << "layer2 = " << (*this)[tri2[0]].GetLayer() << endl;
for (k = 1; k <= 3; k++)
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 0ee52fa64bf5599f3bb6c52ee252c031d5e997b2..48880602553c0123c459b8c8894d1a381d1cb91b 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -2226,6 +2226,7 @@ elementary surfaces should be oriented consistently (using negative
identification numbers to specify reverse orientation). (Surface loops are
used to create volumes: see @ref{Volumes}.)
+
@item Compound Surface ( @var{expression} ) = @{ @var{expression-list} @} < Boundary @{ @{ @var{expression-list} @}, @{ @var{expression-list} @}, @{ @var{expression-list} @}, @{ @var{expression-list} @} @} > < Harmonic | Conformal | Harmonic_NoSplit | Conformal_NoSplit >;
Creates a compound surface from several elementary surfaces. When
meshed, a compound surface will be reparametrized as a single surface,
@@ -2233,14 +2234,12 @@ whose mesh can thus cross internal boundaries. Compound surfaces are
mostly useful for remeshing discrete models; see ``J.-F. Remacle,
C. Geuzaine, G. Compere and E. Marchandise, @emph{High Quality Surface
Remeshing Using Harmonic Maps}, International Journal for Numerical
-Methods in Engineering, 2009'' for details. The @var{expression} inside
-the parentheses is the compound surface's identification number; the
+Methods in Engineering, 2009'' for details as well as the wiki for more
+examples. The @var{expression} inside the parentheses is the compound surface's identification number; the
mandatory @var{expression-list} on the right hand side contains the
identification number of the elementary surfaces that should be
reparametrized as a single surface.
-@c TODO document the optional arguments.
-
@item Physical Surface ( @var{expression} | @var{char-expression} ) = @{ @var{expression-list} @};
Creates a physical surface. The @var{expression} inside the parentheses
is the physical surface's identification number (if a