diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index e693d950565e0f672c1b1538a266b825e01ad344..01f8e88eb86a5a198a3ebbddf1fc0e1f6df8a371 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -622,7 +622,7 @@ void GModel::getMeshVerticesForPhysicalGroup(int dim, int num, std::vector<MVert
MElement *GModel::getMeshElementByTag(int num)
{
- int n = (*_newElemNumbers)(num);
+ int n = (_newElemNumbers) ? (*_newElemNumbers)(num) : num;
if(_elementVectorCache.empty() && _elementMapCache.empty()){
Msg::Debug("Rebuilding mesh element cache");
_elementVectorCache.clear();
@@ -637,14 +637,16 @@ MElement *GModel::getMeshElementByTag(int num)
for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
- _elementVectorCache[(*_newElemNumbers)(e->getNum())] = e;
+ int ni = (_newElemNumbers) ? (*_newElemNumbers)(e->getNum()) : e->getNum();
+ _elementVectorCache[ni] = e;
}
}
else{
for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
- _elementMapCache[(*_newElemNumbers)(e->getNum())] = e;
+ int ni = (_newElemNumbers) ? (*_newElemNumbers)(e->getNum()) : e->getNum();
+ _elementMapCache[ni] = e;
}
}
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index e1d8f5d555034c1908b018edf5440dd1848976bb..47946a3669e81177108edeb0158883393d92ddf8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -364,6 +364,29 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3])
}
}
+void MElement::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
+{
+ if(!getParent()) return;
+ SPoint3 p;
+ getParent()->pnt(u, v, w, p);
+ double xyz[3] = {p.x(), p.y(), p.z()};
+ double uvwE[3];
+ xyz2uvw(xyz, uvwE);
+ u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
+}
+
+void MElement::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
+{
+ if(!getParent()) return;
+ SPoint3 p;
+ double uvwE[3] = {u, v, w};
+ pnt(u, v, w, p);
+ double xyz[3] = {p.x(), p.y(), p.z()};
+ double uvwP[3];
+ getParent()->xyz2uvw(xyz, uvwP);
+ u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
+}
+
double MElement::interpolate(double val[], double u, double v, double w, int stride,
int order)
{
@@ -883,12 +906,12 @@ MElement *MElement::copy(int &num, std::map<int, MVertex*> &vertexMap,
}
}
MElementFactory factory;
- MElement *newEl = factory.create(eType, vmv, ++num, _partition);
+ MElement *newEl = factory.create(eType, vmv, getNum(), _partition);
if(eParent && !getDomain(0) && !getDomain(1)) {
std::map<MElement*, MElement*>::iterator it = newParents.find(eParent);
MElement *newParent;
if(it == newParents.end()) {
- newParent = eParent->copy(num, vertexMap, newParents, newDomains);
+ newParent = eParent->copy(++num, vertexMap, newParents, newDomains);
newParents[eParent] = newParent;
}
else
@@ -901,7 +924,7 @@ MElement *MElement::copy(int &num, std::map<int, MVertex*> &vertexMap,
std::map<MElement*, MElement*>::iterator it = newDomains.find(dom);
MElement *newDom;
if(it == newDomains.end()) {
- newDom = dom->copy(num, vertexMap, newParents, newDomains);
+ newDom = dom->copy(++num, vertexMap, newParents, newDomains);
newDomains[dom] = newDom;
}
else
diff --git a/Geo/MElement.h b/Geo/MElement.h
index d13f59b206792236ad22b9535bf44cf95cc79568..ff1165a37cb9199054514bdc5ad5ecba675dfd7f 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -229,6 +229,10 @@ class MElement
// invert the parametrisation
virtual void xyz2uvw(double xyz[3], double uvw[3]);
+ // move point between parent and element parametric spaces
+ virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
+ virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
+
// test if a point, given in parametric coordinates, belongs to the
// element
virtual bool isInside(double u, double v, double w) = 0;
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index eb4267ecf0c6eb7cb4127d40dfb8b83008714893..e38e0ad9f99babaaf73ad09eaab17cbbc69c2d25 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -341,6 +341,12 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
+ *npts = 0;
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQTPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+
double uvw[3][3];
for(int j = 0; j < 3; j++) {
double xyz[3] = {_v[j]->x(), _v[j]->y(), _v[j]->z()};
@@ -350,26 +356,34 @@ void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
MVertex v2(uvw[2][0], uvw[2][1], uvw[2][2]);
MTriangle tt(&v0, &v1, &v2);
- tt.getIntegrationPoints(pOrder, npts, pts);
+ tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
double jac[3][3];
- for(int ip = 0; ip < (*npts); ip++){
- const double u = pts[ip]->pt[0];
- const double v = pts[ip]->pt[1];
- const double w = pts[ip]->pt[2];
- const double weight = pts[ip]->weight;
+ 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];
+ const double weight = ptsi[ip].weight;
const double detJ = tt.getJacobian(u, v, w, jac);
SPoint3 p; tt.pnt(u, v, w, p);
- (*pts[ip]).pt[0] = p.x();
- (*pts[ip]).pt[1] = p.y();
- (*pts[ip]).pt[2] = p.z();
- (*pts[ip]).weight = detJ * weight;
+ _intpt[ip].pt[0] = p.x();
+ _intpt[ip].pt[1] = p.y();
+ _intpt[ip].pt[2] = p.z();
+ _intpt[ip].weight = weight;//detJ * weight;
}
+ *npts = nptsi;
+ *pts = _intpt;
}
//----------------------------------- MPolygonBorder ------------------------------
void MPolygonBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
+ *npts = 0;
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQTPts(pOrder) * _parts.size()];
+ int nptsi;
+ IntPt *ptsi;
+
std::vector<MVertex*> verts;
for(unsigned int i = 0; i < _parts.size(); i++) {
double uvw[3][3];
@@ -383,20 +397,22 @@ void MPolygonBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
verts.push_back(new MVertex(uvw[2][0], uvw[2][1], uvw[2][2]));
}
MPolygon pp(verts);
- pp.getIntegrationPoints(pOrder, npts, pts);
+ pp.getIntegrationPoints(pOrder, &nptsi, &ptsi);
double jac[3][3];
- for(int ip = 0; ip < (*npts); ip++){
- const double u = pts[ip]->pt[0];
- const double v = pts[ip]->pt[1];
- const double w = pts[ip]->pt[2];
- const double weight = pts[ip]->weight;
+ 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];
+ const double weight = ptsi[ip].weight;
const double detJ = pp.getJacobian(u, v, w, jac);
SPoint3 p; pp.pnt(u, v, w, p);
- (*pts[ip]).pt[0] = p.x();
- (*pts[ip]).pt[1] = p.y();
- (*pts[ip]).pt[2] = p.z();
- (*pts[ip]).weight = detJ * weight;
+ _intpt[ip].pt[0] = p.x();
+ _intpt[ip].pt[1] = p.y();
+ _intpt[ip].pt[2] = p.z();
+ _intpt[ip].weight = weight;//detJ * weight;
}
+ *npts = nptsi;
+ *pts = _intpt;
for(unsigned int i = 0; i < verts.size(); i++)
delete verts[i];
}
@@ -405,6 +421,12 @@ void MPolygonBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
+ *npts = 0;
+ if(_intpt) delete [] _intpt;
+ _intpt = new IntPt[getNGQLPts(pOrder)];
+ int nptsi;
+ IntPt *ptsi;
+
double uvw[2][3];
for(int j = 0; j < 2; j++) {
double xyz[3] = {_v[j]->x(), _v[j]->y(), _v[j]->z()};
@@ -413,20 +435,22 @@ void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
MLine ll(&v0, &v1);
- ll.getIntegrationPoints(pOrder, npts, pts);
+ ll.getIntegrationPoints(pOrder, &nptsi, &ptsi);
double jac[3][3];
- for(int ip = 0; ip < (*npts); ip++){
- const double u = pts[ip]->pt[0];
- const double v = pts[ip]->pt[1];
- const double w = pts[ip]->pt[2];
- const double weight = pts[ip]->weight;
+ 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];
+ const double weight = ptsi[ip].weight;
const double detJ = ll.getJacobian(u, v, w, jac);
SPoint3 p; ll.pnt(u, v, w, p);
- (*pts[ip]).pt[0] = p.x();
- (*pts[ip]).pt[1] = p.y();
- (*pts[ip]).pt[2] = p.z();
- (*pts[ip]).weight = detJ * weight;
+ _intpt[ip].pt[0] = p.x();
+ _intpt[ip].pt[1] = p.y();
+ _intpt[ip].pt[2] = p.z();
+ _intpt[ip].weight = weight;//detJ * weight;
}
+ *npts = nptsi;
+ *pts = _intpt;
}
//---------------------------------------- CutMesh ----------------------------
@@ -1120,7 +1144,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
newPhysTags[d][0] = gm->getMaxPhysicalNumber(d); //max value at [dim][0]
}
- int numEle = 0; //element number increment
+ int numEle = gm->getNumMeshElements(); //element number increment
std::map<MElement*, MElement*> newParents; //map<oldParent, newParent>
std::map<MElement*, MElement*> newDomains; //map<oldDomain, newDomain>
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 40c0529c1003bd938b5265268864d7b98c8357c2..359fb2aabfab3d064b5ce997ef6c893e81363f3c 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -320,16 +320,17 @@ class MLineChild : public MLine {
class MTriangleBorder : public MTriangle {
protected:
MElement* _domains[2];
+ IntPt *_intpt;
public:
MTriangleBorder(MVertex *v0, MVertex *v1, MVertex *v2, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MTriangle(v0, v1, v2, num, part)
+ : MTriangle(v0, v1, v2, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
MTriangleBorder(std::vector<MVertex*> v, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MTriangle(v, num, part)
+ : MTriangle(v, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
@@ -342,38 +343,24 @@ class MTriangleBorder : public MTriangle {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_TRI_B; }
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
- {
- return getParent()->getFunctionSpace(order);
- }
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
- {
- getParent()->getShapeFunctions(u, v, w, s, o);
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
- {
- getParent()->getGradShapeFunctions(u, v, w, s, o);
- }
- virtual void xyz2uvw(double xyz[3], double uvw[3])
- {
- getParent()->xyz2uvw(xyz,uvw);
- }
+ // the integration points of the MTriangleBorder are in the parent element space
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
class MPolygonBorder : public MPolygon {
protected:
MElement* _domains[2];
+ IntPt *_intpt;
public:
MPolygonBorder(std::vector<MTriangle*> v, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MPolygon(v, num, part)
+ : MPolygon(v, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
MPolygonBorder(std::vector<MVertex*> v, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MPolygon(v, num, part)
+ : MPolygon(v, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
@@ -386,38 +373,24 @@ class MPolygonBorder : public MPolygon {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_POLYG_B; }
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
- {
- return getParent()->getFunctionSpace(order);
- }
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
- {
- getParent()->getShapeFunctions(u, v, w, s, o);
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
- {
- getParent()->getGradShapeFunctions(u, v, w, s, o);
- }
- virtual void xyz2uvw(double xyz[3], double uvw[3])
- {
- getParent()->xyz2uvw(xyz,uvw);
- }
+ // the integration points of the MPolygonBorder are in the parent element space
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};
class MLineBorder : public MLine {
protected:
MElement* _domains[2];
+ IntPt *_intpt;
public:
MLineBorder(MVertex *v0, MVertex *v1, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MLine(v0, v1, num, part)
+ : MLine(v0, v1, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
MLineBorder(std::vector<MVertex*> v, int num = 0, int part = 0,
MElement* d1 = NULL, MElement* d2 = NULL)
- : MLine(v, num, part)
+ : MLine(v, num, part), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
@@ -430,23 +403,7 @@ class MLineBorder : public MLine {
return NULL;
}
virtual int getTypeForMSH() const { return MSH_LIN_B; }
- virtual const polynomialBasis* getFunctionSpace(int order=-1) const
- {
- if (this->getParent()) return getParent()->getFunctionSpace(order);
- else return NULL;
- }
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
- {
- getParent()->getShapeFunctions(u, v, w, s, o);
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
- {
- getParent()->getGradShapeFunctions(u, v, w, s, o);
- }
- virtual void xyz2uvw(double xyz[3], double uvw[3])
- {
- getParent()->xyz2uvw(xyz,uvw);
- }
+ // the integration points of the MLineBorder are in the parent element space
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
};