Commit 048f7106 authored by akbib's avatar akbib Committed by akbib@FU-BERLIN.DE
Browse files

use dune-localfunctions instead of dune-disc shapefunctions and use the symmetrictensor class

[[Imported from SVN: r10475]]
parent 155d0c5f
......@@ -5,8 +5,11 @@
#include <dune/istl/bvector.hh>
#include <dune/istl/bdmatrix.hh>
#include <dune/istl/solvers.hh>
#include <dune/disc/miscoperators/massmatrix.hh>
#include <dune/disc/operators/p1operator.hh>
#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
#include <dune/ag-common/assemblers/operatorassembler.hh>
#include <dune/ag-common/assemblers/localassemblers/massassembler.hh>
#include <dune/ag-common/symmetrictensor.hh>
#include <dune/ag-common/lumpmatrix.hh>
//#define LUMP_MATRIX
......@@ -18,9 +21,6 @@ class ZienkiewiczZhuEstimator
enum {dim = GridType::dimension};
//! The engineers' way of writing a symmetric second-order tensor
typedef Dune::FieldVector<double, (dim+1)*dim/2> SymmTensor;
public:
ZienkiewiczZhuEstimator(const GridType& grid, double E, double nu)
......@@ -41,17 +41,23 @@ public:
// Create mass matrix
// ////////////////////////////////////////////////////////////////////////
typedef BCRSMatrix<FieldMatrix<double,SymmTensor::size, SymmTensor::size> > MassMatrixType;
MassMatrixLocalStiffness<typename GridType::LeafGridView,double,SymmTensor::size> localMassMatrix;
LeafP1OperatorAssembler<GridType,double,SymmTensor::size> massMatrix(*grid_);
LeafP1Function<GridType,double,SymmTensor::size> u0(*grid_),f0(*grid_);
massMatrix.assemble(localMassMatrix,u0,f0);
typedef FieldMatrix<double,SymmetricTensor<dim>::size, SymmetricTensor<dim>::size> BlockType;
typedef BCRSMatrix<BlockType> MassMatrixType;
typedef P1NodalBasis<typename GridType::LeafGridView,double> P1Basis;
P1Basis p1Basis(grid_->leafView());
MassMatrixType massMatrix;
OperatorAssembler<P1Basis,P1Basis> globalAssembler(p1Basis,p1Basis);
MassAssembler<GridType, typename P1Basis::LocalFiniteElement, typename P1Basis::LocalFiniteElement,BlockType> localMassAssembler;
globalAssembler.assemble(localMassAssembler,massMatrix);
// /////////////////////////////////////////////////////////////
// Lump and invert mass matrix
// /////////////////////////////////////////////////////////////
#ifdef LUMP_MATRIX
BDMatrix<FieldMatrix<double,1,1> > invMassMatrix(baseSet.size());
/* BDMatrix<FieldMatrix<double,1,1> > invMassMatrix(baseSet.size());
invMassMatrix = 0;
for (int i=0; i<baseSet.size(); i++) {
......@@ -60,33 +66,40 @@ public:
invMassMatrix[i][i] += massMatrix[i][j];
}
*/
BDMatrix<BlockType> invMassMatrix;
invMassMatrix = 0;
lumpMatrix(massMatrix,invMassMatrix);
invMassMatrix.invert();
#endif
BlockVector<SymmTensor> unscaledP1Stress(indexSet.size(dim));
BlockVector<SymmetricTensor<dim>> unscaledP1Stress(indexSet.size(dim));
unscaledP1Stress = 0;
LeafElementIterator eIt = grid_->template leafbegin<0>();
LeafElementIterator eEndIt = grid_->template leafend<0>();
typedef Dune::PQkLocalFiniteElementCache<typename GridType::ctype, double, dim, 1> FiniteElementCache;
FiniteElementCache cache;
for (; eIt != eEndIt; ++eIt) {
// Element type
GeometryType gt = eIt->geometry().type();
// First order Lagrange shape functions
const typename LagrangeShapeFunctionSetContainer<double,double,dim>::value_type& baseSet
= LagrangeShapeFunctions<double,double,dim>::general(gt,1);
const typename FiniteElementCache::FiniteElementType& finiteElement = cache.get(gt);
// /////////////////////////////////////////
// Extract the solution on this element
// /////////////////////////////////////////
VectorType localSolution(baseSet.size());
VectorType localSolution(finiteElement.localBasis().size());
for (int i=0; i<baseSet.size(); i++)
localSolution[i] = x[indexSet.template subIndex<dim>(*eIt,i)];
for (int i=0; i<finiteElement.localBasis().size(); i++)
localSolution[i] = x[indexSet.subIndex(*eIt,i,dim)];
/** \todo Determine correct quadrature order */
int quadOrder = (gt.isSimplex()) ? 2 : dim;
......@@ -110,18 +123,18 @@ public:
// /////////////////////////////////////////////
// Compute p0 stress at this quadrature node
// /////////////////////////////////////////////
SymmTensor p0Strain(0);
SymmetricTensor<dim> p0Strain(0.0);
for (int i=0; i<baseSet.size(); i++) {
// evaluate gradients at Gauss points
FieldVector<double,dim> grad, temp;
for (int j=0; j<dim; j++)
temp[j] = baseSet[i].evaluateDerivative(0,j,quadPos);
grad = 0;
jac.umv(temp,grad); // transform gradient to global coordinates
// evaluate gradients at Gauss points
std::vector<FieldMatrix<double,1,dim> > temp;
FieldVector<double,dim> grad;
finiteElement.localBasis().evaluateJacobian(quadPos,temp);
for (int i=0; i<finiteElement.localBasis().size(); i++) {
grad = 0;
jac.umv(temp[i][0],grad); // transform gradient to global coordinates
for (int j=0; j<dim; j++) {
......@@ -130,7 +143,7 @@ public:
Grad[j] = grad;
/* Computes the linear strain tensor from the deformation gradient*/
SymmTensor scaledStrain;
SymmetricTensor<dim> scaledStrain;
computeStrain(Grad,scaledStrain);
scaledStrain *= localSolution[i][j];
......@@ -141,16 +154,19 @@ public:
}
// compute stress
SymmTensor p0Stress = hookeTimesStrain(p0Strain);
SymmetricTensor<dim> p0Stress = hookeTimesStrain(p0Strain);
std::vector<FieldVector<double,1>>values;
finiteElement.localBasis().evaluateFunction(quadPos,values);
// loop over test function number
for (int row=0; row<baseSet.size(); row++) {
for (int row=0; row<finiteElement.localBasis().size(); row++) {
int globalRow = indexSet.template subIndex<dim>(*eIt,row);
int globalRow = indexSet.subIndex(*eIt,row,dim);
for (int rcomp=0; rcomp<SymmTensor::size; rcomp++) {
for (int rcomp=0; rcomp<SymmetricTensor<dim>::size; rcomp++) {
unscaledP1Stress[globalRow][rcomp] += p0Stress[rcomp]*baseSet[row].evaluateFunction(0,quadPos) * factor;
unscaledP1Stress[globalRow][rcomp] += p0Stress[rcomp]*values[row] * factor;
}
......@@ -164,24 +180,28 @@ public:
//
// /////////////////////////////////////////////////////////////
BlockVector<SymmTensor> elementP1Stress(indexSet.size(dim));
BlockVector<SymmetricTensor<dim>> elementP1Stress(indexSet.size(dim));
elementP1Stress = 0;
#ifdef LUMP_MATRIX
elementP1Stress = unscaledP1Stress;
/*elementP1Stress = unscaledP1Stress;
for (int i=0; i<baseSet.size(); i++)
elementP1Stress *= invMassMatrix[i][i];
*/
BlockVector<SymmetricTensor<dim>> tmp(indexSet.size(dim));
tmp = unscaledP1Stress;
invMassMatrix.mv(tmp,elementP1Stress);
#else
// Technicality: turn the matrix into a linear operator
MatrixAdapter<MassMatrixType, BlockVector<SymmTensor>, BlockVector<SymmTensor> > op(*massMatrix);
MatrixAdapter<MassMatrixType, BlockVector<SymmetricTensor<dim>>, BlockVector<SymmetricTensor<dim>> > op(massMatrix);
// A preconditioner
SeqILU0<MassMatrixType, BlockVector<SymmTensor>, BlockVector<SymmTensor> > ilu0(*massMatrix,1.0);
SeqILU0<MassMatrixType, BlockVector<SymmetricTensor<dim>>, BlockVector<SymmetricTensor<dim>> > ilu0(massMatrix,1.0);
// A preconditioned conjugate-gradient solver
int numIt = 100;
CGSolver<BlockVector<SymmTensor> > cg(op,ilu0,1E-16,numIt,0);
CGSolver<BlockVector<SymmetricTensor<dim>> > cg(op,ilu0,1E-16,numIt,0);
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
......@@ -198,18 +218,17 @@ public:
// Element type
GeometryType gt = eIt->geometry().type();
// First order Lagrange shape functions
const typename LagrangeShapeFunctionSetContainer<double,double,dim>::value_type& baseSet
= LagrangeShapeFunctions<double,double,dim>::general(gt,1);
// First order finite element
const typename FiniteElementCache::FiniteElementType& finiteElement = cache.get(gt);
// /////////////////////////////////////////
// Extract the solution on this element
// /////////////////////////////////////////
VectorType localSolution(baseSet.size());
VectorType localSolution(finiteElement.localBasis().size());
for (int i=0; i<baseSet.size(); i++)
localSolution[i] = x[indexSet.template subIndex<dim>(*eIt,i)];
for (int i=0; i<finiteElement.localBasis().size(); i++)
localSolution[i] = x[indexSet.subIndex(*eIt,i,dim)];
/** \todo Determine correct quadrature order */
int quadOrder = (gt.isSimplex()) ? 2 : dim;
......@@ -236,18 +255,21 @@ public:
double detjac = eIt->geometry().integrationElement(quadPos);
// Stresses at this quadrature point
SymmTensor p0Stress(0), p1Stress(0);
SymmetricTensor<dim> p0Stress(0.0), p1Stress(0.0);
// evaluate gradients at Gauss points
std::vector<FieldMatrix<double,1,dim>>temp;
std::vector<FieldVector<double,1>>values;
FieldVector<double,dim> grad;
finiteElement.localBasis().evaluateJacobian(quadPos,temp);
finiteElement.localBasis().evaluateFunction(quadPos,values);
// loop over test function number
for (int i=0; i<baseSet.size(); i++) {
for (int i=0; i<finiteElement.localBasis().size(); i++) {
// evaluate gradient of the current shape function at the Gauss point
FieldVector<double,dim> grad, temp;
for (int l=0; l<dim; l++)
temp[l] = baseSet[i].evaluateDerivative(0,l,quadPos);
grad = 0;
jac.umv(temp,grad); // transform gradient to global coordinates
jac.umv(temp[i][0],grad); // transform gradient to global coordinates
for (int j=0; j<dim; j++) {
......@@ -256,7 +278,7 @@ public:
Grad[j] = grad;
/* Computes the linear strain tensor from the deformation gradient*/
SymmTensor shapeFunctionStrain;
SymmetricTensor<dim> shapeFunctionStrain;
computeStrain(Grad,shapeFunctionStrain);
shapeFunctionStrain *= localSolution[i][j];
......@@ -267,8 +289,8 @@ public:
}
SymmTensor scaledStress = elementP1Stress[i];
scaledStress *= baseSet[i].evaluateFunction(0,quadPos);
SymmetricTensor<dim> scaledStress = elementP1Stress[i];
scaledStress *= values[i];
p1Stress += scaledStress;
}
......@@ -285,94 +307,35 @@ public:
}
private:
// computes the linear strain from the deformation gradient
void computeStrain(const Dune::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(Dune::NotImplemented, "No elasticity assembler for " << dim << "-dimensional problems");
}
SymmTensor hookeTimesStrain(const SymmTensor& strain) const {
if (dim==3) {
// compute Hooke Tensor
Dune::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
Dune::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(Dune::NotImplemented, "No elasticity assembler for " << dim << "-dimensional problems");
}
void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad, SymmetricTensor<dim>& strain) const
{
for (int i=0; i<dim ; ++i)
{
strain(i,i) = grad[i][i];
for (int j=i+1; j<dim; ++j)
strain(i,j) = 0.5*(grad[i][j] + grad[j][i]);
}
}
SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
SymmetricTensor<dim> stress = strain;
stress *= E_/(1+nu_);
// Compute the trace of the strain
double trace = 0;
for (int i=0; i<dim; i++)
trace += strain[i];
double f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * trace;
for (int i=0; i<dim; i++)
stress[i] += f;
return stress;
}
// /////////////////////////////////
// Data members
......
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