From 7c7c983e102e097112d99e12e55ac319e71ae95d Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 14 Jul 2011 10:14:53 +0000
Subject: [PATCH] Add import gmshCommom.i
---
Geo/GEdgeCompound.cpp | 6 +--
Geo/GFaceCompound.cpp | 40 ++++++++++++-----
Geo/OCCFace.cpp | 1 +
Geo/discreteEdge.cpp | 2 +-
Mesh/BackgroundMesh.cpp | 7 +--
Mesh/meshGEdge.cpp | 2 +
benchmarks/curvature/Torus.geo | 8 ++--
benchmarks/curvature/TorusRemesh.geo | 23 +++++-----
benchmarks/curvature/occtorus.brep | 64 ++++++++++++++++++++++++++++
benchmarks/curvature/torus.py | 20 +++++++++
gmshpy/gmshGeo.i | 1 +
11 files changed, 140 insertions(+), 34 deletions(-)
create mode 100644 benchmarks/curvature/occtorus.brep
create mode 100644 benchmarks/curvature/torus.py
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index e951da4cd3..64970ab931 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -160,9 +160,9 @@ void GEdgeCompound::orderEdges()
int GEdgeCompound::minimumMeshSegments() const
{
- int N = 0;
- for (unsigned int i = 0; i < _compound.size(); i++)
- N +=_compound[i]->minimumMeshSegments();
+ // int N = 0;
+ // for (unsigned int i = 0; i < _compound.size(); i++)
+ // N +=_compound[i]->minimumMeshSegments();
return 3;
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 173ebadda8..346950716b 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -671,6 +671,7 @@ bool GFaceCompound::parametrize() const
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
_rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
+ double t0 = Cpu();
//_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(_rbf->getXYZ(), _rbf->getN(), Oper);
@@ -679,8 +680,11 @@ bool GFaceCompound::parametrize() const
//_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
_rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
-
- printStuff();
+ double t1 = Cpu();
+
+ printf("Time =%g \n", t1-t0);
+
+ //printStuff();
}
@@ -952,6 +956,9 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
}
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+ //EMI FIX
+ //if ((*it)->tag() == 3) _mapping = CONFORMAL;
+
if(!(*it)){
Msg::Error("Incorrect face in compound surface %d\n", tag);
Msg::Exit(1);
@@ -1066,7 +1073,6 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
dofManager<double> myAssembler(_lsys);
if(_type == UNITCIRCLE){
- printf("unit circle \n");
for(unsigned int i = 0; i < _ordered.size(); i++){
MVertex *v = _ordered[i];
const double theta = 2 * M_PI * _coords[i];
@@ -1369,7 +1375,6 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
if(!oct) parametrize();
-
if(trivial())
{
return (*(_compound.begin()))->curvatureMax(param);
@@ -1462,12 +1467,10 @@ GPoint GFaceCompound::point(double par1, double par2) const
if(!oct) parametrize();
-
-
double U,V;
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
- getTriangle (par1, par2, <, U,V);
+ getTriangle (par1, par2, <, U,V);
SPoint3 p(3, 3, 0);
if(!lt && _mapping != RBF){
GPoint gp (p.x(),p.y(),p.z(),this);
@@ -1566,9 +1569,9 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
if(!lt && _mapping != RBF)
return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
else if (!lt && _mapping == RBF){
- double x, y, z;
- SVector3 dXdu, dXdv ;
- bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
+ double x, y, z;
+ SVector3 dXdu, dXdv ;
+ bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
return Pair<SVector3, SVector3>(dXdu,dXdv);
}
@@ -1700,7 +1703,22 @@ void GFaceCompound::getTriangle(double u, double v,
double &_u, double &_v) const
{
double uv[3] = {u, v, 0};
- *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
+ // if (_mapping == RBF){
+ // std::list<void*> l;
+ // Octree_SearchAll(uv, oct, &l);
+ // if (l.size() > 1 || l.size() == 0){
+ // GFaceCompoundTriangle *gfct = NULL;
+ // *lt = gfct;
+ // return;
+ // }
+ // else{
+ // std::list<void*>::iterator it = l.begin();
+ // *lt = (GFaceCompoundTriangle*)(*it);
+ // }
+ // }
+
+ *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
+
// if(!(*lt)) {
// for(int i=0;i<nbT;i++){
// if(GFaceCompoundInEle (&_gfct[i],uv)){
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 7ad8dd3254..238ffac0cf 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -233,6 +233,7 @@ double OCCFace::curvatureMax(const SPoint2 ¶m) const
return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
}
+
double OCCFace::curvatures(const SPoint2 ¶m,
SVector3 *dirMax,
SVector3 *dirMin,
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 7ef8239de2..4ed94e0d98 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -501,7 +501,6 @@ double discreteEdge::curvature(double par) const
double c0, c1;
-
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() )
@@ -522,6 +521,7 @@ double discreteEdge::curvature(double par) const
double cv = (1-tLoc)*c0 + tLoc*c1;
+ //printf("curv edge =%g \n", cv);
return cv;
}
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 369d107f27..c439a0a63d 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -303,7 +303,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
double l2 = MAX_LC;
if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
l2 = LC_MVertex_PNTS(ge, U, V);
-
+
// lc from curvature
double l3 = MAX_LC;
if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
@@ -314,15 +314,16 @@ double BGM_MeshSize(GEntity *ge, double U, double V,
FieldManager *fields = ge->model()->getFields();
if(fields->background_field > 0){
Field *f = fields->get(fields->background_field);
- // printf("field %p %s %d %p\n",f,f->getName(),fields->size(), ge->model());
+ //printf("field %p %s %d %p\n",f,f->getName(),fields->size(), ge->model());
if(f) l4 = (*f)(X, Y, Z, ge);
+ //printf("X Y Z =%g %g %g L4=%g L3=%g L2=%g L1=%g\n", X, Y, Z, l4, l3, l2, l1);
}
// take the minimum, then constrain by lcMin and lcMax
double lc = std::min(std::min(std::min(l1, l2), l3), l4);
lc = std::max(lc, CTX::instance()->mesh.lcMin);
lc = std::min(lc, CTX::instance()->mesh.lcMax);
-
+
if(lc <= 0.){
Msg::Error("Wrong mesh element size lc = %g (lcmin = %g, lcmax = %g)",
lc, CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 733bc67a30..20e0e3752d 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -39,8 +39,10 @@ static double F_Lc(GEdge *ge, double t)
SVector3 der = ge->firstDer(t);
//const double d = norm(der);
//return d / lc_here;
+
SMetric3 metric(1. / (lc_here * lc_here));
double lSquared = dot(der, metric, der);
+
return sqrt(lSquared);
}
diff --git a/benchmarks/curvature/Torus.geo b/benchmarks/curvature/Torus.geo
index 5e6adbe6e2..b38ec2eabb 100644
--- a/benchmarks/curvature/Torus.geo
+++ b/benchmarks/curvature/Torus.geo
@@ -1,11 +1,11 @@
-//Mesh.CharacteristicLengthFactor=0.3;
-Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
+
+Mesh.CharacteristicLengthFactor = 0.6;
+Mesh.CharacteristicLengthFromCurvature = 1; //-clcurv
Mesh.CharacteristicLengthMin = 0.1; //-clmin
Mesh.CharacteristicLengthMax = 2.5; //-clmax
Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
-
Mesh.MinimumCirclePoints=15; //default=7
-Mesh.CharacteristicLengthExtendFromBoundary=0;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;
lc = 0.1;
Point(1) = {2.0,0.0,0.0,lc};
diff --git a/benchmarks/curvature/TorusRemesh.geo b/benchmarks/curvature/TorusRemesh.geo
index b1c1fb203e..53e0d06c9b 100644
--- a/benchmarks/curvature/TorusRemesh.geo
+++ b/benchmarks/curvature/TorusRemesh.geo
@@ -1,27 +1,26 @@
// ======================================================================
// The following settings were copied from benchmarks/step/linkrods.geo
// ======================================================================
-//Mesh.RemeshAlgorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal
-//Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
-Mesh.CharacteristicLengthFactor=0.3;
+
+Mesh.RemeshAlgorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal
+Mesh.RemeshParametrization= 1; //(0) harmonic (1) conformal
+Mesh.RemeshAlgorithm = 1;
+
+Mesh.CharacteristicLengthFactor=0.2;
+Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature=1; //-clcurv
Mesh.CharacteristicLengthMin = 0.2; //-clmin
-Mesh.CharacteristicLengthMax = 1.5; //-clmax
-Mesh.LcIntegrationPrecision=1.e-5; //-epslc1d
+Mesh.CharacteristicLengthMax = 3.0; //-clmax
+Mesh.LcIntegrationPrecision = 1.e-5; //-epslc1d
Mesh.MinimumCirclePoints=15; //default=7
Mesh.CharacteristicLengthExtendFromBoundary=0;
-// ======================================================================
-
//Merge "TorusRemeshedBAMG.stl";
-// Merge "SimpleTorus.msh";
+//Merge "SimpleTorus.msh";
Merge "TorusInput.stl";
-//CreateTopology;
-
-Mesh.RemeshAlgorithm=6; //1=MeshAdapt, 5=Delaunay, 6=Frontal
-Mesh.RemeshParametrization=1; //(0) harmonic (1) conformal
+CreateTopology;
Compound Surface(20) = {1};
//Compound Surface(100) = {15,19,23,27};
diff --git a/benchmarks/curvature/occtorus.brep b/benchmarks/curvature/occtorus.brep
new file mode 100644
index 0000000000..c97da515a0
--- /dev/null
+++ b/benchmarks/curvature/occtorus.brep
@@ -0,0 +1,64 @@
+
+CASCADE Topology V1, (c) Matra-Datavision
+Locations 0
+Curve2ds 4
+1 0 0 1 0
+1 0 6.2831853071795862 1 0
+1 6.2831853071795862 -0 0 1
+1 0 -0 0 1
+Curves 2
+2 -2.4492935982947064e-16 0 0 1 0 0 -0 0 1 0 -1 0 4
+2 0 7.3478807948841188e-16 3 0 1 -2.4492935982947064e-16 0 2.4492935982947064e-16 1 1 0 0 1
+Polygon3D 0
+PolygonOnTriangulations 0
+Surfaces 1
+5 0 0 0 1 0 0 -0 0 1 0 -1 0 3 1
+Triangulations 0
+
+TShapes 8
+Ve
+1e-07
+-2.44929359829471e-16 9.79717439317883e-16 4
+0 0
+
+0101101
+*
+Ed
+ 1e-07 1 1 0
+1 1 0 0 6.28318530717959
+3 1 2CN 1 0 0 6.28318530717959
+0
+
+0101000
++8 0 -8 0 *
+Ed
+ 1e-07 1 1 0
+1 2 0 0 6.28318530717959
+3 3 4CN 1 0 0 6.28318530717959
+0
+
+0101000
++8 0 -8 0 *
+Wi
+
+0101000
+-7 0 +6 0 +7 0 -6 0 *
+Fa
+0 1e-07 1 0
+
+0111000
++5 0 *
+Sh
+
+0101100
++4 0 *
+So
+
+0100100
++3 0 *
+Co
+
+1100000
++2 0 *
+
++1 0
\ No newline at end of file
diff --git a/benchmarks/curvature/torus.py b/benchmarks/curvature/torus.py
new file mode 100644
index 0000000000..701a4959f1
--- /dev/null
+++ b/benchmarks/curvature/torus.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+
+from gmshpy import *
+
+lc = 0.8
+GmshSetOption('Mesh', 'CharacteristicLengthFactor', lc)
+
+R = 1 ;
+g = GModel();
+g.setFactory('OpenCASCADE')
+
+#g.addBlock([-R,-R,-R],[R,R,R]);
+g.addTorus([0,0,0],[1,0,0], 3.0, 1.0);
+
+#g.mesh(2);
+#g.save("occtorus.msh");
+g.save("occtorus.brep");
+
+#FlGui.instance();
+#FlGui.run();
diff --git a/gmshpy/gmshGeo.i b/gmshpy/gmshGeo.i
index 21607ce3b4..2c4bb975d6 100644
--- a/gmshpy/gmshGeo.i
+++ b/gmshpy/gmshGeo.i
@@ -4,6 +4,7 @@
%include std_string.i
%include std_list.i
%include std_vector.i
+%import "gmshCommon.i"
%{
#include "GmshConfig.h"
--
GitLab