diff --git a/CMakeLists.txt b/CMakeLists.txt
index 78caf910c30e3705cde6e606f62f494923d9a4cc..5d3feb2718a2b770c24f55cfd70739298b5ff502 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,7 +17,9 @@ endif(DEFINED CMAKE_BUILD_TYPE)
project(gmsh CXX C)
option(ENABLE_ACIS "Enable ACIS geometrical models" ON)
-option(ENABLE_ANN "Enable ANN to compute Approximate Nearest Neighbors" ON)
+option(ENABLE_ANN "Enable ANN to compute Approximate Nearest
+Neighbors" ON)
+option(ENABLE_RBF "Enable RBF project" OFF)
option(ENABLE_BAMG "Enable Bamg mesh generator" ON)
option(ENABLE_BFGS "Enable BFGS" ON)
option(ENABLE_BLAS_LAPACK "Use BLAS/Lapack for basic linear algebra" ON)
@@ -454,6 +456,11 @@ if(ENABLE_ANN)
set_config_option(HAVE_ANN "Ann")
endif(ENABLE_ANN)
+if(ENABLE_RBF)
+ include_directories(projects/rbf)
+ set_config_option(HAVE_RBF "RBF")
+endif(ENABLE_RBF)
+
if(ENABLE_CHACO)
add_subdirectory(contrib/Chaco)
include_directories(contrib/Chaco/main)
diff --git a/Common/GmshConfig.h.in b/Common/GmshConfig.h.in
index 8997d6166b192d610d5905acb387ee6eb4bffdf6..4b0eef98b7b02a2c9f7cd420b73872ce068b99be 100644
--- a/Common/GmshConfig.h.in
+++ b/Common/GmshConfig.h.in
@@ -52,6 +52,7 @@
#cmakedefine HAVE_SOLVER
#cmakedefine HAVE_TAUCS
#cmakedefine HAVE_TETGEN
+#cmakedefine HAVE_RBF
#define GMSH_CONFIG_OPTIONS "${GMSH_CONFIG_OPTIONS}"
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index c8a0e169425b2c907853ee0461dbb78e9eae6cf4..9a3b33bf2f67e988e92b80634d44ff4d1667ae15 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -40,8 +40,9 @@
#include "discreteFace.h"
#include "eigenSolver.h"
#include "multiscaleLaplace.h"
-//#include "rbf.h"
-
+#ifdef HAVE_RBF
+#include "rbf.h"
+#endif
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
@@ -520,9 +521,7 @@ bool GFaceCompound::checkOverlap() const
// the mapped triangles have the same normal orientation
bool GFaceCompound::checkOrientation(int iter) const
{
- //Only check orientation for stl files (1 patch)
- //if(_compound.size() > 1.0) return true;
-
+
std::list<GFace*>::const_iterator it = _compound.begin();
double a_old = 0.0, a_new=0.0;
bool oriented = true;
@@ -560,8 +559,8 @@ bool GFaceCompound::checkOrientation(int iter) const
return checkOrientation(iter+1);
}
else if (oriented && iter < iterMax){
- //Msg::Info("Parametrization is bijective (no flips)");
- //printStuff();
+ Msg::Info("Parametrization is bijective (no flips)");
+ printStuff();
}
return oriented;
@@ -631,7 +630,6 @@ bool GFaceCompound::parametrize() const
// Laplace parametrization
if (_mapping == HARMONIC){
- printf("Parametrizing surface %d with 'harmonic map' \n", tag());
Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag());
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
@@ -662,21 +660,22 @@ bool GFaceCompound::parametrize() const
}
// Radial-Basis Function parametrization
else if (_mapping == RBF){
- printf("Parametrizing surface %d with 'RBF' \n", tag());
Msg::Debug("Parametrizing surface %d with 'RBF''", tag());
}
buildOct();
-
- if (!checkOrientation(0)){
- Msg::Info("### Flipping: parametrization switched to convex combination map");
- coordinates.clear();
- Octree_Delete(oct);
- fillNeumannBCS();
- parametrize(ITERU,CONVEXCOMBINATION);
- parametrize(ITERV,CONVEXCOMBINATION);
- checkOrientation(0);
- buildOct();
+
+ if (_mapping != RBF){
+ if (!checkOrientation(0) ){
+ Msg::Info("### Flipping: parametrization switched to convex combination map");
+ coordinates.clear();
+ Octree_Delete(oct);
+ fillNeumannBCS();
+ parametrize(ITERU,CONVEXCOMBINATION);
+ parametrize(ITERV,CONVEXCOMBINATION);
+ checkOrientation(0);
+ buildOct();
+ }
}
double AR = checkAspectRatio();
@@ -688,12 +687,14 @@ bool GFaceCompound::parametrize() const
return paramOK;
}
-void GFaceCompound::setParam(std::map<MVertex*, SPoint3> rbf_param) {//, rbf *myRBF){
+#ifdef HAVE_RBF
+void GFaceCompound::setParam(std::map<MVertex*, SPoint3> rbf_param, rbf *myRBF){
- //_rbf = myRBF;
+ _rbf = myRBF;
coordinates = rbf_param;
-
+
}
+#endif
void GFaceCompound::getBoundingEdges()
{
@@ -714,11 +715,6 @@ void GFaceCompound::getBoundingEdges()
l_edges.push_back(*it);
(*it)->addFace(this);
}
- it = _V0.begin();
- for( ; it != _V0.end() ; ++it){
- l_edges.push_back(*it);
- (*it)->addFace(this);
- }
std::list<GEdge*> loop;
computeALoop(_unique, loop);
while(!_unique.empty()) computeALoop(_unique,loop);
@@ -920,11 +916,9 @@ 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, std::list<GEdge*> &V0,
- std::list<GEdge*> &U1, std::list<GEdge*> &V1,
linearSystem<double> *lsys, typeOfMapping mpg,
int allowPartition)
- : GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
+ : GFace(m, tag), _compound(compound), oct(0),
_lsys(lsys),_mapping(mpg), _allowPartition(allowPartition)
{
ONE = new simpleFunction<double>(1.0);
@@ -952,9 +946,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
}
getBoundingEdges();
- if(!_U0.size()) _type = UNITCIRCLE;
- else if(!_V1.size()) _type = UNITCIRCLE;
- else _type = SQUARE;
+ _type = UNITCIRCLE;
nbSplit = 0;
fillTris.clear();
@@ -1082,28 +1074,6 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));
}
}
- else if(_type == SQUARE){
- if(step == ITERU){
- std::list<GEdge*>::const_iterator it = _U0.begin();
- for( ; it != _U0.end(); ++it){
- fixEdgeToValue(*it, 0.0, myAssembler);
- }
- it = _U1.begin();
- for( ; it != _U1.end(); ++it){
- fixEdgeToValue(*it, 1.0, myAssembler);
- }
- }
- else if(step == ITERV){
- std::list<GEdge*>::const_iterator it = _V0.begin();
- for( ; it != _V0.end(); ++it){
- fixEdgeToValue(*it, 0.0, myAssembler);
- }
- it = _V1.begin();
- for( ; it != _V1.end(); ++it){
- fixEdgeToValue(*it, 1.0, myAssembler);
- }
- }
- }
else{
Msg::Error("Unknown type of parametrization");
return;
@@ -1173,229 +1143,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
_lsys->clear();
}
-void GFaceCompound::computeThetaDerivatives (MVertex *prev, MVertex *curr, MVertex *next,
- double &dTdu1, double &dTdv1,
- double &dTdu2, double &dTdv2,
- double &dTdu3, double &dTdv3) const
-{
- SPoint2 p1 = getCoordinates(prev);
- SPoint2 p2 = getCoordinates(curr);
- SPoint2 p3 = getCoordinates(next);
- SVector3 va(p2.x()-p1.x(), p2.y()-p1.y(),0.);
- SVector3 vb(p3.x()-p2.x(), p3.y()-p2.y(), 0.);
- SVector3 vc(p1.x()-p3.x(), p1.y()-p3.y(), 0.);
- double a = norm(va), b=norm(vb), c=norm(vc);
- SVector3 n = crossprod(va, vb);
-
- double sign = 1.;
- if (n.z() < 0.0) sign = -1.;
-
- //- compute grad_uv theta
- //dTheta/du1 sign*acos((a^2+b^2 -c^2)/(2*a*b))
- //a=sqrt((u2-u1)^2+(v2-v1)^2))
- //b=sqrt((u3-u2)^2+(v3-v2)^2))
- //c=sqrt((u1-u3)^2+(v1-v3)^2))
- double u1 = p1.x();double v1 = p1.y();
- double u2 = p2.x();double v2 = p2.y();
- double u3 = p3.x();double v3 = p3.y();
-
- dTdu1 = sign*(v1*v3*u1-v1*v3*u2-v3*v2*u1+v3*v2*u2-v2*v1*u1+2*u3*v2*v1+u2*SQU(v1)-u3*SQU(v2)-u3*SQU(v1)+SQU(v2)*u1-v2*v1*u2)/pow(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1), 3./2.)/sqrt(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2))/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
- dTdv1 = sign*(v1*u1*u3-SQU(u1)*v3+v1*SQU(u2)-u2*u3*v1-v1*u1*u2-u2*u1*v2+u3*u2*v2-u1*u3*v2+2*v3*u2*u1-v3*SQU(u2)+v2*SQU(u1))/pow(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1),3./2.)/sqrt(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2))/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
- dTdu2= sign*(u1*SQU(v1)*SQU(v2)+u1*SQU(v1)*SQU(v3)+v2*CUB(v1)*u3-u2*v2*CUB(v1)-v3*CUB(v2)*u3+CUB(u1)*SQU(v3)+CUB(u1)*SQU(v2)+5*u2*v2*v1*SQU(u3)-3*SQU(u2)*v2*v1*u3+u2*v2*v1*SQU(v3)-2*u2*SQU(v2)*v1*v3+u1*v2*v1*SQU(u3)+3*u1*v2*v1*SQU(u2)+u1*v2*v1*SQU(v3)+u1*SQU(v2)*v1*v3+4*u3*u2*u1*SQU(v2)-u2*v3*v2*SQU(u3)+3*SQU(u2)*v3*v2*u3+u1*v3*v2*SQU(u3)-3*u1*v3*v2*SQU(u2)-u1*v1*v3*SQU(u3)+v2*v1*u3*SQU(u1)-u2*v2*v1*SQU(u1)+v3*SQU(v2)*u3*v1+v3*v2*u3*SQU(v1)-3*u3*v1*v3*SQU(u2)-u3*v1*v3*SQU(u1)+u2*v1*v3*SQU(u1)-4*u1*v2*v1*u3*u2-u1*CUB(v2)*v1-u2*CUB(v3)*v2+u2*SQU(v3)*SQU(v2)+u1*CUB(v3)*v2-2*u1*SQU(v3)*SQU(v2)+u1*v3*CUB(v2)-u1*v1*CUB(v3)-SQU(v2)*u3*SQU(u1)+CUB(v2)*u3*v1-2*SQU(v2)*u3*SQU(v1)+3*u3*SQU(u2)*SQU(v1)+4*u1*v1*v3*u3*u2-3*u1*v1*v3*SQU(u2)+SQU(u1)*u3*v3*v2+u2*SQU(v1)*v3*v2+5*u2*SQU(u1)*v3*v2-2*u1*SQU(v1)*v3*v2-4*u3*u2*u1*v3*v2-2*u2*SQU(v2)*SQU(u1)+u2*SQU(v2)*SQU(v1)-2*SQU(u3)*u2*SQU(v2)-3*SQU(u3)*u2*SQU(v1)-u1*SQU(u3)*SQU(v2)+3*SQU(u2)*u1*SQU(v3)+SQU(v1)*CUB(u3)-2*u2*SQU(v1)*SQU(v3)-3*u2*SQU(u1)*SQU(v3)+u2*v1*v3*SQU(u3)-2*u3*v2*v1*SQU(v3)-2*CUB(u1)*v3*v2-CUB(u2)*SQU(v3)-SQU(v1)*CUB(u2)-u3*CUB(v1)*v3+2*CUB(u2)*v1*v3+u2*CUB(v1)*v3+u2*v1*CUB(v3)-2*v2*v1*CUB(u3)+u3*SQU(v1)*SQU(v3)+u3*SQU(v2)*SQU(v3)+SQU(v2)*CUB(u3))/pow(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1), 3./2.)/pow(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2),3./2.)/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
- dTdv2= -sign*(CUB(u3)*u2*v2-u1*CUB(u3)*v2-2*u1*u3*CUB(v2)+3*u1*u3*SQU(v2)*v1-u1*u3*v2*SQU(v1)+v3*u1*u3*SQU(v1)-4*v3*u1*u3*v2*v1+v1*u1*u3*SQU(v3)+u3*u2*v2*SQU(v3)-u2*u3*v1*SQU(v3)+v3*CUB(u1)*u3-u2*CUB(u3)*v1+u2*CUB(u1)*v2+v1*u1*CUB(u3)+v1*SQU(u2)*SQU(v3)+v1*CUB(u2)*u1-v1*CUB(u2)*u3+2*v1*SQU(u2)*SQU(u3)-CUB(u2)*u1*v3-v2*SQU(u1)*u3*u2+2*v2*SQU(u2)*u1*u3-v2*u2*u1*SQU(u3)-u1*u3*v2*SQU(v3)+3*u1*u3*SQU(v2)*v3+2*SQU(u2)*SQU(u1)*v3-u2*CUB(u1)*v3+4*u2*u1*v3*v2*v1-CUB(v1)*SQU(u3)-CUB(v1)*SQU(u2)+2*CUB(v1)*u3*u2+2*v2*SQU(v1)*SQU(u2)+3*v2*SQU(v1)*SQU(u3)-3*v1*SQU(v2)*SQU(u3)-v1*SQU(u1)*SQU(u2)-v1*SQU(u1)*SQU(u3)+3*SQU(u1)*SQU(v3)*v2+2*u2*u1*CUB(v3)+2*SQU(u2)*SQU(v3)*v2+v3*SQU(v1)*SQU(u2)-3*v3*SQU(u1)*SQU(v2)-5*u2*u1*SQU(v3)*v2+u2*u1*v2*SQU(v1)-SQU(u1)*CUB(v3)-SQU(u2)*CUB(v3)+4*v3*v2*v1*u3*u2-v3*SQU(u1)*SQU(u3)+v3*CUB(u2)*u3-v3*SQU(u2)*SQU(u3)-SQU(u2)*v2*SQU(u3)+3*u2*u1*v3*SQU(v2)-u2*u1*v3*SQU(v1)-4*v3*v2*v1*SQU(u2)-5*v2*SQU(v1)*u3*u2+3*v1*SQU(v2)*u3*u2-CUB(u1)*u3*v2-v2*SQU(u1)*SQU(u2)+2*v2*SQU(u1)*SQU(u3)+SQU(u1)*CUB(v2)+CUB(v2)*SQU(u3)+2*v1*SQU(u1)*u3*u2-3*v1*u2*u1*SQU(v2)-v1*u2*u1*SQU(v3)-v1*SQU(u2)*u1*u3-v1*u2*u1*SQU(u3)-v3*SQU(v1)*u3*u2-3*v3*SQU(v2)*u3*u2-v3*SQU(u1)*u3*u2-v3*SQU(u2)*u1*u3+2*v3*u2*u1*SQU(u3))/pow(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1),3./2.)/pow(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2), 3./2.)/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
- dTdu3= -sign*(-SQU(v2)*u3-2*u1*v3*v2-u2*SQU(v3)+u2*v3*v2+v2*v1*u3+u2*v1*v3+u1*SQU(v3)-u2*v2*v1+u1*SQU(v2)+v3*v2*u3-u3*v1*v3)/sqrt(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/pow(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2), 3./2.)/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
- dTdv3= sign*(v3*u1*u3-v1*SQU(u3)+2*u2*u3*v1-v1*SQU(u2)+u2*u1*v2-u2*u1*v3-u3*u2*v2-u3*u2*v3-u1*u3*v2+SQU(u2)*v3+v2*SQU(u3))/sqrt(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/pow(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2), 3./2.)/sqrt(SQU(u2*v3-u1*v3+u1*v2-u2*v1-v2*u3+v1*u3)/(SQU(u2)-2*u2*u1+SQU(u1)+SQU(v2)-2*v2*v1+SQU(v1))/(SQU(u3)-2*u3*u2+SQU(u2)+SQU(v3)-2*v3*v2+SQU(v2)));
-
-}
-
-bool GFaceCompound::parametrize_conformal_nonLinear() const
-{
- bool converged = false;
-
- //--create dofManager
- linearSystem<double>* lsysNL;
- //lsysNL = new linearSystemFull<double>;
- //lsysNL = new linearSystemPETSc<double>;
- lsysNL = new linearSystemCSRTaucs<double>;
-
- dofManager<double> myAssembler(lsysNL);
- //--- first compute mapping harmonic
- //parametrize(ITERU,HARMONIC);
- //parametrize(ITERV,HARMONIC);
- printStuff(100);
-
- //---order boundary vertices
- // std::vector<MVertex*> ordered;
- // std::vector<double> coords;
- // bool success = orderVertices(_U0, ordered, coords);
- // int nb = ordered.size();
-
- //--fix vertex for du=0 and dv=0
- MVertex *v1 = _ordered[0];
- MVertex *v2 = _ordered[1]; //(int)ceil((double)_ordered.size()/2.)];
- myAssembler.fixVertex(v1, 0, 1, 0.);
- myAssembler.fixVertex(v1, 0, 2, 0.);
- myAssembler.fixVertex(v2, 0, 1, 0.);
- myAssembler.fixVertex(v2, 0, 2, 0.);
-
- //--Assemble linear system
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- myAssembler.numberVertex(v, 0, 1);
- myAssembler.numberVertex(v, 0, 2);
- }
- MVertex *lag = new MVertex(0.,0.,0.);
- myAssembler.numberVertex(lag, 0, 3);//ghost vertex for lagrange multiplier
-
- //--- newton Loop
- int nbNewton = 5;
- double lambda = 1.e6;
- double fac = 1.0; //1.e7;
- for (int iNewton = 0; iNewton < nbNewton; iNewton++){
-
- //-- assemble conformal matrix
- std::vector<MElement *> allElems;
- laplaceTerm laplace1(model(), 1, ONE, &coordinates);
- laplaceTerm laplace2(model(), 2, ONE, &coordinates);
- crossConfTerm cross12(model(), 1, 2, ONE, &coordinates);
- crossConfTerm cross21(model(), 2, 1, MONE, &coordinates);
- std::list<GFace*>::const_iterator it = _compound.begin();
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- SElement se((*it)->triangles[i]);
- laplace1.addToMatrix(myAssembler, &se);
- laplace2.addToMatrix(myAssembler, &se);
- cross12.addToMatrix(myAssembler, &se);
- cross21.addToMatrix(myAssembler, &se);
- allElems.push_back((MElement*)((*it)->triangles[i]));
- }
- }
-
- groupOfElements gr(allElems);
- laplace1.addToRightHandSide(myAssembler, gr);
- laplace2.addToRightHandSide(myAssembler, gr);
- cross12.addToRightHandSide(myAssembler, gr);
- cross21.addToRightHandSide(myAssembler, gr);
-
- //-- compute all boundary angles
- double sumTheta = 0.0;
- int nb = _ordered.size();
- for (int i=0; i< nb; i++){
- MVertex *prev = (i!=0) ? _ordered[i-1] : _ordered[nb-1];
- MVertex *curr = _ordered[i];
- MVertex *next = (i+1!=nb) ? _ordered[i+1] : _ordered[0];
- SPoint2 p1 = getCoordinates(prev);
- SPoint2 p2 = getCoordinates(curr);
- SPoint2 p3 = getCoordinates(next);
- SVector3 va(p2.x()-p1.x(), p2.y()-p1.y(),0.);
- SVector3 vb(p3.x()-p2.x(), p3.y()-p2.y(), 0.);
- SVector3 vc(p1.x()-p3.x(), p1.y()-p3.y(), 0.);
- double a = norm(va), b=norm(vb), c=norm(vc);
- SVector3 n = crossprod(va, vb);
-
- double theta = acos((a*a+b*b-c*c)/(2*a*b)); //in rad
- double sign = 1.;
- if (n.z() < 0.0) {sign = -1.; theta = 2*M_PI-theta;}
- sumTheta +=theta;
-
- //- compute grad_uv theta
- //dTheta/du1 sign*acos((a^2+b^2 -c^2)/(2*a*b))
- //a=sqrt((u2-u1)^2+(v2-v1)^2))
- //b=sqrt((u3-u2)^2+(v3-v2)^2))
- //c=sqrt((u1-u3)^2+(v1-v3)^2))
- double dTdu1, dTdv1, dTdu2, dTdv2, dTdu3, dTdv3;
- computeThetaDerivatives(prev, curr, next,
- dTdu1, dTdv1, dTdu2, dTdv2, dTdu3, dTdv3);
-
- //- assemble constraint terms
- myAssembler.assemble(lag, 0, 3, prev, 0, 1, -fac*dTdu1);
- myAssembler.assemble(lag, 0, 3, prev, 0, 2, -fac*dTdv1);
- myAssembler.assemble(lag, 0, 3, curr, 0, 1, -fac*dTdu2);
- myAssembler.assemble(lag, 0, 3, curr, 0, 2, -fac*dTdv2);
- myAssembler.assemble(lag, 0, 3, next, 0, 1, -fac*dTdu3);
- myAssembler.assemble(lag, 0, 3, next, 0, 2, -fac*dTdv3);
-
- myAssembler.assemble(prev, 0, 1, lag, 0, 3, -fac*dTdu1);
- myAssembler.assemble(prev, 0, 2, lag, 0, 3, -fac*dTdv1);
- myAssembler.assemble(curr, 0, 1, lag, 0, 3, -fac*dTdu2);
- myAssembler.assemble(curr, 0, 2, lag, 0, 3, -fac*dTdv2);
- myAssembler.assemble(next, 0, 1, lag, 0, 3, -fac*dTdu3);
- myAssembler.assemble(next, 0, 2, lag, 0, 3, -fac*dTdv3);
-
- myAssembler.assemble(prev, 0, 1, lambda*fac*dTdu1);
- myAssembler.assemble(prev, 0, 2, lambda*fac*dTdv1);
- myAssembler.assemble(curr, 0, 1, lambda*fac*dTdu2);
- myAssembler.assemble(curr, 0, 2, lambda*fac*dTdv2);
- myAssembler.assemble(next, 0, 1, lambda*fac*dTdu3);
- myAssembler.assemble(next, 0, 2, lambda*fac*dTdv3);
-
- }
-
- //--- compute constraint
- double G = sumTheta - (nb-2)*M_PI;
- printf("**NL** Sum of angles G = %g \n", G );
- myAssembler.assemble(lag, 0, 3, lag, 0, 3, 1.e-7);
- myAssembler.assemble(lag, 0, 3, fac*G);
-
- //--- solve linear system
- Msg::Debug("Assembly done");
- lsysNL->systemSolve();
- Msg::Debug("System solved");
-
- //-- update newton
- //U=U+DU
- //lambda = lambda +dLambda
- double meandu = 0.0;
- double meandv = 0.0;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- double du,dv;
- myAssembler.getDofValue(v, 0, 1, du);
- myAssembler.getDofValue(v, 0, 2, dv);
- meandu +=du; meandv +=dv;
- std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
- double oldu = itf->second.x(), oldv = itf->second.y();
- if(itf != coordinates.end()){
- itf->second[0]= oldu+du;
- itf->second[1]= oldv+dv;
- }
- }
- meandu /=allNodes.size();
- meandv /=allNodes.size();
- double dLambda;
- myAssembler.getDofValue(lag, 0, 3, dLambda);
- lambda = lambda+dLambda;
-
- lsysNL->zeroMatrix();
- lsysNL->zeroRightHandSide();
-
- //--priting
- printStuff(iNewton);
-
- //-- exit newton criteria
- bool noOverlap = checkOverlap();
- printf("**NL** ---- iNewton %d --- \n", iNewton);
- printf("**NL** System solved: du=%g, dv=%g dL=%g \n", meandu, meandv, dLambda);
- if (noOverlap) {
- if (checkOrientation(0)){
- converged = true;
- break;
- }
- }
-
- }//end Newton
-
- lsysNL->clear();
- //exit(1);
-
- return converged;
-}
bool GFaceCompound::parametrize_conformal_spectral() const
{
@@ -1678,13 +1426,17 @@ GPoint GFaceCompound::point(double par1, double par2) const
if(!oct) parametrize();
+#ifdef HAVE_RBF
if (_mapping == RBF){
- double x, y,z;
- //_rbf->UVStoXYZ(par1, par2,x,y,z);
+ double x, y, z;
+ SVector3 dXdu, dXdv;
+ _rbf->UVStoXYZ(par1, par2,x,y,z, dXdu, dXdv);
//printf("XYZ =%g %g %g \n", x,y,z);
//exit(1);
- //return GPoint(x,y,z);
+ return GPoint(x,y,z);
}
+#endif
+
double U,V;
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
@@ -1785,6 +1537,15 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
if(trivial())
return (*(_compound.begin()))->firstDer(param);
+#ifdef HAVE_RBF
+ if (_mapping == RBF){
+ double x, y, z;
+ SVector3 dXdu, dXdv ;
+ _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
+ return Pair<SVector3, SVector3>(dXdu,dXdv);
+ }
+#endif
+
double U, V;
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
@@ -1944,6 +1705,7 @@ void GFaceCompound::getTriangle(double u, double v,
void GFaceCompound::buildOct() const
{
+
SBoundingBox3d bb;
int count = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -2009,11 +1771,9 @@ void GFaceCompound::buildOct() const
bool GFaceCompound::checkTopology() const
{
- // FIXME!!! I think those things are wrong with cross-patch reparametrization
- //if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
-
- //TODO: smthg to exit here for lloyd remeshing
+ if (_mapping == RBF) return true;
+
//fixme tristan
//return true;
@@ -2274,7 +2034,7 @@ int GFaceCompound::genusGeom() const
void GFaceCompound::printStuff(int iNewton) const
{
- //if( !CTX::instance()->mesh.saveAll) return;
+ if( !CTX::instance()->mesh.saveAll) return;
std::list<GFace*>::const_iterator it = _compound.begin();
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 7e79b65bf964fb2614b3e33469281be910b5452a..bd1d33cd62eb170e385879c8bebf746641d36e99 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -19,8 +19,9 @@
#include "GEdgeCompound.h"
#include "meshGFaceOptimize.h"
#include "linearSystem.h"
-//#include "rbf.h"
-
+#ifdef HAVE_RBF
+#include "rbf.h"
+#endif
#define AR_MAX 5 //maximal geometrical aspect ratio
/*
@@ -51,21 +52,24 @@ class GFaceCompoundTriangle {
};
class Octree;
-//class rbf;
-
+#ifdef HAVE_RBF
+class rbf;
+#endif
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 {UNITCIRCLE, SQUARE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
- void setParam(std::map<MVertex*, SPoint3> rbf_param); //, rbf *myRBF);
+#ifdef HAVE_RBF
+ void setParam(std::map<MVertex*, SPoint3> rbf_param, rbf *myRBF);
+#endif
protected:
- //rbf *_rbf;
+ rbf *_rbf;
simpleFunction<double> *ONE;
simpleFunction<double> *MONE;
std::list<GFace*> _compound;
- std::list<GEdge*> _U0, _U1, _V0, _V1;
+ std::list<GEdge*> _U0;
std::list<std::list<GEdge*> > _interior_loops;
mutable int nbT;
mutable GFaceCompoundTriangle *_gfct;
@@ -85,7 +89,6 @@ class GFaceCompound : public GFace {
void parametrize(iterationStep, typeOfMapping) const;
bool parametrize_conformal() const;
bool parametrize_conformal_spectral() const;
- bool parametrize_conformal_nonLinear() const;
void compute_distance() const;
bool checkOrientation(int iter) const;
bool checkOverlap() const;
@@ -107,14 +110,9 @@ class GFaceCompound : public GFace {
SOrientedBoundingBox obb_boundEdges(const std::list<GEdge* > &elist) const;
void fillNeumannBCS() const;
/* double sumAngles(std::vector<MVertex*> ordered) const; */
- void computeThetaDerivatives (MVertex *prev, MVertex *curr, MVertex *next,
- double &dTdu1, double &dTdv1,
- double &dTdu2, double &dTdv2,
- double &dTdu3, double &dTdv3) const;
+
public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
- std::list<GEdge*> &U0, std::list<GEdge*> &U1,
- std::list<GEdge*> &V0, std::list<GEdge*> &V1,
linearSystem<double>* lsys =0,
typeOfMapping typ = HARMONIC, int allowPartition=1);
virtual ~GFaceCompound();
@@ -152,8 +150,6 @@ class GFaceCompound : public GFace {
public:
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4} typeOfMapping;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
- std::list<GEdge*> &U0, std::list<GEdge*> &U1,
- std::list<GEdge*> &V0, std::list<GEdge*> &V1,
linearSystem<double>* lsys =0,
typeOfMapping typ = HARMONIC, int allowPartition=1)
: GFace(m, tag)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index c7bd8591ebc63c2db06cec9b6a5ba9c33d8cc9f7..81950486d8b1b4288ec2089acb5e6d4b442096ef 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1321,7 +1321,7 @@ void GModel::makeDiscreteFacesSimplyConnected()
Msg::Debug("Done making discrete faces simply connected");
}
-void GModel::createTopologyFromMesh()
+void GModel::createTopologyFromMesh(int ignoreHoles)
{
Msg::StatusBar(2, true, "Creating topology from mesh...");
double t1 = Cpu();
@@ -1343,7 +1343,7 @@ void GModel::createTopologyFromMesh()
for(fiter it = firstFace(); it != lastFace(); it++)
if((*it)->geomType() == GEntity::DiscreteSurface)
discFaces.push_back((discreteFace*) *it);
- createTopologyFromFaces(discFaces);
+ createTopologyFromFaces(discFaces, ignoreHoles);
// create old format (necessary for boundary layers)
exportDiscreteGEOInternals();
@@ -1508,7 +1508,7 @@ void GModel::createTopologyFromRegions(std::vector<discreteRegion*> &discRegions
Msg::Debug("Done creating topology from regions");
}
-void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int onlyOneBoundary)
+void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int ignoreHoles)
{
Msg::Debug("Creating topology from faces...");
@@ -1598,12 +1598,22 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces, int
int nbBounds = connectedSurfaceBoundaries(myEdges, boundaries);
//EMI RBF fix
- if (onlyOneBoundary){
- nbBounds = 1;
+ if (ignoreHoles){
+ int index = 0;
+ int boundSize = 0;
+ for (int ib = 0; ib < nbBounds; ib++){
+ if (boundaries[ib].size() > boundSize){
+ boundSize = boundaries[ib].size() ;
+ index = ib;
+ }
+ }
+ std::vector<std::vector<MEdge> > new_boundaries;
+ new_boundaries.push_back(boundaries[index]);
+ boundaries = new_boundaries;
}
// create new discrete edges
- for (int ib = 0; ib < nbBounds; ib++){
+ for (int ib = 0; ib < boundaries.size(); ib++){
int numE = getMaxElementaryNumber(1) + 1;
discreteEdge *e = new discreteEdge(this, numE, 0, 0);
add(e);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 329834195be4b53f0c86785c8a44a2c41c03c130..fdad5cffd5a2c06ddc9b9188cd11fea70079a99b 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -349,9 +349,9 @@ class GModel
int removeDuplicateMeshVertices(double tolerance);
// create topology from mesh
- void createTopologyFromMesh();
+ void createTopologyFromMesh(int ignoreHoles=0);
void createTopologyFromRegions(std::vector<discreteRegion*> &discRegions);
- void createTopologyFromFaces(std::vector<discreteFace*> &pFaces, int onlyOneBoundary=0);
+ void createTopologyFromFaces(std::vector<discreteFace*> &pFaces, int ignoreHoles=0);
void makeDiscreteRegionsSimplyConnected();
void makeDiscreteFacesSimplyConnected();
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 673902aa8309fb8105ff94bb90477e8a38d26dc2..ad54023c7be0b27a02565e62599416d566fc47ba 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -255,7 +255,6 @@ GEdge *OCCFactory::addCircleArc(GModel *gm, const arcCreationMethod &method,
gp_Pnt aP3(end->x(), end->y(), end->z());
TopoDS_Edge occEdge;
-
OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
@@ -851,7 +850,7 @@ GFace *OCCFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
}
}
TopoDS_Wire myWire = wire_maker.Wire();
- if (i)myWire.Reverse();
+ //if (i)myWire.Reverse();
aGenerator.Add (myWire);
}
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 553fb0df0dcab227396e3910e39e4e67c4d1f463..e152ec7439f6482fe63f40fa867900deeb4e1477 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -175,17 +175,16 @@ int GModel::importGEOInternals()
GFace *gf = getFaceByTag(s->compound[j]);
if(gf) comp.push_back(gf);
}
- std::list<GEdge*> b[4];
+ std::list<GEdge*> U0;
for(int j = 0; j < 4; j++){
for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
- if(ge) b[j].push_back(ge);
+ if(ge) U0.push_back(ge);
}
}
int param = CTX::instance()->mesh.remeshParam;
int algo = CTX::instance()->mesh.remeshAlgo;
- f = new GFaceCompound(this, s->Num, comp,
- b[0], b[1], b[2], b[3], 0,
+ f = new GFaceCompound(this, s->Num, comp, 0,
(param == 0) ? GFaceCompound::HARMONIC :
GFaceCompound::CONFORMAL, algo);
f->meshAttributes.recombine = s->Recombine;
@@ -193,6 +192,29 @@ int GModel::importGEOInternals()
f->meshAttributes.Method = s->Method;
f->meshAttributes.extrude = s->Extrude;
add(f);
+
+ if(s->EmbeddedCurves){
+ for(int i = 0; i < List_Nbr(s->EmbeddedCurves); i++){
+ Curve *c;
+ List_Read(s->EmbeddedCurves, i, &c);
+ GEdge *e = getEdgeByTag(abs(c->Num));
+ if(e)
+ f->addEmbeddedEdge(e);
+ else
+ Msg::Error("Unknown curve %d", c->Num);
+ }
+ }
+ if(s->EmbeddedPoints){
+ for(int i = 0; i < List_Nbr(s->EmbeddedPoints); i++){
+ Vertex *v;
+ List_Read(s->EmbeddedPoints, i, &v);
+ GVertex *gv = getVertexByTag(v->Num);
+ if(gv)
+ f->addEmbeddedVertex(gv);
+ else
+ Msg::Error("Unknown point %d", v->Num);
+ }
+ }
}
else if(!f){
f = new gmshFace(this, s);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 55b78834c367157fbb325310c929f3784a294817..83aa26877eb40ede8fc7baa949abcdd56018d4f5 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -690,8 +690,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version = 1.0;
// if there are no physicals we save all the elements
- if(noPhysicalGroups()) saveAll = true;
-
+ if(noPhysicalGroups()) saveAll = true;
+
// get the number of vertices and index the vertices in a continuous
// sequence
int numVertices = indexMeshVertices(saveAll, saveSinglePartition);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4da46d9e240529f1e1fbe2060d653757d12c0d41..c2095c079d54fe290bd83520fd0c1dcb1e4c0e75 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -395,6 +395,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
bool repairSelfIntersecting1dMesh,
bool debug = true)
{
+
BDS_GeomEntity CLASS_F(1, 2);
BDS_GeomEntity CLASS_EXTERIOR(1, 3);
std::map<BDS_Point*, MVertex*> recoverMap;
@@ -1604,7 +1605,6 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
{
bool isMeshed = false;
#if defined(HAVE_SOLVER)
-
bool correctTopo = gf->checkTopology();
if (!correctTopo && gf->allowPartition()){
partitionAndRemesh((GFaceCompound*) gf);
@@ -1694,7 +1694,6 @@ void partitionAndRemesh(GFaceCompound *gf)
}
//Parametrize Compound surfaces
- std::list<GEdge*> b[4];
std::set<MVertex*> allNod;
for (int i=0; i < NF; i++){
std::list<GFace*> f_compound;
@@ -1704,8 +1703,7 @@ void partitionAndRemesh(GFaceCompound *gf)
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound,
- b[0], b[1], b[2], b[3], 0,
+ GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, 0,
gf->getTypeOfMapping());
gfc->meshAttributes.recombine = gf->meshAttributes.recombine;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 5751b845fbfcbe93ed25af19864cd924869e934f..18b4b799db30e226d1ca50b72cc9ff2c37d15d92 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -174,6 +174,7 @@ void circumCenterMetric(MTriangle *base, const double *metric,
void buildMetric(GFace *gf, double *uv, double *metric)
{
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
+
metric[0] = dot(der.first(), der.first());
metric[1] = dot(der.second(), der.first());
metric[2] = dot(der.second(), der.second());
@@ -185,8 +186,9 @@ void buildMetric(GFace *gf, double *uv, double *metric)
void buildMetric(GFace *gf, double *uv, SMetric3 & m, double *metric)
{
+
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
-
+
SVector3 x1(m(0,0) * der.first().x() +
m(1,0) * der.first().y() +
m(2,0) * der.first().z(),
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 048c15c3f246689c8ef812d10d32b8b8ac4ed98f..366bae5b86c633cf5f5c0027f32b1c4e6bdce84b 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -3,11 +3,12 @@
//Field[1].NNodesByEdge = 10;
//Background Field = 1;
-Field[1] = MathEval;
-Field[1].F = "1.0"; //0.1*x+0.1";
-Background Field = 1;
+//Field[1] = MathEval;
+//Field[1].F = "1.0"; //0.1*x+0.1";
+//Background Field = 1;
-//Mesh.CharacteristicLengthFactor=1.0;
+Mesh.Algorithm=1;
+Mesh.CharacteristicLengthFactor=1.5;
lc=0.1;
Point(1) = {0, 0, 0}; //,lc};
@@ -24,11 +25,13 @@ Plane Surface(10) = {5};
//----------------------
-//Compound Line(10)={1,2,3,4};
-//Compound Surface(100)={10};
+Compound Line(10)={1,2,3,4};
+Compound Surface(100)={10};
+
+Line {1,2,3,4} In Surface{100};
-Physical Surface(100)={10};
-Physical Line(200)={1,2,3,4};
+//Physical Surface(100)={10};
+//Physical Line(200)={1,2,3,4};
diff --git a/benchmarks/3d/sph2.geo b/benchmarks/3d/sph2.geo
index 74ff0f2d12d782c3b6992df1b7639033767d897f..61dbe66357bc8bfdcc16960aa1f1f852fe4a4186 100644
--- a/benchmarks/3d/sph2.geo
+++ b/benchmarks/3d/sph2.geo
@@ -3,7 +3,9 @@ dext = 1;
dint = 0.2;
angle = Pi/4;
nbdiv = 2;
-//Mesh.ElementOrder = 2;
+Mesh.ElementOrder = 2;
+Mesh.SecondOrderLinear = 0;
+
Point(1) = {0,0,0,lc};
Point(2) = {dint,0,0,lc};
diff --git a/benchmarks/curvature/ComputeCurvature.py b/benchmarks/curvature/ComputeCurvature.py
index 53704e844ea168bfc2ff1370e2105606d3c5f9bc..3ade0a48816996ff3f0a7f744e91ad042af2d9f6 100644
--- a/benchmarks/curvature/ComputeCurvature.py
+++ b/benchmarks/curvature/ComputeCurvature.py
@@ -1,4 +1,3 @@
-from dgpy import *
from gmshpy import *
import os
import sys
diff --git a/benchmarks/step/capot.geo b/benchmarks/step/capot.geo
index c49a17d9d40362669f239f0dfea4ae84c2eb8ac6..d1e2eda5c15dab652e10f02624e4fd0964a8a662 100644
--- a/benchmarks/step/capot.geo
+++ b/benchmarks/step/capot.geo
@@ -1,13 +1,10 @@
Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
Mesh.RemeshAlgorithm=0; //(0) nosplit (1) automatic (2) split metis
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg)
-Mesh.RecombinationAlgorithm = 1;
+Mesh.CharacteristicLengthFactor=0.2;
+//Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+//Mesh.RecombinationAlgorithm = 1;
-
-
-
-//Mesh.CharacteristicLengthFactor=0.2;
Merge "capot.brep";
Compound Line(1000) = {47,50};