diff --git a/CMakeLists.txt b/CMakeLists.txt
index f7cf85fa1075d1073c12df50872d0aee6d461f20..7abf939ed1ea98fb1569026e3c471a09c12ae30e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -69,7 +69,8 @@ set(GMSH_SHORT_LICENSE "GNU General Public License")
set(GMSH_API
${CMAKE_CURRENT_BINARY_DIR}/Common/GmshConfig.h
${CMAKE_CURRENT_BINARY_DIR}/Common/GmshVersion.h
- Common/Gmsh.h Common/Context.h Common/GmshDefines.h Common/GmshMessage.h Common/VertexArray.h
+ Common/Gmsh.h Common/Context.h Common/GmshDefines.h Common/GmshMessage.h
+ Common/VertexArray.h Common/Octree.h Common/OctreeInternals.h
Numeric/Numeric.h Numeric/Gauss.h Numeric/polynomialBasis.h
Numeric/JacobianBasis.h Numeric/fullMatrix.h
Numeric/simpleFunction.h Numeric/cartesian.h
@@ -83,18 +84,20 @@ set(GMSH_API
Geo/discreteFace.h Geo/discreteRegion.h Geo/SPoint2.h Geo/SPoint3.h
Geo/SVector3.h Geo/STensor3.h Geo/SBoundingBox3d.h Geo/Pair.h Geo/Range.h
Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h
- Geo/Homology.h Geo/partitionEdge.h Geo/CGNSOptions.h
+ Geo/Homology.h Geo/partitionEdge.h Geo/CGNSOptions.h Geo/gmshLevelset.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h Mesh/meshPartition.h
Mesh/meshGFaceDelaunayInsertion.h Mesh/Levy3D.h Mesh/meshPartitionOptions.h
+ Numeric/mathEvaluator.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Solver/crossConfTerm.h Solver/orthogonalTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
Solver/linearSystemFull.h Solver/elasticitySolver.h
- Post/PView.h Post/PViewData.h Plugin/PluginManager.h
+ Post/PView.h Post/PViewData.h Plugin/PluginManager.h Post/OctreePost.h
Graphics/drawContext.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
contrib/kbipack/gmp_blas.h contrib/kbipack/mpz.h
- contrib/DiscreteIntegration/DILevelset.h)
+ contrib/DiscreteIntegration/Integration3D.h
+ contrib/MathEx/mathex.h)
execute_process(COMMAND date "+%Y%m%d" OUTPUT_VARIABLE DATE
OUTPUT_STRIP_TRAILING_WHITESPACE)
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index a651e1e14169d5fe9864d30f7dac7a291b145496..8ff7dc399288f3e56562864fad0d0de49180ba5b 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -213,12 +213,14 @@ static void getDomains(int dom1Num, int dom2Num, int type,
MElement *doms[])
{
int domNums[2] = {dom1Num, dom2Num};
+ int nbD = (dom1Num == 0 || dom2Num == 0) ? 1 : 2;
+ if(dom1Num == 0) domNums[0] = dom2Num;
switch(type){
case MSH_LIN_B :
- getElementsByNum(domNums, elements[8], false, doms, 2);
+ getElementsByNum(domNums, elements[8], false, doms, nbD);
return;
case MSH_TRI_B : case MSH_POLYG_B :
- getElementsByNum(domNums, elements[9], false, doms, 2);
+ getElementsByNum(domNums, elements[9], false, doms, nbD);
return;
}
}
@@ -461,11 +463,11 @@ int GModel::readMSH(const std::string &name)
if(dom1) {
getDomains(dom1, dom2, type, elements, doms);
if(!doms[0]) Msg::Error("Domain element %d not found for element %d", dom1, num);
- if(!doms[1]) Msg::Error("Domain element %d not found for element %d", dom2, num);
+ if(dom2 && !doms[1]) Msg::Error("Domain element %d not found for element %d", dom2, num);
}
MElement *e = createElementMSH(this, num, type, physical, elementary,
partition, vertices, elements, physicals,
- own, p, doms[1], doms[0]);
+ own, p, doms[0], doms[1]);
for(unsigned int j = 0; j < ghosts.size(); j++)
_ghostCells.insert(std::pair<MElement*, short>(e, ghosts[j]));
if(numElements > 100000)
@@ -702,8 +704,8 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
version = 1.0;
// if there are no physicals we save all the elements
- if(noPhysicalGroups()) saveAll = true;
-
+ if(noPhysicalGroups()) saveAll = true;
+
// get the number of vertices and index the vertices in a continuous
// sequence
int numVertices = indexMeshVertices(saveAll, saveSinglePartition);
@@ -767,7 +769,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
_elementIndexCache.clear();
//parents
- if ( !CTX::instance()->mesh.saveTri){
+ if (!CTX::instance()->mesh.saveTri){
for(eiter it = firstEdge(); it != lastEdge(); ++it)
for(unsigned int i = 0; i < (*it)->lines.size(); i++)
if((*it)->lines[i]->ownsParent())
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 3f2a3ca74463c25615973062421cccdc2f973873..ea20fce584a34703b7179d6267784d7e368d9b9b 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -769,7 +769,8 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
- p1 = new MPolyhedron(poly[0], ++numEle, ePart, own, parent);
+ int n = (e->getParent()) ? e->getNum() : ++numEle;
+ p1 = new MPolyhedron(poly[0], n, ePart, own, parent);
own = false;
int reg = getElementaryTag(-1, elementary, newElemTags[3]);
std::vector<int> phys;
@@ -778,7 +779,8 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
assignPhysicals(GM, phys, reg, 3, physicals);
}
if(poly[1].size()) {
- p2 = new MPolyhedron(poly[1], ++numEle, ePart, own, parent);
+ int n = (e->getParent() && poly[0].size() == 0) ? e->getNum() : ++numEle;
+ p2 = new MPolyhedron(poly[1], n, ePart, own, parent);
elements[9][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 3, physicals);
}
@@ -843,7 +845,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
MTriangle *tri;
- if(p1 || p2) tri = new MTriangleBorder(mv[0], mv[1], mv[2], ++numEle, ePart, p1, p2);
+ if(p1 || p2){
+ if(!p1) tri = new MTriangleBorder(mv[0], mv[1], mv[2], ++numEle, ePart, p2, p1);
+ else tri = new MTriangleBorder(mv[0], mv[1], mv[2], ++numEle, ePart, p1, p2);
+ }
else tri = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
int lsT = triangles[i]->lsTag();
int c = elements[2].count(lsT) + elements[3].count(lsT) + elements[8].count(lsT);
@@ -927,10 +932,11 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
+ int n = (e->getParent()) ? e->getNum() : ++numEle;
if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
- p1 = new MPolygonBorder(poly[0], ++numEle, ePart, own, parent,
+ p1 = new MPolygonBorder(poly[0], n, ePart, own, parent,
copy->getDomain(0), copy->getDomain(1));
- else p1 = new MPolygon(poly[0], ++numEle, ePart, own, parent);
+ else p1 = new MPolygon(poly[0], n, ePart, own, parent);
own = false;
int reg = getElementaryTag(-1, elementary, newElemTags[2]);
std::vector<int> phys;
@@ -942,10 +948,11 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
borders[1].insert(std::pair<MElement*, MElement*>(p1->getDomain(i), p1));
}
if(poly[1].size()) {
+ int n = (e->getParent() && poly[0].size() == 0) ? e->getNum() : ++numEle;
if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
- p2 = new MPolygonBorder(poly[1], ++numEle, ePart, own, parent,
+ p2 = new MPolygonBorder(poly[1], n, ePart, own, parent,
copy->getDomain(0), copy->getDomain(1));
- else p2 = new MPolygon(poly[1], ++numEle, ePart, own, parent);
+ else p2 = new MPolygon(poly[1], n, ePart, own, parent);
elements[8][elementary].push_back(p2);
assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
for(int i = 0; i < 2; i++)
@@ -1011,7 +1018,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
MLine *lin;
- if(p1 || p2) lin = new MLineBorder(mv[0], mv[1], ++numEle, ePart, p1, p2);
+ if(p1 || p2){
+ if(!p1) lin = new MLineBorder(mv[0], mv[1], ++numEle, ePart, p2, p1);
+ else lin = new MLineBorder(mv[0], mv[1], ++numEle, ePart, p1, p2);
+ }
else lin = new MLine(mv[0], mv[1], ++numEle, ePart);
int lsL = lines[i]->lsTag();
int c = elements[1].count(lsL);
@@ -1058,8 +1068,9 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
MLine *ml;
- if(eType != MSH_LIN_B) ml = new MLineChild(mv[0], mv[1], ++numEle, ePart, own, parent);
- else ml = new MLineBorder(mv[0], mv[1], ++numEle, ePart,
+ int n = (e->getParent() && i == nbL) ? e->getNum() : ++numEle;
+ if(eType != MSH_LIN_B) ml = new MLineChild(mv[0], mv[1], n, ePart, own, parent);
+ else ml = new MLineBorder(mv[0], mv[1], n, ePart,
copy->getDomain(0), copy->getDomain(1));
own = false;
int reg = getElementaryTag(lines[i]->lsTag(), elementary, newElemTags[1]);
@@ -1160,6 +1171,18 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
}
}
+ int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements(); //element number increment
+ for(unsigned int i = 0; i < gmEntities.size(); i++) {
+ for(unsigned int j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
+ MElement *e = gmEntities[i]->getMeshElement(j);
+ if(e->getNum() > numEle)
+ numEle = e->getNum();
+ if(e->getParent())
+ if(e->getParent()->getNum() > numEle)
+ numEle = e->getParent()->getNum();
+ }
+ }
+
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++){
@@ -1167,7 +1190,6 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
newPhysTags[d][0] = gm->getMaxPhysicalNumber(d); //max value at [dim][0]
}
- int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements(); //element number increment
std::map<MElement*, MElement*> newParents; //map<oldParent, newParent>
std::map<MElement*, MElement*> newDomains; //map<oldDomain, newDomain>