diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index f77eeb18a37fa925e27d7b28ea047a7b54a00ae2..beaf1615064deb7b9fd6bf1cd3d080326393a16c 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -884,9 +884,9 @@ void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
double t1 = Cpu();
Msg::StatusBar(2, true, "(C) Done Computing Curvature (%g s)", t1-t0);
- writeToMshFile("curvature.msh");
+ //writeToMshFile("curvature.msh");
writeToPosFile("curvature.pos");
- writeToVtkFile("curvature.vtk");
+ //writeToVtkFile("curvature.vtk");
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1e2680c76f23f65925c9729b3938696ba4614698..59d308f5b1e350615c1c9f83963596fd4c81b501 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -804,21 +804,21 @@ void GFaceCompound::convexBoundary(double nTot) const
}
- // char name[256];
- // sprintf(name, "myBC-%d.pos", kk);
- // FILE * f2 = fopen(name,"w");
- // fprintf(f2, "View \"\"{\n");
- // for (int i = 0; i< oVert.size()-1; i++){
- // SPoint3 uv0 = coordinates[oVert[i]];
- // SPoint3 uv1 = coordinates[oVert[i+1]];
- // fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- // uv0.x(),uv0.y(), uv0.z(),
- // uv1.x(),uv1.y(), uv1.z(),
- // (double)i, (double)i+1);
- // }
- // fprintf(f2,"};\n");
- // fclose(f2);
- // kk++;
+ char name[256];
+ sprintf(name, "myBC-%d.pos", kk);
+ FILE * f2 = fopen(name,"w");
+ fprintf(f2, "View \"\"{\n");
+ for (int i = 0; i< oVert.size()-1; i++){
+ SPoint3 uv0 = coordinates[oVert[i]];
+ SPoint3 uv1 = coordinates[oVert[i+1]];
+ fprintf(f2, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+ uv0.x(),uv0.y(), uv0.z(),
+ uv1.x(),uv1.y(), uv1.z(),
+ (double)i, (double)i+1);
+ }
+ fprintf(f2,"};\n");
+ fclose(f2);
+ kk++;
}
@@ -889,16 +889,17 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
- bool success = orderVertices(_U0, _ordered, _coords);
- if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
-
+ if (_type != SQUARE){
+ bool success = orderVertices(_U0, _ordered, _coords);
+ if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
+ }
+
// Convex parametrization
if (_mapping == CONVEX){
Msg::Info("Parametrizing surface %d with 'convex map'", tag());
fillNeumannBCS();
parametrize(ITERU,CONVEX);
parametrize(ITERV,CONVEX);
- printStuff(11);
if (_type==MEANPLANE){
checkOrientation(0, true);
// printStuff(22);
@@ -906,45 +907,45 @@ bool GFaceCompound::parametrize() const
// parametrize(ITERU,CONVEX);
// parametrize(ITERV,CONVEX);
// checkOrientation(0, true);
- // printStuff(33);
}
- printStuff(44);
}
// Laplace parametrization
else if (_mapping == HARMONIC){
Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
- parametrize(ITERV,HARMONIC);
- printStuff(111);
+ parametrize(ITERV,HARMONIC);
if (_type == MEANPLANE) checkOrientation(0, true);
- printStuff(222);
}
// Conformal map parametrization
else if (_mapping == CONFORMAL){
-
+
fillNeumannBCS();
std::vector<MVertex *> vert;
bool oriented, hasOverlap;
if (_type == SPECTRAL){
Msg::Info("Parametrizing surface %d with 'spectral conformal map'", tag());
- hasOverlap = parametrize_conformal_spectral();
+ parametrize_conformal_spectral();
}
else if (_type == FE){
Msg::Info("Parametrizing surface %d with 'FE conformal map'", tag());
parametrize_conformal(0, NULL, NULL);
}
+ printStuff(55);
oriented = checkOrientation(0);
+ printStuff(66);
if (!oriented) oriented = checkOrientation(0, true);
+ printStuff(77);
if (_type==SPECTRAL && (!oriented || checkOverlap(vert)) ){
Msg::Warning("!!! parametrization switched to 'FE conformal' map");
parametrize_conformal(0, NULL, NULL);
oriented = checkOrientation(0, true);
}
if (!oriented || checkOverlap(vert)){
- Msg::Warning("$$$ parametrization switched to 'harmonic' map");
- parametrize(ITERU,HARMONIC);
- parametrize(ITERV,HARMONIC);
+ Msg::Warning("$$$ parametrization switched to 'convex' map");
+ _type = UNITCIRCLE;
+ parametrize(ITERU,CONVEX);
+ parametrize(ITERV,CONVEX);
}
}
// Radial-Basis Function parametrization
@@ -953,15 +954,9 @@ bool GFaceCompound::parametrize() const
int variableEps = 0;
int radFunInd = 1; // 1 MQ RBF , 0 GA
double sizeBox = getSizeH();
-
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
_rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
-
- //_rbf->RbfLapSurface_global_CPM_low(_rbf->getXYZ(), _rbf->getN(), Oper);
- //_rbf->RbfLapSurface_local_CPM(true, _rbf->getXYZ(), _rbf->getN(), Oper);
_rbf->RbfLapSurface_global_CPM_high_2(_rbf->getXYZ(), _rbf->getN(), Oper);
- //_rbf->RbfLapSurface_local_CPM(false, _rbf->getXYZ(), _rbf->getN(), Oper);
-
_rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
}
@@ -1725,15 +1720,15 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
// vertices
std::vector<MVertex *> vert;
bool hasOverlap = checkOverlap(vert);
- if ( hasOverlap && iter < 3){
- Msg::Info("Loop FE conformal iter (%d) v1=%d v2=%d", iter,
- vert[0]->getNum(), vert[1]->getNum());
- printStuff(100+iter);
- return hasOverlap = parametrize_conformal(iter+1, vert[0],vert[1]);
- }
- else{
- return hasOverlap;
- }
+ // if ( hasOverlap && iter < 3){
+ // Msg::Info("Loop FE conformal iter (%d) v1=%d v2=%d", iter,
+ // vert[0]->getNum(), vert[1]->getNum());
+ // printStuff(100+iter);
+ // return hasOverlap = parametrize_conformal(iter+1, vert[0],vert[1]);
+ // }
+ // else{
+ // return hasOverlap;
+ // }
}
void GFaceCompound::computeNormals(std::map<MVertex*,SVector3> &normals) const
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index bd87dc2215800f91cbfa8f3ed1d81ce3b0378bc4..29f4a8b202f000b538db2957b0325f5f17c6654b 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -553,9 +553,6 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
}
fields->setBackgroundField(metric);
- // int id = fields->newId();
- // (*fields)[id] = new meshMetric(this, technique, f, parameters);
- // fields->background_field = id;
opt_mesh_lc_integration_precision(0, GMSH_SET, 1.e-4);
opt_mesh_algo2d(0, GMSH_SET, 7.0); //bamg
@@ -572,7 +569,7 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
char name[256];
sprintf(name, "meshAdapt-%d.msh", ITER);
writeMSH(name);
- //metric->exportInfo(name);
+ metric->exportInfo(name);
if (ITER++ >= niter) break;
if (ITER > 3 && fabs((double)(nbElems - nbElemsOld)) < 0.01 * nbElemsOld) break;
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index ab3b5ebf246294d4c4953a80d9b1913ccdf5b6fa..040bafa1b3b0ac91c79d35c1ff50181e4c922746 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -634,7 +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::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE;
+ GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE;
+ //GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE;
//GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL_SPECTRAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 913d36d4567364614c917fd129a58b38a1e22a85..f02382dc20f17839d65e27601f5a6fb62f197fe1 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1181,6 +1181,7 @@ void bowyerWatsonFrontal(GFace *gf)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
double newPoint[2], metric[3];
+ //optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
insertAPoint(gf, AllTris.end(), newPoint, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index b1fa1ad386627ca461a7afa62f260f8a617bf682..dc3d537c460e00f0eaa82610b6c9baf0f1752407 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -109,17 +109,26 @@ void meshMetric::intersectMetrics(){
}
void meshMetric::exportInfo(const char * fileendname){
+
+ printf("printing metric \n");
+ if (_technique == meshMetric::HESSIAN) printf("HESSIANS \n");
+ else if (_technique == meshMetric::LEVELSET) printf("LEVELSETS \n");
+ else if (_technique == meshMetric::EIGENDIRECTIONS) printf("EIGEN \n");
+
if (needMetricUpdate) intersectMetrics();
- std::stringstream sg,sm,sl;
- sg << "meshmetric_gradients" << fileendname;
- sm << "meshmetric_metric" << fileendname;
- sl << "meshmetric_levelset" << fileendname;
+ std::stringstream sg,sm,sl,sh;
+ sg << "meshmetric_gradients_" << fileendname;
+ sm << "meshmetric_metric_" << fileendname;
+ sl << "meshmetric_levelset_" << fileendname;
+ sh << "meshmetric_hessian_" << fileendname;
std::ofstream out_grad(sg.str().c_str());
std::ofstream out_metric(sm.str().c_str());
std::ofstream out_ls(sl.str().c_str());
+ std::ofstream out_hess(sh.str().c_str());
out_grad << "View \"ls_gradient\"{" << std::endl;
out_metric << "View \"metric\"{" << std::endl;
out_ls << "View \"ls\"{" << std::endl;
+ out_hess << "View \"hessian\"{" << std::endl;
std::vector<MElement*>::iterator itelem = _elements.begin();
std::vector<MElement*>::iterator itelemen = _elements.end();
for (;itelem!=itelemen;itelem++){
@@ -127,27 +136,41 @@ void meshMetric::exportInfo(const char * fileendname){
out_metric << "TT(";
out_grad << "VT(";
out_ls << "ST(";
+ out_hess << "ST(";
for ( int i = 0; i < e->getNumVertices(); i++) {
MVertex *ver = e->getVertex(i);
out_metric << ver->x() << "," << ver->y() << "," << ver->z();
out_grad << ver->x() << "," << ver->y() << "," << ver->z();
out_ls << ver->x() << "," << ver->y() << "," << ver->z();
+ out_hess << ver->x() << "," << ver->y() << "," << ver->z();
if (i!=e->getNumVertices()-1){
out_metric << ",";
out_grad << ",";
out_ls << ",";
+ out_hess << ",";
}
else{
out_metric << "){";
out_grad << "){";
out_ls << "){";
+ out_hess << "){";
}
}
for ( int i = 0; i < e->getNumVertices(); i++) {
MVertex *ver = e->getVertex(i);
out_ls << vals[ver];
- if ((i==(e->getNumVertices()-1))) out_ls << "};" << std::endl;
- else out_ls << ",";
+ SVector3 gradudx = dgrads[0][ver];
+ SVector3 gradudy = dgrads[1][ver];
+ SVector3 gradudz = dgrads[2][ver];
+ out_hess << (gradudx(0)+gradudy(1)+gradudz(2));
+ if ((i==(e->getNumVertices()-1))){
+ out_ls << "};" << std::endl;
+ out_hess << "};" << std::endl;
+ }
+ else {
+ out_ls << ",";
+ out_hess << ",";
+ }
for (int k=0;k<3;k++){
out_grad << grads[ver](k);
if ((k==2)&&(i==(e->getNumVertices()-1))) out_grad << "};" << std::endl;
@@ -163,9 +186,13 @@ void meshMetric::exportInfo(const char * fileendname){
out_grad << "};" << std::endl;
out_metric << "};" << std::endl;
out_ls << "};" << std::endl;
+ out_hess << "};" << std::endl;
out_grad.close();
out_metric.close();
out_ls.close();
+ out_hess.close();
+
+ //exit(1);
}
@@ -329,21 +356,7 @@ void meshMetric::computeMetric(){
hfrey(1,0) = hfrey(0,1) = C*gr(1)*gr(0)/(norm) + hessian(1,0)/epsGeom;
hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) + hessian(2,0)/epsGeom;
hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) + hessian(2,1)/epsGeom;
- // hfrey(0,0) += C*gr(0)*gr(0)/norm;
- // hfrey(1,1) += C*gr(1)*gr(1)/norm;
- // hfrey(2,2) += C*gr(2)*gr(2)/norm;
- // hfrey(1,0) = hfrey(0,1) = gr(1)*gr(0)/(norm) ;
- // hfrey(2,0) = hfrey(0,2) = gr(2)*gr(0)/(norm) ;
- // hfrey(2,1) = hfrey(1,2) = gr(2)*gr(1)/(norm) ;
}
- // SMetric3 sss=hessian;
- // sss *= divEps;
- // sss(0,0) += 1/(hmax*hmax);
- // sss(1,1) += 1/(hmax*hmax);
- // sss(2,2) += 1/(hmax*hmax);
- // H = intersection(sss,hfrey);
- // if (dist < _E) H = intersection(sss,hfrey);
- // else H = hfrey;
H = hfrey;
}
else if ((_technique == meshMetric::EIGENDIRECTIONS )||(_technique == meshMetric::EIGENDIRECTIONS_LINEARINTERP_H )){
@@ -396,16 +409,11 @@ void meshMetric::computeMetric(){
std::vector<int> ti_index;
for (int i=0;i<3;i++)
if (i!=grad_index) ti_index.push_back(i);
-// std::cout << "gr tgr t1 t2 dots_tgr_gr_t1_t2 (" << gr(0) << "," << gr(1) << "," << gr(2) << ") (" << ti[grad_index](0) << "," << ti[grad_index](1) << "," << ti[grad_index](2) << ") (" << ti[ti_index[0]](0) << "," << ti[ti_index[0]](1) << "," << ti[ti_index[0]](2) << ") (" << ti[ti_index[1]](0) << "," << ti[ti_index[1]](1) << "," << ti[ti_index[1]](2) << ") " << fabs(dot(ti[grad_index],gr)) << " " << (dot(ti[grad_index],ti[ti_index[0]])) << " " << (dot(ti[grad_index],ti[ti_index[1]])) << std::endl;
-
// finally, creating the metric
std::vector<double> eigenvals;
eigenvals.push_back(std::min(std::max(eigenval_direction,metric_value_hmax),metric_value_hmin));// in gradient direction
eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[0]],metric_value_hmax),metric_value_hmin));
eigenvals.push_back(std::min(std::max(eigenvals_curvature[ti_index[1]],metric_value_hmax),metric_value_hmin));
-// eigenvals.push_back(std::min(std::max(metric_value_hmax,metric_value_hmax),metric_value_hmin));
-// eigenvals.push_back(std::min(std::max(metric_value_hmax,metric_value_hmax),metric_value_hmin));
-
metric = SMetric3(eigenvals[0],eigenvals[1],eigenvals[2],gr,ti[ti_index[0]],ti[ti_index[1]]);
setOfSizes[metricNumber].insert(std::make_pair(ver, std::min(std::min(1/sqrt(eigenvals[0]),1/sqrt(eigenvals[1])),1/sqrt(eigenvals[2]))));
}
@@ -432,23 +440,9 @@ void meshMetric::computeMetric(){
H.eig(V,S);
double lambda1, lambda2, lambda3;
- // if (dist < _E && _technique == meshMetric::FREY){
- // fullMatrix<double> Vhess(3,3);
- // fullVector<double> Shess(3);
- // hessian.eig(Vhess,Shess);
- // double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
- // double lam1 = Shess(0);
- // double lam2 = Shess(1);
- // double lam3 = (_dim == 3)? Shess(2) : 1.;
- // lambda1 = lam1;
- // lambda2 = lam2/lambda1;
- // lambda3 = (_dim == 3)? lam3/lambda1: 1.0;
- // }
- // else{
lambda1 = S(0);
lambda2 = S(1);
lambda3 = (_dim == 3)? S(2) : 1.;
- //}
if (_technique == meshMetric::HESSIAN || (dist < _E && _technique == meshMetric::FREY)){
lambda1 = std::min(std::max(fabs(S(0))/_epsilon,1./(hmax*hmax)),1./(hmin*hmin));
@@ -471,6 +465,8 @@ void meshMetric::computeMetric(){
}
}
+ //exportInfo("EMI");
+
//Adapt epsilon
//-----------------
@@ -553,36 +549,6 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
}
}
-void meshMetric::printMetric(const char* n){
- if (needMetricUpdate) intersectMetrics();
- if (!setOfMetrics.size()){
- std::cout << "meshMetric::printMetric : No metric defined ! " << std::endl;
- throw;
- }
-
- FILE *f = fopen (n,"w");
- fprintf(f,"View \"\"{\n");
-
- //std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin();
- std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
- //for (; it != _nodalMetrics.end(); ++it){
- for (; it != _hessian.end(); ++it){
- MVertex *v = it->first;
- SMetric3 h = it->second;
- double lapl = h(0,0)+h(1,1)+h(2,2);
- //fprintf(f, "SP(%g,%g,%g){%g};\n", it->first->x(),it->first->y(),it->first->z(),lapl);
- fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
- it->first->x(),it->first->y(),it->first->z(),
- h(0,0),h(0,1),h(0,2),
- h(1,0),h(1,1),h(1,2),
- h(2,0),h(2,1),h(2,2)
- );
-
- }
- fprintf(f,"};\n");
- fclose (f);
-}
-
double meshMetric::getLaplacian (MVertex *v) {
MVertex *vNew = _vertexMap[v->getNum()];
std::map<MVertex*, SMetric3 >::const_iterator it = _hessian.find(vNew);
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 262f418716edb9cab590b98dc575a30fee414c81..bd166dcabd0149d6ca256b7fc97240bb6cd406b7 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -96,7 +96,6 @@ class meshMetric: public Field {
virtual double operator() (double x, double y, double z, GEntity *ge=0) ;
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
- void printMetric(const char* n);
// export pos files of fct, fct gradients (fct is the lattest fct passed to meshMetric !!) and resulting metric (intersection of all computed metrics)
void exportInfo(const char *fileendname);
};
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index aaaf4ec24ee406f73d5a85ac82acd6a6519f54b1..a0a5c9b2a6d9cb62850009416f79a8e3c3f981ab 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -876,8 +876,8 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
//Split the mesh in left and right
//or Cut the mesh in left and right
- splitElems(elements);
- //cutElems(elements);
+ //splitElems(elements);
+ cutElems(elements);
}
diff --git a/benchmarks/3d/CubeAniso.geo b/benchmarks/3d/CubeAniso.geo
index 03ca0392a9425701c119d1b492bc456d9ed01690..ea3963d7f85ea3fd47a49fa34e19b171d6249c16 100644
--- a/benchmarks/3d/CubeAniso.geo
+++ b/benchmarks/3d/CubeAniso.geo
@@ -50,7 +50,7 @@ Field[2] = MathEvalAniso;
Field[2].m11 = "1./(0.1)^2";
Field[2].m12 = "0";
Field[2].m13 = "0";
-Field[2].m22 = "y+1/(0.001)^2";
+Field[2].m22 = "y+1/(0.005)^2";
Field[2].m23 = "0";
Field[2].m33 = "1/(0.1)^2";
Background Field = 2;
diff --git a/benchmarks/curvature/test.geo b/benchmarks/curvature/test.geo
index 6135e8eab3b729c79af89d64857bf962c9fc423a..2e8167cd6dda0f70669c9b4770dd9c825dbb2cc2 100644
--- a/benchmarks/curvature/test.geo
+++ b/benchmarks/curvature/test.geo
@@ -4,25 +4,39 @@
Mesh.Algorithm=1; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
//Mesh.CharacteristicLengthFactor= 0.1;
+
+ Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthMin = 0.001;
- Mesh.CharacteristicLengthMax = 10.00;
+ Mesh.CharacteristicLengthMax = 3.00;
Mesh.LcIntegrationPrecision=1.e-5;
Mesh.MinimumCirclePoints=20;
Mesh.CharacteristicLengthExtendFromBoundary=0;
Merge "Sphere.stl";
//Merge "TorusInput.stl";
+
+
//Merge "DistanceFromOrigin.bgm";
//Merge "hop.bgm";
//Background Mesh View[0];
-
//Field[1] = MathEval;
//Field[1].F = "0.1*sqrt(x^2 + y^2 + z^2)";
//Background Field = 1;
CreateTopology;
- Compound Surface(20) = {1};
- Physical Surface("Wall") = {20};
+ll[] = Line "*";
+For j In {0 : #ll[]-1}
+ Compound Line(newl) = ll[j];
+EndFor
+
+ss[] = Surface "*";
+s = news;
+For i In {0 : #ss[]-1}
+ Compound Surface(s+i) = ss[i];
+EndFor
+
+
+Physical Surface(1) = {s : s + #ss[]-1};
diff --git a/benchmarks/stl/BifurcationRemeshCurvature.geo b/benchmarks/stl/BifurcationRemeshCurvature.geo
index 02515267cdbbe4f0ac4c8a4e9ae47ce360af8eb7..fb8dc04fe359e5ceb1fd28ad857bced8e6a4f433 100644
--- a/benchmarks/stl/BifurcationRemeshCurvature.geo
+++ b/benchmarks/stl/BifurcationRemeshCurvature.geo
@@ -1,23 +1,22 @@
-//Mesh.Algorithm=7; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=BAMG
-Mesh.RemeshParametrization=0; //0) harmonic (1) conformal
-Mesh.RemeshAlgorithm=0; //(0) nosplit (1) automatic
+Mesh.Algorithm=7; //1=MeshAdapt, 5=Delaunay, 6=Frontal, 7=BAMG
-//Mesh.CharacteristicLengthFactor=0.35;
-//Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
-//Mesh.CharacteristicLengthMin = 0.02; //-clmin 0.05
-//Mesh.CharacteristicLengthMax = 4.00; //-clmax 4.0
-//Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
-//Mesh.MinimumCirclePoints=30; //default=7
-//Mesh.CharacteristicLengthFromPoints = 0;
-//Mesh.CharacteristicLengthExtendFromBoundary=0;
+Mesh.RemeshParametrization=1; //(0=harmonic_circle, 1=conformal_spectral, 2=rbf, 3=harmonic_plane, 4=convex_circle, 5=convex_plane, 6=harmonic square, 7=spectral_fe
+Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic
+
+Mesh.CharacteristicLengthFactor=0.35;
+Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+Mesh.CharacteristicLengthMin = 0.02; //-clmin 0.05
+Mesh.CharacteristicLengthMax = 4.00; //-clmax 4.0
+Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+Mesh.MinimumCirclePoints=30; //default=7
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthExtendFromBoundary=0;
Merge "bifurcation.stl";
CreateTopology;
-
Compound Surface(20) = {1};
-
Physical Surface("Wall") = {20};