Skip to content
Snippets Groups Projects
Commit d754ae6d authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Implement the energy that Patrizio and Jendrik call W_magic

parent 846823ae
No related branches found
No related tags found
No related merge requests found
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/parametertree.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune::Elasticity {
template<class GridView, class LocalFiniteElement, class field_type=double>
class NeffMagicEnergy
: public Elasticity::LocalEnergy<GridView,
LocalFiniteElement,
std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
{
// grid types
typedef typename GridView::Grid::ctype DT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
enum {gridDim=GridView::dimension};
enum {dim=GridView::dimension};
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
NeffMagicEnergy(const Dune::ParameterTree& parameters)
{}
/** \brief Assemble the energy for a single element */
field_type energy (const Entity& e,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
NeffMagicEnergy<GridView, LocalFiniteElement, field_type>::
energy(const Entity& element,
const LocalFiniteElement& localFiniteElement,
const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
{
assert(element.type() == localFiniteElement.type());
field_type energy = 0;
// store gradients of shape functions and base functions
std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
for (const auto& qp : quad)
{
const DT integrationElement = element.geometry().integrationElement(qp.position());
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
// get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<gridDim; j++)
derivative[j].axpy(localConfiguration[i][j], gradients[i]);
auto det = derivative.determinant();
auto normSquared = derivative.frobenius_norm2();
auto K = normSquared / (2*det);
using std::acosh;
auto energyDensity = K + std::sqrt(K*K-1) - acosh(K) + log(det);
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
......@@ -39,6 +39,7 @@
#include <dune/elasticity/materials/expacoshenergy.hh>
#include <dune/elasticity/materials/nefflogenergy.hh>
#include <dune/elasticity/materials/nefflogenergy2.hh>
#include <dune/elasticity/materials/neffmagicenergy.hh>
#include <dune/elasticity/materials/neffreducedenergy.hh>
// grid dimension
......@@ -271,6 +272,11 @@ int main (int argc, char *argv[]) try
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters);
if (parameterSet.get<std::string>("energy") == "magic")
elasticEnergy = std::make_shared<Elasticity::NeffMagicEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters);
if(!elasticEnergy)
DUNE_THROW(Exception, "Error: Selected energy not available!");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment