From c9bd75417e570cc2c635069b3d59e2b6c5eaf584 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Wed, 10 Aug 2011 13:42:02 +0000
Subject: [PATCH] work on levelset computation
---
Geo/MElementCut.cpp | 26 ++++++++-----------
Mesh/meshGRegionMMG3D.cpp | 5 ++--
contrib/DiscreteIntegration/DILevelset.cpp | 11 +++++---
contrib/DiscreteIntegration/DILevelset.h | 4 +--
contrib/DiscreteIntegration/Integration3D.cpp | 20 +++++++-------
5 files changed, 33 insertions(+), 33 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index ec619d73ca..c589eecc39 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -1126,45 +1126,41 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::vector<const gLevelset *> primitives;
ls->getPrimitivesPO(primitives);
+ int primS = primitives.size();
int numVert = gm->indexMeshVertices(true);
- int nbLs = (cutElem) ? primitives.size() : 1;
+ int nbLs = (cutElem) ? ((primS > 1) ? primS + 1 : 1) : 1;
fullMatrix<double> verticesLs(nbLs, numVert + 1);
//Emi test compute all at once for POINTS (type = 11)
std::vector<MVertex *> vert;
- std::map<MVertex*, double> myMap;
for(unsigned int i = 0; i < gmEntities.size(); i++) {
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
vert.push_back(gmEntities[i]->getMeshVertex(j));
}
}
- if(cutElem){
- for(unsigned int k = 0; k < primitives.size(); k++){
- if (primitives[k]->type() == 11){ //points
- ((gLevelsetPoints*)primitives[k])->computeLS(vert, myMap);
- }
+ for(int k = 0; k < primS; k++){
+ if (primitives[k]->type() == 11){ //points
+ ((gLevelsetPoints*)primitives[k])->computeLS(vert);
}
}
- else{
- if (ls->type() == 11)((gLevelsetPoints*)ls)->computeLS(vert, myMap);
- }
//compute and store levelset values
for(unsigned int i = 0; i < gmEntities.size(); i++) {
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
MVertex *vi = gmEntities[i]->getMeshVertex(j);
if(cutElem){
- for(unsigned int k = 0; k < primitives.size(); k++){
- //std::map<MVertex*,double>::iterator it = myMap.find(vi);
- //verticesLs(k, vi->getIndex()) = it->second;
- verticesLs(k, vi->getIndex()) = (*primitives[k])(vi->x(), vi->y(), vi->z());
+ int k = 0;
+ for(; k < primS; k++){
+ verticesLs(k, vi->getIndex()) = (*primitives[k])(vi->x(), vi->y(), vi->z());
}
+ if(primS > 1)
+ verticesLs(k, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
}
else
verticesLs(0, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
}
}
-
+
std::map<int, int> newElemTags[4]; //map<oldElementary,newElementary>[dim]
std::map<int, int> newPhysTags[4]; //map<oldPhysical,newPhysical>[dim]
for(int d = 0; d < 4; d++){
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index 938fc8f38d..56fd75d03a 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -255,8 +255,9 @@ void refineMeshMMG(GRegion *gr){
mmg3d::MMG_mmg3dlib(opt,mmg,sol);
Msg::Debug("-------- MG3D TERMINATED -------------------");
updateSizes (gr,mmg, sol);
- }
- MMG_saveMesh(mmg,"test.mesh");
+ }
+ char test[] = "test.mesh";
+ MMG_saveMesh(mmg, test);
gr->deleteVertexArrays();
for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
diff --git a/contrib/DiscreteIntegration/DILevelset.cpp b/contrib/DiscreteIntegration/DILevelset.cpp
index 18e308d538..1fbe43fb2d 100644
--- a/contrib/DiscreteIntegration/DILevelset.cpp
+++ b/contrib/DiscreteIntegration/DILevelset.cpp
@@ -307,9 +307,14 @@ gLevelsetPoints::gLevelsetPoints(const gLevelsetPoints &lv) : gLevelsetPrimitive
double gLevelsetPoints::operator()(const double &x, const double &y, const double &z) const{
+ if(mapP.empty()) printf("Levelset Points : call computeLS() before calling operator()\n");
+
SPoint3 sp(x,y,z);
std::map<SPoint3,double>::const_iterator it = mapP.find(sp);
- return it->second;
+ if(it != mapP.end())
+ return it->second;
+ printf("Levelset Points : Point not found\n");
+ return 0;
// fullMatrix<double> xyz_eval(1, 3), surf_eval(1,1);
// xyz_eval(0,0) = x;
@@ -320,9 +325,8 @@ double gLevelsetPoints::operator()(const double &x, const double &y, const doubl
}
-void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert, std::map<MVertex*, double> &myMap){
+void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert){
-
fullMatrix<double> xyz_eval(vert.size(), 3), surf_eval(vert.size(), 1);
for (int i = 0; i< vert.size(); i++){
xyz_eval(i,0) = vert[i]->x();
@@ -331,7 +335,6 @@ void gLevelsetPoints::computeLS(std::vector<MVertex*> &vert, std::map<MVertex*,
}
evalRbfDer(0, 1, points, xyz_eval, surf, surf_eval);
for (int i = 0; i< vert.size(); i++){
- myMap[vert[i]] = surf_eval(i,0);
mapP[SPoint3(vert[i]->x(), vert[i]->y(),vert[i]->z())] = surf_eval(i,0);
}
}
diff --git a/contrib/DiscreteIntegration/DILevelset.h b/contrib/DiscreteIntegration/DILevelset.h
index 7c84f07abc..d4e25ad161 100644
--- a/contrib/DiscreteIntegration/DILevelset.h
+++ b/contrib/DiscreteIntegration/DILevelset.h
@@ -36,7 +36,7 @@ public:
gLevelset() : tag_(-1) {}
gLevelset(const gLevelset &);
virtual ~gLevelset(){}
- virtual gLevelset * clone() const{printf("Error virtual fct called gLevelset::clone()"); return 0;}
+ virtual gLevelset * clone() const{printf("Error virtual fct called gLevelset::clone()"); return 0;}
virtual double operator() (const double &x, const double &y, const double &z) const = 0;
// inline double operator () (const SPoint3 &p) const {return this->operator()(p.x(),p.y(),p.z());}
bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain > 0.;}
@@ -165,7 +165,7 @@ public:
virtual gLevelset * clone() const{return new gLevelsetPoints(*this);}
// return negative value inward and positive value outward
virtual double operator() (const double &x, const double &y, const double &z) const;
- void computeLS(std::vector<MVertex*> &vert, std::map<MVertex*, double> &myMap);
+ void computeLS(std::vector<MVertex*> &vert);
int type() const {return POINTS;}
};
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index dc7403f22b..c25e70f72e 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -2035,14 +2035,14 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
}
else{
- ll_subLines[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
ll->addLs(this, Lsi, iPrim++, nodeLs);
}
}
+ if(iPrim == 1) iPrim--;
+ ll_subLines[0]->addLs(this, RPN.back(), iPrim, nodeLs);
}
for(int l = 0; l < (int)ll_subLines.size(); l++) {
@@ -2190,14 +2190,14 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
}
}
else{
- tt_subTriangles[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
}
}
+ if(iPrim == 1) iPrim--;
+ tt_subTriangles[0]->addLs(this, RPN.back(), iPrim, nodeLs);
}
for(int q = 0; q < (int)tt_subQuads.size(); q++) {
@@ -2373,14 +2373,14 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
}
else {
- qq_subQuads[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
qq->addLs(this, Lsi, iPrim++, nodeLs);
}
}
+ if(iPrim == 1) iPrim--;
+ qq_subQuads[0]->addLs(this, RPN.back(), iPrim, nodeLs);
}
for(int q = 0; q < (int)qq_subQuads.size(); q++) {
@@ -2548,14 +2548,14 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
}
}
else {
- tt_subTetras[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
}
}
+ if(iPrim == 1) iPrim--;
+ tt_subTetras[0]->addLs(this, RPN.back(), iPrim, nodeLs);
}
for(int i = 0; i < (int)QError.size(); i++) {
@@ -2779,14 +2779,14 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
}
else {
- hh_subHexas[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
hh->addLs(this, Lsi, iPrim++, nodeLs);
}
}
+ if(iPrim == 1) iPrim--;
+ hh_subHexas[0]->addLs(this, RPN.back(), iPrim, nodeLs);
}
for(int i = 0; i < (int)QError.size(); i++) {
--
GitLab