From e870360fc4a67ddded57bec04718ede34bc0581d Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Mon, 28 Nov 2011 14:11:29 +0000
Subject: [PATCH] fixed bug cut and split for elems
---
Geo/GFace.cpp | 4 +-
Geo/GFace.h | 2 +-
Geo/MElement.cpp | 4 +-
Geo/MElementCut.cpp | 24 ++--
Geo/gmshLevelset.cpp | 293 +++++++++++++++++++++++++++++++++++++------
Numeric/cartesian.h | 21 +++-
Plugin/Distance.cpp | 20 +--
7 files changed, 300 insertions(+), 68 deletions(-)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index f877f61eb2..9eda7e6c0a 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1017,10 +1017,10 @@ bool GFace::buildSTLTriangulation(bool force)
stl_vertices.clear();
stl_triangles.clear();
}
- else
+ else{
return true;
+ }
}
-
// Build a simple triangulation for surfaces which we know are not
// trimmed
if(geomType() == ParametricSurface || geomType() == ProjectionFace){
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 509eb6573c..35700f9068 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -254,7 +254,7 @@ class GFace : public GEntity
void setCompound(GFaceCompound *gfc) { compound = gfc; }
GFaceCompound *getCompound() const { return compound; }
- // add points (and optionalluy normals) in vectors so that two
+ // add points (and optionally normals) in vectors so that two
// points are at most maxDist apart
bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
std::vector<SVector3> *normals=0);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index da8627e9e8..6f45c85c18 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1093,7 +1093,7 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
if(getNumChildren() == 0) {
for(int i = 0; i < getNumVertices(); i++) {
MVertex *v = getVertex(i);
- int numV = v->getIndex();
+ int numV = v->getNum(); //Index();
if(vertexMap.count(numV))
vmv.push_back(vertexMap[numV]);
else {
@@ -1107,7 +1107,7 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
for(int i = 0; i < getNumChildren(); i++) {
for(int j = 0; j < getChild(i)->getNumVertices(); j++) {
MVertex *v = getChild(i)->getVertex(j);
- int numV = v->getIndex();
+ int numV = v->getNum(); //Index();
if(vertexMap.count(numV))
vmv.push_back(vertexMap[numV]);
else {
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 423e015761..ee41a648ac 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -689,19 +689,19 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
}
//create level set interface (pt in 1D, line in 2D or face in 3D)
- if (splitElem && copy->getDim()==2){
- //printf("*** SPLITELEM \n");
- for(int k = 0; k < copy->getNumEdges(); k++){
- MEdge me = copy->getEdge(k);
+ if (splitElem && e->getDim()==2){
+ for(int k = 0; k < e->getNumEdges(); k++){
+ MEdge me = e->getEdge(k);
MVertex *v0 = me.getVertex(0);
MVertex *v1 = me.getVertex(1);
+ MVertex *v0N = vertexMap[v0->getNum()];
+ MVertex *v1N = vertexMap[v1->getNum()];
double val0 = verticesLs(0, v0->getIndex())-eps;
double val1 = verticesLs(0, v1->getIndex())-eps;
- //printf("splitedge v0= (%g,%g) v1= (%g,%g) val0= %g(%d) val1= %g(%d) nums (%d %d)\n", v0->x(), v0->y(), v1->x(), v1->y(),
- // val0, val1, v0->getIndex(), v1->getIndex(), v0->getNum(), v1->getNum());
+ //printf("splitedge v0= (%g,%g) v1= (%g,%g) val0= %g(%d) val1= %g(%d) nums (%d %d)\n", v0->x(), v0->y(), v1->x(), v1->y(), val0, v0->getIndex(), val1, v1->getIndex(), vertexMap[v0->getNum()]->getIndex(), vertexMap[v1->getNum()]->getIndex());
if ( (val0*val1 > 0.0 && val0 < 0.0 ) ) {
//printf("split edge \n");
- MLine *lin = new MLine(v0, v1);
+ MLine *lin = new MLine(v0N, v1N);
getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
int c = elements[1].count(gLsTag);
int reg = getBorderTag(gLsTag, c, newElemTags[1][0], borderElemTags[0]);
@@ -711,9 +711,9 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
}
}
}
- else if (splitElem && copy->getDim()==3){
- for(int k = 0; k < copy->getNumFaces(); k++){
- MFace mf = copy->getFace(k);
+ else if (splitElem && e->getDim()==3){
+ for(int k = 0; k < e->getNumFaces(); k++){
+ MFace mf = e->getFace(k);
bool sameSign = true;
double val0 = verticesLs(0, mf.getVertex(0)->getIndex());
for (int j= 1; j< mf.getNumVertices(); j++){
@@ -726,7 +726,9 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
int reg = getBorderTag(gLsTag, c, newElemTags[2][0], borderElemTags[1]);
int physTag = (!gePhysicals.size()) ? 0 : getBorderTag(gLsTag, c, newPhysTags[2][0], borderPhysTags[1]);
if(mf.getNumVertices() == 3){
- MTriangle *tri = new MTriangle(mf.getVertex(0), mf.getVertex(1), mf.getVertex(2));
+ MTriangle *tri = new MTriangle(vertexMap[mf.getVertex(0)->getNum()],
+ vertexMap[mf.getVertex(1)->getNum()],
+ vertexMap[mf.getVertex(2)->getNum()]);
elements[2][reg].push_back(tri);
}
else if(mf.getNumVertices() == 4){
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 93cb77cf18..a9e30d497f 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -12,6 +12,98 @@
#include <stack>
#include "fullMatrix.h"
#include "GModel.h"
+#include "MElement.h"
+#include "Numeric.h"
+#include "cartesian.h"
+
+void insertActiveCells(double x, double y, double z, double rmax,
+ cartesianBox<double> &box)
+{
+ int id1 = box.getCellContainingPoint(x - rmax, y - rmax, z - rmax);
+ int id2 = box.getCellContainingPoint(x + rmax, y + rmax, z + rmax);
+ int i1, j1 ,k1;
+ box.getCellIJK(id1, i1, j1, k1);
+ int i2, j2, k2;
+ box.getCellIJK(id2, i2, j2, k2);
+ for (int i = i1; i <= i2; i++)
+ for (int j = j1; j <= j2; j++)
+ for (int k = k1; k <= k2; k++)
+ box.insertActiveCell(box.getCellIndex(i, j, k));
+}
+void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
+{
+ Range<double> t_bounds = ge->parBounds(0);
+ double t_min = t_bounds.low();
+ double t_max = t_bounds.high();
+ double length = ge->length(t_min, t_max, 20);
+ int N = length / sampling;
+ for(int i = 0; i < N; i++) {
+ double t = t_min + (double)i / (double)(N - 1) * (t_max - t_min);
+ GPoint p = ge->point(t);
+ points.push_back(SPoint3(p.x(), p.y(), p.z()));
+ }
+}
+
+int removeBadChildCells(cartesianBox<double> *parent)
+{
+ cartesianBox<double> *child = parent->getChildBox();
+ if(!child) return 0;
+ int I = parent->getNxi();
+ int J = parent->getNeta();
+ int K = parent->getNzeta();
+ for(int i = 0; i < I; i++)
+ for(int j = 0; j < J; j++)
+ for(int k = 0; k < K; k++){
+ // remove children if they do not entirely fill parent, or if
+ // there is no parent on the boundary
+ int idx[8] = {child->getCellIndex(2 * i, 2 * j, 2 * k),
+ child->getCellIndex(2 * i, 2 * j, 2 * k + 1),
+ child->getCellIndex(2 * i, 2 * j + 1, 2 * k),
+ child->getCellIndex(2 * i, 2 * j + 1, 2 * k + 1),
+ child->getCellIndex(2 * i + 1, 2 * j, 2 * k),
+ child->getCellIndex(2 * i + 1, 2 * j, 2 * k + 1),
+ child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k),
+ child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k + 1)};
+ bool exist[8], atLeastOne = false, butNotAll = false;
+ for(int ii = 0; ii < 8; ii++){
+ exist[ii] = child->activeCellExists(idx[ii]);
+ if(exist[ii]) atLeastOne = true;
+ if(!exist[ii]) butNotAll = true;
+ }
+ if(atLeastOne &&
+ (butNotAll ||
+ (i != 0 && !parent->activeCellExists(parent->getCellIndex(i - 1, j, k))) ||
+ (i != I - 1 && !parent->activeCellExists(parent->getCellIndex(i + 1, j, k))) ||
+ (j != 0 && !parent->activeCellExists(parent->getCellIndex(i, j - 1, k))) ||
+ (j != J - 1 && !parent->activeCellExists(parent->getCellIndex(i, j + 1, k))) ||
+ (k != 0 && !parent->activeCellExists(parent->getCellIndex(i, j, k - 1))) ||
+ (k != K - 1 && !parent->activeCellExists(parent->getCellIndex(i, j, k + 1)))))
+ for(int ii = 0; ii < 8; ii++) child->eraseActiveCell(idx[ii]);
+ }
+ return removeBadChildCells(child);
+}
+
+void removeParentCellsWithChildren(cartesianBox<double> *box)
+{
+ if(!box->getChildBox()) return;
+ for(int i = 0; i < box->getNxi(); i++)
+ for(int j = 0; j < box->getNeta(); j++)
+ for(int k = 0; k < box->getNzeta(); k++){
+ if(box->activeCellExists(box->getCellIndex(i, j, k))){
+ cartesianBox<double> *parent = box, *child;
+ int ii = i, jj = j, kk = k;
+ while((child = parent->getChildBox())){
+ ii *= 2; jj *= 2; kk *= 2;
+ if(child->activeCellExists(child->getCellIndex(ii, jj, kk))){
+ box->eraseActiveCell(box->getCellIndex(i, j, k));
+ break;
+ }
+ parent = child;
+ }
+ }
+ }
+ removeParentCellsWithChildren(box->getChildBox());
+}
inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
@@ -495,54 +587,173 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
return 1.;
}
-//EMI TO DO
gLevelsetDistGeom::gLevelsetDistGeom(std::string name, int tag) : gLevelsetPrimitive(tag) {
- _model = new GModel();
- _model->load(name);
- // std::vector<double> distances;
- // distances.clear();
- // for (GModel::fiter fit = _model->firstFace(); fit != _model->lastFace(); fit++){
- // if((*it)->geomType() == GEntity::DiscreteSurface){
- // for(unsigned int k = 0; k < fit->getNumMeshElements(); k++){
- // std::vector<double> iDistances;
- // std::vector<SPoint3> iClosePts;
- // std::vector<double> iDistancesE;
- // MElement *e = fit->getMeshElement(k);
- // MVertex *v1 = e->getVertex(0);
- // MVertex *v2 = e->getVertex(1);
- // SPoint3 p1(v1->x(), v1->y(), v1->z());
- // SPoint3 p2(v2->x(), v2->y(), v2->z());
- // if((e->getNumVertices() == 2 && order==1) || (e->getNumVertices() == 3 && order==2)){
- // signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
- // }
- // else if((e->getNumVertices() == 3 && order == 1) || (e->getNumVertices() == 6 && order==2)){
- // MVertex *v3 = e->getVertex(2);
- // SPoint3 p3 (v3->x(),v3->y(),v3->z());
- // signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
- // }
- // for (unsigned int kk = 0; kk< pts.size(); kk++) {
- // if (std::abs(iDistances[kk]) < distances[kk]){
- // distances[kk] = std::abs(iDistances[kk]);
- // MVertex *v = pt2Vertex[kk];
- // _distance_map[v] = distances[kk];
- // }
- // }
- // }
- // }
- // else{
- // //look in utils_api_demos maincartesian
- // // for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
- // // }
- // }
- // }
+ //_model = new GModel();
+ //_model->load(name);
+
+ GModel *gm = new GModel();
+ gm->load(name);
+
+ double lx = 0.1, ly = 0.1, lz = 0.1;
+ double rmax = 0.05;
+ int levels = 3;
+ int refineCurvedSurfaces = 0;
+
+ double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
+ double rtube = std::max(lx, std::max(ly, lz)) * 2.;
+
+ std::vector<SPoint3> points;
+ std::vector<SPoint3> refinePoints;
+
+ Msg::Info("Filling coarse point cloud on surfaces");
+
+ //FOR CAD
+ for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
+ (*fit)->fillPointCloud(sampling, &points);
+
+ //FOR STL
+ // for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
+ // 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 p1(v1->x(), v1->y(), v1->z());
+ // SPoint3 p2(v2->x(), v2->y(), v2->z());
+ // SPoint3 p3(v3->x(), v3->y(), v3->z());
+ // points.push_back(p1);
+ // points.push_back(p2);
+ // points.push_back(p3);
+ // // std::vector<SPoint3>::iterator it1 = std::find(points.begin(), points.end(), p1);
+ // // if(it1 != points.end()) points.push_back(p1);
+ // // std::vector<SPoint3>::iterator it2 = std::find(points.begin(), points.end(), p2);
+ // // if(it2 != points.end()) points.push_back(p2);
+ // // std::vector<SPoint3>::iterator it3 = std::find(points.begin(), points.end(),p3);
+ // // if(it3 != points.end()) points.push_back(p3);
+ // }
+ // }
+ // }
+ Msg::Info(" %d points in the surface cloud", (int)points.size());
+ if (points.size() == 0) {Msg::Fatal("No points on surfaces \n"); };
+
+
+ // if(levels > 1){
+ // double s = sampling / pow(2., levels - 1);
+ // Msg::Info("Filling refined point cloud on curves and curved surfaces");
+ // for (GModel::eiter eit = gm->firstEdge(); eit != gm->lastEdge(); eit++)
+ // fillPointCloud(*eit, s, refinePoints);
+
+ // // FIXME: refine this by computing e.g. "mean" curvature
+ // if(refineCurvedSurfaces){
+ // for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
+ // if((*fit)->geomType() != GEntity::Plane)
+ // (*fit)->fillPointCloud(2 * s, &refinePoints);
+ // }
+ // Msg::Info(" %d points in the refined cloud", (int)refinePoints.size());
+ // }
+
+ SBoundingBox3d bb;
+ for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
+ for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
+ bb.scale(1.21, 1.21, 1.21);
+ SVector3 range = bb.max() - bb.min();
+ int NX = range.x() / lx;
+ int NY = range.y() / ly;
+ int NZ = range.z() / lz;
+ if(NX < 2) NX = 2;
+ if(NY < 2) NY = 2;
+ if(NZ < 2) NZ = 2;
+
+ Msg::Info(" bounding box min: %g %g %g -- max: %g %g %g",
+ bb.min().x(), bb.min().y(), bb.min().z(),
+ bb.max().x(), bb.max().y(), bb.max().z());
+ Msg::Info(" Nx=%d Ny=%d Nz=%d", NX, NY, NZ);
+
+ cartesianBox<double> box(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);
+
+ Msg::Info("Inserting active cells in the cartesian grid");
+ Msg::Info(" level %d", box.getLevel());
+ for (unsigned int i = 0; i < points.size(); i++)
+ insertActiveCells(points[i].x(), points[i].y(), points[i].z(), rmax, box);
+
+ 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;
+ }
+
+ // remove child cells that do not entirely fill parent cell or for
+ // which there is no parent neighbor; then remove parent cells that
+ // have children
+ Msg::Info("Removing cells to match X-FEM mesh topology constraints");
+ removeBadChildCells(&box);
+ removeParentCellsWithChildren(&box);
+
+ // we generate duplicate nodes at this point so we can easily access
+ // cell values at each level; we will clean up by renumbering after
+ // filtering
+ 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();
+
+ bool decomposeInSimplex = false;
+ box.writeMSH("yeah.msh", decomposeInSimplex);
+ exit(1);
+
}
double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const {
- double dist = 1.0;
+
+ std::vector<SPoint3> nodes;
+ nodes.push_back(SPoint3(x,y,z));
+ double dist = 1.e22;
+
+ for (GModel::fiter fit = _model->firstFace(); fit != _model->lastFace(); fit++){
+ if((*fit)->geomType() == GEntity::DiscreteSurface){
+ for(unsigned int k = 0; k < (*fit)->getNumMeshElements(); k++){
+ std::vector<double> iDistances;
+ std::vector<SPoint3> iClosePts;
+ std::vector<double> iDistancesE;
+ 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 p1(v1->x(), v1->y(), v1->z());
+ SPoint3 p2(v2->x(), v2->y(), v2->z());
+ SPoint3 p3(v3->x(),v3->y(),v3->z());
+ signedDistancesPointsTriangle(iDistances, iClosePts, nodes, p1, p2, p3);
+ }
+ if (std::fabs(iDistances[0]) < std::fabs(dist)) dist = iDistances[0];
+ }
+ }
+ else{
+ printf(" CAD surface \n");
+ //look in utils_api_demos maincartesian
+ // for (GModel::fiter fit = _model->firstFace(); fit != _model->lastFace(); fit++){
+ // for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
+ }
+ }
+ double ana = sqrt((y*y)+(z*z))-0.5;
+ printf("xyz = %g %g %g dist = %g (%g)\n",x, y, z, -dist, ana );
+ //exit(1);
+
return dist;
}
-
#if defined (HAVE_POST)
gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(tag), _viewIndex(index){
if(_viewIndex >= 0 && _viewIndex < (int)PView::list.size()){
diff --git a/Numeric/cartesian.h b/Numeric/cartesian.h
index 0698f7276c..131002ebf0 100644
--- a/Numeric/cartesian.h
+++ b/Numeric/cartesian.h
@@ -45,7 +45,7 @@ class cartesianBox {
// stored a node tag (used for global numbering of the nodes across
// the grid levels)
typename std::map<int, std::pair<scalar, int> > _nodalValues;
- // level of the box (coarset box has highest level; finest box has
+ // level of the box (coarsest box has highest level; finest box has
// level==1)
int _level;
// pointer to a finer (refined by 2) level box, if any
@@ -171,6 +171,25 @@ class cartesianBox {
}
}
}
+ /* double getValueContainingPoint(double x, double y, double z) const */
+ /* { */
+ /* double val; */
+ /* int t = getCellContainingPoint(x, y,z); */
+ /* std::vector<scalar> values; */
+ /* getNodalValuesForCell(t, values); */
+ /* printf("values size =%d \n", values.size()); */
+
+ /* SVector3 DP (x - _X, y - _Y, z - _Z); */
+ /* double xi = dot(DP, _xiAxis); */
+ /* double eta = dot(DP, _etaAxis); */
+ /* double zeta = dot(DP, _zetaAxis); */
+
+ /* double xi_e = xi-values[0]; */
+ /* double eta_e = xi-_dxi/2; */
+ /* double zeta_e = xi-_dxi/2; */
+
+ /* return val; */
+ /* } */
int getCellContainingPoint(double x, double y, double z) const
{
// P = P_0 + xi * _vdx + eta * _vdy + zeta *vdz
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 8e80559364..265c914776 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -200,16 +200,16 @@ PView *GMSH_DistancePlugin::execute(PView *v)
distances.reserve(totNumNodes);
pt2Vertex.reserve(totNumNodes);
- std::map<MVertex*,double> _distanceE_map;
- std::map<MVertex*,int> _isInYarn_map;
- std::vector<int> index;
- std::vector<double> distancesE;
- std::vector<int> isInYarn;
- std::vector<SPoint3> closePts;
- std::vector<double> distances2;
- std::vector<double> distancesE2;
- std::vector<int> isInYarn2;
- std::vector<SPoint3> closePts2;
+ std::map<MVertex*,double> _distanceE_map;
+ std::map<MVertex*,int> _isInYarn_map;
+ std::vector<int> index;
+ std::vector<double> distancesE;
+ std::vector<int> isInYarn;
+ std::vector<SPoint3> closePts;
+ std::vector<double> distances2;
+ std::vector<double> distancesE2;
+ std::vector<int> isInYarn2;
+ std::vector<SPoint3> closePts2;
for (int i = 0; i < totNumNodes; i++) {
distances.push_back(1.e22);
--
GitLab