MechSys: fem/ele/hex8.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_HEX8_H
00023 #define MECHSYS_FEM_HEX8_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/lawrap.h"
00037 
00038 namespace FEM
00039 {
00040 
00041 // Hex8 Constants // {{{
00042 const int HEX8_NNODES=8;
00043 const int HEX8_NINTPTS=8;
00044 const int HEX8_NFACENODES=4;
00045 const int HEX8_NFACEINTPTS=4;
00046 const Element::IntegPoint HEX8_INTPTS[]=
00047 {{ -0.577350,  -0.577350,  -0.577350,  1. }, 
00048 {   0.577350,  -0.577350,  -0.577350,  1. },
00049 {   0.577350,   0.577350,  -0.577350,  1. },
00050 {  -0.577350,   0.577350,  -0.577350,  1. },
00051 {  -0.577350,  -0.577350,   0.577350,  1. },
00052 {   0.577350,  -0.577350,   0.577350,  1. },
00053 {   0.577350,   0.577350,   0.577350,  1. },
00054 {  -0.577350,   0.577350,   0.577350,  1. }};
00055 const Element::IntegPoint HEX8_FACEINTPTS[]=
00056 {{ -0.577350,  -0.577350,  0.0,  1. }, 
00057 {   0.577350,  -0.577350,  0.0,  1. },
00058 {   0.577350,   0.577350,  0.0,  1. },
00059 {  -0.577350,   0.577350,  0.0,  1. }};
00060 // }}}
00061 
00062 class Hex8 : public virtual Element
00063 {
00064 public:
00065     // Constructor
00066     Hex8()
00067     {
00068         // Setup nodes number
00069                _n_nodes = HEX8_NNODES;
00070              _n_int_pts = HEX8_NINTPTS;
00071           _n_face_nodes = HEX8_NFACENODES;
00072         _n_face_int_pts = HEX8_NFACEINTPTS;
00073 
00074         // Allocate nodes (connectivity)
00075         _connects.resize(_n_nodes);
00076 
00077         // Setup pointer to the array of Integration Points
00078         _a_int_pts = HEX8_INTPTS;
00079 
00080         //_models.resize(_n_nodes);
00081     }
00082     // Destructor
00083     virtual ~Hex8() {}
00084     // Methods
00085     void         Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape)       const;
00086     void        Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs)      const;
00087     REAL BoundDistance(REAL r, REAL s, REAL t) const;
00088     int VTKCellType() const { return 12; } // VTK_HEXAHEDRON
00089 private:
00090 protected:
00091     void _distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects,
00092                                          REAL                 const   FaceValue   ,
00093                                          LinAlg::Vector<REAL>       & NodalValues ) const;
00094     void  _face_shape(REAL r, REAL s,         LinAlg::Vector<REAL> & AFaceShape)  const;
00095     void _face_derivs(REAL r, REAL s,         LinAlg::Matrix<REAL> & AFaceDerivs) const;
00096 }; // class Hex8
00097 
00098 
00100 
00101 
00102 inline void Hex8::_distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects, // {{{   // In:  Array of ptrs to face nodes
00103                                                   REAL                 const   FaceValue   ,          // In:  A value applied on a face to be converted to nodes
00104                                                   LinAlg::Vector<REAL>       & NodalValues ) const    // Out: The resultant nodal values
00105 {
00106     //Dimensioning NodalValues
00107     NodalValues.Resize(_n_face_nodes);
00108     NodalValues.SetValues(0.0);
00109     LinAlg::Matrix<REAL> J;                           // Jacobian matrix. size = 2 x 3
00110     LinAlg::Vector<REAL> face_shape(_n_face_nodes); // Shape functions of a face. size = _n_face_nodes
00111     //Integration along the face
00112     for (int i=0; i<_n_face_int_pts; i++)
00113     {
00114         REAL r = HEX8_FACEINTPTS[i].r;
00115         REAL s = HEX8_FACEINTPTS[i].s;
00116         REAL w = HEX8_FACEINTPTS[i].w;
00117         _face_shape(r, s, face_shape);
00118         _face_jacobian(FaceConnects, r, s, J);
00119         //REAL det_J = J.Det();
00120         //LinAlg::Axpy(FaceValue*det_J*w, face_shape, NodalValues); // NodalValues += FaceValue*face_shape*det_J*w
00121         NodalValues += FaceValue*face_shape*det(J)*w;
00122     }
00123 }
00124 // }}}
00125 
00126 inline void Hex8::_face_shape(REAL r, REAL s, LinAlg::Vector<REAL> & FaceShape) const // {{{
00127 {
00128     /*           s
00129      *           ^
00130      *           |             
00131      *         3 @-----------@ 2
00132      *           |           |
00133      *           |           |
00134      *           |           |
00135      *           |           |
00136      *           |           |
00137      *           @-----------@-> r
00138      *           0           1
00139      */
00140     FaceShape.Resize(4);
00141     FaceShape(0) = 0.25*(1.0-r-s+r*s);
00142     FaceShape(1) = 0.25*(1.0+r-s-r*s);
00143     FaceShape(2) = 0.25*(1.0+r+s+r*s);
00144     FaceShape(3) = 0.25*(1.0-r+s-r*s);
00145 } // }}}
00146 
00147 inline void Hex8::_face_derivs(REAL r, REAL s, LinAlg::Matrix<REAL> & FaceDerivs) const // {{{
00148 {
00149     /*           _     _ T
00150      *          |  dNi  |
00151      * Derivs = |  ---  |   , where cj = r, s
00152      *          |_ dcj _|
00153      *
00154      * Derivs(j,i), j=>local coordinate and i=>shape function
00155      */
00156     FaceDerivs.Resize(2,4);
00157     FaceDerivs(0,0) = 0.25*(-1.0+s);   FaceDerivs(1,0) = 0.25*(-1.0+r);
00158     FaceDerivs(0,1) = 0.25*(+1.0-s);   FaceDerivs(1,1) = 0.25*(-1.0-r);
00159     FaceDerivs(0,2) = 0.25*(+1.0+s);   FaceDerivs(1,2) = 0.25*(+1.0+r);
00160     FaceDerivs(0,3) = 0.25*(-1.0-s);   FaceDerivs(1,3) = 0.25*(+1.0-r);
00161 } // }}}
00162 
00163 inline void Hex8::Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape) const // {{{
00164 {
00165 
00166     /*                    t
00167      *                    ^
00168      *                    |     
00169      *                   4                7
00170      *                    @________________@
00171      *                  ,'|              ,'|
00172      *                ,'  |            ,'  |
00173      *              ,'    |          ,'    |
00174      *        5   ,'      |     6  ,'      |
00175      *          @'===============@'        |
00176      *          |         |      |         |
00177      *          |         |      |         |            
00178      *          |       0 @_____ | ________@ 3 --> s
00179      *          |       ,'       |       ,' 
00180      *          |     ,'         |     ,' 
00181      *          |   ,'           |   ,' 
00182      *          | ,'             | ,' 
00183      *          @________________@'
00184      *         1                  2 
00185      *      ,'
00186      *    |_
00187      *   r
00188      */
00189     Shape.Resize(8);
00190     Shape(0) = 0.125*(1.0-r-s+r*s-t+s*t+r*t-r*s*t);
00191     Shape(1) = 0.125*(1.0+r-s-r*s-t+s*t-r*t+r*s*t);
00192     Shape(2) = 0.125*(1.0+r+s+r*s-t-s*t-r*t-r*s*t);
00193     Shape(3) = 0.125*(1.0-r+s-r*s-t-s*t+r*t+r*s*t);
00194     Shape(4) = 0.125*(1.0-r-s+r*s+t-s*t-r*t+r*s*t);
00195     Shape(5) = 0.125*(1.0+r-s-r*s+t-s*t+r*t-r*s*t);
00196     Shape(6) = 0.125*(1.0+r+s+r*s+t+s*t+r*t+r*s*t);
00197     Shape(7) = 0.125*(1.0-r+s-r*s+t+s*t-r*t-r*s*t);
00198 } // }}}
00199 
00200 inline void Hex8::Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs) const // {{{
00201 {
00202     /*           _     _ T
00203      *          |  dNi  |
00204      * Derivs = |  ---  |   , where cj = r, s
00205      *          |_ dcj _|
00206      *
00207      * Derivs(j,i), j=>local coordinate and i=>shape function
00208      */
00209     Derivs.Resize(3,8);
00210     Derivs(0,0) = 0.125*(-1.0+s+t-s*t);   Derivs(1,0)=0.125*(-1.0+r+t-r*t);   Derivs(2,0)=0.125*(-1.0+r+s-r*s);
00211     Derivs(0,1) = 0.125*(+1.0-s-t+s*t);   Derivs(1,1)=0.125*(-1.0-r+t+r*t);   Derivs(2,1)=0.125*(-1.0-r+s+r*s);
00212     Derivs(0,2) = 0.125*(+1.0+s-t-s*t);   Derivs(1,2)=0.125*(+1.0+r-t-r*t);   Derivs(2,2)=0.125*(-1.0-r-s-r*s);
00213     Derivs(0,3) = 0.125*(-1.0-s+t+s*t);   Derivs(1,3)=0.125*(+1.0-r-t+r*t);   Derivs(2,3)=0.125*(-1.0+r-s+r*s);
00214     Derivs(0,4) = 0.125*(-1.0+s-t+s*t);   Derivs(1,4)=0.125*(-1.0+r-t+r*t);   Derivs(2,4)=0.125*(+1.0-r-s+r*s);
00215     Derivs(0,5) = 0.125*(+1.0-s+t-s*t);   Derivs(1,5)=0.125*(-1.0-r-t-r*t);   Derivs(2,5)=0.125*(+1.0+r-s-r*s);
00216     Derivs(0,6) = 0.125*(+1.0+s+t+s*t);   Derivs(1,6)=0.125*(+1.0+r+t+r*t);   Derivs(2,6)=0.125*(+1.0+r+s+r*s);
00217     Derivs(0,7) = 0.125*(-1.0-s-t-s*t);   Derivs(1,7)=0.125*(+1.0-r+t-r*t);   Derivs(2,7)=0.125*(+1.0-r+s-r*s);
00218 } // }}}
00219 
00220 inline REAL Hex8::BoundDistance(REAL r, REAL s, REAL t) const //{{{
00221 {
00222     return std::min(std::min( 1-fabs(r) , 1-fabs(s) ), 1-fabs(t)) ;
00223 } //}}}
00224 
00225 }; // namespace FEM
00226 
00227 #endif // MECHSYS_FEM_HEX8_H
00228 
00229 // vim:fdm=marker