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

change the structure of the localogdenassembler and ogdenassembler:

localogdenassembler is now derived from localassembler.hh and
ogdenassembler is now derived from operatorassembler.hh                                    

[[Imported from SVN: r10472]]
parent 5df2ae7a
......@@ -5,49 +5,51 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
#include <dune/ag-common/functionspacebases/p1nodalbasis.hh>
#include <dune/grid/common/quadraturerules.hh>
template <class GridView, class RT>
void OgdenMaterialLocalStiffness<GridView, RT>::
assemble(const Entity&entity,
const Dune::BlockVector<Dune::FieldVector<double, dim> >& localSolution,
int k)
template <class GridType, class TrialLocalFE, class AnsatzLocalFE>
void OgdenMaterialLocalStiffness<GridType,TrialLocalFE,AnsatzLocalFE>::
assemble(const ElementPointer& 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
{
Dune::GeometryType gt = entity.type();
typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
int rows = tFE.localBasis().size();
int cols = aFE.localBasis().size();
const typename Dune::LagrangeShapeFunctionSetContainer<DT,RT,dim>::value_type& sfs
= Dune::LagrangeShapeFunctions<DT,RT,dim>::general(gt,k);
localMatrix.setSize(rows,cols);
localMatrix = 0.0;
localRhs.resize(rows);
localRhs = 0;
/* ndof is the number of vectors of the element */
int ndof = sfs.size();
int ndof = tFE.localBasis().size();
int ncomp = dim;
this->setcurrentsize(sfs.size());
// Initialize element matrix and local rhs
for (int i=0; i<sfs.size(); i++) {
this->b[i] = 0;
this->bctype[i][0] = Dune::BoundaryConditions::neumann;
for (int j=0; j<sfs.size(); j++)
this->A[i][j] = 0;
}
/** \todo Is this correct/sufficient? */
int polOrd = 3;
// Get quadrature rule
const Dune::QuadratureRule<double, dim>& quad = Dune::QuadratureRules<double, dim>::rule(gt, polOrd);
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<DT,dim> quadPos = quad[ip].position();
const Dune::FieldVector<double,dim> quadPos = quad[ip].position();
// calc Jacobian inverse before integration element is evaluated
const Dune::FieldMatrix<DT,dim,dim>& inv = entity.geometry().jacobianInverseTransposed(quadPos);
const DT integrationElement = entity.geometry().integrationElement(quadPos);
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;
......@@ -55,20 +57,18 @@ assemble(const Entity&entity,
/**********************************************/
/* compute gradients of the shape functions */
/**********************************************/
Dune::FieldVector<DT, dim> shape_grads[sfs.size()];
for (int dof=0; dof<ndof; dof++) {
//baseSet.jacobian(dof, quadPos, shape_grads[dof]);
for (int i=0; i<dim; i++)
shape_grads[dof][i] = sfs[dof].evaluateDerivative(0, i, quadPos);
// get gradients of shape functions
tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
// compute gradients of base functions
for (int i=0; i<gradients.size(); ++i) {
// multiply with jacobian inverse
Dune::FieldVector<double,dim> tmp(0);
inv.umv(shape_grads[dof], tmp);
shape_grads[dof] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
}
// transform gradients
gradients[i] = 0.0;
inv.umv(referenceGradients[i][0], gradients[i]);
}
/****************************************************/
/* The deformation gradients of the shape functions */
......@@ -79,7 +79,7 @@ assemble(const Entity&entity,
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) ? shape_grads[dof][j] : 0;
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);
......@@ -98,7 +98,7 @@ assemble(const Entity&entity,
parDer[i][j] = 0;
for (int k=0; k<ndof; k++)
parDer[i][j] += localSolution[k][i]*shape_grads[k][j];
parDer[i][j] += localSolution[k][i]*gradients[k][j];
}
......@@ -111,14 +111,16 @@ assemble(const Entity&entity,
/********************************************************/
Dune::FieldMatrix<double,dim,dim> fu;
Dune::OgdenAssembler<typename GridView::Grid>::Ogden_du_tmpl(parDer, d_, lambda_, mu_, 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<typename GridView::Grid>::Ogden_dudu_tmpl(parDer, d_, lambda_, mu_, fuu);
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
......@@ -140,7 +142,7 @@ assemble(const Entity&entity,
// Second inner tensor product.
for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++)
this->A[i/dim][j/dim][i%dim][j%dim] += nablaPhi[i][k][l] * tmp[k][l] * weight;
localMatrix[i/dim][j/dim][i%dim][j%dim] += nablaPhi[i][k][l] * tmp[k][l] * weight;
}
......@@ -148,7 +150,7 @@ assemble(const Entity&entity,
/** \todo The volume forces part is missing! */
for (int k=0; k<dim; k++)
for (int l=0; l<dim; l++)
this->b[i/dim][i%dim] -= nablaPhi[i][k][l] * fu[k][l] * weight;
localRhs[i/dim][i%dim] -= nablaPhi[i][k][l] * fu[k][l] * weight;
}
......@@ -183,129 +185,93 @@ assemble(const Entity&entity,
}
template <class GridType>
void Dune::OgdenAssembler<GridType>::
getNeighborsPerVertex(MatrixIndexSet& nb) const
{
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
int n = grid_->size(grid_->maxLevel(), dim);
nb.resize(n, n);
ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
for (; it!=endit; ++it) {
for (int i=0; i<it->template count<dim>(); i++) {
for (int j=0; j<it->template count<dim>(); j++) {
int iIdx = indexSet.template subIndex<dim>(*it,i);
int jIdx = indexSet.template subIndex<dim>(*it,j);
nb.add(iIdx, jIdx);
}
}
}
}
template <class GridType>
void Dune::OgdenAssembler<GridType>::
assembleMatrix(const BlockVector<FieldVector<double, dim> >& sol,
BlockVector<FieldVector<double, dim> >& rhs)
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)
{
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
//use method from the base class to compute indices of the matrix
int rows = this->tBasis_.size();
int cols = this->aBasis_.size();
MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
Dune::MatrixIndexSet indices(rows, cols);
/** \todo Reuse matrix, the occupation pattern won't change! */
if (matrix_)
delete(matrix_);
matrix_ = new BCRSMatrix<MatrixBlock>;
// Create local assembler
OgdenMaterialLocalStiffness<typename GridView::Grid,typename TrialBasis::LocalFiniteElement,typename AnsatzBasis::LocalFiniteElement> localAssembler(lambda_, mu_, d_);
addIndices(localAssembler, indices, false);
neighborsPerVertex.exportIdx(*matrix_);
*matrix_ = 0;
indices.exportIdx(globalMatrix);
globalMatrix=0.0;
// Create local assembler
OgdenMaterialLocalStiffness<typename GridType::LeafGridView,double> localStiffness(lambda_, mu_, d_);
typedef typename OgdenMaterialLocalStiffness<typename GridView::Grid,typename TrialBasis::LocalFiniteElement,typename AnsatzBasis::LocalFiniteElement>::LocalMatrix LocalMatrix;
ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
ElementIterator it = this->tBasis_.getGridView().template begin<0>();
ElementIterator endIt = this->tBasis_.getGridView().template end<0>();
for( ; it != endit; ++it ) {
//const BaseFunctionSetType& baseSet = functionSpace_.getBaseFunctionSet( *it );
const LagrangeShapeFunctionSet<double, double, dim> & baseSet
= Dune::LagrangeShapeFunctions<double, double, dim>::general(it->type(), elementOrder);
const int numOfBaseFct = baseSet.size();
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(numOfBaseFct);
BlockVector<FieldVector<double, dim> > localSolution(aFE.localBasis().size()), localRhs(tFE.localBasis().size());
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.template subIndex<dim>(*it,i)];
for (int i=0; i<aFE.localBasis().size(); i++)
localSolution[i] = sol[this->aBasis_.index(*it,i)];
// setup matrix
localStiffness.template assemble(*it, localSolution);
// 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<numOfBaseFct; i++) {
int row = indexSet.template subIndex<dim>(*it,i);
for (int j=0; j<numOfBaseFct; j++ ) {
int col = indexSet.template subIndex<dim>(*it,j);
(*matrix_)[row][col] += localStiffness.mat(i,j);
}
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<numOfBaseFct; i++)
rhs[indexSet.template subIndex<dim>(*it,i)] += localStiffness.rhs(i);
for (int i=0; i<tFE.localBasis().size(); i++)
rhs[this->tBasis_.index(*it, i)] += localRhs[i];
}
}
template <class GridType>
double Dune::OgdenAssembler<GridType>::
template <class TrialBasis,class AnsatzBasis>
double Dune::OgdenAssembler<TrialBasis,AnsatzBasis>::
computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
{
const int maxlevel = grid_->maxLevel();
double energy = 0;
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(maxlevel);
if (sol.size()!=grid_->size(maxlevel, dim))
if (sol.size()!=this->aBasis_.size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
ElementIterator it = grid_->template lbegin<0>(maxlevel);
ElementIterator endIt = grid_->template lend<0>(maxlevel);
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 LagrangeShapeFunctionSet<double, double, dim> & baseSet = Dune::LagrangeShapeFunctions<double, double, dim>::general(it->geometry().type(), elementOrder);
const int numOfBaseFct = baseSet.size();
const typename AnsatzBasis::LocalFiniteElement& localFE = this->aBasis_.getLocalFiniteElement(*it);
FieldVector<double, dim> localSolution[localFE.localBasis().size()];
FieldVector<double, dim> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.template subIndex<dim>(*it,i)];
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);
......@@ -323,18 +289,17 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
// ///////////////////////////////////////
// Compute deformation gradient
// ///////////////////////////////////////
std::vector<FieldVector<double, dim> > shape_grads(numOfBaseFct);
std::vector<FieldMatrix<double,1, dim> > shape_grads;
localFE.localBasis().evaluateJacobian(quadPos,shape_grads);
for (int dof=0; dof<numOfBaseFct; dof++) {
for (int dof=0; dof<localFE.localBasis().size(); dof++) {
for (int i=0; i<dim; i++)
shape_grads[dof][i] = baseSet[dof].evaluateDerivative(0, i, quadPos);
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
// multiply with jacobian inverse
FieldVector<double, dim> tmp(0);
inv.umv(shape_grads[dof], tmp);
shape_grads[dof] = tmp;
inv.umv(shape_grads[dof][0], tmp);
shape_grads[dof][0] = tmp;
//std::cout << "Gradient " << dof << ": " << shape_grads[dof] << std::endl;
}
......@@ -345,8 +310,8 @@ computeEnergy(const BlockVector<FieldVector<double, dim> >& sol) const
defGrad[i][j] = 0;
for (int k=0; k<numOfBaseFct; k++)
defGrad[i][j] += localSolution[k][i]*shape_grads[k][j];
for (int k=0; k<localFE.localBasis().size(); k++)
defGrad[i][j] += localSolution[k][i]*shape_grads[k][0][j];
}
......
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