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)
Expand all 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
d13f34469065230dd86e63733e4dc8d1f775acfb
# Cleanup after having moved everything into the Dune::GFE namespace
c2eb8dacb5885363c30c2b7373c479ed9c66ac12
.gitlab-ci.yml
View file @ 48b9358e
......@@ -7,9 +7,10 @@ before_script: &before
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-vtk.git
- duneci-install-module https://gitlab.dune-project.org/fufem/dune-fufem.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/iwr/dune-gmsh4
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-gmsh4
# Tests with the 2.9 release. That's the one in Debian bookworm
#--------------------------------------------------------------------
dune:2.9 gcc:
variables:
DUNECI_BRANCH: releases/2.9
......@@ -18,6 +19,76 @@ dune:2.9 gcc:
- *before
script: duneci-standard-test
dune:2.9 dune-elasticity gcc:
variables:
DUNECI_BRANCH: releases/2.9
image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Tests with the 2.10 release. That's the one in Debian trixie
#--------------------------------------------------------------------
dune:2.10 gcc-10:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
script: duneci-standard-test
dune:2.10 dune-elasticity gcc-10:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:2.10 dune-curvedgrid dune-foamgrid gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Full tests again, but this time with clang.
dune:2.10 dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-clang-11-20
variables:
DUNECI_BRANCH: releases/2.10
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Tests with the git master branch
# Some, but not all optional dependencies are made available.
#------------------------------------------------------------------
dune:git gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
......@@ -30,24 +101,44 @@ dune:git clang-11 C++20:
- *before
script: duneci-standard-test
dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid gcc-10 C++20:
dune:git dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:git dune-curvedgrid dune-foamgrid gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
script: duneci-standard-test
# Tests with all optional dependencies
#-------------------------------------------
dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity gcc-10 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
dune:git dune-parmg dune-curvedgeometry dune-curvedgrid dune-foamgrid clang-11 C++20:
dune:git dune-parmg dune-curvedgrid dune-foamgrid dune-elasticity clang-11 C++20:
image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-11-20
before_script:
- *before
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/paraphase/dune-parmg.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git
- duneci-install-module https://gitlab-ci-token:${CI_JOB_TOKEN}@gitlab.mn.tu-dresden.de/iwr/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgeometry.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-curvedgrid.git
- duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git
- duneci-install-module https://gitlab.mn.tu-dresden.de/ag-sander/dune/dune-elasticity.git
script: duneci-standard-test
# Check for spelling mistakes in text
......
CHANGELOG.md
View file @ 48b9358e
# Master
- The `SimoFoxEnergy` class (formerly known as `SimoFoxEnergyLocalStiffness`)
does not accept a surface load density anymore. Use the
`NeumannEnergy` class for assembling surface loads.
- The class `LocalIntegralEnergy` now does not construct the local interpolation
rule itself anymore. Rather, it expects to be given one in its constructor.
This will allow interpolation rules with state.
- The same holds for `CosseratRodEnergy`. Unfortunately, that move
implies that `CosseratRodEnergy` gains an extra template argument
which specifies the GFE function used to implement the reference
configuration.
- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
now take a `dune-functions` basis as their first arguments, instead
of a `LocalFiniteElement`. In addition, they need to be given an object
of this type in the constructor. This is the basis to be used for
interpolation weights (for `LocalGeodesicFEFunction`) and for
interpolating in the embedding space (for `LocalProjectedFEFunction`).
- The entire code is now in the `Dune::GFE` namespace.
- Remove the files `localgfetestfunctionbasis.hh` and `localtangentfefunction.hh`.
They were never used for anything at all, I and currently wouldn't know
what to use them for, either.
- All implementations of functions have been move to a `functions` subdirectory.
- The `LocalDensity` class now has a template parameter `ElementOrIntersection`,
which replaces the parameter `Position`. As the name says, this parameter
has to be a grid element, or a grid intersection. As it turned out,
just depending on the position type of the integration domain was not enough.
- The `LocalDensity` class now has a `bind` method, which binds the density
to a given element or intersection. This is necessary, for example, to in turn
bind coefficient functions that the density may own.
# Release 2.10 (2024-12-29)
- The module `dune-elasticity` has been downgraded from a required
dependency to an optional dependency. It is currently only needed
by the `film-on-substrate.cc` program.
- `cosserat-rod.cc` (formerly `rod3d.cc`) now reads the reference configuration
and the initial iterate from the Python file.
- The file `rod3d.cc` has been renamed to `cosserat-rod.cc`,
to better reflect what it does.
- In the `BulkCosseratDensity` class: Make the type of curvature
tensor controllable at run-time (previously, ugly preprocessor switches
where used). Among its parameters `BulkCosseratDensity` now expects
a string `curvatureType`, which has to take one of the values `norm`,
`curl`, or `wryness`.
- Building `dune-gfe` now requires Dune 2.9 or newer. Older configurations
have not gotten CI-tested for quite a while, anyway.
......
CMakeLists.txt
View file @ 48b9358e
......@@ -17,6 +17,13 @@ include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
# Create a module library
if(dune-common_VERSION VERSION_GREATER_EQUAL 2.9.0)
dune_add_library(dunegfe INTERFACE EXPORT_NAME GFE NAMESPACE Dune::
LINK_LIBRARIES ${DUNE_LIBS})
endif()
dune_enable_all_packages()
add_subdirectory("src")
add_subdirectory("dune")
......
dune.module
View file @ 48b9358e
......@@ -4,8 +4,8 @@
#Name of the module
Module: dune-gfe
Version: svn
Version: 2.11-git
Maintainer: oliver.sander@tu-dresden.de
#depending on
Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-elasticity (>=2.9) dune-gmsh4
Suggests: dune-foamgrid dune-parmg dune-vtk dune-curvedgrid
Depends: dune-grid(>=2.9) dune-uggrid dune-istl dune-localfunctions dune-functions (>=2.9) dune-solvers dune-fufem dune-gmsh4 dune-vtk
Suggests: dune-foamgrid dune-parmg dune-curvedgrid dune-elasticity
dune/gfe/CMakeLists.txt
View file @ 48b9358e
......@@ -4,20 +4,12 @@ install(FILES
cosseratstrain.hh
cosseratvtkreader.hh
cosseratvtkwriter.hh
embeddedglobalgfefunction.hh
filereader.hh
geodesicdifference.hh
geodesicfefunctionadaptor.hh
gfedifferencenormsquared.hh
globalgfefunction.hh
gramschmidtsolver.hh
interpolationderivatives.hh
linearalgebra.hh
localgeodesicfefunction.hh
localgfetestfunctionbasis.hh
localprojectedfefunction.hh
localquickanddirtyfefunction.hh
localtangentfefunction.hh
maxnormtrustregion.hh
mixedriemanniantrsolver.cc
mixedriemanniantrsolver.hh
......@@ -42,3 +34,4 @@ install(FILES
add_subdirectory("spaces")
add_subdirectory("assemblers")
add_subdirectory("densities")
add_subdirectory("functions")
dune/gfe/assemblers/CMakeLists.txt
View file @ 48b9358e
#install headers
install(FILES
cosseratenergystiffness.hh
cosseratrodenergy.hh
geodesicfeassembler.hh
geodesicfeassemblerwrapper.hh
......
dune/gfe/assemblers/cosseratenergystiffness.hh deleted 100644 → 0
View file @ 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 @@
#define DUNE_GFE_COSSERATRODENERGY_HH
#include <array>
#include <memory>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
......@@ -18,9 +19,10 @@
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
namespace Dune::GFE {
namespace Dune::GFE
{
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
class CosseratRodEnergy
: public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
{
......@@ -42,6 +44,14 @@ namespace Dune::GFE {
public:
/** \brief The current configuration whose energy is being computed
*/
const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
/** \brief The GFE function that implements the stress-free reference configuration
*/
const std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration_;
/** \brief The stress-free configuration
The number type cannot be RT, because RT become `adouble` when
......@@ -67,9 +77,25 @@ namespace Dune::GFE {
GridView gridView_;
//! Constructor
CosseratRodEnergy(const GridView& gridView,
CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
: localGFEFunction_(localGFEFunction)
, localReferenceConfiguration_(localReferenceConfiguration)
, K_(K),
A_(A),
gridView_(gridView)
{}
//! Constructor
CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
ReferenceInterpolationRule&& localReferenceConfiguration,
const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
: K_(K),
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
, localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
, K_(K),
A_(A),
gridView_(gridView)
{}
......@@ -80,9 +106,39 @@ namespace Dune::GFE {
\param E Young's modulus
\param nu Poisson number
*/
CosseratRodEnergy(const GridView& gridView,
CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
const GridView& gridView,
double A, double J1, double J2, double E, double nu)
: gridView_(gridView)
: localGFEFunction_(localGFEFunction)
, localReferenceConfiguration_(localReferenceConfiguration)
, gridView_(gridView)
{
// shear modulus
double G = E/(2+2*nu);
K_[0] = E * J1;
K_[1] = E * J2;
K_[2] = G * (J1 + J2);
A_[0] = G * A;
A_[1] = G * A;
A_[2] = E * A;
}
/** \brief Constructor setting shape constants and material parameters
\param A The rod section area
\param J1, J2 The geometric moments (Flächenträgheitsmomente)
\param E Young's modulus
\param nu Poisson number
*/
CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
ReferenceInterpolationRule&& localReferenceConfiguration,
const GridView& gridView,
double A, double J1, double J2, double E, double nu)
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
, localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
, gridView_(gridView)
{
// shear modulus
double G = E/(2+2*nu);
......@@ -172,21 +228,18 @@ namespace Dune::GFE {
};
template<class Basis, class LocalInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localCoefficients) const
{
const auto& localFiniteElement = localView.tree().finiteElement();
LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
const auto& element = localView.element();
localGFEFunction_->bind(element, localCoefficients);
RT energy = 0;
std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::other;
InactiveLocalInterpolationRule localReferenceConfiguration(localFiniteElement, localReferenceCoefficients);
localReferenceConfiguration_->bind(element, localReferenceCoefficients);
// ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate
......@@ -205,10 +258,10 @@ namespace Dune::GFE {
double weight = shearingQuad[pt].weight() * integrationElement;
auto strain = getStrain(localConfiguration, element, quadPos);
auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
for (int i=0; i<3; i++)
energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
......@@ -225,10 +278,10 @@ namespace Dune::GFE {
double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
auto strain = getStrain(localConfiguration, element, quadPos);
auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
// Part II: the bending and twisting energy
for (int i=0; i<3; i++)
......@@ -240,9 +293,9 @@ namespace Dune::GFE {
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class ReboundLocalInterpolationRule>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const ReboundLocalInterpolationRule& localInterpolation,
const Entity& element,
const FieldVector<double,1>& pos) const
......@@ -285,9 +338,9 @@ namespace Dune::GFE {
return strain;
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class Number>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element,
const FieldVector<double, 1>& pos) const
......@@ -308,8 +361,8 @@ namespace Dune::GFE {
return stress;
}
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const
{
......@@ -363,8 +416,8 @@ namespace Dune::GFE {
}
}
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const
{
......@@ -386,9 +439,9 @@ namespace Dune::GFE {
}
}
template<class Basis, class LocalInterpolationRule, class RT>
template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class PatchGridView>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const
{
......
dune/gfe/assemblers/discretekirchhoffbendingenergy.hh 0 → 100644
View file @ 48b9358e
#ifndef DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
#define DUNE_GFE_ASSEMBLERS_DISCRETEKIRCHHOFFBENDINGENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/bendingisometryhelper.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
/**
* \brief Assemble the discrete Kirchhoff bending energy for a single element.
*
* The energy is essentially the harmonic energy of a
* discrete Jacobian operator which is represented as a
* linear combination in a P2 finite element space.
* The gradient of that linear combination serves as an
* discrete Hessian in the energy.
*/
namespace Dune::GFE
{
template<class Basis, class DiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
class DiscreteKirchhoffBendingEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
typedef typename DiscreteKirchhoffFunction::DerivativeType DerivativeType;
// Grid dimension
constexpr static int gridDim = GridView::dimension;
// P2-LocalFiniteElement used to represent the discrete Jacobian on the current element.
typedef typename Dune::LagrangeSimplexLocalFiniteElement<DT, double, gridDim, 2> P2LocalFiniteElement;
public:
DiscreteKirchhoffBendingEnergy(DiscreteKirchhoffFunction &localFunction)
: localFunction_(localFunction)
{}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(NotImplemented, "!");
}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localConfiguration) const override
{
const auto element = localView.element();
int quadOrder = 3;
const auto &quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
/**
bind the underlying Hermite basis to the current element.
Note: The DiscreteKirchhoffFunction object stores a vector of global
and local basis coefficients that are accessed for evaluation
methods (linear combinations of hermite basis function and local coefficients).
Therefore the local configuration is passed to the bind method as well since
these coefficient vectors need to be updated with the new local configuration
previous to evaluation.
*/
localFunction_.bind(element,localConfiguration);
RT bendingEnergy = 0;
for (auto&& quadPoint : quad)
{
auto discreteHessian= localFunction_.evaluateDiscreteHessian(quadPoint.position());
auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
bendingEnergy += weight * discreteHessian.frobenius_norm2();
}
return 0.5 * bendingEnergy;
}
private:
DiscreteKirchhoffFunction& localFunction_;
P2LocalFiniteElement lagrangeLFE_;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/forceenergy.hh 0 → 100644
View file @ 48b9358e
#ifndef DUNE_GFE_FORCEENERGY_HH
#define DUNE_GFE_FORCEENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/bendingisometryhelper.hh>
/**
* \brief Assemble the energy of a force term (force times deformation) for a single element.
*/
namespace Dune::GFE
{
template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
class ForceEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
// Grid dimension
constexpr static int gridDim = GridView::dimension;
public:
ForceEnergy(LocalDiscreteKirchhoffFunction &localFunction,
LocalForce &localForce)
: localFunction_(localFunction),
localForce_(localForce)
{}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(NotImplemented, "!");
}
/** \brief Assemble the energy for a single element */
virtual RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localConfiguration) const override
{
RT energy = 0;
int quadOrder = 3;
const auto &quad = QuadratureRules<double, gridDim>::rule(GeometryTypes::simplex(gridDim), quadOrder);
const auto element = localView.element();
/** bind the local force to the current element */
localForce_.bind(element);
/**
bind the underlying Hermite basis to the current element.
Note: The DiscreteKirchhoffFunction object stores a vector of global
and local basis coefficients that are accessed for evaluation
methods (linear comibnations of hermite basis function and local coefficients).
Therefore the local configuration is passed to the bind method as well since
these coefficient vectors need to be updated with the new local configuration
previous to evaluation.
*/
localFunction_.bind(element,localConfiguration);
/**
* @brief Compute contribution of the force Term
*
* Integrate the scalar product of the force
* with the local deformation function 'localFunction_'
* over the current element.
*/
RT forceTermEnergy = 0;
for (auto&& quadPoint : quad)
{
auto deformationValue = localFunction_.evaluate(quadPoint.position());
auto forceValue = localForce_(quadPoint.position());
auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
}
return forceTermEnergy;
}
private:
LocalDiscreteKirchhoffFunction& localFunction_;
LocalForce& localForce_;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/geodesicfeassembler.hh
View file @ 48b9358e
......@@ -10,269 +10,273 @@
#include <dune/solvers/common/wrapownshare.hh>
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
using field_type = typename TargetSpace::field_type;
typedef typename Basis::GridView GridView;
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
namespace Dune::GFE
{
//! Dimension of the grid.
constexpr static int gridDim = GridView::dimension;
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
using field_type = typename TargetSpace::field_type;
typedef typename Basis::GridView GridView;
using LocalStiffness = LocalGeodesicFEStiffness<Basis, TargetSpace>;
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
//! Dimension of the grid.
constexpr static int gridDim = GridView::dimension;
//! Dimension of a tangent space
constexpr static int blocksize = TargetSpace::TangentVector::dimension;
protected:
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
//! The global basis
const Basis basis_;
//! The local stiffness operator
std::shared_ptr<LocalStiffness> localStiffness_;
protected:
public:
//! The global basis
const Basis basis_;
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const Basis& basis)
: basis_(basis)
{}
//! The local stiffness operator
std::shared_ptr<LocalStiffness> localStiffness_;
/** \brief Constructor for a given grid */
template <class LocalStiffnessT>
GeodesicFEAssembler(const Basis& basis,
LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
public:
/** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
template <class LocalStiffnessT>
void setLocalStiffness(LocalStiffnessT&& localStiffness) {
localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
}
/** \brief Get the local stiffness operator. */
const LocalStiffness& getLocalStiffness() const {
return *localStiffness_;
}
/** \brief Constructor for a given grid */
GeodesicFEAssembler(const Basis& basis)
: basis_(basis)
{}
/** \brief Get the local stiffness operator. */
LocalStiffness& getLocalStiffness() {
return *localStiffness_;
}
/** \brief Constructor for a given grid */
template <class LocalStiffnessT>
GeodesicFEAssembler(const Basis& basis,
LocalStiffnessT&& localStiffness)
: basis_(basis),
localStiffness_(Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness)))
{}
/** \brief Get the basis. */
const Basis& getBasis() const {
return basis_;
}
/** \brief Set the local stiffness assembler. This can be a temporary, l-value or shared pointer. */
template <class LocalStiffnessT>
void setLocalStiffness(LocalStiffnessT&& localStiffness) {
localStiffness_ = Dune::Solvers::wrap_own_share<LocalStiffness>(std::forward<LocalStiffnessT>(localStiffness));
}
/** \brief Assemble the tangent stiffness matrix and the functional gradient together
*
* This is more efficient than computing them separately, because you need the gradient
* anyway to compute the Riemannian Hessian.
*/
virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const;
/** \brief Get the local stiffness operator. */
const LocalStiffness& getLocalStiffness() const {
return *localStiffness_;
}
/** \brief Assemble the gradient */
virtual void assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Get the local stiffness operator. */
LocalStiffness& getLocalStiffness() {
return *localStiffness_;
}
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
/** \brief Get the basis. */
const Basis& getBasis() const {
return basis_;
}
//protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
/** \brief Assemble the tangent stiffness matrix and the functional gradient together
*
* This is more efficient than computing them separately, because you need the gradient
* anyway to compute the Riemannian Hessian.
*/
virtual void assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const;
}; // end class
/** \brief Assemble the gradient */
virtual void assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
//protected:
void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
auto n = basis_.size();
}; // end class
nb.resize(n, n);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
// Bind the local FE basis view to the current element
localView.bind(element);
auto n = basis_.size();
const auto& lfe = localView.tree().finiteElement();
nb.resize(n, n);
for (size_t i=0; i<lfe.size(); i++) {
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t j=0; j<lfe.size(); j++) {
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
const auto& lfe = localView.tree().finiteElement();
nb.add(iIdx, jIdx);
for (size_t i=0; i<lfe.size(); i++) {
}
for (size_t j=0; j<lfe.size(); j++) {
}
auto iIdx = localView.index(i);
auto jIdx = localView.index(j);
}
nb.add(iIdx, jIdx);
}
}
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
{
if (computeOccupationPattern) {
}
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
neighborsPerVertex.exportIdx(hessian);
}
}
hessian = 0;
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > &gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
{
if (computeOccupationPattern) {
gradient.resize(sol.size());
gradient = 0;
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
neighborsPerVertex.exportIdx(hessian);
// A view on the FE basis on a single element
auto localView = basis_.localView();
}
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
hessian = 0;
const int numOfBaseFct = localView.tree().size();
gradient.resize(sol.size());
gradient = 0;
// Extract local solution
std::vector<TargetSpace> localSolution(numOfBaseFct);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[localView.index(i)];
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
std::vector<double> localGradient(numOfBaseFct*blocksize);
Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
const int numOfBaseFct = localView.tree().size();
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
// Extract local solution
std::vector<TargetSpace> localSolution(numOfBaseFct);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[localView.index(i)];
auto row = localView.index(i);
std::vector<double> localGradient(numOfBaseFct*blocksize);
Dune::Matrix<double> localHessian(numOfBaseFct*blocksize, numOfBaseFct*blocksize);
for (int j=0; j<numOfBaseFct; j++ ) {
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
auto col = localView.index(j);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
for (std::size_t ii=0; ii<blocksize; ii++)
for (std::size_t jj=0; jj<blocksize; jj++)
hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
auto row = localView.index(i);
}
}
for (int j=0; j<numOfBaseFct; j++ ) {
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
for (int j=0; j<blocksize; j++)
gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
auto col = localView.index(j);
}
for (std::size_t ii=0; ii<blocksize; ii++)
for (std::size_t jj=0; jj<blocksize; jj++)
hessian[row][col][ii][jj] += localHessian[i*blocksize+ii][j*blocksize+jj];
}
}
}
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
for (int j=0; j<blocksize; j++)
gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
grad.resize(sol.size());
grad = 0;
}
// A view on the FE basis on a single element
auto localView = basis_.localView();
}
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
localView.bind(element);
if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
// The number of degrees of freedom of the current element
const auto nDofs = localView.tree().size();
grad.resize(sol.size());
grad = 0;
// Extract local solution
std::vector<TargetSpace> localSolution(nDofs);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)];
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
// Assemble local gradient
std::vector<double> localGradient(nDofs*blocksize);
// The number of degrees of freedom of the current element
const auto nDofs = localView.tree().size();
localStiffness_->assembleGradient(localView, localSolution, localGradient);
// Extract local solution
std::vector<TargetSpace> localSolution(nDofs);
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<blocksize; j++)
grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)];
}
// Assemble local gradient
std::vector<double> localGradient(nDofs*blocksize);
localStiffness_->assembleGradient(localView, localSolution, localGradient);
template <class Basis, class TargetSpace>
double GeodesicFEAssembler<Basis, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
double energy = 0;
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
for (size_t j=0; j<blocksize; j++)
grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
if (sol.size() != basis_.size())
DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
}
// A view on the FE basis on a single element
auto localView = basis_.localView();
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
template <class Basis, class TargetSpace>
double GeodesicFEAssembler<Basis, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
localView.bind(element);
double energy = 0;
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
if (sol.size() != basis_.size())
DUNE_THROW(Dune::Exception, "Coefficient vector doesn't match the function space basis!");
std::vector<TargetSpace> localSolution(nDofs);
// A view on the FE basis on a single element
auto localView = basis_.localView();
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)[0]];
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
localView.bind(element);
energy += localStiffness_->energy(localView, localSolution);
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
}
std::vector<TargetSpace> localSolution(nDofs);
return energy;
for (size_t i=0; i<nDofs; i++)
localSolution[i] = sol[localView.index(i)[0]];
}
energy += localStiffness_->energy(localView, localSolution);
}
return energy;
}
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
View file @ 48b9358e
......@@ -5,7 +5,8 @@
#include <dune/common/tuplevector.hh>
#include <dune/gfe/spaces/productmanifold.hh>
namespace Dune::GFE {
namespace Dune::GFE
{
/** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
......@@ -13,9 +14,8 @@ namespace Dune::GFE {
*/
template <class Basis, class ScalarBasis, class TargetSpace>
class
GeodesicFEAssemblerWrapper {
class GeodesicFEAssemblerWrapper
{
using MixedSpace0 = std::tuple_element_t<0,TargetSpace>;
using MixedSpace1 = std::tuple_element_t<1,TargetSpace>;
......@@ -72,7 +72,8 @@ namespace Dune::GFE {
auto splitVector(const std::vector<TargetSpace>& sol) const;
std::unique_ptr<MatrixType> hessianMixed_;
}; // end class
} //end namespace
} // namespace Dune::GFE
template <class Basis, class ScalarBasis, class TargetSpace>
......
dune/gfe/assemblers/l2distancesquaredenergy.hh
View file @ 48b9358e
......@@ -5,86 +5,94 @@
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/globalgfefunction.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/functions/globalgfefunction.hh>
#include <dune/gfe/functions/localgeodesicfefunction.hh>
template<class Basis, class TargetSpace>
class L2DistanceSquaredEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
namespace Dune::GFE
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
// some other sizes
constexpr static int gridDim = GridView::dimension;
template<class Basis, class TargetSpace>
class L2DistanceSquaredEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
using LocalInterpolationRule = LocalGeodesicFEFunction<gridDim, DT, typename Basis::LocalView::Tree::FiniteElement, typename TargetSpace::template rebind<double>::other>;
// some other sizes
constexpr static int gridDim = GridView::dimension;
public:
using LocalInterpolationRule = LocalGeodesicFEFunction<Basis, typename TargetSpace::template rebind<double>::other>;
// This is the function that we are computing the L2-distance to
std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
public:
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
// This is the function that we are computing the L2-distance to
std::shared_ptr<Dune::GFE::GlobalGFEFunction<Basis, LocalInterpolationRule, typename TargetSpace::template rebind<double>::other > > origin_;
{
RT energy = 0;
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const override
const auto& localFiniteElement = localView.tree().finiteElement();
typedef LocalGeodesicFEFunction<gridDim, double, decltype(localFiniteElement), TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localFiniteElement,localSolution);
{
RT energy = 0;
const auto element = localView.element();
auto localOrigin = localFunction(*origin_);
localOrigin.bind(element);
typedef LocalGeodesicFEFunction<Basis, TargetSpace> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(localView.globalBasis());
// Just guessing an appropriate quadrature order
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
const auto element = localView.element();
localGeodesicFEFunction.bind(element,localSolution);
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
auto localOrigin = localFunction(*origin_);
localOrigin.bind(element);
for (size_t pt=0; pt<quad.size(); pt++)
{
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
// Just guessing an appropriate quadrature order
const auto& localFiniteElement = localView.tree().finiteElement();
auto quadOrder = localFiniteElement.localBasis().order() * 2 * gridDim;
const auto& quad = Dune::QuadratureRules<double, gridDim>::rule(localFiniteElement.type(), quadOrder);
const double integrationElement = element.geometry().integrationElement(quadPos);
for (size_t pt=0; pt<quad.size(); pt++)
{
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
auto weight = quad[pt].weight() * integrationElement;
const double integrationElement = element.geometry().integrationElement(quadPos);
// The function value
auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originValue = localOrigin(quadPos);
auto weight = quad[pt].weight() * integrationElement;
// The derivative of the 'origin' function
// First: as function defined on the reference element
auto originReferenceDerivative = derivative(localOrigin)(quadPos);
// The function value
auto value = localGeodesicFEFunction.evaluate(quadPos);
auto originValue = localOrigin(quadPos);
// The derivative of the function defined on the actual element
auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
// The derivative of the 'origin' function
// First: as function defined on the reference element
auto originReferenceDerivative = derivative(localOrigin)(quadPos);
double weightFactor = originDerivative.frobenius_norm();
// un-comment the following line to switch off the weight factor
//weightFactor = 1.0;
// The derivative of the function defined on the actual element
auto originDerivative = originReferenceDerivative * element.geometry().jacobianInverse(quadPos);
// Add the local energy density
energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
double weightFactor = originDerivative.frobenius_norm();
// un-comment the following line to switch off the weight factor
//weightFactor = 1.0;
// Add the local energy density
energy += weight * weightFactor * TargetSpace::distanceSquared(originValue, value);
}
return energy;
}
return energy;
}
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
virtual RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
}
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localenergy.hh
View file @ 48b9358e
......@@ -7,9 +7,9 @@
#include <dune/gfe/spaces/productmanifold.hh>
namespace Dune {
namespace GFE:: Impl
namespace Dune::GFE
{
namespace Impl
{
/** \brief A class exporting container types for coefficient sets
*
......@@ -34,43 +34,40 @@ namespace Dune {
using Coefficients = std::vector<ProductManifold<Factors...> >;
using CompositeCoefficients = TupleVector<std::vector<Factors>... >;
};
}
namespace GFE {
} // namespace Impl
/** \brief Base class for energies defined by integrating over one grid element */
template<class Basis, class TargetSpace>
class LocalEnergy
{
public:
using RT = typename TargetSpace::ctype;
/** \brief Compute the energy
*
* \param localView Local view specifying the current element and the FE space there
* \param coefficients The coefficients of a FE function on the current element
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** \brief Base class for energies defined by integrating over one grid element */
template<class Basis, class TargetSpace>
class LocalEnergy
{
public:
using RT = typename TargetSpace::ctype;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
/** \brief Compute the energy
*
* \param localView Local view specifying the current element and the FE space there
* \param coefficients The coefficients of a FE function on the current element
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients) const = 0;
/** Empty virtual default destructor
*
* To allow proper destruction of derived classes through a base class pointer
*/
virtual ~LocalEnergy() = default;
/** \brief ProductManifolds: Compute the energy from coefficients in separate containers
* for each factor
*/
virtual RT
energy (const typename Basis::LocalView& localView,
const typename Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const = 0;
};
/** Empty virtual default destructor
*
* To allow proper destruction of derived classes through a base class pointer
*/
virtual ~LocalEnergy() = default;
} // namespace GFE
};
} // namespace Dune
} // namespace Dune::GFE
#endif // DUNE_GFE_LOCALENERGY_HH
dune/gfe/assemblers/localfirstordermodel.hh
View file @ 48b9358e
......@@ -20,6 +20,6 @@ namespace Dune::GFE
};
} // namespace Dune
} // namespace Dune::GFE
#endif // DUNE_GFE_LOCALFIRSTORDERMODEL_HH
dune/gfe/assemblers/localgeodesicfestiffness.hh
View file @ 48b9358e
......@@ -64,29 +64,30 @@ namespace Dune::GFE
// Type of the local Hessian
using CompositeHessian = FieldMatrix<Matrix<double>,2,2>;
};
}
}
} // namespace Impl
template<class Basis, class TargetSpace>
class LocalGeodesicFEStiffness
: public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
{
public:
/** \brief Assemble the local gradient and stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
template<class Basis, class TargetSpace>
class LocalGeodesicFEStiffness
: public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
{
public:
/** \brief Assemble the local gradient and stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
std::vector<double>& localGradient,
typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
std::vector<double>& localGradient,
typename GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
} // namespace Dune::GFE
#endif
dune/gfe/assemblers/localintegralenergy.hh
View file @ 48b9358e
......@@ -10,36 +10,8 @@
#include <dune/gfe/densities/localdensity.hh>
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
namespace Dune::GFE
{
/** \brief An energy given as an integral over a density
*
* \tparam Basis The scalar finite element basis used to construct the interpolation rule
......@@ -52,6 +24,7 @@ namespace Dune::GFE {
{
using LocalView = typename Basis::LocalView;
using GridView = typename LocalView::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using DT = typename GridView::Grid::ctype;
using RT = typename GFE::LocalEnergy<Basis,TargetSpace>::RT;
......@@ -59,10 +32,28 @@ namespace Dune::GFE {
public:
/** \brief Constructor with a Dune::GFE::LocalDensity
/** \brief Constructor from a GFE function as a shared pointer
*
* \param localGFEFunction The geometric finite element function that
* the density will be evaluated on
* \param ld The density that will be integrated
*/
LocalIntegralEnergy(const std::shared_ptr<GFE::LocalDensity<FieldVector<DT,gridDim>,TargetSpace> >& ld)
: localDensity_(ld)
LocalIntegralEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
: localGFEFunction_(localGFEFunction),
localDensity_(ld)
{}
/** \brief Constructor from a GFE function r-value reference
*
* \param localGFEFunction The geometric finite element function that
* the density will be evaluated on
* \param ld The density that will be integrated
*/
LocalIntegralEnergy(LocalInterpolationRule&& localGFEFunction,
const std::shared_ptr<GFE::LocalDensity<Element,TargetSpace> >& ld)
: localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction))),
localDensity_(ld)
{}
private:
......@@ -75,34 +66,20 @@ namespace Dune::GFE {
if constexpr (Basis::LocalView::Tree::isLeaf || Basis::LocalView::Tree::isPower)
{
#if DUNE_VERSION_GTE(DUNE_LOCALFUNCTIONS, 2, 10)
// Get an appropriate scalar local finite element, to construct the interpolation rule with
// TODO: This is not a good design, for several reasons:
// * The interpolation rule could want to have state beyond what we know here.
// It should therefore be constructed outside of the LocalIntegralEnergy class
// * I don't really see why Basis should be allowed to be scalar-valued
// to begin with, but a lot of code still currently does that.
auto lfeGetter = [&localView]()
{
if constexpr (Basis::LocalView::Tree::isPower)
return localView.tree().child(0).finiteElement();
else
return localView.tree().finiteElement();
};
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);
// Bind GFE function to the current element
const auto& element = localView.element();
localGFEFunction_->bind(element,localConfiguration);
// Bind density to the element
localDensity_->bind(element);
// Get a suitable quadrature rule
const auto localGFEOrder = localGFEFunction_->localFiniteElement().localBasis().order();
int quadOrder = (element.type().isSimplex())
? (localGFEOrder-1) * 2
: (localGFEOrder * gridDim - 1) * 2;
const auto& quad = QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
for (auto&& qp : quad)
{
......@@ -117,13 +94,13 @@ namespace Dune::GFE {
{
if (localDensity_->dependsOnDerivative())
{
auto [value, derivative] = localInterpolationRule.evaluateValueAndDerivative(quadPos);
auto [value, derivative] = localGFEFunction_->evaluateValueAndDerivative(quadPos);
derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),derivative);
}
else
{
auto value = localInterpolationRule.evaluate(quadPos);
const auto value = localGFEFunction_->evaluate(quadPos);
typename LocalInterpolationRule::DerivativeType dummyDerivative;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,value.globalCoordinates(),dummyDerivative);
}
......@@ -133,7 +110,7 @@ namespace Dune::GFE {
if (localDensity_->dependsOnDerivative())
{
typename TargetSpace::CoordinateType dummyValue;
auto derivative = localInterpolationRule.evaluateDerivative(quadPos);
auto derivative = localGFEFunction_->evaluateDerivative(quadPos);
derivative = derivative * geometryJacobianInverse;
energy += qp.weight() * integrationElement * (*localDensity_)(quadPos,dummyValue,derivative);
}
......@@ -158,27 +135,26 @@ namespace Dune::GFE {
{
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,
"LocalIntegralEnergy only implemented for product spaces with two factors!");
const auto& element = localView.element();
using namespace Indices;
// composite Basis: grab the finite element of the first child
const auto& localFiniteElement0 = localView.tree().child(_0,0).finiteElement();
const auto& localFiniteElement1 = localView.tree().child(_1,0).finiteElement();
using LocalGFEFunctionType0 = typename std::tuple_element<0, LocalInterpolationRule>::type;
using LocalGFEFunctionType1 = typename std::tuple_element<1, LocalInterpolationRule>::type;
// Bind GFE functions to the current element
const auto& element = localView.element();
std::get<0>(*localGFEFunction_).bind(element, coefficients[_0]);
std::get<1>(*localGFEFunction_).bind(element, coefficients[_1]);
LocalGFEFunctionType0 localGFEFunction0(localFiniteElement0, coefficients[_0]);
LocalGFEFunctionType1 localGFEFunction1(localFiniteElement1, coefficients[_1]);
// Bind density to the element
localDensity_->bind(element);
int quadOrder = (element.type().isSimplex()) ? localFiniteElement0.localBasis().order()
: localFiniteElement0.localBasis().order() * gridDim;
// Get a suitable quadrature rule
const auto feOrder = localView.tree().child(_0,0).finiteElement().localBasis().order();
int quadOrder = (element.type().isSimplex()) ? feOrder : feOrder * gridDim;
const auto& quad = QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
......@@ -199,19 +175,19 @@ namespace Dune::GFE {
// A value is needed either directly, or for computing the derivative.
TargetSpace value;
if (localDensity_->dependsOnValue(0) || localDensity_->dependsOnDerivative(0))
value[_0] = localGFEFunction0.evaluate(quadPos);
value[_0] = std::get<0>(*localGFEFunction_).evaluate(quadPos);
if (localDensity_->dependsOnValue(1) || localDensity_->dependsOnDerivative(1))
value[_1] = localGFEFunction1.evaluate(quadPos);
value[_1] = std::get<1>(*localGFEFunction_).evaluate(quadPos);
// Compute the derivatives of the interpolation function factors
typename LocalGFEFunctionType0::DerivativeType derivative0;
typename LocalGFEFunctionType1::DerivativeType derivative1;
typename std::tuple_element_t<0, LocalInterpolationRule>::DerivativeType derivative0;
typename std::tuple_element_t<1, LocalInterpolationRule>::DerivativeType derivative1;
if (localDensity_->dependsOnDerivative(0))
derivative0 = localGFEFunction0.evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
derivative0 = std::get<0>(*localGFEFunction_).evaluateDerivative(quadPos,value[_0]) * geometryJacobianInverse;
if (localDensity_->dependsOnDerivative(1))
derivative1 = localGFEFunction1.evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
derivative1 = std::get<1>(*localGFEFunction_).evaluateDerivative(quadPos,value[_1]) * geometryJacobianInverse;
// Copy the two derivatives into a joint matrix object
// TODO: I am not sure about this. May the densities should get the
......@@ -230,13 +206,18 @@ namespace Dune::GFE {
}
}
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;
}
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
......