diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 7c61c3bbdbdd18a8eada74592879a3e30f0b026b..778c0128e1d23b30585b2d1d57905ef51840aad5 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -50,6 +50,7 @@ void Print_Usage(const char *name)
Msg::Direct(" bdf, p3d, cgns, med)");
Msg::Direct(" -bin Use binary format when available");
Msg::Direct(" -parametric Save vertices with their parametric coordinates");
+ Msg::Direct(" -numsubedges Set the number of subdivisions when displaying high order elements");
Msg::Direct(" -algo string Select mesh algorithm (de, del2d, frontal, iso, netgen, tetgen)");
Msg::Direct(" -smooth int Set number of mesh smoothing steps");
Msg::Direct(" -optimize[_netgen] Optimize quality of tetrahedral elements");
@@ -417,6 +418,13 @@ void Get_Options(int argc, char *argv[])
else
Msg::Fatal("Missing number");
}
+ else if(!strcmp(argv[i] + 1, "numsubedges")) {
+ i++;
+ if(argv[i] != NULL)
+ opt_mesh_num_sub_edges(0, GMSH_SET, atof(argv[i++]));
+ else
+ Msg::Fatal("Missing number");
+ }
else if(!strcmp(argv[i] + 1, "c1")) {
i++;
opt_mesh_c1(0, GMSH_SET, 1);
diff --git a/Common/Context.h b/Common/Context.h
index bbc82317292b17fdc7e6821a5226e1e518c15ecc..54e35733ca80e286b83e4d11cf7f090085ee6baf 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -155,6 +155,7 @@ class Context_T {
int points, lines, triangles, quadrangles, tetrahedra, hexahedra, prisms, pyramids;
int surfaces_edges, surfaces_faces, volumes_edges, volumes_faces;
int points_num, lines_num, surfaces_num, volumes_num;
+ int num_sub_edges;
double label_frequency;
int point_type; // flat or 3D
double point_size, line_width;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index ba972ab05302173aac732775f2093850b2a71813..5b08c7023d47d2d3dd75e1587ddbe4910aad620c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1013,6 +1013,8 @@ StringXNumber MeshOptions_Number[] = {
"Number of triangles in the current mesh (read-only)" },
{ F|O, "Normals" , opt_mesh_normals , 0.0 ,
"Display size of normal vectors (in pixels)" },
+ { F|O, "NumSubEdges" , opt_mesh_num_sub_edges , 2. ,
+ "Number of edge subdivisions when displaying high order elements" },
{ F|O, "Optimize" , opt_mesh_optimize , 0. ,
"Optimize the mesh to improve the quality of tetrahedral elements" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 0ffc03d2253aabdae8980768f503fa2986c89942..18410864462b8ae9642884f061295cd6363e2e2b 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4443,6 +4443,13 @@ double opt_mesh_normals(OPT_ARGS_NUM)
return CTX.mesh.normals;
}
+double opt_mesh_num_sub_edges(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.num_sub_edges =(int) val;
+ return CTX.mesh.num_sub_edges;
+}
+
double opt_mesh_tangents(OPT_ARGS_NUM)
{
if(action & GMSH_SET) {
diff --git a/Common/Options.h b/Common/Options.h
index d8f535d5812fe794692523ef5815f34fc8b334b6..b8e71c9b2f4119c537e9171083acb54d679584b9 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -421,6 +421,7 @@ double opt_mesh_optimize(OPT_ARGS_NUM);
double opt_mesh_optimize_netgen(OPT_ARGS_NUM);
double opt_mesh_refine_steps(OPT_ARGS_NUM);
double opt_mesh_normals(OPT_ARGS_NUM);
+double opt_mesh_num_sub_edges(OPT_ARGS_NUM);
double opt_mesh_tangents(OPT_ARGS_NUM);
double opt_mesh_explode(OPT_ARGS_NUM);
double opt_mesh_scaling_factor(OPT_ARGS_NUM);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e5146f327b768b0925b9ae97c924b0fca87177a4..c96dfd1c18cf0a010e7357a2ba5c3026d6be957f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -10,6 +10,37 @@
#include "Octree.h"
#include "gmshLinearSystemGmm.h"
+class gmshGradientBasedDiffusivity : public gmshFunction
+{
+ MElement *_current;
+ int _iComp;
+ mutable std::map<MVertex*,SPoint2> _coordinates;
+public:
+ gmshGradientBasedDiffusivity (std::map<MVertex*,SPoint2> &coordinates)
+ : _current (0), _iComp(-1),_coordinates(coordinates){}
+ void setCurrent (MElement *current){_current=current;}
+ void setComponent (int iComp){_iComp=iComp;}
+ virtual double operator () (double x, double y, double z) const {
+ if (_iComp == -1)return 1.0;
+ double xyz[3]={x,y,z},uvw[3];
+ _current->xyz2uvw(xyz,uvw);
+ double value1[256];
+ double value2[256];
+ for (int i=0;i<_current->getNumVertices();i++){
+ const SPoint2 p = _coordinates[_current->getVertex(i)];
+ value1[i] = p[0];
+ value2[i] = p[1];
+ }
+ double g1[3],g2[3];
+ _current->interpolateGrad(value1,uvw[0],uvw[1],uvw[2],g1);
+ _current->interpolateGrad(value2,uvw[0],uvw[1],uvw[2],g2);
+ // printf("%g %g %g -- %g %g %g \n",g[0],g[1],g[2],value[0],value[1],value[2]);
+ return sqrt(g1[0]*g1[0]+g1[1]*g1[1]+g1[2]*g1[2] +
+ g2[0]*g2[0]+g2[1]*g2[1]+g2[2]*g2[2]);
+ }
+};
+
+
static void fixEdgeToValue (GEdge *ed,
double value,
gmshAssembler &myAssembler)
@@ -23,21 +54,16 @@ static void fixEdgeToValue (GEdge *ed,
}
}
-static void fixEdgeToValueX (GEdge *ed,
- gmshAssembler &myAssembler)
-{
- for (unsigned int i = 0 ; i < ed->lines.size(); i++){
- myAssembler.fixVertex(ed->lines[i]->getVertex(0), 0, 1, ed->lines[i]->getVertex(0)->x());
- myAssembler.fixVertex(ed->lines[i]->getVertex(1), 0, 1, ed->lines[i]->getVertex(1)->x());
- }
-}
-
-static void fixEdgeToValueY (GEdge *ed,
- gmshAssembler &myAssembler)
+void GFaceCompound::parametrize() const
{
- for (unsigned int i = 0 ; i < ed->lines.size(); i++){
- myAssembler.fixVertex(ed->lines[i]->getVertex(0), 0, 1, ed->lines[i]->getVertex(0)->y());
- myAssembler.fixVertex(ed->lines[i]->getVertex(1), 0, 1, ed->lines[i]->getVertex(1)->y());
+ if (!oct){
+ coordinates.clear();
+ parametrize(0,0);
+ parametrize(1,0);
+ // parametrize(0,1);
+ // parametrize(1,1);
+ computeNormals () ;
+ buildOct();
}
}
@@ -47,7 +73,6 @@ void GFaceCompound::getBoundingEdges ()
if (_U0.size() && !_U1.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
for ( ; it != _U0.end() ; ++it){
- // printf("edge %d\n",(*it)->tag());
l_edges.push_back(*it);
(*it)->addFace(this);
}
@@ -55,16 +80,26 @@ void GFaceCompound::getBoundingEdges ()
}
std::set<GEdge*> _unique;
+ std::multiset<GEdge*> _touched;
std::list<GFace*> :: iterator it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
std::list<GEdge*> ed = (*it)->edges();
std::list<GEdge*> :: iterator ite = ed.begin();
for ( ; ite != ed.end() ; ++ite){
- std::set<GEdge*>::iterator itf = _unique.find(*ite);
- if (itf == _unique.end())_unique.insert(*ite);
- else _unique.erase(itf);
+ _touched.insert(*ite);
}
}
+ it = _compound.begin();
+ for ( ; it != _compound.end() ; ++it){
+ std::list<GEdge*> ed = (*it)->edges();
+ std::list<GEdge*> :: iterator ite = ed.begin();
+ for ( ; ite != ed.end() ; ++ite){
+ if (!(*ite)->degenerate(0) && _touched.count(*ite) == 1)_unique.insert(*ite);
+ }
+ }
+
+
+
std::set<GEdge*>::iterator itf = _unique.begin();
for ( ; itf != _unique.end() ; ++itf){
l_edges.push_back(*itf);
@@ -80,7 +115,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &V1) :
GFace(m,tag),_compound(compound),_U0(U0),_U1(U1),_V0(V0),_V1(V1), oct(0)
{
- //printf("%d %d %d %d \n",_U0.size(),_U1.size(),_V0.size(),_V1.size());
getBoundingEdges();
}
@@ -92,23 +126,104 @@ GFaceCompound::~GFaceCompound()
}
}
-void GFaceCompound::parametrize (bool _isU) const
+static bool orderVertices (const std::list<GEdge*> &e ,
+ std::vector<MVertex*> &l,
+ std::vector<double> &coord){
+ l.clear();
+ coord.clear();
+ std::list<GEdge*> :: const_iterator it = e.begin();
+ std::list<MLine*> temp;
+ double tot_length = 0;
+ for ( ; it != e.end() ; ++it ){
+ // printf("GLine %d\n",(*it)->tag());
+ for ( int i = 0 ; i < (*it)->lines.size(); i++ ){
+ temp.push_back((*it)->lines[i]);
+ MVertex *v0 = (*it)->lines[i]->getVertex(0);
+ MVertex *v1 = (*it)->lines[i]->getVertex(1);
+ // printf("UNSORTED %d -> %d \n",v0->getNum(),v1->getNum());
+ const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()) );
+ tot_length += length;
+ }
+ }
+
+ MVertex *first_v = (*temp.begin())->getVertex(0);
+ MVertex *current_v = (*temp.begin())->getVertex(1);
+
+ l.push_back(first_v);
+ coord.push_back(0.0);
+ temp.erase(temp.begin());
+
+
+ // printf("SORTED %d -> %d \n",first_v->getNum(),current_v->getNum());
+
+ while (temp.size()){
+ bool found = 0;
+ for (std::list<MLine*>::iterator itl = temp.begin() ; itl !=temp.end() ; ++itl){
+ MLine *ll = *itl;
+ MVertex *v0 = ll->getVertex(0);
+ MVertex *v1 = ll->getVertex(1);
+ if (v0 == current_v){
+ found = true;
+ l.push_back(current_v);
+ current_v = v1;
+ temp.erase(itl);
+ const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()) );
+ coord.push_back(coord[coord.size()-1] + length/tot_length);
+ // printf("SORTED %d -> %d (%d,%g)\n",v0->getNum(),v1->getNum(),temp.size(),coord[coord.size()-1]);
+ break;
+ }
+ else if (v1 == current_v){
+ found = true;
+ l.push_back(current_v);
+ current_v = v0;
+ temp.erase(itl);
+ const double length = sqrt ((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()) );
+ coord.push_back(coord[coord.size()-1] + length/tot_length);
+ // printf("SORTED %d -> %d (%d,%g)\n",v1->getNum(),v0->getNum(),temp.size(),coord[coord.size()-1]);
+ break;
+ }
+ }
+ if (!found)return false;
+ }
+ // printf(" DONE\n");
+
+ return true;
+}
+
+
+void GFaceCompound::parametrize (bool _isU, int ITER) const
{
+ Msg::Info("Parametrizing Surface %d coordinate %d (ITER %d)", tag(), _isU, ITER);
+
gmshLinearSystemGmm lsys;
+ lsys.setPrec(1.e-10);
+ lsys.setNoisy(2);
gmshAssembler myAssembler(&lsys);
-
- if (!_U0.size()){
- std::list<GEdge*> :: const_iterator it = l_edges.begin();
- for ( ; it != l_edges.end() ; ++it){
- if (_isU)fixEdgeToValueX (*it, myAssembler);
- else fixEdgeToValueY (*it, myAssembler);
- }
- }
- else if (!_U1.size()){
- std::list<GEdge*> :: const_iterator it = _U0.begin();
- for ( ; it != _U0.end() ; ++it){
- if (_isU)fixEdgeToValueX (*it, myAssembler);
- else fixEdgeToValueY (*it, myAssembler);
+ gmshGradientBasedDiffusivity diffusivity (coordinates);
+ if (ITER > 0) diffusivity.setComponent(_isU);
+
+ if (!_U1.size()){
+ // maps the boundary onto a circle
+ std::vector<MVertex*> ordered;
+ std::vector<double> coords;
+ Msg::Info("%d edges on the contour", l_edges.size());
+ for (std::list<GEdge*>::const_iterator it = l_edges.begin();
+ it !=l_edges.end();++it)printf("%d ",(*it)->tag());
+ printf("\n");
+ bool success = orderVertices (l_edges, ordered, coords);
+ if (!success)throw;
+ for (int i=0; i<ordered.size();i++){
+ MVertex *v = ordered[i];
+ const double theta = 2*M_PI*coords[i];
+ // printf("fixing %d to %g\n",v->getIndex(),theta);
+ if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta));
+ else myAssembler.fixVertex(v, 0, 1, sin(theta));
}
}
else{
@@ -133,7 +248,7 @@ void GFaceCompound::parametrize (bool _isU) const
}
}
}
-
+
std::list<GFace*> :: const_iterator it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -145,13 +260,19 @@ void GFaceCompound::parametrize (bool _isU) const
}
// printf("%d %d %d %d\n",_U0.size(),_U1.size(),_V0.size(),_V1.size());
// printf("creating term %d dofs numbered %d fixed\n", myAssembler.sizeOfR(),myAssembler.sizeOfF());
- gmshLaplaceTerm laplace (model(), 1.0, 1);
+ gmshLaplaceTerm laplace (model(), &diffusivity, 1);
it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
- laplace.addToMatrix (myAssembler, *it);
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ diffusivity.setCurrent(t);
+ laplace.addToMatrix (myAssembler, t);
+ }
}
+ // printf("assembly done\n");
lsys.systemSolve();
-
+ // printf("system souleved\n");
+
it = _compound.begin();
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -171,23 +292,77 @@ void GFaceCompound::parametrize (bool _isU) const
uu[j] = itf->second[0];
vv[j] = itf->second[1];
}
- }
+ }
}
}
+}
+
+
+void GFaceCompound::computeNormals () const
+{
+ _normals.clear();
+ std::list<GFace*>::const_iterator it = _compound.begin();
+ double J[3][3];
+ for ( ; it != _compound.end() ; ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ double det = t->getJacobian(0,0,0,J);
+ // SVector3 n (J[2][0],J[2][1],J[2][2]);
+ SVector3 d1(J[0][0], J[0][1], J[0][2]);
+ SVector3 d2(J[1][0], J[1][1], J[1][2]);
+ SVector3 n = crossprod(d1, d2);
+ n.normalize();
+ for (int j=0;j<3;j++){
+ std::map<MVertex*,SVector3>::iterator itn = _normals.find(t->getVertex(j));
+ if (itn == _normals.end())_normals[t->getVertex(j)] = n;
+ else itn->second += n;
+ }
+ }
+ }
+ std::map<MVertex*,SVector3>::iterator itn = _normals.begin();
+ for ( ; itn != _normals.end() ; ++itn) itn->second.normalize();
+}
+
+double GFaceCompound::curvature(const SPoint2 ¶m) const
+{
+ parametrize();
+ double U,V;
+ GFaceCompoundTriangle *lt;
+ getTriangle (param.x(),param.y(), <, U,V);
+ if (!lt){
+ return 0.0;
+ }
+ return curvature(lt->t);
+}
+double GFaceCompound::curvature(MTriangle *t) const
+{
+ SVector3 n1 = _normals[t->getVertex(0)];
+ SVector3 n2 = _normals[t->getVertex(1)];
+ SVector3 n3 = _normals[t->getVertex(2)];
+ double val[9] = {n1.x(),n2.x(),n3.x(),
+ n1.y(),n2.y(),n3.y(),
+ n1.z(),n2.z(),n3.z()};
+ return fabs(t->interpolateDiv (val,0,0,0.));
}
+
GPoint GFaceCompound::point(double par1, double par2) const
{
parametrize();
- double U, V;
- MTriangle *t;
- getTriangle (par1, par2, &t, U,V);
- SPoint3 p(0, 0, 0);
- if (!t) return GPoint(p.x(), p.y(), p.z(), this);
- t->pnt(U, V, 0, p);
- double par[2] = {par1, par2};
- return GPoint(p.x(), p.y(), p.z(), this, par);
+ double U,V;
+ GFaceCompoundTriangle *lt;
+ getTriangle (par1, par2, <, U,V);
+ SPoint3 p(0,0,0);
+ if (!lt){
+ printf("coucou\n");
+ return GPoint(p.x(),p.y(),p.z(),this);
+ }
+ p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
+ // lt->t->pnt(U,V,0,p);
+ double par[2] = {par1,par2};
+ // printf("coucou %g %g -> %g %g %g\n",par1,par2,p.x(),p.y(),p.z());
+ return GPoint(p.x(),p.y(),p.z(),this,par);
}
/*
@@ -197,13 +372,53 @@ GPoint GFaceCompound::point(double par1, double par2) const
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
parametrize();
- double U, V;
- MTriangle *t;
- getTriangle(param.x(), param.y(), &t, U,V);
- double jac[3][3];
- t->getJacobian(U, V, 0, jac);
- return Pair<SVector3, SVector3>(SVector3(jac[0][0], jac[0][1], jac[0][2]),
- SVector3(jac[1][0], jac[1][1], jac[1][2]));
+ double U,V,UDU,UDV,VDU,VDV;
+ GFaceCompoundTriangle *lt;
+ // GFaceCompoundTriangle *ltdu;
+ // GFaceCompoundTriangle *ltdv;
+ getTriangle (param.x(), param.y(), <, U,V);
+ // getTriangle (param.x()+1.e-4, param.y(), <du, UDU,VDU);
+ // getTriangle (param.x(), param.y()+1.e-4, <dv, UDV,VDV);
+
+ if (!lt)
+ return Pair<SVector3, SVector3>(SVector3(1,0,0),SVector3(0,1,0));
+
+ double mat[2][2] = {{lt->p2.x()-lt->p1.x(),lt->p3.x()-lt->p1.x()},
+ {lt->p2.y()-lt->p1.y(),lt->p3.y()-lt->p1.y()}};
+ double inv[2][2];
+ inv2x2(mat,inv);
+
+ // MVertex *v0 = lt->t->getVertex(0);
+ // MVertex *v1 = lt->t->getVertex(1);
+ // MVertex *v2 = lt->t->getVertex(2);
+
+ // SPoint3 p(0,0,0);
+ // SPoint3 pdu(0,0,0);
+ // SPoint3 pdv(0,0,0);
+
+ // lt->t->pnt(U,V,0,p);
+ // ltdu->t->pnt(UDU,VDU,0,pdu);
+ // ltdv->t->pnt(UDV,VDV,0,pdv);
+
+
+ SVector3 dXdxi (lt->v2 - lt->v1);
+ SVector3 dXdeta (lt->v3 - lt->v1);
+
+ // return Pair<SVector3, SVector3>(dXdxi*inv[0][0]+dXdeta*inv[1][0],
+ // dXdxi*inv[0][1]+dXdeta*inv[1][1]);
+
+ SVector3 dXdu (dXdxi*inv[0][0]+dXdeta*inv[1][0]);
+ SVector3 dXdv (dXdxi*(inv[0][1])+dXdeta*(inv[1][1]));
+
+ // SVector3 dXduFD (SVector3(pdu-p)*1.e4);
+ // SVector3 dXdvFD (SVector3(pdv-p)*1.e4);
+
+// printf("FD %g %g %g / %g %g %g AN %g %g %g / %g %g %g \n",
+// dXdu.x(),dXdu.y(),dXdu.z(),
+// dXdv.x(),dXdv.y(),dXdv.z(),
+// dXduFD.x(),dXduFD.y(),dXduFD.z(),
+// dXdvFD.x(),dXdvFD.y(),dXdvFD.z());
+ return Pair<SVector3, SVector3>(dXdu,dXdv);
}
static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
@@ -215,6 +430,13 @@ static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
mmax[1] = std::max(std::max(t->p1.y(), t->p2.y()), t->p3.y());
mmin[2] = -1;
mmax[2] = 1;
+
+ const double dx = mmax[0] - mmin[0];
+ const double dy = mmax[1] - mmin[1];
+ mmin[0] -= .02*dx;
+ mmin[1] -= .02*dy;
+ mmax[0] += .02*dx;
+ mmax[1] += .02*dy;
}
static void GFaceCompoundCentroid(void *a, double*c)
@@ -229,7 +451,7 @@ static int GFaceCompoundInEle(void *a, double*c)
{
GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a;
double M[2][2],R[2],X[2];
- const double eps = 1.e-6;
+ const double eps = 1.e-8;
const SPoint2 p0 = t->p1;
const SPoint2 p1 = t->p2;
const SPoint2 p2 = t->p3;
@@ -246,22 +468,20 @@ static int GFaceCompoundInEle(void *a, double*c)
return 0;
}
-void GFaceCompound::getTriangle(double u, double v,
- MTriangle **t,
- double &_u, double &_v) const
+void GFaceCompound::getTriangle (double u, double v,
+ GFaceCompoundTriangle **lt,
+ double &_u, double &_v) const
{
parametrize();
double uv[3] = {u,v,0};
- GFaceCompoundTriangle *tt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
- if (!tt){*t = 0; return;}
+ *lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
+ if (!(*lt)){return;}
- *t = tt->t;
double M[2][2],X[2],R[2];
-// const double eps = 1.e-6;
- const SPoint2 p0 = tt->p1;
- const SPoint2 p1 = tt->p2;
- const SPoint2 p2 = tt->p3;
+ const SPoint2 p0 = (*lt)->p1;
+ const SPoint2 p1 = (*lt)->p2;
+ const SPoint2 p2 = (*lt)->p3;
M[0][0] = p1.x() - p0.x();
M[0][1] = p2.x() - p0.x();
M[1][0] = p1.y() - p0.y();
@@ -272,40 +492,6 @@ void GFaceCompound::getTriangle(double u, double v,
_u = X[0];
_v = X[1];
return;
-
-// std::list<GFace*> :: const_iterator it = _compound.begin();
-// for ( ; it != _compound.end() ; ++it){
-// for ( int i=0;i<(*it)->triangles.size() ; ++i){
-// *t = (*it)->triangles[i];
-// std::map<MVertex*,SPoint2>::const_iterator it0 =
-// coordinates.find((*t)->getVertex(0));
-// std::map<MVertex*,SPoint2>::const_iterator it1 =
-// coordinates.find((*t)->getVertex(1));
-// std::map<MVertex*,SPoint2>::const_iterator it2 =
-// coordinates.find((*t)->getVertex(2));
-// const SPoint2 p0 = it0->second;
-// const SPoint2 p1 = it1->second;
-// const SPoint2 p2 = it2->second;
-// M[0][0] = p1.x() - p0.x();
-// M[0][1] = p2.x() - p0.x();
-// M[1][0] = p1.y() - p0.y();
-// M[1][1] = p2.y() - p0.y();
-// R[0] = (u - p0.x());
-// R[1] = (v - p0.y());
-// sys2x2(M,R,X);
-// // printf("T (%g %g,%g %g,%g %g) %g %g vs %g %g\n",p0.x(),p0.y(),p1.x(),p1.y(),p2.x(),p2.y(),u,v,X[0],X[1]);
-// if (X[0] > -eps &&
-// X[1] > -eps &&
-// 1.-X[0]-X[1] > -eps){
-// _u = X[0];
-// _v = X[1];
-// // printf("OK -- %g %g\n",_u,_v);
-// return;
-// }
-// }
-// }
-// printf("NOT OK -- %g %g\n",u,v);
-// *t = 0;
}
void GFaceCompound::buildOct() const
@@ -313,6 +499,7 @@ void GFaceCompound::buildOct() const
SBoundingBox3d bb;
int count = 0;
std::list<GFace*> :: const_iterator it = _compound.begin();
+ // printf("building octree %d coords \n",coordinates.size());
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
@@ -340,12 +527,14 @@ void GFaceCompound::buildOct() const
FILE * uvz = fopen ("UVZ.pos","w");
FILE * xyzu = fopen("XYZU.pos","w");
FILE * xyzv = fopen("XYZV.pos","w");
+ FILE * xyzc = fopen("XYZC.pos","w");
fprintf(uvx,"View \"\"{\n");
fprintf(uvy,"View \"\"{\n");
fprintf(uvz,"View \"\"{\n");
fprintf(xyzu,"View \"\"{\n");
fprintf(xyzv,"View \"\"{\n");
+ fprintf(xyzc,"View \"\"{\n");
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -360,6 +549,9 @@ void GFaceCompound::buildOct() const
_gfct[count].p2 = it1->second;
_gfct[count].p3 = it2->second;
_gfct[count].t = t;
+ _gfct[count].v1 = SPoint3(t->getVertex(0)->x(),t->getVertex(0)->y(),t->getVertex(0)->z());
+ _gfct[count].v2 = SPoint3(t->getVertex(1)->x(),t->getVertex(1)->y(),t->getVertex(1)->z());
+ _gfct[count].v3 = SPoint3(t->getVertex(2)->x(),t->getVertex(2)->y(),t->getVertex(2)->z());
fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(),
t->getVertex(0)->y(),
@@ -382,6 +574,18 @@ void GFaceCompound::buildOct() const
t->getVertex(2)->y(),
t->getVertex(2)->z(),
it0->second.y(),it1->second.y(),it2->second.y());
+ const double K = fabs(curvature (t));
+ fprintf(xyzc,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ t->getVertex(0)->x(),
+ t->getVertex(0)->y(),
+ t->getVertex(0)->z(),
+ t->getVertex(1)->x(),
+ t->getVertex(1)->y(),
+ t->getVertex(1)->z(),
+ t->getVertex(2)->x(),
+ t->getVertex(2)->y(),
+ t->getVertex(2)->z(),K,K,K);
+
fprintf(uvx,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(),it0->second.y(),0.0,
it1->second.x(),it1->second.y(),0.0,
@@ -412,6 +616,8 @@ void GFaceCompound::buildOct() const
fclose(xyzu);
fprintf(xyzv,"};\n");
fclose(xyzv);
+ fprintf(xyzc,"};\n");
+ fclose(xyzc);
Octree_Arrange(oct);
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 7bc2698bb7b165d5e752e1d3b92829a20c95b876..d7db76c451ccb6a7f24e086ac39275ceadf02a8d 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -28,6 +28,7 @@ generator of gmsh!
typedef struct {
SPoint2 p1, p2, p3;
+ SPoint3 v1, v2, v3;
MTriangle *t;
} GFaceCompoundTriangle;
@@ -40,19 +41,15 @@ class GFaceCompound : public GFace {
mutable GFaceCompoundTriangle *_gfct;
mutable Octree *oct;
mutable std::map<MVertex*,SPoint2> coordinates;
+ mutable std::map<MVertex*,SVector3> _normals;
void buildOct() const ;
- void parametrize(bool) const ;
- void parametrize() const
- {
- if (!oct){
- parametrize(0);
- parametrize(1);
- buildOct();
- }
- }
+ void parametrize(bool,int) const ;
+ void parametrize() const ;
+ void computeNormals () const;
void getBoundingEdges();
- void getTriangle(double u, double v, MTriangle **t, double &_u, double &_v) const;
- public:
+ void getTriangle(double u, double v, GFaceCompoundTriangle **lt, double &_u, double &_v) const;
+ virtual double curvature(MTriangle *t) const;
+public:
GFaceCompound(GModel *m, int tag,
std::list<GFace*> &compound,
std::list<GEdge*> &U0,
@@ -68,6 +65,8 @@ class GFaceCompound : public GFace {
void * getNativePtr() const { return 0; }
SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; }
virtual bool buildRepresentationCross(){ return false; }
+ virtual double curvature(const SPoint2 ¶m) const;
+
};
#endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d37ff5779fe6da6026addac05b023789639559da..7cb9df82d6fe4d45179f160102d1e733e37a2577 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -813,20 +813,20 @@ const gmshFunctionSpace* MTriangle::getFunctionSpace(int o) const
return 0;
}
-#define NUM_SUB_EDGES (_order)
+#define NUM_SUB_EDGES (CTX.mesh.num_sub_edges)
int MTriangleN::getNumEdgesRep(){ return 3 * NUM_SUB_EDGES; }
+int MTriangle6::getNumEdgesRep(){ return 3 * NUM_SUB_EDGES; }
-void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+static void _myGetEdgeRep(MTriangle *t, int num, double *x, double *y, double *z, SVector3 *n,
+ int numSubEdges)
{
- int numSubEdges = NUM_SUB_EDGES;
-
- n[0] = n[1] = n[2] = getFace(0).normal();
+ n[0] = n[1] = n[2] = t->getFace(0).normal();
if (num < numSubEdges){
SPoint3 pnt1, pnt2;
- pnt((double)num / numSubEdges, 0., 0.,pnt1);
- pnt((double)(num + 1) / numSubEdges, 0., 0, pnt2);
+ t->pnt((double)num / numSubEdges, 0., 0.,pnt1);
+ t->pnt((double)(num + 1) / numSubEdges, 0., 0, pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -835,8 +835,8 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
if (num < 2 * numSubEdges){
SPoint3 pnt1, pnt2;
num -= numSubEdges;
- pnt(1. - (double)num / numSubEdges, (double)num / numSubEdges, 0, pnt1);
- pnt(1. - (double)(num + 1) / numSubEdges, (double)(num + 1) / numSubEdges, 0, pnt2);
+ t->pnt(1. - (double)num / numSubEdges, (double)num / numSubEdges, 0, pnt1);
+ t->pnt(1. - (double)(num + 1) / numSubEdges, (double)(num + 1) / numSubEdges, 0, pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -845,19 +845,27 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
{
SPoint3 pnt1, pnt2;
num -= 2 * numSubEdges;
- pnt(0, (double)num / numSubEdges, 0,pnt1);
- pnt(0, (double)(num + 1) / numSubEdges, 0,pnt2);
+ t->pnt(0, (double)num / numSubEdges, 0,pnt1);
+ t->pnt(0, (double)(num + 1) / numSubEdges, 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
}
}
+void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n){
+ _myGetEdgeRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n){
+ _myGetEdgeRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+
+int MTriangle6::getNumFacesRep(){ return NUM_SUB_EDGES * NUM_SUB_EDGES; }
int MTriangleN::getNumFacesRep(){ return NUM_SUB_EDGES * NUM_SUB_EDGES; }
-void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+static void _myGetFaceRep(MTriangle *t, int num, double *x, double *y, double *z, SVector3 *n,
+ int numSubEdges)
{
- int numSubEdges = NUM_SUB_EDGES;
// on the first layer, we have (numSubEdges-1) * 2 + 1 triangles
// on the second layer, we have (numSubEdges-2) * 2 + 1 triangles
@@ -879,20 +887,20 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
SPoint3 pnt1, pnt2, pnt3;
double J1[3][3], J2[3][3], J3[3][3];
if (ix % 2 == 0){
- pnt(ix / 2 * d, iy * d, 0, pnt1);
- pnt((ix / 2 + 1) * d, iy * d, 0, pnt2);
- pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
- getJacobian(ix / 2 * d, iy * d, 0, J1);
- getJacobian((ix / 2 + 1) * d, iy * d, 0, J2);
- getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
+ t->pnt(ix / 2 * d, iy * d, 0, pnt1);
+ t->pnt((ix / 2 + 1) * d, iy * d, 0, pnt2);
+ t->pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
+ t->getJacobian(ix / 2 * d, iy * d, 0, J1);
+ t->getJacobian((ix / 2 + 1) * d, iy * d, 0, J2);
+ t->getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
}
else{
- pnt((ix / 2 + 1) * d, iy * d, 0, pnt1);
- pnt((ix / 2 + 1) * d, (iy + 1) * d, 0, pnt2);
- pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
- getJacobian((ix / 2 + 1) * d, iy * d, 0, J1);
- getJacobian((ix / 2 + 1) * d, (iy + 1) * d, 0, J2);
- getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
+ t->pnt((ix / 2 + 1) * d, iy * d, 0, pnt1);
+ t->pnt((ix / 2 + 1) * d, (iy + 1) * d, 0, pnt2);
+ t->pnt(ix / 2 * d, (iy + 1) * d, 0, pnt3);
+ t->getJacobian((ix / 2 + 1) * d, iy * d, 0, J1);
+ t->getJacobian((ix / 2 + 1) * d, (iy + 1) * d, 0, J2);
+ t->getJacobian(ix / 2 * d, (iy + 1) * d, 0, J3);
}
{
SVector3 d1(J1[0][0], J1[0][1], J1[0][2]);
@@ -918,6 +926,14 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
}
+void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){
+ _myGetFaceRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+void MTriangle6::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n){
+ _myGetFaceRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+
+
void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
#if !defined(HAVE_GMSH_EMBEDDED)
@@ -1036,12 +1052,11 @@ const gmshFunctionSpace* MTetrahedron::getFunctionSpace(int o) const
return 0;
}
+int MTetrahedron10::getNumEdgesRep(){ return 6 * NUM_SUB_EDGES; }
int MTetrahedronN::getNumEdgesRep(){ return 6 * NUM_SUB_EDGES; }
-void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+static void _myGetEdgeRep(MTetrahedron *tet, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges)
{
- int numSubEdges = NUM_SUB_EDGES;
-
static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
static int ed [6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
int iEdge = num / numSubEdges;
@@ -1049,30 +1064,39 @@ void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector
int iVertex1 = ed [iEdge][0];
int iVertex2 = ed [iEdge][1];
- double t1 = (double)iSubEdge / (double)numSubEdges;
- double u1 = pp[iVertex1][0] * (1. - t1) + pp[iVertex2][0] * t1;
- double v1 = pp[iVertex1][1] * (1. - t1) + pp[iVertex2][1] * t1;
- double w1 = pp[iVertex1][2] * (1. - t1) + pp[iVertex2][2] * t1;
+ double t1 = (double) iSubEdge / (double) numSubEdges;
+ double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
- double t2 = (double)(iSubEdge + 1) / (double)numSubEdges;
- double u2 = pp[iVertex1][0] * (1. - t2) + pp[iVertex2][0] * t2;
- double v2 = pp[iVertex1][1] * (1. - t2) + pp[iVertex2][1] * t2;
- double w2 = pp[iVertex1][2] * (1. - t2) + pp[iVertex2][2] * t2;
+ double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+ double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
SPoint3 pnt1, pnt2;
- pnt(u1, v1, w1, pnt1);
- pnt(u2, v2, w2, pnt2);
+ tet->pnt(u1,v1,w1,pnt1);
+ tet->pnt(u2,v2,w2,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
}
-int MTetrahedronN::getNumFacesRep(){ return 4 * NUM_SUB_EDGES * NUM_SUB_EDGES; }
+void MTetrahedron10::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetEdgeRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
-void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
- int numSubEdges = NUM_SUB_EDGES;
+ _myGetEdgeRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+int MTetrahedronN::getNumFacesRep(){ return 4 * NUM_SUB_EDGES * NUM_SUB_EDGES; }
+int MTetrahedron10::getNumFacesRep(){ return 4 * NUM_SUB_EDGES * NUM_SUB_EDGES; }
+
+void _myGetFaceRep (MTetrahedron *tet, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges )
+{
static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
static int fak [4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}};
int iFace = num / (numSubEdges * numSubEdges);
@@ -1134,9 +1158,9 @@ void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector
double W2 = pp[iVertex1][2] * (1.-u2-v2) + pp[iVertex2][2] * u2 + pp[iVertex3][2] * v2;
double W3 = pp[iVertex1][2] * (1.-u3-v3) + pp[iVertex2][2] * u3 + pp[iVertex3][2] * v3;
- pnt(U1,V1,W1,pnt1);
- pnt(U2,V2,W2,pnt2);
- pnt(U3,V3,W3,pnt3);
+ tet->pnt(U1,V1,W1,pnt1);
+ tet->pnt(U2,V2,W2,pnt2);
+ tet->pnt(U3,V3,W3,pnt3);
x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
@@ -1173,6 +1197,15 @@ void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector
}
}
+void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
+
+void MTetrahedron10::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this,num,x,y,z,n,NUM_SUB_EDGES);
+}
void MTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 97a4543c65b88bf1be960cea34c0ea65b168b2e4..4eba7945cda37b4411a8f455188f106703168267 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -641,31 +641,16 @@ class MTriangle6 : public MTriangle {
return getVertex(map[num]);
}
virtual int getNumEdgeVertices() const { return 3; }
- virtual int getNumEdgesRep(){ return 6; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[6][2] = {
- {0, 3}, {3, 1},
- {1, 4}, {4, 2},
- {2, 5}, {5, 0}
- };
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MTriangle::_getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 4; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[4][3] = {
- {0, 3, 5}, {1, 4, 3}, {2, 5, 4}, {3, 4, 5}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+ virtual int getNumFacesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(6);
@@ -1306,38 +1291,19 @@ class MTetrahedron10 : public MTetrahedron {
return getVertex(map[num]);
}
virtual int getNumEdgeVertices() const { return 6; }
- virtual int getNumEdgesRep(){ return 12; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[12][2] = {
- {0, 4}, {4, 1},
- {1, 5}, {5, 2},
- {2, 6}, {6, 0},
- {3, 7}, {7, 0},
- {3, 8}, {8, 2},
- {3, 9}, {9, 1}
- };
- static const int f[12] = {0, 0, 0, 1, 2, 3};
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
- }
+
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumEdgesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
+ virtual int getNumFacesRep();
+
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MTetrahedron::_getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 16; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[16][3] = {
- {0, 6, 4}, {2, 5, 6}, {1, 4, 5}, {6, 5, 4},
- {0, 4, 7}, {1, 9, 4}, {3, 7, 9}, {4, 9, 7},
- {0, 7, 6}, {3, 8, 7}, {2, 6, 8}, {7, 8, 6},
- {3, 9, 8}, {1, 5, 9}, {2, 8, 5}, {9, 5, 8}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(6);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 50214ebd5fec65033000adfdc4e8bb55083f57f2..9659fbcfa012cc4bd324add8fdf1bc228042fef4 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -30,6 +30,8 @@ void MVertex::writeMSH(FILE *fp, bool binary, bool saveParametric, double scalin
{
if(_index < 0) return; // negative index vertices are never saved
+ // printf("tag = %d index %d pos = %g %g %g\n",onWhat()->tag(),_index,x(),y(),z());
+
int myDim = 0, myTag = 0;
if(saveParametric){
if(onWhat()){
@@ -294,7 +296,8 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
{
- if (gf->geomType() == GEntity::CompoundSurface){
+ if (gf->geomType() == GEntity::CompoundSurface &&
+ v->onWhat()->dim() < 2){
GFaceCompound *gfc = (GFaceCompound*) gf;
param = gfc->getCoordinates(v);
return true;
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index ca27b04d5aa582485f5983cc2d0ffa6a23805ed6..0ef33b11566259f4f1f95ba8971f2154879e74d9 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -119,13 +119,17 @@ static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
GEdge *ged = (GEdge *)ge;
GVertex *v1 = ged->getBeginVertex();
GVertex *v2 = ged->getEndVertex();
- Range<double> range = ged->parBounds(0);
- double a = (U - range.low()) / (range.high() - range.low());
- double lc = (1 - a) * v1->prescribedMeshSizeAtVertex() +
- (a) * v2->prescribedMeshSizeAtVertex() ;
- // FIXME we might want to remove this to make all lc treatment consistent
- if(lc >= MAX_LC) return CTX.lc / 10.;
- return lc;
+ if (v1 && v2){
+ Range<double> range = ged->parBounds(0);
+ double a = (U - range.low()) / (range.high() - range.low());
+ double lc = (1 - a) * v1->prescribedMeshSizeAtVertex() +
+ (a) * v2->prescribedMeshSizeAtVertex() ;
+ // FIXME we might want to remove this to make all lc treatment consistent
+ if(lc >= MAX_LC) return CTX.lc / 10.;
+ return lc;
+ }
+ else
+ return MAX_LC;
}
default:
return MAX_LC;
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 1ed7bad20fc74689f946d24731e70e87b1fdfe81..730f5f2b02d6169510fab47630b65905fabe3cdb 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -407,7 +407,7 @@ static void Mesh2D(GModel *m)
}
}
- gmshCollapseSmallEdges (*m);
+ // gmshCollapseSmallEdges (*m);
double t2 = Cpu();
CTX.mesh_timer[1] = t2 - t1;
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 6f2ddec5bc7a3a4605f168f2a32a7e93f9569779..839253efb3dba1fbdd246ffe4ec30b6b967aa335 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -439,7 +439,8 @@ void localHarmonicMapping(GModel *gm,
gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
gmshAssembler myAssembler(lsys);
- gmshLaplaceTerm Laplace (gm,1.0,0);
+ gmshFunction f(1.0);
+ gmshLaplaceTerm Laplace (gm,&f,0);
myAssembler.fixVertex ( n1 , 0 , 0 , -1.0);
myAssembler.fixVertex ( n2 , 0 , 0 , -1.0);
@@ -453,7 +454,7 @@ void localHarmonicMapping(GModel *gm,
gmshLinearSystemGmm *lsys1 = new gmshLinearSystemGmm;
gmshAssembler myAssembler1(lsys1);
- gmshLaplaceTerm Laplace1 (gm,1.0,1);
+ gmshLaplaceTerm Laplace1 (gm,&f,1);
myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
myAssembler1.fixVertex ( n3 , 0 , 1 , -1.0);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 3c448a3f3f466afe6ba70bd91ab27ab98b60674e..f6eb2c1a32b950550cfa56435e9ee157341d33a5 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -29,6 +29,59 @@
extern Context_T CTX;
+void computeEdgeLoops(const GFace *gf,
+ std::vector<MVertex*> &all_mvertices,
+ std::vector<int> &indices)
+{
+ std::list<GEdge*> edges = gf->edges();
+ std::list<int> ori = gf->orientations();
+ std::list<GEdge*>::iterator it = edges.begin();
+ std::list<int>::iterator ito = ori.begin();
+
+ indices.push_back(0);
+ GVertex *start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ GVertex *v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ all_mvertices.push_back(start->mesh_vertices[0]);
+ if (*ito == 1)
+ for (unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
+ all_mvertices.push_back((*it)->mesh_vertices[i]);
+ else
+ for (int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
+ all_mvertices.push_back((*it)->mesh_vertices[i]);
+
+ GVertex *v_start = start;
+ while(1){
+ ++it;
+ ++ito;
+ if(v_end == start){
+ indices.push_back(all_mvertices.size());
+ if(it == edges.end ()) break;
+ start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ v_start = start;
+ }
+ else{
+ if(it == edges.end ()){
+ Msg::Error("Something wrong in edge loop computation");
+ return;
+ }
+ v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ if(v_start != v_end){
+ Msg::Error("Something wrong in edge loop computation");
+ return;
+ }
+ v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
+ }
+ all_mvertices.push_back(v_start->mesh_vertices[0]);
+ if(*ito == 1)
+ for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
+ all_mvertices.push_back((*it)->mesh_vertices[i]);
+ else
+ for (int i = (*it)->mesh_vertices.size()-1; i >= 0; i--)
+ all_mvertices.push_back((*it)->mesh_vertices[i]);
+ }
+}
+
void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
double c[3];
@@ -149,59 +202,6 @@ static bool AlgoDelaunay2D(GFace *gf)
return false;
}
-void computeEdgeLoops(const GFace *gf,
- std::vector<MVertex*> &all_mvertices,
- std::vector<int> &indices)
-{
- std::list<GEdge*> edges = gf->edges();
- std::list<int> ori = gf->orientations();
- std::list<GEdge*>::iterator it = edges.begin();
- std::list<int>::iterator ito = ori.begin();
-
- indices.push_back(0);
- GVertex *start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- GVertex *v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- all_mvertices.push_back(start->mesh_vertices[0]);
- if (*ito == 1)
- for (unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- all_mvertices.push_back((*it)->mesh_vertices[i]);
- else
- for (int i = (*it)->mesh_vertices.size() - 1; i >= 0; i--)
- all_mvertices.push_back((*it)->mesh_vertices[i]);
-
- GVertex *v_start = start;
- while(1){
- ++it;
- ++ito;
- if(v_end == start){
- indices.push_back(all_mvertices.size());
- if(it == edges.end ()) break;
- start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- v_start = start;
- }
- else{
- if(it == edges.end ()){
- Msg::Error("Something wrong in edge loop computation");
- return;
- }
- v_start = ((*ito) == 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- if(v_start != v_end){
- Msg::Error("Something wrong in edge loop computation");
- return;
- }
- v_end = ((*ito) != 1) ? (*it)->getBeginVertex() : (*it)->getEndVertex();
- }
- all_mvertices.push_back(v_start->mesh_vertices[0]);
- if(*ito == 1)
- for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
- all_mvertices.push_back((*it)->mesh_vertices[i]);
- else
- for (int i = (*it)->mesh_vertices.size()-1; i >= 0; i--)
- all_mvertices.push_back((*it)->mesh_vertices[i]);
- }
-}
-
void computeElementShapes(GFace *gf, double &worst, double &avg, double &best,
int &nT, int &greaterThan)
{
@@ -230,6 +230,51 @@ static bool recover_medge(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
g = m->get_geom(ge->tag(),1);
}
+ for (unsigned int i = 0; i < ge->lines.size(); i++){
+ MVertex *vstart = ge->lines[i]->getVertex(0);
+ MVertex *vend = ge->lines[i]->getVertex(1);
+ if (pass_ == 1) e2r->insert(EdgeToRecover(vstart->getIndex(), vend->getIndex(), ge));
+ else{
+ BDS_Edge *e = m->recover_edge(vstart->getIndex(), vend->getIndex(), e2r, not_recovered);
+ if (e) e->g = g;
+ else {
+ // Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
+ // vstart->x(), vstart->y(), vend->x(), vend->y(), i,
+ // ge->mesh_vertices.size());
+ // return false;
+ }
+ }
+ }
+
+ if (pass_ == 2 && ge->getBeginVertex()){
+ MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
+ MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
+ BDS_Point *pstart = m->find_point(vstart->getIndex());
+ BDS_Point *pend = m->find_point(vend->getIndex());
+ if(!pstart->g){
+ m->add_geom (vstart->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vstart->getIndex(), 0);
+ pstart->g = g0;
+ }
+ if(!pend->g){
+ m->add_geom (vend->getIndex(), 0);
+ BDS_GeomEntity *g0 = m->get_geom(vend->getIndex(), 0);
+ pend->g = g0;
+ }
+ }
+
+ return true;
+}
+
+static bool recover_medge_old(BDS_Mesh *m, GEdge *ge, std::set<EdgeToRecover> *e2r,
+ std::set<EdgeToRecover> *not_recovered, int pass_)
+{
+ BDS_GeomEntity *g = 0;
+ if (pass_ == 2){
+ m->add_geom (ge->tag(), 1);
+ g = m->get_geom(ge->tag(),1);
+ }
+
if(ge->mesh_vertices.size() == 0){
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
@@ -342,12 +387,6 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
all_vertices.insert((*it)->lines[i]->getVertex(0));
all_vertices.insert((*it)->lines[i]->getVertex(1));
}
- // all_vertices.insert((*it)->mesh_vertices.begin(),
- // (*it)->mesh_vertices.end());
- // all_vertices.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
- // (*it)->getBeginVertex()->mesh_vertices.end());
- // all_vertices.insert((*it)->getEndVertex()->mesh_vertices.begin(),
- // (*it)->getEndVertex()->mesh_vertices.end());
}
++it;
}
@@ -1318,9 +1357,6 @@ void meshGFace::operator() (GFace *gf)
// compute loops on the fly (indices indicate start and end points
// of a loop; loops are not yet oriented)
Msg::Debug("Computing edge loops");
- // std::vector<MVertex*> points;
- // std::vector<int> indices;
- // computeEdgeLoops(gf, points, indices);
Msg::Debug("Generating the mesh");
if(noseam(gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index afa2da7ed02bff5d0c9f00fa85b5d99398e932fe..b3315b208ab3f9193dcfdbe70aa9f3fd3c0d82c3 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -572,6 +572,9 @@ static void insertAPoint(GFace *gf,
GPoint p = gf->point(center[0], center[1]);
// printf("the new point is %g %g\n",p.x(),p.y());
MVertex *v = new MFaceVertex(p.x(), p.y(), p.z(), gf, center[0], center[1]);
+
+ // printf("%g %g -> %g %g %g\n",center[0], center[1],v->x(), v->y(), v->z());
+
v->setNum(Us.size());
double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] +
uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
@@ -591,8 +594,9 @@ static void insertAPoint(GFace *gf,
AllTris.insert(worst);
delete v;
}
- else
+ else {
gf->mesh_vertices.push_back(v);
+ }
}
else {
Msg::Debug("Point %g %g is outside (%g %g , %g %g , %g %g) (metric %g %g %g)",
@@ -652,11 +656,11 @@ void gmshBowyerWatson(GFace *gf)
circumCenterMetric(worst->tri(), metric, Us, Vs, center, r2);
insertAPoint(gf,AllTris.begin(),center,metric,Us,Vs,vSizes,vSizesBGM,AllTris);
}
-// if(ITER % 1000== 0){
-// char name[245];
-// sprintf(name,"del2d%d-ITER%d.pos",gf->tag(),ITER);
-// _printTris (name, AllTris, Us,Vs,false);
-// }
+ // if(ITER % 10== 0){
+ // char name[245];
+ // sprintf(name,"del2d%d-ITER%d.pos",gf->tag(),ITER);
+ // _printTris (name, AllTris, Us,Vs,false);
+ // }
}
transferDataStructure(gf, AllTris);
}
@@ -794,7 +798,7 @@ void gmshBowyerWatsonFrontal(GFace *gf){
char name[245];
sprintf(name,"frontal%d.pos",gf->tag());
- //_printTris (name, AllTris, Us,Vs);
+ _printTris (name, AllTris, Us,Vs);
transferDataStructure(gf, AllTris);
}
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index d09753228a2e424940f14273ffb6aa2af8430864..1410d342654a555d378c8e8c901e0ecb4ade3318 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -26,6 +26,8 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, Double_Matrix & m) const
const double w = GP[i].pt[2];
const double weight = GP[i].weight;
const double detJ = e->getJacobian (u,v,w,jac);
+ SPoint3 p; e->pnt(u,v,w,p);
+ const double _diff = (*_diffusivity)(p.x(),p.y(),p.z());
inv3x3 ( jac, invjac) ;
e->getGradShapeFunctions (u,v,w,grads);
for (int j=0;j<nbNodes;j++){
@@ -37,7 +39,7 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, Double_Matrix & m) const
for (int k=0;k<=j;k++){
m(j,k) += (Grads[j][0]*Grads[k][0]+
Grads[j][1]*Grads[k][1]+
- Grads[j][2]*Grads[k][2]) * weight * detJ * _diffusivity*0.5;
+ Grads[j][2]*Grads[k][2]) * weight * detJ * _diff*0.5;
}
}
}
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index 34a5fee87c377601a0c4094e1fa0b46e63be13ed..ae6444cf2093d8c61500c933ea09c0d0d6f3d537 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -7,14 +7,15 @@
#define _GMSH_LAPLACE_H_
#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "GmshMatrix.h"
class gmshLaplaceTerm : public gmshNodalFemTerm {
- const double _diffusivity;
const int _iField ;
+ const gmshFunction * _diffusivity;
protected:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
@@ -25,10 +26,9 @@ class gmshLaplaceTerm : public gmshNodalFemTerm {
*iFieldR = _iField;
}
public:
- gmshLaplaceTerm(GModel *gm, double diffusivity = 1.0, int iField = 0) :
- gmshNodalFemTerm(gm), _diffusivity(diffusivity), _iField(iField){}
+ gmshLaplaceTerm(GModel *gm, gmshFunction *diffusivity, int iField = 0) :
+ gmshNodalFemTerm(gm),_diffusivity(diffusivity), _iField(iField){}
virtual void elementMatrix(MElement *e, Double_Matrix &m) const;
- double getDiffusivity () const { return _diffusivity; }
};
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index cb06605f853e2438cdfb417f783c72f8fd16bb89..6f27a41781e435c693d29754fa24cea763469a66 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/Numeric/gmshLinearSystemGmm.h
@@ -18,8 +18,11 @@
class gmshLinearSystemGmm : public gmshLinearSystem {
gmm::row_matrix<gmm::wsvector<double> > *_a;
std::vector<double> *_b, *_x;
+ double _prec;
+ int _noisy;
+ int _gmres;
public :
- gmshLinearSystemGmm () : _a(0), _b(0), _x(0) {}
+ gmshLinearSystemGmm () : _a(0), _b(0), _x(0), _prec(1.e-8), _noisy(0), _gmres(0) {}
virtual bool isAllocated () const {return _a != 0;}
virtual void allocate (int _nbRows)
{
@@ -64,14 +67,17 @@ public :
{
for (unsigned int i = 0; i < _b->size(); i++) (*_b)[i] = 0;
}
+ void setPrec(double p){_prec=p;}
+ void setNoisy(int n){_noisy=n;}
+ void setGmres(int n){_gmres=n;}
virtual int systemSolve ()
{
// gmm::ilutp_precond< gmm::row_matrix< gmm::rsvector<double> > > P(*_a, 10,0.);
gmm::ildltt_precond< gmm::row_matrix< gmm::wsvector<double> > > P(*_a, 2, 1.e-10);
- gmm::iteration iter(1E-8); // defines an iteration object, with a max residu of 1E-8
- //iter.set_noisy(2);
- //gmm::gmres(*_a, *_x, *_b, P, 100, iter); // execute the GMRES algorithm
- gmm::cg(*_a, *_x, *_b, P, iter); // execute the CG algorithm
+ gmm::iteration iter(_prec); // defines an iteration object, with a max residu of 1E-8
+ iter.set_noisy(_noisy);
+ if (_gmres)gmm::gmres(*_a, *_x, *_b, P, 100, iter); // execute the GMRES algorithm
+ else gmm::cg(*_a, *_x, *_b, P, iter); // execute the CG algorithm
return 1;
}
};
diff --git a/benchmarks/3d/Sphere.geo b/benchmarks/3d/Sphere.geo
index a05b93386d81b79ac3d52017754f32c4119a0a68..b48c3ef31aa6b4df8b2d2a607a3592326e1aaf51 100644
--- a/benchmarks/3d/Sphere.geo
+++ b/benchmarks/3d/Sphere.geo
@@ -1,15 +1,25 @@
-lc = .2;
-Point(1) = {0.0,0.0,0.0,lc};
-Point(2) = {1,0.0,0.0,lc};
-Point(3) = {0,1,0.0,lc};
+// Gmsh project created on Fri Feb 1 12:17:44 2008
+
+Delete All;
+
+RSphere = 1.;
+lcSphere = .25;
+
+RDom = 10;
+lcDom = 1.;
+
+
+Point(1) = {0,0,0,lcSphere};
+Point(2) = {RSphere,0,0,lcSphere};
+Point(3) = {0,RSphere,0,lcSphere};
Circle(1) = {2,1,3};
-Point(4) = {-1,0,0.0,lc};
-Point(5) = {0,-1,0.0,lc};
+Point(4) = {-RSphere,0,0.0,lcSphere};
+Point(5) = {0,-RSphere,0.0,lcSphere};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};
-Point(6) = {0,0,-1,lc};
-Point(7) = {0,0,1,lc};
+Point(6) = {0,0,-RSphere,lcSphere};
+Point(7) = {0,0,RSphere,lcSphere};
Circle(5) = {3,1,6};
Circle(6) = {6,1,5};
Circle(7) = {5,1,7};
@@ -18,7 +28,7 @@ Circle(9) = {2,1,7};
Circle(10) = {7,1,4};
Circle(11) = {4,1,6};
Circle(12) = {6,1,2};
-Line Loop(13) = {-2,-8,10};
+Line Loop(13) = {2,8,-10};
Ruled Surface(14) = {13};
Line Loop(15) = {10,3,7};
Ruled Surface(16) = {15};
@@ -26,17 +36,61 @@ Line Loop(17) = {-8,-9,1};
Ruled Surface(18) = {17};
Line Loop(19) = {-11,-2,5};
Ruled Surface(20) = {19};
-Line Loop(21) = {5,12,1};
+Line Loop(21) = {-5,-12,-1};
Ruled Surface(22) = {21};
Line Loop(23) = {-3,11,6};
Ruled Surface(24) = {23};
Line Loop(25) = {-7,4,9};
Ruled Surface(26) = {25};
-Line Loop(27) = {4,-12,6};
+Line Loop(27) = {-4,12,-6};
Ruled Surface(28) = {27};
+Surface Loop(29) = {28,26,16,14,20,24,22,18};
+
+Point(101) = {0,0,0,lcDom};
+Point(102) = {RDom,0,0,lcDom};
+Point(103) = {0,RDom,0,lcDom};
+Circle(101) = {102,101,103};
+Point(104) = {-RDom,0,0.0,lcDom};
+Point(105) = {0,-RDom,0.0,lcDom};
+Circle(102) = {103,101,104};
+Circle(103) = {104,101,105};
+Circle(104) = {105,101,102};
+Point(106) = {0,0,-RDom,lcDom};
+Point(107) = {0,0,RDom,lcDom};
+Circle(105) = {103,101,106};
+Circle(106) = {106,101,105};
+Circle(107) = {105,101,107};
+Circle(108) = {107,101,103};
+Circle(109) = {102,101,107};
+Circle(110) = {107,101,104};
+Circle(111) = {104,101,106};
+Circle(112) = {106,101,102};
+Line Loop(113) = {102,108,-110};
+Ruled Surface(114) = {113};
+Line Loop(115) = {110,103,107};
+Ruled Surface(116) = {115};
+Line Loop(117) = {-108,-109,101};
+Ruled Surface(118) = {117};
+Line Loop(119) = {-111,-102,105};
+Ruled Surface(120) = {119};
+Line Loop(121) = {-105,-112,-101};
+Ruled Surface(122) = {121};
+Line Loop(123) = {-103,111,106};
+Ruled Surface(124) = {123};
+Line Loop(125) = {-107,104,109};
+Ruled Surface(126) = {125};
+Line Loop(127) = {-104,112,-106};
+Ruled Surface(128) = {127};
+Surface Loop(129) = {128,126,116,114,120,124,122,118};
+
+Volume(200) = {129,29};
+
+
+// Wall ID = 2100
+Physical Surface(2100) = {28,26,16,14,20,24,22,18};
+// Farfield CBC ID = 2222
+Physical Surface(2222) = {128,126,116,114,120,124,122,118};
+// Interior ID = 1000
+Physical Volume(1000) = {200};
-Surface Loop(29) = {-28,26,16,-14,20,24,-22,18};
-Volume(30) = {29};
-Physical Volume(1) = {0,30};
-Physical Surface(2) = {14:28:2};