Commit 1ac3ef4f authored by Oliver Sander's avatar Oliver Sander Committed by sander@FU-BERLIN.DE
Browse files

file now in dune module 'cellelast'

[[Imported from SVN: r10487]]
parent 1eac3d6e
#include <dune/localfunctions/lagrange/pqkfactory.hh>
* @file
* @brief Compute local stiffness matrix for conforming finite elements for linear elasticity equation
* @author Oliver Sander
namespace Dune
/** @addtogroup DISC
* @{
* @brief compute local stiffness matrix for the linear elasticity operator
//! A class for computing local stiffness matrices
/*! A class for computing local stiffness matrix for the
linear elasticity equation
- div \sigma = f in Omega
u = g on Gamma1; j*n = J on Gamma2.
Uses conforming finite elements with the Lagrange shape functions.
It should work for all dimensions and element types.
All the numbering is with respect to the reference element and the
Lagrange shape functions
Template parameters are:
- GridType a DUNE grid type
- RT type used for return values
template<class GridType, class RT>
class UniformLinearElasticityLocalStiffness
: public LocalStiffness<UniformLinearElasticityLocalStiffness<GridType,RT>,GridType,RT,GridType::dimension>
// grid types
typedef typename GridType::ctype DT;
typedef typename GridType::Traits::template Codim<0>::Entity Entity;
typedef typename GridType::template Codim<0>::EntityPointer EEntityPointer;
// some other sizes
enum {dim=GridType::dimension};
//! The engineers' way of writing a symmetric second-order tensor
typedef FieldVector<double, (dim+1)*dim/2> SymmTensor;
// 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 FixedArray<BoundaryConditions::Flags,m> BCBlockType; // componentwise boundary conditions
enum {SIZE = LocalStiffness<UniformLinearElasticityLocalStiffness,GridType,RT,m>::SIZE};
// /////////////////////////////////
// The material parameters
// /////////////////////////////////
// Young's modulus
double E_;
// Poisson ratio
double nu_;
//! Constructor
UniformLinearElasticityLocalStiffness (double E, double nu, bool procBoundaryAsDirichlet_=true)
: E_(E), nu_(nu), isComputed_(false)
this->currentsize_ = 0;
procBoundaryAsDirichlet = procBoundaryAsDirichlet_;
// For the time being: all boundary conditions are homogeneous Neumann
// This means no boundary condition handling is done at all
for (int i=0; i<SIZE; i++)
this->bctype[i] = BoundaryConditions::neumann;
/** \brief Set material parameters */
void setEandNu(double E, double nu)
E_ = E;
nu_ = nu;
//! assemble local stiffness matrix for given element and order
/*! On exit the following things have been done:
- The stiffness matrix for the given entity and polynomial degree has been assembled and is
accessible with the mat() method.
- The boundary conditions have been evaluated and are accessible with the bc() method
- The right hand side has been assembled. It contains either the value of the essential boundary
condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
@param[in] e a codim 0 entity reference
@param[in] k order of Lagrange basis
template<typename TypeTag>
void assemble (const Entity& e, int k=1)
if (!isComputed_) {
// extract some important parameters
GeometryType gt = e.geometry().type();
//chache of local finite elements
typedef Dune::PQkLocalFiniteElementCache<typename GridType::ctype, RT, GridType::dimension, k> FiniteElementCache;
FiniteElementCache0 cache;
const typename FiniteElementCache::FiniteElementType& finiteElement = cache.get(gt);
this->currentsize_ = finiteElement.localBasis();
// clear assemble data
for (int i=0; i<finiteElement.localBasis(); i++)
this->b[i] = 0;
this->bctype[i][0] = BoundaryConditions::neumann;
for (int j=0; j<finiteElement.localBasis(); j++)
privateA[i][j] = 0;
// Compute suitable quadrature order
int p = (gt.isSimplex())
? 2*(k-1)
: 2*k*dim;
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
std::vector<Dune::FieldMatrix<dim,1> > referenceGradients(finiteElement.localBasis().size());
std::vector<Dune::FieldVector<double,dim> > grad(finiteElement.localBasis().size());
for (int i=0; i<finiteElement.localBasis.size(); i++) {
grad[i] = 0;
jac.umv(referenceGradients[i][0],grad[i]); // transform gradient to global ooordinates
// /////////////////////////////////////////////
// Compute strain for all shape functions
// /////////////////////////////////////////////
std::vector<FixedArray<SymmTensor,dim> > strain(finiteElement.localBasis().size());
for (int i=0; i<finiteElement.localBasis().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*/
// loop over test function number
for (int row=0; row<finiteElement.localBasis().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++)
privateA[row][col][rcomp][ccomp] += stress*strain[col][ccomp] * weight * detjac;
isComputed_ = true;
// /////////////////////////////////////////////////////////////////////////////
// Hand over the locally computed (or kept) matrix to the global assembler
// /////////////////////////////////////////////////////////////////////////////
for (int row=0; row<this->currentsize_; row++) {
for (int col=0; col<=row; col++) {
// Complete the symmetric matrix by copying the lower left triangular matrix
// to the upper right one
for (int rcomp=0; rcomp<dim; rcomp++)
for (int ccomp=0; ccomp<dim; ccomp++) {
this->A[row][col][rcomp][ccomp] = privateA[row][col][rcomp][ccomp];
this->A[col][row][ccomp][rcomp] = privateA[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");
//! should allow to assmble boundary conditions only
template<typename Tag>
void assembleBoundaryCondition (const Entity& e, int k=1)
/** If this variable is false, then the local element matrix is computed and stored.
Then variable is then set to true. If it is true, then the kept value is simply
handed over to the calling assembler. That saves lots of computing time when
we know that all element matrices are identical. */
bool isComputed_;
MBlockType privateA[SIZE][SIZE];
// parameters given in constructor
bool procBoundaryAsDirichlet;
/** @} */
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment