diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 363ea7ad1ac917264d166513ca557682ebcf1f1d..227b996a0d9a979fe2f81e5d2c0d3193c9cc9c4c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -998,7 +998,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, 3=harmonicplane, 4=convex)" },
+ "Remeshing using discrete parametrization (0=harmonic_circle, 1=conformal, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square" },
{ F|O, "RefineSteps" , opt_mesh_refine_steps , 10 ,
"Number of refinement steps in the MeshAdapt-based 2D algorithms" },
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 4338608e1554f64cd3c2be56df3aa4269521b364..7b37ae4ed62e7759b0fa184c340710595ec7a3b6 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -644,11 +644,12 @@ bool GFaceCompound::parametrize() const
if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
// Convex parametrization
- if (_mapping == CONVEXCOMBINATION){
+ if (_mapping == CONVEX){
Msg::Info("Parametrizing surface %d with 'convex map'", tag());
fillNeumannBCS();
- parametrize(ITERU,CONVEXCOMBINATION);
- parametrize(ITERV,CONVEXCOMBINATION);
+ parametrize(ITERU,CONVEX);
+ parametrize(ITERV,CONVEX);
+ printStuff(111);
}
// Laplace parametrization
else if (_mapping == HARMONIC){
@@ -707,8 +708,8 @@ bool GFaceCompound::parametrize() const
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
- parametrize(ITERU,CONVEXCOMBINATION);
- parametrize(ITERV,CONVEXCOMBINATION);
+ parametrize(ITERU,CONVEX);
+ parametrize(ITERV,CONVEX);
checkOrientation(0);
buildOct();
}
@@ -895,11 +896,11 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
- std::list<GEdge*> &U0, typeOfMapping mpg,
+ std::list<GEdge*> &U0, typeOfCompound toc,
int allowPartition,
linearSystem<double>* lsys)
: GFace(m, tag), _compound(compound), oct(0), _U0(U0),
- _mapping(mpg), _allowPartition(allowPartition), _lsys(lsys)
+ _toc(toc), _allowPartition(allowPartition), _lsys(lsys)
{
ONE = new simpleFunction<double>(1.0);
MONE = new simpleFunction<double>(-1.0);
@@ -912,9 +913,16 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
}
getBoundingEdges();
+
+ _mapping = HARMONIC;
_type = UNITCIRCLE;
- if (_mapping == HARMONICPLANE){
- _mapping = HARMONIC;
+ if (toc == HARMONIC_CIRCLE) _mapping = HARMONIC;
+ else if(toc == CONFORMAL) _mapping = CONFORMAL;
+ else if(toc == RBF) _mapping = RBF;
+ else if (toc == HARMONIC_PLANE) _type = MEANPLANE;
+ else if (toc == CONVEX_CIRCLE) _mapping = CONVEX;
+ else if (toc == CONVEX_PLANE){
+ _mapping = CONVEX;
_type = MEANPLANE;
}
@@ -922,14 +930,14 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
fillTris.clear();
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
- std::list<GEdge*> &U0, std::list<GEdge*> &V0,
- std::list<GEdge*> &U1, std::list<GEdge*> &V1,
- typeOfMapping mpg,
- int allowPartition,
- linearSystem<double>* lsys)
+ std::list<GEdge*> &U0, std::list<GEdge*> &V0,
+ std::list<GEdge*> &U1, std::list<GEdge*> &V1,
+ typeOfCompound toc,
+ int allowPartition,
+ linearSystem<double>* lsys)
: GFace(m, tag), _compound(compound), oct(0),
_U0(U0), _V0(V0), _U1(U1), _V1(V1),
- _mapping(mpg), _allowPartition(allowPartition), _lsys(lsys)
+ _toc(toc), _allowPartition(allowPartition), _lsys(lsys)
{
ONE = new simpleFunction<double>(1.0);
MONE = new simpleFunction<double>(-1.0);
@@ -942,12 +950,19 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
}
getBoundingEdges();
+
+ _mapping = HARMONIC;
_type = UNITCIRCLE;
- if (_mapping == HARMONICPLANE){
- _mapping = HARMONIC;
+ if (toc == HARMONIC_CIRCLE) _mapping = HARMONIC;
+ else if(toc == CONFORMAL) _mapping = CONFORMAL;
+ else if(toc == RBF) _mapping = RBF;
+ else if (toc == HARMONIC_PLANE) _type = MEANPLANE;
+ else if (toc == CONVEX_CIRCLE) _mapping = CONVEX;
+ else if (toc == CONVEX_PLANE){
+ _mapping = CONVEX;
_type = MEANPLANE;
}
- else if(_mapping == HARMONIC && _U0.size() && _V0.size() && _U1.size() && _V1.size()) {
+ else if(toc == HARMONIC_SQUARE && _U0.size() && _V0.size() && _U1.size() && _V1.size()) {
_type = SQUARE;
}
@@ -1053,20 +1068,34 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
+
linearSystem<double> *lsys = 0;
if (_lsys) lsys = _lsys;
else{
-#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
+ if (tom==HARMONIC){
+#if defined(HAVE_TAUCS)
+ lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
lsys = new linearSystemPETSc<double>;
#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
lsysb->setGmres(1);
lsys = lsysb;
-#elif defined(HAVE_TAUCS)
- lsys = new linearSystemCSRTaucs<double>;
#else
lsys = new linearSystemFull<double>;
#endif
+ }
+ else if (tom==CONVEX){
+#if defined(HAVE_PETSC)
+ lsys = new linearSystemPETSc<double>;
+#elif (defined(HAVE_GMM)
+ linearSystemGmm<double> *lsysb = new linearSystemGmm<double>;
+ lsysb->setGmres(1);
+ lsys = lsysb;
+#else
+ lsys = new linearSystemFull<double>;
+#endif
+ }
}
dofManager<double> myAssembler(lsys);
@@ -1154,9 +1183,9 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
femTerm<double> *mapping;
if (tom == HARMONIC)
mapping = new laplaceTerm(0, 1, ONE);
- else
+ else
mapping = new convexCombinationTerm(0, 1, ONE);
-
+
it = _compound.begin();
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 7684deb6e12b33b937e91503b241b1333dc597c3..17a093b7f9cbf94c9d9e51537ae6d98282b10c78 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -26,17 +26,14 @@ template <class scalar> class simpleFunction;
/*
A GFaceCompound is a model face that is the compound of model faces.
-The user selects 4 corners to the compound, which are mmapped to
-points (0,0), (1,0), (1,1) and (0,1) of the compound parametric plane.
-
It is assumed that all the faces of the compound have been meshed
-first. We use this discretization to solve 2 elliptic problems on the
-compound. Those 2 problems enables to compute the parametric
+first. We use this discretization to solve elliptic problems on the
+compound. Those problems enable to compute the parametric
coordinates of the mesh points. The parametrization of the compound
consist in a triangulation in the (u,v) space with parameter values at
nodes.
-The compound can therefore be re-meshed using any surface mesh
+The compound can therefore be (re)-meshed using any surface mesh
generator of gmsh!
*/
@@ -56,7 +53,9 @@ class GRbf;
class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
- typedef enum {HARMONIC=1,CONFORMAL=2, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
+ typedef enum {HARMONIC_CIRCLE=0, CONFORMAL_SPECTRAL=1, RADIAL_BASIS=2, HARMONIC_PLANE=3,
+ CONVEX_CIRCLE=4,CONVEX_PLANE=5, HARMONIC_SQUARE=6} typeOfCompound;
+ typedef enum {HARMONIC=0,CONFORMAL=1, RBF=2, CONVEX=3} typeOfMapping;
typedef enum {UNITCIRCLE, MEANPLANE, SQUARE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
@@ -108,20 +107,20 @@ class GFaceCompound : public GFace {
public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
- std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC,
+ std::list<GEdge*> &U0, typeOfCompound typ = HARMONIC_CIRCLE,
int allowPartition=1,
linearSystem<double>* lsys =0);
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
std::list<GEdge*> &U1, std::list<GEdge*> &V1,
- typeOfMapping typ = HARMONIC,
+ typeOfCompound typ = HARMONIC_CIRCLE,
int allowPartition=1,
linearSystem<double>* lsys =0);
virtual ~GFaceCompound();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
virtual GPoint point(double par1, double par2) const;
- typeOfMapping getTypeOfMapping() { return _mapping;}
+ typeOfCompound getTypeOfCompound() { return _toc;}
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
@@ -143,18 +142,22 @@ class GFaceCompound : public GFace {
int allowPartition() const{ return _allowPartition; }
void setType(typeOfIsomorphism type){ _type=type;}
private:
- typeOfIsomorphism _type;
+ mutable typeOfCompound _toc;
mutable typeOfMapping _mapping;
+ typeOfIsomorphism _type;
int _allowPartition;
};
#else
+//define empty class ifndef HAVE_SOLVER
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, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
+ typedef enum {HARMONIC_CIRCLE=0, CONFORMAL_SPECTRAL=1, RADIAL_BASIS=2, HARMONIC_PLANE=3,
+ CONVEX_CIRCLE=4,CONVEX_PLANE=5, HARMONIC_SQUARE=6} typeOfCompound;
+ typedef enum {HARMONIC=0,CONFORMAL=1, RBF=2, CONVEX=3} typeOfMapping;
typedef enum {UNITCIRCLE, MEANPLANE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC,
@@ -167,7 +170,7 @@ class GFaceCompound : public GFace {
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &V0,
std::list<GEdge*> &U1, std::list<GEdge*> &V1,
- typeOfMapping typ = HARMONIC,
+ typeOfCompound typ = HARMONIC_CIRCLE,
int allowPartition=1,
linearSystem<double>* lsys =0)
: GFace(m, tag)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index ce1a2e9d45809943d8ed1b825562107049e09bac..3d443e936e809f8212ceb000875da5c1fd6cabf6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2066,7 +2066,7 @@ void GModel::save(std::string fileName)
GModel::setCurrent(temp);
}
-GFaceCompound* GModel::addCompoundFace(std::vector<GFace*> faces, int typeP, int typeS)
+GFaceCompound* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int typeS)
{
int num = getMaxElementaryNumber(2) + 1;
@@ -2074,9 +2074,13 @@ GFaceCompound* GModel::addCompoundFace(std::vector<GFace*> faces, int typeP, int
std::list<GFace*> comp(faces.begin(), faces.end());
std::list<GEdge*> U0;
- GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
- if (typeP == 1) typ = GFaceCompound::CONFORMAL;
- if (typeP == 2) typ = GFaceCompound::RBF;
+ GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
+ if (param == 1) typ = GFaceCompound::CONFORMAL_SPECTRAL;
+ if (param == 2) typ = GFaceCompound::RADIAL_BASIS;
+ if (param == 3) typ = GFaceCompound::HARMONIC_PLANE;
+ if (param == 4) typ = GFaceCompound::CONVEX_CIRCLE;
+ if (param == 5) typ = GFaceCompound::CONVEX_PLANE;
+ if (param == 6) typ = GFaceCompound::HARMONIC_SQUARE;
GFaceCompound *gfc = new GFaceCompound(this, num, comp, U0, typ, typeS);
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index f3b8fab3fad80280068ac500d231a85f92ec8113..b78cd64679b9ff41378d177c4f76a1144e94bccc 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -191,11 +191,14 @@ 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;
- if (param == 3) typ = GFaceCompound::HARMONICPLANE;
- if (param == 4) typ = GFaceCompound::CONVEXCOMBINATION;
+
+ GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_CIRCLE;
+ if (param == 1) typ = GFaceCompound::CONFORMAL_SPECTRAL;
+ if (param == 2) typ = GFaceCompound::RADIAL_BASIS;
+ if (param == 3) typ = GFaceCompound::HARMONIC_PLANE;
+ if (param == 4) typ = GFaceCompound::CONVEX_CIRCLE;
+ if (param == 5) typ = GFaceCompound::CONVEX_PLANE;
+ if (param == 6) typ = GFaceCompound::HARMONIC_SQUARE;
int algo = CTX::instance()->mesh.remeshAlgo;
f = new GFaceCompound(this, s->Num, comp, b[0], b[1], b[2], b[3], typ, algo);
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 67f524934745798842800d8f25cb631b5c60fbb2..ab3b5ebf246294d4c4953a80d9b1913ccdf5b6fa 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -634,8 +634,8 @@ void Centerline::createSplitCompounds(){
int num_gfc = NF+i+1;
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONICPLANE;
- //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL;
+ GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE;
+ //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL_SPECTRAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 22f503a820749f56a2372258eb82ed0cc556d637..d1a9771d51ff1c455e60436f7b6145ee849cf7b0 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1952,7 +1952,7 @@ void partitionAndRemesh(GFaceCompound *gf)
num_gfc, pf->tag());
GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, U0,
- gf->getTypeOfMapping());
+ gf->getTypeOfCompound());
gfc->meshAttributes.recombine = gf->meshAttributes.recombine;
gf->model()->add(gfc);
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index d999a2182c1e894c1ef23ba213775c8127d322ff..cb64fb86c6f8039930e5aff619e465b872446c07 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -42,31 +42,32 @@ class convexCombinationTerm : public femTerm<double> {
MElement *e = se->getMeshElement();
m.setAll(0.);
const int nbSF = e->getNumShapeFunctions();
- const double _diff = 1.0;
for (int j = 0; j < nbSF; j++){
+ MVertex *vj = e->getShapeFunctionNode(j);
+ double diag = 0.0;
for (int k = 0; k < nbSF; k++){
- m(j, k) = -1. * _diff;
+ MVertex *vk = e->getShapeFunctionNode(k);
+ MVertex *vl = NULL;
+ SVector3 a, b;
+ double tan = 0.0;
+ if ( k!=j && 3-j-k >= 0 && 3-j-k <=2) {
+ vl = e->getShapeFunctionNode(3-j-k);
+ a = SVector3(vk->x()-vj->x(),vk->y()-vj->y(),vk->z()-vj->z());
+ b = SVector3(vl->x()-vj->x(),vl->y()-vj->y(),vl->z()-vj->z());
+ double angle = myacos(dot(a,b)/(norm(a)*norm(b)));
+ tan = sin(angle*0.5)/cos(angle*0.5);
+ }
+ double length = vj->distance(vk);
+ m(j, k) = -tan/length;
+ //m(j, k) = -1;
}
- m(j, j) = (nbSF - 1) * _diff;
+ for (int k = 0; k < nbSF; k++){
+ if (k != j) diag += (-m(j,k));
+ }
+ m(j, j) = diag;
}
}
-//diag matrix
-/* virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const */
-/* { */
-
-/* MElement *e = se->getMeshElement(); */
-/* m.setAll(0.); */
-/* const int nbSF = e->getNumShapeFunctions(); */
-/* for (int j = 0; j < nbSF; j++){ */
-/* for (int k = 0; k < nbSF; k++) m(j,k) = 0.0; */
-/* MVertex *v = e->getShapeFunctionNode(j); */
-/* if( v->onWhat()->dim() < 2) m(j,j) = (nbSF - 1) ; */
-/* else m(j,j) = 0.0; */
-/* } */
-/* } */
-
-
};
#endif