diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 285642b4926fc564b6282a35617453bcbbb77f2d..866d52238506717bbd26c04cdd4c09ba88d66342 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -2,6 +2,9 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Gaetan Bricteux
#include <stdlib.h>
#include <sstream>
@@ -653,21 +656,21 @@ 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){
- for(int k = 0; k < copy->getNumEdges(); k++){
+ for(int k = 0; k < copy->getNumEdges(); k++){
MEdge me = copy->getEdge(k);
MVertex *v0 = me.getVertex(0);
MVertex *v1 = me.getVertex(1);
double val0 = verticesLs(0, v0->getIndex());
double val1 = verticesLs(0, v1->getIndex());
if (val0*val1 > 0.0 && val0 < 0.0) {
- MLine *lin = new MLine(v0, v1);
- getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
- int c = elements[1].count(gLsTag);
+ MLine *lin = new MLine(v0, v1);
+ getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
+ int c = elements[1].count(gLsTag);
int reg = getBorderTag(gLsTag, c, newElemTags[1][0], borderElemTags[0]);
int physTag = (!gePhysicals.size()) ? 0 : getBorderTag(gLsTag, c, newPhysTags[1][0], borderPhysTags[0]);
- elements[1][reg].push_back(lin);
- if(physTag) assignLsPhysical(GM, reg, 1, physicals, physTag, gLsTag);
- }
+ elements[1][reg].push_back(lin);
+ if(physTag) assignLsPhysical(GM, reg, 1, physicals, physTag, gLsTag);
+ }
}
}
else if (splitElem && copy->getDim()==3){
@@ -676,27 +679,26 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
bool sameSign = true;
double val0 = verticesLs(0, mf.getVertex(0)->getIndex());
for (int j= 1; j< mf.getNumVertices(); j++){
- double valj = verticesLs(0, mf.getVertex(j)->getIndex());
- if (val0*valj < 0.0){ sameSign = false; break;}
+ double valj = verticesLs(0, mf.getVertex(j)->getIndex());
+ if (val0*valj < 0.0){ sameSign = false; break;}
}
if (sameSign && val0 < 0.0) {
- getPhysicalTag(-1, gePhysicals, newPhysTags[2]);
- int c = elements[2].count(gLsTag) + elements[3].count(gLsTag) + elements[8].count(gLsTag);
+ getPhysicalTag(-1, gePhysicals, newPhysTags[2]);
+ int c = elements[2].count(gLsTag) + elements[3].count(gLsTag) + elements[8].count(gLsTag);
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));
+ if(mf.getNumVertices() == 3){
+ MTriangle *tri = new MTriangle(mf.getVertex(0), mf.getVertex(1), mf.getVertex(2));
elements[2][reg].push_back(tri);
- }
+ }
else if(mf.getNumVertices() == 4){
- MQuadrangle *quad = new MQuadrangle(mf.getVertex(0), mf.getVertex(1), mf.getVertex(2), mf.getVertex(3));
+ MQuadrangle *quad = new MQuadrangle(mf.getVertex(0), mf.getVertex(1), mf.getVertex(2), mf.getVertex(3));
elements[3][reg].push_back(quad);
- }
+ }
if(physTag) assignLsPhysical(GM, reg, 2, physicals, physTag, gLsTag);
}
- }
+ }
}
-
}
#if defined(HAVE_DINTEGRATION)
@@ -759,7 +761,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
MElement *copy = e->copy(vertexMap, newParents, newDomains);
MElement *parent = eParent ? copy->getParent() : copy;
- double **nodeLs = new double*[e->getNumPrimaryVertices()];
+ double **nodeLs = new double*[e->getNumVertices()];
switch (eType) {
case TYPE_TET :
@@ -774,7 +776,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
T.setPolynomialOrder(recur+1);
- for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
+ for(int i = 0; i < e->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, recur, nodeLs);
}
@@ -788,7 +791,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(6)->x(), e->getVertex(6)->y(), e->getVertex(6)->z(),
e->getVertex(7)->x(), e->getVertex(7)->y(), e->getVertex(7)->z());
H.setPolynomialOrder(recur+1);
- for(int i = 0; i < 8; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
+ for(int i = 0; i < e->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = H.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
hexas, tetras, quads, triangles, lines, recur, nodeLs);
}
@@ -797,7 +801,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- T1.setPolynomialOrder(recur+1);
+ //T1.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
nodeLs[3] = &verticesLs(0, e->getVertex(5)->getIndex());
bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
@@ -806,7 +810,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- T2.setPolynomialOrder(recur+1);
+ //T2.setPolynomialOrder(recur+1);
nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
nodeLs[1] = &verticesLs(0, e->getVertex(4)->getIndex());
nodeLs[2] = &verticesLs(0, e->getVertex(1)->getIndex());
@@ -817,7 +821,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- T3.setPolynomialOrder(recur+1);
+ //T3.setPolynomialOrder(recur+1);
for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i+2)->getIndex());
bool iC3 = T3.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, recur, nodeLs);
@@ -828,7 +832,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
- T1.setPolynomialOrder(recur+1);
+ //T1.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
nodeLs[3] = &verticesLs(0, e->getVertex(4)->getIndex());
bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
@@ -837,7 +841,7 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
- T2.setPolynomialOrder(recur+1);
+ //T2.setPolynomialOrder(recur+1);
nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i+1)->getIndex());
bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
@@ -852,7 +856,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
t->getVertex(3)->x(), t->getVertex(3)->y(), t->getVertex(3)->z());
Tet.setPolynomialOrder(recur+1);
- for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
+ for(int i = 0; i < t->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
bool iC = Tet.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, recur, nodeLs);
isCut = (isCut || iC);
@@ -997,7 +1002,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
T.setPolynomialOrder(recur+1);
- for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
+ for(int i = 0; i < e->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, recur, nodeLs);
}
@@ -1007,7 +1013,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
Q.setPolynomialOrder(recur+1);
- for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
+ for(int i = 0; i < e->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = Q.cut(RPN, ipV, ipS, cp, integOrder,integOrder,integOrder,
quads, triangles, lines, recur, nodeLs);
}
@@ -1018,7 +1025,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());
Tri.setPolynomialOrder(recur+1);
- for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
+ for(int i = 0; i < t->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
bool iC = Tri.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, recur, nodeLs);
isCut = (isCut || iC);
@@ -1203,7 +1211,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
DI_Line L(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z());
L.setPolynomialOrder(recur+1);
- for(int i = 0; i < 2; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
+ for(int i = 0; i < e->getNumVertices(); i++)
+ nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = L.cut(RPN, ipV, cp, integOrder, lines, recur, nodeLs);
if(isCut) {
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 7c0a7c7f5d73dbe244db835268782a4af812508c..ae33a5ef76df3f3d5d44f5c445cdffddba1e9837 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -2,6 +2,9 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Gaetan Bricteux
#ifndef _MELEMENTCUT_H_
#define _MELEMENTCUT_H_
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 4641caf25705aed9a48ee006655fd406f0805180..ccd5058f1d90bedb08e7deecb6171bcacc6193eb 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1,3 +1,12 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Gaetan Bricteux
+//
+
#include "gmshLevelset.h"
#include <queue>
#include <stack>
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index de05cd67edf988883c5cb9f846c74bc599144b6f..7cdd4b85706ebcfc0b4034f9cc654f43157d488e 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -2,6 +2,9 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Gaetan Bricteux
#ifndef _GMSH_LEVELSET_H_
#define _GMSH_LEVELSET_H_
@@ -122,7 +125,7 @@ protected:
public:
gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag=1) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
virtual double operator () (const double x, const double y, const double z) const
- {return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;}
+ {return (sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - abs(r)) * r / abs(r);}
int type() const {return SPHERE;}
};
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 28374f7c0b81c1ec6c18394033b8cd8a617ee38e..18b2caf491d553d135483d7afdc16629baa8c3fc 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -1399,7 +1399,7 @@ const polynomialBasis *polynomialBases::find(int tag)
case MSH_LIN_1 : F.parentType = TYPE_LIN; F.order = 0; F.serendip = false; break;
case MSH_LIN_2 : F.parentType = TYPE_LIN; F.order = 1; F.serendip = false; break;
case MSH_LIN_3 : F.parentType = TYPE_LIN; F.order = 2; F.serendip = false; break;
- case MSH_LIN_4 : printf("line 4 \n"); F.parentType = TYPE_LIN; F.order = 3; F.serendip = false; break;
+ case MSH_LIN_4 : F.parentType = TYPE_LIN; F.order = 3; F.serendip = false; break;
case MSH_LIN_5 : F.parentType = TYPE_LIN; F.order = 4; F.serendip = false; break;
case MSH_LIN_6 : F.parentType = TYPE_LIN; F.order = 5; F.serendip = false; break;
case MSH_LIN_7 : F.parentType = TYPE_LIN; F.order = 6; F.serendip = false; break;
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 87a64b4c1a5ea6c716e22ae5321a3927ca20e744..b3bd8ca6a5682646379a193eb491aa74b16fef59 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -93,7 +93,7 @@ class polynomialBasis
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
int numFaces;
- inline int getNumShapeFunctions() const {return coefficients.size1();};
+ inline int getNumShapeFunctions() const {return coefficients.size1();}
// for a given face/edge, with both a sign and a rotation, give an
// ordered list of nodes on this face/edge
inline const std::vector<int> &getClosure(int id) const
diff --git a/benchmarks/levelset/square.geo b/benchmarks/levelset/square.geo
index 413dc38e2d855defb1a20ebe1cb56b9aa89bacb8..fe6e8719bffa3bf1122f256f3dcc37259d7e2a73 100644
--- a/benchmarks/levelset/square.geo
+++ b/benchmarks/levelset/square.geo
@@ -9,6 +9,7 @@
Mesh.Algorithm=1;
Mesh.CharacteristicLengthFactor=1.5;
+Mesh.ElementOrder=2;
lc=0.1;
Point(1) = {0.1, 0.1, 0.1}; //,lc};
@@ -25,12 +26,10 @@ Plane Surface(10) = {5};
Mesh 2;
Levelset Plane (1) = {0,1,0,-7};
+Levelset Sphere (2) = {{0,0,0},-10};
//Levelset Point (1) = {{0.1, 2, 0},{3,2,0},{9,2,0},{5,2,0}, {0.1, 2.2, -1},{3,2.5,-1},{9,2,-1},{5,2,-1}, {0.2, 2, 1},{3,2,1},{9,2,1},{5,2,1}};
-
-Levelset CutMeshTri {1};
+Levelset Intersection (3) = {1,2};
+Levelset CutMeshTri {3};
Print "square_cut.msh";
-
-
-
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index a4ad62b753431f1526a7487a67d2c1cdfbacc3eb..a9f8ed9c36855ab477e3874636c596a4ddb3466b 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -1,3 +1,5 @@
+// Author : Gaetan Bricteux
+
#ifndef INTEGRATION3D_CC
#define INTEGRATION3D_CC
@@ -580,11 +582,10 @@ bool DI_PointLessThan::operator()(const DI_Point *p1, const DI_Point *p2) const
// DI_Element methods -----------------------------------------------------------------------------
DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.getPolynomialOrder()),
integral_(cp.integral()) {
- //printf("constructeur de copie d'element %d : ",cp.type()); cout << &cp << "->" << this << endl;
pts_ = new DI_Point [cp.nbVert()];
for(int i = 0; i < cp.nbVert(); i++)
pts_[i] = DI_Point(*cp.pt(i));
- if(cp.getPolynomialOrder() > 1) {
+ if(cp.nbMid()) {
mid_ = new DI_Point [cp.nbMid()];
for(int i = 0; i < cp.nbMid(); i++)
mid_[i] = DI_Point(*cp.mid(i));
@@ -603,12 +604,13 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
for(int i = 0; i < nbVert(); i++) {
pts_[i] = DI_Point(*rhs.pt(i));
}
- if(rhs.nbMid() > 0) {
- delete mid_;
+ if(rhs.nbMid()) {
+ if(mid_) delete [] mid_;
mid_ = new DI_Point[rhs.nbMid()];
+ for(int i = 0; i < rhs.nbMid(); i++)
+ mid_[i] = DI_Point(*(rhs.mid(i)));
}
- for(int i = 0; i < rhs.nbMid(); i++)
- mid_[i] = DI_Point(*rhs.mid(i));
+ else mid_ = NULL;
polOrder_ = rhs.getPolynomialOrder();
integral_ = rhs.integral();
lsTag_ = rhs.lsTag();
@@ -632,7 +634,10 @@ void DI_Element::getGradShapeFunctions(double u, double v, double w, double s[][
}
void DI_Element::setPolynomialOrder (int o) {
if(polOrder_ == o) return;
- if (mid_) delete [] mid_;
+ if (mid_){
+ delete [] mid_;
+ mid_ = NULL;
+ }
polOrder_ = o;
switch (o) {
case 1 :
@@ -649,13 +654,15 @@ void DI_Element::setPolynomialOrder (int o) {
}
return;
default :
- return;
- //printf("Order %d line function space not implemented ", o); print();
+ printf("Order %d function space not implemented ", o); print();
}
}
void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
if(polOrder_ == o) return;
- if (mid_) delete [] mid_;
+ if (mid_){
+ delete [] mid_;
+ mid_ = NULL;
+ }
polOrder_ = o;
switch (o) {
case 1 :
@@ -672,38 +679,25 @@ void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vect
}
return;
default :
- printf("Order %d line function space not implemented ", o); print();
+ printf("Order %d function space not implemented ", o); print();
}
}
void DI_Element::addLs (const double *ls) {
- for(int i = 0; i < nbVert(); i++)
- pts_[i].addLs(ls[i]);
- for(int i = 0; i < nbMid(); i++)
- mid_[i].addLs(ls[nbVert()+i]);
+ for(int i = 0; i < nbVert() + nbMid(); i++)
+ pt(i)->addLs(ls[i]);
}
void DI_Element::addLs (const DI_Element *e) {
if(e->sizeLs() < 1) return;
- for(int i = 0; i < nbVert(); i++)
- pts_[i].addLs(e);
- for(int i = 0; i < nbMid(); i++)
- mid_[i].addLs(e);
+ for(int i = 0; i < nbVert() + nbMid(); i++)
+ pt(i)->addLs(e);
}
void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls) {
if(type() != e->type()) {
printf("Error : addLs with element of different type\n");
}
- for(int j = 0; j < nbVert(); ++j) {
- double ls = (*Ls)(e->x(j), e->y(j), e->z(j));
- pts_[j].addLs(ls);
- }
- for(int j = 0; j < nbMid(); ++j) {
- std::vector<int> s(nbVert()); int n;
- e->midV(j, &s[0], n);
- double xc = 0, yc = 0, zc = 0;
- for(int k = 0; k < n; k++){
- xc += e->x(s[k]); yc += e->y(s[k]); zc += e->z(s[k]); }
- double ls = (*Ls)(xc/n, yc/n, zc/n);
- mid_[j].addLs(ls);
+ for(int i = 0; i < nbVert() + nbMid(); ++i) {
+ double ls = (*Ls)(e->x(i), e->y(i), e->z(i));
+ pt(i)->addLs(ls);
}
}
void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls, int iLs, double **nodeLs) {
@@ -711,31 +705,18 @@ void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls, int iLs, doubl
addLs(e, Ls);
return;
}
- for(int i = 0; i < nbVert(); i++)
- pts_[i].addLs(nodeLs[i][iLs]);
- for(int j = 0; j < nbMid(); ++j) {
- std::vector<int> s(nbVert()); int n;
- e->midV(j, &s[0], n);
- double xc = 0, yc = 0, zc = 0;
- for(int k = 0; k < n; k++){
- xc += e->x(s[k]); yc += e->y(s[k]); zc += e->z(s[k]); }
- double ls = (*Ls)(xc/n, yc/n, zc/n);
- mid_[j].addLs(ls);
- }
+ for(int i = 0; i < nbVert() + nbMid(); i++)
+ pt(i)->addLs(nodeLs[i][iLs]);
}
void DI_Element::chooseLs (const gLevelset *Lsi) {
if(sizeLs() < 2)
printf("chooseLs with element ls size < 2 : typeEl=%d\n", type());
- for(int i = 0; i < nbVert(); i++)
- pts_[i].chooseLs(Lsi);
- for(int i = 0; i < nbMid(); i++)
- mid_[i].chooseLs(Lsi);
+ for(int i = 0; i < nbVert() + nbMid(); i++)
+ pt(i)->chooseLs(Lsi);
}
void DI_Element::clearLs() {
- for(int i = 0; i < nbVert(); i++)
- pts_[i].clearLs();
- for(int i = 0; i < nbMid(); i++)
- mid_[i].clearLs();
+ for(int i = 0; i < nbVert() + nbMid(); i++)
+ pt(i)->clearLs();
}
bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
const DI_Element *e, const std::vector<gLevelset *> &RPNi) {
@@ -822,13 +803,9 @@ void DI_Element::mappingCP (DI_CuttingPoint *cp) const {
}
void DI_Element::mappingEl (DI_Element *el) const {
double s[3];
- for (int i = 0; i < el->nbVert(); i++) {
- evalC(el->x(i), el->y(i), el->z(i), s);
- el->pts_[i].move(s[0], s[1], s[2]);
- }
- for(int i = el->nbVert(); i < el->nbVert() + el->nbMid(); i++) {
+ for (int i = 0; i < el->nbVert() + el->nbMid(); i++) {
evalC(el->x(i), el->y(i), el->z(i), s);
- el->mid_[i - el->nbVert()].move(s[0], s[1], s[2]);
+ el->pt(i)->move(s[0], s[1], s[2]);
}
el->computeIntegral();
}
@@ -873,19 +850,20 @@ void DI_Element::evalC (const double u, const double v, const double w, double *
int nbV = nbVert() + nbMid();
std::vector<double> s(nbV);
ev[0] = 0; ev[1] = 0; ev[2] = 0;
- getShapeFunctions (u, v, w, &s[0], order);
+ getShapeFunctions (u, v, w, &s[0], order);
for(int i = 0; i < nbV; i++){
- ev[0] += x(i) * s[i];
- ev[1] += y(i) * s[i];
- ev[2] += z(i) * s[i];
+ ev[0] += x(i) * s[i];
+ ev[1] += y(i) * s[i];
+ ev[2] += z(i) * s[i];
}
}
double DI_Element::evalLs (const double u, const double v, const double w, int iLs, int order) const{
+ int nbV = nbVert() + nbMid();
if(iLs == -1) iLs = sizeLs() - 1; //last ls value
double vls = 0;
- std::vector<double> s(nbVert() + nbMid());
+ std::vector<double> s(nbV);
getShapeFunctions (u, v, w, &s[0], order);
- for(int i = 0; i < nbVert() + nbMid(); i++)
+ for(int i = 0; i < nbV; i++)
vls += ls(i, iLs) * s[i];
return adjustLs(vls);
}
@@ -905,11 +883,12 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
s[i] = new double[3];//s[i] = new double[getDim()];
}*/
typedef double d3[3];
- d3 *s=new d3[nbVert()+ nbMid()];
+ int nbV = nbVert() + nbMid();
+ d3 *s = new d3[nbV];
getGradShapeFunctions(u, v, w, s);
switch(getDim()){
case 3 :
- for(int i = 0; i < nbVert() + nbMid(); i++) {
+ for(int i = 0; i < nbV; i++) {
J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
J[1][0] += x(i) * s[i][1]; J[1][1] += y(i) * s[i][1]; J[1][2] += z(i) * s[i][1];
J[2][0] += x(i) * s[i][2]; J[2][1] += y(i) * s[i][2]; J[2][2] += z(i) * s[i][2];
@@ -918,7 +897,7 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
delete [] s;
return det3(J[0][0], J[0][1], J[0][2], J[1][0], J[1][1], J[1][2], J[2][0], J[2][1], J[2][2]);
case 2 :
- for(int i = 0; i < nbVert() + nbMid(); i++) {
+ for(int i = 0; i < nbV; i++) {
J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
J[1][0] += x(i) * s[i][1]; J[1][1] += y(i) * s[i][1]; J[1][2] += z(i) * s[i][1];
// delete [] s[i];
@@ -928,7 +907,7 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
sq2(J[0][2] * J[1][0] - J[0][0] * J[1][2]) +
sq2(J[0][1] * J[1][2] - J[0][2] * J[1][1])); // allways > 0 => does'nt allow to check if twist !!!!
case 1 :
- for(int i = 0; i < nbVert() + nbMid(); i++) {
+ for(int i = 0; i < nbV; i++) {
J[0][0] += x(i) * s[i][0]; J[0][1] += y(i) * s[i][0]; J[0][2] += z(i) * s[i][0];
// delete [] s[i];
}
@@ -1110,7 +1089,7 @@ const polynomialBasis* DI_Triangle::getFunctionSpace(int o) const
case 9: return polynomialBases::find(MSH_TRI_55);
case 10: return polynomialBases::find(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
- }
+ }
}
void DI_Triangle::computeIntegral() {
integral_ = TriSurf(pt(0), pt(1), pt(2));
@@ -1805,7 +1784,6 @@ bool DI_Line::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoin
DI_Point::Container &cp, const int polynomialOrder,
std::vector<DI_Line *> &lines, int recurLevel, double **nodeLs) const
{
- //printf(" "); print();
std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Line *ll = new DI_Line(-1, 0, 0, 1, 0, 0, 2.); //reference line
@@ -1902,7 +1880,6 @@ bool DI_Triangle::cut (std::vector<gLevelset *> &RPN, std::vector<DI_Integration
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
{
- //printf(" ");print();
std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Triangle *tt = new DI_Triangle(0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5); //reference triangle
@@ -1920,7 +1897,7 @@ bool DI_Triangle::cut (std::vector<gLevelset *> &RPN, std::vector<DI_Integration
int iPrim = 0;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
- gLevelset *Lsi = RPN[l]; //printf("LS(0,0)=%g LS(1,1)=%g\n",(*Lsi)(0,0,0),(*Lsi)(1,1,0));
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
@@ -2081,7 +2058,6 @@ bool DI_Quad::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoin
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
{
- //printf(" "); printls();
std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Quad *qq = new DI_Quad(-1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 4.); // reference quad
@@ -2275,7 +2251,6 @@ bool DI_Tetra::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoi
std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
std::vector<DI_Triangle *> &surfTriangles, int recurLevel, double **nodeLs) const
{
- //printf(" ");print();
std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Tetra *tt = new DI_Tetra(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1./6.); // reference tetra
@@ -2448,7 +2423,6 @@ bool DI_Hexa::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoin
std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
std::vector<DI_Line *> &frontLines, int recurLevel, double **nodeLs) const
{
- //printf(" "); print();
std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
// reference hexa
@@ -2586,7 +2560,6 @@ bool DI_Hexa::cut (std::vector<gLevelset *> &RPN, std::vector<DI_IntegrationPoin
for(int p = 0; p < (int)hh_cp.size(); p++)
hh_cp[p]->chooseLs(Lsi);
}
- //printf("Lsi=%d sizeLs=%d nbH=%d nbT=%d nbQ=%d nbTr=%d sizeLs(hh_cp)=%d\n",Lsi->type(),(hh_subHexas.size()>0)?hh_subHexas[0].sizeLs():hh_subTetras[0].sizeLs(), hh_subHexas.size(),hh_subTetras.size(),hh_surfQuads.size(),hh_surfTriangles.size(),(hh_cp.size() > 0)?hh_cp[0].sizeLs():0);
}
}
else {
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 0e9ab75693097a6145d2e656b81fffb30a9bddd6..3e5f22156cf2bce5505498a0ee5988c9f85039cd 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -1,3 +1,5 @@
+// Author : Gaetan Bricteux
+
#ifndef _INTEGRATION3D_H_
#define _INTEGRATION3D_H_
@@ -322,17 +324,17 @@ class DI_Element
void getCuttingPoints (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_CuttingPoint*> &cp) const;
// return the ith point
- virtual DI_Point* pt (int i) const = 0;
+ virtual DI_Point* pt(int i) const {return (i < nbVert()) ? &(pts_[i]) : &(mid_[i - nbVert()]);}
// return the ith middle point
- inline DI_Point* mid (int i) const {return &mid_[i];}
+ inline DI_Point* mid(int i) const {return mid_ ? &(mid_[i]) : NULL;}
// return the coordinates of the ith point
- virtual double x (int i) const = 0;
- virtual double y (int i) const = 0;
- virtual double z (int i) const = 0;
+ virtual double x(int i) const {return pt(i)->x();}
+ virtual double y(int i) const {return pt(i)->y();}
+ virtual double z(int i) const {return pt(i)->z();}
// return the last levelset value of the ith point
- virtual double ls (int i) const = 0;
+ virtual double ls(int i) const {return pt(i)->ls();}
// return the jth levelset value of the ith point
- virtual double ls (int i, int j) const = 0;
+ virtual double ls(int i, int j) const {return pt(i)->ls(j);}
// return the number of levelset values of the points
inline int sizeLs() const {return pts_[0].sizeLs();}
// return the interpolating nodal shape functions evaluated at point (u,v,w)
@@ -405,27 +407,16 @@ class DI_Line : public DI_Element
inline int type() const {return DI_LIN;}
inline int getDim() const {return 1;}
inline int nbVert() const {return 2;}
- inline int nbMid() const {
- switch(polOrder_) {
- case 1 : return 0;
- case 2 : return 1;
- default : return 0;
- }
- }
+ inline int nbMid() const {return polOrder_ - 1;}
inline int nbEdg() const {return 1;}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 2.;}
- inline DI_Point* pt (int i) const {return i < 2 ? &pts_[i] : &mid_[i - 2];}
- inline double x (int i) const {return i < 2 ? pts_[i].x() : mid_[i - 2].x();}
- inline double y (int i) const {return i < 2 ? pts_[i].y() : mid_[i - 2].y();}
- inline double z (int i) const {return i < 2 ? pts_[i].z() : mid_[i - 2].z();}
- inline double ls (int i) const {return i < 2 ? pts_[i].ls() : mid_[i - 2].ls();}
- inline double ls (int i, int j) const {return i < 2 ? pts_[i].ls(j) : mid_[i - 2].ls(j);}
void getRefIntegrationPoints (const int polynomialOrder,
std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const {
- s1 = 0; s2 = 1;}
+ s1 = 0; s2 = 1;
+ }
void midV (const int e, int *s, int &n) const {
s[0] = 0; s[1] = 1; n = 2;
}
@@ -482,22 +473,12 @@ class DI_Triangle : public DI_Element
inline int getDim() const {return 2;}
inline int nbVert() const {return 3;}
inline int nbMid() const {
- switch(polOrder_) {
- case 1 : return 0;
- case 2 : return 3;
- default : return 0;
- }
+ return 0.5 * (polOrder_ - 1) * (polOrder_ + 4);
}
inline int nbEdg() const {return 3;}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 0.5;}
- inline DI_Point* pt (int i) const {return i < 3 ? &pts_[i] : &mid_[i - 3];}
- inline double x (int i) const {return i < 3 ? pts_[i].x() : mid_[i - 3].x();}
- inline double y (int i) const {return i < 3 ? pts_[i].y() : mid_[i - 3].y();}
- inline double z (int i) const {return i < 3 ? pts_[i].z() : mid_[i - 3].z();}
- inline double ls (int i) const {return i < 3 ? pts_[i].ls() : mid_[i - 3].ls();}
- inline double ls (int i, int j) const {return i < 3 ? pts_[i].ls(j) : mid_[i - 3].ls(j);}
void getRefIntegrationPoints (const int polynomialOrder,
std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const {
@@ -576,22 +557,12 @@ class DI_Quad : public DI_Element
inline int getDim() const {return 2;}
inline int nbVert() const {return 4;}
inline int nbMid() const {
- switch(polOrder_) {
- case 1 : return 0;
- case 2 : return 5;
- default : return 0;
- }
+ return (polOrder_ - 1) * (polOrder_ + 3);
}
inline int nbEdg() const {return 4;}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 4.;}
- inline DI_Point* pt (int i) const {return i < 4 ? &pts_[i] : &mid_[i - 4];}
- inline double x (int i) const {return i < 4 ? pts_[i].x() : mid_[i - 4].x();}
- inline double y (int i) const {return i < 4 ? pts_[i].y() : mid_[i - 4].y();}
- inline double z (int i) const {return i < 4 ? pts_[i].z() : mid_[i - 4].z();}
- inline double ls (int i) const {return i < 4 ? pts_[i].ls() : mid_[i - 4].ls();}
- inline double ls (int i, int j) const {return i < 4 ? pts_[i].ls(j) : mid_[i - 4].ls(j);}
void getRefIntegrationPoints (const int polynomialOrder,
std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const{
@@ -668,22 +639,12 @@ class DI_Tetra : public DI_Element
inline int getDim() const {return 3;}
inline int nbVert() const {return 4;}
inline int nbMid() const {
- switch(polOrder_) {
- case 1 : return 0;
- case 2 : return 6;
- default : return 0;
- }
+ return (polOrder_ + 1) * (polOrder_ + 2) * (polOrder_ + 3) / 6 - 4;
}
inline int nbEdg() const {return 6;}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 1. / 6.;}
- inline DI_Point* pt (int i) const {return i < 4 ? &pts_[i] : &mid_[i - 4];}
- inline double x (int i) const {return i < 4 ? pts_[i].x() : mid_[i - 4].x();}
- inline double y (int i) const {return i < 4 ? pts_[i].y() : mid_[i - 4].y();}
- inline double z (int i) const {return i < 4 ? pts_[i].z() : mid_[i - 4].z();}
- inline double ls (int i) const {return i < 4 ? pts_[i].ls() : mid_[i - 4].ls();}
- inline double ls (int i, int j) const {return i < 4 ? pts_[i].ls(j) : mid_[i - 4].ls(j);}
inline void getRefIntegrationPoints (const int polynomialOrder,
std::vector<DI_IntegrationPoint *> &ip) const;
inline void vert(const int edge, int &s1, int &s2) const {
@@ -780,22 +741,12 @@ class DI_Hexa : public DI_Element
inline int getDim() const {return 3;}
inline int nbVert() const {return 8;}
inline int nbMid() const {
- switch(polOrder_) {
- case 1 : return 0;
- case 2 : return 19;
- default : return 0;
- }
+ return (polOrder_ + 1) * (polOrder_ + 1) * (polOrder_ + 1) - 8;
}
inline int nbEdg() const {return 12;}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
void computeIntegral();
inline double refIntegral() const {return 8.;}
- inline DI_Point* pt (int i) const {return i < 8 ? &pts_[i] : &mid_[i - 8];}
- inline double x (int i) const {return i < 8 ? pts_[i].x() : mid_[i - 8].x();}
- inline double y (int i) const {return i < 8 ? pts_[i].y() : mid_[i - 8].y();}
- inline double z (int i) const {return i < 8 ? pts_[i].z() : mid_[i - 8].z();}
- inline double ls (int i) const {return i < 8 ? pts_[i].ls() : mid_[i - 8].ls();}
- inline double ls (int i, int j) const {return i < 8 ? pts_[i].ls(j) : mid_[i - 8].ls(j);}
void getRefIntegrationPoints (const int polynomialOrder,
std::vector<DI_IntegrationPoint *> &ip) const;
inline void vert(const int edge, int &s1, int &s2) const {
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index 39a897128c6c2d6f1919587a0b79f3ec1fc1c2a3..0c0565265ae1e708785d70395d4480cb33274dde 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -1,3 +1,5 @@
+//Author : Gaetan Bricteux
+
#ifndef RECURCUT_CC
#define RECURCUT_CC
@@ -232,13 +234,13 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<gLevel
elem->addLs(e, Lsi, iPrim++, nodeLs);
for(unsigned int i = 0; i < cp.size(); i++)
cp[i]->addLs(elem);
- if (re->super) re->el->addLs(elem);
+ if(re->super) re->el->addLs(elem);
re->el->getCuttingPoints(elem, RPNi, cp);
}
else {
for(unsigned int p = 0; p < cp.size(); p++)
cp[p]->chooseLs(Lsi);
- if (re->super) re->el->chooseLs(Lsi);
+ if(re->super) re->el->chooseLs(Lsi);
}
}
re->root()->el->clearLs();
diff --git a/contrib/DiscreteIntegration/recurCut.h b/contrib/DiscreteIntegration/recurCut.h
index 569ec66ce90af1777aeba5b6bcd5e64e48319abd..0bb3799810bff3ba114fdda8c3608d12517c3602 100644
--- a/contrib/DiscreteIntegration/recurCut.h
+++ b/contrib/DiscreteIntegration/recurCut.h
@@ -1,3 +1,5 @@
+//Author : Gaetan Bricteux
+
#ifndef _RECUR_H_
#define _RECUR_H_