MechSys: fem/ele/hex20.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_HEX20_H
00023 #define MECHSYS_FEM_HEX20_H
00024 
00025 #ifndef REAL
00026     #define REAL double
00027 #endif
00028 
00029 #include <algorithm> 
00030 #include "fem/ele/element.h"
00031 #include "linalg/vector.h"
00032 #include "linalg/matrix.h"
00033 #include "linalg/lawrap.h"
00034 
00035 
00036 namespace FEM
00037 {
00038 
00039 // Hex20 Constants // {{{
00040 const int HEX20_NNODES=20;
00041 const int HEX20_NINTPTS=8;
00042 const int HEX20_NFACENODES=8;
00043 const int HEX20_NFACEINTPTS=4;
00044 const Element::IntegPoint HEX20_INTPTS[]=
00045 {{ -0.577350,  -0.577350,  -0.577350,  1. }, 
00046 {   0.577350,  -0.577350,  -0.577350,  1. },
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 const Element::IntegPoint HEX20_FACEINTPTS[]=
00054 {{ -0.577350,  -0.577350,  0.0,  1. }, 
00055 {   0.577350,  -0.577350,  0.0,  1. },
00056 {   0.577350,   0.577350,  0.0,  1. },
00057 {  -0.577350,   0.577350,  0.0,  1. }};
00058 // }}}
00059 
00060 class Hex20 : public virtual Element
00061 {
00062 public:
00063     // Constructor
00064     Hex20()
00065     {
00066         // Setup nodes number
00067                _n_nodes = HEX20_NNODES;
00068              _n_int_pts = HEX20_NINTPTS;
00069           _n_face_nodes = HEX20_NFACENODES;
00070         _n_face_int_pts = HEX20_NFACEINTPTS;
00071 
00072         // Allocate nodes (connectivity)
00073         _connects.resize(_n_nodes);
00074 
00075         // Setup pointer to the array of Integration Points
00076         _a_int_pts = HEX20_INTPTS;
00077 
00078         //_models.resize(_n_nodes);
00079     }
00080     // Destructor
00081     virtual ~Hex20() {}
00082     void       Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape)       const;
00083     void      Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs)      const;
00084     REAL    BoundDistance(REAL r, REAL s, REAL t) const;
00085     // Methods
00086     int VTKCellType() const { return 25; } // VTK_QUADRATIC_HEXAHEDRON
00087 private:
00088 protected:
00089     void _distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects,
00090                                          REAL                 const   FaceValue   ,
00091                                          LinAlg::Vector<REAL>       & NodalValues ) const;
00092     void  _face_shape(REAL r, REAL s,         LinAlg::Vector<REAL> & AFaceShape)  const;
00093     void _face_derivs(REAL r, REAL s,         LinAlg::Matrix<REAL> & AFaceDerivs) const;
00094 }; // class Hex20
00095 
00096 
00098 
00099 
00100 inline void Hex20::_distrib_val_to_face_nodal_vals(Array<Node*>         const & FaceConnects, // {{{   // In:  Array of ptrs to face nodes
00101                                                   REAL                 const   FaceValue   ,          // In:  A value applied on a face to be converted to nodes
00102                                                   LinAlg::Vector<REAL>       & NodalValues ) const    // Out: The resultant nodal values
00103 {
00104     //Dimensioning NodalValues
00105     NodalValues.Resize(_n_face_nodes);
00106     NodalValues.SetValues(0.0);
00107     LinAlg::Matrix<REAL> J;                           // Jacobian matrix. size = 2 x 3
00108     LinAlg::Vector<REAL> face_shape(_n_face_nodes); // Shape functions of a face. size = _n_face_nodes
00109     //Integration along the face
00110     for (int i=0; i<_n_face_int_pts; i++)
00111     {
00112         REAL r = HEX20_FACEINTPTS[i].r;
00113         REAL s = HEX20_FACEINTPTS[i].s;
00114         REAL w = HEX20_FACEINTPTS[i].w;
00115         _face_shape(r, s, face_shape);
00116         _face_jacobian(FaceConnects, r, s, J);
00117         //REAL det_J = J.Det();
00118         NodalValues += FaceValue*face_shape*det(J)*w;
00119         //LinAlg::Axpy(FaceValue*det_J*w, face_shape, NodalValues); // NodalValues += FaceValue*face_shape*det_J*w
00120     }
00121 }
00122 // }}}
00123 
00124 inline void Hex20::_face_shape(REAL r, REAL s, LinAlg::Vector<REAL> & FaceShape) const // {{{
00125 {
00126     /*                 ^ s
00127      *                 |
00128      *                 5       
00129      *         6 @-----------@ 4
00130      *           |           |
00131      *           |           |
00132      *         7 |           | 3  -> r
00133      *           |           |
00134      *           |           |
00135      *           @-----------@
00136      *           0     1     2
00137      */
00138     FaceShape.Resize(8);
00139     // revisar estas funçoes de forma
00140     FaceShape(0) = 0.25 * (1-r) * (1-s) * (-1-r-s);    
00141     FaceShape(1) =  0.5 * (1-r) * (1+r) * (1-s);
00142     FaceShape(2) = 0.25 * (1+r) * (1-s) * (-1+r-s);
00143     FaceShape(3) =  0.5 * (1+r) * (1+s) * (1-s);      
00144     FaceShape(4) = 0.25 * (1+r) * (1+s) * (-1+r+s);  
00145     FaceShape(5) =  0.5 * (1-r) * (1+r) * (1+s);
00146     FaceShape(6) = 0.25 * (1-r) * (1+s) * (-1-r+s);
00147     FaceShape(7) =  0.5 * (1-r) * (1-s) * (1+s);      
00148 } // }}}
00149 
00150 inline void Hex20::_face_derivs(REAL r, REAL s, LinAlg::Matrix<REAL> & FaceDerivs) const // {{{
00151 {
00152     /*           _     _ T
00153      *          |  dNi  |
00154      * Derivs = |  ---  |   , where cj = r, s
00155      *          |_ dcj _|
00156      *
00157      * Derivs(j,i), j=>local coordinate and i=>shape function
00158      */
00159     FaceDerivs.Resize(2,8);
00160 
00161     FaceDerivs(0,0) = 0.25 * (s + 2*r - 2*r*s - s*s);
00162     FaceDerivs(0,1) = -r + r*s;
00163     FaceDerivs(0,2) = 0.25 * (-s + 2*r - 2*r*s + s*s);
00164     FaceDerivs(0,3) = 0.5 * (1 - s*s);
00165     FaceDerivs(0,4) = 0.25 * (s + 2*r + 2*r*s + s*s);
00166     FaceDerivs(0,5) = -r - r*s;
00167     FaceDerivs(0,6) = 0.25 * (-s + 2*r + r*s - s*s);
00168     FaceDerivs(0,7) = 0.5 * (-1 + s*s);
00169     
00170     FaceDerivs(1,0) = 0.25 * (r + 2*s - r*r - 2*r*s);
00171     FaceDerivs(1,1) = 0.5  * (r*r - 1);
00172     FaceDerivs(1,2) = 0.25 * (-r + 2*s - r*r + 2*r*s);
00173     FaceDerivs(1,3) = -s - r*s;
00174     FaceDerivs(1,4) = 0.25 * (r + 2*s + r*r + 2*r*s);
00175     FaceDerivs(1,5) = 0.5  * (-r*r + 1);
00176     FaceDerivs(1,6) = 0.25 * (-r + 2*s + r*r - 2*r*s);
00177     FaceDerivs(1,7) = -s + r*s;
00178 
00179 } // }}}
00180 
00181 inline void  Hex20::Shape(REAL r, REAL s, REAL t, LinAlg::Vector<REAL> & Shape) const // {{{
00182 {
00183 
00184     /*                    t
00185      *                    ^
00186      *                    |     
00187      *                  12       19       18
00188      *                    @________________@
00189      *                  ,'|              ,'|
00190      *             13 ,'  |         17 ,'  |
00191      *              ,'    |          ,'    |
00192      *        14  ,'      | 8   16 ,'      |11
00193      *          @'===============@'        |
00194      *          |      15 |      |         |
00195      *          |         |      |7        |            
00196      *          |       0 @_____ | ________@ 6 --> s
00197      *        9 |       ,'       | 10    ,' 
00198      *          |   1 ,'         |     ,' 
00199      *          |   ,'           |   ,' 5
00200      *          | ,'             | ,' 
00201      *          @________________@'
00202      *         2        3         4 
00203      *      ,'
00204      *    |_
00205      *   r
00206      */
00207     Shape.Resize(20);
00208 
00209 // coords 1..8
00210     Shape( 0) = 1/8.*(1-r)  *(1-s)  *(1-t)  *(-r-s-t-2);
00211     Shape( 1) = 1/4.*(1-r*r)*(1-s)  *(1-t)             ;
00212     Shape( 2) = 1/8.*(1+r)  *(1-s)  *(1-t)  *( r-s-t-2);
00213     Shape( 3) = 1/4.*(1+r)  *(1-s*s)*(1-t)             ;
00214     Shape( 4) = 1/8.*(1+r)  *(1+s)  *(1-t)  *( r+s-t-2);
00215     Shape( 5) = 1/4.*(1-r*r)*(1+s)  *(1-t)             ;
00216     Shape( 6) = 1/8.*(1-r)  *(1+s)  *(1-t)  *(-r+s-t-2);
00217     Shape( 7) = 1/4.*(1-r)  *(1-s*s)*(1-t)             ;
00218 // coords 9..12
00219     Shape( 8) = 1/4.*(1-r)  *(1-s)  *(1-t*t)           ;
00220     Shape( 9) = 1/4.*(1+r)  *(1-s)  *(1-t*t)           ;
00221     Shape(10) = 1/4.*(1+r)  *(1+s)  *(1-t*t)           ;
00222     Shape(11) = 1/4.*(1-r)  *(1+s)  *(1-t*t)           ;
00223 // coords 13..20
00224     Shape(12) = 1/8.*(1-r)  *(1-s)  *(1+t)  *(-r-s+t-2);
00225     Shape(13) = 1/4.*(1-r*r)*(1-s)  *(1+t)             ;
00226     Shape(14) = 1/8.*(1+r)  *(1-s)  *(1+t)  *( r-s+t-2);
00227     Shape(15) = 1/4.*(1+r)  *(1-s*s)*(1+t)             ;
00228     Shape(16) = 1/8.*(1+r)  *(1+s)  *(1+t)  *( r+s+t-2);
00229     Shape(17) = 1/4.*(1-r*r)*(1+s)  *(1+t)             ;
00230     Shape(18) = 1/8.*(1-r)  *(1+s)  *(1+t)  *(-r+s+t-2);
00231     Shape(19) = 1/4.*(1-r)  *(1-s*s)*(1+t)             ;
00232     
00233 } // }}}
00234 
00235 inline void Hex20::Derivs(REAL r, REAL s, REAL t, LinAlg::Matrix<REAL> & Derivs) const // {{{
00236 {
00237     /*           _     _ T
00238      *          |  dNi  |
00239      * Derivs = |  ---  |   , where cj = r, s
00240      *          |_ dcj _|
00241      *
00242      * Derivs(j,i), j=>local coordinate and i=>shape function
00243      */
00244 
00245     Derivs.Resize(3,20);
00246     
00247     //Drrivatrs with rrsprtt to r
00248     Derivs(0, 0)= -.125  *(1-s)  *(1-t)  *(-r-s-t-2)-0.125*(1-r)*(1-s)*(1-t);
00249     Derivs(0, 1)= -.5  *r*(1-s)  *(1-t);
00250     Derivs(0, 2)=  .125  *(1-s)  *(1-t)  *( r-s-t-2)+0.125*(1+r)*(1-s)*(1-t);
00251     Derivs(0, 3)=  .25   *(1-s*s)*(1-t);
00252     Derivs(0, 4)=  .125  *(1+s)  *(1-t)  *( r+s-t-2)+0.125*(1+r)*(1+s)*(1-t);
00253     Derivs(0, 5)= -.5  *r*(1+s)  *(1-t);
00254     Derivs(0, 6)= -.125  *(1+s)  *(1-t)  *(-r+s-t-2)-0.125*(1-r)*(1+s)*(1-t);
00255     Derivs(0, 7)= -.25   *(1-s*s)*(1-t);
00256 
00257     Derivs(0, 8)= -.25   *(1-s)  *(1-t*t);
00258     Derivs(0, 9)=  .25   *(1-s)  *(1-t*t);
00259     Derivs(0,10)=  .25   *(1+s)  *(1-t*t);
00260     Derivs(0,11)= -.25   *(1+s)  *(1-t*t);
00261 
00262     Derivs(0,12)= -.125  *(1-s)  *(1+t)  *(-r-s+t-2)-0.125*(1-r)*(1-s)*(1+t);
00263     Derivs(0,13)= -.5  *r*(1-s)  *(1+t);
00264     Derivs(0,14)=  .125  *(1-s)  *(1+t)  *( r-s+t-2)+0.125*(1+r)*(1-s)*(1+t);
00265     Derivs(0,15)=  .25   *(1-s*s)*(1+t);
00266     Derivs(0,16)=  .125  *(1+s)  *(1+t)  *( r+s+t-2)+0.125*(1+r)*(1+s)*(1+t);
00267     Derivs(0,17)= -.5  *r*(1+s)  *(1+t);
00268     Derivs(0,18)= -.125  *(1+s)  *(1+t)  *(-r+s+t-2)-0.125*(1-r)*(1+s)*(1+t);
00269     Derivs(0,19)= -.25   *(1-s*s)*(1+t);
00270 
00271     //Drrivatrs with rrsprtt to s
00272     Derivs(1, 0)= -.125  *(1-r)  *(1-t)  *(-r-s-t-2)-0.125*(1-r)*(1-s)*(1-t);
00273     Derivs(1, 1)= -.25   *(1-r*r)*(1-t)  ;
00274     Derivs(1, 2)= -.125  *(1+r)  *(1-t)  *( r-s-t-2)-0.125*(1+r)*(1-s)*(1-t);
00275     Derivs(1, 3)= -.5  *s*(1+r)  *(1-t)  ;
00276     Derivs(1, 4)=  .125  *(1+r)  *(1-t)  *( r+s-t-2)+0.125*(1+r)*(1+s)*(1-t);
00277     Derivs(1, 5)=  .25   *(1-r*r)*(1-t)  ;
00278     Derivs(1, 6)=  .125  *(1-r)  *(1-t)  *(-r+s-t-2)+0.125*(1-r)*(1+s)*(1-t);
00279     Derivs(1, 7)= -.5  *s*(1-r)  *(1-t)  ;
00280 
00281     Derivs(1, 8)= -.25   *(1-r)  *(1-t*t);
00282     Derivs(1, 9)= -.25   *(1+r)  *(1-t*t);
00283     Derivs(1,10)=  .25   *(1+r)  *(1-t*t);
00284     Derivs(1,11)=  .25   *(1-r)  *(1-t*t);
00285 
00286     Derivs(1,12)= -.125  *(1-r)  *(1+t)  *(-r-s+t-2)-0.125*(1-r)*(1-s)*(1+t);
00287     Derivs(1,13)= -.25   *(1-r*r)*(1+t)  ;
00288     Derivs(1,14)= -.125  *(1+r)  *(1+t)  *( r-s+t-2)-0.125*(1+r)*(1-s)*(1+t);
00289     Derivs(1,15)= -.5  *s*(1+r)  *(1+t)  ;
00290     Derivs(1,16)=  .125  *(1+r)  *(1+t)  *( r+s+t-2)+0.125*(1+r)*(1+s)*(1+t);
00291     Derivs(1,17)=  .25   *(1-r*r)*(1+t)  ;
00292     Derivs(1,18)=  .125  *(1-r)  *(1+t)  *(-r+s+t-2)+0.125*(1-r)*(1+s)*(1+t);
00293     Derivs(1,19)= -.5  *s*(1-r)  *(1+t)  ;
00294 
00295     //Drrivatrs with rrsprtt to t
00296     Derivs(2, 0)= -.125  *(1-r)  *(1-s)  *(-r-s-t-2)-0.125*(1-r)*(1-s)*(1-t);
00297     Derivs(2, 1)= -.25   *(1-r*r)*(1-s)  ;
00298     Derivs(2, 2)= -.125  *(1+r)  *(1-s)  *( r-s-t-2)-0.125*(1+r)*(1-s)*(1-t);
00299     Derivs(2, 3)= -.25   *(1+r)  *(1-s*s);
00300     Derivs(2, 4)= -.125  *(1+r)  *(1+s)  *( r+s-t-2)-0.125*(1+r)*(1+s)*(1-t);
00301     Derivs(2, 5)= -.25   *(1-r*r)*(1+s)  ;
00302     Derivs(2, 6)= -.125  *(1-r)  *(1+s)  *(-r+s-t-2)-0.125*(1-r)*(1+s)*(1-t);
00303     Derivs(2, 7)= -.25   *(1-r)  *(1-s*s);
00304 
00305     Derivs(2, 8)= -.5  *t*(1-r)  *(1-s)  ;
00306     Derivs(2, 9)= -.5  *t*(1+r)  *(1-s)  ;
00307     Derivs(2,10)= -.5  *t*(1+r)  *(1+s)  ;
00308     Derivs(2,11)= -.5  *t*(1-r)  *(1+s)  ;
00309 
00310     Derivs(2,12)=  .125  *(1-r)  *(1-s)  *(-r-s+t-2)+0.125*(1-r)*(1-s)*(1+t);
00311     Derivs(2,13)=  .25   *(1-r*r)*(1-s)  ;
00312     Derivs(2,14)=  .125  *(1+r)  *(1-s)  *( r-s+t-2)+0.125*(1+r)*(1-s)*(1+t);
00313     Derivs(2,15)=  .25   *(1+r)  *(1-s*s);
00314     Derivs(2,16)=  .125  *(1+r)  *(1+s)  *( r+s+t-2)+0.125*(1+r)*(1+s)*(1+t);
00315     Derivs(2,17)=  .25   *(1-r*r)*(1+s)  ;
00316     Derivs(2,18)=  .125  *(1-r)  *(1+s)  *(-r+s+t-2)+0.125*(1-r)*(1+s)*(1+t);
00317     Derivs(2,19)=  .25   *(1-r)  *(1-s*s);
00318 
00319 } // }}}
00320     
00321 inline REAL Hex20::BoundDistance(REAL r, REAL s, REAL t) const //{{{
00322 {
00323     return std::min(std::min( 1-fabs(r) , 1-fabs(s) ), 1-fabs(t)) ;
00324 } //}}}
00325 
00326 }; // namespace FEM
00327 
00328 #endif // MECHSYS_FEM_HEX20_H
00329 
00330 // vim:fdm=marker