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 (90)
Hide whitespace changes
Inline Side-by-side
Showing
with 1645 additions and 1684 deletions
.git-blame-ignore-revs
View file @ 48b9358e
# First use of uncrustify, with the vanilla config from dune-project.org # First use of uncrustify, with the vanilla config from dune-project.org
d13f34469065230dd86e63733e4dc8d1f775acfb 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 ...@@ -7,9 +7,10 @@ before_script: &before
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-vtk.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.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.dune-project.org/extensions/dune-gmsh4
- duneci-install-module https://gitlab.mn.tu-dresden.de/iwr/dune-gmsh4
# Tests with the 2.9 release. That's the one in Debian bookworm
#--------------------------------------------------------------------
dune:2.9 gcc: dune:2.9 gcc:
variables: variables:
DUNECI_BRANCH: releases/2.9 DUNECI_BRANCH: releases/2.9
...@@ -18,6 +19,76 @@ dune:2.9 gcc: ...@@ -18,6 +19,76 @@ dune:2.9 gcc:
- *before - *before
script: duneci-standard-test 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: dune:git gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20 image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script: before_script:
...@@ -30,24 +101,44 @@ dune:git clang-11 C++20: ...@@ -30,24 +101,44 @@ dune:git clang-11 C++20:
- *before - *before
script: duneci-standard-test 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 image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script: before_script:
- *before - *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/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.dune-project.org/extensions/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-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.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 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 image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-11-20
before_script: before_script:
- *before - *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/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.dune-project.org/extensions/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-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.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 script: duneci-standard-test
# Check for spelling mistakes in text # Check for spelling mistakes in text
......
CHANGELOG.md
View file @ 48b9358e
# Master # 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 - Building `dune-gfe` now requires Dune 2.9 or newer. Older configurations
have not gotten CI-tested for quite a while, anyway. have not gotten CI-tested for quite a while, anyway.
......
CMakeLists.txt
View file @ 48b9358e
...@@ -17,6 +17,13 @@ include(DuneMacros) ...@@ -17,6 +17,13 @@ include(DuneMacros)
# start a dune project with information from dune.module # start a dune project with information from dune.module
dune_project() 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() dune_enable_all_packages()
add_subdirectory("src") add_subdirectory("src")
add_subdirectory("dune") add_subdirectory("dune")
......
dune.module
View file @ 48b9358e
...@@ -4,8 +4,8 @@ ...@@ -4,8 +4,8 @@
#Name of the module #Name of the module
Module: dune-gfe Module: dune-gfe
Version: svn Version: 2.11-git
Maintainer: oliver.sander@tu-dresden.de Maintainer: oliver.sander@tu-dresden.de
#depending on #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 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-vtk dune-curvedgrid Suggests: dune-foamgrid dune-parmg dune-curvedgrid dune-elasticity
dune/gfe/CMakeLists.txt
View file @ 48b9358e
...@@ -4,20 +4,12 @@ install(FILES ...@@ -4,20 +4,12 @@ install(FILES
cosseratstrain.hh cosseratstrain.hh
cosseratvtkreader.hh cosseratvtkreader.hh
cosseratvtkwriter.hh cosseratvtkwriter.hh
embeddedglobalgfefunction.hh
filereader.hh filereader.hh
geodesicdifference.hh geodesicdifference.hh
geodesicfefunctionadaptor.hh geodesicfefunctionadaptor.hh
gfedifferencenormsquared.hh gfedifferencenormsquared.hh
globalgfefunction.hh
gramschmidtsolver.hh gramschmidtsolver.hh
interpolationderivatives.hh
linearalgebra.hh linearalgebra.hh
localgeodesicfefunction.hh
localgfetestfunctionbasis.hh
localprojectedfefunction.hh
localquickanddirtyfefunction.hh
localtangentfefunction.hh
maxnormtrustregion.hh maxnormtrustregion.hh
mixedriemanniantrsolver.cc mixedriemanniantrsolver.cc
mixedriemanniantrsolver.hh mixedriemanniantrsolver.hh
...@@ -42,3 +34,4 @@ install(FILES ...@@ -42,3 +34,4 @@ install(FILES
add_subdirectory("spaces") add_subdirectory("spaces")
add_subdirectory("assemblers") add_subdirectory("assemblers")
add_subdirectory("densities") add_subdirectory("densities")
add_subdirectory("functions")
dune/gfe/assemblers/CMakeLists.txt
View file @ 48b9358e
#install headers #install headers
install(FILES install(FILES
cosseratenergystiffness.hh
cosseratrodenergy.hh cosseratrodenergy.hh
geodesicfeassembler.hh geodesicfeassembler.hh
geodesicfeassemblerwrapper.hh geodesicfeassemblerwrapper.hh
......
dune/gfe/assemblers/cosseratenergystiffness.hh deleted 100644 → 0
View file @ 2c964dc9
#ifndef COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
#define COSSERAT_ENERGY_LOCAL_STIFFNESS_HH
#define DONT_USE_CURL
//#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
#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>
#include <dune/gfe/densities/bulkcosseratdensity.hh>
#include <dune/gfe/densities/planarcosseratshelldensity.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>
/** \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;
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters)
: bulkDensity_(parameters),
planarCosseratShellDensity_(parameters)
{}
/** \brief Constructor with material parameters and external loads
* \param parameters The material parameters
* \param neumannBoundary The part of the boundary where a surface load is applied
* \param neumannFunction Surface load density
* \param volumeLoad Volume load density
*/
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)
: bulkDensity_(parameters),
planarCosseratShellDensity_(parameters),
neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction),
volumeLoad_(volumeLoad)
{}
/** \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;
/** The energy density if the grid is three-dimensional */
Dune::GFE::BulkCosseratDensity<Dune::FieldVector<double,3>,field_type> bulkDensity_;
/** The energy density if the grid is two-dimensional */
Dune::GFE::PlanarCosseratShellDensity<Dune::FieldVector<double,2>,field_type> planarCosseratShellDensity_;
/** \brief The Neumann boundary */
const BoundaryPatch<GridView>* neumannBoundary_ = nullptr;
/** \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 constexpr (gridDim==2) {
energy += weight * planarCosseratShellDensity_(quadPos,
value.globalCoordinates(),
derivative);
} else if constexpr (gridDim==3) {
energy += weight * bulkDensity_(quadPos,
value.globalCoordinates(),
derivative);
} 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);
TargetSpace value;
value[_0] = deformationValue;
value[_1] = orientationValue;
// 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]);
// The derivative of both factors together
typename Dune::FieldMatrix<field_type,deformationDerivative.N() + orientationDerivative.N(), gridDim> derivative(0);
size_t rowCount = 0;
for (auto&& row : deformationDerivative)
derivative[rowCount++] = row;
for (auto&& row : orientationDerivative)
derivative[rowCount++] = row;
/////////////////////////////////////////////////////////
// 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 constexpr (gridDim==2) {
energy += weight * planarCosseratShellDensity_(quadPos,
value.globalCoordinates(),
derivative);
} else if constexpr (gridDim==3) {
energy += weight * bulkDensity_(quadPos,
value.globalCoordinates(),
derivative);
} 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 @@ ...@@ -2,6 +2,7 @@
#define DUNE_GFE_COSSERATRODENERGY_HH #define DUNE_GFE_COSSERATRODENERGY_HH
#include <array> #include <array>
#include <memory>
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/version.hh> #include <dune/common/version.hh>
...@@ -18,9 +19,10 @@ ...@@ -18,9 +19,10 @@
#include <dune/gfe/spaces/realtuple.hh> #include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.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 class CosseratRodEnergy
: public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > > : public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
{ {
...@@ -42,6 +44,14 @@ namespace Dune::GFE { ...@@ -42,6 +44,14 @@ namespace Dune::GFE {
public: 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 /** \brief The stress-free configuration
The number type cannot be RT, because RT become `adouble` when The number type cannot be RT, because RT become `adouble` when
...@@ -67,9 +77,25 @@ namespace Dune::GFE { ...@@ -67,9 +77,25 @@ namespace Dune::GFE {
GridView gridView_; GridView gridView_;
//! Constructor //! 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) 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), A_(A),
gridView_(gridView) gridView_(gridView)
{} {}
...@@ -80,9 +106,39 @@ namespace Dune::GFE { ...@@ -80,9 +106,39 @@ namespace Dune::GFE {
\param E Young's modulus \param E Young's modulus
\param nu Poisson number \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) 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 // shear modulus
double G = E/(2+2*nu); double G = E/(2+2*nu);
...@@ -172,21 +228,18 @@ namespace Dune::GFE { ...@@ -172,21 +228,18 @@ namespace Dune::GFE {
}; };
template<class Basis, class LocalInterpolationRule, class RT> template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>:: RT CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
energy(const typename Basis::LocalView& localView, energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localCoefficients) const const std::vector<TargetSpace>& localCoefficients) const
{ {
const auto& localFiniteElement = localView.tree().finiteElement();
LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
const auto& element = localView.element(); const auto& element = localView.element();
localGFEFunction_->bind(element, localCoefficients);
RT energy = 0; RT energy = 0;
std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView); 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; localReferenceConfiguration_->bind(element, localReferenceCoefficients);
InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
// /////////////////////////////////////////////////////////////////////////////// // ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate // The following two loops are a reduced integration scheme. We integrate
...@@ -205,10 +258,10 @@ namespace Dune::GFE { ...@@ -205,10 +258,10 @@ namespace Dune::GFE {
double weight = shearingQuad[pt].weight() * integrationElement; double weight = shearingQuad[pt].weight() * integrationElement;
auto strain = getStrain(localConfiguration, element, quadPos); auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain // The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos); auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
for (int i=0; i<3; i++) for (int i=0; i<3; i++)
energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]); energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
...@@ -225,10 +278,10 @@ namespace Dune::GFE { ...@@ -225,10 +278,10 @@ namespace Dune::GFE {
double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos); double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
auto strain = getStrain(localConfiguration, element, quadPos); auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain // The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos); auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
// Part II: the bending and twisting energy // Part II: the bending and twisting energy
for (int i=0; i<3; i++) for (int i=0; i<3; i++)
...@@ -240,9 +293,9 @@ namespace Dune::GFE { ...@@ -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> template <class ReboundLocalInterpolationRule>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>:: auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const ReboundLocalInterpolationRule& localInterpolation, getStrain(const ReboundLocalInterpolationRule& localInterpolation,
const Entity& element, const Entity& element,
const FieldVector<double,1>& pos) const const FieldVector<double,1>& pos) const
...@@ -285,9 +338,9 @@ namespace Dune::GFE { ...@@ -285,9 +338,9 @@ namespace Dune::GFE {
return strain; return strain;
} }
template<class Basis, class LocalInterpolationRule, class RT> template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class Number> 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, getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element, const Entity& element,
const FieldVector<double, 1>& pos) const const FieldVector<double, 1>& pos) const
...@@ -308,8 +361,8 @@ namespace Dune::GFE { ...@@ -308,8 +361,8 @@ namespace Dune::GFE {
return stress; return stress;
} }
template<class Basis, class LocalInterpolationRule, class RT> template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>:: void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol, getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const BlockVector<FieldVector<double, blocksize> >& strain) const
{ {
...@@ -363,8 +416,8 @@ namespace Dune::GFE { ...@@ -363,8 +416,8 @@ namespace Dune::GFE {
} }
} }
template<class Basis, class LocalInterpolationRule, class RT> template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>:: void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol, getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const BlockVector<FieldVector<double, blocksize> >& stress) const
{ {
...@@ -386,9 +439,9 @@ namespace Dune::GFE { ...@@ -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> template <class PatchGridView>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>:: auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary, getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const 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 @@ ...@@ -10,269 +10,273 @@
#include <dune/solvers/common/wrapownshare.hh> #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; namespace Dune::GFE
typedef typename Basis::GridView GridView; {
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
//! Dimension of the grid. /** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
constexpr static int gridDim = GridView::dimension; */
template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
//! Dimension of a tangent space using field_type = typename TargetSpace::field_type;
constexpr static int blocksize = TargetSpace::TangentVector::dimension; typedef typename Basis::GridView GridView;
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
//! //! Dimension of the grid.
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock; 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 protected:
std::shared_ptr<LocalStiffness> localStiffness_;
public: //! The global basis
const Basis basis_;
/** \brief Constructor for a given grid */ //! The local stiffness operator
GeodesicFEAssembler(const Basis& basis) std::shared_ptr<LocalStiffness> localStiffness_;
: basis_(basis)
{}
/** \brief Constructor for a given grid */ public:
template <class LocalStiffnessT>
GeodesicFEAssembler(const Basis& basis,
LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
/** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */ /** \brief Constructor for a given grid */
template <class LocalStiffnessT> GeodesicFEAssembler(const Basis& basis)
void setLocalStiffness(LocalStiffnessT&& localStiffness) { : basis_(basis)
localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)); {}
}
/** \brief Get the local stiffness operator. */
const LocalStiffness& getLocalStiffness() const {
return *localStiffness_;
}
/** \brief Get the local stiffness operator. */ /** \brief Constructor for a given grid */
LocalStiffness& getLocalStiffness() { template <class LocalStiffnessT>
return *localStiffness_; GeodesicFEAssembler(const Basis& basis,
} LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
/** \brief Get the basis. */ /** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
const Basis& getBasis() const { template <class LocalStiffnessT>
return basis_; 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 /** \brief Get the local stiffness operator. */
* const LocalStiffness& getLocalStiffness() const {
* This is more efficient than computing them separately, because you need the gradient return *localStiffness_;
* 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 Assemble the gradient */ /** \brief Get the local stiffness operator. */
virtual void assembleGradient(const std::vector<TargetSpace>& sol, LocalStiffness& getLocalStiffness() {
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const; return *localStiffness_;
}
/** \brief Compute the energy of a deformation state */ /** \brief Get the basis. */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const; const Basis& getBasis() const {
return basis_;
}
//protected: /** \brief Assemble the tangent stiffness matrix and the functional gradient together
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const; *
* 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> }; // end class
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
auto n = basis_.size();
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 auto n = basis_.size();
localView.bind(element);
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); const auto& lfe = localView.tree().finiteElement();
auto jIdx = localView.index(j);
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()); Dune::MatrixIndexSet neighborsPerVertex;
gradient = 0; 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)) hessian = 0;
{
localView.bind(element);
const int numOfBaseFct = localView.tree().size(); gradient.resize(sol.size());
gradient = 0;
// Extract local solution // A view on the FE basis on a single element
std::vector<TargetSpace> localSolution(numOfBaseFct); auto localView = basis_.localView();
for (int i=0; i<numOfBaseFct; i++) for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
localSolution[i] = sol[localView.index(i)]; {
localView.bind(element);
std::vector<double> localGradient(numOfBaseFct*blocksize); const int numOfBaseFct = localView.tree().size();
Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
// setup local matrix and gradient // Extract local solution
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian); 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++) auto row = localView.index(i);
for (std::size_t jj=0; jj<blocksize; jj++)
hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
} for (int j=0; j<numOfBaseFct; j++ ) {
}
// Add local gradient to global gradient auto col = localView.index(j);
for (int i=0; i<numOfBaseFct; i++)
for (int j=0; j<blocksize; j++)
gradient[localView.index(i)][j] += localGradient[i*blocksize + 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> // Add local gradient to global gradient
void GeodesicFEAssembler<Basis,TargetSpace>:: for (int i=0; i<numOfBaseFct; i++)
assembleGradient(const std::vector<TargetSpace>& sol, for (int j=0; j<blocksize; j++)
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
{
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size()); }
grad = 0;
// A view on the FE basis on a single element }
auto localView = basis_.localView();
// Loop over all elements template <class Basis, class TargetSpace>
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior)) 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 grad.resize(sol.size());
const auto nDofs = localView.tree().size(); grad = 0;
// Extract local solution // A view on the FE basis on a single element
std::vector<TargetSpace> localSolution(nDofs); auto localView = basis_.localView();
for (size_t i=0; i<nDofs; i++) // Loop over all elements
localSolution[i] = sol[localView.index(i)]; for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
// Assemble local gradient // The number of degrees of freedom of the current element
std::vector<double> localGradient(nDofs*blocksize); 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 i=0; i<nDofs; i++) localSolution[i] = sol[localView.index(i)];
for (size_t j=0; j<blocksize; j++)
grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
} // Assemble local gradient
std::vector<double> localGradient(nDofs*blocksize);
localStiffness_->assembleGradient(localView, localSolution, localGradient);
template <class Basis, class TargetSpace> // Add to global gradient
double GeodesicFEAssembler<Basis, TargetSpace>:: for (size_t i=0; i<nDofs; i++)
computeEnergy(const std::vector<TargetSpace>& sol) const for (size_t j=0; j<blocksize; j++)
{ grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
double energy = 0; }
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 template <class Basis, class TargetSpace>
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior)) 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 if (sol.size() != basis_.size())
size_t nDofs = localView.tree().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++) // Loop over all elements
localSolution[i] = sol[localView.index(i)[0]]; 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 #endif
dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
View file @ 48b9358e
...@@ -5,7 +5,8 @@ ...@@ -5,7 +5,8 @@
#include <dune/common/tuplevector.hh> #include <dune/common/tuplevector.hh>
#include <dune/gfe/spaces/productmanifold.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 /** \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 { ...@@ -13,9 +14,8 @@ namespace Dune::GFE {
*/ */
template <class Basis, class ScalarBasis, class TargetSpace> template <class Basis, class ScalarBasis, class TargetSpace>
class class GeodesicFEAssemblerWrapper
GeodesicFEAssemblerWrapper { {
using MixedSpace0 = std::tuple_element_t<0,TargetSpace>; using MixedSpace0 = std::tuple_element_t<0,TargetSpace>;
using MixedSpace1 = std::tuple_element_t<1,TargetSpace>; using MixedSpace1 = std::tuple_element_t<1,TargetSpace>;
...@@ -72,7 +72,8 @@ namespace Dune::GFE { ...@@ -72,7 +72,8 @@ namespace Dune::GFE {
auto splitVector(const std::vector<TargetSpace>& sol) const; auto splitVector(const std::vector<TargetSpace>& sol) const;
std::unique_ptr<MatrixType> hessianMixed_; std::unique_ptr<MatrixType> hessianMixed_;
}; // end class }; // end class
} //end namespace
} // namespace Dune::GFE
template <class Basis, class ScalarBasis, class TargetSpace> template <class Basis, class ScalarBasis, class TargetSpace>
......
dune/gfe/assemblers/l2distancesquaredenergy.hh
View file @ 48b9358e
...@@ -5,86 +5,94 @@ ...@@ -5,86 +5,94 @@
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/assemblers/localenergy.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 namespace Dune::GFE
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{ {
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
// some other sizes template<class Basis, class TargetSpace>
constexpr static int gridDim = GridView::dimension; 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 public:
std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
/** \brief Assemble the energy for a single element */ // This is the function that we are computing the L2-distance to
RT energy (const typename Basis::LocalView& localView, std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
const std::vector<TargetSpace>& localSolution) const override
{ /** \brief Assemble the energy for a single element */
RT energy = 0; 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; RT energy = 0;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
const auto element = localView.element(); typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
auto localOrigin = localFunction(*origin_); LocalGFEFunctionType localGeodesicFEFunction(localView.globalBasis());
localOrigin.bind(element);
// Just guessing an appropriate quadrature order const auto element = localView.element();
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim; 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++) // Just guessing an appropriate quadrature order
{ const auto& localFiniteElement = localView.tree().finiteElement();
// Local position of the quadrature point auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
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 weight = quad[pt].weight() * integrationElement;
auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originValue = localOrigin(quadPos);
// The derivative of the 'origin' function // The function value
// First: as function defined on the reference element auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originReferenceDerivative = derivative(localOrigin)(quadPos); auto originValue = localOrigin(quadPos);
// The derivative of the function defined on the actual element // The derivative of the 'origin' function
auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos); // First: as function defined on the reference element
auto originReferenceDerivative = derivative(localOrigin)(quadPos);
double weightFactor = originDerivative.frobenius_norm(); // The derivative of the function defined on the actual element
// un-comment the following line to switch off the weight factor auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
//weightFactor = 1.0;
// Add the local energy density double weightFactor = originDerivative.frobenius_norm();
energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value); // 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, } // namespace Dune::GFE
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
#endif #endif
dune/gfe/assemblers/localenergy.hh
View file @ 48b9358e
...@@ -7,9 +7,9 @@ ...@@ -7,9 +7,9 @@
#include <dune/gfe/spaces/productmanifold.hh> #include <dune/gfe/spaces/productmanifold.hh>
namespace Dune { namespace Dune::GFE
{
namespace GFE:: Impl namespace Impl
{ {
/** \brief A class exporting container types for coefficient sets /** \brief A class exporting container types for coefficient sets
* *
...@@ -34,43 +34,40 @@ namespace Dune { ...@@ -34,43 +34,40 @@ namespace Dune {
using Coefficients = std::vector<ProductManifold<Factors...> >; using Coefficients = std::vector<ProductManifold<Factors...> >;
using CompositeCoefficients = TupleVector<std::vector<Factors>... >; using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
}; };
} } // namespace Impl
namespace GFE {
/** \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 /** \brief Base class for energies defined by integrating over one grid element */
* template<class Basis, class TargetSpace>
* \param localView Local view specifying the current element and the FE space there class LocalEnergy
* \param coefficients The coefficients of a FE function on the current element {
*/ public:
virtual RT using RT = typename TargetSpace::ctype;
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers /** \brief Compute the energy
* for each factor *
*/ * \param localView Local view specifying the current element and the FE space there
virtual RT * \param coefficients The coefficients of a FE function on the current element
energy (const typename Basis::LocalView& localView, */
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0; virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** Empty virtual default destructor /** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* * for each factor
* To allow proper destruction of derived classes through a base class pointer */
*/ virtual RT
virtual ~LocalEnergy() = default; 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 #endif // DUNE_GFE_LOCALENERGY_HH
dune/gfe/assemblers/localfirstordermodel.hh
View file @ 48b9358e
...@@ -20,6 +20,6 @@ namespace Dune::GFE ...@@ -20,6 +20,6 @@ namespace Dune::GFE
}; };
} // namespace Dune } // namespace Dune::GFE
#endif // DUNE_GFE_LOCALFIRSTORDERMODEL_HH #endif // DUNE_GFE_LOCALFIRSTORDERMODEL_HH
dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
View file @ 48b9358e
...@@ -14,145 +14,83 @@ ...@@ -14,145 +14,83 @@
#include <dune/gfe/assemblers/localgeodesicfestiffness.hh> #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
#include <dune/gfe/spaces/productmanifold.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_; namespace Dune::GFE
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
{ {
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; /** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
} */
template<class Basis, class TargetSpace>
class LocalGeodesicFEADOLCStiffness
template <class Basis, class TargetSpace> : public LocalGeodesicFEStiffness<Basis,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)
{ {
// TODO: The following fails in the clang-11 C++20 CI jobs // grid types
//static_assert(localConfiguration.size()==1); typedef typename Basis::GridView GridView;
return energy(localView, localConfiguration[_0]); typedef typename GridView::ctype DT;
} typedef typename TargetSpace::ctype RT;
else 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; double pureEnergy;
int rank = localView.globalBasis().gridView().comm().rank();
std::vector<ATargetSpace> localASolution(localSolution.size());
trace_on(rank); trace_on(rank);
adouble energy = 0; adouble energy = 0;
...@@ -167,374 +105,196 @@ energy(const typename Basis::LocalView& localView, ...@@ -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 // 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.) // (Presumably because several independent variables use the same memory location.)
std::vector<typename ATargetSpace0::CoordinateType> aRaw0(localConfiguration[_0].size()); std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localConfiguration[_0].size(); i++) { for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace0::CoordinateType raw = localConfiguration[_0][i].globalCoordinates(); typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++) for (size_t j=0; j<raw.size(); j++)
aRaw0[i][j] <<= raw[j]; aRaw[i][j] <<= raw[j];
localAConfiguration[_0][i] = aRaw0[i]; // may contain a projection onto M -- needs to be done in adouble localASolution[i] = aRaw[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
} }
try { try {
energy = localEnergy_->energy(localView,localASolution);
energy = localEnergy_->energy(localView, localAConfiguration);
} catch (Dune::Exception &e) { } catch (Dune::Exception &e) {
trace_off(); trace_off();
throw e; throw e;
} }
energy >>= pureEnergy; energy >>= pureEnergy;
trace_off(); trace_off();
return pureEnergy; 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; template <class Basis, class TargetSpace>
MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data())); typename LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::RT
if (rc < 0) LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!"); 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++) int rank = Dune::MPIHelper::getCommunication().rank();
embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
// Make v the null vector again if constexpr (not Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::isProductManifold)
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++)
{ {
for (int i=0; i<n; i++) // TODO: The following fails in the clang-11 C++20 CI jobs
tangent[i][j] = 0.0; //static_assert(localConfiguration.size()==1);
return energy(localView, localConfiguration[_0]);
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); else
{
// Copy Hessian into Dune data type using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
for(size_t i=0; i<nDoubles; i++) using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
for (int j=0; j<nDirections; j++)
embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j]; using ATargetSpace0 = typename TargetSpace0::template rebind<adouble>::other;
using ATargetSpace1 = typename TargetSpace1::template rebind<adouble>::other;
for(size_t i=0; i<nDoubles; i++) {
free(rawHessian[i]); Dune::TupleVector<std::vector<ATargetSpace0>, std::vector<ATargetSpace1> > localAConfiguration;
free(tangent[i]); localAConfiguration[_0].resize(localConfiguration[_0].size());
} localAConfiguration[_1].resize(localConfiguration[_1].size());
}
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded double pureEnergy;
// isometrically in a Euclidean space. trace_on(rank);
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
// at the Riemannian Hessian". adouble energy = 0;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector; // The following loop is not quite intuitive: we copy the localSolution into an
typedef typename TargetSpace::TangentVector TangentVector; // array of FieldVector<double>, go from there to FieldVector<adouble> and
// only then to ATargetSpace.
localHessian.setSize(nDofs*blocksize,nDofs*blocksize); // Rationale: The constructor/assignment-from-vector of TargetSpace frequently
// contains a projection onto the manifold from the surrounding Euclidean space.
for (size_t col=0; col<nDofs; col++) { // 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 energy = localEnergy_->energy(localView, localAConfiguration);
for (size_t row=0; row<nDofs; row++) {
TangentVector semiEmbeddedProduct;
embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize; subRow++) } catch (Dune::Exception &e) {
localHessian[row*blocksize+subRow][col*blocksize+subCol] = semiEmbeddedProduct[subRow]; 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 template <class Basis, class TargetSpace>
std::vector<typename TargetSpace::EmbeddedTangentVector> projectedGradient(localSolution.size()); void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
for (size_t i=0; i<localSolution.size(); i++) assembleGradient(const typename Basis::LocalView& localView,
projectedGradient[i] = localSolution[i].projectOntoNormalSpace(localEmbeddedGradient[i]); 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]; // Compute gradient
EmbeddedTangentVector tmp1 = localSolution[row].weingarten(z,projectedGradient[row]); std::vector<double> g(nDoubles);
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
TangentVector tmp2; // Copy into Dune type
orthonormalFrame[row].mv(tmp1,tmp2); std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
for (size_t subCol=0; subCol<blocksize; subCol++) idx=0;
localHessian[row*blocksize+subRow][row*blocksize+subCol] += tmp2[subCol]; 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. // 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 // Compute the gradient. It is needed to transform the Hessian
// into the correct coordinates. // 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 // Compute the actual gradient
size_t nDofs0 = localConfiguration[_0].size(); constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
size_t nDofs1 = localConfiguration[_1].size(); size_t nDofs = localSolution.size();
size_t nDoubles = nDofs0*embeddedBlocksize0 + nDofs1*embeddedBlocksize1; size_t nDoubles = nDofs*embeddedBlocksize;
std::vector<double> xp(nDoubles); std::vector<double> xp(nDoubles);
int idx=0; int idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++) for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize0; j++) for (size_t j=0; j<embeddedBlocksize; j++)
xp[idx++] = localConfiguration[_0][i].globalCoordinates()[j]; xp[idx++] = localSolution[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 // Compute gradient
std::vector<double> g(nDoubles); 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 // Copy into Dune type
std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration[_0].size()); std::vector<typename TargetSpace::EmbeddedTangentVector> localEmbeddedGradient(localSolution.size());
std::vector<typename TargetSpace1::EmbeddedTangentVector> localEmbeddedGradient1(localConfiguration[_1].size());
idx=0; idx=0;
for (size_t i=0; i<localConfiguration[_0].size(); i++) { for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<embeddedBlocksize0; j++) for (size_t j=0; j<embeddedBlocksize; j++)
localEmbeddedGradient0[i][j] = g[idx++]; localEmbeddedGradient[i][j] = g[idx++];
// Express gradient in local coordinate system // Express gradient in local coordinate system
typename TargetSpace0::TangentVector tmp; for (size_t i=0; i<nDofs; i++) {
localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],tmp); typename TargetSpace::TangentVector tmp;
localSolution[i].orthonormalFrame().mv(localEmbeddedGradient[i],tmp);
for (size_t j=0; j<blocksize0; j++) for (size_t j=0; j<blocksize; j++)
localGradient[i*blocksize0+j] = tmp[j]; localGradient[i*blocksize+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];
} }
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
...@@ -564,81 +324,40 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView, ...@@ -564,81 +324,40 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// this reinterpretation allows to reduce the number of calls to ADOL-C. // this reinterpretation allows to reduce the number of calls to ADOL-C.
// In my Cosserat shell tests this reduced assembly time by about 10%. // 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++) for (size_t i=0; i<nDofs; i++)
orthonormalFrame0[i] = localConfiguration[_0][i].orthonormalFrame(); orthonormalFrame[i] = localSolution[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,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
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_) { if (adolcScalarMode_) {
std::vector<double> v(nDoubles); std::vector<double> v(nDoubles);
std::vector<double> w(nDoubles); std::vector<double> w(nDoubles);
std::fill(v.begin(), v.end(), 0.0); std::fill(v.begin(), v.end(), 0.0);
size_t nDoubles0 = nDofs0*embeddedBlocksize0; // nDoubles = nDoubles0 + nDoubles1 for (size_t i=0; i<nDofs; i++)
size_t nDoubles1 = nDofs1*embeddedBlocksize1; for (int ii=0; ii<blocksize; ii++)
{
std::fill(v.begin(), v.end(), 0.0); // Evaluate Hessian in the direction of each vector of the orthonormal frame
for (size_t k=0; k<embeddedBlocksize; k++)
for (size_t i=0; i<nDofs0 + nDofs1; i++) { v[i*embeddedBlocksize + k] = orthonormalFrame[i][ii][k];
// 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 int rc= 3;
embeddedHessian10[i1][j0/embeddedBlocksize0][ii1][j0%embeddedBlocksize0] = w[j0]; 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 for (size_t j=0; j<nDoubles; j++)
embeddedHessian11[i1][j1/embeddedBlocksize1][ii1][j1%embeddedBlocksize1] = w[nDoubles0 + j1]; embeddedHessian[i][j/embeddedBlocksize][ii][j%embeddedBlocksize] = w[j];
// Make v the null vector again // Make v the null vector again
std::fill(&v[nDoubles0 + i1*embeddedBlocksize1], &v[nDoubles0 + (i1+1)*embeddedBlocksize1], 0.0); 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 n = nDoubles;
int nDirections = nDofs0 * blocksize0 + nDofs1 * blocksize1; int nDirections = nDofs * blocksize;
double* tangent[nDoubles]; double* tangent[nDoubles];
for(size_t i=0; i<nDoubles; i++) for(size_t i=0; i<nDoubles; i++)
tangent[i] = (double*)malloc(nDirections*sizeof(double)); tangent[i] = (double*)malloc(nDirections*sizeof(double));
...@@ -647,182 +366,469 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView, ...@@ -647,182 +366,469 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
for(size_t i=0; i<nDoubles; i++) for(size_t i=0; i<nDoubles; i++)
rawHessian[i] = (double*)malloc(nDirections*sizeof(double)); rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
// Initialize directions field with zeros
for (int j=0; j<nDirections; j++) for (int j=0; j<nDirections; j++)
{
for (int i=0; i<n; i++) for (int i=0; i<n; i++)
tangent[i][j] = 0.0; tangent[i][j] = 0.0;
for (size_t j=0; j<nDofs0*blocksize0; j++) for (int i=0; i<embeddedBlocksize; i++)
for (int i=0; i<embeddedBlocksize0; i++) tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][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); hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type // Copy Hessian into Dune data type
size_t offset0 = nDofs0*embeddedBlocksize0; for(size_t i=0; i<nDoubles; i++)
size_t offset1 = nDofs0*blocksize0; for (int j=0; j<nDirections; j++)
embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
// 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(size_t i=0; i<nDoubles; i++) {
free(rawHessian[i]); free(rawHessian[i]);
free(tangent[i]); free(tangent[i]);
} }
} }
// From this, compute the Hessian with respect to the manifold (which we assume here is embedded // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
// isometrically in a Euclidean space. // isometrically in a Euclidean space.
// For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look // For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
// at the Riemannian Hessian". // 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 col=0; col<nDofs; col++) {
{
for (size_t subCol=0; subCol<blocksize0; subCol++) for (size_t subCol=0; subCol<blocksize; subCol++) {
{
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol]; EmbeddedTangentVector z = orthonormalFrame[col][subCol];
// P_x \partial^2 f z // P_x \partial^2 f z
for (size_t row=0; row<nDofs0; row++) for (size_t row=0; row<nDofs; row++) {
{ TangentVector semiEmbeddedProduct;
typename TargetSpace0::TangentVector semiEmbeddedProduct; embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++) for (int subRow=0; subRow<blocksize; subRow++)
localHessian[_0][_0][row*blocksize0+subRow][col*blocksize0 + subCol] = semiEmbeddedProduct[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++) // TODO: The following fails in the clang-11 C++20 CI jobs
{ //static_assert(localConfiguration.size()==1);
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol]; return assembleGradientAndHessian(localView,
localConfiguration[_0],
localGradient,
localHessian[_0][_0]);
}
else
{
int rank = Dune::MPIHelper::getCommunication().rank();
// P_x \partial^2 f z // Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
for (size_t row=0; row<nDofs0; row++) 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; typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[col][subCol];
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++) // P_x \partial^2 f z
localHessian[_0][_1][row*blocksize0+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow]; 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 col=0; col<nDofs1; col++)
{
for (size_t subCol=0; subCol<blocksize0; subCol++)
{ {
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 // P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++) { for (size_t row=0; row<nDofs0; row++)
typename TargetSpace1::TangentVector semiEmbeddedProduct; {
embeddedHessian10[row][col].mv(z,semiEmbeddedProduct); typename TargetSpace0::TangentVector semiEmbeddedProduct;
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++) for (int subRow=0; subRow<blocksize0; subRow++)
localHessian[_1][_0][row*blocksize1+subRow][col*blocksize0+subCol] = semiEmbeddedProduct[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 col=0; col<nDofs0; col++)
{
for (size_t subCol=0; subCol<blocksize1; subCol++)
{ {
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 // P_x \partial^2 f z
for (size_t row=0; row<nDofs1; row++) 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; typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[col][subCol];
embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++) // P_x \partial^2 f z
localHessian[_1][_1][row*blocksize1+subRow][col*blocksize1+subCol] = semiEmbeddedProduct[subRow]; 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 // Further correction due to non-planar configuration space
// + \mathfrak{A}_x(z,P^\orth_x \partial f) // + \mathfrak{A}_x(z,P^\orth_x \partial f)
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
// Project embedded gradient onto normal space // Project embedded gradient onto normal space
std::vector<typename TargetSpace0::EmbeddedTangentVector> projectedGradient0(nDofs0); std::vector<typename TargetSpace0::EmbeddedTangentVector> projectedGradient0(nDofs0);
for (size_t i=0; i<nDofs0; i++) for (size_t i=0; i<nDofs0; i++)
projectedGradient0[i] = localConfiguration[_0][i].projectOntoNormalSpace(localEmbeddedGradient0[i]); projectedGradient0[i] = localConfiguration[_0][i].projectOntoNormalSpace(localEmbeddedGradient0[i]);
std::vector<typename TargetSpace1::EmbeddedTangentVector> projectedGradient1(nDofs1); std::vector<typename TargetSpace1::EmbeddedTangentVector> projectedGradient1(nDofs1);
for (size_t i=0; i<nDofs1; i++) for (size_t i=0; i<nDofs1; i++)
projectedGradient1[i] = localConfiguration[_1][i].projectOntoNormalSpace(localEmbeddedGradient1[i]); projectedGradient1[i] = localConfiguration[_1][i].projectOntoNormalSpace(localEmbeddedGradient1[i]);
// The Weingarten map has only diagonal entries // The Weingarten map has only diagonal entries
for (size_t row=0; row<nDofs0; row++) for (size_t row=0; row<nDofs0; row++)
{
for (size_t subRow=0; subRow<blocksize0; subRow++)
{ {
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[row][subRow]; for (size_t subRow=0; subRow<blocksize0; subRow++)
typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration[_0][row].weingarten(z,projectedGradient0[row]); {
typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[row][subRow];
typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration[_0][row].weingarten(z,projectedGradient0[row]);
typename TargetSpace0::TangentVector tmp2; typename TargetSpace0::TangentVector tmp2;
orthonormalFrame0[row].mv(tmp1,tmp2); orthonormalFrame0[row].mv(tmp1,tmp2);
for (size_t subCol=0; subCol<blocksize0; subCol++) for (size_t subCol=0; subCol<blocksize0; subCol++)
localHessian[_0][_0][row*blocksize0+subRow][row*blocksize0+subCol] += tmp2[subCol]; localHessian[_0][_0][row*blocksize0+subRow][row*blocksize0+subCol] += tmp2[subCol];
}
} }
}
for (size_t row=0; row<nDofs1; row++) for (size_t row=0; row<nDofs1; row++)
{
for (size_t subRow=0; subRow<blocksize1; subRow++)
{ {
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[row][subRow]; for (size_t subRow=0; subRow<blocksize1; subRow++)
typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration[_1][row].weingarten(z,projectedGradient1[row]); {
typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[row][subRow];
typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration[_1][row].weingarten(z,projectedGradient1[row]);
typename TargetSpace1::TangentVector tmp2; typename TargetSpace1::TangentVector tmp2;
orthonormalFrame1[row].mv(tmp1,tmp2); orthonormalFrame1[row].mv(tmp1,tmp2);
for (size_t subCol=0; subCol<blocksize1; subCol++) for (size_t subCol=0; subCol<blocksize1; subCol++)
localHessian[_1][_1][row*blocksize1+subRow][row*blocksize1+subCol] += tmp2[subCol]; localHessian[_1][_1][row*blocksize1+subRow][row*blocksize1+subCol] += tmp2[subCol];
}
} }
} }
} }
}
} // namespace Dune::GFE
#endif #endif
dune/gfe/assemblers/localgeodesicfefdstiffness.hh
View file @ 48b9358e
...@@ -6,272 +6,278 @@ ...@@ -6,272 +6,278 @@
#include <dune/gfe/assemblers/localgeodesicfestiffness.hh> #include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
/** \brief Assembles energy gradient and Hessian with finite difference approximations
*/ namespace Dune::GFE
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;
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 typedef typename TargetSpace::template rebind<field_type>::other ATargetSpace;
constexpr static int gridDim = GridView::dimension;
public: // some other sizes
constexpr static int gridDim = GridView::dimension;
//! Dimension of a tangent space public:
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
//! Dimension of the embedding space //! Dimension of a tangent space
constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension; constexpr static int blocksize = TargetSpace::TangentVector::dimension;
LocalGeodesicFEFDStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy) //! Dimension of the embedding space
: localEnergy_(energy) constexpr static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
{}
/** \brief Compute the energy at the current configuration */ LocalGeodesicFEFDStiffness(const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* energy)
virtual RT energy (const typename Basis::LocalView& localView, : localEnergy_(energy)
const std::vector<TargetSpace>& localSolution) const override {}
{
return localEnergy_->energy(localView,localSolution);
}
RT energy (const typename Basis::LocalView& localView, /** \brief Compute the energy at the current configuration */
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override virtual RT energy (const typename Basis::LocalView& localView,
{ const std::vector<TargetSpace>& localSolution) const override
return localEnergy_->energy(localView,coefficients); {
} 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 */ /** \brief Assemble the element gradient of the energy functional
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<field_type>& gradient) const override;
/** \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 /** \brief Assemble the local tangent matrix and gradient at the current position
*/
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 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 /** \brief Assemble the local tangent matrix and gradient at the current position
*/
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_; 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
{
// /////////////////////////////////////////////////////////// template <class Basis, class TargetSpace, class field_type>
// Compute gradient by finite-difference approximation 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()); field_type eps = 1e-6;
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
}
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; localGradient.resize(localSolution.size());
std::vector<ATargetSpace> backwardSolution = localASolution;
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 for (size_t i=0; i<localSolution.size(); i++) {
const Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> B = localSolution[i].orthonormalFrame();
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]; for (int j=0; j<blocksize; j++) {
forwardCorrection *= eps;
typename ATargetSpace::EmbeddedTangentVector backwardCorrection = B[j]; typename ATargetSpace::EmbeddedTangentVector forwardCorrection = B[j];
backwardCorrection *= -eps; forwardCorrection *= eps;
forwardSolution[i] = ATargetSpace::exp(localASolution[i], forwardCorrection); typename ATargetSpace::EmbeddedTangentVector backwardCorrection = B[j];
backwardSolution[i] = ATargetSpace::exp(localASolution[i], backwardCorrection); 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 #ifdef MULTIPRECISION
localGradient[i*blocksize+j] = foo.template convert_to<double>(); localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else #else
localGradient[i*blocksize+j] = foo; localGradient[i*blocksize+j] = foo;
#endif #endif
} }
forwardSolution[i] = localASolution[i]; forwardSolution[i] = localASolution[i];
backwardSolution[i] = localASolution[i]; backwardSolution[i] = localASolution[i];
} }
} }
///////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////
// Compute gradient and Hessian together // Compute gradient and Hessian together
// To compute the Hessian we need to compute the gradient anyway, so we may // To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time. // as well return it. This saves assembly time.
///////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////
template <class Basis, class TargetSpace, class field_type> template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>:: void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradientAndHessian(const typename Basis::LocalView& localView, assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution, const std::vector<TargetSpace>& localSolution,
std::vector<field_type>& localGradient, std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{ {
// Number of degrees of freedom for this element // Number of degrees of freedom for this element
size_t nDofs = localSolution.size(); size_t nDofs = localSolution.size();
// Clear assemble data // Clear assemble data
localHessian.setSize(nDofs*blocksize, nDofs*blocksize); localHessian.setSize(nDofs*blocksize, nDofs*blocksize);
localHessian = 0; localHessian = 0;
#ifdef MULTIPRECISION #ifdef MULTIPRECISION
const field_type eps = 1e-10; const field_type eps = 1e-10;
#else #else
const field_type eps = 1e-4; const field_type eps = 1e-4;
#endif #endif
std::vector<ATargetSpace> localASolution(localSolution.size()); std::vector<ATargetSpace> localASolution(localSolution.size());
std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size()); std::vector<typename ATargetSpace::CoordinateType> aRaw(localSolution.size());
for (size_t i=0; i<localSolution.size(); i++) { for (size_t i=0; i<localSolution.size(); i++) {
typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates(); typename TargetSpace::CoordinateType raw = localSolution[i].globalCoordinates();
for (size_t j=0; j<raw.size(); j++) for (size_t j=0; j<raw.size(); j++)
aRaw[i][j] = raw[j]; aRaw[i][j] = raw[j];
localASolution[i] = aRaw[i]; localASolution[i] = aRaw[i];
} }
std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size()); std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
for (size_t i=0; i<B.size(); i++) for (size_t i=0; i<B.size(); i++)
B[i] = localSolution[i].orthonormalFrame(); B[i] = localSolution[i].orthonormalFrame();
// Precompute negative energy at the current configuration // Precompute negative energy at the current configuration
// (negative because that is how we need it as part of the 2nd-order fd formula) // (negative because that is how we need it as part of the 2nd-order fd formula)
field_type centerValue = -localEnergy_->energy(localView, localASolution); field_type centerValue = -localEnergy_->energy(localView, localASolution);
// Precompute energy infinitesimal corrections in the directions of the local basis vectors // 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> > forwardEnergy(nDofs);
std::vector<std::array<field_type,blocksize> > backwardEnergy(nDofs); std::vector<std::array<field_type,blocksize> > backwardEnergy(nDofs);
//#pragma omp parallel for schedule (dynamic) //#pragma omp parallel for schedule (dynamic)
for (size_t i=0; i<localSolution.size(); i++) { for (size_t i=0; i<localSolution.size(); i++) {
for (size_t i2=0; i2<blocksize; i2++) { for (size_t i2=0; i2<blocksize; i2++) {
typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2]; typename ATargetSpace::EmbeddedTangentVector epsXi = B[i][i2];
epsXi *= eps; epsXi *= eps;
typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi; typename ATargetSpace::EmbeddedTangentVector minusEpsXi = epsXi;
minusEpsXi *= -1; minusEpsXi *= -1;
std::vector<ATargetSpace> forwardSolution = localASolution; std::vector<ATargetSpace> forwardSolution = localASolution;
std::vector<ATargetSpace> backwardSolution = localASolution; std::vector<ATargetSpace> backwardSolution = localASolution;
forwardSolution[i] = ATargetSpace::exp(localASolution[i],epsXi); forwardSolution[i] = ATargetSpace::exp(localASolution[i],epsXi);
backwardSolution[i] = ATargetSpace::exp(localASolution[i],minusEpsXi); backwardSolution[i] = ATargetSpace::exp(localASolution[i],minusEpsXi);
forwardEnergy[i][i2] = localEnergy_->energy(localView, forwardSolution); forwardEnergy[i][i2] = localEnergy_->energy(localView, forwardSolution);
backwardEnergy[i][i2] = localEnergy_->energy(localView, backwardSolution); 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 (size_t i=0; i<localSolution.size(); i++)
for (int j=0; j<blocksize; j++) for (int j=0; j<blocksize; j++)
{ {
field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps); field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
#ifdef MULTIPRECISION #ifdef MULTIPRECISION
localGradient[i*blocksize+j] = foo.template convert_to<double>(); localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else #else
localGradient[i*blocksize+j] = foo; localGradient[i*blocksize+j] = foo;
#endif #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 #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 #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 #endif
}
} }
} }
} }
} }
}
} // namespace Dune::GFE
#endif #endif
dune/gfe/assemblers/localgeodesicfestiffness.hh
View file @ 48b9358e
...@@ -64,29 +64,30 @@ namespace Dune::GFE ...@@ -64,29 +64,30 @@ namespace Dune::GFE
// Type of the local Hessian // Type of the local Hessian
using CompositeHessian = FieldMatrix<Matrix<double>,2,2>; using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
}; };
} } // namespace Impl
}
template<class Basis, class TargetSpace> template<class Basis, class TargetSpace>
class LocalGeodesicFEStiffness class LocalGeodesicFEStiffness
: public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace> : public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
{ {
public: public:
/** \brief Assemble the local gradient and stiffness matrix at the current position /** \brief Assemble the local gradient and stiffness matrix at the current position
*/ */
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients, const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient, std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0; typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version /** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/ */
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView, virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients, const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient, std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0; typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
}; };
} // namespace Dune::GFE
#endif #endif
dune/gfe/assemblers/localintegralenergy.hh
View file @ 48b9358e
...@@ -10,36 +10,8 @@ ...@@ -10,36 +10,8 @@
#include <dune/gfe/densities/localdensity.hh> #include <dune/gfe/densities/localdensity.hh>
namespace Dune::GFE { namespace Dune::GFE
{
#if ! DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10)
namespace Impl
{
template <class Basis>
class LocalFiniteElementFactory
{
public:
static auto get(const typename Basis::LocalView& localView)
-> decltype(localView.tree().child(0).finiteElement())
{
return localView.tree().child(0).finiteElement();
}
};
/** \brief Specialize for scalar bases, here we cannot call tree().child() */
template <class GridView, int order>
class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order> >
{
public:
static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView)
-> decltype(localView.tree().finiteElement())
{
return localView.tree().finiteElement();
}
};
}
#endif
/** \brief An energy given as an integral over a density /** \brief An energy given as an integral over a density
* *
* \tparam Basis The scalar finite element basis used to construct the interpolation rule * \tparam Basis The scalar finite element basis used to construct the interpolation rule
...@@ -52,6 +24,7 @@ namespace Dune::GFE { ...@@ -52,6 +24,7 @@ namespace Dune::GFE {
{ {
using LocalView = typename Basis::LocalView; using LocalView = typename Basis::LocalView;
using GridView = typename LocalView::GridView; using GridView = typename LocalView::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using DT = typename GridView::Grid::ctype; using DT = typename GridView::Grid::ctype;
using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT; using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT;
...@@ -59,10 +32,28 @@ namespace Dune::GFE { ...@@ -59,10 +32,28 @@ namespace Dune::GFE {
public: 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(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld) LocalIntegralEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
: localDensity_(ld) 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(LocalInterpolationRule&& localGFEFunction,
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction))),
localDensity_(ld)
{} {}
private: private:
...@@ -75,34 +66,20 @@ namespace Dune::GFE { ...@@ -75,34 +66,20 @@ namespace Dune::GFE {
if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower) if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
{ {
#if DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10) // Bind GFE function to the current element
// Get an appropriate scalar local finite element, to construct the interpolation rule with const auto& element = localView.element();
// TODO: This is not a good design, for several reasons: localGFEFunction_->bind(element,localConfiguration);
// * The interpolation rule could want to have state beyond what we know here.
// It should therefore be constructed outside of the LocalIntegralEnergy class // Bind density to the element
// * I don't really see why Basis should be allowed to be scalar-valued localDensity_->bind(element);
// to begin with, but a lot of code still currently does that.
auto lfeGetter = [&localView]() // Get a suitable quadrature rule
{ const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
if constexpr (Basis::LocalView::Tree::isPower) int quadOrder = (element.type().isSimplex())
return localView.tree().child(0).finiteElement(); ? (localGFEOrder-1) * 2
else : (localGFEOrder * gridDim - 1) * 2;
return localView.tree().finiteElement();
}; const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
const auto& localFiniteElement = lfeGetter();
#else
const auto& localFiniteElement = Impl::LocalFiniteElementFactory<Basis>::get(localView);
#endif
LocalInterpolationRule localInterpolationRule(localFiniteElement,localConfiguration);
int quadOrder = (localFiniteElement.type().isSimplex())
? (localFiniteElement.localBasis().order()-1) * 2
: (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
const auto element = localView.element();
const auto& quad = QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
for (auto&& qp : quad) for (auto&& qp : quad)
{ {
...@@ -117,13 +94,13 @@ namespace Dune::GFE { ...@@ -117,13 +94,13 @@ namespace Dune::GFE {
{ {
if (localDensity_->dependsOnDerivative()) if (localDensity_->dependsOnDerivative())
{ {
auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos); auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
derivative = derivative * geometryJacobianInverse; derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative); energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
} }
else else
{ {
auto value = localInterpolationRule.evaluate(quadPos); const auto value = localGFEFunction_->evaluate(quadPos);
typename LocalInterpolationRule::DerivativeType dummyDerivative; typename LocalInterpolationRule::DerivativeType dummyDerivative;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative); energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
} }
...@@ -133,7 +110,7 @@ namespace Dune::GFE { ...@@ -133,7 +110,7 @@ namespace Dune::GFE {
if (localDensity_->dependsOnDerivative()) if (localDensity_->dependsOnDerivative())
{ {
typename TargetSpace::CoordinateType dummyValue; typename TargetSpace::CoordinateType dummyValue;
auto derivative = localInterpolationRule.evaluateDerivative(quadPos); auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
derivative = derivative * geometryJacobianInverse; derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative); energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
} }
...@@ -158,27 +135,26 @@ namespace Dune::GFE { ...@@ -158,27 +135,26 @@ namespace Dune::GFE {
{ {
RT energy = 0; RT energy = 0;
if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold && Basis::LocalView::Tree::isComposite) if constexpr (Impl::LocalEnergyTypes<TargetSpace>::isProductManifold
&& Basis::LocalView::Tree::isComposite
&& gridDim==GridView::dimensionworld) // TODO: Implement the case gridDim!=dimworld
{ {
static_assert(TargetSpace::size() == 2, static_assert(TargetSpace::size() == 2,
"LocalIntegralEnergy only implemented for product spaces with two factors!"); "LocalIntegralEnergy only implemented for product spaces with two factors!");
const auto& element = localView.element();
using namespace Indices; using namespace Indices;
// composite Basis: grab the finite element of the first child // Bind GFE functions to the current element
const auto& localFiniteElement0 = localView.tree().child(_0,0).finiteElement(); const auto& element = localView.element();
const auto& localFiniteElement1 = localView.tree().child(_1,0).finiteElement(); std::get<0>(*localGFEFunction_).bind(element, coefficients[_0]);
std::get<1>(*localGFEFunction_).bind(element, coefficients[_1]);
using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
LocalGFEFunctionType0 localGFEFunction0(localFiniteElement0, coefficients[_0]); // Bind density to the element
LocalGFEFunctionType1 localGFEFunction1(localFiniteElement1, coefficients[_1]); localDensity_->bind(element);
int quadOrder = (element.type().isSimplex()) ? localFiniteElement0.localBasis().order() // Get a suitable quadrature rule
: localFiniteElement0.localBasis().order() * gridDim; 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); const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
...@@ -199,19 +175,19 @@ namespace Dune::GFE { ...@@ -199,19 +175,19 @@ namespace Dune::GFE {
// A value is needed either directly, or for computing the derivative. // A value is needed either directly, or for computing the derivative.
TargetSpace value; TargetSpace value;
if (localDensity_->dependsOnValue(0) || localDensity_->dependsOnDerivative(0)) 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)) 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 // Compute the derivatives of the interpolation function factors
typename LocalGFEFunctionType0::DerivativeType derivative0; typename std::tuple_element_t<0, LocalInterpolationRule>::DerivativeType derivative0;
typename LocalGFEFunctionType1::DerivativeType derivative1; typename std::tuple_element_t<1, LocalInterpolationRule>::DerivativeType derivative1;
if (localDensity_->dependsOnDerivative(0)) 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)) 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 // Copy the two derivatives into a joint matrix object
// TODO: I am not sure about this. May the densities should get the // TODO: I am not sure about this. May the densities should get the
...@@ -230,13 +206,18 @@ namespace Dune::GFE { ...@@ -230,13 +206,18 @@ namespace Dune::GFE {
} }
} }
else else
DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis"); DUNE_THROW(Dune::NotImplemented, "Non-product manifold or non-composite basis or gridDim!=dimworld");
return energy; return energy;
} }
protected: 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 } // namespace Dune::GFE
......