Commit 40a6b7a0 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

obsolete

[[Imported from SVN: r10415]]
parent 31cd1a80
#ifndef VISCOELASTICITY_ASSEMBLER_HH
#define VISCOELASTICITY_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>
/**
* @file
* @brief Compute local stiffness matrix for conforming finite elements for viscous part of viscoelasticity equation
*/
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 viscous equation stress = N * strainrate N viscosity tensor
- 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
\tparam GridView a DUNE gridview type
\tparam RT type used for return values
*/
template<class GridView, class RT>
class LinearViscosityLocalStiffness
: 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};
//! The engineers' way of writing a symmetric second-order tensor
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
// /////////////////////////////////
// Shear Viscosity
double nu_shear_;
// Bulk Viscosity
double nu_bulk_;
//! Default Constructor
LinearViscosityLocalStiffness ()
{}
//! Constructor
LinearViscosityLocalStiffness (double nu_shear, double nu_bulk, bool procBoundaryAsDirichlet_=true)
: nu_shear_(nu_shear), nu_bulk_(nu_bulk)
{}
/** \brief Set material parameters */
void setNu_shearandNu_bulk(double nu_shear, double nu_bulk)
{
nu_shear_ = nu_shear;
nu_bulk_ = nu_bulk;
}
//! 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
*/
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);
setcurrentsize(sfs.size());
this->A = 0;
// clear assemble data
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;
}
// Compute suitable quadrature order
int p = (gt.isSimplex())
? 2*(k-1)
: 2*(k*dim-1);
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++) {
for (int l=0; l<dim; l++)
temp[l] = sfs[i].evaluateDerivative(0,l,local);
grad[i] = 0;
jac.umv(temp,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 viscous_stress = viscoustensorTimesStrain(strain[row][rcomp]);
for (int col=0; col<=row; col++) {
for (int ccomp=0; ccomp<dim; ccomp++) {
// test this->A[col][row][ccomp][rcomp] += viscous_stress*strain[col][ccomp] * weight * detjac;
this->A[row][col][rcomp][ccomp] += viscous_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];
}
}
}
protected:
// 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 viscosity assembler for " << dim << "-dimensional problems");
}
SymmTensor viscoustensorTimesStrain(const SymmTensor& strain) const {
double lambda = nu_bulk_-(2/3)*nu_shear_; //second viscosityparameter
if (dim==3) {
// compute ViscosityTensor
FieldMatrix<double, 6, 6> viscousTensor;
viscousTensor = 0;
viscousTensor[0][0] = 2*nu_shear_ + lambda;
viscousTensor[0][1] = lambda;
viscousTensor[0][2] = lambda;
viscousTensor[1][0] = lambda;
viscousTensor[1][1] = 2*nu_shear_ + lambda;
viscousTensor[1][2] = lambda;
viscousTensor[2][0] = lambda;
viscousTensor[2][1] = lambda;
viscousTensor[2][2] = 2*nu_shear_ + lambda;
viscousTensor[3][3] = 2*nu_shear_;
viscousTensor[4][4] = 2*nu_shear_;
viscousTensor[5][5] = 2*nu_shear_;
// compute viscous_stress
SymmTensor viscous_stress;
viscous_stress = 0;
viscousTensor.umv(strain, viscous_stress);
return viscous_stress;
} else if (dim==2) {
// compute ViscosityTensor
FieldMatrix<double, 3, 3> viscousTensor;
viscousTensor = 0;
viscousTensor[0][0] = 2*nu_shear_ + lambda;
viscousTensor[0][1] = lambda;
viscousTensor[1][0] = lambda;
viscousTensor[1][1] = 2*nu_shear_ + lambda;
viscousTensor[2][2] = 2*nu_shear_;
// compute viscous_stress
SymmTensor viscous_stress;
viscous_stress = 0;
viscousTensor.umv(strain, viscous_stress);
return viscous_stress;
} else
DUNE_THROW(NotImplemented, "No viscosity assembler for " << dim << "-dimensional problems");
}
//! should allow to assmble boundary conditions only
void assembleBoundaryCondition (const Entity& e, int k=1)
{
}
};
/** @} */
}
#endif
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