diff --git a/Geo/GModel.h b/Geo/GModel.h
index 6bf05e1046fcfe8002c171df132157df964c38b3..7039a68b7a069f75fc62d885a1f537c1d7179098 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -364,7 +364,7 @@ class GModel
// mesh the model
int mesh(int dimension);
- // adapt the mesh anisotropically using a metric that is computed from a function f(x,y,z).
+ // adapt the mesh anisotropically using a metric that is computed from a scalar function f(x,y,z).
// One can either
// For all algorithms
// parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index b110564b82a08c21d24fe1ea82709db501e03685..7a2ca832ba8b90ca87ab73ec81deb944c64a4103 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -302,6 +302,12 @@ gLevelset::gLevelset(const gLevelset &lv)
tag_ = lv.tag_;
}
+gLevelsetPlane::gLevelsetPlane(const std::vector<double> &pt, const std::vector<double> &norm, int tag) : gLevelsetPrimitive(tag) {
+ a = norm[0];
+ b = norm[1];
+ c = norm[2];
+ d = -a * pt[0] - b * pt[1] - c * pt[2];
+}
gLevelsetPlane::gLevelsetPlane(const double * pt, const double *norm, int tag) : gLevelsetPrimitive(tag) {
a = norm[0];
b = norm[1];
@@ -322,7 +328,6 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv
}
//level set defined by points (RBF interpolation)
-
fullMatrix<double> gLevelsetPoints::generateRbfMat(int p, int index,
const fullMatrix<double> &nodes1,
const fullMatrix<double> &nodes2) const {
@@ -675,13 +680,14 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
return 1.;
}
-gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int tag) : gLevelsetPrimitive(tag) {
+gLevelsetDistGeom::gLevelsetDistGeom(std::string box, std::string geom, int tag) : gLevelsetPrimitive(tag), _box(NULL) {
- GModel *gmg = new GModel();
- gmg->load(geomBox+std::string(".msh"));
- GModel *gm = new GModel();
- gm->load(name);
+ modelBox = new GModel();
+ modelBox->load(box+std::string(".msh"));
+ modelGeom = new GModel();
+ modelGeom->load(geom);
+ //EMI FIXME THIS
double lx = 0.5, ly = 0.5, lz = 0.5;
int levels = 3;
int refineCurvedSurfaces = 0;
@@ -693,14 +699,14 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int
//FILLING POINTS FROM GEOMBOX
std::vector<SPoint3> points;
std::vector<SPoint3> refinePoints;
- for(GModel::viter vit = gmg->firstVertex(); vit != gmg->lastVertex(); vit++){
+ for(GModel::viter vit = modelBox->firstVertex(); vit != modelBox->lastVertex(); vit++){
for(unsigned int k = 0; k < (*vit)->getNumMeshVertices(); k++){
MVertex *v = (*vit)->getMeshVertex(k);
SPoint3 p(v->x(), v->y(), v->z());
points.push_back(p);
}
}
- for (GModel::eiter eit = gmg->firstEdge(); eit != gmg->lastEdge(); eit++){
+ for (GModel::eiter eit = modelBox->firstEdge(); eit != modelBox->lastEdge(); eit++){
for(unsigned int k = 0; k < (*eit)->getNumMeshVertices(); k++){
MVertex *ve = (*eit)->getMeshVertex(k);
SPoint3 pe(ve->x(), ve->y(), ve->z());
@@ -709,7 +715,7 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int
}
//FILLING POINTS FROM STL
- for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
+ for (GModel::fiter fit = modelGeom->firstFace(); fit != modelGeom->lastFace(); fit++){
for(unsigned int k = 0; k < (*fit)->getNumMeshVertices(); k++){
MVertex *vf = (*fit)->getMeshVertex(k);
SPoint3 pf(vf->x(), vf->y(), vf->z());
@@ -718,18 +724,18 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int
for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
MElement *e = (*fit)->getMeshElement(k);
if (e->getType() == TYPE_TRI){
- MVertex *v1 = e->getVertex(0);
- MVertex *v2 = e->getVertex(1);
- MVertex *v3 = e->getVertex(2);
- SPoint3 cg( (v1->x()+v2->x()+v3->x())/3.,
- (v1->y()+v2->y()+v3->y())/3.,
- (v1->z()+v2->z()+v3->z())/3.);
- refinePoints.push_back(cg);
+ MVertex *v1 = e->getVertex(0);
+ MVertex *v2 = e->getVertex(1);
+ MVertex *v3 = e->getVertex(2);
+ SPoint3 cg( (v1->x()+v2->x()+v3->x())/3.,
+ (v1->y()+v2->y()+v3->y())/3.,
+ (v1->z()+v2->z()+v3->z())/3.);
+ refinePoints.push_back(cg);
}
}
}
//FOR CAD
- //for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
+ //for (GModel::fiter fit = modelGeom->firstFace(); fit != modelGeom->lastFace(); fit++)
// (*fit)->fillPointCloud(sampling, &points);
if (points.size() == 0) {Msg::Fatal("No points on surfaces \n"); };
@@ -751,48 +757,46 @@ gLevelsetDistGeom::gLevelsetDistGeom(std::string geomBox, std::string name, int
// bb.max().x(), bb.max().y(), bb.max().z());
// Msg::Info(" Nx=%d Ny=%d Nz=%d", NX, NY, NZ);
- _box = new cartesianBox<double>(bb.min().x(), bb.min().y(), bb.min().z(),
- SVector3(range.x(), 0, 0),
- SVector3(0, range.y(), 0),
- SVector3(0, 0, range.z()),
- NX, NY, NZ, levels);
- for (int i = 0; i < NX; i++)
- for (int j = 0; j < NY; j++)
- for (int k = 0; k < NZ; k++)
- _box->insertActiveCell(_box->getCellIndex(i, j, k));
-
- cartesianBox<double> *parent = _box, *child;
- while((child = parent->getChildBox())){
- //Msg::Info(" level %d ", child->getLevel());
- for(unsigned int i = 0; i < refinePoints.size(); i++)
- insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
- rtube / pow(2., (levels - child->getLevel())), *child);
- parent = child;
- }
-
- //Msg::Info("Removing cells to match mesh topology constraints");
- removeBadChildCells(_box);
- removeParentCellsWithChildren(_box);
-
- //Msg::Info("Initializing nodal values in the cartesian grid");
- _box->createNodalValues();
-
- //Msg::Info("Computing levelset on the cartesian grid");
- computeLevelset(gm, *_box);
-
- //Msg::Info("Renumbering mesh vertices across levels");
- _box->renumberNodes();
+ // _box = new cartesianBox<double>(bb.min().x(), bb.min().y(), bb.min().z(),
+ // SVector3(range.x(), 0, 0),
+ // SVector3(0, range.y(), 0),
+ // SVector3(0, 0, range.z()),
+ // NX, NY, NZ, levels);
+ // for (int i = 0; i < NX; i++)
+ // for (int j = 0; j < NY; j++)
+ // for (int k = 0; k < NZ; k++)
+ // _box->insertActiveCell(_box->getCellIndex(i, j, k));
+
+ // cartesianBox<double> *parent = _box, *child;
+ // while((child = parent->getChildBox())){
+ // //Msg::Info(" level %d ", child->getLevel());
+ // for(unsigned int i = 0; i < refinePoints.size(); i++)
+ // insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
+ // rtube / pow(2., (levels - child->getLevel())), *child);
+ // parent = child;
+ // }
+
+ // //Msg::Info("Removing cells to match mesh topology constraints");
+ // removeBadChildCells(_box);
+ // removeParentCellsWithChildren(_box);
+
+ // //Msg::Info("Initializing nodal values in the cartesian grid");
+ // _box->createNodalValues();
+
+ // //Msg::Info("Computing levelset on the cartesian grid");
+ // computeLevelset(gm, *_box);
+
+ // //Msg::Info("Renumbering mesh vertices across levels");
+ // _box->renumberNodes();
- _box->writeMSH("yeah.msh", false);
+ // _box->writeMSH("yeah.msh", false);
- delete gm;
- delete gmg;
-
+
}
double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const {
- double dist = _box->getValueContainingPoint(x,y,z);
+ double dist = 0.0 ; //_box->getValueContainingPoint(x,y,z);
return dist;
}
@@ -920,6 +924,19 @@ gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *p
gLevelsetBox::gLevelsetBox(const gLevelsetBox &lv) : gLevelsetImproved(lv){}
+gLevelsetCylinder::gLevelsetCylinder(const std::vector<double> &pt, const std::vector<double> &dir, const double &R, const double &H, int tag) : gLevelsetImproved() {
+ double pt1[3]={pt[0], pt[1], pt[2]};
+ double dir1[3] = {dir[0], dir[1], dir[2]};
+ double dir2[3] = {-dir1[0], -dir1[1], -dir1[2]};
+ double n[3]; norm(dir1, n);
+ double pt2[3] = {pt1[0] + H * n[0], pt1[1] + H * n[1], pt1[2] + H * n[2]};
+ std::vector<gLevelset *> p;
+ p.push_back(new gLevelsetGenCylinder(pt1, dir1, R, tag++));
+ p.push_back(new gLevelsetPlane(pt1, dir2, tag++));
+ p.push_back(new gLevelsetPlane(pt2, dir1, tag));
+ Ls = new gLevelsetIntersection(p);
+}
+
gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const double &R, const double &H, int tag) : gLevelsetImproved() {
double dir2[3] = {-dir[0], -dir[1], -dir[2]};
double n[3]; norm(dir, n);
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 7ce2a06cfda44e94154069d8d57e861e5527d0f9..838115cb60c0b73dbafc43656e711ec7d1d0285a 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -46,6 +46,7 @@
#define UNION 13
#define INTER 14
#define CRACK 15
+#define DISK 16
class gLevelset : public simpleFunction<double>
{
@@ -143,6 +144,7 @@ public:
// define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
gLevelsetPlane (const double _a, const double _b, const double _c, const double _d, int tag=1) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
// define the plane passing through the point pt and with outward normal norm
+ gLevelsetPlane (const std::vector<double> &pt, const std::vector<double> &norm, int tag=1);
gLevelsetPlane (const double *pt, const double *norm, int tag=1);
// define the plane passing through the 3 points pt1,pt2,pt3 and with outward normal (pt1,pt2)x(pt1,pt3)
gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int tag=1);
@@ -155,6 +157,8 @@ public:
int type() const {return PLANE;}
};
+
+
class gLevelsetPoints : public gLevelsetPrimitive
{
protected:
@@ -276,9 +280,14 @@ public:
class gLevelsetDistGeom: public gLevelsetPrimitive
{
cartesianBox<double> *_box;
+ GModel *modelBox, *modelGeom;
public:
gLevelsetDistGeom(std::string g, std::string f, int tag=1);
- ~gLevelsetDistGeom(){if (_box) delete _box; }
+ ~gLevelsetDistGeom(){
+ delete modelBox;
+ delete modelGeom;
+ if (_box) delete _box;
+ }
double operator () (const double x, const double y, const double z) const;
int type() const { return UNKNOWN; }
};
@@ -404,12 +413,13 @@ public:
class gLevelsetCrack : public gLevelsetTools
{
public:
- gLevelsetCrack (std::vector<gLevelset *> p) {
+ gLevelsetCrack (std::vector<gLevelset *> p, bool delC=false) {
if (p.size() != 2)
printf("Error : gLevelsetCrack needs 2 levelsets\n");
children.push_back(p[0]);
children.push_back(new gLevelsetReverse(p[0]));
if(p[1]) children.push_back(p[1]);
+ _delChildren = delC;
}
double choose (double d1, double d2) const {
return (d1 > d2) ? d1 : d2; // greater of d1 and d2
@@ -417,6 +427,7 @@ public:
int type2() const {return CRACK;}
};
+
// --------------------------------------------------------------------------------------------------------------
// IMPROVED LEVELSET
class gLevelsetImproved : public gLevelset
@@ -476,6 +487,7 @@ public:
// tags of the faces are : exterior face : tag+0
// plane face including pt : tag+1
// plane face opposite to pt : tag+2
+ gLevelsetCylinder (const std::vector<double> &pt, const std::vector<double> &dir, const double &R, const double &H, int tag=1);
gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int tag=1);
// create a cylinder : pt is the point in the middle of the cylinder base,
// dir is the direction of the cylinder axis,
@@ -529,5 +541,32 @@ public:
int type() const {return CONROD;}
};
+/* class gLevelsetDisk : public gLevelsetTools */
+/* { */
+/* public: */
+/* // define the disk of given radius centered at a point pt and with outward normal norm */
+/* gLevelsetDisk (std::vector<double> pt, std::vector<double> dir, const double R, bool delC=false) { */
+/* double *ptP, *normP; */
+/* ptP[0] = pt[0]; ptP[1] = pt[1]; ptP[2] = pt[2]; */
+/* normP[0] = dir[0];normP[1] = dir[1]; normP[2] = dir[2]; */
+/* gLevelsetPlane *plane = new gLevelsetPlane(ptP, normP); */
+/* children.push_back(plane); */
+/* children.push_back(new gLevelsetReverse(plane)); */
+/* double H = 4.*R; */
+/* double *ptC, *normC; */
+/* double val = sqrt(dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]); */
+/* normC[0] = dir[0]/val; normC[1] = dir[1]/val; normC[2] = dir[2]/val; */
+/* ptC[0] = pt[0]-2.*R*normC[0]; */
+/* ptC[1] = pt[1]-2.*R*normC[1]; */
+/* ptC[2] = pt[2]-2.*R*normC[2]; */
+/* gLevelsetCylinder *cyl = new gLevelsetCylinder(ptC, normC, R, H); */
+/* children.push_back(cyl); */
+/* _delChildren = delC; */
+/* } */
+/* double choose (double d1, double d2) const { */
+/* return (d1 > d2) ? d1 : d2; // greater of d1 and d2 */
+/* } */
+/* int type2() const {return DISK;} */
+/* }; */
#endif
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 876f3a7e301be420d82d04f6064ea5cd757978f2..ed2b0353f818cbb3e95c4c4b536014016e97e521 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -723,7 +723,7 @@ double backgroundMesh::operator() (double u, double v, double w) const
double uv2[3];
MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("cannot find %g %g %g", u, v, w);
+ Msg::Error("BGM octree: cannot find %g %g %g", u, v, w);
return -1000.0;//0.4;
}
e->xyz2uvw(uv, uv2);
@@ -739,7 +739,7 @@ double backgroundMesh::getAngle(double u, double v, double w) const
double uv2[3];
MElement *e = _octree->find(u, v, w, 2, true);
if (!e) {
- Msg::Error("cannot find %g %g %g", u, v, w);
+ Msg::Error("BGM octree : cannot find %g %g %g", u, v, w);
return 0.0;
}
e->xyz2uvw(uv, uv2);
diff --git a/contrib/DiscreteIntegration/recurCut.h b/contrib/DiscreteIntegration/recurCut.h
index 0bb3799810bff3ba114fdda8c3608d12517c3602..8633ed879ce385e3f8503530b64c23e31f1af637 100644
--- a/contrib/DiscreteIntegration/recurCut.h
+++ b/contrib/DiscreteIntegration/recurCut.h
@@ -3,8 +3,8 @@
#ifndef _RECUR_H_
#define _RECUR_H_
-#include "gmshLevelset.h"
#include "Integration3D.h"
+class gLevelset;
class RecurElement
{