diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index 4c28872ce9fa63aa9668d7dabbd8d4730cd0235e..b06f5660f59eab26f047973f4d493cd7b13aebc1 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -26,6 +26,8 @@
#include "convexCombinationTerm.h" // #FIXME
#include "qualityMeasuresJacobian.h" // #temporary?
+#include "OptHomRun.h"
+#include "MeshQualityOptimizer.h"
static inline void functionShapes(int p, double Xi[2], double* phi)
{
@@ -286,13 +288,17 @@ discreteDiskFace::discreteDiskFace(int id, GFace *gf, triangulation* diskTriangu
if (!checkOrientationUV()){
Msg::Info("discreteDiskFace:: parametrization is not one-to-one; fixing "
"the discrete system.");
- parametrize(true);
+
+ //parametrize(true);
+ //buildOct(CAD);
+ optimize();
buildOct(CAD);
- }
-
+
+ }
putOnView(true,false);
printParamMesh();
if(!checkOrientationUV()) Msg::Fatal("discreteDiskFace:: failing to fix the discrete system");
+ else Msg::Info("Parameterization done :-)");
}
/*end BUILDER*/
@@ -497,6 +503,87 @@ bool discreteDiskFace::checkOrientationUV()
}
}
+
+void discreteDiskFace::optimize(){
+
+ // parameters for mesh optimization
+ // -- high order
+ OptHomParameters optParams;
+ optParams.dim = 2;
+ optParams.optPassMax = 10;
+ optParams.TMAX = 300;
+ optParams.fixBndNodes = true;
+ optParams.strategy = 1;
+ // -- linear
+ MeshQualOptParameters opt;
+ opt.dim = 2;
+ opt.fixBndNodes = true;
+
+ //creating the "geometry" of the parametrization, and its corresponding mesh
+ // -- generation of parametric nodes
+ std::map<SPoint3,MVertex*> sp2mv;
+ std::vector<MElement*> paramTriangles;
+ for(std::map<MVertex*,SPoint3>::iterator it=coordinates.begin(); it!= coordinates.end(); ++it)
+ sp2mv[it->second] = new MVertex(it->second.x(),it->second.y(),0.);
+ // -- generation of parametric triangles
+ paramTriangles.resize(discrete_triangles.size() - geoTriangulation->fillingHoles.size());
+ for(unsigned int i=0; i<discrete_triangles.size() -geoTriangulation->fillingHoles.size(); i++){
+ discreteDiskFaceTriangle* ct = &_ddft[i];
+ std::vector<MVertex*> mv;
+ mv.resize(ct->tri->getNumVertices());
+ for (int j=0; j<ct->tri->getNumVertices(); j++)
+ mv[j] = sp2mv[ct->p[j]];
+ if(_order==1)
+ paramTriangles[i] = new MTriangle(mv);
+ else
+ paramTriangles[i] = new MTriangle6(mv);
+ }
+ // -- generation of the parametric topology #mark what about the GFace for the GModel ?
+ std::map<int,std::vector<MElement*> > e2e;
+ e2e[0] = paramTriangles;
+ std::vector<int> v;
+ v.push_back(0);
+ std::map<int,std::vector<int> > e2p;
+ e2p[0] = v;
+ GModel* paramDisk = GModel::createGModel(e2e,e2p);
+ discreteVertex dv(paramDisk,NEWPOINT());
+ dv.mesh_vertices.push_back(sp2mv[coordinates[_U0[0]]]);
+ sp2mv[coordinates[_U0[0]]]->setEntity(&dv);
+ paramDisk->add(&dv);
+ discreteEdge de(paramDisk,NEWLINE(),&dv,&dv);
+ paramDisk->add(&de);
+ for(unsigned int i=1; i<_U0.size(); i++){
+ sp2mv[coordinates[_U0[i]]]->setEntity(&de);
+ de.mesh_vertices.push_back(sp2mv[coordinates[_U0[i]]]);
+ de.lines.push_back(new MLine(sp2mv[coordinates[_U0[i-1]]],sp2mv[coordinates[_U0[i]]]));
+ }
+ de.createGeometry();
+
+
+ // optimization
+ if(_order >1)
+ HighOrderMeshOptimizer(paramDisk, optParams);
+ else
+ MeshQualityOptimizer(paramDisk,opt);
+
+ // update the parametrization
+ paramTriangles = e2e[0];
+ for(unsigned int i=0; i< paramTriangles.size(); i++){
+ discreteDiskFaceTriangle* ct = &_ddft[i];
+ MElement* tri = paramTriangles[i];
+ for(int j=0; j<ct->tri->getNumVertices(); j++){
+ MVertex* mv = tri->getVertex(j);
+ SPoint3 p(mv->x(), mv->y(), 0);
+ coordinates[ct->tri->getVertex(j)] = p;
+ ct->p[j] = p;
+ }
+ }
+
+ // cleaning
+ delete paramDisk;
+
+}
+
// (u;v) |-> < (x;y;z); GFace; (u;v) >
GPoint discreteDiskFace::point(const double par1, const double par2) const
{
@@ -1010,7 +1097,7 @@ double triangulation::geodesicDistance () {
addTo (v2t, tri[i]->getVertex(2),tri[i]);
}
- // printf("computing geodesic distance with %d triangles and %d vertcie\n",
+ // printf("computing geodesic distance with %d triangles and %d vertices\n",
// tri.size(),v2t.size());
std::map<MVertex*,double> Fixed;
diff --git a/Geo/discreteDiskFace.h b/Geo/discreteDiskFace.h
index d5e394904502fceadcf6bf2921b9e559d6bff6f0..aa4908def9a7f21d9966566c2a754b5f9ca4a6f2 100644
--- a/Geo/discreteDiskFace.h
+++ b/Geo/discreteDiskFace.h
@@ -22,20 +22,7 @@
#include "meshGFaceOptimize.h"
#include "PView.h"
#include "robustPredicates.h"
-
-inline double tri3Darea(MVertex* mv0, MVertex* mv1, MVertex* mv2)
-{
-
- SVector3 v1(mv1->x()-mv0->x(),mv1->y()-mv0->y(),mv1->z()-mv0->z());
- SVector3 v2(mv2->x()-mv0->x(),mv2->y()-mv0->y(),mv2->z()-mv0->z());
-
- SVector3 n(v1.y()*v2.z()-v2.y()*v1.z(),v2.x()*v1.z()-v1.x()*v2.z(),v1.x()*v2.y()-v2.x()*v1.y());
-
-
- return .5*n.norm();
-
-}
-
+#include "MLine.h"
inline int nodeLocalNum(MElement* e, MVertex* v)
{
@@ -231,6 +218,7 @@ class discreteDiskFace : public GFace {
bool parametrize(bool one2one = false) const;
void putOnView(bool,bool);
bool checkOrientationUV();
+ void optimize();
void printParamMesh();
public:
@@ -274,6 +262,7 @@ class discreteDiskFace : public GFace {
mutable discreteDiskFaceTriangle *_ddft;
mutable Octree *oct;
mutable std::vector<double> _coords;
+
};
#endif
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 776a0a8dddece283f6d778ef2c2ce85afcc97d56..1e5f01ae02b2ebc2cae5be547d051bbda510f939 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -111,7 +111,7 @@ void discreteFace::secondDer(const SPoint2 ¶m,
void discreteFace::createGeometry()
{
checkAndFixOrientation();
- int order = 1;
+ int order = 2;
int nPart = 2;
double eta = 2*3.14/7;
@@ -154,12 +154,12 @@ void discreteFace::createGeometry()
for(unsigned int i=0; i<toParam.size(); i++){
//printf("MAP(%d) : aspect ratio = %12.5E\n",i,toParam[i]->aspectRatio());
- char name[256];
+ //char name[256];
//sprintf(name,"map%d.pos",i);
- toParam[i]->print(name,i);
+ //toParam[i]->print(name,i);
fillHoles(toParam[i]);
//sprintf(name,"mapFilled%d.pos",i);
- toParam[i]->print(name,i);
+ //toParam[i]->print(name,i);
}
for(unsigned int i=0; i<toParam.size(); i++){
std::vector<MElement*> mytri = toParam[i]->tri;
@@ -167,6 +167,7 @@ void discreteFace::createGeometry()
df->replaceEdges(toParam[i]->my_GEdges);
_atlas.push_back(df);
}
+
#endif
}
@@ -215,7 +216,6 @@ void discreteFace::gatherMeshes()
void discreteFace::mesh(bool verbose)
{
-
#if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
if (!CTX::instance()->meshDiscrete) return;
for (unsigned int i=0;i<_atlas.size();i++){