diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 521560a86e2ed4de4f1dee8c3fdf2ce7ac813aa3..0794ef839646d5da0773c5775a9dffbdbb3748a2 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -996,7 +996,7 @@ StringXNumber MeshOptions_Number[] = {
{ 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 , 0 ,
- "Remsh Parametrization (0=harmonic, 1=conformal, 2=rbf harmonic)" },
+ "Remsh Parametrization (0=harmonic, 1=conformal, 2=rbf harmonic, 3=harmonicplane, 4=convex)" },
{ F|O, "RefineSteps" , opt_mesh_refine_steps , 10 ,
"Number of refinement steps in the MeshAdapt-based 2D algorithms" },
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 541ee8436883cb0c6312bb12632576caa10f4284..16ca3dafb5099b334fa8c9c5c01f9aaf98a624c9 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -555,88 +555,6 @@ end:
}
}
-void computeMeanPlane(const std::vector<SPoint3> &points, mean_plane &meanPlane)
-{
- double xm = 0., ym = 0., zm = 0.;
- int ndata = points.size();
- int na = 3;
- for(int i = 0; i < ndata; i++) {
- xm += points[i].x();
- ym += points[i].y();
- zm += points[i].z();
- }
- xm /= (double)ndata;
- ym /= (double)ndata;
- zm /= (double)ndata;
-
- fullMatrix<double> U(ndata, na), V(na, na);
- fullVector<double> sigma(na);
- for(int i = 0; i < ndata; i++) {
- U(i, 0) = points[i].x() - xm;
- U(i, 1) = points[i].y() - ym;
- U(i, 2) = points[i].z() - zm;
- }
- U.svd(V, sigma);
- double res[4], svd[3];
- svd[0] = sigma(0);
- svd[1] = sigma(1);
- svd[2] = sigma(2);
- int min;
- if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
- min = 0;
- else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
- min = 1;
- else
- min = 2;
- res[0] = V(0, min);
- res[1] = V(1, min);
- res[2] = V(2, min);
- norme(res);
-
- double ex[3], t1[3], t2[3];
-
- ex[0] = ex[1] = ex[2] = 0.0;
- if(res[0] == 0.)
- ex[0] = 1.0;
- else if(res[1] == 0.)
- ex[1] = 1.0;
- else
- ex[2] = 1.0;
-
- prodve(res, ex, t1);
- norme(t1);
- prodve(t1, res, t2);
- norme(t2);
-
- res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
-
- for(int i = 0; i < 3; i++)
- meanPlane.plan[0][i] = t1[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[1][i] = t2[i];
- for(int i = 0; i < 3; i++)
- meanPlane.plan[2][i] = res[i];
-
- meanPlane.a = res[0];
- meanPlane.b = res[1];
- meanPlane.c = res[2];
- meanPlane.d = res[3];
-
- meanPlane.x = meanPlane.y = meanPlane.z = 0.;
- if(fabs(meanPlane.a) >= fabs(meanPlane.b) &&
- fabs(meanPlane.a) >= fabs(meanPlane.c) ){
- meanPlane.x = meanPlane.d / meanPlane.a;
- }
- else if(fabs(meanPlane.b) >= fabs(meanPlane.a) &&
- fabs(meanPlane.b) >= fabs(meanPlane.c)){
- meanPlane.y = meanPlane.d / meanPlane.b;
- }
- else{
- meanPlane.z = meanPlane.d / meanPlane.c;
- }
-}
-
-
void GFace::getMeanPlaneData(double VX[3], double VY[3],
double &x, double &y, double &z) const
{
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 9a0eeab3adcfce331adc87a05001808af8a1cfd0..c99a2ed3714cf3033ec16a6436a4ffcbc4ffd680 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -650,26 +650,15 @@ bool GFaceCompound::parametrize() const
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
printStuff(111);
- checkOrientation(0, true);
+ if (_type == MEANPLANE) checkOrientation(0, true);
printStuff(222);
}
- // Multiscale Laplace parametrization
- else if (_mapping == MULTISCALE){
- std::vector<MElement*> _elements;
- for(std::list<GFace*>::const_iterator itt = _compound.begin();
- itt != _compound.end(); ++itt)
- for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
- _elements.push_back((*itt)->triangles[i]);
- multiscaleLaplace multiLaplace(_elements, coordinates);
- }
// Conformal map parametrization
else if (_mapping == CONFORMAL){
Msg::Info("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- // bool hasOverlap = parametrize_conformal_spectral();
std::vector<MVertex *> vert;
bool oriented, hasOverlap;
- //hasOverlap = parametrize_conformal(0,NULL,NULL);
hasOverlap = parametrize_conformal_spectral();
printStuff(11);
if (hasOverlap) oriented = checkOrientation(0);
@@ -781,13 +770,11 @@ void GFaceCompound::getBoundingEdges()
for(std::list<std::list<GEdge*> >::iterator it = _interior_loops.begin();
it != _interior_loops.end(); it++){
double size = getSizeBB(*it);
- printf("size(%d) = %g\n",(*(it->begin()))->tag(),size);
if (size > maxSize) {
_U0 = *it;
maxSize = size;
}
}
- printf("maxSize(%d) = %g\n",tag(),maxSize);
}
}
@@ -980,7 +967,10 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
getBoundingEdges();
_type = UNITCIRCLE;
- //_type = MEANPLANE;
+ if (_mapping == HARMONICPLANE){
+ _mapping = HARMONIC;
+ _type = MEANPLANE;
+ }
nbSplit = 0;
fillTris.clear();
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 6f92c9691dce87ee55bbf6c054824d653425e28d..e050e7a10dd3f3ec3bf60c555d2b9eab53cd276c 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -56,7 +56,7 @@ class GRbf;
class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
- typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4, RBF=5} typeOfMapping;
+ typedef enum {HARMONIC=1,CONFORMAL=2, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
@@ -147,7 +147,7 @@ template<class scalar> class linearSystem;
class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
- typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4, RBF=5} typeOfMapping;
+ typedef enum {HARMONIC=1,CONFORMAL=2, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC, int allowPartition=1)
@@ -165,7 +165,6 @@ class GFaceCompound : public GFace {
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{}
virtual SPoint2 getCoordinates(MVertex *v) const { return SPoint2(); }
void parametrize() const {}
- int getNbSplit() const { return 0; }
};
#endif
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index c5c06cb65d35b8e5e76347ea5488052c1b565cc0..4780eb6743de85d2f1765fb92800ffa140d1f45e 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -13,6 +13,7 @@
#include "gmshRegion.h"
#include "MLine.h"
#include "GModel.h"
+#include "Numeric.h"
GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
@@ -872,8 +873,6 @@ GFace *OCCFactory::addFace(GModel *gm, std::vector<GEdge *> edges,
return gm->_occ_internals->addFaceToModel(gm, TopoDS::Face(aResult));
}
-extern void computeMeanPlane(const std::vector<SPoint3> &points, mean_plane &meanPlane);
-
GFace *OCCFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > wires)
{
@@ -892,7 +891,7 @@ GFace *OCCFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
}
mean_plane meanPlane;
- computeMeanPlane(points, meanPlane);
+ computeMeanPlaneSimple(points, meanPlane);
gp_Pln aPlane (meanPlane.a,meanPlane.b,meanPlane.c,meanPlane.d);
BRepBuilderAPI_MakeFace aGenerator (aPlane);
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index c8453df8833ec40009eea055c1e883a6cd620264..985df492a52dbad0009ce89268807a023ab7044f 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -190,9 +190,12 @@ int GModel::importGEOInternals()
}
}
int param = CTX::instance()->mesh.remeshParam;
- GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
- if (param == 1) typ = GFaceCompound::CONFORMAL;
- if (param == 2) typ = GFaceCompound::RBF;
+ GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
+ if (param == 1) typ = GFaceCompound::CONFORMAL;
+ if (param == 2) typ = GFaceCompound::RBF;
+ if (param == 3) typ = GFaceCompound::HARMONICPLANE;
+ if (param == 4) typ = GFaceCompound::CONVEXCOMBINATION;
+
int algo = CTX::instance()->mesh.remeshAlgo;
f = new GFaceCompound(this, s->Num, comp, U0, typ, algo);