From aedf3e88d2862bc0e408040a266604ffbada7829 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 21 Oct 2011 10:45:26 +0000
Subject: [PATCH] quad cut for gLevelset
---
Geo/MElementCut.cpp | 51 ++++++++++---------
contrib/DiscreteIntegration/Integration3D.cpp | 10 ++--
2 files changed, 33 insertions(+), 28 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index ec300ed0ff..676cee9c34 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -950,30 +950,31 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
else
poly[1].push_back(mt);
}
- // //if quads
- // for (unsigned int i = nbQ; i < quads.size(); i++){
- // MVertex *mv[4] = {NULL, NULL, NULL, NULL};
- // for(int j = 0; j < 4; j++){
- // int numV = getElementVertexNum(triangles[i]->pt(j), e);
- // if(numV == -1) {
- // MVertex *newv = new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
- // std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
- // mv[j] = *(it.first);
- // if (!it.second) newv->deleteLast();
- // }
- // else {
- // std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
- // if(it == vertexMap.end()) {
- // mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
- // quads[i]->z(j), 0, numV);
- // vertexMap[numV] = mv[j];
- // }
- // else mv[j] = it->second;
- // }
- // }
- // MQuadrangle *mq = new MQuadrangle(mv[0], mv[1], mv[2], mv[3], 0, 0);
- // elements[3][reg].push_back(mq);
- // }
+ //if quads
+ for (unsigned int i = nbQ; i < quads.size(); i++){
+ MVertex *mv[4] = {NULL, NULL, NULL, NULL};
+ for(int j = 0; j < 4; j++){
+ int numV = getElementVertexNum(quads[i]->pt(j), e);
+ if(numV == -1) {
+ MVertex *newv = new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
+ std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+ mv[j] = *(it.first);
+ if (!it.second) newv->deleteLast();
+ }
+ else {
+ std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
+ if(it == vertexMap.end()) {
+ mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
+ quads[i]->z(j), 0, numV);
+ vertexMap[numV] = mv[j];
+ }
+ else mv[j] = it->second;
+ }
+ }
+ MQuadrangle *mq = new MQuadrangle(mv[0], mv[1], mv[2], mv[3], 0, 0);
+ int reg = getElementaryTag(quads[i]->lsTag(), elementary, newElemTags[2]);
+ elements[3][reg].push_back(mq);
+ }
bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
@@ -1337,7 +1338,9 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
}
for(DI_Point::Container::iterator it = cp.begin(); it != cp.end(); it++) delete *it;
+ printf("coucou delet dim =%d\n", gmEntities[i]->dim());
for(unsigned int k = 0; k < lines.size(); k++) delete lines[k];
+ printf("coucou delet \n");
for(unsigned int k = 0; k < triangles.size(); k++) delete triangles[k];
for(unsigned int k = 0; k < quads.size(); k++) delete quads[k];
for(unsigned int k = 0; k < tetras.size(); k++) delete tetras[k];
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index b8d68b5b64..76fb9977c6 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -5,7 +5,7 @@
#include "recurCut.h"
#include "Gauss.h"
-#define ZERO_LS_TOL 1.e-7
+#define ZERO_LS_TOL 1.e-3
#define EQUALITY_TOL 1.e-15
// cross product
@@ -729,7 +729,7 @@ bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
double sideLength = distance(pt(s1), pt(s2));
if (dist / sideLength < 1.e-5) {
if(order0 == 1) setPolynomialOrder(1);
- printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quad edge\n", dist, sideLength, dist/sideLength);
+ printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quadratic edge\n", dist, sideLength, dist/sideLength);
return true; // do not add the quadratic edge if xm is very close from the middle of the edge
}
DI_Point *p0 = &mid_[edge];
@@ -1447,7 +1447,8 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<gLevelset *>
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
{
- bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
+ //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
+ bool quadT = false;
int LStag = RPNi.back()->getTag();
int nbZe = 0, ze[3];
@@ -1565,7 +1566,8 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &R
std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
std::vector<DI_CuttingPoint *> &cuttingPoints, std::vector<DI_QualError *> &QError) const
{
- bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
+ //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
+ bool quadT = false;
int tag = RPNi.back()->getTag();
int nbZe = 0, ze[4];
--
GitLab