diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 874b15a4d25e809c33454429269a9e7dd0ca1091..8770d1b851ed929e06e3cddd3ca0907946cfc867 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -198,7 +198,12 @@
// Additional types
#define MSH_PYR_1 132
-#define MSH_NUM_TYPE 132
+#define MSH_PNT_SUB 133
+#define MSH_LIN_SUB 134
+#define MSH_TRI_SUB 135
+#define MSH_TET_SUB 136
+
+#define MSH_NUM_TYPE 136
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 17466e538b478fe9b6692012ab3ae94ce84100b8..275e064e0d6e4b528339ee6e99486789eba426d1 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -39,7 +39,7 @@ set(SRC
MFace.cpp
MElement.cpp MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
- MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
+ MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp MSubElement.cpp
MZone.cpp MZoneBoundary.cpp
Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp Chain.cpp
Curvature.cpp
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index a66072fa4cd23fc9a4b86c789994b4a2dd91f6ea..3ba8a89fdc2b6bc84fb24af6c65491ea86db368b 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -816,41 +816,61 @@ static int getNumElementsMSH(GModel *m, bool saveAll, int saveSinglePartition)
int n = 0;
for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it)
+ {
n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
- for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
+ if(!CTX::instance()->mesh.saveTri){
+ for(unsigned int i = 0; i < (*it)->points.size(); i++)
+ if((*it)->points[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
+ }
+ }
+ for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
+ {
n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
if(!CTX::instance()->mesh.saveTri){
for(unsigned int i = 0; i < (*it)->lines.size(); i++)
- if((*it)->lines[i]->ownsParent())
- n += (saveAll ? 1 : (*it)->physicals.size());
+ if((*it)->lines[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
}
}
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ {
n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
if(CTX::instance()->mesh.saveTri){
- for(unsigned int i = 0; i < (*it)->polygons.size(); i++){
+ for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
+ {
int nbC = (*it)->polygons[i]->getNumChildren()-1;
- n += (saveAll ? nbC : nbC * (*it)->physicals.size());
+ n += (saveAll ? nbC : nbC * (*it)->physicals.size());
}
}
- else{
+ else
+ {
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
+ if((*it)->triangles[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
- if((*it)->polygons[i]->ownsParent())
- n += (saveAll ? 1 : (*it)->physicals.size());
+ if((*it)->polygons[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
}
}
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
+ {
n += getNumElementsMSH(*it, saveAll, saveSinglePartition);
if(CTX::instance()->mesh.saveTri){
- for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++){
+ for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++)
+ {
int nbC = (*it)->polyhedra[i]->getNumChildren()-1;
- n += (saveAll ? nbC : nbC * (*it)->physicals.size());
+ n += (saveAll ? nbC : nbC * (*it)->physicals.size());
}
}
- else{
+ else
+ {
+ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
+ if((*it)->tetrahedra[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++)
- if((*it)->polyhedra[i]->ownsParent())
- n += (saveAll ? 1 : (*it)->physicals.size());
+ if((*it)->polyhedra[i]->ownsParent())
+ n += (saveAll ? 1 : (*it)->physicals.size());
}
}
return n;
@@ -939,11 +959,26 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
//parents
if (!CTX::instance()->mesh.saveTri){
+ for(viter it = firstVertex(); it != lastVertex(); ++it)
+ for(unsigned int i = 0; i < (*it)->points.size(); i++)
+ if((*it)->points[i]->ownsParent())
+ writeElementMSH(fp, this, (*it)->points[i]->getParent(),
+ saveAll, version, binary, num, (*it)->tag(), (*it)->physicals);
for(eiter it = firstEdge(); it != lastEdge(); ++it)
for(unsigned int i = 0; i < (*it)->lines.size(); i++)
if((*it)->lines[i]->ownsParent())
writeElementMSH(fp, this, (*it)->lines[i]->getParent(),
saveAll, version, binary, num, (*it)->tag(), (*it)->physicals);
+ for(fiter it = firstFace(); it != lastFace(); ++it)
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
+ if((*it)->triangles[i]->ownsParent())
+ writeElementMSH(fp, this, (*it)->triangles[i]->getParent(),
+ saveAll, version, binary, num, (*it)->tag(), (*it)->physicals);
+ for(riter it = firstRegion(); it != lastRegion(); ++it)
+ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++)
+ if((*it)->tetrahedra[i]->ownsParent())
+ writeElementMSH(fp, this, (*it)->tetrahedra[i]->getParent(),
+ saveAll, version, binary, num, (*it)->tag(), (*it)->physicals);
for(fiter it = firstFace(); it != lastFace(); ++it)
for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
if((*it)->polygons[i]->ownsParent())
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index ca89017d999ce28202f23ecb2b1998a87e3f3dc7..6190fd8c989fbbc8fc1fc5dc8a753009807ba3ec 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -17,6 +17,7 @@
#include "MPrism.h"
#include "MPyramid.h"
#include "MElementCut.h"
+#include "MSubElement.h"
#include "GEntity.h"
#include "StringUtils.h"
#include "Numeric.h"
@@ -734,18 +735,30 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
}
return;
}
- if(type == MSH_LIN_B || type == MSH_LIN_C){
- MLine *l = new MLine(getVertex(0), getVertex(1));
- l->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
- delete l;
+ if(type == MSH_TET_SUB){
+ MTetrahedron *tt = new MTetrahedron(getVertex(0), getVertex(1), getVertex(2), getVertex(3));
+ tt->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
+ delete tt;
return;
}
- if(type == MSH_TRI_B){
+ if(type == MSH_TRI_B || type == MSH_TRI_SUB){
MTriangle *t = new MTriangle(getVertex(0), getVertex(1), getVertex(2));
t->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
delete t;
return;
}
+ if(type == MSH_LIN_B || type == MSH_LIN_C || type == MSH_LIN_SUB){
+ MLine *l = new MLine(getVertex(0), getVertex(1));
+ l->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
+ delete l;
+ return;
+ }
+ if(type == MSH_PNT_SUB){
+ MPoint *p = new MPoint(getVertex(0));
+ p->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
+ delete p;
+ return;
+ }
}
if(!binary){
@@ -1171,6 +1184,10 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8;
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
+ case MSH_PNT_SUB: if(name) *name = "Point Xfem"; return 1;
+ case MSH_LIN_SUB: if(name) *name = "Line Xfem"; return 2;
+ case MSH_TRI_SUB: if(name) *name = "Triangle Xfem"; return 3;
+ case MSH_TET_SUB: if(name) *name = "Tetrahedron Xfem";return 4;
default:
Msg::Error("Unknown type of element %d", typeMSH);
if(name) *name = "Unknown";
@@ -1332,6 +1349,10 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_HEX_512: return new MHexahedronN(v, 7, num, part);
case MSH_HEX_729: return new MHexahedronN(v, 8, num, part);
case MSH_HEX_1000:return new MHexahedronN(v, 9, num, part);
+ case MSH_PNT_SUB:return new MSubPoint(v, num, part, owner, parent);
+ case MSH_LIN_SUB:return new MSubLine(v, num, part, owner, parent);
+ case MSH_TRI_SUB:return new MSubTriangle(v, num, part, owner, parent);
+ case MSH_TET_SUB:return new MSubTetrahedron(v, num, part, owner, parent);
default: return 0;
}
}
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b2f83347f8f911e32658e22e7b542d872c119cc8
--- /dev/null
+++ b/Geo/MSubElement.cpp
@@ -0,0 +1,341 @@
+#include "MSubElement.h"
+
+
+
+//--------------------------------------------------------------------------
+// MSubTetrahedron
+//--------------------------------------------------------------------------
+
+MSubTetrahedron::~MSubTetrahedron()
+{
+ if(_owner)
+ delete _orig;
+}
+
+const polynomialBasis* MSubTetrahedron::getFunctionSpace(int order) const
+{
+ if(_orig) return _orig->getFunctionSpace(order);
+ return 0;
+}
+const JacobianBasis* MSubTetrahedron::getJacobianFuncSpace(int order) const
+{
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+}
+void MSubTetrahedron::getShapeFunctions(double u, double v, double w, double s[], int o)
+{
+ if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+}
+void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+{
+ if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+}
+void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+{
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+}
+
+bool MSubTetrahedron::isInside(double u, double v, double w)
+{
+ if(!_orig) return false;
+ double uvw[3] = {u, v, w};
+ double v_uvw[4][3];
+ for(int i=0; i<4; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MVertex v3(v_uvw[3][0], v_uvw[3][1], v_uvw[3][2]);
+ MTetrahedron t(&v0, &v1, &v2, &v3);
+ double ksi[3];
+ t.xyz2uvw(uvw, ksi);
+ if(t.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
+}
+void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ *npts=0;
+ if(_intpt) delete [] _intpt;
+ if(!_orig) return;
+ _intpt = new IntPt[getNGQTetPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+
+ double v_uvw[4][3];
+ for(int i=0; i<4; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MVertex v3(v_uvw[3][0], v_uvw[3][1], v_uvw[3][2]);
+ MTetrahedron tt(&v0, &v1, &v2, &v3);
+ tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
+ double jac[3][3];
+ for(int ip=0; ip<nptsi; ++ip)
+ {
+ const double u = ptsi[ip].pt[0];
+ const double v = ptsi[ip].pt[1];
+ const double w = ptsi[ip].pt[2];
+ tt.getJacobian(u, v, w, jac);
+ SPoint3 p; tt.pnt(u, v, w, p);
+ _intpt[*npts + ip].pt[0] = p.x();
+ _intpt[*npts + ip].pt[1] = p.y();
+ _intpt[*npts + ip].pt[2] = p.z();
+ _intpt[*npts + ip].weight = ptsi[ip].weight;
+ }
+ *npts = nptsi;
+ *pts = _intpt;
+}
+
+
+//--------------------------------------------------------------------------
+// MSubTriangle
+//--------------------------------------------------------------------------
+
+MSubTriangle::~MSubTriangle()
+{
+ if(_owner)
+ delete _orig;
+}
+
+const polynomialBasis* MSubTriangle::getFunctionSpace(int order) const
+{
+ if(_orig) return _orig->getFunctionSpace(order);
+ return 0;
+}
+const JacobianBasis* MSubTriangle::getJacobianFuncSpace(int order) const
+{
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+}
+void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[], int o)
+{
+ if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+}
+void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+{
+ if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+}
+void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+{
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+}
+
+bool MSubTriangle::isInside(double u, double v, double w)
+{
+ if(!_orig) return false;
+ double uvw[3] = {u, v, w};
+ double v_uvw[3][3];
+ for(int i=0; i<3; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MTriangle t(&v0, &v1, &v2);
+ double ksi[3];
+ t.xyz2uvw(uvw, ksi);
+ if(t.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
+}
+void MSubTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ *npts=0;
+ if(_intpt) delete [] _intpt;
+ if(!_orig) return;
+ _intpt = new IntPt[getNGQTPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+
+ double v_uvw[3][3];
+ for(int i=0; i<3; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MVertex v2(v_uvw[2][0], v_uvw[2][1], v_uvw[2][2]);
+ MTriangle t(&v0, &v1, &v2);
+ t.getIntegrationPoints(pOrder, &nptsi, &ptsi);
+ double jac[3][3];
+ for(int ip=0; ip<nptsi; ++ip)
+ {
+ const double u = ptsi[ip].pt[0];
+ const double v = ptsi[ip].pt[1];
+ const double w = ptsi[ip].pt[2];
+ t.getJacobian(u, v, w, jac);
+ SPoint3 p; t.pnt(u, v, w, p);
+ _intpt[*npts + ip].pt[0] = p.x();
+ _intpt[*npts + ip].pt[1] = p.y();
+ _intpt[*npts + ip].pt[2] = p.z();
+ _intpt[*npts + ip].weight = ptsi[ip].weight;
+ }
+ *npts = nptsi;
+ *pts = _intpt;
+}
+
+
+//--------------------------------------------------------------------------
+// MSubLine
+//--------------------------------------------------------------------------
+
+MSubLine::~MSubLine()
+{
+ if(_owner)
+ delete _orig;
+}
+
+const polynomialBasis* MSubLine::getFunctionSpace(int order) const
+{
+ if(_orig) return _orig->getFunctionSpace(order);
+ return 0;
+}
+const JacobianBasis* MSubLine::getJacobianFuncSpace(int order) const
+{
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+}
+void MSubLine::getShapeFunctions(double u, double v, double w, double s[], int o)
+{
+ if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+}
+void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+{
+ if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+}
+void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+{
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+}
+
+bool MSubLine::isInside(double u, double v, double w)
+{
+ if(!_orig) return false;
+ double uvw[3] = {u, v, w};
+ double v_uvw[2][3];
+ for(int i=0; i<2; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MLine l(&v0, &v1);
+ double ksi[3];
+ l.xyz2uvw(uvw, ksi);
+ if(l.isInside(ksi[0], ksi[1], ksi[2]))
+ return true;
+ return false;
+}
+void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ *npts=0;
+ if(_intpt) delete [] _intpt;
+ if(!_orig) return;
+ _intpt = new IntPt[getNGQLPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+ double v_uvw[2][3];
+ for(int i=0; i<2; ++i)
+ {
+ MVertex *vi = getVertex(i);
+ double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
+ _orig->xyz2uvw(v_xyz, v_uvw[i]);
+ }
+ MVertex v0(v_uvw[0][0], v_uvw[0][1], v_uvw[0][2]);
+ MVertex v1(v_uvw[1][0], v_uvw[1][1], v_uvw[1][2]);
+ MLine l(&v0, &v1);
+ l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
+ for(int ip=0; ip<nptsi; ++ip)
+ {
+ const double u = ptsi[ip].pt[0];
+ const double v = ptsi[ip].pt[1];
+ const double w = ptsi[ip].pt[2];
+ SPoint3 p; l.pnt(u, v, w, p);
+ _intpt[*npts + ip].pt[0] = p.x();
+ _intpt[*npts + ip].pt[1] = p.y();
+ _intpt[*npts + ip].pt[2] = p.z();
+ _intpt[*npts + ip].weight = ptsi[ip].weight;
+ }
+ *npts = nptsi;
+ *pts = _intpt;
+}
+
+//--------------------------------------------------------------------------
+// MSubPoint
+//--------------------------------------------------------------------------
+
+MSubPoint::~MSubPoint()
+{
+ if(_owner)
+ delete _orig;
+}
+
+const polynomialBasis* MSubPoint::getFunctionSpace(int order) const
+{
+ if(_orig) return _orig->getFunctionSpace(order);
+ return 0;
+}
+const JacobianBasis* MSubPoint::getJacobianFuncSpace(int order) const
+{
+ if(_orig) return _orig->getJacobianFuncSpace(order);
+ return 0;
+}
+void MSubPoint::getShapeFunctions(double u, double v, double w, double s[], int o)
+{
+ if(_orig) _orig->getShapeFunctions(u, v, w, s, o);
+}
+void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+{
+ if(_orig) _orig->getGradShapeFunctions(u, v, w, s, o);
+}
+void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o)
+{
+ if(_orig) _orig->getHessShapeFunctions(u, v, w, s, o);
+}
+
+bool MSubPoint::isInside(double u, double v, double w)
+{
+ if(!_orig) return false;
+ MVertex *v0 = getVertex(0);
+ double v_xyz[3] = {v0->x(), v0->y(), v0->z()};
+ double v_uvw[3];
+ _orig->xyz2uvw(v_xyz, v_uvw);
+
+ double p_xyz[3] = {u, v, w};
+ double p_uvw[3];
+ _orig->xyz2uvw(p_xyz, p_uvw);
+
+ double d_xyz[3] = {p_uvw[0]-v_uvw[0], p_uvw[1]-v_uvw[1], p_uvw[2]-v_uvw[2]};
+ double tol = _isInsideTolerance;
+
+ if (d_xyz[0]*d_xyz[0]+d_xyz[1]*d_xyz[1]+d_xyz[2]*d_xyz[2]<tol*tol)
+ return true;
+ return false;
+}
+void MSubPoint::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+{
+ static IntPt GQL[1];
+ GQL[0].pt[0] = 0;
+ GQL[0].pt[1] = 0;
+ GQL[0].pt[2] = 0;
+ GQL[0].weight = 1;
+ *npts = 1;
+ *pts = GQL;
+}
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
new file mode 100644
index 0000000000000000000000000000000000000000..6bf5646e666eb5a2a1778d695cbddae648a90053
--- /dev/null
+++ b/Geo/MSubElement.h
@@ -0,0 +1,144 @@
+//
+// C++ Interface: MSubElement
+//
+// Description:
+//
+//
+// Authors: <Frederic Duboeuf>, (C) 2012
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _MSUBELEMENT_H_
+#define _MSUBELEMENT_H_
+
+#include "GmshMessage.h"
+#include "MElement.h"
+#include "MTetrahedron.h"
+#include "MTriangle.h"
+#include "MLine.h"
+#include "MPoint.h"
+
+
+class MSubTetrahedron : public MTetrahedron
+{
+ protected:
+ bool _owner;
+ MElement* _orig;
+ std::vector<MElement*> _parents;
+ IntPt *_intpt;
+ public:
+ MSubTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MTetrahedron(v0, v1, v2, v3, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ MSubTetrahedron(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MTetrahedron(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ ~MSubTetrahedron();
+ virtual int getTypeForMSH() const { return MSH_TET_SUB; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+ // the parametric coordinates of the LineChildren are
+ // the coordinates in the local parent element.
+ virtual bool isInside(double u, double v, double w);
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual bool ownsParent() const { return _owner; }
+ virtual std::vector<MElement*> getMultiParents() const { return _parents; }
+ virtual void setMultiParent(std::vector<MElement*> &parents, bool owner = false) { _parents = parents; _orig = _parents[0]; _owner = owner; }
+};
+
+class MSubTriangle : public MTriangle
+{
+ protected:
+ bool _owner;
+ MElement* _orig;
+ std::vector<MElement*> _parents;
+ IntPt *_intpt;
+ public:
+ MSubTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ MSubTriangle(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MTriangle(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ ~MSubTriangle();
+ virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+ // the parametric coordinates of the LineChildren are
+ // the coordinates in the local parent element.
+ virtual bool isInside(double u, double v, double w);
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual bool ownsParent() const { return _owner; }
+ virtual std::vector<MElement*> getMultiParents() const { return _parents; }
+ virtual void setMultiParent(std::vector<MElement*> &parents, bool owner = false) { _parents = parents; _orig = _parents[0]; _owner = owner; }
+};
+
+class MSubLine : public MLine
+{
+ protected:
+ bool _owner;
+ MElement* _orig;
+ std::vector<MElement*> _parents;
+ IntPt *_intpt;
+ public:
+ MSubLine(MVertex *v0, MVertex *v1, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ MSubLine(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MLine(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ ~MSubLine();
+ virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+ // the parametric coordinates of the LineChildren are
+ // the coordinates in the local parent element.
+ virtual bool isInside(double u, double v, double w);
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual bool ownsParent() const { return _owner; }
+ virtual std::vector<MElement*> getMultiParents() const { return _parents; }
+ virtual void setMultiParent(std::vector<MElement*> &parents, bool owner = false) { _parents = parents; _orig = _parents[0]; _owner = owner; }
+};
+
+class MSubPoint : public MPoint
+{
+ protected:
+ bool _owner;
+ MElement* _orig;
+ std::vector<MElement*> _parents;
+ IntPt *_intpt;
+ public:
+ MSubPoint(MVertex *v0, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MPoint(v0, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ MSubPoint(std::vector<MVertex*> v, int num=0, int part=0, bool owner=false, MElement* orig=NULL)
+ : MPoint(v, num, part), _owner(owner), _orig(orig), _intpt(0) {}
+ ~MSubPoint();
+ virtual int getTypeForMSH() const { return MSH_PNT_SUB; }
+ virtual const polynomialBasis* getFunctionSpace(int order=-1) const;
+ virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o);
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o);
+ virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int o);
+ // the parametric coordinates of the PointChildren are
+ // the coordinates in the local parent element.
+ virtual bool isInside(double u, double v, double w);
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual MElement *getParent() const { return _orig; }
+ virtual void setParent(MElement *p, bool owner = false) { _orig = p; _owner = owner; }
+ virtual bool ownsParent() const { return _owner; }
+ virtual std::vector<MElement*> getMultiParents() const { return _parents; }
+ virtual void setMultiParent(std::vector<MElement*> &parents, bool owner = false) { _parents = parents; _orig = _parents[0]; _owner = owner; }
+};
+
+#endif // _MSUBELEMENT_H_
\ No newline at end of file