Skip to content
Snippets Groups Projects
Commit 0fa525b5 authored by Gaetan Bricteux's avatar Gaetan Bricteux
Browse files

keep parent number after cut + use shape functions of child

parent 040cff7c
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GModel.cpp
+ 5
− 3
View file @ 0fa525b5
......@@ -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;
}
}
}
......
Geo/MElement.cpp
+ 26
− 3
View file @ 0fa525b5
......@@ -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
......
Geo/MElement.h
+ 4
− 0
View file @ 0fa525b5
......@@ -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;
......
Geo/MElementCut.cpp
+ 55
− 31
View file @ 0fa525b5
......@@ -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>
......
Geo/MElementCut.h
+ 12
− 55
View file @ 0fa525b5
......@@ -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);
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment