• Home
• About DUNE
  • What is DUNE
  • News
  • Gallery
  • Publications
• Download
  • Core Modules (Stable)
    • Release Notes (v1.1)
    • Release Notes (v1.0)
  • Recent Changes
  • External Modules
  • Dune Logo
• Documentation
  • Installation
  • Using DUNE
  • Class Documentation
• Community
  • People
  • Mailinglists
  • FAQ
• Development
  • Developing DUNE
  • Unstable Modules
    • Download via SVN
    • Class Documentation
  • Bug Tracking System
  • Automated Test Environment

linearelasticityassembler.hh

Go to the documentation of this file.

#ifndef DUNE_P1_LINEAR_ELASTICITY_ASSEMBLER_HH
#define DUNE_P1_LINEAR_ELASTICITY_ASSEMBLER_HH

#include<map>
#include<iostream>
#include<iomanip>
#include<fstream>
#include<vector>
#include<sstream>

#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/common/exceptions.hh>
#include<dune/grid/common/grid.hh>
#include<dune/grid/common/referenceelements.hh>
#include<dune/istl/operators.hh>
#include<dune/istl/bvector.hh>

#include<dune/grid/common/quadraturerules.hh>
#include<dune/disc/operators/localstiffness.hh>
#include<dune/disc/shapefunctions/lagrangeshapefunctions.hh>
#include<dune/disc/operators/boundaryconditions.hh>

namespace Dune
{


  template<class GridView, class RT>
  class LinearElasticityLocalStiffness 
    : public LinearLocalStiffness<GridView,RT,GridView::dimension>
  {
        // grid types
      typedef typename GridView::Grid::ctype DT;
        typedef typename GridView::template Codim<0>::Entity Entity;

        // some other sizes
        enum {dim=GridView::dimension};

      typedef FieldVector<double, (dim+1)*dim/2> SymmTensor;

  public:
        // define the number of components of your system, this is used outside
        // to allocate the correct size of (dense) blocks with a FieldMatrix
        enum {m=dim};

        // types for matrics, vectors and boundary conditions
        typedef FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
        typedef FieldVector<RT,m> VBlockType;   // one entry in the global vectors
        typedef array<BoundaryConditions::Flags,m> BCBlockType;     // componentwise boundary conditions
        // /////////////////////////////////
        //   The material parameters
        // /////////////////////////////////
        
        // Young's modulus
        double E_;

        // Poisson ratio
        double nu_;

      typedef std::map<std::pair<GeometryType, int>, std::vector<FieldVector<DT,dim> > > CacheType;
      CacheType gradientCache_;
      
      LinearElasticityLocalStiffness ()
      {}

      LinearElasticityLocalStiffness (double E, double nu, bool procBoundaryAsDirichlet_=true)
          : E_(E), nu_(nu)
      {
          procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
      }

      void setEandNu(double E, double nu)
      {
          E_  = E;
          nu_ = nu;
      }
      

      void assemble (const Entity& e, int k=1)
      {
          // extract some important parameters
          GeometryType gt = e.type();
          const typename LagrangeShapeFunctionSetContainer<DT,RT,dim>::value_type& sfs
              = LagrangeShapeFunctions<DT,RT,dim>::general(gt,k);
          
          // clear assemble data
          setcurrentsize(sfs.size());

          for (int i=0; i<sfs.size(); i++)
              {
                  this->b[i] = 0;
                  for (int j=0; j<this->bctype[i].size(); j++)
                      this->bctype[i][j] = BoundaryConditions::neumann;
                  for (int j=0; j<sfs.size(); j++) 
                      this->A[i][j] = 0;
              }
          
          // Compute suitable quadrature order
          int p = (gt.isSimplex()) 
              ? 2*(k-1)
              : 2*(k*dim-1);

          const std::vector<FieldVector<DT,dim> >* shapeGradients = NULL;

          typename CacheType::iterator shapeFunctionGradients = gradientCache_.find(std::make_pair(gt,p));

          if (shapeFunctionGradients == gradientCache_.end()) {

              std::vector<FieldVector<DT,dim> > newGradientSet(sfs.size() * QuadratureRules<DT,dim>::rule(gt,p).size());
              
              for (size_t g=0; g<Dune::QuadratureRules<DT,dim>::rule(gt,p).size(); ++g) {

                  // pos of integration point
                  const Dune::FieldVector<DT,dim>& quadPos = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].position();

                  for (int i=0; i<sfs.size(); i++) {
                  
                      for (int l=0; l<dim; l++) 
                          newGradientSet[sfs.size()*g + i][l] = sfs[i].evaluateDerivative(0,l,quadPos);
                      
                  }

              }

              gradientCache_.insert(std::make_pair(std::make_pair(gt,p), newGradientSet));

          }

          shapeGradients = &gradientCache_.find(std::make_pair(gt,p))->second;


          for (size_t g=0; g<Dune::QuadratureRules<DT,dim>::rule(gt,p).size(); ++g) {
              
              // pos of integration point
              const Dune::FieldVector<DT,dim>& local = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].position();
              
              // eval jacobian inverse
              const Dune::FieldMatrix<DT,dim,dim>& jac = e.geometry().jacobianInverseTransposed(local);

              // weight of quadrature point
              double weight = Dune::QuadratureRules<DT,dim>::rule(gt,p)[g].weight();
              
              // determinant of jacobian
              DT detjac = e.geometry().integrationElement(local);

              // evaluate gradients at Gauss points
              Dune::FieldVector<DT,dim> grad[sfs.size()], temp, gv;
              for (int i=0; i<sfs.size(); i++) {
                  grad[i] = 0;
                  jac.umv((*shapeGradients)[sfs.size()*g + i],grad[i]); // transform gradient to global coordinates
              }
              
              // /////////////////////////////////////////////
              //   Compute strain for all shape functions
              // /////////////////////////////////////////////
              std::vector<array<SymmTensor,dim> > strain(sfs.size());
              
              for (int i=0; i<sfs.size(); i++) 
                  for (int k=0; k<dim; k++) {
                      
                      // The deformation gradient of the shape function
                      FieldMatrix<double, dim, dim> deformationGradient(0);
                      deformationGradient[k] = grad[i];
                 
                      /* Computes the linear strain tensor from the deformation gradient*/
                      computeStrain(deformationGradient,strain[i][k]);
                      
                  }
              
              // loop over test function number
              for (int row=0; row<sfs.size(); row++) {
                  
                  for (int rcomp=0; rcomp<dim; rcomp++) {
                      
                      SymmTensor stress = hookeTimesStrain(strain[row][rcomp]);
                      
                      for (int col=0; col<=row; col++) {
                          
                          for (int ccomp=0; ccomp<dim; ccomp++) {
                              
                              this->A[row][col][rcomp][ccomp] += stress*strain[col][ccomp] * weight * detjac;

                          }
                      }
                  }    

              }
              
          }

          for (int row=0; row<sfs.size(); row++)
              for (int col=0; col<=row; col++) {
                  
                  // Complete the symmetric matrix by copying the lower left triangular matrix
                  // to the upper right one
                  if (row!=col) {
                      for (int rcomp=0; rcomp<dim; rcomp++)
                          for (int ccomp=0; ccomp<dim; ccomp++)
                              this->A[col][row][ccomp][rcomp] = this->A[row][col][rcomp][ccomp];
                  }
              }

      }

      // computes the linear strain from the deformation gradient
      void computeStrain(const FieldMatrix<double, dim, dim>& grad, 
                         SymmTensor& strain) const
      {
          
          if (dim==2) {
              
              strain[0] = grad[0][0];
              strain[1] = grad[1][1];
              strain[2] = grad[0][1] + grad[1][0];
              
          } else if (dim==3) {
              
              strain[0] = grad[0][0];
              strain[1] = grad[1][1];
              strain[2] = grad[2][2];
              strain[3] = grad[0][1] + grad[1][0];
              strain[4] = grad[0][2] + grad[2][0];
              strain[5] = grad[2][1] + grad[1][2];
              
          } else
                DUNE_THROW(NotImplemented, "No elasticity assembler for " << dim << "-dimensional problems");
      }
      
      
      SymmTensor hookeTimesStrain(const SymmTensor& strain) const {
          
          if (dim==3) {
              
              // compute Hooke Tensor
              FieldMatrix<double, 6, 6> hookeTensor;
              
              hookeTensor = 0;
              
              hookeTensor[0][0] = 1-nu_;
              hookeTensor[0][1] = nu_;
              hookeTensor[0][2] = nu_;
              
              hookeTensor[1][0] = nu_;
              hookeTensor[1][1] = 1-nu_;
              hookeTensor[1][2] = nu_;
              
              hookeTensor[2][0] = nu_;
              hookeTensor[2][1] = nu_;
              hookeTensor[2][2] = 1-nu_;
              
              hookeTensor[3][3] = 0.5 -nu_;
              hookeTensor[4][4] = 0.5 -nu_;
              hookeTensor[5][5] = 0.5 -nu_;
              
              hookeTensor *= (E_/(1+nu_)/(1-2*nu_));
              
              // compute stress
              SymmTensor stress;
              stress = 0;
              hookeTensor.umv(strain, stress);
              
              return stress;
              
          } else if (dim==2) {
              
              // compute Hooke Tensor
              FieldMatrix<double, 3, 3> hookeTensor;
              
              hookeTensor = 0;
              
              hookeTensor[0][0] = 1-nu_;
              hookeTensor[0][1] = nu_;
              
              hookeTensor[1][0] = nu_;
              hookeTensor[1][1] = 1-nu_;
              
              hookeTensor[2][2] = 0.5 -nu_;
              
              hookeTensor *= (E_/(1+nu_)/(1-2*nu_));
              
              // compute stress
              SymmTensor stress;
              stress = 0;
              hookeTensor.umv(strain, stress);
              
              return stress;
              
          } else
              DUNE_THROW(NotImplemented, "No elasticity assembler for " << dim << "-dimensional problems");
          
      }

        void assembleBoundaryCondition (const Entity& e, int k=1)
        {
        }

  private:

        // parameters given in constructor
        bool procBoundaryAsDirichlet;

  };

}
#endif

Generated on 6 Jan 2009 with Doxygen (ver 1.5.1) [logfile].