diff --git a/Geo/MElement.h b/Geo/MElement.h
index 6f98c42918f07210126c84662aaea07cfcf49520..46107c3bb2f4b35ffc56877aa770ae19488d2a9e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -216,8 +216,7 @@ class MElement
int order=-1);
// integration routines
- virtual int getNumIntegrationPointsToAllocate(int pOrder) const { return 0; }
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
Msg::Error("No integration points defined for this type of element");
}
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index fdecf9a6705d07ca149e953da4a9ac42591b7f23..7e20a17989dd748f41b7b7d5d90f1697953fd4d1 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -96,10 +96,12 @@ bool MPolyhedron::isInside(double u, double v, double w)
return false;
}
-void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = 0;
double jac[3][3];
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQTetPts(pOrder) * _parts.size()];
for(unsigned int i = 0; i < _parts.size(); i++) {
int nptsi;
IntPt *ptsi;
@@ -124,13 +126,14 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
const double weight = ptsi[ip].weight;
const double detJ = tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
- (*pts[*npts + ip]).pt[0] = p.x();
- (*pts[*npts + ip]).pt[1] = p.y();
- (*pts[*npts + ip]).pt[2] = p.z();
- (*pts[*npts + ip]).weight = detJ * weight;
+ _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 = detJ * weight;
}
*npts += nptsi;
}
+ *pts = _intpt;
}
void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
@@ -275,10 +278,12 @@ bool MPolygon::isInside(double u, double v, double w)
return false;
}
-void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = 0;
double jac[3][3];
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQTPts(pOrder) * _parts.size()];
for(unsigned int i = 0; i < _parts.size(); i++) {
int nptsi;
IntPt *ptsi;
@@ -300,13 +305,14 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
const double weight = ptsi[ip].weight;
const double detJ = tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
- (*pts[*npts + ip]).pt[0] = p.x();
- (*pts[*npts + ip]).pt[1] = p.y();
- (*pts[*npts + ip]).pt[2] = p.z();
- (*pts[*npts + ip]).weight = detJ * weight;
+ _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 = detJ * weight;
}
*npts += nptsi;
}
+ *pts = _intpt;
}
void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
@@ -362,7 +368,7 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
//----------------------------------- MTriangleBorder ------------------------------
-void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
double uvw[3][3];
for(int j = 0; j < 3; j++) {
@@ -391,7 +397,7 @@ void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) c
//-------------------------------------- MLineBorder ------------------------------
-void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
double uvw[2][3];
for(int j = 0; j < 2; j++) {
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 7027e5d9afb947190848298cb8aaf7718c5b6b0d..fa37382c5a5d12df1e5202e392f494840078f6e6 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -19,6 +19,7 @@ class MPolyhedron : public MElement {
protected:
bool _owner;
MElement* _orig;
+ IntPt *_intpt;
std::vector<MTetrahedron*> _parts;
std::vector<MVertex*> _vertices;
std::vector<MVertex*> _innerVertices;
@@ -28,7 +29,7 @@ class MPolyhedron : public MElement {
public:
MPolyhedron(std::vector<MVertex*> v, int num = 0, int part = 0,
bool owner = false, MElement* orig = NULL)
- : MElement(num, part), _owner(owner), _orig(orig)
+ : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
{
if(v.size() % 4){
Msg::Error("Got %d vertices for polyhedron", (int)v.size());
@@ -40,7 +41,7 @@ class MPolyhedron : public MElement {
}
MPolyhedron(std::vector<MTetrahedron*> vT, int num = 0, int part = 0,
bool owner = false, MElement* orig = NULL)
- : MElement(num, part), _owner(owner), _orig(orig)
+ : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
{
for(unsigned int i = 0; i < vT.size(); i++)
_parts.push_back(vT[i]);
@@ -52,6 +53,7 @@ class MPolyhedron : public MElement {
delete _orig;
for(unsigned int i = 0; i < _parts.size(); i++)
delete _parts[i];
+ if(_intpt) delete [] _intpt;
}
virtual int getDim() { return 3; }
virtual int getNumVertices() const { return _vertices.size() + _innerVertices.size(); }
@@ -120,11 +122,7 @@ class MPolyhedron : public MElement {
_orig->xyz2uvw(xyz,uvw);
}
virtual bool isInside(double u, double v, double w);
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
- virtual int getNumIntegrationPointsToAllocate (int pOrder)
- {
- return _parts.size() * getNGQTetPts(pOrder);
- }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
int num=0, int elementary=1, int physical=1, int parentNum=0);
virtual MElement *getParent() const { return _orig; }
@@ -138,6 +136,7 @@ class MPolygon : public MElement {
protected:
bool _owner;
MElement* _orig;
+ IntPt *_intpt;
std::vector<MTriangle*> _parts;
std::vector<MVertex*> _vertices;
std::vector<MVertex*> _innerVertices;
@@ -146,7 +145,7 @@ class MPolygon : public MElement {
public:
MPolygon(std::vector<MVertex*> v, int num = 0, int part = 0,
bool owner = false, MElement* orig = NULL)
- : MElement(num, part), _owner(owner), _orig(orig)
+ : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
{
for(unsigned int i = 0; i < v.size() / 3; i++)
_parts.push_back(new MTriangle(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]));
@@ -154,7 +153,7 @@ class MPolygon : public MElement {
}
MPolygon(std::vector<MTriangle*> vT, int num = 0, int part = 0,
bool owner = false, MElement* orig = NULL)
- : MElement(num, part), _owner(owner), _orig(orig)
+ : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
{
for(unsigned int i = 0; i < vT.size(); i++){
MTriangle *t = (MTriangle*) vT[i];
@@ -168,6 +167,7 @@ class MPolygon : public MElement {
delete _orig;
for(unsigned int i = 0; i < _parts.size(); i++)
delete _parts[i];
+ if(_intpt) delete [] _intpt;
}
virtual int getDim(){ return 2; }
virtual int getNumVertices() const { return _vertices.size() + _innerVertices.size(); }
@@ -224,11 +224,7 @@ class MPolygon : public MElement {
_orig->xyz2uvw(xyz,uvw);
}
virtual bool isInside(double u, double v, double w);
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
- virtual int getNumIntegrationPointsToAllocate (int pOrder)
- {
- return _parts.size() * getNGQTPts(pOrder);
- }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
int num=0, int elementary=1, int physical=1, int parentNum=0);
virtual MElement *getParent() const { return _orig; }
@@ -271,11 +267,7 @@ class MTriangleBorder : public MTriangle {
{
getParent()->xyz2uvw(xyz,uvw);
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
- virtual int getNumIntegrationPointsToAllocate (int pOrder)
- {
- return getNGQTPts(pOrder);
- }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
class MLineBorder : public MLine {
@@ -311,11 +303,7 @@ class MLineBorder : public MLine {
{
getParent()->xyz2uvw(xyz,uvw);
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
- virtual int getNumIntegrationPointsToAllocate (int pOrder)
- {
- return pOrder / 2 + 1;
- }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
GModel *buildCutMesh(GModel *gm, gLevelset *ls,
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index f7bbf706dcd667175a378421c114d5032da391fd..8c53b9438c6ceac03466a2461b54c74f52b9577c 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -21,7 +21,7 @@ int MHexahedron::getVolumeSign()
return sign(det3x3(mat));
}
-void MHexahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MHexahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQHPts(pOrder);
*pts = getGQHPts(pOrder);
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 31a59a7cca40d09a9b8a1533011bdcc12ec4d830..350a23e23bcea6aebbbca4afc10174c8552f5d07 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -171,7 +171,7 @@ class MHexahedron : public MElement {
return false;
return true;
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
private:
int edges_hexa(const int edge, const int vert) const
{
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 16c1a7a6dbc5096b5d76048e72d4c073b497c4da..cdd91eb117e667e31bc0270a87901eb1c59271bc 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -23,7 +23,7 @@ const functionSpace* MLine::getFunctionSpace(int o) const
return 0;
}
-void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
static IntPt GQL[100];
double *t, *w;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 265c4662ebed76f85722c37ef36aca0e8a387bc8..c78fe34652fa691f3a1ab38a0f7209d8f3d3a66e 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -71,7 +71,7 @@ class MLine : public MElement {
return false;
return true;
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
/*
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 6880925b401ec8499927e9251c4d0d0c85f0e431..3daf1324c1f5ad2b7a5b98cf18d0fae633725fed 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -153,7 +153,7 @@ void MQuadrangle9::getFaceRep(int num, double *x, double *y, double *z, SVector3
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
-void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQQPts(pOrder);
*pts = getGQQPts(pOrder);
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 62603daf10fed9431ee9147855e979053b49f9a3..d6457574b0787d6b698561b9d928599ff3970f65 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -127,7 +127,7 @@ class MQuadrangle : public MElement {
return false;
return true;
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
private:
int edges_quad(const int edge, const int vert) const
{
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index fb4d18f7918bb11f16d0ac1391b6dd511e79cf52..e4dcfaadcc29da91f26f45e5227770292f3ca3fb 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -249,7 +249,7 @@ void MTetrahedron10::getFaceRep(int num, double *x, double *y, double *z, SVecto
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
-void MTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQTetPts(pOrder);
*pts = getGQTetPts(pOrder);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 0d0c2745071069a4e935e47a44798cc29eea6085..dec9e6e20a47ae64b3b8118880fc5c8aa2ed4ac3 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -138,7 +138,7 @@ class MTetrahedron : public MElement {
return false;
return true;
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual SPoint3 circumcenter();
private:
int edges_tetra(const int edge, const int vert) const
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 1a8ecf699fc9c7dd1fb9e561f7c13563f3f3b57f..4d66d6a1264f837da131a257261f65bc16d92c68 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -193,7 +193,7 @@ void MTriangle6::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
-void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
*npts = getNGQTPts(pOrder);
*pts = getGQTPts(pOrder);
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index c63c328c04f89a37d797f73b86bac785951a19df..e59c1b159a91510807db523518289e8e8e5d6bb7 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -129,7 +129,7 @@ class MTriangle : public MElement {
return false;
return true;
}
- virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual SPoint3 circumcenter();
private:
int edges_tri(const int edge, const int vert) const