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 (47)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1489 additions and 446 deletions
.gitlab-ci.yml
View file @ 2c68ec15
......@@ -8,8 +8,8 @@ before_script: &before
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-fufem.git
- duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-elasticity.git
# ADOLC's comparison operators return `int` instead of `bool` which confuses
# template meta-programming in dune-common...
# ADOL-C's comparison operators prior to version 2.7.2 return `int` instead of `bool`,
# which confuses template meta-programming in dune-common...
.patch-dune-common: &patch-dune-common
- |-
patch -d /duneci/modules/dune-common -p1 <<EOT
......@@ -34,11 +34,12 @@ dune:2.6 gcc:
variables:
DUNECI_BRANCH: releases/2.6-1
dune:2.6 clang:
image: registry.dune-project.org/docker/ci/dune:2.6-debian-10-clang-7-libcpp-17
dune:2.7 gcc:
image: registry.dune-project.org/docker/ci/dune:2.7-debian-10-gcc-8-17
before_script:
- *patch-dune-common
- *before
script: duneci-standard-test
variables:
DUNECI_BRANCH: releases/2.6-1
dune:git gcc:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
......@@ -48,24 +49,26 @@ dune:git gcc:
script: duneci-standard-test
dune:git clang:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
before_script:
- *patch-dune-common
- *before
script: duneci-standard-test
dune:git parmg gcc:
dune:git parmg dune-vtk gcc:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
before_script:
- *patch-dune-common
- *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-vtk.git
script: duneci-standard-test
dune:git parmg clang:
image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
dune:git parmg dune-vtk clang:
image: registry.dune-project.org/docker/ci/dune:git-ubuntu-20.04-clang-10-20
before_script:
- *patch-dune-common
- *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-vtk.git
script: duneci-standard-test
dune.module
View file @ 2c68ec15
......@@ -8,4 +8,4 @@ Version: svn
Maintainer: oliver.sander@tu-dresden.de
#depending on
Depends: dune-grid dune-uggrid dune-istl dune-localfunctions dune-functions dune-solvers dune-fufem dune-elasticity
Suggests: dune-foamgrid dune-parmg
Suggests: dune-foamgrid dune-parmg dune-vtk
dune/gfe/cosseratenergystiffness.hh
View file @ 2c68ec15
......@@ -37,9 +37,9 @@ class LocalFiniteElementFactory
public:
static auto get(const typename Basis::LocalView& localView,
std::integral_constant<std::size_t, i> iType)
-> decltype(localView.tree().child(iType).finiteElement())
-> decltype(localView.tree().child(iType,0).finiteElement())
{
return localView.tree().child(iType).finiteElement();
return localView.tree().child(iType,0).finiteElement();
}
};
......@@ -74,37 +74,6 @@ class CosseratEnergyLocalStiffness
enum {gridDim=GridView::dimension};
enum {dimworld=GridView::dimensionworld};
/** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
{
Dune::FieldMatrix<field_type,dim,dim> result;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
result[i][j] = 0.5 * (A[i][j] + A[j][i]);
return result;
}
/** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
{
Dune::FieldMatrix<field_type,dim,dim> result;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
result[i][j] = 0.5 * (A[i][j] - A[j][i]);
return result;
}
/** \brief Return the square of the trace of a matrix */
template <int N>
static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
{
field_type trace = 0;
for (int i=0; i<N; i++)
trace += A[i][i];
return trace*trace;
}
/** \brief Compute the (row-wise) curl of a matrix R \f$
\param DR The partial derivatives of the matrix R
*/
......@@ -128,8 +97,8 @@ public:
*/
CosseratEnergyLocalStiffness(const Dune::ParameterTree& parameters,
const BoundaryPatch<GridView>* neumannBoundary,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction)
{
......@@ -171,9 +140,9 @@ public:
for (int i=0; i<dim; i++)
UMinus1[i][i] -= 1;
return mu_ * sym(UMinus1).frobenius_norm2()
+ mu_c_ * skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * traceSquared(sym(UMinus1));
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * Dune::GFE::traceSquared(Dune::GFE::sym(UMinus1));
}
/** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
......@@ -219,7 +188,7 @@ public:
RT detU = U.determinant();
return mu_ * sym(UMinus1).frobenius_norm2() + mu_c_ * skew(UMinus1).frobenius_norm2()
return mu_ * Dune::GFE::sym(UMinus1).frobenius_norm2() + mu_c_ * Dune::GFE::skew(UMinus1).frobenius_norm2()
+ (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
}
......@@ -275,9 +244,9 @@ public:
for (int k=0; k<3; k++)
RT_DR3[i][j] += R[k][i] * DR[k][2][j];
return mu_ * sym(RT_DR3).frobenius_norm2()
+ mu_c_ * skew(RT_DR3).frobenius_norm2()
+ mu_*lambda_/(2*mu_+lambda_) * traceSquared(RT_DR3);
return mu_ * Dune::GFE::sym(RT_DR3).frobenius_norm2()
+ mu_c_ * Dune::GFE::skew(RT_DR3).frobenius_norm2()
+ mu_*lambda_/(2*mu_+lambda_) * Dune::GFE::traceSquared(RT_DR3);
}
/** \brief The shell thickness */
......@@ -302,10 +271,10 @@ public:
const BoundaryPatch<GridView>* neumannBoundary_;
/** \brief The function implementing the Neumann data */
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction_;
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
/** \brief The function implementing a volume load */
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
};
template <class Basis, int dim, class field_type>
......@@ -398,12 +367,7 @@ energy(const typename Basis::LocalView& localView,
continue;
// Value of the volume load density at the current position
Dune::FieldVector<double,3> volumeLoadDensity;
if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
else
volumeLoad_->evaluate(element.geometry().global(quad[pt].position()), volumeLoadDensity);
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++)
......@@ -437,12 +401,7 @@ energy(const typename Basis::LocalView& localView,
RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
Dune::FieldVector<double,3> neumannValue;
if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
else
neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
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++)
......@@ -581,12 +540,7 @@ energy(const typename Basis::LocalView& localView,
RealTuple<field_type,dim> deformationValue = localDeformationGFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
Dune::FieldVector<double,3> neumannValue;
if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
else
neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
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++)
......
dune/gfe/geodesicfeassemblerwrapper.hh 0 → 100644
View file @ 2c68ec15
#ifndef GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
#define GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
#include <dune/gfe/mixedgfeassembler.hh>
#include <dune/common/tuplevector.hh>
namespace Dune::GFE {
/** \brief A wrapper that wraps a MixedGFEAssembler into an assembler that does not distinguish between the two finite element spaces
It reimplements the same methods as these two and transfers the structure of the gradient and the Hessian
*/
template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
class
GeodesicFEAssemblerWrapper {
typedef typename Basis::GridView GridView;
//! Dimension of the grid.
enum { gridDim = GridView::dimension };
//! Dimension of the tangent space
enum { blocksize = TargetSpace::TangentVector::dimension };
enum { blocksize0 = MixedSpace0::TangentVector::dimension };
enum { blocksize1 = MixedSpace1::TangentVector::dimension };
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
typedef typename MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>::MatrixType MatrixType;
protected:
MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler_;
public:
const ScalarBasis& basis_;
/** \brief Constructor for a given grid */
GeodesicFEAssemblerWrapper(MixedGFEAssembler<Basis, MixedSpace0, MixedSpace1>* mixedAssembler, ScalarBasis& basis)
: mixedAssembler_(mixedAssembler),
basis_(basis)
{
hessianMixed_ = std::make_unique<MatrixType>();
// The two spaces from the mixed version need to have the same embeddedDim as the Target Space
assert(MixedSpace0::embeddedDim + MixedSpace1::embeddedDim == TargetSpace::embeddedDim);
assert(blocksize0 + blocksize1 == blocksize);
}
/** \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<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern=true) const;
/** \brief Compute the energy of a deformation state */
virtual double computeEnergy(const std::vector<TargetSpace>& sol) const;
/** \brief Get the occupation structure of the Hessian */
virtual void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
private:
Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> splitVector(const std::vector<TargetSpace>& sol) const;
std::unique_ptr<MatrixType> hessianMixed_;
}; // end class
} //end namespace
template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
splitVector(const std::vector<TargetSpace>& sol) const {
using namespace Indices;
// Split the solution into the Deformation and the Rotational part
auto n = basis_.size();
assert (sol.size() == n);
Dune::TupleVector<std::vector<MixedSpace0>,std::vector<MixedSpace1>> solutionSplit;
solutionSplit[_0].resize(n);
solutionSplit[_1].resize(n);
for (int i = 0; i < n; i++) {
solutionSplit[_0][i] = sol[i].r; // Deformation part
solutionSplit[_1][i] = sol[i].q; // Rotational part
}
return solutionSplit;
}
template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
void Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
auto n = basis_.size();
nb.resize(n, n);
//Retrieve occupation structure for the mixed space and convert it to the non-mixed space
//In the mixed space, each index set is for one part of the composite basis
//So: nb00 is for the displacement part, nb11 is for the rotation part and nb01 and nb10 (where nb01^T = nb10) is for the coupling
Dune::MatrixIndexSet nb00, nb01, nb10, nb11;
mixedAssembler_->getMatrixPattern(nb00, nb01, nb10, nb11);
auto size0 = nb00.rows();
auto size1 = nb11.rows();
assert(size0 == size1);
assert(size0 == n);
//After checking if the sizes match, we can just copy over the occupation pattern
//as all occupation patterns work on the same basis combination, so they are all equal
nb = nb00;
}
template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
void Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
Dune::BCRSMatrix<MatrixBlock>& hessian,
bool computeOccupationPattern) const
{
using namespace Dune::TypeTree::Indices;
auto n = basis_.size();
// Get a split up version of the input
auto solutionSplit = splitVector(sol);
// Define the Matrix and the Gradient in Block Stucture
Dune::BlockVector<Dune::FieldVector<double, blocksize0> > gradient0(n);
Dune::BlockVector<Dune::FieldVector<double, blocksize1> > gradient1(n);
if (computeOccupationPattern) {
Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_0]);
neighborsPerVertex.exportIdx((*hessianMixed_)[_0][_1]);
neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_0]);
neighborsPerVertex.exportIdx((*hessianMixed_)[_1][_1]);
neighborsPerVertex.exportIdx(hessian);
}
*hessianMixed_ = 0;
mixedAssembler_->assembleGradientAndHessian(solutionSplit[_0], solutionSplit[_1], gradient0, gradient1, *hessianMixed_, false);
// Transform gradient and hessian back to the non-mixed structure
hessian = 0;
gradient.resize(n);
gradient = 0;
for (int i = 0; i < n; i++) {
for (int j = 0; j < blocksize0 + blocksize1; j++)
gradient[i][j] = j < blocksize0 ? gradient0[i][j] : gradient1[i][j - blocksize0];
}
// All 4 matrices are nxn;
// Each FieldMatrix in (*hessianMixed_)[_0][_0] is blocksize0 x blocksize0
// Each FieldMatrix in (*hessianMixed_)[_1][_0] is blocksize1 x blocksize0
// Each FieldMatrix in (*hessianMixed_)[_0][_1] is blocksize0 x blocksize1
// Each FieldMatrix in (*hessianMixed_)[_1][_1] is blocksize1 x blocksize1
// The hessian will then be a nxn matrix where each FieldMatrix is (blocksize0+blocksize1)x(blocksize0+blocksize1)
for (size_t l = 0; l < blocksize0; l++) {
// Extract Upper Left Block of mixed matrix
for (auto rowIt = (*hessianMixed_)[_0][_0].begin(); rowIt != (*hessianMixed_)[_0][_0].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for(size_t k = 0; k < blocksize0; k++)
hessian[rowIt.index()][colIt.index()][k][l] = (*colIt)[k][l];
// Extract Lower Left Block of mixed matrix
for (auto rowIt = (*hessianMixed_)[_1][_0].begin(); rowIt != (*hessianMixed_)[_1][_0].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for(size_t k = 0; k < blocksize1; k++)
hessian[rowIt.index()][colIt.index()][k + blocksize0][l] = (*colIt)[k][l];
}
for (size_t l = 0; l < blocksize1; l++) {
// Extract Upper Right Block of mixed matrix
for (auto rowIt = (*hessianMixed_)[_0][_1].begin(); rowIt != (*hessianMixed_)[_0][_1].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for(size_t k = 0; k < blocksize0; k++)
hessian[rowIt.index()][colIt.index()][k][l + blocksize0] = (*colIt)[k][l];
// Extract Lower Right Block of mixed matrix
for (auto rowIt = (*hessianMixed_)[_1][_1].begin(); rowIt != (*hessianMixed_)[_1][_1].end(); ++rowIt)
for (auto colIt = rowIt->begin(); colIt != rowIt->end(); ++colIt)
for(size_t k = 0; k < blocksize1; k++)
hessian[rowIt.index()][colIt.index()][k + blocksize0][l + blocksize0] = (*colIt)[k][l];
}
}
template <class Basis, class ScalarBasis, class TargetSpace, class MixedSpace0, class MixedSpace1>
double Dune::GFE::GeodesicFEAssemblerWrapper<Basis, ScalarBasis, TargetSpace,MixedSpace0,MixedSpace1>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
using namespace Dune::TypeTree::Indices;
auto solutionSplit = splitVector(sol);
return mixedAssembler_->computeEnergy(solutionSplit[_0], solutionSplit[_1]);
}
#endif //GLOBAL_GEODESIC_FE_ASSEMBLERWRAPPER_HH
\ No newline at end of file
dune/gfe/linearalgebra.hh
View file @ 2c68ec15
#ifndef LINEAR_ALGEBRA_HH
#define LINEAR_ALGEBRA_HH
#ifndef DUNE_GFE_LINEARALGEBRA_HH
#define DUNE_GFE_LINEARALGEBRA_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
///////////////////////////////////////////////////////////////////////////////////////////
// Various vector-space matrix methods
///////////////////////////////////////////////////////////////////////////////////////////
#if ADOLC_ADOUBLE_H
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
template< int m, int n, int p >
......@@ -171,4 +175,123 @@ Dune::FieldVector<K,m> operator/ ( const Dune::FieldVector<K, m> &A, const K& s)
}
#endif
///////////////////////////////////////////////////////////////////////////////////////////
// Various other matrix methods
///////////////////////////////////////////////////////////////////////////////////////////
namespace Dune {
namespace GFE {
/** \brief Return the trace of a matrix */
template <class T, int n>
static T trace(const FieldMatrix<T,n,n>& A)
{
T trace = 0;
for (int i=0; i<n; i++)
trace += A[i][i];
return trace;
}
/** \brief Return the square of the trace of a matrix */
template <class T, int n>
static T traceSquared(const FieldMatrix<T,n,n>& A)
{
T trace = 0;
for (int i=0; i<n; i++)
trace += A[i][i];
return trace*trace;
}
/** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
template <class T, int n>
static FieldMatrix<T,n,n> sym(const FieldMatrix<T,n,n>& A)
{
FieldMatrix<T,n,n> result;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result[i][j] = 0.5 * (A[i][j] + A[j][i]);
return result;
}
/** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
template <class T, int n>
static FieldMatrix<T,n,n> skew(const FieldMatrix<T,n,n>& A)
{
FieldMatrix<T,n,n> result;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result[i][j] = 0.5 * (A[i][j] - A[j][i]);
return result;
}
/** \brief Compute the deviator of a matrix A */
template <class T, int n>
static FieldMatrix<T,n,n> dev(const FieldMatrix<T,n,n>& A)
{
FieldMatrix<T,n,n> result = A;
auto t = trace(A);
for (int i=0; i<n; i++)
result[i][i] -= t / n;
return result;
}
/** \brief Return the transposed matrix */
template <class T, int n>
static FieldMatrix<T,n,n> transpose(const FieldMatrix<T,n,n>& A)
{
FieldMatrix<T,n,n> result;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result[i][j] = A[j][i];
return result;
}
/** \brief The Frobenius (i.e., componentwise) product of two matrices */
template <class T, int n>
static T frobeniusProduct(const FieldMatrix<T,n,n>& A, const FieldMatrix<T,n,n>& B)
{
T result(0.0);
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result += A[i][j] * B[i][j];
return result;
}
template <class T, int n>
static auto dyadicProduct(const FieldVector<T,n>& A, const FieldVector<T,n>& B)
{
FieldMatrix<T,n,n> result;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result[i][j] = A[i]*B[j];
return result;
}
#if ADOLC_ADOUBLE_H
template <int n>
static auto dyadicProduct(const FieldVector<adouble,n>& A, const FieldVector<double,n>& B)
-> FieldMatrix<adouble,n,n>
{
FieldMatrix<adouble,n,n> result;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
result[i][j] = A[i]*B[j];
return result;
}
#endif
}
}
#endif
dune/gfe/localenergy.hh
View file @ 2c68ec15
......@@ -8,19 +8,20 @@ namespace Dune {
namespace GFE {
/** \brief Base class for energies defined by integrating over one grid element */
template<class Basis, class TargetSpace>
template<class Basis, class... TargetSpaces>
class LocalEnergy
{
public:
using RT = typename std::common_type<typename TargetSpaces::ctype...>::type;
/** \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 typename TargetSpace::ctype
virtual RT
energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const = 0;
const std::vector<TargetSpaces>&... localSolution) const = 0;
/** Empty virtual default destructor
*
......
dune/gfe/localgeodesicfeadolcstiffness.hh
View file @ 2c68ec15
......@@ -106,13 +106,13 @@ energy(const typename Basis::LocalView& localView,
try {
energy = localEnergy_->energy(localView,localASolution);
} catch (Dune::Exception &e) {
trace_off(rank);
trace_off();
throw e;
}
energy >>= pureEnergy;
trace_off(rank);
trace_off();
#if 0
size_t tape_stats[STAT_SIZE];
tapestats(rank,tape_stats); // reading of tape statistics
......
dune/gfe/localintegralenergy.hh 0 → 100644
View file @ 2c68ec15
#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH
#define DUNE_GFE_LOCALINTEGRALENERGY_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/gfe/localenergy.hh>
#include <dune/gfe/realtuple.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/elasticity/materials/localdensity.hh>
namespace Dune::GFE {
/**
\brief Assembles the elastic energy for a single element integrating the localdensity over one element.
This class works similarly to the class Dune::Elasticity::LocalIntegralEnergy, where Dune::Elasticity::LocalIntegralEnergy extends
Dune::Elasticity::LocalEnergy and Dune::GFE::LocalIntegralEnergy extends Dune::GFE::LocalEnergy.
*/
template<class Basis, class... TargetSpaces>
class LocalIntegralEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
{
using LocalView = typename Basis::LocalView;
using GridView = typename LocalView::GridView;
using DT = typename GridView::Grid::ctype;
using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
enum {gridDim=GridView::dimension};
public:
/** \brief Constructor with a local energy density
*/
LocalIntegralEnergy(const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,RT,DT>>& ld)
: localDensity_(ld)
{}
/** \brief Assemble the energy for a single element */
RT energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpaces>&... localSolutions) const
{
static_assert(sizeof...(TargetSpaces) > 0, "LocalIntegralEnergy needs at least one TargetSpace!");
using TargetSpace = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
static_assert( (std::is_same<TargetSpace, RigidBodyMotion<RT,gridDim>>::value)
or (std::is_same<TargetSpace, RealTuple<RT,gridDim>>::value), "The first TargetSpace of LocalIntegralEnergy needs to be RigidBodyMotion or RealTuple!" );
const std::vector<TargetSpace>& localSolution = std::get<0>(std::forward_as_tuple(localSolutions...));
using namespace Dune::Indices;
// powerBasis: grab the finite element of the first child
const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
const auto& element = localView.element();
RT energy = 0;
// store gradients of shape functions and base functions
std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
: localFiniteElement.localBasis().order() * gridDim;
const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
for (const auto& qp : quad)
{
const DT integrationElement = element.geometry().integrationElement(qp.position());
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
// Global position
auto x = element.geometry().global(qp.position());
// Get gradients of shape functions
localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
// compute gradients of base functions
for (size_t i=0; i<gradients.size(); ++i)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
// Deformation gradient
FieldMatrix<RT,gridDim,gridDim> deformationGradient(0);
for (size_t i=0; i<gradients.size(); i++)
for (int j=0; j<gridDim; j++)
deformationGradient[j].axpy(localSolution[i][j], gradients[i]);
// Integrate the energy density
energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
}
return energy;
}
protected:
const std::shared_ptr<Dune::Elasticity::LocalDensity<gridDim,RT,DT>> localDensity_ = nullptr;
};
} // namespace Dune::GFE
#endif //#ifndef DUNE_GFE_LOCALINTEGRALENERGY_HH
dune/gfe/mixedgfeassembler.hh
View file @ 2c68ec15
......@@ -31,11 +31,9 @@ class MixedGFEAssembler {
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize0> > MatrixBlock10;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > MatrixBlock11;
public:
typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixBlock00,MatrixBlock01>,
Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > MatrixType;
protected:
public:
const Basis basis_;
MixedLocalGeodesicFEStiffness<Basis,
......@@ -195,26 +193,38 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
#endif
using namespace Dune::TypeTree::Indices;
const int nDofs0 = localView.tree().child(_0).finiteElement().size();
const int nDofs1 = localView.tree().child(_1).finiteElement().size();
const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
// This loop reads out the pattern for a local matrix; in each element, we have localView.size() degrees of freedom; from the composite and powerbasis layers
// nDofs0 are the degrees of freedom for *one* subspacebasis of the power basis of the displacement part;
// nDofs1 are the degrees of freedom for *one* subspacebasis of the power basis of the rotational part
// this is why the indices (_0,0) and (_1,0) are used: _0 takes the whole displacement part and _1 the whole rotational part; and 0 the first subspacebasis respectively
// Extract local solution
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (int i=0; i<nDofs0+nDofs1; i++)
{
int localIndexI = 0;
if (i < nDofs0) {
auto& node = localView.tree().child(_0,0);
localIndexI = node.localIndex(i);
} else {
auto& node = localView.tree().child(_1,0);
localIndexI = node.localIndex(i-nDofs0);
}
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
auto multiIndex = localIndexSet.index(localIndexI);
#else
if (localView.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localView.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
#endif
auto multiIndex = localView.index(localIndexI);
#endif
//CompositeBasis number is contained in multiIndex[0], the Subspacebasis is contained in multiIndex[2]
//multiIndex[1] contains the actual index
if (multiIndex[0] == 0)
localConfiguration0[i] = configuration0[multiIndex[1]];
else if (multiIndex[0] == 1)
localConfiguration1[i-nDofs0] = configuration1[multiIndex[1]];
}
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
......@@ -228,18 +238,34 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0+nDofs1; i++)
{
int localIndexRow = 0;
if (i < nDofs0) {
auto& node = localView.tree().child(_0,0);
localIndexRow = node.localIndex(i);
} else {
auto& node = localView.tree().child(_1,0);
localIndexRow = node.localIndex(i-nDofs0);
}
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto row = localIndexSet.index(i);
auto row = localIndexSet.index(localIndexRow);
#else
auto row = localView.index(i);
auto row = localView.index(localIndexRow);
#endif
for (int j=0; j<nDofs0+nDofs1; j++ )
{
int localIndexCol = 0;
if (j < nDofs0) {
auto& node = localView.tree().child(_0,0);
localIndexCol = node.localIndex(j);
} else {
auto& node = localView.tree().child(_1,0);
localIndexCol = node.localIndex(j-nDofs0);
}
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto col = localIndexSet.index(j);
auto col = localIndexSet.index(localIndexCol);
#else
auto col = localView.index(j);
auto col = localView.index(localIndexCol);
#endif
if (row[0]==0 and col[0]==0)
......@@ -256,17 +282,10 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
}
// Add local gradient to global gradient
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
gradient0[localIndexSet.index(i)[1]] += localGradient0[i];
else
gradient1[localIndexSet.index(i)[1]] += localGradient1[i-nDofs0];
#else
if (localView.index(i)[0] == 0)
gradient0[localView.index(i)[1]] += localGradient0[i];
if (row[0] == 0)
gradient0[row[1]] += localGradient0[i];
else
gradient1[localView.index(i)[1]] += localGradient1[i-nDofs0];
#endif
gradient1[row[1]] += localGradient1[i-nDofs0];
}
}
......@@ -343,25 +362,33 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
// Number of degrees of freedom on this element
using namespace Dune::TypeTree::Indices;
const int nDofs0 = localView.tree().child(_0).finiteElement().size();
const int nDofs1 = localView.tree().child(_1).finiteElement().size();
const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
std::vector<TargetSpace0> localConfiguration0(nDofs0);
std::vector<TargetSpace1> localConfiguration1(nDofs1);
for (int i=0; i<nDofs0+nDofs1; i++)
{
int localIndexI = 0;
if (i < nDofs0) {
auto& node = localView.tree().child(_0,0);
localIndexI = node.localIndex(i);
} else {
auto& node = localView.tree().child(_1,0);
localIndexI = node.localIndex(i-nDofs0);
}
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
auto multiIndex = localIndexSet.index(localIndexI);
#else
if (localView.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localView.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
#endif
auto multiIndex = localView.index(localIndexI);
#endif
// The CompositeBasis number is contained in multiIndex[0]
// multiIndex[1] contains the actual index
if (multiIndex[0] == 0)
localConfiguration0[i] = configuration0[multiIndex[1]];
else if (multiIndex[0] == 1)
localConfiguration1[i-nDofs0] = configuration1[multiIndex[1]];
}
energy += localStiffness_->energy(localView,
......
dune/gfe/mixedlocalgfeadolcstiffness.hh
View file @ 2c68ec15
......@@ -86,12 +86,13 @@ energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace0>& localConfiguration0,
const std::vector<TargetSpace1>& localConfiguration1) const
{
int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
double pureEnergy;
std::vector<ATargetSpace0> localAConfiguration0(localConfiguration0.size());
std::vector<ATargetSpace1> localAConfiguration1(localConfiguration1.size());
trace_on(1);
trace_on(rank);
adouble energy = 0;
......@@ -122,10 +123,14 @@ energy(const typename Basis::LocalView& localView,
}
using namespace Dune::TypeTree::Indices;
energy = localEnergy_->energy(localView,
try {
energy = localEnergy_->energy(localView,
localAConfiguration0,
localAConfiguration1);
} catch (Dune::Exception &e) {
trace_off();
throw e;
}
energy >>= pureEnergy;
trace_off();
......@@ -197,6 +202,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
std::vector<typename TargetSpace0::TangentVector>& localGradient0,
std::vector<typename TargetSpace1::TangentVector>& localGradient1)
{
int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localConfiguration0, localConfiguration1);
......@@ -222,7 +228,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// Compute gradient
std::vector<double> g(nDoubles);
gradient(1,nDoubles,xp.data(),g.data()); // gradient evaluation
gradient(rank,nDoubles,xp.data(),g.data()); // gradient evaluation
// Copy into Dune type
std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration0.size());
......@@ -338,7 +344,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
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(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
// Copy Hessian into Dune data type
size_t offset0 = nDofs0*embeddedBlocksize0;
......
dune/gfe/mixedriemanniantrsolver.cc
View file @ 2c68ec15
#include "omp.h"
#include <dune/common/bitsetvector.hh>
#include <dune/common/timer.hh>
#include <dune/common/parallel/mpihelper.hh>
......@@ -387,6 +385,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
corr_global[_0].resize(rhs_global[_0].size());
corr_global[_1].resize(rhs_global[_1].size());
corr_global = 0;
bool solvedByInnerSolver = true;
if (rank==0)
{
......@@ -407,35 +406,46 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
mmgStep1->preprocess();
///////////////////////////////
// Solve !
///////////////////////////////
std::cout << "Solve quadratic problem..." << std::endl;
double oldEnergy = 0;
Dune::Timer solutionTimer;
for (int ii=0; ii<innerIterations_; ii++)
{
residual[_0] = rhs_global[_0];
stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
mmgStep0->setRhs(residual[_0]);
mmgStep0->iterate();
residual[_1] = rhs_global[_1];
stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
mmgStep1->setRhs(residual[_1]);
mmgStep1->iterate();
// Compute energy
CorrectionType tmp(corr_global);
stiffnessMatrix.mv(corr_global,tmp);
double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
if (energy > oldEnergy)
//DUNE_THROW(Dune::Exception, "energy increase!");
std::cout << "Warning: energy increase!" << std::endl;
oldEnergy = energy;
int ii = 0;
try {
for (; ii<innerIterations_; ii++) {
residual[_0] = rhs_global[_0];
stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
mmgStep0->setRhs(residual[_0]);
mmgStep0->iterate();
residual[_1] = rhs_global[_1];
stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
mmgStep1->setRhs(residual[_1]);
mmgStep1->iterate();
// Compute energy
CorrectionType tmp(corr_global);
stiffnessMatrix.mv(corr_global,tmp);
double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
if (energy > oldEnergy)
//DUNE_THROW(Dune::Exception, "energy increase!");
std::cout << "Warning: energy increase!" << std::endl;
oldEnergy = energy;
if (corr_global.infinity_norm() < innerTolerance_)
break;
}
} catch (Dune::Exception &e) {
std::cerr << "Error while solving: " << e << std::endl;
solvedByInnerSolver = false;
corr_global = 0;
}
std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds and " << ii << " steps." << std::endl;
//std::cout << "Correction: " << std::endl << corr_global << std::endl;
......@@ -447,73 +457,88 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
//corr = vectorComm.scatter(corr_global);
corr = corr_global;
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
if (corr_global.infinity_norm() < tolerance_) {
if (verbosity_ == NumProc::FULL and rank==0)
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
break;
}
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
double energy = 0;
double modelDecrease = 0;
SolutionType newIterate = x_;
for (size_t j=0; j<newIterate[_0].size(); j++)
newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
for (size_t j=0; j<newIterate[_1].size(); j++)
newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
double energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
energy = mpiHelper.getCollectiveCommunication().sum(energy);
// compute the model decrease
// It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
// Note that rhs = -g
CorrectionType tmp(corr);
hessianMatrix_->mv(corr,tmp);
double modelDecrease = rhs*corr - 0.5 * (corr*tmp);
modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
double relativeModelDecrease = modelDecrease / std::fabs(energy);
if (verbosity_ == NumProc::FULL and rank==0) {
std::cout << "Absolute model decrease: " << modelDecrease
<< ", functional decrease: " << oldEnergy - energy << std::endl;
std::cout << "Relative model decrease: " << relativeModelDecrease
<< ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
}
if (solvedByInnerSolver) {
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
if (corr_global.infinity_norm() < tolerance_) {
if (verbosity_ == NumProc::FULL and rank==0)
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
break;
}
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
for (size_t j=0; j<newIterate[_0].size(); j++)
newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
for (size_t j=0; j<newIterate[_1].size(); j++)
newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
try {
energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
} catch (Dune::Exception &e) {
std::cerr << "Error while computing the energy of the new iterate: " << e << std::endl;
std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
newIterate = x_;
solvedByInnerSolver = false;
energy = oldEnergy;
}
if (solvedByInnerSolver) {
energy = mpiHelper.getCollectiveCommunication().sum(energy);
// compute the model decrease
// It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
// Note that rhs = -g
CorrectionType tmp(corr);
hessianMatrix_->mv(corr,tmp);
modelDecrease = rhs*corr - 0.5 * (corr*tmp);
modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
double relativeModelDecrease = modelDecrease / std::fabs(energy);
if (verbosity_ == NumProc::FULL and rank==0) {
std::cout << "Absolute model decrease: " << modelDecrease
<< ", functional decrease: " << oldEnergy - energy << std::endl;
std::cout << "Relative model decrease: " << relativeModelDecrease
<< ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
}
assert(modelDecrease >= 0);
assert(modelDecrease >= 0);
if (energy >= oldEnergy and rank==0) {
if (this->verbosity_ == NumProc::FULL)
printf("Richtung ist keine Abstiegsrichtung!\n");
}
if (energy >= oldEnergy and rank==0) {
if (this->verbosity_ == NumProc::FULL)
printf("Richtung ist keine Abstiegsrichtung!\n");
}
if (energy >= oldEnergy &&
(std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "Suspecting rounding problems" << std::endl;
if (energy >= oldEnergy &&
(std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
if (this->verbosity_ == NumProc::FULL and rank==0)
std::cout << "Suspecting rounding problems" << std::endl;
if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
if (this->verbosity_ != NumProc::QUIET and rank==0)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
x_ = newIterate;
break;
}
x_ = newIterate;
break;
}
}
}
// //////////////////////////////////////////////
// Check for acceptance of the step
// //////////////////////////////////////////////
if ( (oldEnergy-energy) / modelDecrease > 0.9) {
if ( solvedByInnerSolver && (oldEnergy-energy) / modelDecrease > 0.9) {
// very successful iteration
x_ = newIterate;
......@@ -525,7 +550,7 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
recomputeGradientHessian = true;
} else if ( (oldEnergy-energy) / modelDecrease > 0.01
} else if (solvedByInnerSolver && (oldEnergy-energy) / modelDecrease > 0.01
|| std::abs(oldEnergy-energy) < 1e-12) {
// successful iteration
x_ = newIterate;
......
dune/gfe/neumannenergy.hh 0 → 100644
View file @ 2c68ec15
#ifndef DUNE_GFE_NEUMANNENERGY_HH
#define DUNE_GFE_NEUMANNENERGY_HH
#include <dune/geometry/quadraturerules.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
namespace Dune::GFE {
/**
\brief Assembles the Neumann energy for a single element on the Neumann Boundary using the Neumann Function.
This class works similarly to the class Dune::Elasticity::NeumannEnergy, where Dune::Elasticity::NeumannEnergy extends
Dune::Elasticity::LocalEnergy and Dune::GFE::NeumannEnergy extends Dune::GFE::LocalEnergy.
*/
template<class Basis, class... TargetSpaces>
class NeumannEnergy
: public Dune::GFE::LocalEnergy<Basis,TargetSpaces...>
{
using LocalView = typename Basis::LocalView;
using GridView = typename LocalView::GridView;
using DT = typename GridView::Grid::ctype;
using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpaces...>::RT;
enum {dim=GridView::dimension};
public:
/** \brief Constructor with a set of material parameters
* \param parameters The material parameters
*/
NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView>>& neumannBoundary,
std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<DT,dim>)> neumannFunction)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction)
{}
/** \brief Assemble the energy for a single element */
RT energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpaces>&... localSolutions) const
{
static_assert(sizeof...(TargetSpaces) > 0, "NeumannEnergy needs at least one TargetSpace!");
using namespace Dune::Indices;
using TargetSpace = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
const std::vector<TargetSpace>& localSolution = std::get<0>(std::forward_as_tuple(localSolutions...));
const auto& localFiniteElement = localView.tree().child(_0,0).finiteElement();
const auto& element = localView.element();
RT energy = 0;
for (auto&& intersection : intersections(neumannBoundary_->gridView(), element)) {
if (not neumannBoundary_ or not neumannBoundary_->contains(intersection))
continue;
int quadOrder = localFiniteElement.localBasis().order();
const auto& quad = Dune::QuadratureRules<DT, dim-1>::rule(intersection.type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<DT,dim>& quadPos = intersection.geometryInInside().global(quad[pt].position());
const auto integrationElement = intersection.geometry().integrationElement(quad[pt].position());
// The value of the local function
std::vector<Dune::FieldVector<DT,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
Dune::FieldVector<RT,dim> value(0);
for (size_t i=0; i<localFiniteElement.size(); i++)
for (int j=0; j<dim; j++)
value[j] += shapeFunctionValues[i] * localSolution[i][j];
// Value of the Neumann data at the current position
auto neumannValue = neumannFunction_( intersection.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[i]) * quad[pt].weight() * integrationElement;
}
}
return energy;
}
private:
/** \brief The Neumann boundary */
const std::shared_ptr<BoundaryPatch<GridView>> neumannBoundary_;
/** \brief The function implementing the Neumann data */
std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<DT,dim>)> neumannFunction_;
};
} // namespace Dune::GFE
#endif //#ifndef DUNE_GFE_NEUMANNENERGY_HH
dune/gfe/nonplanarcosseratshellenergy.hh
View file @ 2c68ec15
......@@ -34,119 +34,6 @@ class NonplanarCosseratShellEnergy
enum {gridDim=GridView::dimension};
enum {dimworld=GridView::dimensionworld};
/** \brief Compute the symmetric part of a matrix A, i.e. \f$ \frac 12 (A + A^T) \f$ */
static Dune::FieldMatrix<field_type,dim,dim> sym(const Dune::FieldMatrix<field_type,dim,dim>& A)
{
Dune::FieldMatrix<field_type,dim,dim> result;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
result[i][j] = 0.5 * (A[i][j] + A[j][i]);
return result;
}
/** \brief Compute the antisymmetric part of a matrix A, i.e. \f$ \frac 12 (A - A^T) \f$ */
static Dune::FieldMatrix<field_type,dim,dim> skew(const Dune::FieldMatrix<field_type,dim,dim>& A)
{
Dune::FieldMatrix<field_type,dim,dim> result;
for (int i=0; i<dim; i++)
for (int j=0; j<dim; j++)
result[i][j] = 0.5 * (A[i][j] - A[j][i]);
return result;
}
/** \brief Compute the deviator of a matrix A */
static Dune::FieldMatrix<field_type,dim,dim> dev(const Dune::FieldMatrix<field_type,dim,dim>& A)
{
Dune::FieldMatrix<field_type,dim,dim> result = A;
auto t = trace(A);
for (int i=0; i<dim; i++)
result[i][i] -= t / dim;
return result;
}
/** \brief Return the trace of a matrix */
template <class T, int N>
static T trace(const Dune::FieldMatrix<T,N,N>& A)
{
T trace = 0;
for (int i=0; i<N; i++)
trace += A[i][i];
return trace;
}
/** \brief Return the square of the trace of a matrix */
template <int N>
static field_type traceSquared(const Dune::FieldMatrix<field_type,N,N>& A)
{
field_type trace = 0;
for (int i=0; i<N; i++)
trace += A[i][i];
return trace*trace;
}
template <class T, int N>
static T frobeniusProduct(const Dune::FieldMatrix<T,N,N>& A, const Dune::FieldMatrix<T,N,N>& B)
{
T result(0.0);
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result += A[i][j] * B[i][j];
return result;
}
template <int N>
static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<adouble,N>& B)
-> Dune::FieldMatrix<adouble,N,N>
{
Dune::FieldMatrix<adouble,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = A[i]*B[j];
return result;
}
template <int N>
static auto dyadicProduct(const Dune::FieldVector<double,N>& A, const Dune::FieldVector<double,N>& B)
-> Dune::FieldMatrix<double,N,N>
{
Dune::FieldMatrix<double,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = A[i]*B[j];
return result;
}
template <int N>
static auto dyadicProduct(const Dune::FieldVector<adouble,N>& A, const Dune::FieldVector<double,N>& B)
-> Dune::FieldMatrix<adouble,N,N>
{
Dune::FieldMatrix<adouble,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = A[i]*B[j];
return result;
}
template <class T, int N>
static Dune::FieldMatrix<T,N,N> transpose(const Dune::FieldMatrix<T,N,N>& A)
{
Dune::FieldMatrix<T,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
result[i][j] = A[j][i];
return result;
}
public:
/** \brief Constructor with a set of material parameters
......@@ -155,8 +42,8 @@ public:
NonplanarCosseratShellEnergy(const Dune::ParameterTree& parameters,
const std::vector<UnitVector<double,3> >& vertexNormals,
const BoundaryPatch<GridView>* neumannBoundary,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction,
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad)
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction,
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad)
: vertexNormals_(vertexNormals),
neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction),
......@@ -192,19 +79,20 @@ public:
RT W_mixt(const Dune::FieldMatrix<field_type,3,3>& S, const Dune::FieldMatrix<field_type,3,3>& T) const
{
return mu_ * frobeniusProduct(sym(S), sym(T))
+ mu_c_ * frobeniusProduct(skew(S), skew(T))
+ lambda_ * mu_ / (lambda_ + 2*mu_) * trace(S) * trace(T);
return mu_ * Dune::GFE::frobeniusProduct(Dune::GFE::sym(S), Dune::GFE::sym(T))
+ mu_c_ * Dune::GFE::frobeniusProduct(Dune::GFE::skew(S), Dune::GFE::skew(T))
+ lambda_ * mu_ / (lambda_ + 2*mu_) * Dune::GFE::trace(S) * Dune::GFE::trace(T);
}
RT W_mp(const Dune::FieldMatrix<field_type,3,3>& S) const
{
return mu_ * sym(S).frobenius_norm2() + mu_c_ * skew(S).frobenius_norm2() + lambda_ * 0.5 * traceSquared(S);
return mu_ * Dune::GFE::sym(S).frobenius_norm2() + mu_c_ * Dune::GFE::skew(S).frobenius_norm2() + lambda_ * 0.5 * Dune::GFE::traceSquared(S);
}
RT W_curv(const Dune::FieldMatrix<field_type,3,3>& S) const
{
return mu_ * L_c_ * L_c_ * (b1_ * dev(sym(S)).frobenius_norm2() + b2_ * skew(S).frobenius_norm2() + b3_ * traceSquared(S));
return mu_ * L_c_ * L_c_ * (b1_ * Dune::GFE::dev(Dune::GFE::sym(S)).frobenius_norm2()
+ b2_ * Dune::GFE::skew(S).frobenius_norm2() + b3_ * Dune::GFE::traceSquared(S));
}
/** \brief The shell thickness */
......@@ -229,10 +117,10 @@ public:
const BoundaryPatch<GridView>* neumannBoundary_;
/** \brief The function implementing the Neumann data */
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > neumannFunction_;
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> neumannFunction_;
/** \brief The function implementing a volume load */
const std::shared_ptr<Dune::VirtualFunction<Dune::FieldVector<double,dimworld>, Dune::FieldVector<double,3> > > volumeLoad_;
const std::function<Dune::FieldVector<double,3>(Dune::FieldVector<double,dimworld>)> volumeLoad_;
};
template <class Basis, int dim, class field_type>
......@@ -297,7 +185,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<field_type,dim,dim> R;
value.q.matrix(R);
auto RT = transpose(R);
auto RT = Dune::GFE::transpose(R);
Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
......@@ -328,7 +216,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldMatrix<double,3,3> a(0);
for (int alpha=0; alpha<gridDim; alpha++)
a += dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
a += Dune::GFE::dyadicProduct(aCovariant[alpha], aContravariant[alpha]);
auto a00 = aCovariant[0] * aCovariant[0];
auto a01 = aCovariant[0] * aCovariant[1];
......@@ -342,7 +230,7 @@ energy(const typename Basis::LocalView& localView,
for (int alpha=0; alpha<2; alpha++)
for (int beta=0; beta<2; beta++)
c += sqrt(aScalar) * eps[alpha][beta] * dyadicProduct(aContravariant[alpha], aContravariant[beta]);
c += aScalar * eps[alpha][beta] * Dune::GFE::dyadicProduct(aContravariant[alpha], aContravariant[beta]);
// Second fundamental form
// The derivative of the normal field
......@@ -354,14 +242,14 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldVector<double,3> vec;
for (int i=0; i<3; i++)
vec[i] = normalDerivative[i][alpha];
b -= dyadicProduct(vec, aContravariant[alpha]);
b -= Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
// Gauss curvature
auto K = b.determinant();
// Mean curvatue
auto H = 0.5 * trace(b);
auto H = 0.5 * Dune::GFE::trace(b);
//////////////////////////////////////////////////////////
// Strain tensors
......@@ -375,7 +263,7 @@ energy(const typename Basis::LocalView& localView,
Dune::FieldVector<field_type,3> vec;
for (int i=0; i<3; i++)
vec[i] = derivative[i][alpha];
grad_s_m += dyadicProduct(vec, aContravariant[alpha]);
grad_s_m += Dune::GFE::dyadicProduct(vec, aContravariant[alpha]);
}
Ee = RT * grad_s_m - a;
......@@ -389,7 +277,7 @@ energy(const typename Basis::LocalView& localView,
for (int j=0; j<3; j++)
tmp[i][j] = DR[i][j][alpha];
auto tmp2 = RT * tmp;
Ke += dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
Ke += Dune::GFE::dyadicProduct(SkewMatrix<field_type,3>(tmp2).axial(), aContravariant[alpha]);
}
//////////////////////////////////////////////////////////
......@@ -418,12 +306,7 @@ energy(const typename Basis::LocalView& localView,
continue;
// Value of the volume load density at the current position
Dune::FieldVector<double,3> volumeLoadDensity;
if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_))
std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(volumeLoad_)->evaluateLocal(element, quadPos, volumeLoadDensity);
else
volumeLoad_->evaluate(geometry.global(quad[pt].position()), volumeLoadDensity);
Dune::FieldVector<double,3> volumeLoadDensity = volumeLoad_(geometry.global(quad[pt].position()));
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
......@@ -456,12 +339,7 @@ energy(const typename Basis::LocalView& localView,
RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
Dune::FieldVector<double,3> neumannValue;
if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_))
std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,Dune::FieldVector<double,3> > >(neumannFunction_)->evaluateLocal(element, quadPos, neumannValue);
else
neumannFunction_->evaluate(it.geometry().global(quad[pt].position()), neumannValue);
Dune::FieldVector<double,3> 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++)
......
dune/gfe/parallel/globalmapper.hh
View file @ 2c68ec15
......@@ -28,12 +28,11 @@ namespace Dune {
Master m = Dune::ParMG::entityMasterRank(basis.gridView(), [&](int, int codim) -> bool {
return dofmap.codimSet().test(codim);
});
DofMaster dofmaster = Dune::ParMG::dofMaster(basis, m);
coarseGlobalDof_ = globalDof(basis,m, dofmap);
globalDof_ = globalDof(basis,m, dofmap);
// total number of degrees of freedom
size_ = coarseGlobalDof_.size;
size_ = globalDof_.size;
}
#else
GlobalMapper(const Basis& basis)
......@@ -46,7 +45,7 @@ namespace Dune {
/** \brief Given a local index, retrieve its index globally unique over all processes. */
Index index(const int& localIndex) const {
#if HAVE_DUNE_PARMG
return coarseGlobalDof_.globalDof[localIndex];
return globalDof_.globalDof[localIndex];
#else
return localIndex;
#endif
......@@ -58,7 +57,7 @@ namespace Dune {
}
#if HAVE_DUNE_PARMG
ParMG::GlobalDof coarseGlobalDof_;
ParMG::GlobalDof globalDof_;
#endif
std::size_t size_;
......
dune/gfe/realtuple.hh
View file @ 2c68ec15
......@@ -66,6 +66,17 @@ public:
return *this;
}
RealTuple& operator-=(const Dune::FieldVector<T,N>& other) {
data_ -= other;
return *this;
}
template <class T2>
RealTuple& operator-=(const RealTuple<T2,N>& other) {
data_ -= other.data_;
return *this;
}
/** \brief Assigment from RealTuple with different type -- used for automatic differentiation with ADOL-C */
template <class T2>
RealTuple& operator <<= (const RealTuple<T2,N>& other) {
......@@ -74,6 +85,11 @@ public:
return *this;
}
/** \brief Const random-access operator*/
T operator[] (const size_t indexVariable ) const {
return data_[indexVariable];
}
/** \brief Rebind the RealTuple to another coordinate type */
template<class U>
struct rebind
......
dune/gfe/riemannianpnsolver.cc 0 → 100644
View file @ 2c68ec15
This diff is collapsed.
dune/gfe/riemannianpnsolver.hh 0 → 100644
View file @ 2c68ec15
#ifndef RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
#define RIEMANNIAN_PROXIMAL_NEWTON_SOLVER_HH
#include <vector>
#include <dune/common/bitsetvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/cholmodsolver.hh>
#include <dune/gfe/periodic1dpq1nodalbasis.hh>
#include "riemanniantrsolver.hh"
#include <dune/grid/utility/globalindexset.hh>
#include <dune/gfe/parallel/globalmapper.hh>
#include <dune/gfe/parallel/globalp1mapper.hh>
#include <dune/gfe/parallel/globalp2mapper.hh>
#include <dune/gfe/parallel/p2mapper.hh>
/** \brief Riemannian proximal-newton solver for geodesic finite-element problems */
template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
class RiemannianProximalNewtonSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
{
typedef typename Basis::GridView::Grid GridType;
const static int blocksize = TargetSpace::TangentVector::dimension;
const static int gridDim = GridType::dimension;
// Centralize the field type here
typedef double field_type;
// Some types that I need
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
typedef std::vector<TargetSpace> SolutionType;
#if HAVE_MPI
typedef typename MapperFactory<Basis>::GlobalMapper GlobalMapper;
typedef typename MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianProximalNewtonSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
{
std::size_t finalIteration;
field_type finalEnergy;
};
public:
RiemannianProximalNewtonSolver()
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
{}
/** \brief Set up the solver using a choldmod solver as the inner solver */
void setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
int maxProximalNewtonSteps,
double initialRegularization,
bool instrumented);
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
innerSolver_->ignoreNodes_ = ignoreNodes_;
}
void solve();
void setInitialSolution(const SolutionType& x) DUNE_DEPRECATED {
x_ = x;
}
void setInitialIterate(const SolutionType& x) {
x_ = x;
}
SolutionType getSol() const {return x_;}
const Statistics& getStatistics() const {return statistics_;}
protected:
#if HAVE_MPI
std::unique_ptr<GlobalMapper> globalMapper_;
#endif
/** \brief The grid */
const GridType* grid_;
/** \brief The solution vector */
SolutionType x_;
/** \brief The initial regularization parameter for the proximal newton step */
double initialRegularization_;
double tolerance_;
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
const Assembler* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<typename Dune::Solvers::CholmodSolver<MatrixType,CorrectionType>> innerSolver_;
/** \brief The Dirichlet nodes */
const Dune::BitSetVector<blocksize>* ignoreNodes_;
/** \brief The norm used to measure convergence for statistics*/
std::shared_ptr<H1SemiNorm<CorrectionType> > h1SemiNorm_;
/** \brief An L2-norm, really. The H1SemiNorm class is badly named */
std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
/** \brief If set to true we log convergence speed and other stuff */
bool instrumented_;
/** \brief Store information about solver runs for unit testing */
Statistics statistics_;
#if DUNE_VERSION_LT(DUNE_GEOMETRY, 2, 7)
std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > A;
std::shared_ptr<Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > > massMatrix;
#endif
};
#include "riemannianpnsolver.cc"
#endif
dune/gfe/riemanniantrsolver.cc
View file @ 2c68ec15
......@@ -28,10 +28,10 @@
#include <dune/gfe/parallel/vectorcommunicator.hh>
#endif
template <class Basis, class TargetSpace>
void RiemannianTrustRegionSolver<Basis,TargetSpace>::
template <class Basis, class TargetSpace, class Assembler>
void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
setup(const GridType& grid,
const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
......@@ -306,8 +306,8 @@ setup(const GridType& grid,
}
template <class Basis, class TargetSpace>
void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
template <class Basis, class TargetSpace, class Assembler>
void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
{
int rank = grid_->comm().rank();
......@@ -403,6 +403,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
*hessianMatrix_,
i==0 // assemble occupation pattern only for the first call
);
std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
rhs *= -1; // The right hand side is the _negative_ gradient
......@@ -418,11 +419,12 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
if ((*mgStep->ignoreNodes_)[j][k]) // global Dirichlet nodes set
gradient[j][k] = 0;
if (this->verbosity_ == Solver::FULL and rank==0) {
std::cout << "Gradient operator norm: " << l2Norm_->operator()(gradient) << std::endl;
std::cout << "Gradient norm: " << gradient.two_norm() << std::endl;
}
if (this->verbosity_ == Solver::FULL and rank==0)
std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
if (this->verbosity_ == Solver::FULL)
std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
std::cout << "Oveall assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
totalAssemblyTime += gradientTimer.elapsed();
// Transfer matrix data
......@@ -460,15 +462,14 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
} catch (Dune::Exception &e) {
std::cerr << "Error while solving: " << e << std::endl;
solved = false;
corr_global = 0;
}
std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
totalSolverTime += solutionTimer.elapsed();
if (mgStep && solved)
if (mgStep && solved) {
corr_global = mgStep->getSol();
//std::cout << "Correction: " << std::endl << corr_global << std::endl;
std::cout << "Two norm of the correction: " << corr_global.two_norm() << std::endl;
}
}
// Distribute solution
......@@ -476,7 +477,13 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
std::cout << "Transfer solution back to root process ..." << std::endl;
#if HAVE_MPI
corr = vectorComm.scatter(corr_global);
solved = grid_->comm().min(solved);
if (solved) {
corr = vectorComm.scatter(corr_global);
} else {
corr_global = 0;
corr = 0;
}
#else
corr = corr_global;
#endif
......@@ -582,11 +589,14 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
} catch (Dune::Exception &e) {
std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
newIterate = x_;
solved = false;
energy = oldEnergy;
}
if (solved) {
solved = grid_->comm().min(solved);
if (!solved) {
newIterate = x_;
energy = oldEnergy;
} else {
energy = grid_->comm().sum(energy);
// compute the model decrease
......
dune/gfe/riemanniantrsolver.hh
View file @ 2c68ec15
......@@ -69,7 +69,7 @@ struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,3> >
};
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
template <class Basis, class TargetSpace>
template <class Basis, class TargetSpace, class Assembler = GeodesicFEAssembler<Basis,TargetSpace>>
class RiemannianTrustRegionSolver
: public IterativeSolver<std::vector<TargetSpace>,
Dune::BitSetVector<TargetSpace::TangentVector::dimension> >
......@@ -115,7 +115,7 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const GeodesicFEAssembler<Basis, TargetSpace>* assembler,
const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
......@@ -188,7 +188,7 @@ protected:
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
const GeodesicFEAssembler<Basis, TargetSpace>* assembler_;
const Assembler* assembler_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
......
dune/gfe/sumcosseratenergy.hh deleted 100644 → 0
View file @ 0d02ceed
#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH
#define DUNE_GFE_SUMCOSSERATENERGY_HH
#include <dune/elasticity/assemblers/localfestiffness.hh>
#include <dune/gfe/localenergy.hh>
namespace Dune {
namespace GFE {
template<class Basis, class TargetSpace, class field_type=double>
class SumCosseratEnergy
: public GFE::LocalEnergy<Basis, TargetSpace>
{
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype ctype;
typedef typename Basis::LocalView::Tree::FiniteElement LocalFiniteElement;
typedef typename GridView::ctype DT;
typedef typename TargetSpace::ctype RT;
enum {dim=GridView::dimension};
public:
/** \brief Constructor with elastic energy and cosserat energy
* \param elasticEnergy The elastic energy
* \param cosseratEnergy The cosserat energy
*/
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
SumCosseratEnergy(std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
SumCosseratEnergy(std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#else
SumCosseratEnergy(std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy,
#endif
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace>> cosseratEnergy)
: elasticEnergy_(elasticEnergy),
cosseratEnergy_(cosseratEnergy)
{}
/** \brief Assemble the energy for a single element */
RT energy (const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution) const
{
auto&& element = localView.element();
auto&& localFiniteElement = localView.tree().finiteElement();
std::vector<Dune::FieldVector<field_type,dim>> localConfiguration(localSolution.size());
for (int i = 0; i < localSolution.size(); i++) {
localConfiguration[i] = localSolution[i].r;
}
return elasticEnergy_->energy(element, localFiniteElement, localConfiguration) + cosseratEnergy_->energy(localView, localSolution);
}
private:
#if DUNE_VERSION_LT(DUNE_ELASTICITY, 2, 7)
std::shared_ptr<LocalFEStiffness<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#elif DUNE_VERSION_GTE(DUNE_ELASTICITY, 2, 8)
std::shared_ptr<Dune::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#else
std::shared_ptr<Elasticity::LocalEnergy<GridView,LocalFiniteElement,std::vector<Dune::FieldVector<field_type,dim> > > > elasticEnergy_;
#endif
std::shared_ptr<GFE::LocalEnergy<Basis, TargetSpace> > cosseratEnergy_;
};
} // namespace GFE
} // namespace Dune
#endif //#ifndef DUNE_GFE_SUMCOSSERATENERGY_HH