00001
00029 #ifndef TETRAHEDRON
00030 #define TETRAHEDRON
00031
00032 #include <algorithm>
00033 #include <cassert>
00034
00035 #include "ElementGeometry.h"
00036 #include "Triangle.h"
00037
00060 #define TET_SPD 3
00061
00062 class Tetrahedron: public AbstractGeom<TET_SPD>
00063 {
00064 private:
00065 static const double ParamCoord[];
00066 static const size_t FaceNodes[];
00067
00068
00069 public:
00070
00071 Tetrahedron (const std::vector<double>& globalCoordVec, const std::vector<GlobalNodalIndex>& connectivity)
00072 :AbstractGeom<TET_SPD> (globalCoordVec, connectivity) {
00073 assert (connectivity.size () == 4);
00074 }
00075
00076
00077 Tetrahedron(const Tetrahedron & that) : AbstractGeom<TET_SPD>(that) {
00078 }
00079
00080 virtual Tetrahedron* clone() const {
00081 return new Tetrahedron(*this);
00082 }
00083
00084
00086 inline const size_t getNumVertices() const { return 4; }
00087
00088
00090 inline const std::string getPolytopeName() const { return "TETRAHEDRON"; }
00091
00093 inline const size_t getParametricDimension() const { return 3; }
00094
00096 inline const size_t getEmbeddingDimension() const { return 3; }
00097
00099 inline size_t getNumFaces() const { return 4; }
00100
00104 void map(const double *X, double *Y) const {
00105 const size_t sd = AbstractGeom<TET_SPD>::getSpatialDimension ();
00106
00107 for(size_t i=0; i<sd; i++)
00108 Y[i] =
00109 X[0]*AbstractGeom<TET_SPD>::getCoordinate(0,i) +
00110 X[1]*AbstractGeom<TET_SPD>::getCoordinate(1,i) +
00111 X[2]*AbstractGeom<TET_SPD>::getCoordinate(3,i) +
00112 (1.0-X[0]-X[1]-X[2])*AbstractGeom<TET_SPD>::getCoordinate(2,i);
00113 }
00114
00120 void dMap(const double *X, double *DY, double &Jac) const {
00121 const size_t sd = AbstractGeom<TET_SPD>::getSpatialDimension ();
00122
00123 for(size_t i=0; i<sd; i++)
00124 {
00125 DY[i] = AbstractGeom<TET_SPD>::getCoordinate(0,i)-AbstractGeom<TET_SPD>::getCoordinate(2,i);
00126 DY[sd*1+i] = AbstractGeom<TET_SPD>::getCoordinate(1,i)-AbstractGeom<TET_SPD>::getCoordinate(2,i);
00127 DY[sd*2+i] = AbstractGeom<TET_SPD>::getCoordinate(3,i)-AbstractGeom<TET_SPD>::getCoordinate(2,i);
00128 }
00129
00130 Jac = 0.;
00131 for(size_t i=0; i<sd; i++)
00132 {
00133 size_t i1 = (i+1)%sd;
00134 size_t i2 = (i+2)%sd;
00135 Jac += DY[i]*(DY[1*sd+i1]*DY[2*sd+i2] - DY[2*sd+i1]*DY[1*sd+i2]);
00136 }
00137
00138 Jac = fabs(Jac);
00139 }
00140
00145 Triangle<3> * getFaceGeometry(size_t e) const {
00146
00147 if(e>=0 && e<=3) {
00148 std::vector<GlobalNodalIndex> conn(FaceNodes + 3*e, FaceNodes + 3*e + 2);
00149
00150 return new Triangle<TET_SPD> (AbstractGeom<TET_SPD>::getGlobalCoordVec (), conn);
00151
00152
00153
00154
00155 }
00156 else {
00157 std::cout<<"\nTetrahedron::getFaceGeometry() : Request for invalid face. Quitting...\n";
00158 exit(1);
00159 }
00160 }
00161
00163 const double getInRadius(void) const {
00164 double a[3],b[3],c[3],o,t[3],d[3],t1[3],t2[3], t3[3];
00165 for(size_t i=0;i<3;i++) {
00166 o = AbstractGeom<TET_SPD>::getCoordinate(0,i);
00167 a[i]=AbstractGeom<TET_SPD>::getCoordinate(1,i) - o;
00168 b[i]=AbstractGeom<TET_SPD>::getCoordinate(2,i) - o;
00169 c[i]=AbstractGeom<TET_SPD>::getCoordinate(3,i) - o;
00170 }
00171 cross(b,c,t);
00172 cross(c,a,d);
00173 d[0]+=t[0];
00174 d[1]+=t[1];
00175 d[2]+=t[2];
00176 cross(a,b,t);
00177 d[0]+=t[0];
00178 d[1]+=t[1];
00179 d[2]+=t[2];
00180 cross(b,c,t);
00181 cross(b,c,t1);
00182 cross(c,a,t2);
00183 cross(a,b,t3);
00184 double rv = (dot(a,t)/(mag(t1)+mag(t2)+mag(t3)+ mag(d)));
00185
00186 return(rv);
00187 }
00188
00190 const double getOutRadius(void) const {
00191 double x[4],y[4],z[4],r2[4],ones[4];
00192 double M11, M12, M13, M14, M15;
00193 double **a;
00194 a = new double*[4];
00195
00196 for(size_t i = 0; i < 4; i++) {
00197 x[i] = AbstractGeom<TET_SPD>::getCoordinate(i,0);
00198 y[i] = AbstractGeom<TET_SPD>::getCoordinate(i,1);
00199 z[i] = AbstractGeom<TET_SPD>::getCoordinate(i,2);
00200 r2[i] = x[i]*x[i] + y[i]*y[i] + z[i]*z[i];
00201 ones[i] = 1.0;
00202 }
00203 a[0] = x;
00204 a[1] = y;
00205 a[2] = z;
00206 a[3] = ones;
00207 M11 = determinant(a,4);
00208 a[0] = r2;
00209 M12 = determinant(a,4);
00210 a[1] = x;
00211 M13 = determinant(a,4);
00212 a[2] = y;
00213 M14 = determinant(a,4);
00214 a[3] = z;
00215 M15 = determinant(a,4);
00216
00217 double x0,y0,z0;
00218 x0 = 0.5 * M12/M11;
00219 y0 = -0.5 * M13/M11;
00220 z0 = 0.5 * M14/M11;
00221
00222 delete[] a;
00223 return(sqrt(x0*x0 + y0*y0 + z0*z0 - M15/M11));
00224 }
00225
00226
00231 virtual void computeNormal (size_t e, std::vector<double>& vNormal) const {
00232 const size_t sd = getSpatialDimension();
00233
00234 size_t n0, n1, n2;
00235
00236 n0 = FaceNodes[3*e]; n1 = FaceNodes[3*e+1]; n2 = FaceNodes[3*e+2];
00237
00238
00239 double p0[sd], p1[sd], p2[sd];
00240
00241 map(&ParamCoord[3*n0], p0);
00242 map(&ParamCoord[3*n1], p1);
00243 map(&ParamCoord[3*n2], p2);
00244
00245 double L01[sd];
00246 double L02[sd];
00247 for(size_t k=0; k<sd; k++)
00248 {
00249 L01[k] = p1[k]-p0[k];
00250 L02[k] = p2[k]-p0[k];
00251 }
00252
00253 if(vNormal.size() < sd) vNormal.resize(sd);
00254
00255 double vnorm2=0.;
00256 for(size_t k=0; k<sd; k++)
00257 {
00258 size_t n1 = (k+1)%sd;
00259 size_t n2 = (k+2)%sd;
00260 vNormal[k] = L01[n1]*L02[n2] - L01[n2]*L02[n1];
00261 vnorm2 += vNormal[k]*vNormal[k];
00262 }
00263
00264 for(size_t k=0; k<sd; k++)
00265 vNormal[k] /= sqrt(vnorm2);
00266
00267 }
00268
00269 protected:
00272 static double mag(const double *a) {
00273 double rv = sqrt(dot(a,a));
00274 return(rv);
00275 };
00276
00280 static double dot(const double *a,const double *b) {
00281 return(a[0]*b[0] + a[1]*b[1] + a[2]*b[2]);
00282 };
00283
00287 static void cross(const double *a,const double *b, double *rv) {
00288 rv[0] = a[1]*b[2] - a[2]*b[1];
00289 rv[1] = a[2]*b[0] - a[0]*b[2];
00290 rv[2] = a[0]*b[1] - a[1]*b[0];
00291 };
00292
00296 static double determinant(double **a,size_t n) {
00297 size_t i,j,j1,j2;
00298 double det = 0;
00299 double **m = NULL;
00300
00301 if (n < 1) {
00302
00303 } else if (n == 1) {
00304 det = a[0][0];
00305 } else if (n == 2) {
00306 det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
00307 } else {
00308 det = 0;
00309 for (j1=0;j1<n;j1++) {
00310 m = new double*[n-1];
00311 for (i=0;i<n-1;i++) {
00312 m[i] = new double[n-1];
00313 }
00314 for (i=1;i<n;i++) {
00315 j2 = 0;
00316 for (j=0;j<n;j++) {
00317 if (j == j1) {
00318 continue;
00319 }
00320 m[i-1][j2] = a[i][j];
00321 j2++;
00322 }
00323 }
00324 det += pow(-1.0,1.0+j1+1.0) * a[0][j1] * determinant(m,n-1);
00325 for (i=0;i<n-1;i++) {
00326 delete[] m[i];
00327 }
00328 delete[] m;
00329 }
00330 }
00331 return(det);
00332 }
00333
00334
00335 };
00336
00337 #endif
