Select Git revision
ogdenassembler.cc
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
ogdenassembler.cc 13.72 KiB
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/geometry/quadraturerules.hh>
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
void OgdenMaterialLocalStiffness<GridType,TrialLocalFE,AnsatzLocalFE>::
assemble(const Element& element,
const Dune::BlockVector<Dune::FieldVector<double,dim> >& localSolution,
LocalMatrix& localMatrix,
Dune::BlockVector<Dune::FieldVector<double,dim> >& localRhs,
const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
{
typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
int rows = tFE.localBasis().size();
int cols = aFE.localBasis().size();
localMatrix.setSize(rows,cols);
localMatrix = 0.0;
localRhs.resize(rows);
localRhs = 0;
/* ndof is the number of vectors of the element */
int ndof = tFE.localBasis().size();
int ncomp = dim;
/** \todo Is this correct/sufficient? */
int polOrd = 3;
// Get quadrature rule
const Dune::QuadratureRule<double, dim>& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), polOrd);
// store gradients of shape functions and base functions
std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
std::vector<Dune::FieldVector<double,dim> > gradients(tFE.localBasis().size());
/* Loop over all integration points */
for (int ip=0; ip<quad.size(); ip++) {
// Local position of the quadrature point
const Dune::FieldVector<double,dim> quadPos = quad[ip].position();
// calc Jacobian inverse before integration element is evaluated
const Dune::FieldMatrix<double,dim,dim>& inv = element.geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = element.geometry().integrationElement(quadPos);
/* Compute the weight of the current integration point */
double weight = quad[ip].weight() * integrationElement;
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
// get gradients of shape functions
tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradients of base functions
for (int i=0; i<gradients.size(); ++i) {
// transform gradients
gradients[i] = 0.0;
inv.umv(referenceGradients[i][0], gradients[i]);
}
/****************************************************/
/* The deformation gradients of the shape functions */
/****************************************************/
Dune::FieldMatrix<double,dim,dim> nablaPhi[100*dim];
for (int dof=0; dof<ndof; dof++)
for (int comp=0; comp<ncomp; comp++)
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
nablaPhi[dof*ncomp+comp][i][j] = (i==comp) ? gradients[dof][j] : 0;
// for (int i=0; i<ndof*ncomp; i++) {
// printf("Deformation gradient %d\n", i);
// std::cout << nablaPhi[i] << std::endl;
// }
/****************************************************/
// the partial derivatives of u
// in formulas: parDer(i,j) = $\partial u_i / \partial x_j$
/****************************************************/
Dune::FieldMatrix<double, dim, dim> parDer;
for (int i=0; i<ncomp; i++) {
for (int j=0; j<dim; j++) {
parDer[i][j] = 0;
for (int k=0; k<ndof; k++)
parDer[i][j] += localSolution[k][i]*gradients[k][j];
}
}
//std::cout << "parDer:\n" << parDer << std::endl;
/********************************************************/
/* Compute first derivative of Ogden Energy functional */
/********************************************************/
Dune::FieldMatrix<double,dim,dim> fu;
//it doesn't matter which basis we choose here since the function Ogden_du_tmpl doesn't depend on it
Dune::OgdenAssembler<P1NodalBasis<typename GridType::LeafGridView,double>,P1NodalBasis<typename GridType::LeafGridView,double> >::Ogden_du_tmpl(parDer, d_, lambda_, mu_, fu);
/********************************************************/
/* Compute second derivative of Ogden Energy functional */
/********************************************************/
double fuu[dim][dim][dim][dim];
Dune::OgdenAssembler<P1NodalBasis<typename GridType::LeafGridView,double>,P1NodalBasis<typename GridType::LeafGridView,double> >::Ogden_dudu_tmpl(parDer, d_, lambda_, mu_, fuu);
// loop over all entries of the element stiffness matrix
// and add the contribution of the current quadrature point
for (int i=0; i<ndof*ncomp; i++) {
for (int j=0; j<ndof*ncomp; j++) {
// Compute \nabla \phi_i * fuu * \nabla \phi_j
// First tensor product:
Dune::FieldMatrix<double,dim,dim> tmp(0);
for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++)
for (int m=0; m<dim; m++)
for (int n=0; n<dim; n++)
tmp[k][l] += fuu[k][l][m][n]*nablaPhi[j][m][n];
// Second inner tensor product.
for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++)
localMatrix[i/dim][j/dim][i%dim][j%dim] += nablaPhi[i][k][l] * tmp[k][l] * weight;
}
/* add the nonlinear part to the right hand side */
/** \todo The volume forces part is missing! */
for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++)
localRhs[i/dim][i%dim] -= nablaPhi[i][k][l] * fu[k][l] * weight;
}
}
#if 0
static int eleme = 0;
printf("********* Element %d **********\n", eleme++);
for (int row=0; row<matSize; row++) {
for (int rcomp=0; rcomp<dim; rcomp++) {
for (int col=0; col<matSize; col++) {
for (int ccomp=0; ccomp<dim; ccomp++)
std::cout << mat[row][col][rcomp][ccomp] << " ";
std::cout << " ";
}
std::cout << std::endl;
}
std::cout << std::endl;
}
exit(0);
#endif
}
template <class TrialBasis,class AnsatzBasis>
void Dune::OgdenAssembler<TrialBasis,AnsatzBasis>::
assembleProblem(BCRSMatrix<FieldMatrix<double,dim,dim> >& globalMatrix,
const BlockVector<FieldVector<double, dim> >& sol,
BlockVector<FieldVector<double, dim> >& rhs)
{
//use method from the base class to compute indices of the matrix
int rows = this->tBasis_.size();
int cols = this->aBasis_.size();
Dune::MatrixIndexSet indices(rows, cols);
// Create local assembler
OgdenMaterialLocalStiffness<typename GridView::Grid,typename TrialBasis::LocalFiniteElement,typename AnsatzBasis::LocalFiniteElement> localAssembler(lambda_, mu_, d_);
this->addIndices(localAssembler, indices, false);
indices.exportIdx(globalMatrix);
globalMatrix=0.0;
typedef typename OgdenMaterialLocalStiffness<typename GridView::Grid,typename TrialBasis::LocalFiniteElement,typename AnsatzBasis::LocalFiniteElement>::LocalMatrix LocalMatrix;
ElementIterator it = this->tBasis_.getGridView().template begin<0>();
ElementIterator endIt = this->tBasis_.getGridView().template end<0>();
for( ; it != endIt; ++it ) {
// get shape functions
const typename TrialBasis::LocalFiniteElement& tFE = this->tBasis_.getLocalFiniteElement(*it);
const typename AnsatzBasis::LocalFiniteElement& aFE = this->aBasis_.getLocalFiniteElement(*it);
LocalMatrix localA(tFE.localBasis().size(), aFE.localBasis().size());
// Extract local solution
BlockVector<FieldVector<double, dim> > localSolution(aFE.localBasis().size()), localRhs(tFE.localBasis().size());
for (int i=0; i<aFE.localBasis().size(); i++)
localSolution[i] = sol[this->aBasis_.index(*it,i)];
// assemble local matrix and rhs
localAssembler.assemble(*it, localSolution, localA, localRhs, tFE, aFE);
// Add element matrix to global stiffness matrix
for(int i=0; i<tFE.localBasis().size(); i++) {
int rowIndex = this->tBasis_.index(*it, i);
for (int j=0; j<aFE.localBasis().size(); j++ ) {
int colIndex = this->aBasis_.index(*it, j);
globalMatrix[rowIndex][colIndex] += localA[i][j];
}
}
// Add elements of local rhs to global rhs
for (int i=0; i<tFE.localBasis().size(); i++)
rhs[this->tBasis_.index(*it, i)] += localRhs[i];
}
}
template <class TrialBasis,class AnsatzBasis>
double Dune::OgdenAssembler<TrialBasis,AnsatzBasis>::
computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
{
double energy = 0;
if (sol.size()!=this->aBasis_.size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
ElementIterator it = this->tBasis_.getGridView().template begin<0>();
ElementIterator endIt = this->tBasis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
// Extract local solution on this element
const typename AnsatzBasis::LocalFiniteElement& localFE = this->aBasis_.getLocalFiniteElement(*it);
FieldVector<double, dim> localSolution[localFE.localBasis().size()];
for (int i=0; i<localFE.localBasis().size(); i++)
localSolution[i] = sol[this->aBasis_.index(*it, i)];
// Get quadrature rule
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(it->type(), polOrd);
for (int pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const FieldVector<double,dim>& quadPos = quad[pt].position();
const MatrixBlock& inv = it->geometry().jacobianInverseTransposed(quadPos);
const double integrationElement = it->geometry().integrationElement(quadPos);
double weight = quad[pt].weight() * integrationElement;
// ///////////////////////////////////////
// Compute deformation gradient
// ///////////////////////////////////////
std::vector<FieldMatrix<double,1, dim> > shape_grads;
localFE.localBasis().evaluateJacobian(quadPos,shape_grads);
for (int dof=0; dof<localFE.localBasis().size(); dof++) {
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
// multiply with jacobian inverse
FieldVector<double, dim> tmp(0);
inv.umv(shape_grads[dof][0], tmp);
shape_grads[dof][0] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
}
FieldMatrix<double, dim, dim> defGrad;
for (int i=0; i<dim; i++) {
for (int j=0; j<dim; j++) {
defGrad[i][j] = 0;
for (int k=0; k<localFE.localBasis().size(); k++)
defGrad[i][j] += localSolution[k][i]*shape_grads[k][0][j];
}
}
// ///////////////////////////////////////
// Compute nonlinear strain tensor
// ///////////////////////////////////////
FieldMatrix<double,dim,dim> E;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
E[i][j] = E_val(defGrad, i, j);
// ///////////////////////////////////////
// Sum up result
// ///////////////////////////////////////
double trE = 0;
for (int i=0; i<dim; i++)
trE += E[i][i];
FieldMatrix<double,dim,dim> EE = E;
EE.leftmultiply(E);
double trEE = 0;
for (int i=0; i<dim; i++)
trEE += EE[i][i];
double a = -d_ * Gamma_x(1);
double b = (lambda_ - d_ * (Gamma_x(1)+Gamma_xx(1)))/2;
double c = mu_ + d_ * Gamma_x(1);
// printf("trE: %g, trEE: %g\n", trE, trEE);
// printf("weight: %g\n", weight);
energy += weight * (a * trE + b*trE*trE + c * trEE + d_ * Gamma(det_val(defGrad)));
#if 0
if (det_val(defGrad) <= 0) {
printf("bad element: %d\n", it->index());
printf("vertices: %d %d %d %d\n", it->template subIndex<3>(0),
it->template subIndex<3>(1),
it->template subIndex<3>(2),
it->template subIndex<3>(3));
}
#endif
// printf("det %g\n", det_val(defGrad));
// printf("energy %g\n", energy);
}
}
return energy;
}