MechSys: fem/ele/tet10.h Source File

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 #ifndef MECHSYS_FEM_TET10_H
00023 #define MECHSYS_FEM_TET10_H
00024 
00025 #ifdef HAVE_CONFIG_H
00026   #include "config.h"
00027 #else
00028   #ifndef REAL
00029     #define REAL double
00030   #endif
00031 #endif
00032 
00033 #include "fem/ele/element.h"
00034 #include "linalg/vector.h"
00035 #include "linalg/matrix.h"
00036 #include "linalg/laexpr.h"
00037 
00038 namespace FEM
00039 {
00040 
00041 class Debug;
00042 
00043 // Tet10 Constants {{{
00044 /* IPS[i,j] i=>ip_number, j=>r,s,t,w
00045  * ip=0 {{ r0, s0, t0, w0 },
00046  * ip=1  { r1, s1, t1, w1 },
00047  * ip=2  { r2, s2, t2, w2 }}
00048  */
00049 const int TET10_NNODES=10;
00050 const int TET10_NINTPTS=4;
00051 const int TET10_NFACENODES=6;
00052 const int TET10_NFACEINTPTS=3;
00053 const Element::IntegPoint TET10_INTPTS[]=
00054 {{ 0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.04166666666666667 },
00055  { 0.1381966011250105, 0.5854101966249685, 0.1381966011250105, 0.04166666666666667 },
00056  { 0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.04166666666666667 },
00057  { 0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.04166666666666667 }};
00058 const Element::IntegPoint TET10_FACEINTPTS[]=
00059 {{ 0.166666666666667,  0.166666666666667,  0.0,  0.166666666666667 }, 
00060  { 0.666666666666667,  0.166666666666667,  0.0,  0.166666666666667 },
00061  { 0.166666666666667,  0.666666666666667,  0.0,  0.166666666666667 }};
00062 // }}}
00063 
00064 class Tet10 : public virtual Element
00065 {
00066     friend class FEM::Debug;
00067 public:
00068     // Constructor
00069     Tet10()
00070     {
00071         // Setup nodes number
00072                _n_nodes = TET10_NNODES;
00073              _n_int_pts = TET10_NINTPTS;
00074           _n_face_nodes = TET10_NFACENODES;
00075         _n_face_int_pts = TET10_NFACEINTPTS;
00076 
00077         // Allocate nodes (connectivity)
00078         _connects.resize(_n_nodes);
00079 
00080         // Setup pointer to the array of Integration Points
00081         _a_int_pts = TET10_INTPTS;
00082     }
00083     // Destructor
00084     virtual ~Tet10() {}
00085     // Methods
00086     void Shape      (REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape)       const;
00087     void Derivs     (REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs)      const;
00088     REAL BoundDistance(REAL r, REAL s, REAL t) const;
00089     int VTKCellType() const { return 24; } // VTK_QUADRATIC_TETRA
00090 private:
00091 protected:
00092     void _distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects,
00093                                          REAL                 const   FaceValue   ,
00094                                          LinAlg::Vector<REAL>       & NodalValues ) const;
00095     void _face_shape (REAL r, REAL s,         LinAlg::Vector<REAL> & AFaceShape)  const;
00096     void _face_derivs(REAL r, REAL s,         LinAlg::Matrix<REAL> & AFaceDerivs) const;
00097 }; // class Tet10
00098 
00099 
00101 
00102 
00103 inline void Tet10::_distrib_val_to_face_nodal_vals(Array<Node*>  const & FaceConnects, // {{{   // In:  Array of ptrs to face nodes
00104                                                   REAL           const   FaceValue   ,          // In:  A value applied on a face to be converted to nodes
00105                                                   LinAlg::Vector<REAL> & NodalValues ) const    // Out: The resultant nodal values
00106 {
00107     //Dimensioning NodalValues
00108     NodalValues.Resize(_n_face_nodes);
00109     NodalValues.SetValues(0.0);
00110     LinAlg::Matrix<REAL> J;                           // Jacobian matrix. size = 2 x 3
00111     LinAlg::Vector<REAL> face_shape(_n_face_nodes); // Shape functions of a face. size = _n_face_nodes
00112     //Integration along the face
00113     for (int i=0; i<_n_face_int_pts; i++)
00114     {
00115         REAL r = TET10_FACEINTPTS[i].r;
00116         REAL s = TET10_FACEINTPTS[i].s;
00117         REAL w = TET10_FACEINTPTS[i].w;
00118         _face_shape(r, s, face_shape);
00119         _face_jacobian(FaceConnects, r, s, J);
00120         //REAL det_J = J.Det();
00121         NodalValues += FaceValue*face_shape*det(J)*w;
00122         //LinAlg::Axpy(FaceValue*det_J*w, face_shape, NodalValues); // NodalValues += FaceValue*face_shape*det_J*W
00123     }
00124 } // }}}
00125 
00126 inline void Tet10::_face_shape(REAL r, REAL s, LinAlg::Vector<REAL> & FaceShape) const // {{{
00127 {
00128     /*           s
00129      *           ^
00130      *           |
00131      *         2 @.
00132      *           | '.
00133      *           |   '. 4
00134      *         5 @     @.  
00135      *           |       '.
00136      *           |         '.
00137      *           @-----@-----@-> r
00138      *           0     3     1
00139      */
00140     FaceShape.Resize(6);
00141 
00142     REAL u = 1.0 - r - s;
00143 
00144     FaceShape(0) = u*(2.0*u - 1.0);
00145     FaceShape(1) = r*(2.0*r - 1.0);
00146     FaceShape(2) = s*(2.0*s - 1.0);
00147     FaceShape(3) = 4.0 * r * u;
00148     FaceShape(4) = 4.0 * r * s;
00149     FaceShape(5) = 4.0 * s * u;
00150 } // }}}
00151 
00152 inline void Tet10::_face_derivs(REAL r, REAL s, LinAlg::Matrix<REAL> & FaceDerivs) const // {{{
00153 {
00154     /*           _     _ T
00155      *          |  dNi  |
00156      * Derivs = |  ---  |   , where cj = r, s
00157      *          |_ dcj _|
00158      *
00159      * Derivs(j,i), j=>local coordinate and i=>shape function
00160      */
00161     FaceDerivs.Resize(2,6);
00162 
00163     // r-derivatives: dN0/dr to dN5/dr
00164     FaceDerivs(0,0) =  4.0*r + 4.0*s - 3.0;
00165     FaceDerivs(0,1) =  4.0*r - 1.0;
00166     FaceDerivs(0,2) =  0.0;
00167     FaceDerivs(0,3) =  4.0 - 8.0*r - 4.0*s;
00168     FaceDerivs(0,4) =  4.0*s;
00169     FaceDerivs(0,5) = -4.0*s;
00170     
00171     // s-derivatives: dN0/ds to dN5/ds
00172     FaceDerivs(1,0) =  4.0*r + 4.0*s - 3.0;
00173     FaceDerivs(1,1) =  0.0;
00174     FaceDerivs(1,2) =  4.0*s - 1.0;
00175     FaceDerivs(1,3) = -4.0*r;
00176     FaceDerivs(1,4) =  4.0*r;
00177     FaceDerivs(1,5) =  4.0 - 8.0*s - 4.0*r;
00178 } // }}}
00179 
00180 inline void  Tet10::Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape) const // {{{
00181 {
00182     /*                         
00183      *                         t
00184      *                         |
00185      *                         |                                           
00186      *                         | 3                                         
00187      *                         @,                                          
00188      *                        /|`                                          
00189      *                        ||  `,                                       
00190      *                       / |    ',                                     
00191      *                       | |      \                                    
00192      *                      /  |       `.                                  
00193      *                      |  |         `,  9                             
00194      *                     /   @ 7         `@                              
00195      *                     |   |             \                             
00196      *                    /    |              `.                           
00197      *                    |    |                ',                         
00198      *                 8 @     |                  \                        
00199      *                   |     @.,,_       6       `.                      
00200      *                  |     / 0   ``'-.,,@_        `.                    
00201      *                  |    /              ``''-.,,_  ', 2                
00202      *                 |    /                        ``'@.,,,              
00203      *                 |   '                       ,.-``     ``''- s       
00204      *                |  ,@ 4                 _,-'`                        
00205      *                ' /                 ,.'`                             
00206      *               | /             _.@``                                 
00207      *               '/          ,-'`   5                                  
00208      *              |/      ,.-``                                          
00209      *              /  _,-``                                               
00210      *            .@ '`                                                    
00211      *           / 1                                                       
00212      *          /                                                          
00213      *         /                                                           
00214      *        r                                                            
00215      */
00216     Shape.Resize(10);
00217 
00218     REAL u = 1.0 - r - s - t;
00219     
00220     // corners
00221     Shape(0) = u*(2.0*u - 1.0);
00222     Shape(1) = r*(2.0*r - 1.0);
00223     Shape(2) = s*(2.0*s - 1.0);
00224     Shape(3) = t*(2.0*t - 1.0);
00225     
00226     // midedge
00227     Shape(4) = 4.0 * u * r;
00228     Shape(5) = 4.0 * r * s;
00229     Shape(6) = 4.0 * s * u;
00230     Shape(7) = 4.0 * u * t;
00231     Shape(8) = 4.0 * r * t;
00232     Shape(9) = 4.0 * s * t;
00233 } // }}}
00234 
00235 inline void Tet10::Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs) const // {{{
00236 {
00237     /*           _     _ T
00238      *          |  dNi  |
00239      * Derivs = |  ---  |   , where cj = r, s, t
00240      *          |_ dcj _|
00241      *
00242      * Derivs(j,i), j=>local coordinate and i=>shape function
00243      */
00244     Derivs.Resize(3,10);
00245 
00246     // r-derivatives: dN0/dr to dN9/dr
00247     Derivs(0,0) =  4.0*(r + s + t) - 3.0;
00248     Derivs(0,1) =  4.0*r - 1.0;
00249     Derivs(0,2) =  0.0;
00250     Derivs(0,3) =  0.0;
00251     Derivs(0,4) =  4.0 - 8.0*r - 4.0*s - 4.0*t;
00252     Derivs(0,5) =  4.0*s;
00253     Derivs(0,6) = -4.0*s;
00254     Derivs(0,7) = -4.0*t;
00255     Derivs(0,8) =  4.0*t;
00256     Derivs(0,9) =  0.0;
00257     
00258     // s-derivatives: dN0/ds to dN9/ds
00259     Derivs(1,0) =  4.0*(r + s + t) - 3.0;
00260     Derivs(1,1) =  0.0;
00261     Derivs(1,2) =  4.0*s - 1.0;
00262     Derivs(1,3) =  0.0;
00263     Derivs(1,4) = -4.0*r;
00264     Derivs(1,5) =  4.0*r;
00265     Derivs(1,6) =  4.0 - 4.0*r - 8.0*s - 4.0*t;
00266     Derivs(1,7) = -4.0*t;
00267     Derivs(1,8) =  0.0;
00268     Derivs(1,9) =  4.0*t;
00269     
00270     // t-derivatives: dN0/dt to dN9/dt
00271     Derivs(2,0) =  4.0*(r + s + t) - 3.0;
00272     Derivs(2,1) =  0.0;
00273     Derivs(2,2) =  0.0;
00274     Derivs(2,3) =  4.0*t - 1.0;
00275     Derivs(2,4) = -4.0*r;
00276     Derivs(2,5) =  0.0;
00277     Derivs(2,6) = -4.0*s;
00278     Derivs(2,7) =  4.0 - 4.0*r - 4.0*s - 8.0*t;
00279     Derivs(2,8) =  4.0*r;
00280     Derivs(2,9) =  4.0*s;
00281 } // }}}
00282 
00283 inline REAL Tet10::BoundDistance(REAL r, REAL s, REAL t) const //{{{
00284 {
00285     return std::min(std::min(std::min(r, s), t), 1-r-s-t) ;
00286 } //}}}
00287 
00288 }; // namespace FEM
00289 
00290 #endif // MECHSYS_FEM_TET10_H
00291 
00292 // vim:fdm=marker