diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 78b3d77db3a7ef4b1d9d1e899a46410f9e05a2be..ef0fe4bcdd36ae056414a5ed5e02eb2ae77436f2 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -829,7 +829,7 @@ bool GFaceCompound::parametrize() const
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
- printStuff(111);
+ printStuff(111);
if (_type == MEANPLANE) checkOrientation(0, true);
printStuff(222);
}
@@ -2074,9 +2074,6 @@ void GFaceCompound::buildOct() const
bool GFaceCompound::checkTopology() const
{
if (_mapping == RBF) return true;
-
- //fixme tristan
- //return true;
bool correctTopo = true;
if(allNodes.empty()) buildAllNodes();
@@ -2498,9 +2495,9 @@ void GFaceCompound::printStuff(int iNewton) const
// useful for mesh generators ----------------------------------------
// Intersection of a circle and a plane
-
GPoint GFaceCompound::intersectionWithCircle (const SVector3 &n1, const SVector3 &n2, const SVector3 &p, const double &d, double uv[2]) const
{
+
SVector3 n = crossprod(n1,n2);
n.normalize();
for (int i=-1;i<nbT;i++){
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index b13a884e1ae79f11c131892a9332a36b9cade51a..aaaf4ec24ee406f73d5a85ac82acd6a6519f54b1 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -15,6 +15,7 @@
#include "laplaceTerm.h"
#include "convexCombinationTerm.h"
#include "linearSystemCSR.h"
+#include "linearSystemPETSc.h"
#include "robustPredicates.h"
#include "meshGFaceOptimize.h"
#include "GFaceCompound.h"
@@ -832,47 +833,7 @@ static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoi
return;
}
-//--------------------------------------------------------------
-static bool checkOrientation(std::vector<MElement *> &elements,
- std::map<MVertex*,SPoint2> &solution, int iter) {
-
- double a_old = 0, a_new;
- bool oriented = true;
- int count = 0;
- for(unsigned int i = 0; i < elements.size(); ++i){
- MElement *t = elements[i];
- SPoint2 v1 = solution[t->getVertex(0)];
- SPoint2 v2 = solution[t->getVertex(1)];
- SPoint2 v3 = solution[t->getVertex(2)];
- double p0[2] = {v1[0],v1[1]};
- double p1[2] = {v2[0],v2[1]};
- double p2[2] = {v3[0],v3[1]};
- a_new = robustPredicates::orient2d(p0, p1, p2);
- if(count == 0) a_old=a_new;
- if(a_new*a_old < 0.){
- oriented = false;
- break;
- }
- else{
- a_old = a_new;
- }
- count++;
- }
-
- int iterMax = 1;
- if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Debug("Parametrization is NOT 1 to 1 : applying cavity checks.");
- Msg::Debug("Cavity Check - iter %d -",iter);
- one2OneMap(elements, solution);
- return checkOrientation(elements, solution, iter+1);
- }
- else if (iter < iterMax){
- Msg::Debug("Parametrization is 1 to 1 :-)");
- }
- return oriented;
-
-}
//--------------------------------------------------------------
multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
@@ -962,11 +923,7 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
//Parametrize level
std::map<MVertex*,SPoint2> solution;
- parametrize_method(level, allNodes, solution, HARMONIC);
- if (!checkOrientation(level.elements, solution, 0)){
- Msg::Debug("Parametrization switched to convex combination");
- parametrize_method(level, allNodes, solution, CONVEXCOMBINATION);
- }
+ parametrize_method(level, allNodes, solution);
//Compute the bbox of the parametric space
SBoundingBox3d bbox;
@@ -983,7 +940,7 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
- if (local_size < 1.e-5*global_size) //1.e-6
+ if (local_size < 1.e-6*global_size)
tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -1042,28 +999,6 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
regions.clear();
connectedRegions (tooSmall,regions);
- //ensure that small elements are circular patches
- // for (int i=0;i< regions.size() ; i++){
- // SPoint2 c = cents[i];
- // double rad = rads[i];
- // for (std::vector<MElement*>::iterator it = regions[i].begin(); it != regions[i].end(); ++it){
- // MElement *e = *it;
- // SPoint2 p0 = solution[e->getVertex(0)];
- // SPoint2 p1 = solution[e->getVertex(1)];
- // SPoint2 p2 = solution[e->getVertex(2)];
- // SPoint2 ec = (.3*(p0.x()+p1.x()+p2.x()), .3*(p0.y()+p1.y()+p2.y()));
- // double dist = sqrt((ec.x()-c.x())*(ec.x()-c.x()) + (ec.y()-c.y())*(ec.y()-c.y()));
- // std::vector<MElement*>::iterator itp;
- // if (dist > 0.5*rad ) {
- // goodSize.push_back(e);
- // itp = it;
- // it++;
- // regions[i].erase(itp);
- // }
- // }
- // keepConnected(regions[i], goodSize);
- //}
-
level.elements.clear();
level.elements = goodSize;
@@ -1132,13 +1067,12 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
- std::map<MVertex*,SPoint2> &solution,
- typeOfMapping tom)
+ std::map<MVertex*,SPoint2> &solution)
{
linearSystem<double> *_lsys;
-#if defined(HAVE_TAUCS)
- _lsys = new linearSystemCSRTaucs<double>;
+#if defined(HAVE_PETSC)
+ _lsys = new linearSystemPETSc<double>;
#elif defined(HAVE_GMM)
linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
_lsysb->setGmres(1);
@@ -1167,12 +1101,7 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
// assemble
femTerm<double> *mapping;
- if (tom == HARMONIC){
- mapping = new laplaceTerm(0, 1, &ONE);
- }
- else if (tom == CONVEXCOMBINATION){
- mapping = new convexCombinationTerm(0, 1, &ONE);
- }
+ mapping = new convexCombinationTerm(0, 1, &ONE);
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index b4be14fa30179c7c93a4529f0861642878bce317..751a839ca497fb27d5065065d351981fa54c519b 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -34,8 +34,7 @@ public:
std::map<MVertex*, SPoint3> &allCoordinates);
void cutElems (std::vector<MElement *> &elements);
void splitElems (std::vector<MElement *> &elements);
- typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
-
+
multiscaleLaplaceLevel* root;
void fillCoordinates (multiscaleLaplaceLevel & level,
std::map<MVertex*, SPoint3> &allCoordinates,
@@ -44,8 +43,7 @@ public:
void parametrize (multiscaleLaplaceLevel &);
void parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
- std::map<MVertex*,SPoint2> &solution,
- typeOfMapping tom);
+ std::map<MVertex*,SPoint2> &solution);
};
diff --git a/benchmarks/stl/artery.geo b/benchmarks/stl/artery.geo
index 2f9d95a1c45080e559ee79176c839ebffd577ee7..130ca3de1a31eeaa52a64376dd3af13ecf359cc0 100644
--- a/benchmarks/stl/artery.geo
+++ b/benchmarks/stl/artery.geo
@@ -1,7 +1,9 @@
-Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
-Mesh.RemeshAlgorithm=2; //(0) nosplit (1) automatic (2) split metis
+Mesh.Algorithm = 6; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
-Mesh.CharacteristicLengthFactor=0.05;
+Mesh.RemeshParametrization=0; //(0=harmonic_circle, 1=conformal, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square"
+Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
+
+Mesh.CharacteristicLengthFactor=0.1;
Merge "artery.stl";
CreateTopology;
diff --git a/benchmarks/stl/skull.geo b/benchmarks/stl/skull.geo
index 3dee9774285c953f8b8daff51314ec760be9f468..25230c5d64030f83a2e55b571f617c9d9da43400 100644
--- a/benchmarks/stl/skull.geo
+++ b/benchmarks/stl/skull.geo
@@ -1,4 +1,6 @@
-Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
+Mesh.Algorithm = 6; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+
+Mesh.RemeshParametrization=0; //(0=harmonic_circle, 1=conformal, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square"
Mesh.RemeshAlgorithm=2; //(0) nosplit (1) automatic (2) split metis
Mesh.CharacteristicLengthFactor=0.2;