diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 248f98f1895baee7386309bc969a6dab7251a759..adcf36203dcafa8f518a080cadbb026cb0e88fcb 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.59 2008-03-18 19:30:14 geuzaine Exp $
+// $Id: MElement.cpp,v 1.60 2008-03-19 21:22:36 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -133,27 +133,20 @@ double MTetrahedron::etaShapeMeasure()
}
double MTetrahedron::getVolume()
-{
+{
double mat[3][3];
getMat(mat);
return det3x3(mat) / 6.;
}
-bool MTetrahedron::invertmapping(double *p, double *uvw, double tol)
+void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3])
{
- double mat[3][3];
- double b[3], dum;
+ double mat[3][3], b[3], det;
getMat(mat);
- b[0] = p[0] - getVertex(0)->x();
- b[1] = p[1] - getVertex(0)->y();
- b[2] = p[2] - getVertex(0)->z();
- sys3x3(mat, b, uvw, &dum);
- if(uvw[0] >= -tol && uvw[1] >= -tol && uvw[2] >= -tol &&
- uvw[0] <= 1. + tol && uvw[1] <= 1. + tol && uvw[2] <= 1. + tol &&
- 1. - uvw[0] - uvw[1] - uvw[2] > -tol) {
- return true;
- }
- return false;
+ b[0] = xyz[0] - getVertex(0)->x();
+ b[1] = xyz[1] - getVertex(0)->y();
+ b[2] = xyz[2] - getVertex(0)->z();
+ sys3x3(mat, b, uvw, &det);
}
int MHexahedron::getVolumeSign()
@@ -224,6 +217,172 @@ std::string MElement::getInfoString()
return std::string(tmp);
}
+double MElement::getJacobian(double u, double v, double w, double jac[3][3])
+{
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+ double s[3];
+ switch(getDim()){
+ case 3 :
+ for(int i = 0; i < getNumVertices(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ MVertex *p = getVertex(i);
+ jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
+ jac[1][0] += p->x() * s[1]; jac[1][1] += p->y() * s[1]; jac[1][2] += p->z() * s[1];
+ jac[2][0] += p->x() * s[2]; jac[2][1] += p->y() * s[2]; jac[2][2] += p->z() * s[2];
+ }
+ return fabs(jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
+ jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
+ jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
+ case 2 :
+ for(int i = 0; i < getNumVertices(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ MVertex *p = getVertex(i);
+ jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
+ jac[1][0] += p->x() * s[1]; jac[1][1] += p->y() * s[1]; jac[1][2] += p->z() * s[1];
+ }
+ {
+ double a[3], b[3], c[3];
+ a[0]= getVertex(1)->x() - getVertex(0)->x();
+ a[1]= getVertex(1)->y() - getVertex(0)->y();
+ a[2]= getVertex(1)->z() - getVertex(0)->z();
+ b[0]= getVertex(2)->x() - getVertex(0)->x();
+ b[1]= getVertex(2)->y() - getVertex(0)->y();
+ b[2]= getVertex(2)->z() - getVertex(0)->z();
+ prodve(a, b, c);
+ jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ }
+ return sqrt(SQR(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+ SQR(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
+ SQR(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
+ case 1:
+ for(int i = 0; i < getNumVertices(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ MVertex *p = getVertex(i);
+ jac[0][0] += p->x() * s[0]; jac[0][1] += p->y() * s[0]; jac[0][2] += p->z() * s[0];
+ }
+ {
+ double a[3], b[3], c[3];
+ a[0]= getVertex(1)->x() - getVertex(0)->x();
+ a[1]= getVertex(1)->y() - getVertex(0)->y();
+ a[2]= getVertex(1)->z() - getVertex(0)->z();
+ if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
+ (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
+ b[0] = a[1]; b[1] = -a[0]; b[2] = 0.;
+ }
+ else {
+ b[0] = 0.; b[1] = a[2]; b[2] = -a[1];
+ }
+ prodve(a, b, c);
+ jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
+ jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ }
+ return sqrt(SQR(jac[0][0]) + SQR(jac[0][1]) + SQR(jac[0][2]));
+ default:
+ return 1.;
+ }
+}
+
+void MElement::xyz2uvw(double xyz[3], double uvw[3])
+{
+ // general Newton routine for the nonlinear case (more efficient
+ // routines are implemented for simplices, where the basis functions
+ // are linear)
+ uvw[0] = uvw[1] = uvw[2] = 0.;
+ int iter = 1, maxiter = 20;
+ double error = 1., tol = 1.e-6;
+
+ while (error > tol && iter < maxiter){
+ double jac[3][3];
+ if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
+ double xn = 0., yn = 0., zn = 0.;
+ for (int i = 0; i < getNumVertices(); i++) {
+ double s;
+ getShapeFunction(i, uvw[0], uvw[1], uvw[2], s);
+ MVertex *v = getVertex(i);
+ xn += v->x() * s;
+ yn += v->y() * s;
+ zn += v->z() * s;
+ }
+ double inv[3][3];
+ inv3x3(jac, inv);
+ double un = uvw[0] +
+ inv[0][0] * (xyz[0] - xn) + inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
+ double vn = uvw[1] +
+ inv[0][1] * (xyz[0] - xn) + inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn) ;
+ double wn = uvw[2] +
+ inv[0][2] * (xyz[0] - xn) + inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn) ;
+ error = sqrt(SQR(un - uvw[0]) + SQR(vn - uvw[1]) + SQR(wn - uvw[2]));
+ uvw[0] = un;
+ uvw[1] = vn;
+ uvw[2] = wn;
+ iter++ ;
+ }
+}
+
+double MElement::interpolate(double val[], double u, double v, double w, int stride)
+{
+ double sum = 0;
+ int j = 0;
+ for(int i = 0; i < getNumVertices(); i++){
+ double s;
+ getShapeFunction(i, u, v, w, s);
+ sum += val[j] * s;
+ j += stride;
+ }
+ return sum;
+}
+
+void MElement::interpolateGrad(double val[], double u, double v, double w, double f[3],
+ int stride, double invjac[3][3])
+{
+ double dfdu[3] = {0., 0., 0.};
+ int j = 0;
+ for(int i = 0; i < getNumVertices(); i++){
+ double s[3];
+ getGradShapeFunction(i, u, v, w, s);
+ dfdu[0] += val[j] * s[0];
+ dfdu[1] += val[j] * s[1];
+ dfdu[2] += val[j] * s[2];
+ j += stride;
+ }
+ if(invjac){
+ matvec(invjac, dfdu, f);
+ }
+ else{
+ double jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ matvec(inv, dfdu, f);
+ }
+}
+
+void MElement::interpolateCurl(double val[], double u, double v, double w, double f[3],
+ int stride)
+{
+ double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ interpolateGrad(&val[0], u, v, w, fx, stride, inv);
+ interpolateGrad(&val[1], u, v, w, fy, stride, inv);
+ interpolateGrad(&val[2], u, v, w, fz, stride, inv);
+ f[0] = fz[1] - fy[2];
+ f[1] = -(fz[0] - fx[2]);
+ f[2] = fy[0] - fx[1];
+}
+
+double MElement::interpolateDiv(double val[], double u, double v, double w, int stride)
+{
+ double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ interpolateGrad(&val[0], u, v, w, fx, stride, inv);
+ interpolateGrad(&val[1], u, v, w, fy, stride, inv);
+ interpolateGrad(&val[2], u, v, w, fz, stride, inv);
+ return fx[0] + fy[1] + fz[2];
+}
+
void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
int elementary, int physical)
{
@@ -534,35 +693,6 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
#endif
}
-void MTriangleN::jac(double uu, double vv , double j[2][3])
-{
- MTriangle::jac(_order, &(*(_vs.begin())), uu, vv, j);
-}
-
-void MTriangleN::pnt(double uu, double vv, SPoint3 &p)
-{
- MTriangle::pnt(_order, &(*(_vs.begin())), uu, vv, p);
-}
-
-void MTriangle6::jac(double uu, double vv , double j[2][3])
-{
- MTriangle::jac(2, _vs, uu, vv, j);
-}
-
-void MTriangle6::pnt(double uu, double vv, SPoint3 &p){
- MTriangle::pnt(2, _vs, uu, vv, p);
-}
-
-void MTriangle::jac(double uu, double vv, double j[2][3])
-{
- jac(1, 0, uu, vv, j);
-}
-
-void MTriangle::pnt(double uu, double vv, SPoint3 &p)
-{
- MTriangle::pnt(1, 0, uu, vv, p);
-}
-
int MTriangle6::getNumEdgesRep(){ return 3 * 6; }
void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 27de345dcf0c027b00902857d33d7a3505f2f354..222070a36b2d4f7f3982e716c77a7af7b7be0e0a 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -142,6 +142,37 @@ class MElement
// Returns an information string for the element
virtual std::string getInfoString();
+ // Returns the interpolating nodal shape function associated with
+ // node num, evaluated at point (u,v,w) in parametric coordinates
+ virtual void getShapeFunction(int num, double u, double v, double w,
+ double &s) = 0;
+
+ // Returns the gradient of of the nodal shape function associated
+ // with node num, evaluated at point (u,v,w) in parametric
+ // coordinates
+ virtual void getGradShapeFunction(int num, double u, double v, double w,
+ double s[3]) = 0;
+
+ // Returns the Jacobian of the element evaluated at point (u,v,w) in
+ // parametric coordinates
+ double getJacobian(double u, double v, double w, double jac[3][3]);
+
+ // Inverts the parametrisation
+ virtual void xyz2uvw(double xyz[3], double uvw[3]);
+
+ // Tests if a point, given in parametric coordinates, belongs to the
+ // element
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8) = 0;
+
+ // Interpolate the given nodal data (resp. its gradient, curl and
+ // divergence) at point (u,v,w) in parametric coordinates
+ double interpolate(double val[], double u, double v, double w, int stride=1);
+ void interpolateGrad(double val[], double u, double v, double w, double f[3],
+ int stride=1, double invjac[3][3]=0);
+ void interpolateCurl(double val[], double u, double v, double w, double f[3],
+ int stride=3);
+ double interpolateDiv(double val[], double u, double v, double w, int stride=3);
+
// IO routines
virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
int num=0, int elementary=1, int physical=1);
@@ -217,6 +248,28 @@ class MLine : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 0.5 * (1. - u); break;
+ case 1 : s = 0.5 * (1. + u); break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.5; s[1] = 0.; s[2] = 0.; break;
+ case 1 : s[0] = 0.5; s[1] = 0.; s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < -(1. + tol) || u > (1. + tol))
+ return false;
+ return true;
+ }
};
class MLine3 : public MLine {
@@ -274,7 +327,7 @@ class MLineN : public MLine {
{
for(unsigned int i = 2; i < v.size(); i++)
_vs.push_back(v[i]);
- for(unsigned int i = 0 ; i < _vs.size(); i++)
+ for(unsigned int i = 0; i < _vs.size(); i++)
_vs[i]->setPolynomialOrder(_vs.size() + 1);
}
~MLineN(){}
@@ -304,8 +357,6 @@ class MTriangle : public MElement {
protected:
MVertex *_v[3];
public :
- virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][3]);
- virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
: MElement(num, part)
{
@@ -318,13 +369,6 @@ class MTriangle : public MElement {
}
~MTriangle(){}
virtual int getDim(){ return 2; }
- void getMat(double mat[2][2])
- {
- mat[0][0] = _v[1]->x() - _v[0]->x();
- mat[0][1] = _v[2]->x() - _v[0]->x();
- mat[1][0] = _v[1]->y() - _v[0]->y();
- mat[1][1] = _v[2]->y() - _v[0]->y();
- }
virtual double gammaShapeMeasure();
virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
@@ -369,8 +413,40 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual void jac(double u, double v, double j[2][3]);
- virtual void pnt(double u, double v, SPoint3&);
+ virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][3]);
+ virtual void jac(double u, double v, double j[2][3])
+ {
+ jac(1, 0, u, v, j);
+ }
+ virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
+ virtual void pnt(double u, double v, SPoint3 &p)
+ {
+ pnt(1, 0, u, v, p);
+ }
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num){
+ case 0 : s = 1. - u - v; break;
+ case 1 : s = u ; break;
+ case 2 : s = v; break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -1.; s[1] = -1.; s[2] = 0.; break;
+ case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
+ case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
+ return false;
+ return true;
+ }
};
class MTriangle6 : public MTriangle {
@@ -429,8 +505,14 @@ class MTriangle6 : public MTriangle {
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
- virtual void jac ( double u, double v , double j[2][3]) ;
- virtual void pnt(double u, double v, SPoint3&);
+ virtual void jac(double u, double v, double j[2][3])
+ {
+ MTriangle::jac(2, _vs, u, v, j);
+ }
+ virtual void pnt(double u, double v, SPoint3 &p)
+ {
+ MTriangle::pnt(2, _vs, u, v, p);
+ }
};
class MTriangleN : public MTriangle {
@@ -461,7 +543,7 @@ class MTriangleN : public MTriangle {
}
~MTriangleN(){}
virtual int getPolynomialOrder(){ return _order; }
- virtual int getNumVertices(){ return 3 + _vs.size() ; }
+ virtual int getNumVertices(){ return 3 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
virtual int getNumFaceVertices()
{
@@ -502,8 +584,14 @@ class MTriangleN : public MTriangle {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
- virtual void jac(double u, double v, double jac[2][3]);
- virtual void pnt(double u, double v, SPoint3&);
+ virtual void jac(double u, double v, double j[2][3])
+ {
+ MTriangle::jac(_order, &(*(_vs.begin())), u, v, j);
+ }
+ virtual void pnt(double u, double v, SPoint3 &p)
+ {
+ MTriangle::pnt(_order, &(*(_vs.begin())), u, v, p);
+ }
};
class MQuadrangle : public MElement {
@@ -560,6 +648,32 @@ class MQuadrangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
}
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 0.25 * (1. - u) * (1. - v); break;
+ case 1 : s = 0.25 * (1. + u) * (1. - v); break;
+ case 2 : s = 0.25 * (1. + u) * (1. + v); break;
+ case 3 : s = 0.25 * (1. - u) * (1. + v); break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.25 * (1. - v); s[1] = -0.25 * (1. - u); s[2] = 0.; break;
+ case 1 : s[0] = 0.25 * (1. - v); s[1] = -0.25 * (1. + u); s[2] = 0.; break;
+ case 2 : s[0] = 0.25 * (1. + v); s[1] = 0.25 * (1. + u); s[2] = 0.; break;
+ case 3 : s[0] = -0.25 * (1. + v); s[1] = 0.25 * (1. - u); s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol))
+ return false;
+ return true;
+ }
};
class MQuadrangle8 : public MQuadrangle {
@@ -744,7 +858,7 @@ class MTetrahedron : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- void getMat(double mat[3][3])
+ void getMat(double mat[3][3])
{
mat[0][0] = _v[1]->x() - _v[0]->x();
mat[0][1] = _v[2]->x() - _v[0]->x();
@@ -760,8 +874,33 @@ class MTetrahedron : public MElement {
virtual int getVolumeSign(){ return (getVolume() >= 0) ? 1 : -1; }
virtual double gammaShapeMeasure();
virtual double etaShapeMeasure();
- // returns true if the point lies inside the tet
- bool invertmapping(double *p, double *uvw, double tol = 1.e-8);
+ void xyz2uvw(double xyz[3], double uvw[3]);
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 1. - u - v - w; break;
+ case 1 : s = u ; break;
+ case 2 : s = v ; break;
+ case 3 : s = w; break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -1.; s[1] = -1.; s[2] = -1.; break;
+ case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
+ case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
+ case 3 : s[0] = 0.; s[1] = 0.; s[2] = 1.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < (-tol) || v < (-tol) || w < (-tol) || u > ((1. + tol) - v - w))
+ return false;
+ return true;
+ }
};
class MTetrahedron10 : public MTetrahedron {
@@ -923,6 +1062,57 @@ class MHexahedron : public MElement {
tmp = _v[4]; _v[4] = _v[6]; _v[6] = tmp;
}
virtual int getVolumeSign();
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = (1. - u) * (1. - v) * (1. - w) * 0.125; break;
+ case 1 : s = (1. + u) * (1. - v) * (1. - w) * 0.125; break;
+ case 2 : s = (1. + u) * (1. + v) * (1. - w) * 0.125; break;
+ case 3 : s = (1. - u) * (1. + v) * (1. - w) * 0.125; break;
+ case 4 : s = (1. - u) * (1. - v) * (1. + w) * 0.125; break;
+ case 5 : s = (1. + u) * (1. - v) * (1. + w) * 0.125; break;
+ case 6 : s = (1. + u) * (1. + v) * (1. + w) * 0.125; break;
+ case 7 : s = (1. - u) * (1. + v) * (1. + w) * 0.125; break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.125 * (1. - v) * (1. - w);
+ s[1] = -0.125 * (1. - u) * (1. - w);
+ s[2] = -0.125 * (1. - u) * (1. - v); break;
+ case 1 : s[0] = 0.125 * (1. - v) * (1. - w);
+ s[1] = -0.125 * (1. + u) * (1. - w);
+ s[2] = -0.125 * (1. + u) * (1. - v); break;
+ case 2 : s[0] = 0.125 * (1. + v) * (1. - w);
+ s[1] = 0.125 * (1. + u) * (1. - w);
+ s[2] = -0.125 * (1. + u) * (1. + v); break;
+ case 3 : s[0] = -0.125 * (1. + v) * (1. - w);
+ s[1] = 0.125 * (1. - u) * (1. - w);
+ s[2] = -0.125 * (1. - u) * (1. + v); break;
+ case 4 : s[0] = -0.125 * (1. - v) * (1. + w);
+ s[1] = -0.125 * (1. - u) * (1. + w);
+ s[2] = 0.125 * (1. - u) * (1. - v); break;
+ case 5 : s[0] = 0.125 * (1. - v) * (1. + w);
+ s[1] = -0.125 * (1. + u) * (1. + w);
+ s[2] = 0.125 * (1. + u) * (1. - v); break;
+ case 6 : s[0] = 0.125 * (1. + v) * (1. + w);
+ s[1] = 0.125 * (1. + u) * (1. + w);
+ s[2] = 0.125 * (1. + u) * (1. + v); break;
+ case 7 : s[0] = -0.125 * (1. + v) * (1. + w);
+ s[1] = 0.125 * (1. - u) * (1. + w);
+ s[2] = 0.125 * (1. - u) * (1. + v); break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < -(1. + tol) || v < -(1. + tol) || w < -(1. + tol) ||
+ u > (1. + tol) || v > (1. + tol) || w > (1. + tol))
+ return false;
+ return true;
+ }
};
class MHexahedron20 : public MHexahedron {
@@ -1200,6 +1390,49 @@ class MPrism : public MElement {
tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
}
virtual int getVolumeSign();
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = (1. - u - v) * (1. - w) * 0.5; break;
+ case 1 : s = u * (1. - w) * 0.5; break;
+ case 2 : s = v * (1. - w) * 0.5; break;
+ case 3 : s = (1. - u - v) * (1. + w) * 0.5; break;
+ case 4 : s = u * (1. + w) * 0.5; break;
+ case 5 : s = v * (1. + w) * 0.5; break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.5 * (1. - w) ;
+ s[1] = -0.5 * (1. - w) ;
+ s[2] = -0.5 * (1. - u - v); break;
+ case 1 : s[0] = 0.5 * (1. - w) ;
+ s[1] = 0. ;
+ s[2] = -0.5 * u ; break;
+ case 2 : s[0] = 0. ;
+ s[1] = 0.5 * (1. - w) ;
+ s[2] = -0.5 * v ; break;
+ case 3 : s[0] = -0.5 * (1. + w) ;
+ s[1] = -0.5 * (1. + w) ;
+ s[2] = 0.5 * (1. - u - v); break;
+ case 4 : s[0] = 0.5 * (1. + w) ;
+ s[1] = 0. ;
+ s[2] = 0.5 * u ; break;
+ case 5 : s[0] = 0. ;
+ s[1] = 0.5 * (1. + w) ;
+ s[2] = 0.5 * v ; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(w > (1. + tol) || w < -(1. + tol) || u < (1. + tol)
+ || v < (1. + tol) || u > ((1. + tol) - v))
+ return false;
+ return true;
+ }
};
class MPrism15 : public MPrism {
@@ -1439,6 +1672,55 @@ class MPyramid : public MElement {
MVertex *tmp = _v[0]; _v[0] = _v[2]; _v[2] = tmp;
}
virtual int getVolumeSign();
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ double r;
+ if(w != 1. && num != 4) r = u * v * w / (1. - w);
+ else r = 0.;
+ switch(num) {
+ case 0 : s = 0.25 * ((1. - u) * (1. - v) - w + r); break;
+ case 1 : s = 0.25 * ((1. + u) * (1. - v) - w - r); break;
+ case 2 : s = 0.25 * ((1. + u) * (1. + v) - w + r); break;
+ case 3 : s = 0.25 * ((1. - u) * (1. + v) - w - r); break;
+ case 4 : s = w; break;
+ default : s = 0.; break;
+ }
+ }
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ if(w == 1. && num != 4) {
+ s[0] = 0.25;
+ s[1] = 0.25;
+ s[2] = -0.25;
+ }
+ else{
+ switch(num) {
+ case 0 : s[0] = 0.25 * (-(1. - v) + v * w / (1. - w));
+ s[1] = 0.25 * (-(1. - u) + u * w / (1. - w));
+ s[2] = 0.25 * (-1. + u * v / (1. - w) / (1. - w)); break;
+ case 1 : s[0] = 0.25 * ( (1. - v) + v * w / (1. - w));
+ s[1] = 0.25 * (-(1. + u) + u * w / (1. - w));
+ s[2] = 0.25 * (-1. + u * v / (1. - w) / (1. - w)); break;
+ case 2 : s[0] = 0.25 * ( (1. + v) + v * w / (1. - w));
+ s[1] = 0.25 * ( (1. + u) + u * w / (1. - w));
+ s[2] = 0.25 * (-1. + u * v / (1. - w) / (1. - w)); break;
+ case 3 : s[0] = 0.25 * (-(1. + v) + v * w / (1. - w));
+ s[1] = 0.25 * ( (1. - u) + u * w / (1. - w));
+ s[2] = 0.25 * (-1. + u * v / (1. - w) / (1. - w)); break;
+ case 4 : s[0] = 0.;
+ s[1] = 0.;
+ s[2] = 1.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ }
+ virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+ {
+ if(u < (w - (1. + tol)) || u > ((1. + tol) - w) || v < (w - (1. + tol)) ||
+ v > ((1. + tol) - w) || w < (-tol) || w > (1. + tol))
+ return false;
+ return true;
+ }
};
class MPyramid13 : public MPyramid {
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 82df613ff9003101164d1df4a0d8925d4e53d89a..f903157c0237451014e0e20ead27493ac1e867b8 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.41 2008-03-18 19:30:14 geuzaine Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.42 2008-03-19 21:22:36 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -811,8 +811,8 @@ void insertVerticesInRegion (GRegion *gr)
double center[3];
worst->circumcenter(center);
double uvw[3];
- bool inside = worst->tet()->invertmapping(center, uvw);
- if(inside){
+ worst->tet()->xyz2uvw(center, uvw);
+ if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setNum(NUM++);
double lc1 =
@@ -820,8 +820,8 @@ void insertVerticesInRegion (GRegion *gr)
uvw[0] * vSizes[worst->tet()->getVertex(1)->getNum()] +
uvw[1] * vSizes[worst->tet()->getVertex(2)->getNum()] +
uvw[2] * vSizes[worst->tet()->getVertex(3)->getNum()];
- double lc = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
- // double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
+ double lc = BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]);
+ // double lc = std::min(lc1, BGM_MeshSize(gr, 0, 0, center[0], center[1], center[2]));
vSizes.push_back(lc1);
vSizesBGM.push_back(lc);
// compute mesh spacing there