Skip to content
Snippets Groups Projects
Commit 48b43a53 authored by Éric Béchet's avatar Éric Béchet
Browse files

No commit message

No commit message
parent 36eae5e3
Branches
Tags
No related merge requests found
......@@ -161,6 +161,10 @@ class MElement
virtual int getNumChildren() const { return 0; }
virtual MElement *getChild(int i) const { return NULL; }
virtual bool ownsParent() const { return false; }
// get base element in case of MSubElement
virtual const MElement *getBaseElement() const { return this; }
virtual MElement *getBaseElement() { return this; }
// get and set domain for borders
virtual MElement *getDomain(int i) const { return NULL; }
virtual void setDomain (MElement *e, int i) { }
......@@ -239,10 +243,10 @@ class MElement
int order=-1);
// return the Jacobian of the element evaluated at point (u,v,w) in
// parametric coordinates
double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
// To be compatible with _vgrads of functionSpace without having to put under fullMatrix form
double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
double getJacobian(double u, double v, double w, double jac[3][3]);
virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
virtual double getJacobian(double u, double v, double w, double jac[3][3]);
inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
double JAC[3][3];
const double detJ = getJacobian (u,v,w,JAC);
......@@ -253,7 +257,7 @@ class MElement
}
return detJ;
}
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
double getJacobianDeterminant(double u, double v, double w)
{
double jac[3][3]; return getJacobian(u, v, w, jac);
......
......@@ -7,12 +7,107 @@
// Frederic Duboeuf
#include "MSubElement.h"
#include "Numeric.h"
// MSubTetrahedron
MSubTetrahedron::~MSubTetrahedron()
{
if(_intpt) delete [] _intpt;
if(_pts) delete [] _pts;
if(_base) delete _base;
}
const nodalBasis* 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 order)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
}
void MSubTetrahedron::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
{
if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
}
void MSubTetrahedron::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
{
if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
}
void MSubTetrahedron::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
{
if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
}
double MSubTetrahedron::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubTetrahedron::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubTetrahedron::getJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getJacobian(u, v, w, jac);
return 0;
}
double MSubTetrahedron::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
return 0;
}
int MSubTetrahedron::getNumShapeFunctions()
{
if(_orig) return _orig->getNumShapeFunctions();
return 0;
}
int MSubTetrahedron::getNumPrimaryShapeFunctions()
{
if(_orig) return _orig->getNumPrimaryShapeFunctions();
return 0;
}
MVertex* MSubTetrahedron::getShapeFunctionNode(int i)
{
if(_orig) return _orig->getShapeFunctionNode(i);
return 0;
}
void MSubTetrahedron::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
_orig->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwE[3];
getBaseElement()->xyz2uvw(xyz, uvwE);
u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
}
void MSubTetrahedron::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
getBaseElement()->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwP[3];
_orig->xyz2uvw(xyz, uvwP);
u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
}
bool MSubTetrahedron::isInside(double u, double v, double w)
......@@ -40,44 +135,57 @@ bool MSubTetrahedron::isInside(double u, double v, double w)
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;
if(_pts)
{
if(pOrder==_pOrder)
{
*npts = _npts;
*pts = _pts;
return;
}
else
delete [] _pts;
}
// work in the parametric space of the parent element
// (i) get the parametric coordinates of MSubElement vertices in this space
double v_uvw[4][3];
for(int i=0; i<4; ++i)
_pOrder = pOrder;
if(!_orig)
{
MVertex *vi = getVertex(i);
double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
_orig->xyz2uvw(v_xyz, v_uvw[i]);
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
*npts = _npts;
*pts = _pts;
return;
}
// (ii) build a MElement on the vertices with these coordinates in this space
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);
// (iii) get the integration points of the new element in its parametric space
tt.getIntegrationPoints(pOrder, &nptsi, &ptsi);
// (iv) get the coordinates of these points in the parametric space of parent element
for(int ip=0; ip<nptsi; ++ip)
// work in the parametric space of the parent element
*npts=0;
_pts = new IntPt[getNGQTetPts(pOrder)];
// (i) get the integration points of the base element in its parametric space
IntPt *ptsb;
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
// (ii) get the coordinates of these points in the parametric space of parent element
double u,v,w;
double jac[3][3];
double baseJac, origJac;
for(int i=0; i<_npts; ++i)
{
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
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;
u = ptsb[i].pt[0];
v = ptsb[i].pt[1];
w = ptsb[i].pt[2];
baseJac = getBaseElement()->getJacobian(u, v, w, jac);
movePointFromElementSpaceToParentSpace(u, v, w);
origJac = _orig->getJacobian(u, v, w, jac);
_pts[*npts + i].pt[0] = u;
_pts[*npts + i].pt[1] = v;
_pts[*npts + i].pt[2] = w;
_pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
}
*npts = nptsi;
*pts = _intpt;
*npts = _npts;
*pts = _pts;
}
......@@ -85,7 +193,148 @@ void MSubTetrahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
MSubTriangle::~MSubTriangle()
{
if(_intpt) delete [] _intpt;
if(_pts) delete [] _pts;
if(_base) delete _base;
}
const nodalBasis* 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 order)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
}
void MSubTriangle::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
{
if(!_orig)
return;
if (_orig->getDim()==getDim())
return _orig->getGradShapeFunctions(u, v, w, s, order);
int nsf = getNumShapeFunctions();
double gradsuvw[nsf][3];
_orig->getGradShapeFunctions(u, v, w, gradsuvw, order);
// work in the parametric space of the parent element
double jac[3][3];
double invjac[3][3];
_orig->getJacobian(u, v, w, jac);
inv3x3(jac, invjac);
MEdge edge[2];
edge[0]=getBaseElement()->getEdge(0);
edge[1]=getBaseElement()->getEdge(1);
SVector3 tang[2];
tang[0]=edge[0].tangent();
tang[1]=edge[1].tangent();
SVector3 vect = crossprod(tang[0],tang[1]);
tang[1] = crossprod(vect,tang[0]);
double gradxyz[3];
double projgradxyz[3];
for (int i=0; i<nsf; ++i)
{
// (i) get the cartesian coordinates of the gradient
gradxyz[0] = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
gradxyz[1] = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
gradxyz[2] = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
// (ii) projection of the gradient on edges in the cartesian space
SVector3 grad(&gradxyz[0]);
double prodscal[2];
prodscal[0]=dot(tang[0],grad);
prodscal[1]=dot(tang[1],grad);
projgradxyz[0] = prodscal[0]*tang[0].x() + prodscal[1]*tang[1].x();
projgradxyz[1] = prodscal[0]*tang[0].y() + prodscal[1]*tang[1].y();
projgradxyz[2] = prodscal[0]*tang[0].z() + prodscal[1]*tang[1].z();
// (iii) get the parametric coordinates of the projection in the parametric space of the parent element
s[i][0] = jac[0][0] * projgradxyz[0] + jac[0][1] * projgradxyz[1] + jac[0][2] * projgradxyz[2];
s[i][1] = jac[1][0] * projgradxyz[0] + jac[1][1] * projgradxyz[1] + jac[1][2] * projgradxyz[2];
s[i][2] = jac[2][0] * projgradxyz[0] + jac[2][1] * projgradxyz[1] + jac[2][2] * projgradxyz[2];
}
}
void MSubTriangle::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
{
if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
}
void MSubTriangle::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
{
if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
}
double MSubTriangle::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubTriangle::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubTriangle::getJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getJacobian(u, v, w, jac);
return 0;
}
double MSubTriangle::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
return 0;
}
int MSubTriangle::getNumShapeFunctions()
{
if(_orig) return _orig->getNumShapeFunctions();
return 0;
}
int MSubTriangle::getNumPrimaryShapeFunctions()
{
if(_orig) return _orig->getNumPrimaryShapeFunctions();
return 0;
}
MVertex* MSubTriangle::getShapeFunctionNode(int i)
{
if(_orig) return _orig->getShapeFunctionNode(i);
return 0;
}
void MSubTriangle::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
_orig->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwE[3];
getBaseElement()->xyz2uvw(xyz, uvwE);
u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
}
void MSubTriangle::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
getBaseElement()->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwP[3];
_orig->xyz2uvw(xyz, uvwP);
u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
}
bool MSubTriangle::isInside(double u, double v, double w)
......@@ -112,51 +361,195 @@ bool MSubTriangle::isInside(double u, double v, double w)
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;
if(_pts)
{
if(pOrder==_pOrder)
{
*npts = _npts;
*pts = _pts;
return;
}
else
delete [] _pts;
}
// work in the parametric space of the parent element
// (i) get the parametric coordinates of MSubElement vertices in this space
double v_uvw[3][3];
for(int i=0; i<3; ++i)
_pOrder = pOrder;
if(!_orig)
{
MVertex *vi = getVertex(i);
double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
_orig->xyz2uvw(v_xyz, v_uvw[i]);
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
*npts = _npts;
*pts = _pts;
return;
}
// (ii) build a MElement on the vertices with these coordinates in this space
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);
// (iii) get the integration points of the new element in its parametric space
t.getIntegrationPoints(pOrder, &nptsi, &ptsi);
// (iv) get the coordinates of these points in the parametric space of parent element
for(int ip=0; ip<nptsi; ++ip)
// work in the parametric space of the parent element
*npts=0;
_pts = new IntPt[getNGQTetPts(pOrder)];
// (i) get the integration points of the base element in its parametric space
IntPt *ptsb;
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
// (ii) get the coordinates of these points in the parametric space of parent element
double u,v,w;
double jac[3][3];
double baseJac, origJac;
for(int i=0; i<_npts; ++i)
{
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
const double w = ptsi[ip].pt[2];
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;
u = ptsb[i].pt[0];
v = ptsb[i].pt[1];
w = ptsb[i].pt[2];
baseJac = getBaseElement()->getJacobian(u, v, w, jac);
movePointFromElementSpaceToParentSpace(u, v, w);
origJac = _orig->getJacobian(u, v, w, jac);
_pts[*npts + i].pt[0] = u;
_pts[*npts + i].pt[1] = v;
_pts[*npts + i].pt[2] = w;
_pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
}
*npts = nptsi;
*pts = _intpt;
*npts = _npts;
*pts = _pts;
}
// MSubLine
MSubLine::~MSubLine()
{
if(_intpt) delete [] _intpt;
if(_pts) delete [] _pts;
if(_base) delete _base;
}
const nodalBasis* 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 order)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
}
void MSubLine::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
{
if(!_orig)
return;
if (_orig->getDim()==getDim())
return _orig->getGradShapeFunctions(u, v, w, s, order);
int nsf = _orig->getNumShapeFunctions();
double gradsuvw[nsf][3];
_orig->getGradShapeFunctions(u, v, w, gradsuvw, order);
double jac[3][3];
double invjac[3][3];
_orig->getJacobian(u, v, w, jac);
inv3x3(jac, invjac);
MEdge edge = getBaseElement()->getEdge(0);
SVector3 tang = edge.tangent();
double gradxyz[3];
double projgradxyz[3];
for (int i=0; i<nsf; ++i)
{
// (i) get the cartesian coordinates of the gradient
gradxyz[0] = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
gradxyz[1] = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
gradxyz[2] = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
// (ii) projection of the gradient on edges in the cartesian space
SVector3 grad(&gradxyz[0]);
double prodscal = dot(tang,grad);
projgradxyz[0] = prodscal * tang.x();
projgradxyz[1] = prodscal * tang.y();
projgradxyz[2] = prodscal * tang.z();
// (iii) get the parametric coordinates of the projection in the parametric space of the parent element
s[i][0] = jac[0][0] * projgradxyz[0] + jac[0][1] * projgradxyz[1] + jac[0][2] * projgradxyz[2];
s[i][1] = jac[1][0] * projgradxyz[0] + jac[1][1] * projgradxyz[1] + jac[1][2] * projgradxyz[2];
s[i][2] = jac[2][0] * projgradxyz[0] + jac[2][1] * projgradxyz[1] + jac[2][2] * projgradxyz[2];
}
}
void MSubLine::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
{
if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
}
void MSubLine::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
{
if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
}
double MSubLine::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubLine::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubLine::getJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getJacobian(u, v, w, jac);
return 0;
}
double MSubLine::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
return 0;
}
int MSubLine::getNumShapeFunctions()
{
if(_orig) return _orig->getNumShapeFunctions();
return 0;
}
int MSubLine::getNumPrimaryShapeFunctions()
{
if(_orig) return _orig->getNumPrimaryShapeFunctions();
return 0;
}
MVertex* MSubLine::getShapeFunctionNode(int i)
{
if(_orig) return _orig->getShapeFunctionNode(i);
return 0;
}
void MSubLine::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
_orig->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwE[3];
getBaseElement()->xyz2uvw(xyz, uvwE);
u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
}
void MSubLine::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
getBaseElement()->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwP[3];
_orig->xyz2uvw(xyz, uvwP);
u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
}
bool MSubLine::isInside(double u, double v, double w)
......@@ -182,42 +575,57 @@ bool MSubLine::isInside(double u, double v, double w)
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;
if(_pts)
{
if(pOrder==_pOrder)
{
*npts = _npts;
*pts = _pts;
return;
}
else
delete [] _pts;
}
// work in the parametric space of the parent element
// (i) get the parametric coordinates of MSubElement vertices in this space
double v_uvw[2][3];
for(int i=0; i<2; ++i)
_pOrder = pOrder;
if(!_orig)
{
MVertex *vi = getVertex(i);
double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
_orig->xyz2uvw(v_xyz, v_uvw[i]);
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &_pts);
*npts = _npts;
*pts = _pts;
return;
}
// (ii) build a MElement on the vertices with these coordinates in this space
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);
// (iii) get the integration points of the new element in its parametric space
l.getIntegrationPoints(pOrder, &nptsi, &ptsi);
// (iv) get the coordinates of these points in the parametric space of parent element
for(int ip=0; ip<nptsi; ++ip)
// work in the parametric space of the parent element
*npts=0;
_pts = new IntPt[getNGQTetPts(pOrder)];
// (i) get the integration points of the base element in its parametric space
IntPt *ptsb;
getBaseElement()->getIntegrationPoints(pOrder, &_npts, &ptsb);
// (ii) get the coordinates of these points in the parametric space of parent element
double u,v,w;
double jac[3][3];
double baseJac, origJac;
for(int i=0; i<_npts; ++i)
{
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;
u = ptsb[i].pt[0];
v = ptsb[i].pt[1];
w = ptsb[i].pt[2];
baseJac = getBaseElement()->getJacobian(u, v, w, jac);
movePointFromElementSpaceToParentSpace(u, v, w);
origJac = _orig->getJacobian(u, v, w, jac);
_pts[*npts + i].pt[0] = u;
_pts[*npts + i].pt[1] = v;
_pts[*npts + i].pt[2] = w;
_pts[*npts + i].weight = ptsb[i].weight * baseJac/origJac;
}
*npts = nptsi;
*pts = _intpt;
*npts = _npts;
*pts = _pts;
}
......@@ -225,9 +633,102 @@ void MSubLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
MSubPoint::~MSubPoint()
{
if(_intpt) delete [] _intpt;
if(_pts) delete [] _pts;
if(_base) delete _base;
}
const nodalBasis* 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 order)
{
if(_orig) _orig->getShapeFunctions(u, v, w, s, order);
}
void MSubPoint::getGradShapeFunctions(double u, double v, double w, double s[][3], int order)
{
if(_orig) _orig->getGradShapeFunctions(u, v, w, s, order);
}
void MSubPoint::getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order)
{
if(_orig) _orig->getHessShapeFunctions(u, v, w, s, order);
}
void MSubPoint::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order)
{
if(_orig) _orig->getThirdDerivativeShapeFunctions(u, v, w, s, order);
}
double MSubPoint::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubPoint::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
{
if(_orig) return _orig->getJacobian(gsf, jac);
return 0;
}
double MSubPoint::getJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getJacobian(u, v, w, jac);
return 0;
}
double MSubPoint::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
if(_orig) return _orig->getPrimaryJacobian(u, v, w, jac);
return 0;
}
int MSubPoint::getNumShapeFunctions()
{
if(_orig) return _orig->getNumShapeFunctions();
return 0;
}
int MSubPoint::getNumPrimaryShapeFunctions()
{
if(_orig) return _orig->getNumPrimaryShapeFunctions();
return 0;
}
MVertex* MSubPoint::getShapeFunctionNode(int i)
{
if(_orig) return _orig->getShapeFunctionNode(i);
return 0;
}
void MSubPoint::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
_orig->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwE[3];
getBaseElement()->xyz2uvw(xyz, uvwE);
u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
}
void MSubPoint::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
{
if(!_orig) return;
SPoint3 p;
getBaseElement()->pnt(u, v, w, p);
double xyz[3] = {p.x(), p.y(), p.z()};
double uvwP[3];
_orig->xyz2uvw(xyz, uvwP);
u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
}
bool MSubPoint::isInside(double u, double v, double w)
{
......@@ -248,20 +749,23 @@ bool MSubPoint::isInside(double u, double v, double w)
void MSubPoint::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
{
// invariable regardless of the order
if(!_intpt)
if(!_pts)
{
if(!_orig) return;
_intpt = new IntPt[1];
_pts = new IntPt[1];
// work in the parametric space of the parent element
MVertex *vi = getVertex(0);
double v_xyz[3] = {vi->x(), vi->y(), vi->z()};
double v_uvw[3];
_orig->xyz2uvw(v_xyz, v_uvw);
_intpt[0].pt[0] = v_uvw[0];
_intpt[0].pt[1] = v_uvw[1];
_intpt[0].pt[2] = v_uvw[2];
_intpt[0].weight = 1;
double jac[3][3];
double origJac = _orig->getJacobian(v_uvw[0], v_uvw[1], v_uvw[2], jac);
_pts[0].pt[0] = v_uvw[0];
_pts[0].pt[1] = v_uvw[1];
_pts[0].pt[2] = v_uvw[2];
_pts[0].weight = 1./origJac;
}
*npts = 1;
*pts = _intpt;
*pts = _pts;
}
......@@ -23,22 +23,43 @@ class MSubTetrahedron : public MTetrahedron
bool _owner;
MElement* _orig;
std::vector<MElement*> _parents;
IntPt *_intpt;
mutable MElement* _base;
int _pOrder;
int _npts;
IntPt *_pts;
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(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), _base(0), _pOrder(-1), _npts(0), _pts(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), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
MSubTetrahedron(const MTetrahedron &tet, bool owner=false, MElement* orig=NULL)
: MTetrahedron(tet), _owner(owner), _orig(orig), _intpt(0) {}
: MTetrahedron(tet), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
~MSubTetrahedron();
virtual int getTypeForMSH() const { return MSH_TET_SUB; }
virtual const nodalBasis* getFunctionSpace(int order=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
// the parametric coordinates are the coordinates in the local parent element
virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
virtual double getJacobian(double u, double v, double w, double jac[3][3]);
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
virtual int getNumShapeFunctions();
virtual int getNumPrimaryShapeFunctions();
virtual MVertex* getShapeFunctionNode(int i);
virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
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; }
......@@ -47,6 +68,9 @@ class MSubTetrahedron : public MTetrahedron
{
_parents = parents; _orig = _parents[0]; _owner = owner;
}
virtual const MElement *getBaseElement() const {if(!_base) _base=new MTetrahedron(*this); return _base;}
virtual MElement *getBaseElement() {if(!_base) _base=new MTetrahedron(*this); return _base;}
};
// A sub triangle, contained in another element
......@@ -56,21 +80,43 @@ class MSubTriangle : public MTriangle
bool _owner;
MElement* _orig;
std::vector<MElement*> _parents;
IntPt *_intpt;
mutable MElement* _base;
int _pOrder;
int _npts;
IntPt *_pts;
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) {}
: MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(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), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
MSubTriangle(const MTriangle &tri, bool owner=false, MElement* orig=NULL)
: MTriangle(tri), _owner(owner), _orig(orig), _intpt(0) {}
: MTriangle(tri), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
~MSubTriangle();
virtual int getTypeForMSH() const { return MSH_TRI_SUB; }
virtual const nodalBasis* getFunctionSpace(int order=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
// the parametric coordinates are the coordinates in the local parent element
virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
virtual double getJacobian(double u, double v, double w, double jac[3][3]);
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
virtual int getNumShapeFunctions();
virtual int getNumPrimaryShapeFunctions();
virtual MVertex* getShapeFunctionNode(int i);
virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
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; }
......@@ -79,6 +125,9 @@ class MSubTriangle : public MTriangle
{
_parents = parents; _orig = _parents[0]; _owner = owner;
}
virtual const MElement *getBaseElement() const {if(!_base) _base=new MTriangle(*this); return _base;}
virtual MElement *getBaseElement() {if(!_base) _base=new MTriangle(*this); return _base;}
};
// A sub line, contained in another element
......@@ -88,22 +137,43 @@ class MSubLine : public MLine
bool _owner;
MElement* _orig;
std::vector<MElement*> _parents;
IntPt *_intpt;
mutable MElement* _base;
int _pOrder;
int _npts;
IntPt *_pts;
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(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), _base(0), _pOrder(-1), _npts(0), _pts(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), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
MSubLine(const MLine &lin, bool owner=false, MElement* orig=NULL)
: MLine(lin), _owner(owner), _orig(orig), _intpt(0) {}
: MLine(lin), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
~MSubLine();
virtual int getTypeForMSH() const { return MSH_LIN_SUB; }
virtual const nodalBasis* getFunctionSpace(int order=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
// the parametric coordinates are the coordinates in the local parent element
virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
virtual double getJacobian(double u, double v, double w, double jac[3][3]);
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
virtual int getNumShapeFunctions();
virtual int getNumPrimaryShapeFunctions();
virtual MVertex* getShapeFunctionNode(int i);
virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
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; }
......@@ -112,6 +182,9 @@ class MSubLine : public MLine
{
_parents = parents; _orig = _parents[0]; _owner = owner;
}
virtual const MElement *getBaseElement() const {if(!_base) _base=new MLine(*this); return _base;}
virtual MElement *getBaseElement() {if(!_base) _base=new MLine(*this); return _base;}
};
// A sub point, contained in another element
......@@ -121,21 +194,43 @@ class MSubPoint : public MPoint
bool _owner;
MElement* _orig;
std::vector<MElement*> _parents;
IntPt *_intpt;
mutable MElement* _base;
int _pOrder;
int _npts;
IntPt *_pts;
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(MVertex *v0, int num=0, int part=0,
bool owner=false, MElement* orig=NULL)
: MPoint(v0, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(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), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
MSubPoint(const MPoint &pt, bool owner=false, MElement* orig=NULL)
: MPoint(pt), _owner(owner), _orig(orig), _intpt(0) {}
: MPoint(pt), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0), _pts(0) {}
~MSubPoint();
virtual int getTypeForMSH() const { return MSH_PNT_SUB; }
virtual const nodalBasis* getFunctionSpace(int order=-1) const;
virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
// the parametric coordinates are the coordinates in the local parent element
virtual void getShapeFunctions(double u, double v, double w, double s[], int order=-1);
virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int order=-1);
virtual void getHessShapeFunctions(double u, double v, double w, double s[][3][3], int order=-1);
virtual void getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3], int order=-1);
virtual double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
virtual double getJacobian(const std::vector<SVector3> &gsf, double jac[3][3]);
virtual double getJacobian(double u, double v, double w, double jac[3][3]);
virtual double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
virtual int getNumShapeFunctions();
virtual int getNumPrimaryShapeFunctions();
virtual MVertex* getShapeFunctionNode(int i);
virtual void movePointFromParentSpaceToElementSpace(double &u, double &v, double &w);
virtual void movePointFromElementSpaceToParentSpace(double &u, double &v, double &w);
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; }
......@@ -144,6 +239,9 @@ class MSubPoint : public MPoint
{
_parents = parents; _orig = _parents[0]; _owner = owner;
}
virtual const MElement *getBaseElement() const {if(!_base) _base=new MPoint(*this); return _base;}
virtual MElement *getBaseElement() {if(!_base) _base=new MPoint(*this); return _base;}
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment