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

Target

Select target project
  • osander/dune-gfe
  • lnebel/dune-gfe
  • spraetor/dune-gfe
3 results
Show changes
Commits on Source (109)
9 additional commits have been omitted to prevent performance issues.
Hide whitespace changes
Inline Side-by-side
Showing
with 1650 additions and 1835 deletions
.git-blame-ignore-revs
View file @ 48b9358e
# First use of uncrustify, with the vanilla config from dune-project.org
d13f34469065230dd86e63733e4dc8d1f775acfb
# Cleanup after having moved everything into the Dune::GFE namespace
c2eb8dacb5885363c30c2b7373c479ed9c66ac12
.gitlab-ci.yml
View file @ 48b9358e
......@@ -7,9 +7,10 @@ before_script: &before
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-vtk.git
- duneci-install-module https://gitlab.dune-project.org/fufem/dune-fufem.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/iwr/dune-gmsh4
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-gmsh4
# Tests with the 2.9 release. That's the one in Debian bookworm
#--------------------------------------------------------------------
dune:2.9 gcc:
variables:
DUNECI_BRANCH: releases/2.9
......@@ -18,6 +19,76 @@ dune:2.9 gcc:
- *before
script: duneci-standard-test
dune:2.9 dune-elasticity gcc:
variables:
DUNECI_BRANCH: releases/2.9
image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Tests with the 2.10 release. That's the one in Debian trixie
#--------------------------------------------------------------------
dune:2.10 gcc-10:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
script: duneci-standard-test
dune:2.10 dune-elasticity gcc-10:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:2.10 dune-curvedgrid dune-foamgrid gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Full tests again, but this time with clang.
dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-clang-11-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Tests with the git master branch
# Some, but not all optional dependencies are made available.
#------------------------------------------------------------------
dune:git gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
......@@ -30,24 +101,44 @@ dune:git clang-11 C++20:
- *before
script: duneci-standard-test
dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid gcc-10 C++20:
dune:git dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:git dune-curvedgrid dune-foamgrid gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
# Tests with all optional dependencies
#-------------------------------------------
dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid clang-11 C++20:
dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-11-20
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Check for spelling mistakes in text
......
CHANGELOG.md
View file @ 48b9358e
# Master
- The `SimoFoxEnergy` class (formerly known as `SimoFoxEnergyLocalStiffness`)
does not accept a surface load density anymore. Use the
`NeumannEnergy` class for assembling surface loads.
- The class `LocalIntegralEnergy` now does not construct the local interpolation
rule itself anymore. Rather, it expects to be given one in its constructor.
This will allow interpolation rules with state.
- The same holds for `CosseratRodEnergy`. Unfortunately, that move
implies that `CosseratRodEnergy` gains an extra template argument
which specifies the GFE function used to implement the reference
configuration.
- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
now take a `dune-functions` basis as their first arguments, instead
of a `LocalFiniteElement`. In addition, they need to be given an object
of this type in the constructor. This is the basis to be used for
interpolation weights (for `LocalGeodesicFEFunction`) and for
interpolating in the embedding space (for `LocalProjectedFEFunction`).
- The entire code is now in the `Dune::GFE` namespace.
- Remove the files `localgfetestfunctionbasis.hh` and `localtangentfefunction.hh`.
They were never used for anything at all, I and currently wouldn't know
what to use them for, either.
- All implementations of functions have been move to a `functions` subdirectory.
- The `LocalDensity` class now has a template parameter `ElementOrIntersection`,
which replaces the parameter `Position`. As the name says, this parameter
has to be a grid element, or a grid intersection. As it turned out,
just depending on the position type of the integration domain was not enough.
- The `LocalDensity` class now has a `bind` method, which binds the density
to a given element or intersection. This is necessary, for example, to in turn
bind coefficient functions that the density may own.
# Release 2.10 (2024-12-29)
- The module `dune-elasticity` has been downgraded from a required
dependency to an optional dependency. It is currently only needed
by the `film-on-substrate.cc` program.
- `cosserat-rod.cc` (formerly `rod3d.cc`) now reads the reference configuration
and the initial iterate from the Python file.
- The file `rod3d.cc` has been renamed to `cosserat-rod.cc`,
to better reflect what it does.
- In the `BulkCosseratDensity` class: Make the type of curvature
tensor controllable at run-time (previously, ugly preprocessor switches
where used). Among its parameters `BulkCosseratDensity` now expects
a string `curvatureType`, which has to take one of the values `norm`,
`curl`, or `wryness`.
- Building `dune-gfe` now requires Dune 2.9 or newer. Older configurations
have not gotten CI-tested for quite a while, anyway.
......
CMakeLists.txt
View file @ 48b9358e
......@@ -17,6 +17,13 @@ include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
# Create a module library
if(dune-common_VERSION VERSION_GREATER_EQUAL 2.9.0)
dune_add_library(dunegfe INTERFACE EXPORT_NAME GFE NAMESPACE Dune::
LINK_LIBRARIES ${DUNE_LIBS})
endif()
dune_enable_all_packages()
add_subdirectory("src")
add_subdirectory("dune")
......
dune.module
View file @ 48b9358e
......@@ -4,8 +4,8 @@
#Name of the module
Module: dune-gfe
Version: svn
Version: 2.11-git
Maintainer: oliver.sander@tu-dresden.de
#depending on
Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-elasticity (>=2.9) dune-gmsh4
Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgrid
Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-gmsh4 dune-vtk
Suggests: dune-foamgrid dune-parmg dune-curvedgrid dune-elasticity
dune/gfe/CMakeLists.txt
View file @ 48b9358e
......@@ -4,20 +4,12 @@ install(FILES
cosseratstrain.hh
cosseratvtkreader.hh
cosseratvtkwriter.hh
embeddedglobalgfefunction.hh
filereader.hh
geodesicdifference.hh
geodesicfefunctionadaptor.hh
gfedifferencenormsquared.hh
globalgfefunction.hh
gramschmidtsolver.hh
interpolationderivatives.hh
linearalgebra.hh
localgeodesicfefunction.hh
localgfetestfunctionbasis.hh
localprojectedfefunction.hh
localquickanddirtyfefunction.hh
localtangentfefunction.hh
maxnormtrustregion.hh
mixedriemanniantrsolver.cc
mixedriemanniantrsolver.hh
......@@ -42,3 +34,4 @@ install(FILES
add_subdirectory("spaces")
add_subdirectory("assemblers")
add_subdirectory("densities")
add_subdirectory("functions")
dune/gfe/assemblers/CMakeLists.txt
View file @ 48b9358e
#install headers
install(FILES
cosseratenergystiffness.hh
cosseratrodenergy.hh
geodesicfeassembler.hh
geodesicfeassemblerwrapper.hh
......
dune/gfe/assemblers/cosseratenergystiffness.hh deleted 100644 → 0
View file @ 5db924c4
#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
#define COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/version.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#ifdef PROJECTED_INTERPOLATION
#include <dune/gfe/localprojectedfefunction.hh>
#else
#include <dune/gfe/localgeodesicfefunction.hh>
#endif
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/cosseratstrain.hh>
#include <dune/gfe/spaces/orthogonalmatrix.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
#define DONT_USE_CURL
#define CURVATURE_WITH_WRYNESS //only relevant for gridDim == 3
//#define QUADRATIC_2006 //only relevant for gridDim == 3
// Use the formula for the quadratic membrane energy (first equation of (4.4)) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
// instead of the equation (2.27) of 2019: Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature
//#define QUADRATIC_MEMBRANE_ENERGY //only relevant for gridDim == 2
/** \brief Get LocalFiniteElements from a localView, for different tree depths of the local view
*
* We instantiate the CosseratEnergyLocalStiffness class with two different kinds of Basis:
* A scalar one and a composite one that combines two scalar ones. But code for accessing the
* finite elements in the basis tree only work for one kind of basis, not for the other.
* To allow both kinds of basis in a single class we need this trickery below.
*/
template <class Basis, std::size_t i>
class LocalFiniteElementFactory
{
public:
static auto get(const typename Basis::LocalView& localView,
std::integral_constant<std::size_t, i> iType)
-> decltype(localView.tree().child(iType,0).finiteElement())
{
return localView.tree().child(iType,0).finiteElement();
}
};
/** \brief Specialize for scalar bases, here we cannot call tree().child() */
template <class GridView, int order, std::size_t i>
class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order>,i>
{
public:
static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView,
std::integral_constant<std::size_t, i> iType)
-> decltype(localView.tree().finiteElement())
{
return localView.tree().finiteElement();
}
};
template<class Basis, int dim, class field_type=double>
class CosseratEnergyLocalStiffness
: public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
constexpr static int gridDim = GridView::dimension;
constexpr static int dimworld = GridView::dimensionworld;
/** \brief Compute the (row-wise) curl of a matrix R \f$
\param DR The partial derivatives of the matrix R
*/
static Dune::FieldMatrix<field_type,dim,dim> curl(const Tensor3<field_type,dim,dim,dim>& DR)
{
Dune::FieldMatrix<field_type,dim,dim> result;
for (int i=0; i<dim; i++) {
result[i][0] = DR[i][2][1] - DR[i][1][2];
result[i][1] = DR[i][0][2] - DR[i][2][0];
result[i][2] = DR[i][1][0] - DR[i][0][1];
}
return result;
}
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
const BoundaryPatch<GridView>* neumannBoundary,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{
// The shell thickness // only relevant for dim == 2
thickness_ = parameters.template get<double>("thickness");
// Lame constants
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
// Cosserat couple modulus
mu_c_ = parameters.template get<double>("mu_c");
// Length scale parameter
L_c_ = parameters.template get<double>("L_c");
// Curvature exponent
q_ = parameters.template get<double>("q");
// Shear correction factor // only relevant for dim == 2
kappa_ = parameters.template get<double>("kappa");
// Curvature parameters
b1_ = parameters.template get<double>("b1");
b2_ = parameters.template get<double>("b2");
b3_ = parameters.template get<double>("b3");
}
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override;
RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the first equation of (4.4) in Neff's paper from 2006: A geometrically exact planar Cosserat shell model with microstructure: Existence of minimizers for zero Cosserat couple modulus
* OR: the equation (2.27) of 2019: Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature
*/
RT quadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
#ifdef QUADRATIC_2006
+ (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
#else
+ lambda_/2 * Dune::GFE::traceSquared(UMinus1); // Dune::GFE::traceSquared(UMinus1) = Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1))
#endif
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the second equation of (4.4) in Neff's paper
*/
RT longQuadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
RT result = 0;
// shear-stretch energy
Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
result += mu_ * sym2x2.frobenius_norm2();
// first order drill energy
Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
result += mu_c_ * skew2x2.frobenius_norm2();
// classical transverse shear energy
result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
// elongational stretch energy
result += mu_*lambda_ / (2*mu_ + lambda_) * traceSquared(sym2x2);
return result;
}
/** \brief Energy for large-deformation problems (private communication by Patrizio Neff)
*/
RT nonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
Dune::FieldMatrix<field_type,3,3> UMinus1 = U.matrix();
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
RT detU = U.determinant();
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2() + mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
* the second equation of (4.4) in Neff's paper
*/
RT longNonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
{
RT result = 0;
// shear-stretch energy
Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
result += mu_ * sym2x2.frobenius_norm2();
// first order drill energy
Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
for (int i=0; i<dim-1; i++)
for (int j=0; j<dim-1; j++)
skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
result += mu_c_ * skew2x2.frobenius_norm2();
// classical transverse shear energy
result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
// elongational stretch energy
RT detU = U.determinant();
result += (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
return result;
}
RT curvatureEnergy(const Tensor3<field_type,3,3,gridDim>& DR) const
{
using std::pow;
#ifdef DONT_USE_CURL
return mu_ * pow(L_c_ * L_c_ * DR.frobenius_norm2(),q_/2.0);
#else
return mu_ * pow(L_c_ * L_c_ * curl(DR).frobenius_norm2(),q_/2.0);
#endif
}
RT curvatureWithWryness(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
{
// construct Wryness tensor \Gamma as in "Refined dimensional reduction for isotropic elastic Cosserat shells with initial curvature"
Dune::FieldMatrix<field_type,3,3> dRx1(0); //Derivative of R with respect to x1
Dune::FieldMatrix<field_type,3,3> dRx2(0); //Derivative of R with respect to x2
Dune::FieldMatrix<field_type,3,3> dRx3(0); //Derivative of R with respect to x3, this of course only exists if gridDim = 3 (then dim = 3 also, because dim >= gridDim, always)
for (int i=0; i<3; i++)
for (int j=0; j<3; j++) {
dRx1[i][j] = DR[i][j][0];
dRx2[i][j] = DR[i][j][1];
dRx3[i][j] = gridDim == 3 ? DR[i][j][2] : 0;
}
Dune::FieldMatrix<field_type,3,3> wryness(0);
auto axialVectorx1 = SkewMatrix<field_type,3>(transpose(R)*dRx1).axial();
auto axialVectorx2 = SkewMatrix<field_type,3>(transpose(R)*dRx2).axial();
auto axialVectorx3 = SkewMatrix<field_type,3>(transpose(R)*dRx3).axial();
for (int i=0; i<3; i++) {
wryness[i][0] = axialVectorx1[i];
wryness[i][1] = axialVectorx2[i];
wryness[i][2] = gridDim == 3 ? axialVectorx3[i] : 0;
}
// The choice for the Frobenius norm here is b1=b2=1 and b3 = 1/3
return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(wryness)).frobenius_norm2()
+ b2_ * Dune::GFE::skew(wryness).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(wryness));
}
RT bendingEnergy(const Dune::FieldMatrix<field_type,dim,dim>& R, const Tensor3<field_type,3,3,gridDim>& DR) const
{
// left-multiply the derivative of the third director (in DR[][2][]) with R^T
Dune::FieldMatrix<field_type,3,3> RT_DR3(0);
for (int i=0; i<3; i++)
for (int j=0; j<gridDim; j++)
for (int k=0; k<3; k++)
RT_DR3[i][j] += R[k][i] * DR[k][2][j];
return mu_ * Dune::GFE::sym(RT_DR3).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(RT_DR3).frobenius_norm2()
+ mu_*lambda_/(2*mu_+lambda_) * Dune::GFE::traceSquared(RT_DR3);
}
/** \brief The shell thickness */
double thickness_;
/** \brief Lame constants */
double mu_, lambda_;
/** \brief Cosserat couple modulus, preferably 0 */
double mu_c_;
/** \brief Length scale parameter */
double L_c_;
/** \brief Curvature exponent */
double q_;
/** \brief Shear correction factor */
double kappa_;
/** \brief Curvature parameters */
double b1_, b2_, b3_;
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_;
/** \brief The function implementing the Neumann data */
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
/** \brief The function implementing a volume load */
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
};
template <class Basis, int dim, class field_type>
[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> >>.")]]
typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
CosseratEnergyLocalStiffness<Basis,dim,field_type>::
energy(const typename Basis::LocalView& localView,
const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
{
using namespace Dune::Indices;
RT energy = 0;
auto element = localView.element();
using namespace Dune::Indices;
const auto& localFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
#ifdef PROJECTED_INTERPOLATION
typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
#else
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
#endif
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
const DT integrationElement = element.geometry().integrationElement(quadPos);
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
DT weight = quad[pt].weight() * integrationElement;
// The value of the local function
Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function defined on the reference element
typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
// The derivative of the function defined on the actual element
typename LocalGFEFunctionType::DerivativeType derivative(0);
for (size_t comp=0; comp<referenceDerivative.N(); comp++)
jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
/////////////////////////////////////////////////////////
// compute U, the Cosserat strain
/////////////////////////////////////////////////////////
static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
//
Dune::FieldMatrix<field_type,dim,dim> R;
value[_1].matrix(R);
Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(derivative,R);
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
// Extract the orientation derivative
Dune::FieldMatrix<field_type,4,gridDim> orientationDerivative;
for (size_t i=0; i<4; ++i)
for (size_t j=0; j<gridDim; ++j)
orientationDerivative[i][j] = derivative[3+i][j];
Tensor3<field_type,3,3,gridDim> DR = value[_1].quaternionTangentToMatrixTangent(orientationDerivative);
// Add the local energy density
if (gridDim==2) {
#ifdef QUADRATIC_MEMBRANE_ENERGY
//energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
#else
//energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
energy += weight * thickness_ * longNonquadraticMembraneEnergy(U);
#endif
energy += weight * thickness_ * curvatureEnergy(DR);
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
#ifdef CURVATURE_WITH_WRYNESS
energy += weight * curvatureWithWryness(R,DR);
#else
energy += weight * curvatureEnergy(DR);
#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
///////////////////////////////////////////////////////////
// Volume load contribution
///////////////////////////////////////////////////////////
if (not volumeLoad_)
continue;
// Value of the volume load density at the current position
auto volumeLoadDensity = volumeLoad_(element.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
energy += (volumeLoadDensity[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
//////////////////////////////////////////////////////////////////////////////
// Assemble boundary contributions
//////////////////////////////////////////////////////////////////////////////
if (not neumannFunction_)
return energy;
for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
{
if (not neumannBoundary_ or not neumannBoundary_->contains(it))
continue;
const Dune::QuadratureRule<DT, gridDim-1>& quad
= Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function
Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
energy += (neumannValue[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
}
return energy;
}
template <class Basis, int dim, class field_type>
typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
CosseratEnergyLocalStiffness<Basis,dim,field_type>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
{
auto element = localView.element();
RT energy = 0;
using namespace Dune::Indices;
const auto& deformationLocalFiniteElement = LocalFiniteElementFactory<Basis,0>::get(localView,_0);
const auto& orientationLocalFiniteElement = LocalFiniteElementFactory<Basis,1>::get(localView,_1);
#ifdef PROJECTED_INTERPOLATION
typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
LocalDeformationGFEFunctionType;
#else
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RealTuple<field_type,dim> >
LocalDeformationGFEFunctionType;
#endif
LocalDeformationGFEFunctionType localDeformationGFEFunction(deformationLocalFiniteElement,localConfiguration[_0]);
#ifdef PROJECTED_INTERPOLATION
typedef Dune::GFE::LocalProjectedFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
#else
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(orientationLocalFiniteElement), Rotation<field_type,dim> > LocalOrientationGFEFunctionType;
#endif
LocalOrientationGFEFunctionType localOrientationGFEFunction(orientationLocalFiniteElement,localConfiguration[_1]);
// \todo Implement smarter quadrature rule selection for more efficiency, i.e., less evaluations of the Rotation GFE function
int quadOrder = deformationLocalFiniteElement.localBasis().order() * ((element.type().isSimplex()) ? 1 : gridDim);
const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++)
{
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = quad[pt].position();
const DT integrationElement = element.geometry().integrationElement(quadPos);
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
DT weight = quad[pt].weight() * integrationElement;
// The value of the local deformation
RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
Rotation<field_type,dim> orientationValue = localOrientationGFEFunction.evaluate(quadPos);
// The derivative of the local function defined on the reference element
typename LocalDeformationGFEFunctionType::DerivativeType deformationReferenceDerivative = localDeformationGFEFunction.evaluateDerivative(quadPos,deformationValue);
typename LocalOrientationGFEFunctionType::DerivativeType orientationReferenceDerivative = localOrientationGFEFunction.evaluateDerivative(quadPos,orientationValue);
// The derivative of the function defined on the actual element
typename LocalDeformationGFEFunctionType::DerivativeType deformationDerivative;
typename LocalOrientationGFEFunctionType::DerivativeType orientationDerivative;
for (size_t comp=0; comp<deformationReferenceDerivative.N(); comp++)
jacobianInverseTransposed.mv(deformationReferenceDerivative[comp], deformationDerivative[comp]);
for (size_t comp=0; comp<orientationReferenceDerivative.N(); comp++)
jacobianInverseTransposed.mv(orientationReferenceDerivative[comp], orientationDerivative[comp]);
/////////////////////////////////////////////////////////
// compute U, the Cosserat strain
/////////////////////////////////////////////////////////
static_assert(dim>=gridDim, "Codim of the grid must be nonnegative");
//
Dune::FieldMatrix<field_type,dim,dim> R;
orientationValue.matrix(R);
Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(deformationDerivative,R);
//////////////////////////////////////////////////////////
// Compute the derivative of the rotation
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
Tensor3<field_type,3,3,gridDim> DR = orientationValue.quaternionTangentToMatrixTangent(orientationDerivative);
// Add the local energy density
if (gridDim==2) {
#ifdef QUADRATIC_MEMBRANE_ENERGY
//energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
#else
energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
#endif
energy += weight * thickness_ * curvatureEnergy(DR);
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
} else if (gridDim==3) {
energy += weight * quadraticMembraneEnergy(U);
#ifdef CURVATURE_WITH_WRYNESS
energy += weight * curvatureWithWryness(R,DR);
#else
energy += weight * curvatureEnergy(DR);
#endif
} else
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
///////////////////////////////////////////////////////////
// Volume load contribution
///////////////////////////////////////////////////////////
if (not volumeLoad_)
continue;
// Value of the volume load density at the current position
auto volumeLoadDensity = volumeLoad_(element.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
energy += (volumeLoadDensity[i] * deformationValue[i]) * quad[pt].weight() * integrationElement;
}
//////////////////////////////////////////////////////////////////////////////
// Assemble boundary contributions
//////////////////////////////////////////////////////////////////////////////
if (not neumannFunction_)
return energy;
for (auto&& it : intersections(neumannBoundary_->gridView(),element) )
{
if (not neumannBoundary_ or not neumannBoundary_->contains(it))
continue;
const auto& quad = Dune::QuadratureRules<DT, gridDim-1>::rule(it.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,gridDim>& quadPos = it.geometryInInside().global(quad[pt].position());
const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function
RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
energy += (neumannValue[i] * deformationValue.globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
}
return energy;
}
#endif //#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
dune/gfe/assemblers/cosseratrodenergy.hh
View file @ 48b9358e
......@@ -2,6 +2,7 @@
#define DUNE_GFE_COSSERATRODENERGY_HH
#include <array>
#include <memory>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
......@@ -18,9 +19,10 @@
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
namespace Dune::GFE {
namespace Dune::GFE
{
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
class CosseratRodEnergy
: public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
{
......@@ -42,6 +44,14 @@ namespace Dune::GFE {
public:
/** \brief The current configuration whose energy is being computed
*/
const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
/** \brief The GFE function that implements the stress-free reference configuration
*/
const std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration_;
/** \brief The stress-free configuration
The number type cannot be RT, because RT become `adouble` when
......@@ -67,9 +77,25 @@ namespace Dune::GFE {
GridView gridView_;
//! Constructor
CosseratRodEnergy(const GridView& gridView,
CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
: localGFEFunction_(localGFEFunction)
, localReferenceConfiguration_(localReferenceConfiguration)
, K_(K),
A_(A),
gridView_(gridView)
{}
//! Constructor
CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
ReferenceInterpolationRule&& localReferenceConfiguration,
const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
: K_(K),
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
, localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
, K_(K),
A_(A),
gridView_(gridView)
{}
......@@ -80,9 +106,39 @@ namespace Dune::GFE {
\param E Young's modulus
\param nu Poisson number
*/
CosseratRodEnergy(const GridView& gridView,
CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
const GridView& gridView,
double A, double J1, double J2, double E, double nu)
: gridView_(gridView)
: localGFEFunction_(localGFEFunction)
, localReferenceConfiguration_(localReferenceConfiguration)
, gridView_(gridView)
{
// shear modulus
double G = E/(2+2*nu);
K_[0] = E * J1;
K_[1] = E * J2;
K_[2] = G * (J1 + J2);
A_[0] = G * A;
A_[1] = G * A;
A_[2] = E * A;
}
/** \brief Constructor setting shape constants and material parameters
\param A The rod section area
\param J1, J2 The geometric moments (Flächenträgheitsmomente)
\param E Young's modulus
\param nu Poisson number
*/
CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
ReferenceInterpolationRule&& localReferenceConfiguration,
const GridView& gridView,
double A, double J1, double J2, double E, double nu)
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
, localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
, gridView_(gridView)
{
// shear modulus
double G = E/(2+2*nu);
......@@ -172,21 +228,18 @@ namespace Dune::GFE {
};
template<class Basis, class LocalInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localCoefficients) const
{
const auto& localFiniteElement = localView.tree().finiteElement();
LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
const auto& element = localView.element();
localGFEFunction_->bind(element, localCoefficients);
RT energy = 0;
std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::other;
InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
localReferenceConfiguration_->bind(element, localReferenceCoefficients);
// ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate
......@@ -205,10 +258,10 @@ namespace Dune::GFE {
double weight = shearingQuad[pt].weight() * integrationElement;
auto strain = getStrain(localConfiguration, element, quadPos);
auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
for (int i=0; i<3; i++)
energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
......@@ -225,10 +278,10 @@ namespace Dune::GFE {
double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
auto strain = getStrain(localConfiguration, element, quadPos);
auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
// Part II: the bending and twisting energy
for (int i=0; i<3; i++)
......@@ -240,9 +293,9 @@ namespace Dune::GFE {
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class ReboundLocalInterpolationRule>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const ReboundLocalInterpolationRule& localInterpolation,
const Entity& element,
const FieldVector<double,1>& pos) const
......@@ -285,9 +338,9 @@ namespace Dune::GFE {
return strain;
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class Number>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element,
const FieldVector<double, 1>& pos) const
......@@ -308,8 +361,8 @@ namespace Dune::GFE {
return stress;
}
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const
{
......@@ -363,8 +416,8 @@ namespace Dune::GFE {
}
}
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const
{
......@@ -386,9 +439,9 @@ namespace Dune::GFE {
}
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class PatchGridView>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const
{
......
dune/gfe/assemblers/discretekirchhoffbendingenergy.hh 0 → 100644
View file @ 48b9358e
#ifndef DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
#define DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/bendingisometryhelper.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
/**
* \brief Assemble the discrete Kirchhoff bending energy for a single element.
*
* The energy is essentially the harmonic energy of a
* discrete Jacobian operator which is represented as a
* linear combination in a P2 finite element space.
* The gradient of that linear combination serves as an
* discrete Hessian in the energy.
*/
namespace Dune::GFE
{
template<class Basis, class DiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
class DiscreteKirchhoffBendingEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
typedef typename DiscreteKirchhoffFunction::DerivativeType DerivativeType;
// Grid dimension
constexpr static int gridDim = GridView::dimension;
// P2-LocalFiniteElement used to represent the discrete Jacobian on the current element.
typedef typename Dune::LagrangeSimplexLocalFiniteElement<DT, double, gridDim, 2> P2LocalFiniteElement;
public:
DiscreteKirchhoffBendingEnergy(DiscreteKirchhoffFunction &localFunction)
: localFunction_(localFunction)
{}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(NotImplemented, "!");
}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localConfiguration) const override
{
const auto element = localView.element();
int quadOrder = 3;
const auto &quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
/**
bind the underlying Hermite basis to the current element.
Note: The DiscreteKirchhoffFunction object stores a vector of global
and local basis coefficients that are accessed for evaluation
methods (linear combinations of hermite basis function and local coefficients).
Therefore the local configuration is passed to the bind method as well since
these coefficient vectors need to be updated with the new local configuration
previous to evaluation.
*/
localFunction_.bind(element,localConfiguration);
RT bendingEnergy = 0;
for (auto&& quadPoint : quad)
{
auto discreteHessian= localFunction_.evaluateDiscreteHessian(quadPoint.position());
auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
bendingEnergy += weight * discreteHessian.frobenius_norm2();
}
return 0.5 * bendingEnergy;
}
private:
DiscreteKirchhoffFunction& localFunction_;
P2LocalFiniteElement lagrangeLFE_;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/forceenergy.hh 0 → 100644
View file @ 48b9358e
#ifndef DUNE_GFE_FORCEENERGY_HH
#define DUNE_GFE_FORCEENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/bendingisometryhelper.hh>
/**
* \brief Assemble the energy of a force term (force times deformation) for a single element.
*/
namespace Dune::GFE
{
template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
class ForceEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
// Grid dimension
constexpr static int gridDim = GridView::dimension;
public:
ForceEnergy(LocalDiscreteKirchhoffFunction &localFunction,
LocalForce &localForce)
: localFunction_(localFunction),
localForce_(localForce)
{}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(NotImplemented, "!");
}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localConfiguration) const override
{
RT energy = 0;
int quadOrder = 3;
const auto &quad = QuadratureRules<double, gridDim>::rule(GeometryTypes::simplex(gridDim), quadOrder);
const auto element = localView.element();
/** bind the local force to the current element */
localForce_.bind(element);
/**
bind the underlying Hermite basis to the current element.
Note: The DiscreteKirchhoffFunction object stores a vector of global
and local basis coefficients that are accessed for evaluation
methods (linear comibnations of hermite basis function and local coefficients).
Therefore the local configuration is passed to the bind method as well since
these coefficient vectors need to be updated with the new local configuration
previous to evaluation.
*/
localFunction_.bind(element,localConfiguration);
/**
* @brief Compute contribution of the force Term
*
* Integrate the scalar product of the force
* with the local deformation function 'localFunction_'
* over the current element.
*/
RT forceTermEnergy = 0;
for (auto&& quadPoint : quad)
{
auto deformationValue = localFunction_.evaluate(quadPoint.position());
auto forceValue = localForce_(quadPoint.position());
auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
}
return forceTermEnergy;
}
private:
LocalDiscreteKirchhoffFunction& localFunction_;
LocalForce& localForce_;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/geodesicfeassembler.hh
View file @ 48b9358e
......@@ -10,269 +10,273 @@
#include <dune/solvers/common/wrapownshare.hh>
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
using field_type = typename TargetSpace::field_type;
typedef typename Basis::GridView GridView;
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
namespace Dune::GFE
{
//! Dimension of the grid.
constexpr static int gridDim = GridView::dimension;
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
using field_type = typename TargetSpace::field_type;
typedef typename Basis::GridView GridView;
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
//! Dimension of the grid.
constexpr static int gridDim = GridView::dimension;
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
protected:
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
//! The global basis
const Basis basis_;
//! The local stiffness operator
std::shared_ptr<LocalStiffness> localStiffness_;
protected:
public:
//! The global basis
const Basis basis_;
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const Basis& basis)
: basis_(basis)
{}
//! The local stiffness operator
std::shared_ptr<LocalStiffness> localStiffness_;
/** \brief Constructor for a given grid */
template <class LocalStiffnessT>
GeodesicFEAssembler(const Basis& basis,
LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
public:
/** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
template <class LocalStiffnessT>
void setLocalStiffness(LocalStiffnessT&& localStiffness) {
localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
}
/** \brief Get the local stiffness operator. */
const LocalStiffness& getLocalStiffness() const {
return *localStiffness_;
}
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const Basis& basis)
: basis_(basis)
{}
/** \brief Get the local stiffness operator. */
LocalStiffness& getLocalStiffness() {
return *localStiffness_;
}
/** \brief Constructor for a given grid */
template <class LocalStiffnessT>
GeodesicFEAssembler(const Basis& basis,
LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
/** \brief Get the basis. */
const Basis& getBasis() const {
return basis_;
}
/** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
template <class LocalStiffnessT>
void setLocalStiffness(LocalStiffnessT&& localStiffness) {
localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
}
/** \brief Assemble the tangent stiffness matrix and the functional gradient together
*
* This is more efficient than computing them separately, because you need the gradient
* anyway to compute the Riemannian Hessian.
*/
virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const;
/** \brief Get the local stiffness operator. */
const LocalStiffness& getLocalStiffness() const {
return *localStiffness_;
}
/** \brief Assemble the gradient */
virtual void assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Get the local stiffness operator. */
LocalStiffness& getLocalStiffness() {
return *localStiffness_;
}
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
/** \brief Get the basis. */
const Basis& getBasis() const {
return basis_;
}
//protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
/** \brief Assemble the tangent stiffness matrix and the functional gradient together
*
* This is more efficient than computing them separately, because you need the gradient
* anyway to compute the Riemannian Hessian.
*/
virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const;
}; // end class
/** \brief Assemble the gradient */
virtual void assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
//protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
auto n = basis_.size();
}; // end class
nb.resize(n, n);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
// Bind the local FE basis view to the current element
localView.bind(element);
auto n = basis_.size();
const auto& lfe = localView.tree().finiteElement();
nb.resize(n, n);
for (size_t i=0; i<lfe.size(); i++) {
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t j=0; j<lfe.size(); j++) {
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
const auto& lfe = localView.tree().finiteElement();
nb.add(iIdx, jIdx);
for (size_t i=0; i<lfe.size(); i++) {
}
for (size_t j=0; j<lfe.size(); j++) {
}
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
}
nb.add(iIdx, jIdx);
}
}
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
{
if (computeOccupationPattern) {
}
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
neighborsPerVertex.exportIdx(hessian);
}
}
hessian = 0;
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
{
if (computeOccupationPattern) {
gradient.resize(sol.size());
gradient = 0;
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
neighborsPerVertex.exportIdx(hessian);
// A view on the FE basis on a single element
auto localView = basis_.localView();
}
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
hessian = 0;
const int numOfBaseFct = localView.tree().size();
gradient.resize(sol.size());
gradient = 0;
// Extract local solution
std::vector<TargetSpace> localSolution(numOfBaseFct);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[localView.index(i)];
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
std::vector<double> localGradient(numOfBaseFct*blocksize);
Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
const int numOfBaseFct = localView.tree().size();
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
// Extract local solution
std::vector<TargetSpace> localSolution(numOfBaseFct);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[localView.index(i)];
auto row = localView.index(i);
std::vector<double> localGradient(numOfBaseFct*blocksize);
Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
for (int j=0; j<numOfBaseFct; j++ ) {
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
auto col = localView.index(j);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
for (std::size_t ii=0; ii<blocksize; ii++)
for (std::size_t jj=0; jj<blocksize; jj++)
hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
auto row = localView.index(i);
}
}
for (int j=0; j<numOfBaseFct; j++ ) {
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
for (int j=0; j<blocksize; j++)
gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
auto col = localView.index(j);
}
for (std::size_t ii=0; ii<blocksize; ii++)
for (std::size_t jj=0; jj<blocksize; jj++)
hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
}
}
}
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
for (int j=0; j<blocksize; j++)
gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
grad.resize(sol.size());
grad = 0;
}
// A view on the FE basis on a single element
auto localView = basis_.localView();
}
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
localView.bind(element);
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
// The number of degrees of freedom of the current element
const auto nDofs = localView.tree().size();
grad.resize(sol.size());
grad = 0;
// Extract local solution
std::vector<TargetSpace> localSolution(nDofs);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)];
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
// Assemble local gradient
std::vector<double> localGradient(nDofs*blocksize);
// The number of degrees of freedom of the current element
const auto nDofs = localView.tree().size();
localStiffness_->assembleGradient(localView, localSolution, localGradient);
// Extract local solution
std::vector<TargetSpace> localSolution(nDofs);
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<blocksize; j++)
grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)];
}
// Assemble local gradient
std::vector<double> localGradient(nDofs*blocksize);
localStiffness_->assembleGradient(localView, localSolution, localGradient);
template <class Basis, class TargetSpace>
double GeodesicFEAssembler<Basis, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
double energy = 0;
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<blocksize; j++)
grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
if (sol.size() != basis_.size())
DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
}
// A view on the FE basis on a single element
auto localView = basis_.localView();
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
double GeodesicFEAssembler<Basis, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
localView.bind(element);
double energy = 0;
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
if (sol.size() != basis_.size())
DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
std::vector<TargetSpace> localSolution(nDofs);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)[0]];
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
energy += localStiffness_->energy(localView, localSolution);
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
}
std::vector<TargetSpace> localSolution(nDofs);
return energy;
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)[0]];
}
energy += localStiffness_->energy(localView, localSolution);
}
return energy;
}
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
View file @ 48b9358e
......@@ -5,7 +5,8 @@
#include <dune/common/tuplevector.hh>
#include <dune/gfe/spaces/productmanifold.hh>
namespace Dune::GFE {
namespace Dune::GFE
{
/** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
......@@ -13,9 +14,8 @@ namespace Dune::GFE {
*/
template <class Basis, class ScalarBasis, class TargetSpace>
class
GeodesicFEAssemblerWrapper {
class GeodesicFEAssemblerWrapper
{
using MixedSpace0 = std::tuple_element_t<0,TargetSpace>;
using MixedSpace1 = std::tuple_element_t<1,TargetSpace>;
......@@ -72,7 +72,8 @@ namespace Dune::GFE {
auto splitVector(const std::vector<TargetSpace>& sol) const;
std::unique_ptr<MatrixType> hessianMixed_;
}; // end class
} //end namespace
} // namespace Dune::GFE
template <class Basis, class ScalarBasis, class TargetSpace>
......
dune/gfe/assemblers/l2distancesquaredenergy.hh
View file @ 48b9358e
......@@ -5,86 +5,94 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/functions/globalgfefunction.hh>
#include <dune/gfe/functions/localgeodesicfefunction.hh>
template<class Basis, class TargetSpace>
class L2DistanceSquaredEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
namespace Dune::GFE
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
// some other sizes
constexpr static int gridDim = GridView::dimension;
template<class Basis, class TargetSpace>
class L2DistanceSquaredEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, typename TargetSpace::template rebind<double>::other>;
// some other sizes
constexpr static int gridDim = GridView::dimension;
public:
using LocalInterpolationRule = LocalGeodesicFEFunction<Basis, typename TargetSpace::template rebind<double>::other>;
// This is the function that we are computing the L2-distance to
std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
public:
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
// This is the function that we are computing the L2-distance to
std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
{
RT energy = 0;
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
const auto& localFiniteElement = localView.tree().finiteElement();
typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
{
RT energy = 0;
const auto element = localView.element();
auto localOrigin = localFunction(*origin_);
localOrigin.bind(element);
typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localView.globalBasis());
// Just guessing an appropriate quadrature order
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
const auto element = localView.element();
localGeodesicFEFunction.bind(element,localSolution);
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
auto localOrigin = localFunction(*origin_);
localOrigin.bind(element);
for (size_t pt=0; pt<quad.size(); pt++)
{
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
// Just guessing an appropriate quadrature order
const auto& localFiniteElement = localView.tree().finiteElement();
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
const double integrationElement = element.geometry().integrationElement(quadPos);
for (size_t pt=0; pt<quad.size(); pt++)
{
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
auto weight = quad[pt].weight() * integrationElement;
const double integrationElement = element.geometry().integrationElement(quadPos);
// The function value
auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originValue = localOrigin(quadPos);
auto weight = quad[pt].weight() * integrationElement;
// The derivative of the 'origin' function
// First: as function defined on the reference element
auto originReferenceDerivative = derivative(localOrigin)(quadPos);
// The function value
auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originValue = localOrigin(quadPos);
// The derivative of the function defined on the actual element
auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
// The derivative of the 'origin' function
// First: as function defined on the reference element
auto originReferenceDerivative = derivative(localOrigin)(quadPos);
double weightFactor = originDerivative.frobenius_norm();
// un-comment the following line to switch off the weight factor
//weightFactor = 1.0;
// The derivative of the function defined on the actual element
auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
// Add the local energy density
energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
double weightFactor = originDerivative.frobenius_norm();
// un-comment the following line to switch off the weight factor
//weightFactor = 1.0;
// Add the local energy density
energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
}
return energy;
}
return energy;
}
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localenergy.hh
View file @ 48b9358e
......@@ -7,9 +7,9 @@
#include <dune/gfe/spaces/productmanifold.hh>
namespace Dune {
namespace GFE:: Impl
namespace Dune::GFE
{
namespace Impl
{
/** \brief A class exporting container types for coefficient sets
*
......@@ -34,43 +34,40 @@ namespace Dune {
using Coefficients = std::vector<ProductManifold<Factors...> >;
using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
};
}
namespace GFE {
} // namespace Impl
/** \brief Base class for energies defined by integrating over one grid element */
template<class Basis, class TargetSpace>
class LocalEnergy
{
public:
using RT = typename TargetSpace::ctype;
/** \brief Compute the energy
*
* \param localView Local view specifying the current element and the FE space there
* \param coefficients The coefficients of a FE function on the current element
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** \brief Base class for energies defined by integrating over one grid element */
template<class Basis, class TargetSpace>
class LocalEnergy
{
public:
using RT = typename TargetSpace::ctype;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
/** \brief Compute the energy
*
* \param localView Local view specifying the current element and the FE space there
* \param coefficients The coefficients of a FE function on the current element
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** Empty virtual default destructor
*
* To allow proper destruction of derived classes through a base class pointer
*/
virtual ~LocalEnergy() = default;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
};
/** Empty virtual default destructor
*
* To allow proper destruction of derived classes through a base class pointer
*/
virtual ~LocalEnergy() = default;
} // namespace GFE
};
} // namespace Dune
} // namespace Dune::GFE
#endif // DUNE_GFE_LOCALENERGY_HH
dune/gfe/assemblers/localfirstordermodel.hh
View file @ 48b9358e
......@@ -20,6 +20,6 @@ namespace Dune::GFE
};
} // namespace Dune
} // namespace Dune::GFE
#endif // DUNE_GFE_LOCALFIRSTORDERMODEL_HH
dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
View file @ 48b9358e
......@@ -14,145 +14,83 @@
#include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
#include <dune/gfe/spaces/productmanifold.hh>
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class Basis, class TargetSpace>
class LocalGeodesicFEADOLCStiffness
: public LocalGeodesicFEStiffness<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
constexpr static int gridDim = GridView::dimension;
// The 'active' target spaces, i.e., the number type is replaced by adouble
using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
public:
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
LocalGeodesicFEADOLCStiffness(std::shared_ptr<const Dune::GFE::LocalEnergy<Basis, ATargetSpace> > energy, bool adolcScalarMode = false)
: localEnergy_(energy),
adolcScalarMode_(adolcScalarMode)
{}
/** \brief Compute the energy at the current configuration */
virtual RT energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients) const override;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
/** \brief Assemble the element gradient of the energy functional
*/
virtual void assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& gradient) const override;
/** \brief Assemble the local stiffness matrix at the current position
This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
/** \brief Assemble the local stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override;
std::shared_ptr<const Dune::GFE::LocalEnergy<Basis, ATargetSpace> > localEnergy_;
const bool adolcScalarMode_;
};
template <class Basis, class TargetSpace>
typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT
LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution) const
namespace Dune::GFE
{
double pureEnergy;
int rank = localView.globalBasis().gridView().comm().rank();
std::vector<ATargetSpace> localASolution(localSolution.size());
trace_on(rank);
adouble energy = 0;
// The following loop is not quite intuitive: we copy the localSolution into an
// array of FieldVector<double>, go from there to FieldVector<adouble> and
// only then to ATargetSpace.
// Rationale: The constructor/assignment-from-vector of TargetSpace frequently
// contains a projection onto the manifold from the surrounding Euclidean space.
// ADOL-C needs a function on the whole Euclidean space, hence that projection
// is part of the function and needs to be taped.
// The following variable cannot be declared inside of the loop, or ADOL-C will report wrong results
// (Presumably because several independent variables use the same memory location.)
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] <<= raw[j];
localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble
}
try {
energy = localEnergy_->energy(localView,localASolution);
} catch (Dune::Exception &e) {
trace_off();
throw e;
}
energy >>= pureEnergy;
trace_off();
return pureEnergy;
}
template <class Basis, class TargetSpace>
typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT
LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
{
using namespace Dune::Indices;
int rank = Dune::MPIHelper::getCommunication().rank();
if constexpr (not Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class Basis, class TargetSpace>
class LocalGeodesicFEADOLCStiffness
: public LocalGeodesicFEStiffness<Basis,TargetSpace>
{
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==1);
return energy(localView, localConfiguration[_0]);
}
else
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
constexpr static int gridDim = GridView::dimension;
// The 'active' target spaces, i.e., the number type is replaced by adouble
using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
public:
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
LocalGeodesicFEADOLCStiffness(std::shared_ptr<const Dune::GFE::LocalEnergy<Basis, ATargetSpace> > energy, bool adolcScalarMode = false)
: localEnergy_(energy),
adolcScalarMode_(adolcScalarMode)
{}
/** \brief Compute the energy at the current configuration */
virtual RT energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients) const override;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override;
/** \brief Assemble the element gradient of the energy functional
*/
virtual void assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& gradient) const override;
/** \brief Assemble the local stiffness matrix at the current position
This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
/** \brief Assemble the local stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override;
std::shared_ptr<const Dune::GFE::LocalEnergy<Basis, ATargetSpace> > localEnergy_;
const bool adolcScalarMode_;
};
template <class Basis, class TargetSpace>
typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT
LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution) const
{
using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
using ATargetSpace0 = typename TargetSpace0::template rebind<adouble>::other;
using ATargetSpace1 = typename TargetSpace1::template rebind<adouble>::other;
Dune::TupleVector<std::vector<ATargetSpace0>, std::vector<ATargetSpace1> > localAConfiguration;
localAConfiguration[_0].resize(localConfiguration[_0].size());
localAConfiguration[_1].resize(localConfiguration[_1].size());
double pureEnergy;
int rank = localView.globalBasis().gridView().comm().rank();
std::vector<ATargetSpace> localASolution(localSolution.size());
trace_on(rank);
adouble energy = 0;
......@@ -167,374 +105,196 @@ energy(const typename Basis::LocalView& localView,
// The following variable cannot be declared inside of the loop, or ADOL-C will report wrong results
// (Presumably because several independent variables use the same memory location.)
std::vector<typename ATargetSpace0::CoordinateType> aRaw0(localConfiguration[_0].size());
for (size_t i=0; i<localConfiguration[_0].size(); i++) {
typename TargetSpace0::CoordinateType raw = localConfiguration[_0][i].globalCoordinates();
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw0[i][j] <<= raw[j];
localAConfiguration[_0][i] = aRaw0[i]; // may contain a projection onto M -- needs to be done in adouble
}
std::vector<typename ATargetSpace1::CoordinateType> aRaw1(localConfiguration[_1].size());
for (size_t i=0; i<localConfiguration[_1].size(); i++) {
typename TargetSpace1::CoordinateType raw = localConfiguration[_1][i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw1[i][j] <<= raw[j];
localAConfiguration[_1][i] = aRaw1[i]; // may contain a projection onto M -- needs to be done in adouble
aRaw[i][j] <<= raw[j];
localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble
}
try {
energy = localEnergy_->energy(localView, localAConfiguration);
energy = localEnergy_->energy(localView,localASolution);
} catch (Dune::Exception &e) {
trace_off();
throw e;
}
energy >>= pureEnergy;
trace_off();
return pureEnergy;
}
}
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
std::vector<double>& localGradient) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
int rank = localView.globalBasis().gridView().comm().rank();
// Compute the actual gradient
constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
size_t nDofs = localSolution.size();
size_t nDoubles = nDofs*embeddedBlocksize;
std::vector<double> xp(nDoubles);
int idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
xp[idx++] = localSolution[i].globalCoordinates()[j];
// Compute gradient
std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
// Copy into Dune type
std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
orthonormalFrame.mv(localEmbeddedGradient[i],localGradient[i]);
}
}
// ///////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
int rank = localView.globalBasis().gridView().comm().rank();
/////////////////////////////////////////////////////////////////
// Compute the gradient. It is needed to transform the Hessian
// into the correct coordinates.
/////////////////////////////////////////////////////////////////
// Compute the actual gradient
constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
size_t nDofs = localSolution.size();
size_t nDoubles = nDofs*embeddedBlocksize;
std::vector<double> xp(nDoubles);
int idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
xp[idx++] = localSolution[i].globalCoordinates()[j];
// Compute gradient
std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
// Copy into Dune type
std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
typename TargetSpace::TangentVector tmp;
localSolution[i].orthonormalFrame().mv(localEmbeddedGradient[i],tmp);
for (size_t j=0; j<blocksize; j++)
localGradient[i*blocksize+j] = tmp[j];
}
/////////////////////////////////////////////////////////////////
// Compute Hessian
/////////////////////////////////////////////////////////////////
// We compute the Hessian of the energy functional using the ADOL-C system.
// Since ADOL-C does not know about nonlinear spaces, what we get is actually
// the Hessian of a prolongation of the energy functional into the surrounding
// Euclidean space. To obtain the Riemannian Hessian from this we apply the
// formula described in Absil, Mahoney, Trumpf, "An extrinsic look at the Riemannian Hessian".
// This formula consists of two steps:
// 1) Remove all entries of the Hessian pertaining to the normal space of the
// manifold. In the aforementioned paper this is done by projection onto the
// tangent space. Since we want a matrix that is really smaller (but full rank again),
// we can achieve the same effect by multiplying the embedded Hessian from the left
// and from the right by the matrix of orthonormal frames.
// 2) Add a correction involving the Weingarten map.
//
// This works, and is easy to implement using the ADOL-C "hessian" driver.
// However, here we implement a small shortcut. Computing the embedded Hessian and
// multiplying one side by the orthonormal frame is the same as evaluating the Hessian
// (seen as an operator from R^n to R^n) in the directions of the vectors of the
// orthonormal frame. By luck, ADOL-C can compute the evaluations of the Hessian in
// a given direction directly (in fact, this is also how the "hessian" driver works).
// Since there are less frame vectors than the dimension of the embedding space,
// this reinterpretation allows to reduce the number of calls to ADOL-C.
// In my Cosserat shell tests this reduced assembly time by about 10%.
std::vector<Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> > orthonormalFrame(nDofs);
for (size_t i=0; i<nDofs; i++)
orthonormalFrame[i] = localSolution[i].orthonormalFrame();
Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
if (adolcScalarMode_) {
std::vector<double> v(nDoubles);
std::vector<double> w(nDoubles);
std::fill(v.begin(), v.end(), 0.0);
for (size_t i=0; i<nDofs; i++)
for (int ii=0; ii<blocksize; ii++)
{
// Evaluate Hessian in the direction of each vector of the orthonormal frame
for (size_t k=0; k<embeddedBlocksize; k++)
v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
template <class Basis, class TargetSpace>
typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT
LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
{
using namespace Dune::Indices;
for (size_t j=0; j<nDoubles; j++)
embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
int rank = Dune::MPIHelper::getCommunication().rank();
// Make v the null vector again
std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
}
} else { //ADOL-C vector mode
int n = nDoubles;
int nDirections = nDofs * blocksize;
double* tangent[nDoubles];
for(size_t i=0; i<nDoubles; i++)
tangent[i] = (double*)malloc(nDirections*sizeof(double));
double* rawHessian[nDoubles];
for(size_t i=0; i<nDoubles; i++)
rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
for (int j=0; j<nDirections; j++)
if constexpr (not Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
{
for (int i=0; i<n; i++)
tangent[i][j] = 0.0;
for (int i=0; i<embeddedBlocksize; i++)
tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==1);
return energy(localView, localConfiguration[_0]);
}
hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type
for(size_t i=0; i<nDoubles; i++)
for (int j=0; j<nDirections; j++)
embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
for(size_t i=0; i<nDoubles; i++) {
free(rawHessian[i]);
free(tangent[i]);
}
}
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
// at the Riemannian Hessian".
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
typedef typename TargetSpace::TangentVector TangentVector;
localHessian.setSize(nDofs*blocksize,nDofs*blocksize);
for (size_t col=0; col<nDofs; col++) {
else
{
using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
using ATargetSpace0 = typename TargetSpace0::template rebind<adouble>::other;
using ATargetSpace1 = typename TargetSpace1::template rebind<adouble>::other;
Dune::TupleVector<std::vector<ATargetSpace0>, std::vector<ATargetSpace1> > localAConfiguration;
localAConfiguration[_0].resize(localConfiguration[_0].size());
localAConfiguration[_1].resize(localConfiguration[_1].size());
double pureEnergy;
trace_on(rank);
adouble energy = 0;
// The following loop is not quite intuitive: we copy the localSolution into an
// array of FieldVector<double>, go from there to FieldVector<adouble> and
// only then to ATargetSpace.
// Rationale: The constructor/assignment-from-vector of TargetSpace frequently
// contains a projection onto the manifold from the surrounding Euclidean space.
// ADOL-C needs a function on the whole Euclidean space, hence that projection
// is part of the function and needs to be taped.
// The following variable cannot be declared inside of the loop, or ADOL-C will report wrong results
// (Presumably because several independent variables use the same memory location.)
std::vector<typename ATargetSpace0::CoordinateType> aRaw0(localConfiguration[_0].size());
for (size_t i=0; i<localConfiguration[_0].size(); i++) {
typename TargetSpace0::CoordinateType raw = localConfiguration[_0][i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw0[i][j] <<= raw[j];
localAConfiguration[_0][i] = aRaw0[i]; // may contain a projection onto M -- needs to be done in adouble
}
for (size_t subCol=0; subCol<blocksize; subCol++) {
std::vector<typename ATargetSpace1::CoordinateType> aRaw1(localConfiguration[_1].size());
for (size_t i=0; i<localConfiguration[_1].size(); i++) {
typename TargetSpace1::CoordinateType raw = localConfiguration[_1][i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw1[i][j] <<= raw[j];
localAConfiguration[_1][i] = aRaw1[i]; // may contain a projection onto M -- needs to be done in adouble
}
EmbeddedTangentVector z = orthonormalFrame[col][subCol];
try {
// P_x \partial^2 f z
for (size_t row=0; row<nDofs; row++) {
TangentVector semiEmbeddedProduct;
embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
energy = localEnergy_->energy(localView, localAConfiguration);
for (int subRow=0; subRow<blocksize; subRow++)
localHessian[row*blocksize+subRow][col*blocksize+subCol] = semiEmbeddedProduct[subRow];
} catch (Dune::Exception &e) {
trace_off();
throw e;
}
energy >>= pureEnergy;
}
trace_off();
return pureEnergy;
}
}
//////////////////////////////////////////////////////////////////////////////////////
// Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f)
//////////////////////////////////////////////////////////////////////////////////////
// Project embedded gradient onto normal space
std::vector<typename TargetSpace::EmbeddedTangentVector> projectedGradient(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++)
projectedGradient[i] = localSolution[i].projectOntoNormalSpace(localEmbeddedGradient[i]);
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
std::vector<double>& localGradient) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
for (size_t row=0; row<nDofs; row++) {
int rank = localView.globalBasis().gridView().comm().rank();
for (size_t subRow=0; subRow<blocksize; subRow++) {
// Compute the actual gradient
constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
size_t nDofs = localSolution.size();
size_t nDoubles = nDofs*embeddedBlocksize;
std::vector<double> xp(nDoubles);
int idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
xp[idx++] = localSolution[i].globalCoordinates()[j];
EmbeddedTangentVector z = orthonormalFrame[row][subRow];
EmbeddedTangentVector tmp1 = localSolution[row].weingarten(z,projectedGradient[row]);
// Compute gradient
std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
TangentVector tmp2;
orthonormalFrame[row].mv(tmp1,tmp2);
// Copy into Dune type
std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
for (size_t subCol=0; subCol<blocksize; subCol++)
localHessian[row*blocksize+subRow][row*blocksize+subCol] += tmp2[subCol];
}
idx=0;
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
orthonormalFrame.mv(localEmbeddedGradient[i],localGradient[i]);
}
}
}
/////////////////////////////////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
/////////////////////////////////////////////////////////////////////////////////////
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const
{
using namespace Dune::Indices;
if constexpr (not Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
{
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==1);
return assembleGradientAndHessian(localView,
localConfiguration[_0],
localGradient,
localHessian[_0][_0]);
}
else
{
int rank = Dune::MPIHelper::getCommunication().rank();
// ///////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
energy(localView, localSolution);
int rank = localView.globalBasis().gridView().comm().rank();
/////////////////////////////////////////////////////////////////
// Compute the gradient. It is needed to transform the Hessian
// into the correct coordinates.
/////////////////////////////////////////////////////////////////
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==2, "ADOL-C assembly only implemented for product spaces with two factors");
//! Dimension of the embedding space
using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
constexpr static int embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension;
constexpr static int embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension;
// Compute the actual gradient
size_t nDofs0 = localConfiguration[_0].size();
size_t nDofs1 = localConfiguration[_1].size();
size_t nDoubles = nDofs0*embeddedBlocksize0 + nDofs1*embeddedBlocksize1;
constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
size_t nDofs = localSolution.size();
size_t nDoubles = nDofs*embeddedBlocksize;
std::vector<double> xp(nDoubles);
int idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++)
for (size_t j=0; j<embeddedBlocksize0; j++)
xp[idx++] = localConfiguration[_0][i].globalCoordinates()[j];
for (size_t i=0; i<localConfiguration[_1].size(); i++)
for (size_t j=0; j<embeddedBlocksize1; j++)
xp[idx++] = localConfiguration[_1][i].globalCoordinates()[j];
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
xp[idx++] = localSolution[i].globalCoordinates()[j];
// Compute gradient
std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data());
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
// Copy into Dune type
std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration[_0].size());
std::vector<typename TargetSpace1::EmbeddedTangentVector> localEmbeddedGradient1(localConfiguration[_1].size());
std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++) {
for (size_t j=0; j<embeddedBlocksize0; j++)
localEmbeddedGradient0[i][j] = g[idx++];
// Express gradient in local coordinate system
typename TargetSpace0::TangentVector tmp;
localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],tmp);
for (size_t j=0; j<blocksize0; j++)
localGradient[i*blocksize0+j] = tmp[j];
}
auto offset = localConfiguration[_0].size()*blocksize0;
for (size_t i=0; i<localConfiguration[_1].size(); i++) {
for (size_t j=0; j<embeddedBlocksize1; j++)
localEmbeddedGradient1[i][j] = g[idx++];
// Express gradient in local coordinate system
typename TargetSpace1::TangentVector tmp;
localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],tmp);
for (size_t j=0; j<blocksize1; j++)
localGradient[offset + i*blocksize1 + j] = tmp[j];
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient[i][j] = g[idx++];
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
typename TargetSpace::TangentVector tmp;
localSolution[i].orthonormalFrame().mv(localEmbeddedGradient[i],tmp);
for (size_t j=0; j<blocksize; j++)
localGradient[i*blocksize+j] = tmp[j];
}
/////////////////////////////////////////////////////////////////
......@@ -564,81 +324,40 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// this reinterpretation allows to reduce the number of calls to ADOL-C.
// In my Cosserat shell tests this reduced assembly time by about 10%.
std::vector<Dune::FieldMatrix<RT,blocksize0,embeddedBlocksize0> > orthonormalFrame0(nDofs0);
std::vector<Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> > orthonormalFrame(nDofs);
for (size_t i=0; i<nDofs0; i++)
orthonormalFrame0[i] = localConfiguration[_0][i].orthonormalFrame();
std::vector<Dune::FieldMatrix<RT,blocksize1,embeddedBlocksize1> > orthonormalFrame1(nDofs1);
for (size_t i=0; i<nDofs1; i++)
orthonormalFrame1[i] = localConfiguration[_1][i].orthonormalFrame();
for (size_t i=0; i<nDofs; i++)
orthonormalFrame[i] = localSolution[i].orthonormalFrame();
Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize0> > embeddedHessian00(nDofs0,nDofs0);
Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize1> > embeddedHessian01(nDofs0,nDofs1);
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize0> > embeddedHessian10(nDofs1,nDofs0);
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize1> > embeddedHessian11(nDofs1,nDofs1);
Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
if(adolcScalarMode_) {
if (adolcScalarMode_) {
std::vector<double> v(nDoubles);
std::vector<double> w(nDoubles);
std::fill(v.begin(), v.end(), 0.0);
size_t nDoubles0 = nDofs0*embeddedBlocksize0; // nDoubles = nDoubles0 + nDoubles1
size_t nDoubles1 = nDofs1*embeddedBlocksize1;
std::fill(v.begin(), v.end(), 0.0);
for (size_t i=0; i<nDofs0 + nDofs1; i++) {
// Evaluate Hessian in the direction of each vector of the orthonormal frame
if (i < nDofs0) { //Upper half
auto i0 = i;
for (int ii0=0; ii0<blocksize0; ii0++) {
for (size_t k0=0; k0<embeddedBlocksize0; k0++) {
v[i0*embeddedBlocksize0 + k0] = orthonormalFrame0[i0][ii0][k0];
}
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
for (size_t j0=0; j0<nDoubles0; j0++) //Upper left
embeddedHessian00[i0][j0/embeddedBlocksize0][ii0][j0%embeddedBlocksize0] = w[j0];
for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
embeddedHessian01[i0][j1/embeddedBlocksize1][ii0][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
// Make v the null vector again
std::fill(&v[i0*embeddedBlocksize0], &v[(i0+1)*embeddedBlocksize0], 0.0);
}
} else { // i = nDofs0 ... nDofs0 + nDofs1 //lower half
auto i1 = i - nDofs0;
for (int ii1=0; ii1<blocksize1; ii1++) {
for (size_t k1=0; k1<embeddedBlocksize1; k1++) {
v[nDoubles0 + i1*embeddedBlocksize1 + k1] = orthonormalFrame1[i1][ii1][k1];
}
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
for (size_t i=0; i<nDofs; i++)
for (int ii=0; ii<blocksize; ii++)
{
// Evaluate Hessian in the direction of each vector of the orthonormal frame
for (size_t k=0; k<embeddedBlocksize; k++)
v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
for (size_t j0=0; j0<nDoubles0; j0++) // Upper left
embeddedHessian10[i1][j0/embeddedBlocksize0][ii1][j0%embeddedBlocksize0] = w[j0];
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
embeddedHessian11[i1][j1/embeddedBlocksize1][ii1][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
for (size_t j=0; j<nDoubles; j++)
embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
// Make v the null vector again
std::fill(&v[nDoubles0 + i1*embeddedBlocksize1], &v[nDoubles0 + (i1+1)*embeddedBlocksize1], 0.0);
}
// Make v the null vector again
std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
}
}
} else { //ADOL-C vector mode}
} else { //ADOL-C vector mode
int n = nDoubles;
int nDirections = nDofs0 * blocksize0 + nDofs1 * blocksize1;
int nDirections = nDofs * blocksize;
double* tangent[nDoubles];
for(size_t i=0; i<nDoubles; i++)
tangent[i] = (double*)malloc(nDirections*sizeof(double));
......@@ -647,182 +366,469 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
for(size_t i=0; i<nDoubles; i++)
rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
// Initialize directions field with zeros
for (int j=0; j<nDirections; j++)
{
for (int i=0; i<n; i++)
tangent[i][j] = 0.0;
for (size_t j=0; j<nDofs0*blocksize0; j++)
for (int i=0; i<embeddedBlocksize0; i++)
tangent[(j/blocksize0)*embeddedBlocksize0+i][j] = orthonormalFrame0[j/blocksize0][j%blocksize0][i];
for (size_t j=0; j<nDofs1*blocksize1; j++)
for (int i=0; i<embeddedBlocksize1; i++)
tangent[nDofs0*embeddedBlocksize0 + (j/blocksize1)*embeddedBlocksize1+i][nDofs0*blocksize0 + j] = orthonormalFrame1[j/blocksize1][j%blocksize1][i];
for (int i=0; i<embeddedBlocksize; i++)
tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
}
hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type
size_t offset0 = nDofs0*embeddedBlocksize0;
size_t offset1 = nDofs0*blocksize0;
// upper left block
for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
for (size_t j=0; j<nDofs0*blocksize0; j++)
embeddedHessian00[j/blocksize0][i/embeddedBlocksize0][j%blocksize0][i%embeddedBlocksize0] = rawHessian[i][j];
// upper right block
for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
for (size_t j=0; j<nDofs0*blocksize0; j++)
embeddedHessian01[j/blocksize0][i/embeddedBlocksize1][j%blocksize0][i%embeddedBlocksize1] = rawHessian[offset0+i][j];
// lower left block
for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
for (size_t j=0; j<nDofs1*blocksize1; j++)
embeddedHessian10[j/blocksize1][i/embeddedBlocksize0][j%blocksize1][i%embeddedBlocksize0] = rawHessian[i][offset1+j];
// lower right block
for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
for (size_t j=0; j<nDofs1*blocksize1; j++)
embeddedHessian11[j/blocksize1][i/embeddedBlocksize1][j%blocksize1][i%embeddedBlocksize1] = rawHessian[offset0+i][offset1+j];
for(size_t i=0; i<nDoubles; i++)
for (int j=0; j<nDirections; j++)
embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
for(size_t i=0; i<nDoubles; i++) {
free(rawHessian[i]);
free(tangent[i]);
}
}
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
// at the Riemannian Hessian".
using namespace Dune::Indices;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
typedef typename TargetSpace::TangentVector TangentVector;
localHessian[_0][_0].setSize(nDofs0*blocksize0, nDofs0*blocksize0);
localHessian.setSize(nDofs*blocksize,nDofs*blocksize);
for (size_t col=0; col<nDofs0; col++)
{
for (size_t subCol=0; subCol<blocksize0; subCol++)
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol];
for (size_t col=0; col<nDofs; col++) {
for (size_t subCol=0; subCol<blocksize; subCol++) {
EmbeddedTangentVector z = orthonormalFrame[col][subCol];
// P_x \partial^2 f z
for (size_t row=0; row<nDofs0; row++)
{
typename TargetSpace0::TangentVector semiEmbeddedProduct;
embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
for (size_t row=0; row<nDofs; row++) {
TangentVector semiEmbeddedProduct;
embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
localHessian[_0][_0][row*blocksize0+subRow][col*blocksize0 + subCol] = semiEmbeddedProduct[subRow];
for (int subRow=0; subRow<blocksize; subRow++)
localHessian[row*blocksize+subRow][col*blocksize+subCol] = semiEmbeddedProduct[subRow];
}
}
}
localHessian[_0][_1].setSize(nDofs0*blocksize0, nDofs1*blocksize1);
//////////////////////////////////////////////////////////////////////////////////////
// Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f)
//////////////////////////////////////////////////////////////////////////////////////
// Project embedded gradient onto normal space
std::vector<typename TargetSpace::EmbeddedTangentVector> projectedGradient(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++)
projectedGradient[i] = localSolution[i].projectOntoNormalSpace(localEmbeddedGradient[i]);
for (size_t row=0; row<nDofs; row++) {
for (size_t subRow=0; subRow<blocksize; subRow++) {
EmbeddedTangentVector z = orthonormalFrame[row][subRow];
EmbeddedTangentVector tmp1 = localSolution[row].weingarten(z,projectedGradient[row]);
TangentVector tmp2;
orthonormalFrame[row].mv(tmp1,tmp2);
for (size_t col=0; col<nDofs1; col++)
for (size_t subCol=0; subCol<blocksize; subCol++)
localHessian[row*blocksize+subRow][row*blocksize+subCol] += tmp2[subCol];
}
}
}
/////////////////////////////////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
/////////////////////////////////////////////////////////////////////////////////////
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const
{
using namespace Dune::Indices;
if constexpr (not Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
{
for (size_t subCol=0; subCol<blocksize1; subCol++)
{
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol];
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==1);
return assembleGradientAndHessian(localView,
localConfiguration[_0],
localGradient,
localHessian[_0][_0]);
}
else
{
int rank = Dune::MPIHelper::getCommunication().rank();
// P_x \partial^2 f z
for (size_t row=0; row<nDofs0; row++)
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration);
/////////////////////////////////////////////////////////////////
// Compute the gradient. It is needed to transform the Hessian
// into the correct coordinates.
/////////////////////////////////////////////////////////////////
// TODO: The following fails in the clang-11 C++20 CI jobs
//static_assert(localConfiguration.size()==2, "ADOL-C assembly only implemented for product spaces with two factors");
//! Dimension of the embedding space
using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
constexpr static int embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension;
constexpr static int embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension;
// Compute the actual gradient
size_t nDofs0 = localConfiguration[_0].size();
size_t nDofs1 = localConfiguration[_1].size();
size_t nDoubles = nDofs0*embeddedBlocksize0 + nDofs1*embeddedBlocksize1;
std::vector<double> xp(nDoubles);
int idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++)
for (size_t j=0; j<embeddedBlocksize0; j++)
xp[idx++] = localConfiguration[_0][i].globalCoordinates()[j];
for (size_t i=0; i<localConfiguration[_1].size(); i++)
for (size_t j=0; j<embeddedBlocksize1; j++)
xp[idx++] = localConfiguration[_1][i].globalCoordinates()[j];
// Compute gradient
std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data());
// Copy into Dune type
std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration[_0].size());
std::vector<typename TargetSpace1::EmbeddedTangentVector> localEmbeddedGradient1(localConfiguration[_1].size());
idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++) {
for (size_t j=0; j<embeddedBlocksize0; j++)
localEmbeddedGradient0[i][j] = g[idx++];
// Express gradient in local coordinate system
typename TargetSpace0::TangentVector tmp;
localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],tmp);
for (size_t j=0; j<blocksize0; j++)
localGradient[i*blocksize0+j] = tmp[j];
}
auto offset = localConfiguration[_0].size()*blocksize0;
for (size_t i=0; i<localConfiguration[_1].size(); i++) {
for (size_t j=0; j<embeddedBlocksize1; j++)
localEmbeddedGradient1[i][j] = g[idx++];
// Express gradient in local coordinate system
typename TargetSpace1::TangentVector tmp;
localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],tmp);
for (size_t j=0; j<blocksize1; j++)
localGradient[offset + i*blocksize1 + j] = tmp[j];
}
/////////////////////////////////////////////////////////////////
// Compute Hessian
/////////////////////////////////////////////////////////////////
// We compute the Hessian of the energy functional using the ADOL-C system.
// Since ADOL-C does not know about nonlinear spaces, what we get is actually
// the Hessian of a prolongation of the energy functional into the surrounding
// Euclidean space. To obtain the Riemannian Hessian from this we apply the
// formula described in Absil, Mahoney, Trumpf, "An extrinsic look at the Riemannian Hessian".
// This formula consists of two steps:
// 1) Remove all entries of the Hessian pertaining to the normal space of the
// manifold. In the aforementioned paper this is done by projection onto the
// tangent space. Since we want a matrix that is really smaller (but full rank again),
// we can achieve the same effect by multiplying the embedded Hessian from the left
// and from the right by the matrix of orthonormal frames.
// 2) Add a correction involving the Weingarten map.
//
// This works, and is easy to implement using the ADOL-C "hessian" driver.
// However, here we implement a small shortcut. Computing the embedded Hessian and
// multiplying one side by the orthonormal frame is the same as evaluating the Hessian
// (seen as an operator from R^n to R^n) in the directions of the vectors of the
// orthonormal frame. By luck, ADOL-C can compute the evaluations of the Hessian in
// a given direction directly (in fact, this is also how the "hessian" driver works).
// Since there are less frame vectors than the dimension of the embedding space,
// this reinterpretation allows to reduce the number of calls to ADOL-C.
// In my Cosserat shell tests this reduced assembly time by about 10%.
std::vector<Dune::FieldMatrix<RT,blocksize0,embeddedBlocksize0> > orthonormalFrame0(nDofs0);
for (size_t i=0; i<nDofs0; i++)
orthonormalFrame0[i] = localConfiguration[_0][i].orthonormalFrame();
std::vector<Dune::FieldMatrix<RT,blocksize1,embeddedBlocksize1> > orthonormalFrame1(nDofs1);
for (size_t i=0; i<nDofs1; i++)
orthonormalFrame1[i] = localConfiguration[_1][i].orthonormalFrame();
Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize0> > embeddedHessian00(nDofs0,nDofs0);
Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize1> > embeddedHessian01(nDofs0,nDofs1);
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize0> > embeddedHessian10(nDofs1,nDofs0);
Dune::Matrix<Dune::FieldMatrix<double,blocksize1, embeddedBlocksize1> > embeddedHessian11(nDofs1,nDofs1);
if(adolcScalarMode_) {
std::vector<double> v(nDoubles);
std::vector<double> w(nDoubles);
std::fill(v.begin(), v.end(), 0.0);
size_t nDoubles0 = nDofs0*embeddedBlocksize0; // nDoubles = nDoubles0 + nDoubles1
size_t nDoubles1 = nDofs1*embeddedBlocksize1;
std::fill(v.begin(), v.end(), 0.0);
for (size_t i=0; i<nDofs0 + nDofs1; i++) {
// Evaluate Hessian in the direction of each vector of the orthonormal frame
if (i < nDofs0) { //Upper half
auto i0 = i;
for (int ii0=0; ii0<blocksize0; ii0++) {
for (size_t k0=0; k0<embeddedBlocksize0; k0++) {
v[i0*embeddedBlocksize0 + k0] = orthonormalFrame0[i0][ii0][k0];
}
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
for (size_t j0=0; j0<nDoubles0; j0++) //Upper left
embeddedHessian00[i0][j0/embeddedBlocksize0][ii0][j0%embeddedBlocksize0] = w[j0];
for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
embeddedHessian01[i0][j1/embeddedBlocksize1][ii0][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
// Make v the null vector again
std::fill(&v[i0*embeddedBlocksize0], &v[(i0+1)*embeddedBlocksize0], 0.0);
}
} else { // i = nDofs0 ... nDofs0 + nDofs1 //lower half
auto i1 = i - nDofs0;
for (int ii1=0; ii1<blocksize1; ii1++) {
for (size_t k1=0; k1<embeddedBlocksize1; k1++) {
v[nDoubles0 + i1*embeddedBlocksize1 + k1] = orthonormalFrame1[i1][ii1][k1];
}
int rc= 3;
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
if (rc < 0)
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
for (size_t j0=0; j0<nDoubles0; j0++) // Upper left
embeddedHessian10[i1][j0/embeddedBlocksize0][ii1][j0%embeddedBlocksize0] = w[j0];
for (size_t j1=0; j1<nDoubles1; j1++) //Upper right
embeddedHessian11[i1][j1/embeddedBlocksize1][ii1][j1%embeddedBlocksize1] = w[nDoubles0 + j1];
// Make v the null vector again
std::fill(&v[nDoubles0 + i1*embeddedBlocksize1], &v[nDoubles0 + (i1+1)*embeddedBlocksize1], 0.0);
}
}
}
} else { //ADOL-C vector mode}
int n = nDoubles;
int nDirections = nDofs0 * blocksize0 + nDofs1 * blocksize1;
double* tangent[nDoubles];
for(size_t i=0; i<nDoubles; i++)
tangent[i] = (double*)malloc(nDirections*sizeof(double));
double* rawHessian[nDoubles];
for(size_t i=0; i<nDoubles; i++)
rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
// Initialize directions field with zeros
for (int j=0; j<nDirections; j++)
for (int i=0; i<n; i++)
tangent[i][j] = 0.0;
for (size_t j=0; j<nDofs0*blocksize0; j++)
for (int i=0; i<embeddedBlocksize0; i++)
tangent[(j/blocksize0)*embeddedBlocksize0+i][j] = orthonormalFrame0[j/blocksize0][j%blocksize0][i];
for (size_t j=0; j<nDofs1*blocksize1; j++)
for (int i=0; i<embeddedBlocksize1; i++)
tangent[nDofs0*embeddedBlocksize0 + (j/blocksize1)*embeddedBlocksize1+i][nDofs0*blocksize0 + j] = orthonormalFrame1[j/blocksize1][j%blocksize1][i];
hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type
size_t offset0 = nDofs0*embeddedBlocksize0;
size_t offset1 = nDofs0*blocksize0;
// upper left block
for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
for (size_t j=0; j<nDofs0*blocksize0; j++)
embeddedHessian00[j/blocksize0][i/embeddedBlocksize0][j%blocksize0][i%embeddedBlocksize0] = rawHessian[i][j];
// upper right block
for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
for (size_t j=0; j<nDofs0*blocksize0; j++)
embeddedHessian01[j/blocksize0][i/embeddedBlocksize1][j%blocksize0][i%embeddedBlocksize1] = rawHessian[offset0+i][j];
// lower left block
for(size_t i=0; i<nDofs0*embeddedBlocksize0; i++)
for (size_t j=0; j<nDofs1*blocksize1; j++)
embeddedHessian10[j/blocksize1][i/embeddedBlocksize0][j%blocksize1][i%embeddedBlocksize0] = rawHessian[i][offset1+j];
// lower right block
for(size_t i=0; i<nDofs1*embeddedBlocksize1; i++)
for (size_t j=0; j<nDofs1*blocksize1; j++)
embeddedHessian11[j/blocksize1][i/embeddedBlocksize1][j%blocksize1][i%embeddedBlocksize1] = rawHessian[offset0+i][offset1+j];
for(size_t i=0; i<nDoubles; i++) {
free(rawHessian[i]);
free(tangent[i]);
}
}
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space.
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
// at the Riemannian Hessian".
using namespace Dune::Indices;
localHessian[_0][_0].setSize(nDofs0*blocksize0, nDofs0*blocksize0);
for (size_t col=0; col<nDofs0; col++)
{
for (size_t subCol=0; subCol<blocksize0; subCol++)
{
typename TargetSpace0::TangentVector semiEmbeddedProduct;
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol];
for (int subRow=0; subRow<blocksize0; subRow++)
localHessian[_0][_1][row*blocksize0+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
// P_x \partial^2 f z
for (size_t row=0; row<nDofs0; row++)
{
typename TargetSpace0::TangentVector semiEmbeddedProduct;
embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
localHessian[_0][_0][row*blocksize0+subRow][col*blocksize0 + subCol] = semiEmbeddedProduct[subRow];
}
}
}
}
localHessian[_1][_0].setSize(nDofs1*blocksize1, nDofs0*blocksize0);
localHessian[_0][_1].setSize(nDofs0*blocksize0, nDofs1*blocksize1);
for (size_t col=0; col<nDofs0; col++)
{
for (size_t subCol=0; subCol<blocksize0; subCol++)
for (size_t col=0; col<nDofs1; col++)
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol];
for (size_t subCol=0; subCol<blocksize1; subCol++)
{
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol];
// P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++) {
typename TargetSpace1::TangentVector semiEmbeddedProduct;
embeddedHessian10[row][col].mv(z,semiEmbeddedProduct);
// P_x \partial^2 f z
for (size_t row=0; row<nDofs0; row++)
{
typename TargetSpace0::TangentVector semiEmbeddedProduct;
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
localHessian[_1][_0][row*blocksize1+subRow][col*blocksize0+subCol] = semiEmbeddedProduct[subRow];
for (int subRow=0; subRow<blocksize0; subRow++)
localHessian[_0][_1][row*blocksize0+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
}
}
}
}
localHessian[_1][_1].setSize(nDofs1*blocksize1, nDofs1*blocksize1);
localHessian[_1][_0].setSize(nDofs1*blocksize1, nDofs0*blocksize0);
for (size_t col=0; col<nDofs1; col++)
{
for (size_t subCol=0; subCol<blocksize1; subCol++)
for (size_t col=0; col<nDofs0; col++)
{
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol];
for (size_t subCol=0; subCol<blocksize0; subCol++)
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol];
// P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++)
// P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++) {
typename TargetSpace1::TangentVector semiEmbeddedProduct;
embeddedHessian10[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
localHessian[_1][_0][row*blocksize1+subRow][col*blocksize0+subCol] = semiEmbeddedProduct[subRow];
}
}
}
localHessian[_1][_1].setSize(nDofs1*blocksize1, nDofs1*blocksize1);
for (size_t col=0; col<nDofs1; col++)
{
for (size_t subCol=0; subCol<blocksize1; subCol++)
{
typename TargetSpace1::TangentVector semiEmbeddedProduct;
embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol];
for (int subRow=0; subRow<blocksize1; subRow++)
localHessian[_1][_1][row*blocksize1+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
// P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++)
{
typename TargetSpace1::TangentVector semiEmbeddedProduct;
embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
localHessian[_1][_1][row*blocksize1+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow];
}
}
}
}
//////////////////////////////////////////////////////////////////////////////////////
// Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f)
//////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////
// Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f)
//////////////////////////////////////////////////////////////////////////////////////
// Project embedded gradient onto normal space
std::vector<typename TargetSpace0::EmbeddedTangentVector> projectedGradient0(nDofs0);
for (size_t i=0; i<nDofs0; i++)
projectedGradient0[i] = localConfiguration[_0][i].projectOntoNormalSpace(localEmbeddedGradient0[i]);
// Project embedded gradient onto normal space
std::vector<typename TargetSpace0::EmbeddedTangentVector> projectedGradient0(nDofs0);
for (size_t i=0; i<nDofs0; i++)
projectedGradient0[i] = localConfiguration[_0][i].projectOntoNormalSpace(localEmbeddedGradient0[i]);
std::vector<typename TargetSpace1::EmbeddedTangentVector> projectedGradient1(nDofs1);
for (size_t i=0; i<nDofs1; i++)
projectedGradient1[i] = localConfiguration[_1][i].projectOntoNormalSpace(localEmbeddedGradient1[i]);
std::vector<typename TargetSpace1::EmbeddedTangentVector> projectedGradient1(nDofs1);
for (size_t i=0; i<nDofs1; i++)
projectedGradient1[i] = localConfiguration[_1][i].projectOntoNormalSpace(localEmbeddedGradient1[i]);
// The Weingarten map has only diagonal entries
for (size_t row=0; row<nDofs0; row++)
{
for (size_t subRow=0; subRow<blocksize0; subRow++)
// The Weingarten map has only diagonal entries
for (size_t row=0; row<nDofs0; row++)
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[row][subRow];
typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration[_0][row].weingarten(z,projectedGradient0[row]);
for (size_t subRow=0; subRow<blocksize0; subRow++)
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[row][subRow];
typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration[_0][row].weingarten(z,projectedGradient0[row]);
typename TargetSpace0::TangentVector tmp2;
orthonormalFrame0[row].mv(tmp1,tmp2);
typename TargetSpace0::TangentVector tmp2;
orthonormalFrame0[row].mv(tmp1,tmp2);
for (size_t subCol=0; subCol<blocksize0; subCol++)
localHessian[_0][_0][row*blocksize0+subRow][row*blocksize0+subCol] += tmp2[subCol];
for (size_t subCol=0; subCol<blocksize0; subCol++)
localHessian[_0][_0][row*blocksize0+subRow][row*blocksize0+subCol] += tmp2[subCol];
}
}
}
for (size_t row=0; row<nDofs1; row++)
{
for (size_t subRow=0; subRow<blocksize1; subRow++)
for (size_t row=0; row<nDofs1; row++)
{
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[row][subRow];
typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration[_1][row].weingarten(z,projectedGradient1[row]);
for (size_t subRow=0; subRow<blocksize1; subRow++)
{
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[row][subRow];
typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration[_1][row].weingarten(z,projectedGradient1[row]);
typename TargetSpace1::TangentVector tmp2;
orthonormalFrame1[row].mv(tmp1,tmp2);
typename TargetSpace1::TangentVector tmp2;
orthonormalFrame1[row].mv(tmp1,tmp2);
for (size_t subCol=0; subCol<blocksize1; subCol++)
localHessian[_1][_1][row*blocksize1+subRow][row*blocksize1+subCol] += tmp2[subCol];
for (size_t subCol=0; subCol<blocksize1; subCol++)
localHessian[_1][_1][row*blocksize1+subRow][row*blocksize1+subCol] += tmp2[subCol];
}
}
}
}
}
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localgeodesicfefdstiffness.hh
View file @ 48b9358e
......@@ -6,272 +6,278 @@
#include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
/** \brief Assembles energy gradient and Hessian with finite difference approximations
*/
template<class Basis, class TargetSpace, class field_type=double>
class LocalGeodesicFEFDStiffness
: public LocalGeodesicFEStiffness<Basis,TargetSpace>
namespace Dune::GFE
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
typedef typename TargetSpace::template rebind<field_type>::other ATargetSpace;
/** \brief Assembles energy gradient and Hessian with finite difference approximations
*/
template<class Basis, class TargetSpace, class field_type=double>
class LocalGeodesicFEFDStiffness
: public LocalGeodesicFEStiffness<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
// some other sizes
constexpr static int gridDim = GridView::dimension;
typedef typename TargetSpace::template rebind<field_type>::other ATargetSpace;
public:
// some other sizes
constexpr static int gridDim = GridView::dimension;
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
public:
//! Dimension of the embedding space
constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
LocalGeodesicFEFDStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy)
: localEnergy_(energy)
{}
//! Dimension of the embedding space
constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
/** \brief Compute the energy at the current configuration */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
{
return localEnergy_->energy(localView,localSolution);
}
LocalGeodesicFEFDStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy)
: localEnergy_(energy)
{}
RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
return localEnergy_->energy(localView,coefficients);
}
/** \brief Compute the energy at the current configuration */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
{
return localEnergy_->energy(localView,localSolution);
}
/** \brief Assemble the element gradient of the energy functional
RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
return localEnergy_->energy(localView,coefficients);
}
The default implementation in this class uses a finite difference approximation */
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<field_type>& gradient) const override;
/** \brief Assemble the element gradient of the energy functional
/** \brief Assemble the local tangent matrix and gradient at the current position
The default implementation in this class uses a finite difference approximation */
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<field_type>& gradient) const override;
This implementation uses finite-difference approximations
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
/** \brief Assemble the local tangent matrix and gradient at the current position
/** \brief Assemble the local tangent matrix and gradient at the current position
This implementation uses finite-difference approximations
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
This implementation uses finite-difference approximations
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
/** \brief Assemble the local tangent matrix and gradient at the current position
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
This implementation uses finite-difference approximations
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient) const
{
};
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient) const
{
field_type eps = 1e-6;
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
std::vector<ATargetSpace> localASolution(localSolution.size());
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] = raw[j];
localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble
}
field_type eps = 1e-6;
localGradient.resize(localSolution.size());
std::vector<ATargetSpace> localASolution(localSolution.size());
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] = raw[j];
localASolution[i] = aRaw[i]; // may contain a projection onto M -- needs to be done in adouble
}
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
localGradient.resize(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
// basis vectors of the tangent space of the i-th entry of localSolution
const Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
for (size_t i=0; i<localSolution.size(); i++) {
for (int j=0; j<blocksize; j++) {
// basis vectors of the tangent space of the i-th entry of localSolution
const Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
typename ATargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
forwardCorrection *= eps;
for (int j=0; j<blocksize; j++) {
typename ATargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
backwardCorrection *= -eps;
typename ATargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
forwardCorrection *= eps;
forwardSolution[i] = ATargetSpace::exp(localASolution[i], forwardCorrection);
backwardSolution[i] = ATargetSpace::exp(localASolution[i], backwardCorrection);
typename ATargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
backwardCorrection *= -eps;
field_type foo = (localEnergy_->energy(localView,forwardSolution) - localEnergy_->energy(localView, backwardSolution)) / (2*eps);
forwardSolution[i] = ATargetSpace::exp(localASolution[i], forwardCorrection);
backwardSolution[i] = ATargetSpace::exp(localASolution[i], backwardCorrection);
field_type foo = (localEnergy_->energy(localView,forwardSolution) - localEnergy_->energy(localView, backwardSolution)) / (2*eps);
#ifdef MULTIPRECISION
localGradient[i*blocksize+j] = foo.template convert_to<double>();
localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else
localGradient[i*blocksize+j] = foo;
localGradient[i*blocksize+j] = foo;
#endif
}
}
forwardSolution[i] = localASolution[i];
backwardSolution[i] = localASolution[i];
forwardSolution[i] = localASolution[i];
backwardSolution[i] = localASolution[i];
}
}
}
}
/////////////////////////////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
/////////////////////////////////////////////////////////////////////////////////
template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Number of degrees of freedom for this element
size_t nDofs = localSolution.size();
/////////////////////////////////////////////////////////////////////////////////
// Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
/////////////////////////////////////////////////////////////////////////////////
template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Number of degrees of freedom for this element
size_t nDofs = localSolution.size();
// Clear assemble data
localHessian.setSize(nDofs*blocksize, nDofs*blocksize);
// Clear assemble data
localHessian.setSize(nDofs*blocksize, nDofs*blocksize);
localHessian = 0;
localHessian = 0;
#ifdef MULTIPRECISION
const field_type eps = 1e-10;
const field_type eps = 1e-10;
#else
const field_type eps = 1e-4;
const field_type eps = 1e-4;
#endif
std::vector<ATargetSpace> localASolution(localSolution.size());
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] = raw[j];
localASolution[i] = aRaw[i];
}
std::vector<ATargetSpace> localASolution(localSolution.size());
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] = raw[j];
localASolution[i] = aRaw[i];
}
std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
for (size_t i=0; i<B.size(); i++)
B[i] = localSolution[i].orthonormalFrame();
std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
for (size_t i=0; i<B.size(); i++)
B[i] = localSolution[i].orthonormalFrame();
// Precompute negative energy at the current configuration
// (negative because that is how we need it as part of the 2nd-order fd formula)
field_type centerValue = -localEnergy_->energy(localView, localASolution);
// Precompute negative energy at the current configuration
// (negative because that is how we need it as part of the 2nd-order fd formula)
field_type centerValue = -localEnergy_->energy(localView, localASolution);
// Precompute energy infinitesimal corrections in the directions of the local basis vectors
std::vector<std::array<field_type,blocksize> > forwardEnergy(nDofs);
std::vector<std::array<field_type,blocksize> > backwardEnergy(nDofs);
// Precompute energy infinitesimal corrections in the directions of the local basis vectors
std::vector<std::array<field_type,blocksize> > forwardEnergy(nDofs);
std::vector<std::array<field_type,blocksize> > backwardEnergy(nDofs);
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
minusEpsXi *= -1;
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
minusEpsXi *= -1;
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution;
forwardSolution[i] = ATargetSpace::exp(localASolution[i],epsXi);
backwardSolution[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
forwardSolution[i] = ATargetSpace::exp(localASolution[i],epsXi);
backwardSolution[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
forwardEnergy[i][i2] = localEnergy_->energy(localView, forwardSolution);
backwardEnergy[i][i2] = localEnergy_->energy(localView, backwardSolution);
forwardEnergy[i][i2] = localEnergy_->energy(localView, forwardSolution);
backwardEnergy[i][i2] = localEnergy_->energy(localView, backwardSolution);
}
}
}
}
//////////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
//////////////////////////////////////////////////////////////
localGradient.resize(nDofs*blocksize);
localGradient.resize(nDofs*blocksize);
for (size_t i=0; i<localSolution.size(); i++)
for (int j=0; j<blocksize; j++)
{
field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
for (size_t i=0; i<localSolution.size(); i++)
for (int j=0; j<blocksize; j++)
{
field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
#ifdef MULTIPRECISION
localGradient[i*blocksize+j] = foo.template convert_to<double>();
localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else
localGradient[i*blocksize+j] = foo;
localGradient[i*blocksize+j] = foo;
#endif
}
///////////////////////////////////////////////////////////////////////////
// Compute Riemannian Hesse matrix by finite-difference approximation.
// We loop over the lower left triangular half of the matrix.
// The other half follows from symmetry.
///////////////////////////////////////////////////////////////////////////
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
for (size_t j=0; j<=i; j++) {
for (size_t j2=0; j2<((i==j) ? i2+1 : blocksize); j2++) {
std::vector<ATargetSpace> forwardSolutionXiEta = localASolution;
std::vector<ATargetSpace> backwardSolutionXiEta = localASolution;
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
typename ATargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
if (i==j)
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi+epsEta);
else {
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi);
forwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],epsEta);
}
if (i==j)
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi+minusEpsEta);
else {
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
backwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],minusEpsEta);
}
field_type forwardValue = localEnergy_->energy(localView, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
field_type backwardValue = localEnergy_->energy(localView, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
}
field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
///////////////////////////////////////////////////////////////////////////
// Compute Riemannian Hesse matrix by finite-difference approximation.
// We loop over the lower left triangular half of the matrix.
// The other half follows from symmetry.
///////////////////////////////////////////////////////////////////////////
//#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) {
for (size_t j=0; j<=i; j++) {
for (size_t j2=0; j2<((i==j) ? i2+1 : blocksize); j2++) {
std::vector<ATargetSpace> forwardSolutionXiEta = localASolution;
std::vector<ATargetSpace> backwardSolutionXiEta = localASolution;
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector epsEta = B[j][j2]; epsEta *= eps;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; minusEpsXi *= -1;
typename ATargetSpace::EmbeddedTangentVector minusEpsEta = epsEta; minusEpsEta *= -1;
if (i==j)
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi+epsEta);
else {
forwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],epsXi);
forwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],epsEta);
}
if (i==j)
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi+minusEpsEta);
else {
backwardSolutionXiEta[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
backwardSolutionXiEta[j] = ATargetSpace::exp(localASolution[j],minusEpsEta);
}
field_type forwardValue = localEnergy_->energy(localView, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
field_type backwardValue = localEnergy_->energy(localView, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
#ifdef MULTIPRECISION
localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo.template convert_to<double>();
localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo.template convert_to<double>();
#else
localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo;
localHessian[i*blocksize+i2][j*blocksize+j2] = localHessian[j*blocksize+j2][i*blocksize+i2] = foo;
#endif
}
}
}
}
}
}
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localgeodesicfestiffness.hh
View file @ 48b9358e
......@@ -64,29 +64,30 @@ namespace Dune::GFE
// Type of the local Hessian
using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
};
}
}
} // namespace Impl
template<class Basis, class TargetSpace>
class LocalGeodesicFEStiffness
: public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
{
public:
/** \brief Assemble the local gradient and stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
template<class Basis, class TargetSpace>
class LocalGeodesicFEStiffness
: public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
{
public:
/** \brief Assemble the local gradient and stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localintegralenergy.hh
View file @ 48b9358e
......@@ -2,6 +2,7 @@
#define DUNE_GFE_ASSEMBLERS_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
#include <dune/geometry/quadraturerules.hh>
......@@ -9,8 +10,8 @@
#include <dune/gfe/densities/localdensity.hh>
namespace Dune::GFE {
namespace Dune::GFE
{
/** \brief An energy given as an integral over a density
*
* \tparam Basis The scalar finite element basis used to construct the interpolation rule
......@@ -23,6 +24,7 @@ namespace Dune::GFE {
{
using LocalView = typename Basis::LocalView;
using GridView = typename LocalView::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using DT = typename GridView::Grid::ctype;
using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT;
......@@ -30,10 +32,28 @@ namespace Dune::GFE {
public:
/** \brief Constructor with a Dune::GFE::LocalDensity
/** \brief Constructor from a GFE function as a shared pointer
*
* \param localGFEFunction The geometric finite element function that
* the density will be evaluated on
* \param ld The density that will be integrated
*/
LocalIntegralEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
: localGFEFunction_(localGFEFunction),
localDensity_(ld)
{}
/** \brief Constructor from a GFE function r-value reference
*
* \param localGFEFunction The geometric finite element function that
* the density will be evaluated on
* \param ld The density that will be integrated
*/
LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld)
: localDensity_(ld)
LocalIntegralEnergy(LocalInterpolationRule&& localGFEFunction,
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction))),
localDensity_(ld)
{}
private:
......@@ -44,18 +64,22 @@ namespace Dune::GFE {
{
RT energy = 0;
if constexpr (not Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
{
const auto& localFiniteElement = localView.tree().finiteElement();
LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
// Bind GFE function to the current element
const auto& element = localView.element();
localGFEFunction_->bind(element,localConfiguration);
int quadOrder = (localFiniteElement.type().isSimplex())
? (localFiniteElement.localBasis().order()-1) * 2
: (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
// Bind density to the element
localDensity_->bind(element);
const auto element = localView.element();
// Get a suitable quadrature rule
const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
int quadOrder = (element.type().isSimplex())
? (localGFEOrder-1) * 2
: (localGFEOrder * gridDim - 1) * 2;
const auto& quad = QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
for (auto&& qp : quad)
{
......@@ -70,13 +94,13 @@ namespace Dune::GFE {
{
if (localDensity_->dependsOnDerivative())
{
auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
}
else
{
auto value = localInterpolationRule.evaluate(quadPos);
const auto value = localGFEFunction_->evaluate(quadPos);
typename LocalInterpolationRule::DerivativeType dummyDerivative;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
}
......@@ -86,7 +110,7 @@ namespace Dune::GFE {
if (localDensity_->dependsOnDerivative())
{
typename TargetSpace::CoordinateType dummyValue;
auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
}
......@@ -97,6 +121,11 @@ namespace Dune::GFE {
}
}
}
else
{
// You need a scalar basis or a power basis when calling this method.
std::abort();
}
return energy;
}
......@@ -106,27 +135,26 @@ namespace Dune::GFE {
{
RT energy = 0;
if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold
&& Basis::LocalView::Tree::isComposite
&& gridDim==GridView::dimensionworld) // TODO: Implement the case gridDim!=dimworld
{
static_assert(TargetSpace::size() == 2,
"LocalIntegralEnergy only implemented for product spaces with two factors!");
const auto& element = localView.element();
using namespace Indices;
// composite Basis: grab the finite element of the first child
const auto& localFiniteElement0 = localView.tree().child(_0,0).finiteElement();
const auto& localFiniteElement1 = localView.tree().child(_1,0).finiteElement();
using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
// Bind GFE functions to the current element
const auto& element = localView.element();
std::get<0>(*localGFEFunction_).bind(element, coefficients[_0]);
std::get<1>(*localGFEFunction_).bind(element, coefficients[_1]);
LocalGFEFunctionType0 localGFEFunction0(localFiniteElement0, coefficients[_0]);
LocalGFEFunctionType1 localGFEFunction1(localFiniteElement1, coefficients[_1]);
// Bind density to the element
localDensity_->bind(element);
int quadOrder = (element.type().isSimplex()) ? localFiniteElement0.localBasis().order()
: localFiniteElement0.localBasis().order() * gridDim;
// Get a suitable quadrature rule
const auto feOrder = localView.tree().child(_0,0).finiteElement().localBasis().order();
int quadOrder = (element.type().isSimplex()) ? feOrder : feOrder * gridDim;
const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
......@@ -147,19 +175,19 @@ namespace Dune::GFE {
// A value is needed either directly, or for computing the derivative.
TargetSpace value;
if (localDensity_->dependsOnValue(0) || localDensity_->dependsOnDerivative(0))
value[_0] = localGFEFunction0.evaluate(quadPos);
value[_0] = std::get<0>(*localGFEFunction_).evaluate(quadPos);
if (localDensity_->dependsOnValue(1) || localDensity_->dependsOnDerivative(1))
value[_1] = localGFEFunction1.evaluate(quadPos);
value[_1] = std::get<1>(*localGFEFunction_).evaluate(quadPos);
// Compute the derivatives of the interpolation function factors
typename LocalGFEFunctionType0::DerivativeType derivative0;
typename LocalGFEFunctionType1::DerivativeType derivative1;
typename std::tuple_element_t<0, LocalInterpolationRule>::DerivativeType derivative0;
typename std::tuple_element_t<1, LocalInterpolationRule>::DerivativeType derivative1;
if (localDensity_->dependsOnDerivative(0))
derivative0 = localGFEFunction0.evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
derivative0 = std::get<0>(*localGFEFunction_).evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
if (localDensity_->dependsOnDerivative(1))
derivative1 = localGFEFunction1.evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
derivative1 = std::get<1>(*localGFEFunction_).evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
// Copy the two derivatives into a joint matrix object
// TODO: I am not sure about this. May the densities should get the
......@@ -177,12 +205,19 @@ namespace Dune::GFE {
derivative);
}
}
else
DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis or gridDim!=dimworld");
return energy;
}
protected:
const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> > localDensity_;
// The value and derivative of this function are evaluated at the quadrature points,
// and given to the density.
const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> > localDensity_;
};
} // namespace Dune::GFE
......