Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • enable-linear-elasticity
  • feature/move-to-dune-functions-bases
  • fix/comment
  • fix/hyperdual
  • fix/mooney-rivlin-parameter
  • master
  • mooney-rivlin-zero
  • patrizio-convexity-test
  • relaxed-micromorphic-continuum
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.4-1
  • releases/2.5-1
  • releases/2.6-1
  • releases/2.7
  • releases/2.8
  • subversion->git
19 results

Target

Select target project
  • jonathan.drechsel_at_mailbox.tu-dresden.de/dune-elasticity
  • patrick.jaap_at_tu-dresden.de/dune-elasticity
2 results
Select Git revision
  • feature/autodiff
  • feature/hyperdual-vector-mode
  • finite-strain-plasticity
  • master
  • patrizio-convexity-test
  • releases/2.0-1
  • releases/2.1-1
  • releases/2.2-1
  • releases/2.3-1
  • releases/2.4-1
  • releases/2.5-1
  • releases/2.6-1
  • releases/2.8
  • subversion->git
14 results
Show changes
Commits on Source (73)
Hide whitespace changes
Inline Side-by-side
Showing
with 450 additions and 249 deletions
.gitignore
View file @ 8955c6d5
...@@ -4,3 +4,6 @@ ...@@ -4,3 +4,6 @@
# Default cmake build directory # Default cmake build directory
/build-cmake /build-cmake
# ignore Python files
*.pyc
.gitlab-ci.yml
View file @ 8955c6d5
...@@ -6,6 +6,24 @@ before_script: ...@@ -6,6 +6,24 @@ before_script:
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-fufem.git - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-fufem.git
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
dune:2.8 gcc C++17:
image: registry.dune-project.org/docker/ci/dune:2.8-debian-10-gcc-8-17
script: duneci-standard-test
variables:
DUNECI_BRANCH: releases/2.8
dune:2.8 gcc C++20:
image: registry.dune-project.org/docker/ci/dune:2.8-debian-11-gcc-9-20
script: duneci-standard-test
variables:
DUNECI_BRANCH: releases/2.8
dune:2.8 clang C++17:
image: registry.dune-project.org/docker/ci/dune:2.8-debian-10-clang-7-libcpp-17
script: duneci-standard-test
variables:
DUNECI_BRANCH: releases/2.8
dune:git clang C++17: dune:git clang C++17:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17 image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
script: duneci-standard-test script: duneci-standard-test
......
CHANGELOG.md
View file @ 8955c6d5
# Master (will become release 2.8) # Master (will become release 2.9)
## Deprecations and removals
- The support of `Amiramesh` is dropped in dune-grid and thus the program `linear-elasticity` is removed from CMakeLists.txt
# 2.8 Release
- Introduce class `LocalDensity` for material-specific implementations - Introduce class `LocalDensity` for material-specific implementations
- Introduce class `LocalIntegralEnergy` to work with the densities - Introduce class `LocalIntegralEnergy` to work with the densities
- Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element - Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element
- Introduce class `LocalHyperDualStiffness` and `hyperdual` to calculate gradient and hessian using hyper-dual numbers - Introduce class `LocalHyperDualStiffness` and `hyperdual` to calculate gradient and hessian using hyper-dual numbers
- The class `FEAssembler` now stores a `const` reference to the basis, rather than a value.
This allows to use the assembler with bases that are not copyable
(as, e.g., the `RaviartThomasBasis`, see [dune-functions issue 58](https://gitlab.dune-project.org/staging/dune-functions/-/issues/58)).
## Deprecations and removals ## Deprecations and removals
......
CMakeLists.txt
View file @ 8955c6d5
...@@ -27,7 +27,7 @@ add_subdirectory("dune") ...@@ -27,7 +27,7 @@ add_subdirectory("dune")
add_subdirectory("src") add_subdirectory("src")
add_subdirectory("cmake/modules") add_subdirectory("cmake/modules")
if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG) if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND dune-uggrid_FOUND)
set(programs viscoelast nonlinelast) set(programs viscoelast nonlinelast)
foreach(_program ${programs}) foreach(_program ${programs})
...@@ -40,7 +40,7 @@ if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG) ...@@ -40,7 +40,7 @@ if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG)
endforeach() endforeach()
endif() endif()
if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG AND ADOLC_FOUND) if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND dune-uggrid_FOUND AND ADOLC_FOUND)
add_dune_adolc_flags(nonlinelast) add_dune_adolc_flags(nonlinelast)
endif() endif()
......
dune.module
View file @ 8955c6d5
...@@ -4,8 +4,8 @@ ...@@ -4,8 +4,8 @@
#Name of the module #Name of the module
Module:dune-elasticity Module:dune-elasticity
Version: 2.8-git Version: 2.9-git
Maintainer: oliver.sander@tu-dresden.de, youett@mi.fu-berlin.de Maintainer: oliver.sander@tu-dresden.de, patrick.jaap@tu-dresden.de
#depending on #depending on
Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions
Suggests: dune-uggrid dune-parmg Suggests: dune-uggrid dune-parmg
dune/elasticity/assemblers/feassembler.hh
View file @ 8955c6d5
...@@ -13,6 +13,8 @@ ...@@ -13,6 +13,8 @@
#include <dune/functions/functionspacebases/concepts.hh> #include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/backends/istlvectorbackend.hh> #include <dune/functions/backends/istlvectorbackend.hh>
#include <dune/grid/common/partitionset.hh>
#include "localfestiffness.hh" #include "localfestiffness.hh"
...@@ -31,7 +33,7 @@ class FEAssembler { ...@@ -31,7 +33,7 @@ class FEAssembler {
public: public:
const Basis basis_; const Basis& basis_;
/** \brief Partition type on which to assemble /** \brief Partition type on which to assemble
* *
...@@ -52,6 +54,11 @@ public: ...@@ -52,6 +54,11 @@ public:
localStiffness_(localStiffness) localStiffness_(localStiffness)
{} {}
/** \brief Assemble the functional gradient */
void assembleGradient(const VectorType& configuration,
VectorType& gradient) const;
/** \brief Assemble the tangent stiffness matrix and the functional gradient together /** \brief Assemble the tangent stiffness matrix and the functional gradient together
* *
* This may be more efficient than computing them separately * This may be more efficient than computing them separately
...@@ -69,6 +76,47 @@ public: ...@@ -69,6 +76,47 @@ public:
}; // end class }; // end class
template <class Basis, class VectorType>
void FEAssembler<Basis,VectorType>::
assembleGradient(const VectorType& configuration,
VectorType& gradient) const
{
// create backends for multi index access
auto configurationBackend = Dune::Functions::istlVectorBackend(configuration);
auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
gradientBackend.resize(basis_);
gradient = 0;
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), partitionType))
{
localView.bind(element);
const auto size = localView.size();
// Extract local configuration
std::vector<field_type> localConfiguration(size);
std::vector<field_type> localGradient(size);
for (size_t i=0; i<size; i++)
localConfiguration[i] = configurationBackend[ localView.index(i) ];
// setup local matrix and gradient
localStiffness_->assembleGradient(localView, localConfiguration, localGradient);
for (size_t i=0; i<size; i++)
{
auto row = localView.index(i);
gradientBackend[row] += localGradient[i];
}
}
}
template <class Basis, class VectorType> template <class Basis, class VectorType>
template <class MatrixType> template <class MatrixType>
void FEAssembler<Basis,VectorType>:: void FEAssembler<Basis,VectorType>::
...@@ -117,18 +165,18 @@ assembleGradientAndHessian(const VectorType& configuration, ...@@ -117,18 +165,18 @@ assembleGradientAndHessian(const VectorType& configuration,
Matrix<field_type> localHessian; Matrix<field_type> localHessian;
localHessian.setSize(size, size); localHessian.setSize(size, size);
for (int i=0; i<size; i++) for (size_t i=0; i<size; i++)
localConfiguration[i] = configurationBackend[ localView.index(i) ]; localConfiguration[i] = configurationBackend[ localView.index(i) ];
// setup local matrix and gradient // setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian); localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);
for(int i=0; i<size; i++) for (size_t i=0; i<size; i++)
{ {
auto row = localView.index(i); auto row = localView.index(i);
gradientBackend[row] += localGradient[i]; gradientBackend[row] += localGradient[i];
for (int j=0; j<size; j++ ) for (size_t j=0; j<size; j++ )
{ {
auto col = localView.index(j); auto col = localView.index(j);
// fufem ISTLBackend uses operator() for entry access // fufem ISTLBackend uses operator() for entry access
...@@ -194,8 +242,8 @@ FEAssembler ...@@ -194,8 +242,8 @@ FEAssembler
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock; typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
public: public:
const Basis [[deprecated("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.")]]
DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.") basis_; const Basis basis_;
/** \brief Partition type on which to assemble /** \brief Partition type on which to assemble
* *
......
dune/elasticity/assemblers/localadolcstiffness.hh
View file @ 8955c6d5
...@@ -35,10 +35,20 @@ public: ...@@ -35,10 +35,20 @@ public:
, useHessian2_(useHessian2) , useHessian2_(useHessian2)
{} {}
virtual ~LocalADOLCStiffness() {}
/** \brief Compute the energy at the current configuration */ /** \brief Compute the energy at the current configuration */
virtual double energy (const LocalView& localView, virtual double energy (const LocalView& localView,
const std::vector<double>& localConfiguration) const; const std::vector<double>& localConfiguration) const;
/** \brief Assemble the local gradient at the current position
*
* This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient);
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
* *
* This uses the automatic differentiation toolbox ADOL_C. * This uses the automatic differentiation toolbox ADOL_C.
...@@ -59,7 +69,7 @@ double LocalADOLCStiffness<LocalView>:: ...@@ -59,7 +69,7 @@ double LocalADOLCStiffness<LocalView>::
energy(const LocalView& localView, energy(const LocalView& localView,
const std::vector<double>& localConfiguration) const const std::vector<double>& localConfiguration) const
{ {
int rank = Dune::MPIHelper::getCollectiveCommunication().rank(); int rank = Dune::MPIHelper::getCommunication().rank();
double pureEnergy; double pureEnergy;
std::vector<adouble> localAConfiguration(localConfiguration.size()); std::vector<adouble> localAConfiguration(localConfiguration.size());
...@@ -74,18 +84,43 @@ energy(const LocalView& localView, ...@@ -74,18 +84,43 @@ energy(const LocalView& localView,
try { try {
energy = localEnergy_->energy(localView,localAConfiguration); energy = localEnergy_->energy(localView,localAConfiguration);
} catch (Dune::Exception &e) { } catch (Dune::Exception &e) {
trace_off(rank); trace_off();
throw e; throw;
} }
energy >>= pureEnergy; energy >>= pureEnergy;
trace_off(rank); trace_off();
return pureEnergy; return pureEnergy;
} }
template<class LocalView>
void LocalADOLCStiffness<LocalView>::
assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient
)
{
int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
/////////////////////////////////////////////////////////////////
// Compute the energy gradient
/////////////////////////////////////////////////////////////////
// Compute the actual gradient
size_t nDoubles = localConfiguration.size();
// Compute gradient
localGradient.resize(nDoubles);
gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
}
template<class LocalView> template<class LocalView>
void LocalADOLCStiffness<LocalView>:: void LocalADOLCStiffness<LocalView>::
assembleGradientAndHessian(const LocalView& localView, assembleGradientAndHessian(const LocalView& localView,
...@@ -94,25 +129,20 @@ assembleGradientAndHessian(const LocalView& localView, ...@@ -94,25 +129,20 @@ assembleGradientAndHessian(const LocalView& localView,
Dune::Matrix<double>& localHessian Dune::Matrix<double>& localHessian
) )
{ {
int rank = Dune::MPIHelper::getCollectiveCommunication().rank(); int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute the energy gradient // Compute the energy gradient
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute the actual gradient // this tapes the energy computation
size_t nDoubles = localConfiguration.size(); assembleGradient(localView, localConfiguration, localGradient);
// Compute gradient
localGradient.resize(nDoubles);
gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// Compute Hessian // Compute Hessian
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
size_t nDoubles = localConfiguration.size();
localHessian.setSize(nDoubles,nDoubles); localHessian.setSize(nDoubles,nDoubles);
// allocate raw doubles: ADOL-C requires "**double" arguments // allocate raw doubles: ADOL-C requires "**double" arguments
...@@ -180,7 +210,7 @@ public: ...@@ -180,7 +210,7 @@ public:
LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy) LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
: localEnergy_(energy) : localEnergy_(energy)
{} DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!") {} [[deprecated("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!")]]
/** \brief Compute the energy at the current configuration */ /** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e, virtual RT energy (const Entity& e,
...@@ -208,7 +238,7 @@ energy(const Entity& element, ...@@ -208,7 +238,7 @@ energy(const Entity& element,
const LocalFiniteElement& localFiniteElement, const LocalFiniteElement& localFiniteElement,
const VectorType& localSolution) const const VectorType& localSolution) const
{ {
int rank = Dune::MPIHelper::getCollectiveCommunication().rank(); int rank = Dune::MPIHelper::getCommunication().rank();
double pureEnergy; double pureEnergy;
std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size()); std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size());
...@@ -225,7 +255,7 @@ energy(const Entity& element, ...@@ -225,7 +255,7 @@ energy(const Entity& element,
energy >>= pureEnergy; energy >>= pureEnergy;
trace_off(rank); trace_off();
return pureEnergy; return pureEnergy;
} }
...@@ -244,7 +274,7 @@ assembleGradientAndHessian(const Entity& element, ...@@ -244,7 +274,7 @@ assembleGradientAndHessian(const Entity& element,
const VectorType& localSolution, const VectorType& localSolution,
VectorType& localGradient) VectorType& localGradient)
{ {
int rank = Dune::MPIHelper::getCollectiveCommunication().rank(); int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap. // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(element, localFiniteElement, localSolution); energy(element, localFiniteElement, localSolution);
......
dune/elasticity/assemblers/localfestiffness.hh
View file @ 8955c6d5
...@@ -16,6 +16,13 @@ class LocalFEStiffness ...@@ -16,6 +16,13 @@ class LocalFEStiffness
public: public:
/** \brief Assemble the local gradient at the current position
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<RT>& localConfiguration,
std::vector<RT>& localGradient
) = 0;
/** \brief Assemble the local gradient and tangent matrix at the current position /** \brief Assemble the local gradient and tangent matrix at the current position
*/ */
virtual void assembleGradientAndHessian(const LocalView& localView, virtual void assembleGradientAndHessian(const LocalView& localView,
...@@ -55,8 +62,8 @@ public: ...@@ -55,8 +62,8 @@ public:
VectorType& localGradient) = 0; VectorType& localGradient) = 0;
// assembled data // assembled data
Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > [[deprecated("Use dune-functions powerBases with LocalView concept. See Dune::Elasticity::LocalFEStiffness")]]
DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Dune::Elasticity::LocalFEStiffness") A_; Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
}; };
......
dune/elasticity/assemblers/localhyperdualstiffness.hh
View file @ 8955c6d5
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
namespace Dune::Elasticity { namespace Dune::Elasticity {
/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers /** \brief Assembles energy gradient and Hessian using Hyperdual Numbers
*/ */
template<class LocalView> template<class LocalView>
class LocalHyperDualStiffness class LocalHyperDualStiffness
...@@ -27,6 +27,15 @@ public: ...@@ -27,6 +27,15 @@ public:
virtual double energy (const LocalView& localView, virtual double energy (const LocalView& localView,
const std::vector<double>& localConfiguration) const; const std::vector<double>& localConfiguration) const;
/** \brief Assemble the local gradient at the current position
*
* This uses the automatic differentiation using hyperdual numbers
*/
virtual void assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient);
/** \brief Assemble the local stiffness matrix at the current position /** \brief Assemble the local stiffness matrix at the current position
* *
* This uses the automatic differentiation using hyperdual numbers * This uses the automatic differentiation using hyperdual numbers
...@@ -36,7 +45,7 @@ public: ...@@ -36,7 +45,7 @@ public:
std::vector<double>& localGradient, std::vector<double>& localGradient,
Dune::Matrix<double>& localHessian); Dune::Matrix<double>& localHessian);
std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_; std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_;
}; };
...@@ -50,10 +59,10 @@ energy(const LocalView& localView, ...@@ -50,10 +59,10 @@ energy(const LocalView& localView,
std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size()); std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size());
hyperdual energy; hyperdual energy;
for (size_t i=0; i<localConfiguration.size(); i++) for (size_t i=0; i<localConfiguration.size(); i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]); localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
energy = localEnergy_->energy(localView,localHyperDualConfiguration); energy = localEnergy_->energy(localView,localHyperDualConfiguration);
pureEnergy = energy.real(); pureEnergy = energy.real();
...@@ -62,6 +71,37 @@ energy(const LocalView& localView, ...@@ -62,6 +71,37 @@ energy(const LocalView& localView,
} }
template<class LocalView>
void LocalHyperDualStiffness<LocalView>::
assembleGradient(const LocalView& localView,
const std::vector<double>& localConfiguration,
std::vector<double>& localGradient
)
{
size_t nDoubles = localConfiguration.size();
localGradient.resize(nDoubles);
// Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
for (size_t i=0; i<nDoubles; i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
// Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
// localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
for (size_t i=0; i<nDoubles; i++)
{
localHyperDualConfiguration[i].seteps1(1.0);
hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
localGradient[i] = temp.eps1();
localHyperDualConfiguration[i].seteps1(0.0);
}
}
template<class LocalView> template<class LocalView>
void LocalHyperDualStiffness<LocalView>:: void LocalHyperDualStiffness<LocalView>::
assembleGradientAndHessian(const LocalView& localView, assembleGradientAndHessian(const LocalView& localView,
...@@ -70,7 +110,7 @@ assembleGradientAndHessian(const LocalView& localView, ...@@ -70,7 +110,7 @@ assembleGradientAndHessian(const LocalView& localView,
Dune::Matrix<double>& localHessian Dune::Matrix<double>& localHessian
) )
{ {
size_t nDoubles = localConfiguration.size(); size_t nDoubles = localConfiguration.size();
localGradient.resize(nDoubles); localGradient.resize(nDoubles);
...@@ -79,28 +119,28 @@ assembleGradientAndHessian(const LocalView& localView, ...@@ -79,28 +119,28 @@ assembleGradientAndHessian(const LocalView& localView,
std::vector<hyperdual> localHyperDualConfiguration(nDoubles); std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
for (size_t i=0; i<nDoubles; i++) for (size_t i=0; i<nDoubles; i++)
localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]); localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
// Compute gradient and hessian (symmetry is used) using hyperdual function evaluation // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
// localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
for (size_t i=0; i<nDoubles; i++) for (size_t i=0; i<nDoubles; i++)
{ {
localHyperDualConfiguration[i].seteps1(1.0); localHyperDualConfiguration[i].seteps1(1.0);
localHyperDualConfiguration[i].seteps2(1.0); localHyperDualConfiguration[i].seteps2(1.0);
hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration); hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
localGradient[i] = temp.eps1(); localGradient[i] = temp.eps1();
localHessian[i][i] = temp.eps1eps2(); localHessian[i][i] = temp.eps1eps2();
localHyperDualConfiguration[i].seteps2(0.0); localHyperDualConfiguration[i].seteps2(0.0);
for (size_t j=i+1; j<nDoubles; j++) for (size_t j=i+1; j<nDoubles; j++)
{ {
localHyperDualConfiguration[j].seteps2(1.0); localHyperDualConfiguration[j].seteps2(1.0);
localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2(); localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian
localHyperDualConfiguration[j].seteps2(0.0); localHyperDualConfiguration[j].seteps2(0.0);
} }
localHyperDualConfiguration[i].seteps1(0.0); localHyperDualConfiguration[i].seteps1(0.0);
} }
......
dune/elasticity/common/hyperdual.hh
View file @ 8955c6d5
...@@ -58,63 +58,63 @@ public: ...@@ -58,63 +58,63 @@ public:
//basic manipulation //basic manipulation
hyperdual operator+ () const; hyperdual operator+ () const;
hyperdual operator+ (const hyperdual rhs) const; hyperdual operator+ (const hyperdual& rhs) const;
friend hyperdual operator+ (const double lhs, const hyperdual rhs); friend hyperdual operator+ (const double lhs, const hyperdual& rhs);
hyperdual operator- () const; hyperdual operator- () const;
hyperdual operator- (const hyperdual rhs) const; hyperdual operator- (const hyperdual& rhs) const;
friend hyperdual operator- (const double lhs, const hyperdual rhs); friend hyperdual operator- (const double lhs, const hyperdual& rhs);
hyperdual operator* (const hyperdual rhs)const; hyperdual operator* (const hyperdual& rhs)const;
friend hyperdual operator* (const double lhs, const hyperdual rhs); friend hyperdual operator* (const double lhs, const hyperdual& rhs);
friend hyperdual operator/ (const hyperdual lhs, const hyperdual rhs); friend hyperdual operator/ (const hyperdual& lhs, const hyperdual& rhs);
friend hyperdual operator/ (const double lhs, const hyperdual rhs); friend hyperdual operator/ (const double lhs, const hyperdual& rhs);
friend hyperdual operator/ (const hyperdual lhs, const double rhs); friend hyperdual operator/ (const hyperdual& lhs, const double rhs);
hyperdual& operator+= (hyperdual rhs); hyperdual& operator+= (const hyperdual& rhs);
hyperdual& operator-= (hyperdual rhs); hyperdual& operator-= (const hyperdual& rhs);
hyperdual& operator*= (hyperdual rhs); hyperdual& operator*= (const hyperdual& rhs);
hyperdual& operator*= (double rhs); hyperdual& operator*= (double rhs);
hyperdual& operator/= (hyperdual rhs); hyperdual& operator/= (const hyperdual& rhs);
hyperdual& operator/= (double rhs); hyperdual& operator/= (double rhs);
//math.h functions //math.h functions
friend hyperdual pow (hyperdual x, double a); friend hyperdual pow (const hyperdual& x, double a);
friend hyperdual pow (hyperdual x, hyperdual a); friend hyperdual pow (const hyperdual& x, const hyperdual& a);
friend hyperdual exp(hyperdual x); friend hyperdual exp(const hyperdual& x);
friend hyperdual log(hyperdual x); friend hyperdual log(const hyperdual& x);
friend hyperdual sin(hyperdual x); friend hyperdual sin(const hyperdual& x);
friend hyperdual cos(hyperdual x); friend hyperdual cos(const hyperdual& x);
friend hyperdual tan(hyperdual x); friend hyperdual tan(const hyperdual& x);
friend hyperdual asin(hyperdual x); friend hyperdual asin(const hyperdual& x);
friend hyperdual acos(hyperdual x); friend hyperdual acos(const hyperdual& x);
friend hyperdual atan(hyperdual x); friend hyperdual atan(const hyperdual& x);
friend hyperdual sqrt(hyperdual x); friend hyperdual sqrt(const hyperdual& x);
friend hyperdual fabs(hyperdual x); friend hyperdual fabs(const hyperdual& x);
friend hyperdual abs(hyperdual x); friend hyperdual abs(const hyperdual& x);
friend hyperdual max(hyperdual x1, hyperdual x2); friend hyperdual max(const hyperdual& x1, const hyperdual& x2);
friend hyperdual max(hyperdual x1, double x2); friend hyperdual max(const hyperdual& x1, double x2);
friend hyperdual max(double x1, hyperdual x2); friend hyperdual max(double x1, const hyperdual& x2);
friend hyperdual min(hyperdual x1, hyperdual x2); friend hyperdual min(const hyperdual& x1, const hyperdual& x2);
friend hyperdual min(hyperdual x1, double x2); friend hyperdual min(const hyperdual& x1, double x2);
friend hyperdual min(double x1, hyperdual x2); friend hyperdual min(double x1, const hyperdual& x2);
//comparisons //comparisons
friend bool operator> (hyperdual lhs, hyperdual rhs); friend bool operator> (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator> (double lhs, hyperdual rhs); friend bool operator> (double lhs, const hyperdual& rhs);
friend bool operator> (hyperdual lhs, double rhs); friend bool operator> (const hyperdual& lhs, double rhs);
friend bool operator>= (hyperdual lhs, hyperdual rhs); friend bool operator>= (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator>= (double lhs, hyperdual rhs); friend bool operator>= (double lhs, const hyperdual& rhs);
friend bool operator>= (hyperdual lhs, double rhs); friend bool operator>= (const hyperdual& lhs, double rhs);
friend bool operator< (hyperdual lhs, hyperdual rhs); friend bool operator< (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator< (double lhs, hyperdual rhs); friend bool operator< (double lhs, const hyperdual& rhs);
friend bool operator< (hyperdual lhs, double rhs); friend bool operator< (const hyperdual& lhs, double rhs);
friend bool operator<= (hyperdual lhs, hyperdual rhs); friend bool operator<= (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator<= (double lhs, hyperdual rhs); friend bool operator<= (double lhs, const hyperdual& rhs);
friend bool operator<= (hyperdual lhs, double rhs); friend bool operator<= (const hyperdual& lhs, double rhs);
friend bool operator== (hyperdual lhs, hyperdual rhs); friend bool operator== (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator== (double lhs, hyperdual rhs); friend bool operator== (double lhs, const hyperdual& rhs);
friend bool operator== (hyperdual lhs, double rhs); friend bool operator== (const hyperdual& lhs, double rhs);
friend bool operator!= (hyperdual lhs, hyperdual rhs); friend bool operator!= (const hyperdual& lhs, const hyperdual& rhs);
friend bool operator!= (double lhs, hyperdual rhs); friend bool operator!= (double lhs, const hyperdual& rhs);
friend bool operator!= (hyperdual lhs, double rhs); friend bool operator!= (const hyperdual& lhs, double rhs);
}; };
...@@ -192,7 +192,7 @@ hyperdual hyperdual::operator+ () const ...@@ -192,7 +192,7 @@ hyperdual hyperdual::operator+ () const
{ {
return *this; return *this;
} }
hyperdual hyperdual::operator+ (const hyperdual rhs) const hyperdual hyperdual::operator+ (const hyperdual& rhs) const
{ {
hyperdual temp; hyperdual temp;
temp.f0 = f0 + rhs.f0; temp.f0 = f0 + rhs.f0;
...@@ -201,7 +201,7 @@ hyperdual hyperdual::operator+ (const hyperdual rhs) const ...@@ -201,7 +201,7 @@ hyperdual hyperdual::operator+ (const hyperdual rhs) const
temp.f12 = f12 + rhs.f12; temp.f12 = f12 + rhs.f12;
return temp; return temp;
} }
hyperdual operator+ (const double lhs, const hyperdual rhs) hyperdual operator+ (const double lhs, const hyperdual& rhs)
{ {
hyperdual temp; hyperdual temp;
temp.f0 = lhs + rhs.f0; temp.f0 = lhs + rhs.f0;
...@@ -219,7 +219,7 @@ hyperdual hyperdual::operator- () const ...@@ -219,7 +219,7 @@ hyperdual hyperdual::operator- () const
temp.f12 = -f12; temp.f12 = -f12;
return temp; return temp;
} }
hyperdual hyperdual::operator- (const hyperdual rhs) const hyperdual hyperdual::operator- (const hyperdual& rhs) const
{ {
hyperdual temp; hyperdual temp;
temp.f0 = f0 - rhs.f0; temp.f0 = f0 - rhs.f0;
...@@ -228,7 +228,7 @@ hyperdual hyperdual::operator- (const hyperdual rhs) const ...@@ -228,7 +228,7 @@ hyperdual hyperdual::operator- (const hyperdual rhs) const
temp.f12 = f12 - rhs.f12; temp.f12 = f12 - rhs.f12;
return temp; return temp;
} }
hyperdual operator- (const double lhs, const hyperdual rhs) hyperdual operator- (const double lhs, const hyperdual& rhs)
{ {
hyperdual temp; hyperdual temp;
temp.f0 = lhs - rhs.f0; temp.f0 = lhs - rhs.f0;
...@@ -237,7 +237,7 @@ hyperdual operator- (const double lhs, const hyperdual rhs) ...@@ -237,7 +237,7 @@ hyperdual operator- (const double lhs, const hyperdual rhs)
temp.f12 = -rhs.f12; temp.f12 = -rhs.f12;
return temp; return temp;
} }
hyperdual hyperdual::operator* (const hyperdual rhs) const hyperdual hyperdual::operator* (const hyperdual& rhs) const
{ {
hyperdual temp; hyperdual temp;
temp.f0 = f0*rhs.f0; temp.f0 = f0*rhs.f0;
...@@ -246,7 +246,7 @@ hyperdual hyperdual::operator* (const hyperdual rhs) const ...@@ -246,7 +246,7 @@ hyperdual hyperdual::operator* (const hyperdual rhs) const
temp.f12 = f0*rhs.f12 + f1*rhs.f2 + f2*rhs.f1 + f12*rhs.f0; temp.f12 = f0*rhs.f12 + f1*rhs.f2 + f2*rhs.f1 + f12*rhs.f0;
return temp; return temp;
} }
hyperdual operator* (const double lhs, const hyperdual rhs) hyperdual operator* (const double lhs, const hyperdual& rhs)
{ {
hyperdual temp; hyperdual temp;
temp.f0 = lhs*rhs.f0; temp.f0 = lhs*rhs.f0;
...@@ -255,7 +255,7 @@ hyperdual operator* (const double lhs, const hyperdual rhs) ...@@ -255,7 +255,7 @@ hyperdual operator* (const double lhs, const hyperdual rhs)
temp.f12 = lhs*rhs.f12; temp.f12 = lhs*rhs.f12;
return temp; return temp;
} }
hyperdual operator/ (const hyperdual lhs, const hyperdual rhs) hyperdual operator/ (const hyperdual& lhs, const hyperdual& rhs)
{ {
hyperdual inv; hyperdual inv;
double prod; double prod;
...@@ -267,13 +267,13 @@ hyperdual operator/ (const hyperdual lhs, const hyperdual rhs) ...@@ -267,13 +267,13 @@ hyperdual operator/ (const hyperdual lhs, const hyperdual rhs)
inv.f12 = 2.0*rhs.f1*rhs.f2/(prod*rhs.f0) - rhs.f12/prod; inv.f12 = 2.0*rhs.f1*rhs.f2/(prod*rhs.f0) - rhs.f12/prod;
return lhs*inv; return lhs*inv;
} }
hyperdual operator/ (const double lhs, const hyperdual rhs) hyperdual operator/ (const double lhs, const hyperdual& rhs)
{ {
hyperdual temp,inv; hyperdual temp,inv;
temp = hyperdual(lhs)/rhs; temp = hyperdual(lhs)/rhs;
return temp; return temp;
} }
hyperdual operator/ (const hyperdual lhs, const double rhs) hyperdual operator/ (const hyperdual& lhs, const double rhs)
{ {
hyperdual temp; hyperdual temp;
double inv; double inv;
...@@ -284,7 +284,7 @@ hyperdual operator/ (const hyperdual lhs, const double rhs) ...@@ -284,7 +284,7 @@ hyperdual operator/ (const hyperdual lhs, const double rhs)
temp.f12 = inv*lhs.f12; temp.f12 = inv*lhs.f12;
return temp; return temp;
} }
hyperdual& hyperdual::operator+= (hyperdual rhs) hyperdual& hyperdual::operator+= (const hyperdual& rhs)
{ {
f0 += rhs.f0; f0 += rhs.f0;
f1 += rhs.f1; f1 += rhs.f1;
...@@ -292,7 +292,7 @@ hyperdual& hyperdual::operator+= (hyperdual rhs) ...@@ -292,7 +292,7 @@ hyperdual& hyperdual::operator+= (hyperdual rhs)
f12 += rhs.f12; f12 += rhs.f12;
return *this; return *this;
} }
hyperdual& hyperdual::operator-= (hyperdual rhs) hyperdual& hyperdual::operator-= (const hyperdual& rhs)
{ {
f0 -= rhs.f0; f0 -= rhs.f0;
f1 -= rhs.f1; f1 -= rhs.f1;
...@@ -300,7 +300,7 @@ hyperdual& hyperdual::operator-= (hyperdual rhs) ...@@ -300,7 +300,7 @@ hyperdual& hyperdual::operator-= (hyperdual rhs)
f12 -= rhs.f12; f12 -= rhs.f12;
return *this; return *this;
} }
hyperdual& hyperdual::operator*= (hyperdual rhs) hyperdual& hyperdual::operator*= (const hyperdual& rhs)
{ {
double tf0,tf1,tf2,tf12; double tf0,tf1,tf2,tf12;
tf0=f0; tf0=f0;
...@@ -321,7 +321,7 @@ hyperdual& hyperdual::operator*= (double rhs) ...@@ -321,7 +321,7 @@ hyperdual& hyperdual::operator*= (double rhs)
f12 *= rhs; f12 *= rhs;
return *this; return *this;
} }
hyperdual& hyperdual::operator/= (hyperdual rhs) hyperdual& hyperdual::operator/= (const hyperdual& rhs)
{ {
*this = *this / rhs; *this = *this / rhs;
return *this; return *this;
...@@ -334,7 +334,7 @@ hyperdual& hyperdual::operator/= (double rhs) ...@@ -334,7 +334,7 @@ hyperdual& hyperdual::operator/= (double rhs)
f12 /= rhs; f12 /= rhs;
return *this; return *this;
} }
hyperdual pow (hyperdual x, double a) hyperdual pow (const hyperdual& x, double a)
{ {
hyperdual temp; hyperdual temp;
double deriv,xval,tol; double deriv,xval,tol;
...@@ -356,11 +356,11 @@ hyperdual pow (hyperdual x, double a) ...@@ -356,11 +356,11 @@ hyperdual pow (hyperdual x, double a)
return temp; return temp;
} }
hyperdual pow (hyperdual x, hyperdual a) hyperdual pow (const hyperdual& x, const hyperdual& a)
{ {
return exp(a*log(x)); return exp(a*log(x));
} }
hyperdual exp(hyperdual x) hyperdual exp(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double deriv; double deriv;
...@@ -371,7 +371,7 @@ hyperdual exp(hyperdual x) ...@@ -371,7 +371,7 @@ hyperdual exp(hyperdual x)
temp.f12 = deriv*(x.f12 + x.f1*x.f2); temp.f12 = deriv*(x.f12 + x.f1*x.f2);
return temp; return temp;
} }
hyperdual log(hyperdual x) hyperdual log(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double deriv1,deriv2; double deriv1,deriv2;
...@@ -383,7 +383,7 @@ hyperdual log(hyperdual x) ...@@ -383,7 +383,7 @@ hyperdual log(hyperdual x)
temp.f12 = x.f12/x.f0 - (deriv1*deriv2); temp.f12 = x.f12/x.f0 - (deriv1*deriv2);
return temp; return temp;
} }
hyperdual sin(hyperdual x) hyperdual sin(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv; double funval,deriv;
...@@ -395,7 +395,7 @@ hyperdual sin(hyperdual x) ...@@ -395,7 +395,7 @@ hyperdual sin(hyperdual x)
temp.f12 = deriv*x.f12 - funval*x.f1*x.f2; temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
return temp; return temp;
} }
hyperdual cos(hyperdual x) hyperdual cos(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv; double funval,deriv;
...@@ -407,7 +407,7 @@ hyperdual cos(hyperdual x) ...@@ -407,7 +407,7 @@ hyperdual cos(hyperdual x)
temp.f12 = deriv*x.f12 - funval*x.f1*x.f2; temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
return temp; return temp;
} }
hyperdual tan(hyperdual x) hyperdual tan(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv; double funval,deriv;
...@@ -419,7 +419,7 @@ hyperdual tan(hyperdual x) ...@@ -419,7 +419,7 @@ hyperdual tan(hyperdual x)
temp.f12 = deriv*x.f12 + x.f1*x.f2*(2*funval*deriv); temp.f12 = deriv*x.f12 + x.f1*x.f2*(2*funval*deriv);
return temp; return temp;
} }
hyperdual asin(hyperdual x) hyperdual asin(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv1,deriv; double funval,deriv1,deriv;
...@@ -432,7 +432,7 @@ hyperdual asin(hyperdual x) ...@@ -432,7 +432,7 @@ hyperdual asin(hyperdual x)
temp.f12 = deriv*x.f12 + x.f1*x.f2*(x.f0*std::pow(deriv1,-1.5)); temp.f12 = deriv*x.f12 + x.f1*x.f2*(x.f0*std::pow(deriv1,-1.5));
return temp; return temp;
} }
hyperdual acos(hyperdual x) hyperdual acos(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv1,deriv; double funval,deriv1,deriv;
...@@ -445,7 +445,7 @@ hyperdual acos(hyperdual x) ...@@ -445,7 +445,7 @@ hyperdual acos(hyperdual x)
temp.f12 = deriv*x.f12 + x.f1*x.f2*(-x.f0*std::pow(deriv1,-1.5)); temp.f12 = deriv*x.f12 + x.f1*x.f2*(-x.f0*std::pow(deriv1,-1.5));
return temp; return temp;
} }
hyperdual atan(hyperdual x) hyperdual atan(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
double funval,deriv1,deriv; double funval,deriv1,deriv;
...@@ -458,7 +458,7 @@ hyperdual atan(hyperdual x) ...@@ -458,7 +458,7 @@ hyperdual atan(hyperdual x)
temp.f12 = deriv*x.f12 + x.f1*x.f2*(-2*x.f0/(deriv1*deriv1)); temp.f12 = deriv*x.f12 + x.f1*x.f2*(-2*x.f0/(deriv1*deriv1));
return temp; return temp;
} }
hyperdual sqrt(hyperdual x) hyperdual sqrt(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
temp.f0 = std::sqrt(x.f0); temp.f0 = std::sqrt(x.f0);
...@@ -467,7 +467,7 @@ hyperdual sqrt(hyperdual x) ...@@ -467,7 +467,7 @@ hyperdual sqrt(hyperdual x)
temp.f12 = 0.5*(x.f12 - 2.*temp.f1*temp.f2)/temp.f0; temp.f12 = 0.5*(x.f12 - 2.*temp.f1*temp.f2)/temp.f0;
return temp; return temp;
} }
hyperdual fabs(hyperdual x) hyperdual fabs(const hyperdual& x)
{ {
hyperdual temp; hyperdual temp;
if (x < 0.0) if (x < 0.0)
...@@ -477,11 +477,11 @@ hyperdual fabs(hyperdual x) ...@@ -477,11 +477,11 @@ hyperdual fabs(hyperdual x)
return temp; return temp;
} }
hyperdual abs(hyperdual x) hyperdual abs(const hyperdual& x)
{ {
return fabs(x); return fabs(x);
} }
hyperdual max(hyperdual x1, hyperdual x2) hyperdual max(const hyperdual& x1, const hyperdual& x2)
{ {
hyperdual temp; hyperdual temp;
if (x1>x2) if (x1>x2)
...@@ -490,7 +490,7 @@ hyperdual max(hyperdual x1, hyperdual x2) ...@@ -490,7 +490,7 @@ hyperdual max(hyperdual x1, hyperdual x2)
temp = x2; temp = x2;
return temp; return temp;
} }
hyperdual max(hyperdual x1, double x2) hyperdual max(const hyperdual& x1, double x2)
{ {
hyperdual temp; hyperdual temp;
if (x1>x2) if (x1>x2)
...@@ -499,7 +499,7 @@ hyperdual max(hyperdual x1, double x2) ...@@ -499,7 +499,7 @@ hyperdual max(hyperdual x1, double x2)
temp = x2; temp = x2;
return temp; return temp;
} }
hyperdual max(double x1, hyperdual x2) hyperdual max(double x1, const hyperdual& x2)
{ {
hyperdual temp; hyperdual temp;
if (x1>x2) if (x1>x2)
...@@ -508,7 +508,7 @@ hyperdual max(double x1, hyperdual x2) ...@@ -508,7 +508,7 @@ hyperdual max(double x1, hyperdual x2)
temp = x2; temp = x2;
return temp; return temp;
} }
hyperdual min(hyperdual x1, hyperdual x2) hyperdual min(const hyperdual& x1, const hyperdual& x2)
{ {
hyperdual temp; hyperdual temp;
if (x1<x2) if (x1<x2)
...@@ -517,7 +517,7 @@ hyperdual min(hyperdual x1, hyperdual x2) ...@@ -517,7 +517,7 @@ hyperdual min(hyperdual x1, hyperdual x2)
temp = x2; temp = x2;
return temp; return temp;
} }
hyperdual min(hyperdual x1, double x2) hyperdual min(const hyperdual& x1, double x2)
{ {
hyperdual temp; hyperdual temp;
if (x1<x2) if (x1<x2)
...@@ -526,7 +526,7 @@ hyperdual min(hyperdual x1, double x2) ...@@ -526,7 +526,7 @@ hyperdual min(hyperdual x1, double x2)
temp = x2; temp = x2;
return temp; return temp;
} }
hyperdual min(double x1, hyperdual x2) hyperdual min(double x1, const hyperdual& x2)
{ {
hyperdual temp; hyperdual temp;
if (x1<x2) if (x1<x2)
...@@ -536,79 +536,103 @@ hyperdual min(double x1, hyperdual x2) ...@@ -536,79 +536,103 @@ hyperdual min(double x1, hyperdual x2)
return temp; return temp;
} }
bool operator> (hyperdual lhs, hyperdual rhs) bool operator> (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 > rhs.f0); return (lhs.f0 > rhs.f0);
} }
bool operator> (double lhs, hyperdual rhs) bool operator> (double lhs, const hyperdual& rhs)
{ {
return (lhs > rhs.f0); return (lhs > rhs.f0);
} }
bool operator> (hyperdual lhs, double rhs) bool operator> (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 > rhs); return (lhs.f0 > rhs);
} }
bool operator>= (hyperdual lhs, hyperdual rhs) bool operator>= (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 >= rhs.f0); return (lhs.f0 >= rhs.f0);
} }
bool operator>= (double lhs, hyperdual rhs) bool operator>= (double lhs, const hyperdual& rhs)
{ {
return (lhs >= rhs.f0); return (lhs >= rhs.f0);
} }
bool operator>= (hyperdual lhs, double rhs) bool operator>= (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 >= rhs); return (lhs.f0 >= rhs);
} }
bool operator< (hyperdual lhs, hyperdual rhs) bool operator< (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 < rhs.f0); return (lhs.f0 < rhs.f0);
} }
bool operator< (double lhs, hyperdual rhs) bool operator< (double lhs, const hyperdual& rhs)
{ {
return (lhs < rhs.f0); return (lhs < rhs.f0);
} }
bool operator< (hyperdual lhs, double rhs) bool operator< (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 < rhs); return (lhs.f0 < rhs);
} }
bool operator<= (hyperdual lhs, hyperdual rhs) bool operator<= (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 <= rhs.f0); return (lhs.f0 <= rhs.f0);
} }
bool operator<= (double lhs, hyperdual rhs) bool operator<= (double lhs, const hyperdual& rhs)
{ {
return (lhs <= rhs.f0); return (lhs <= rhs.f0);
} }
bool operator<= (hyperdual lhs, double rhs) bool operator<= (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 <= rhs); return (lhs.f0 <= rhs);
} }
bool operator== (hyperdual lhs, hyperdual rhs) bool operator== (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 == rhs.f0); return (lhs.f0 == rhs.f0);
} }
bool operator== (double lhs, hyperdual rhs) bool operator== (double lhs, const hyperdual& rhs)
{ {
return (lhs == rhs.f0); return (lhs == rhs.f0);
} }
bool operator== (hyperdual lhs, double rhs) bool operator== (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 == rhs); return (lhs.f0 == rhs);
} }
bool operator!= (hyperdual lhs, hyperdual rhs) bool operator!= (const hyperdual& lhs, const hyperdual& rhs)
{ {
return (lhs.f0 != rhs.f0); return (lhs.f0 != rhs.f0);
} }
bool operator!= (double lhs, hyperdual rhs) bool operator!= (double lhs, const hyperdual& rhs)
{ {
return (lhs != rhs.f0); return (lhs != rhs.f0);
} }
bool operator!= (hyperdual lhs, double rhs) bool operator!= (const hyperdual& lhs, double rhs)
{ {
return (lhs.f0 != rhs); return (lhs.f0 != rhs);
} }
bool isnormal( const hyperdual& x )
{
using std::isnormal;
return isnormal(x.real());
}
bool isfinite( const hyperdual& x )
{
using std::isfinite;
return isfinite(x.real());
}
bool isnan( const hyperdual& x )
{
using std::isnan;
return isnan(x.real());
}
bool isinf( const hyperdual& x )
{
using std::isinf;
return isinf(x.real());
}
namespace std namespace std
{ {
......
dune/elasticity/common/trustregionsolver.cc
View file @ 8955c6d5
...@@ -62,7 +62,6 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -62,7 +62,6 @@ setup(const typename BasisType::GridView::Grid& grid,
baseTolerance_ = baseTolerance; baseTolerance_ = baseTolerance;
damping_ = damping; damping_ = damping;
int numLevels = grid_->maxLevel()+1;
const auto dim = VectorType::value_type::dimension; const auto dim = VectorType::value_type::dimension;
#if HAVE_DUNE_PARMG #if HAVE_DUNE_PARMG
...@@ -90,7 +89,7 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -90,7 +89,7 @@ setup(const typename BasisType::GridView::Grid& grid,
std::conditional_t< // do we have a dune-functions basis? std::conditional_t< // do we have a dune-functions basis?
Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(), Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
BasisType, BasisType,
DuneFunctionsBasis<BasisType> > basis(grid.leafGridView()); DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
// //////////////////////////////// // ////////////////////////////////
// Create a multigrid solver // Create a multigrid solver
// //////////////////////////////// // ////////////////////////////////
...@@ -289,10 +288,14 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -289,10 +288,14 @@ setup(const typename BasisType::GridView::Grid& grid,
isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value; isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
} }
int numLevels = grid_->maxLevel()+1;
using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType; using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels); std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
// Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space // Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
// TODO: IIRC that move to P1Basis was only implemented because I didn't have general
// prolongation operators between arbitratry FE spaces.
if (not isP1Basis) if (not isP1Basis)
{ {
TransferOperatorType pkToP1TransferMatrix; TransferOperatorType pkToP1TransferMatrix;
...@@ -325,13 +328,8 @@ setup(const typename BasisType::GridView::Grid& grid, ...@@ -325,13 +328,8 @@ setup(const typename BasisType::GridView::Grid& grid,
for (int i=0; i<numLevels-1; i++) { for (int i=0; i<numLevels-1; i++) {
// Construct the local multigrid transfer matrix // Construct the local multigrid transfer matrix
TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>; transferOperators[i] = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType> >();
newTransferOp->setup(*grid_,i,i+1); transferOperators[i]->setup(*grid_,i,i+1);
auto op = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType>>();
transferOperators[i] = op;
std::shared_ptr<TransferOperatorType> transferOperatorMatrix = std::make_shared<TransferOperatorType>(newTransferOp->getMatrix());
op->setMatrix(transferOperatorMatrix);
} }
mmgStep->setTransferOperators(transferOperators); mmgStep->setTransferOperators(transferOperators);
...@@ -361,14 +359,11 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -361,14 +359,11 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&iterationStep); mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&iterationStep);
} }
BasisType basis(grid_->leafGridView());
#else #else
BasisType basis(grid_->levelGridView(grid_->maxLevel()));
BasisType coarseBasis(grid_->levelGridView(0)); BasisType coarseBasis(grid_->levelGridView(0));
std::vector<BoxConstraint<typename VectorType::field_type, blocksize>> coarseTrustRegionObstacles(coarseBasis.size()); std::vector<BoxConstraint<typename VectorType::field_type, blocksize>> coarseTrustRegionObstacles(coarseBasis.size());
#endif #endif
MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_); MaxNormTrustRegion<blocksize> trustRegion(assembler_->basis_.size(), initialTrustRegionRadius_);
std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles; std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles;
...@@ -556,16 +551,12 @@ void TrustRegionSolver<BasisType,VectorType>::solve() ...@@ -556,16 +551,12 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
if (correctionInfinityNorm < this->tolerance_) { if (correctionInfinityNorm < this->tolerance_) {
if (this->verbosity_ == NumProc::FULL and rank==0) if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl; {
if (correctionInfinityNorm < trustRegion.radius())
if (this->verbosity_ != NumProc::QUIET and rank==0) std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
std::cout << i+1 << " trust-region steps were taken." << std::endl; else
break; std::cout << "TRUST-REGION UNDERFLOW!" << std::endl;
} }
if (trustRegion.radius() < this->tolerance_) {
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "TRUST REGION RADIUS IS TOO SMALL, SMALLER THAN TOLERANCE, STOPPING NOW" << std::endl;
if (this->verbosity_ != NumProc::QUIET and rank==0) if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl; std::cout << i+1 << " trust-region steps were taken." << std::endl;
......
dune/elasticity/common/trustregionsolver.hh
View file @ 8955c6d5
...@@ -26,6 +26,11 @@ ...@@ -26,6 +26,11 @@
#include <dune/elasticity/assemblers/feassembler.hh> #include <dune/elasticity/assemblers/feassembler.hh>
#if HAVE_DUNE_PARMG #if HAVE_DUNE_PARMG
// enable a temporary hack to let parmg work with powerBases
// which is needed for the current implementation of elasticty
#define DUNE_PARMG_POWER_BASIS_HACK
#include <dune/parmg/iterationstep/lambdastep.hh> #include <dune/parmg/iterationstep/lambdastep.hh>
#include <dune/parmg/iterationstep/multigrid.hh> #include <dune/parmg/iterationstep/multigrid.hh>
#include <dune/parmg/parallel/dofmap.hh> #include <dune/parmg/parallel/dofmap.hh>
...@@ -68,14 +73,7 @@ public: ...@@ -68,14 +73,7 @@ public:
TrustRegionSolver() TrustRegionSolver()
: IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL), : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(std::shared_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL) hessianMatrix_(std::shared_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
{ {}
#if HAVE_DUNE_PARMG
// TODO currently, we cannot use PARMG together with dune-functions bases
static_assert( not Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() ,
"Currently, dune-parmg and dune-functions bases cannot be used together."
);
#endif
}
/** \brief Set up the solver using a monotone multigrid method as the inner solver */ /** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid, void setup(const GridType& grid,
......
dune/elasticity/materials/exphenckydensity.hh
View file @ 8955c6d5
...@@ -85,7 +85,7 @@ public: ...@@ -85,7 +85,7 @@ public:
} }
/** \brief Lame constants */ /** \brief Lame constants */
field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_; double mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
}; };
} // namespace Dune::Elasticity } // namespace Dune::Elasticity
......
dune/elasticity/materials/exphenckyenergy.hh
View file @ 8955c6d5
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
namespace Dune { namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double> template<class GridView, class LocalFiniteElement, class field_type=double>
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity") class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity")]]
ExpHenckyEnergy ExpHenckyEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type> : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{ {
......
dune/elasticity/materials/henckyenergy.hh
View file @ 8955c6d5
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
namespace Dune { namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double> template<class GridView, class LocalFiniteElement, class field_type=double>
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity") class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity")]]
HenckyEnergy HenckyEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type> : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{ {
......
dune/elasticity/materials/localintegralenergy.hh
View file @ 8955c6d5
...@@ -2,6 +2,10 @@ ...@@ -2,6 +2,10 @@
#define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH #define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 8)
#include <dune/common/transpose.hh>
#endif
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
...@@ -53,8 +57,9 @@ energy(const LocalView& localView, ...@@ -53,8 +57,9 @@ energy(const LocalView& localView,
field_type energy = 0; field_type energy = 0;
// store gradients of shape functions and base functions // store gradients of shape functions and base functions
std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size()); using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size()); using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
std::vector<Jacobian> jacobians(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order() int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim; : localFiniteElement.localBasis().order() * gridDim;
...@@ -65,23 +70,30 @@ energy(const LocalView& localView, ...@@ -65,23 +70,30 @@ energy(const LocalView& localView,
{ {
const DT integrationElement = element.geometry().integrationElement(qp.position()); const DT integrationElement = element.geometry().integrationElement(qp.position());
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position()); const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
// Global position // Global position
auto x = element.geometry().global(qp.position()); auto x = element.geometry().global(qp.position());
// Get gradients of shape functions // Get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients); localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
// compute gradients of base functions // compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i) for (size_t i=0; i<jacobians.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 8)
jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
#else
{
auto referenceJacobian = jacobians[i];
geometryJacobianIT.mv(referenceJacobian[0], jacobians[i][0]);
}
#endif
// Deformation gradient // Deformation gradient
FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0); FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
for (size_t i=0; i<gradients.size(); i++) for (size_t i=0; i<jacobians.size(); i++)
for (int j=0; j<gridDim; j++) for (int j=0; j<gridDim; j++)
deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]); deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
// Integrate the energy density // Integrate the energy density
energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient); energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
} }
...@@ -124,8 +136,8 @@ public: ...@@ -124,8 +136,8 @@ public:
const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const; const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
protected: protected:
const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>> [[deprecated("Use dune-functions powerBases with LocalView concept. See Elasticity::LocalIntegralEnergy")]]
DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::LocalIntegralEnergy") localDensity_ = nullptr; const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
}; };
......
dune/elasticity/materials/mooneyrivlindensity.hh
View file @ 8955c6d5
...@@ -2,7 +2,6 @@ ...@@ -2,7 +2,6 @@
#define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH #define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/elasticity/materials/localdensity.hh> #include <dune/elasticity/materials/localdensity.hh>
...@@ -25,7 +24,9 @@ public: ...@@ -25,7 +24,9 @@ public:
mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a"); mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b"); mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c"); mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
} else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square") { }
else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square")
{
mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10"); mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01"); mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20"); mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
...@@ -35,8 +36,10 @@ public: ...@@ -35,8 +36,10 @@ public:
mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03"); mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21"); mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12"); mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k"); mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
} else { }
else
{
DUNE_THROW(Exception, "Error: Selected mooneyrivlin implementation (" << mooneyrivlin_energy << ") not available!"); DUNE_THROW(Exception, "Error: Selected mooneyrivlin implementation (" << mooneyrivlin_energy << ") not available!");
} }
} }
...@@ -58,78 +61,96 @@ public: ...@@ -58,78 +61,96 @@ public:
for (int k=0; k<dim; k++) for (int k=0; k<dim; k++)
C[i][j] += gradient[k][i] * gradient[k][j]; C[i][j] += gradient[k][i] * gradient[k][j];
////////////////////////////////////////////////////////// /* The Mooney-Rivlin-Density is given as a function of the eigenvalues of the right Cauchy-Green-Deformation tensor C = F^TF
// Eigenvalues of C = F^TF or the right Cauchy-Green-Deformation tensor B = FF^T.
////////////////////////////////////////////////////////// C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
There are three Mooney-Rivlin-Variants:
// eigenvalues of C, i.e., squared singular values \sigma_i of F ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF^2 + mooneyrivlin_c*(detF)^2 -
Dune::FieldVector<field_type, dim> sigmaSquared; (2*mooneyrivlin_a + 4*mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF),
FMatrixHelp::eigenValues(C, sigmaSquared); where the last term is chosen s.t. W( t*I ) is minimal for t=1
log: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2
field_type normFSquared = gradient.frobenius_norm2(); square: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2
For log and square: I1 and I2 are the first two invariants of C (or B), multiplied with a factor depending on det(F):
I1 = (det(F)^(-2/dim)) * [ first invariant of C ]
= (det(F)^(-2/dim)) * (sum of all eigenvalues of C)
= (det(F)^(-2/dim)) * trace(C)
= (det(F)^(-2/dim)) * trace(F^TF)
= (det(F)^(-2/dim)) * (frobenius_norm(F))^2
I2 = (det(F)^(-4/dim)) * [ second invariant of C ]
= (det(F)^(-4/dim)) * 0.5 * [ trace(C)^2 - tr(C^2) ]
= (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12C_21 - C_13C_31 - C_23C_32]
= (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12^2 - C_13^2 - C_23^2]
*/
field_type frobeniusNormFsquared = gradient.frobenius_norm2();
field_type detF = gradient.determinant(); field_type detF = gradient.determinant();
if (detF < 0) if (detF < 0)
DUNE_THROW( Dune::Exception, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy."); DUNE_THROW( Dune::MathError, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
field_type normFinvSquared = 0;
field_type c2Tilde = 0;
for (int i = 0; i < dim; i++) {
normFinvSquared += 1/sigmaSquared[i];
// compute D, which is the sum of the squared eigenvalues
for (int j = i+1; j < dim; j++)
c2Tilde += sigmaSquared[j]*sigmaSquared[i];
}
field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
// \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
c2Tilde /= pow(detF, 4.0/dim);
field_type c2TildeMinus3 = c2Tilde - 3;
// Add the local energy density // Add the local energy density
field_type strainEnergy = 0; field_type strainEnergy = 0;
if (mooneyrivlin_energy == "ciarlet") if (mooneyrivlin_energy == "ciarlet")
{ {
//To calculate the squared frobeniusnorm of F^(-1), we actually invert F^(-1)
auto gradientInverse = gradient;
gradientInverse.invert();
field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
using std::log; using std::log;
return mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF); return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF
} else { - (2.0*mooneyrivlin_a + 4.0*mooneyrivlin_b + 2.0*mooneyrivlin_c)*log(detF);
strainEnergy = mooneyrivlin_10 * trCTildeMinus3 + }
mooneyrivlin_01 * c2TildeMinus3 + else
mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 + {
mooneyrivlin_02 * c2TildeMinus3 * c2TildeMinus3 + // mooneyrivlin_energy is "log" or "square"
mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3 + field_type a = pow(detF, 2.0/dim);
mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 + field_type invariant1Minus3 = frobeniusNormFsquared/a - 3;
mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3 + field_type secondInvariantOfC = 0;
mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3 * c2TildeMinus3 + for (int i=0; i<dim-1; i++)
mooneyrivlin_03 * c2TildeMinus3 * c2TildeMinus3 * c2TildeMinus3; for (int j=i+1; j<dim; j++)
if (mooneyrivlin_energy == "log") { secondInvariantOfC += C[i][i]*C[j][j] - C[i][j]*C[i][j];
field_type invariant2Minus3 = secondInvariantOfC/(a*a) - 3;
strainEnergy = mooneyrivlin_10 * invariant1Minus3 +
mooneyrivlin_01 * invariant2Minus3 +
mooneyrivlin_20 * invariant1Minus3 * invariant1Minus3 +
mooneyrivlin_02 * invariant2Minus3 * invariant2Minus3 +
mooneyrivlin_11 * invariant1Minus3 * invariant2Minus3 +
mooneyrivlin_30 * invariant1Minus3 * invariant1Minus3 * invariant1Minus3 +
mooneyrivlin_21 * invariant1Minus3 * invariant1Minus3 * invariant2Minus3 +
mooneyrivlin_12 * invariant1Minus3 * invariant2Minus3 * invariant2Minus3 +
mooneyrivlin_03 * invariant2Minus3 * invariant2Minus3 * invariant2Minus3;
if (mooneyrivlin_energy == "log")
{
using std::log; using std::log;
field_type logDetF = log(detF); field_type logDetF = log(detF);
return strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF; return strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF;
} else { //mooneyrivlin_energy is "square" }
else
{
//mooneyrivlin_energy is "square"
field_type detFMinus1 = detF - 1; field_type detFMinus1 = detF - 1;
return strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1; return strainEnergy + 0.5 * mooneyrivlin_k* detFMinus1 * detFMinus1;
} }
} }
} }
/** \brief Lame constants */ /** \brief Lame constants */
field_type mooneyrivlin_a, double mooneyrivlin_a,
mooneyrivlin_b, mooneyrivlin_b,
mooneyrivlin_c, mooneyrivlin_c,
mooneyrivlin_10, mooneyrivlin_10,
mooneyrivlin_01, mooneyrivlin_01,
mooneyrivlin_20, mooneyrivlin_20,
mooneyrivlin_02, mooneyrivlin_02,
mooneyrivlin_11, mooneyrivlin_11,
mooneyrivlin_30, mooneyrivlin_30,
mooneyrivlin_21, mooneyrivlin_21,
mooneyrivlin_12, mooneyrivlin_12,
mooneyrivlin_03, mooneyrivlin_03,
mooneyrivlin_k; mooneyrivlin_k;
std::string mooneyrivlin_energy; std::string mooneyrivlin_energy;
}; };
......
dune/elasticity/materials/mooneyrivlinenergy.hh
View file @ 8955c6d5
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
namespace Dune { namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double> template<class GridView, class LocalFiniteElement, class field_type=double>
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::MooneyRivlinDensity") class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::MooneyRivlinDensity")]]
MooneyRivlinEnergy MooneyRivlinEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type> : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{ {
......
dune/elasticity/materials/neohookedensity.hh
View file @ 8955c6d5
...@@ -76,7 +76,7 @@ public: ...@@ -76,7 +76,7 @@ public:
} }
/** \brief Lame constants */ /** \brief Lame constants */
field_type mu_, lambda_, kappa_; double mu_, lambda_, kappa_;
}; };
} // namespace Dune::Elasticity } // namespace Dune::Elasticity
......
dune/elasticity/materials/neohookeenergy.hh
View file @ 8955c6d5
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
namespace Dune { namespace Dune {
template<class GridView, class LocalFiniteElement, class field_type=double> template<class GridView, class LocalFiniteElement, class field_type=double>
class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::NeoHookeDensity") class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::NeoHookeDensity")]]
NeoHookeEnergy NeoHookeEnergy
: public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type> : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
{ {
......