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

Move all quasiconvexity test densities into a separate file

parent dfa47870
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 440 additions and 1160 deletions
dune/elasticity/materials/CMakeLists.txt
+ 1
− 0
View file @ aaab57cf
...@@ -15,6 +15,7 @@ install(FILES ...@@ -15,6 +15,7 @@ install(FILES
neohookeenergy.hh neohookeenergy.hh
neumannenergy.hh neumannenergy.hh
ogdenmaterial.hh ogdenmaterial.hh
quasiconvexitytestdensities.hh
stvenantkirchhoffdensity.hh stvenantkirchhoffdensity.hh
stvenantkirchhoffenergy.hh stvenantkirchhoffenergy.hh
sumenergy.hh sumenergy.hh
......
dune/elasticity/materials/coshenergy.hh deleted 100644 → 0
+ 0
− 105
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class CosHEnergy final
: 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
*/
CosHEnergy(const Dune::ParameterTree& parameters)
{
eps_ = parameters.get<double>("eps");
K0_ = parameters.get<double>("K0");
}
/** \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;
double eps_;
double K0_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
CosHEnergy<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::cosh;
using std::pow;
auto h2 = cosh(pow( (K-K0_)/eps_, 4)) - 1;
auto energyDensity = h2 + 0;//100*f2;
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_COSHENERGY_HH
dune/elasticity/materials/dacorognaenergy.hh deleted 100644 → 0
+ 0
− 97
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_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 DacorognaEnergy final
: 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
*/
DacorognaEnergy(const Dune::ParameterTree& parameters)
{
gamma_ = parameters.template get<double>("gamma");
}
/** \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;
double gamma_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
DacorognaEnergy<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 energyDensity = normSquared * ( normSquared - gamma_ * det);
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_DACOROGNAENERGY_HH
dune/elasticity/materials/expacoshenergy.hh deleted 100644 → 0
+ 0
− 106
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class ExpACosHEnergy final
: 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
*/
ExpACosHEnergy(const Dune::ParameterTree& parameters)
{
c1_ = parameters.get<double>("c1");
c2_ = parameters.get<double>("c2");
}
/** \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;
double c1_;
double c2_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
ExpACosHEnergy<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 h = exp(c1_ * acosh(K) * acosh(K));
auto f = - c2_ * log(det);
auto energyDensity = h + f;
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH
dune/elasticity/materials/fungenergy.hh deleted 100644 → 0
+ 0
− 121
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_FUNGENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_FUNGENERGY_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 FungEnergy final
: 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
*/
FungEnergy(const Dune::ParameterTree& parameters)
{
c_ = parameters.template get<double>("c");
A_[0] = parameters.template get<double>("A1");
A_[1] = parameters.template get<double>("A2");
A_[2] = parameters.template get<double>("A3");
A_[3] = parameters.template get<double>("A4");
A_[4] = parameters.template get<double>("A5");
A_[5] = parameters.template get<double>("A6");
}
/** \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;
double c_;
std::array<double, 6> A_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
FungEnergy<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]);
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]);
FieldMatrix<field_type,gridDim,gridDim> greenLagrangeStrain(0);
for (int i=0; i<gridDim; i++)
for (int j=0; j<gridDim; j++)
{
greenLagrangeStrain[i][j] = -(i==j);
for (int k=0; k<gridDim; k++)
greenLagrangeStrain[i][j] += derivative[k][i] * derivative[k][j];
}
//auto psi = A1 * E11*E11 + A2 * E22 * E22 + 2*A3*E11*E22 + A4 * E12 * E12 + 2*A5 * E12 * E11 + 2*A6 * E12 * E22;
auto psi = A_[0] * greenLagrangeStrain[0][0] * greenLagrangeStrain[0][0]
+ A_[1] * greenLagrangeStrain[1][1] * greenLagrangeStrain[1][1]
+ 2*A_[2]*greenLagrangeStrain[0][0] * greenLagrangeStrain[1][1]
+ A_[3] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][1]
+ 2*A_[4] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][0]
+ 2*A_[5] * greenLagrangeStrain[0][1] * greenLagrangeStrain[1][1];
using std::exp;
auto energyDensity = 0.5 * c_ * (exp(psi) - 1);
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_FUNGENERGY_HH
dune/elasticity/materials/neffenergy.hh deleted 100644 → 0
+ 0
− 96
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class NeffEnergy final
: 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
*/
NeffEnergy(const Dune::ParameterTree& parameters)
{
scaling_ = parameters.get<double>("scaling");
}
/** \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;
private:
double scaling_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
NeffEnergy<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();
energy += qp.weight() * integrationElement *
((scaling_/(4*std::sqrt(3))) * (normSquared*normSquared/(det*det) -4) + std::log(det) - det);
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFENERGY_HH
dune/elasticity/materials/nefflogenergy.hh deleted 100644 → 0
+ 0
− 94
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class NeffLogEnergy final
: 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
*/
NeffLogEnergy(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
NeffLogEnergy<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 energyDensity = std::log( normSquared / (2*det) + std::sqrt(normSquared*normSquared/(4*det*det) -1) );
energyDensity = energyDensity * energyDensity;
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY_HH
dune/elasticity/materials/nefflogenergy2.hh deleted 100644 → 0
+ 0
− 104
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class NeffLogEnergy2 final
: 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
*/
NeffLogEnergy2(const Dune::ParameterTree& parameters,
double w0=1, double w1=1)
{
c0_ = parameters.template get<double>("c0");
weights_ = {w0,w1};
}
/** \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;
double c0_;
// This is needed to compute the two parts of the energy separately
std::array<double,2> weights_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
NeffLogEnergy2<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 = 0.5 * normSquared / det;
auto energyDensity = weights_[0] * c0_ * K * std::log(K) + weights_[1] * std::log(det);
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFLOGENERGY2_HH
dune/elasticity/materials/neffreducedenergy.hh deleted 100644 → 0
+ 0
− 93
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/elasticity/assemblers/localfestiffness.hh>
namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double>
class NeffReducedEnergy final
: 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
*/
NeffReducedEnergy(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
NeffReducedEnergy<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 energyDensity = normSquared / det;
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFREDUCEDENERGY_HH
dune/elasticity/materials/neffw2energy.hh deleted 100644 → 0
+ 0
− 102
View file @ dfa47870
#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
#define DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_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 NeffW2Energy final
: 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
*/
NeffW2Energy(const Dune::ParameterTree& parameters)
{
c_ = 1.0;
if (parameters.hasKey("c"))
c_ = parameters.template get<double>("c");
}
/** \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;
double c_;
};
template <class GridView, class LocalFiniteElement, class field_type>
field_type
NeffW2Energy<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 + acosh(K) - c_*log(det);
energy += qp.weight() * integrationElement * energyDensity;
}
return energy;
}
} // namespace Dune::Elasticity
#endif //#ifndef DUNE_ELASTICITY_MATERIALS_NEFFW2ENERGY_HH
dune/elasticity/materials/quasiconvexitytestdensities.hh 0 → 100644
+ 402
− 0
View file @ aaab57cf
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ELASTICITY_MATERIAL_QUASICONVEXITYTESTDENSITIES_HH
#define DUNE_ELASTICITY_MATERIAL_QUASICONVEXITYTESTDENSITIES_HH
#include <dune/elasticity/materials/localdensity.hh>
namespace Dune:: Elasticity
{
template<int dim, class field_type, class ctype=double>
struct BurkholderDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
BurkholderDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
p_ = parameters.template get<double>("p");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto lambdaMaxSquared = 0.5 * (normSquared + std::sqrt(normSquared*normSquared - 4*det*det));
auto lambdaMax = std::sqrt(lambdaMaxSquared);
return -0.5 * p_ * det * std::pow(lambdaMax,p_-2) - 0.5*(2-p_) * std::pow(lambdaMax,p_);
}
FieldMatrix<ctype,dim,dim> F0_;
double p_;
};
template<int dim, class field_type, class ctype=double>
struct CosHDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
CosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
eps_ = parameters.get<double>("eps");
K0_ = parameters.get<double>("K0");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::cosh;
using std::pow;
auto h2 = cosh(pow( (K-K0_)/eps_, 4)) - 1;
return h2 + 0;//100*f2;
}
FieldMatrix<ctype,dim,dim> F0_;
double eps_;
double K0_;
};
template<int dim, class field_type, class ctype=double>
struct DacorognaDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
DacorognaDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
gamma_ = parameters.template get<double>("gamma");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
return normSquared * ( normSquared - gamma_ * det);
}
FieldMatrix<ctype,dim,dim> F0_;
double gamma_;
};
template<int dim, class field_type, class ctype=double>
struct ExpACosHDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
ExpACosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c1_ = parameters.get<double>("c1");
c2_ = parameters.get<double>("c2");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::acosh;
return exp(c1_ * acosh(K) * acosh(K)) - c2_ * log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c1_;
double c2_;
};
template<int dim, class field_type, class ctype=double>
struct FungDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
FungDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = parameters.template get<double>("c");
A_[0] = parameters.template get<double>("A1");
A_[1] = parameters.template get<double>("A2");
A_[2] = parameters.template get<double>("A3");
A_[3] = parameters.template get<double>("A4");
A_[4] = parameters.template get<double>("A5");
A_[5] = parameters.template get<double>("A6");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
FieldMatrix<field_type,dim,dim> greenLagrangeStrain(0);
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
{
greenLagrangeStrain[i][j] = -(i==j);
for (int k=0; k<dim; k++)
greenLagrangeStrain[i][j] += F[k][i] * F[k][j];
}
//auto psi = A1 * E11*E11 + A2 * E22 * E22 + 2*A3*E11*E22 + A4 * E12 * E12 + 2*A5 * E12 * E11 + 2*A6 * E12 * E22;
auto psi = A_[0] * greenLagrangeStrain[0][0] * greenLagrangeStrain[0][0]
+ A_[1] * greenLagrangeStrain[1][1] * greenLagrangeStrain[1][1]
+ 2*A_[2]*greenLagrangeStrain[0][0] * greenLagrangeStrain[1][1]
+ A_[3] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][1]
+ 2*A_[4] * greenLagrangeStrain[0][1] * greenLagrangeStrain[0][0]
+ 2*A_[5] * greenLagrangeStrain[0][1] * greenLagrangeStrain[1][1];
using std::exp;
return 0.5 * c_ * (exp(psi) - 1);
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
std::array<double, 6> A_;
};
template<int dim, class field_type, class ctype=double>
struct KPowerDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
p_ = parameters.template get<double>("p");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto detF = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*detF);
return std::pow(K, p_);
}
FieldMatrix<ctype,dim,dim> F0_;
double p_;
};
template<int dim, class field_type, class ctype=double>
struct MagicDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = 1.0;
if (parameters.hasKey("c"))
c_ = parameters.template get<double>("c");
}
/** \brief Implement space-dependent coefficients */
double c(const Dune::FieldVector<ctype,dim>& x) const
{
return c_;
#if 0 // Three circles
bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
#endif
// return (x.two_norm() < 0.9) ? 1 : c_;
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::acosh;
using std::sqrt;
return K + sqrt(K*K-1) - acosh(K) + c(x)*log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
template<int dim, class field_type, class ctype=double>
struct NeffDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
NeffDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
scaling_ = parameters.get<double>("scaling");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
return (scaling_/(4*std::sqrt(3))) * (normSquared*normSquared/(det*det) -4) + std::log(det) - det;
}
private:
FieldMatrix<ctype,dim,dim> F0_;
double scaling_;
};
template<int dim, class field_type, class ctype=double>
struct NeffLogDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
NeffLogDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto energyDensity = std::log( normSquared / (2*det) + std::sqrt(normSquared*normSquared/(4*det*det) -1) );
return energyDensity * energyDensity;
}
FieldMatrix<ctype,dim,dim> F0_;
};
template<int dim, class field_type, class ctype=double>
struct NeffLog2Density final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
NeffLog2Density(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters,
double w0=1, double w1=1)
: F0_(F0)
{
c0_ = parameters.template get<double>("c0");
weights_ = {w0,w1};
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = 0.5 * normSquared / det;
return weights_[0] * c0_ * K * std::log(K) + weights_[1] * std::log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c0_;
// This is needed to compute the two parts of the energy separately
std::array<double,2> weights_;
};
template<int dim, class field_type, class ctype=double>
struct NeffReducedDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
NeffReducedDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
return normSquared / det;
}
FieldMatrix<ctype,dim,dim> F0_;
};
template<int dim, class field_type, class ctype=double>
struct NeffW2Density final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
NeffW2Density(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = 1.0;
if (parameters.hasKey("c"))
c_ = parameters.template get<double>("c");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::log;
using std::acosh;
return K + acosh(K) - c_*log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
template<int dim, class field_type, class ctype=double>
struct VossDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = parameters.template get<double>("c");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto KBold = normSquared / (2*det);
using std::acosh;
auto K = KBold + std::sqrt(KBold*KBold - 1);
auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
+ std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
return h + c_ * std::log( det );
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
} // namespace Dune::Elasticity
#endif
src/quasiconvexity-test.cc
+ 37
− 242
View file @ aaab57cf
...@@ -37,16 +37,8 @@ ...@@ -37,16 +37,8 @@
#include <dune/elasticity/common/trustregionsolver.hh> #include <dune/elasticity/common/trustregionsolver.hh>
#include <dune/elasticity/assemblers/localadolcstiffness.hh> #include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh> #include <dune/elasticity/assemblers/feassembler.hh>
//#include <dune/elasticity/materials/neffenergy.hh>
//#include <dune/elasticity/materials/coshenergy.hh>
//#include <dune/elasticity/materials/dacorognaenergy.hh>
//#include <dune/elasticity/materials/expacoshenergy.hh>
//#include <dune/elasticity/materials/fungenergy.hh>
#include <dune/elasticity/materials/localintegralenergy.hh> #include <dune/elasticity/materials/localintegralenergy.hh>
//#include <dune/elasticity/materials/nefflogenergy.hh> #include <dune/elasticity/materials/quasiconvexitytestdensities.hh>
//#include <dune/elasticity/materials/nefflogenergy2.hh>
//#include <dune/elasticity/materials/neffreducedenergy.hh>
//#include <dune/elasticity/materials/neffw2energy.hh>
// grid dimension // grid dimension
const int dim = 2; const int dim = 2;
...@@ -55,164 +47,6 @@ const int order = 1; ...@@ -55,164 +47,6 @@ const int order = 1;
using namespace Dune; using namespace Dune;
template<int dim, class field_type, class ctype=double>
struct BurkholderDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
BurkholderDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
p_ = parameters.template get<double>("p");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto lambdaMaxSquared = 0.5 * (normSquared + std::sqrt(normSquared*normSquared - 4*det*det));
auto lambdaMax = std::sqrt(lambdaMaxSquared);
return -0.5 * p_ * det * std::pow(lambdaMax,p_-2) - 0.5*(2-p_) * std::pow(lambdaMax,p_);
}
FieldMatrix<ctype,dim,dim> F0_;
double p_;
};
template<int dim, class field_type, class ctype=double>
struct ExpACosHDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
ExpACosHDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c1_ = parameters.get<double>("c1");
c2_ = parameters.get<double>("c2");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::acosh;
return exp(c1_ * acosh(K) * acosh(K)) - c2_ * log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c1_;
double c2_;
};
template<int dim, class field_type, class ctype=double>
struct KPowerDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
p_ = parameters.template get<double>("p");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto detF = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*detF);
return std::pow(K, p_);
}
FieldMatrix<ctype,dim,dim> F0_;
double p_;
};
template<int dim, class field_type, class ctype=double>
struct MagicDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = 1.0;
if (parameters.hasKey("c"))
c_ = parameters.template get<double>("c");
}
/** \brief Implement space-dependent coefficients */
double c(const Dune::FieldVector<ctype,dim>& x) const
{
return c_;
#if 0 // Three circles
bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
#endif
// return (x.two_norm() < 0.9) ? 1 : c_;
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto K = normSquared / (2*det);
using std::acosh;
return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
template<int dim, class field_type, class ctype=double>
struct VossDensity final
: public Elasticity::LocalDensity<dim,field_type,ctype>
{
VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
: F0_(F0)
{
c_ = parameters.template get<double>("c");
}
field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
{
auto F = derivative + F0_;
auto det = F.determinant();
auto normSquared = F.frobenius_norm2();
auto KBold = normSquared / (2*det);
using std::acosh;
auto K = KBold + std::sqrt(KBold*KBold - 1);
auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
+ std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
return h + c_ * std::log( det );
}
FieldMatrix<ctype,dim,dim> F0_;
double c_;
};
template <class Basis, class Vector> template <class Basis, class Vector>
void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldMatrix<double,2,2>& F0, std::string filename) void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldMatrix<double,2,2>& F0, std::string filename)
{ {
...@@ -625,92 +459,53 @@ int main (int argc, char *argv[]) try ...@@ -625,92 +459,53 @@ int main (int argc, char *argv[]) try
// Assembler using ADOL-C // Assembler using ADOL-C
std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl; std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
std::shared_ptr<Elasticity::LocalEnergy<FEBasis::LocalView,adouble> > elasticEnergy; std::shared_ptr<Elasticity::LocalDensity<dim,adouble> > density;
if (parameterSet.get<std::string>("energy") == "kpower") if (parameterSet.get<std::string>("energy") == "burkholder")
{ density = std::make_shared<Elasticity::BurkholderDensity<dim,adouble> >(F0,materialParameters);
auto density = std::make_shared<KPowerDensity<dim,adouble> >(F0, materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView, if (parameterSet.get<std::string>("energy") == "cosh")
adouble> >(density); density = std::make_shared<Elasticity::CosHDensity<dim,adouble> >(F0,materialParameters);
}
if (parameterSet.get<std::string>("energy") == "dacorogna")
density = std::make_shared<Elasticity::DacorognaDensity<dim,adouble> >(F0,materialParameters);
// if (parameterSet.get<std::string>("energy") == "neff")
// elasticEnergy = std::make_shared<NeffEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
// if (parameterSet.get<std::string>("energy") == "nefflog")
// elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
// if (parameterSet.get<std::string>("energy") == "nefflog2")
// elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
// if (parameterSet.get<std::string>("energy") == "neffreduced")
// elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
// if (parameterSet.get<std::string>("energy") == "cosh")
// elasticEnergy = std::make_shared<CosHEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
#if 0
if (parameterSet.get<std::string>("energy") == "expacosh")
elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
FEBasis::LocalView::Tree::FiniteElement,
adouble> >(materialParameters);
#else
if (parameterSet.get<std::string>("energy") == "expacosh") if (parameterSet.get<std::string>("energy") == "expacosh")
{ density = std::make_shared<Elasticity::ExpACosHDensity<dim,adouble> >(F0,materialParameters);
auto density = std::make_shared<ExpACosHDensity<dim,adouble> >(F0,materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView, if (parameterSet.get<std::string>("energy") == "fung")
adouble> >(density); density = std::make_shared<Elasticity::FungDensity<dim,adouble> >(F0,materialParameters);
}
#endif if (parameterSet.get<std::string>("energy") == "kpower")
// if (parameterSet.get<std::string>("energy") == "fung") density = std::make_shared<Elasticity::KPowerDensity<dim,adouble> >(F0, materialParameters);
// elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
if (parameterSet.get<std::string>("energy") == "magic") if (parameterSet.get<std::string>("energy") == "magic")
{ density = std::make_shared<Elasticity::MagicDensity<dim,adouble> >(F0,materialParameters);
auto density = std::make_shared<MagicDensity<dim,adouble> >(F0,materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
adouble> >(density);
}
if (parameterSet.get<std::string>("energy") == "burkholder") if (parameterSet.get<std::string>("energy") == "neff")
{ density = std::make_shared<Elasticity::NeffDensity<dim,adouble> >(F0,materialParameters);
auto density = std::make_shared<BurkholderDensity<dim,adouble> >(F0,materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView, if (parameterSet.get<std::string>("energy") == "nefflog")
adouble> >(density); density = std::make_shared<Elasticity::NeffLogDensity<dim,adouble> >(F0,materialParameters);
}
if (parameterSet.get<std::string>("energy") == "nefflog2")
density = std::make_shared<Elasticity::NeffLog2Density<dim,adouble> >(F0,materialParameters);
if (parameterSet.get<std::string>("energy") == "neffreduced")
density = std::make_shared<Elasticity::NeffReducedDensity<dim,adouble> >(F0,materialParameters);
if (parameterSet.get<std::string>("energy") == "w2")
density = std::make_shared<Elasticity::NeffW2Density<dim,adouble> >(F0,materialParameters);
// if (parameterSet.get<std::string>("energy") == "w2")
// elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
// if (parameterSet.get<std::string>("energy") == "dacorogna")
// elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView,
// FEBasis::LocalView::Tree::FiniteElement,
// adouble> >(materialParameters);
//
if (parameterSet.get<std::string>("energy") == "voss") if (parameterSet.get<std::string>("energy") == "voss")
{ auto density = std::make_shared<Elasticity::VossDensity<dim,adouble> >(F0,materialParameters);
auto density = std::make_shared<VossDensity<dim,adouble> >(F0,materialParameters);
elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
adouble> >(density);
}
if(!elasticEnergy) if(!density)
DUNE_THROW(Exception, "Error: Selected energy not available!"); DUNE_THROW(Exception, "Error: Selected energy not available!");
auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<FEBasis::LocalView,
adouble> >(density);
auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy); auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
Elasticity::FEAssembler<FEBasis,SolutionType> assembler(feBasis, localADOLCStiffness); Elasticity::FEAssembler<FEBasis,SolutionType> assembler(feBasis, localADOLCStiffness);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment