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
Hide whitespace changes
Inline Side-by-side
Showing
with 660 additions and 1785 deletions
dune/gfe/localgfetestfunctionbasis.hh deleted 100644 → 0
View file @ 5db924c4
#ifndef LOCAL_GFE_TEST_FUNCTION_HH
#define LOCAL_GFE_TEST_FUNCTION_HH
#include <array>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/linearalgebra.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localbasis.hh>
// forward declaration
template <class LocalFiniteElement, class TargetSpace>
class LocalGfeTestFunctionBasis;
template <class LocalFiniteELement, class TargetSpace>
class LocalGfeTestFunctionInterpolation;
/** \brief The local gfe test function finite element
*
* \tparam LagrangeLfe - A Lagrangian finite element whose shape functions define the interpolation weights
* \tparam TargetSpace - The manifold the tangent spaces this basis for from belong to.
*/
template <class LagrangeLfe, class TargetSpace>
class LocalGfeTestFunctionFiniteElement
{
typedef typename LagrangeLfe::Traits::LocalBasisType::Traits LagrangeBasisTraits;
typedef LocalGfeTestFunctionBasis<LagrangeLfe, TargetSpace> LocalBasis;
typedef LocalGfeTestFunctionInterpolation<LagrangeLfe, TargetSpace> LocalInterpolation;
public:
//! Traits
typedef Dune::LocalFiniteElementTraits<LocalBasis,typename LagrangeLfe::Traits::LocalCoefficientsType, LocalInterpolation> Traits;
/** Construct local finite element from the base coefficients and Lagrange local finite element.
*
* \param lfe - The Lagrange local finite element.
* \param baseCoeff - The coefficients of the base points the tangent spaces live at.
*/
LocalGfeTestFunctionFiniteElement(const LagrangeLfe& lfe, const std::vector<TargetSpace> baseCoeff) :
baseCoeff_(baseCoeff),
basis_(lfe, baseCoeff_),
coefficients_(lfe.localCoefficients()),
gt_(lfe.type())
{}
/** \brief Get reference to the local basis.*/
const typename Traits::LocalBasisType& localBasis () const
{
return basis_;
}
/** \brief Get reference to the local coefficients. */
const typename Traits::LocalCoefficientsType& localCoefficients () const
{
return coefficients_;
}
/** \brief Get reference to the local interpolation handler. */
const typename Traits::LocalInterpolationType& localInterpolation () const
{
return interpolation_;
}
/** \brief Get the element type this finite element lives on. */
Dune::GeometryType type () const
{
return gt_;
}
/** \brief Get base coefficients. */
const std::vector<TargetSpace>& getBaseCoefficients() const {return baseCoeff_;}
private:
const std::vector<TargetSpace> baseCoeff_;
LocalBasis basis_;
const typename Traits::LocalCoefficientsType coefficients_;
LocalInterpolation interpolation_;
Dune::GeometryType gt_;
};
/** \brief A local basis of the first variations of a given geodesic finite element function.
*
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
* \tparam TargetSpace The manifold that the function takes its values in
*
* Note that the shapefunctions of this local basis are given blockwise. Each dof corresponds to a local basis of
* the tangent space at that dof. Thus the methods return a vector of arrays.
*/
template <class LocalFiniteElement, class TargetSpace>
class LocalGfeTestFunctionBasis
{
typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LagrangeBasisTraits;
static const int dim = LagrangeBasisTraits::dimDomain;
typedef typename LagrangeBasisTraits::DomainFieldType ctype;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
static const int spaceDim = TargetSpace::TangentVector::dimension;
public:
//! The local basis traits
typedef Dune::LocalBasisTraits<ctype, dim, Dune::FieldVector<ctype,dim>,
typename EmbeddedTangentVector::value_type, embeddedDim, std::array<EmbeddedTangentVector,spaceDim>,
std::array<Dune::FieldMatrix<ctype, embeddedDim, dim>,spaceDim> > Traits;
/** \brief Constructor
*/
LocalGfeTestFunctionBasis(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& baseCoefficients)
: localGFEFunction_(localFiniteElement, baseCoefficients)
{}
/** \brief The number of Lagrange points, NOT the number of basis functions */
unsigned int size() const
{
return localGFEFunction_.size();
}
/** \brief Evaluate all shape functions at the given point */
void evaluateFunction(const typename Traits::DomainType& local,
std::vector<typename Traits::RangeType>& out) const;
/** \brief Evaluate the derivatives of all shape functions function */
void evaluateJacobian(const typename Traits::DomainType& local,
std::vector<typename Traits::JacobianType>& out) const;
/** \brief Polynomial order */
unsigned int order() const
{
return localGFEFunction_.localFiniteElement_.order();
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
*/
const LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace> localGFEFunction_;
};
template <class LocalFiniteElement, class TargetSpace>
void LocalGfeTestFunctionBasis<LocalFiniteElement,TargetSpace>::evaluateFunction(const typename Traits::DomainType& local,
std::vector<typename Traits::RangeType>& out) const
{
out.resize(size());
for (size_t i=0; i<size(); i++) {
Dune::FieldMatrix< double, embeddedDim, embeddedDim > derivative;
/** \todo This call internally keeps computing the value of the gfe function at 'local'.
* This is expensive. Eventually we should precompute it once and reused the result. */
localGFEFunction_.evaluateDerivativeOfValueWRTCoefficient (local, i, derivative);
Dune::FieldMatrix<ctype,spaceDim,embeddedDim> basisVectors = localGFEFunction_.coefficient(i).orthonormalFrame();
for (int j=0; j<spaceDim; j++)
derivative.mv(basisVectors[j], out[i][j]);
}
}
template <class LocalFiniteElement, class TargetSpace>
void LocalGfeTestFunctionBasis<LocalFiniteElement,TargetSpace>::evaluateJacobian(const typename Traits::DomainType& local,
std::vector<typename Traits::JacobianType>& out) const
{
out.resize(size());
for (size_t i=0; i<size(); i++) {
/** \todo This call internally keeps computing the value of the gfe function at 'local'.
* This is expensive. Eventually we should precompute it once and reused the result. */
Tensor3< double, embeddedDim, embeddedDim, dim > derivative;
localGFEFunction_.evaluateDerivativeOfGradientWRTCoefficient (local, i, derivative);
Dune::FieldMatrix<ctype,spaceDim,embeddedDim> basisVectors = localGFEFunction_.coefficient(i).orthonormalFrame();
for (int j=0; j<spaceDim; j++) {
out[i][j] = 0;
// Contract the second index of the derivative with the tangent vector at the i-th Lagrange point.
// Add that to the result.
for (int k=0; k<embeddedDim; k++)
for (int l=0; l<embeddedDim; l++)
for (int m=0; m<dim; m++)
out[i][j][k][m] += derivative[k][l][m] * basisVectors[j][l];
}
}
}
template <class LocalFiniteElement, class TargetSpace>
class LocalGfeTestFunctionInterpolation
{};
#endif
dune/gfe/localprojectedfefunction.hh deleted 100644 → 0
View file @ 5db924c4
#ifndef DUNE_GFE_LOCALPROJECTEDFEFUNCTION_HH
#define DUNE_GFE_LOCALPROJECTEDFEFUNCTION_HH
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/gfe/linearalgebra.hh>
#include <dune/gfe/polardecomposition.hh>
#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/rotation.hh>
namespace Dune {
namespace GFE {
/** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold
*
* \tparam dim Dimension of the reference element
* \tparam ctype Type used for coordinates on the reference element
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
* \tparam TargetSpace The manifold that the function takes its values in
* \tparam conforming Geometrical conformity of the functions (omits projections when set to false)
*/
template <int dim, class ctype, class LocalFiniteElement, class TS, bool conforming=true>
class LocalProjectedFEFunction
{
public:
using TargetSpace=TS;
private:
typedef typename TargetSpace::ctype RT;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
public:
/** \brief The type used for derivatives */
typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
/** \brief Constructor
* \param localFiniteElement A Lagrangian finite element that provides the interpolation points
* \param coefficients Values of the function at the Lagrange points
*/
LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& coefficients)
: localFiniteElement_(localFiniteElement),
coefficients_(coefficients)
{
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
/** \brief Rebind the FEFunction to another TargetSpace */
template<class U>
struct rebind
{
using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U,conforming>;
};
/** \brief The number of Lagrange points */
unsigned int size() const
{
return localFiniteElement_.localBasis().size();
}
/** \brief The type of the reference element */
Dune::GeometryType type() const
{
return localFiniteElement_.type();
}
/** \brief The scalar finite element used as the interpolation weights
*
* \note This method was added for InterpolationDerivatives, which needs it
* to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
* number type. This is not optimal, because the localFiniteElement
* really is an implementation detail of LocalGeodesicFEFunction and
* should not be needed just to copy an entire object. Other non-Euclidean
* interpolation rules may not have such a finite element at all.
* Therefore, this method may disappear again eventually.
*/
const LocalFiniteElement& localFiniteElement() const
{
return localFiniteElement_;
}
/** \brief Evaluate the function */
auto evaluate(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Evaluate the derivative of the function */
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
* \param local Local coordinates in the reference element where to evaluate the derivative
* \param q Value of the local gfe function at 'local'. If you provide something wrong here the result will be wrong, too!
*
* \note This method is only usable in the conforming setting, because it requires the caller
* to hand over the interpolation value as a TargetSpace object.
*/
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const;
/** \brief Evaluate the value and the derivative of the interpolation function
*
* \return A std::pair containing the value and the first derivative of the interpolation function.
* If the interpolation is conforming then the first member of the pair will be a TargetSpace.
* Otherwise it will be a RealTuple.
*/
auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Get the i'th base coefficient. */
const TargetSpace& coefficient(int i) const
{
return coefficients_[i];
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
* \todo We really only need the local basis
*/
const LocalFiniteElement& localFiniteElement_;
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
};
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace, bool conforming>
auto LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
typename TargetSpace::CoordinateType c(0);
for (size_t i=0; i<coefficients_.size(); i++)
c.axpy(w[i][0], coefficients_[i].globalCoordinates());
if constexpr (conforming)
return TargetSpace::projectOnto(c);
else
return (RealTuple<RT, TargetSpace::CoordinateType::dimension>)c;
}
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
if constexpr(conforming)
{
// the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
// Actually compute the derivative
return evaluateDerivative(local, q);
}
else {
std::vector<Dune::FieldMatrix<ctype, 1, dim> > wDer;
localFiniteElement_.localBasis().evaluateJacobian(local, wDer);
Dune::FieldMatrix<RT, embeddedDim, dim> derivative(0);
for (size_t i = 0; i < embeddedDim; i++)
for (size_t j = 0; j < dim; j++)
for (size_t k = 0; k < coefficients_.size(); k++)
derivative[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
return derivative;
}
}
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::DerivativeType
LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace& q) const
{
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
typename TargetSpace::CoordinateType embeddedInterpolation(0);
for (size_t i=0; i<coefficients_.size(); i++)
embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());
Dune::FieldMatrix<RT,embeddedDim,dim> derivative(0);
for (size_t i=0; i<embeddedDim; i++)
for (size_t j=0; j<dim; j++)
for (size_t k=0; k<coefficients_.size(); k++)
derivative[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
auto derivativeOfProjection = TargetSpace::derivativeOfProjection(embeddedInterpolation);
return derivativeOfProjection*derivative;
}
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace,bool conforming>
auto
LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace,conforming>::
evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
// Construct the type of the result -- it depends on whether the interpolation
// is conforming or not.
using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;
std::pair<Value,DerivativeType> result;
///////////////////////////////////////////////////////////
// Compute the value of the interpolation function
///////////////////////////////////////////////////////////
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
typename TargetSpace::CoordinateType embeddedInterpolation(0);
for (size_t i=0; i<coefficients_.size(); i++)
embeddedInterpolation.axpy(w[i][0], coefficients_[i].globalCoordinates());
if constexpr (conforming)
result.first = TargetSpace::projectOnto(embeddedInterpolation);
else
result.first = (RealTuple<RT, TargetSpace::CoordinateType::dimension>)embeddedInterpolation;
///////////////////////////////////////////////////////////
// Compute the derivative of the interpolation function
///////////////////////////////////////////////////////////
// Compute the interpolation in the surrounding space
std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
result.second = 0.0;
for (size_t i=0; i<embeddedDim; i++)
for (size_t j=0; j<dim; j++)
for (size_t k=0; k<coefficients_.size(); k++)
result.second[i][j] += wDer[k][0][j] * coefficients_[k].globalCoordinates()[i];
// The derivative of the projection onto the manifold
if constexpr(conforming)
result.second = TargetSpace::derivativeOfProjection(embeddedInterpolation) * result.second;
return result;
}
/** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)
*
* \tparam dim Dimension of the reference element
* \tparam ctype Type used for coordinates on the reference element
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
*/
template <int dim, class ctype, class LocalFiniteElement, class field_type, bool conforming>
class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3>,conforming>
{
public:
typedef Rotation<field_type,3> TargetSpace;
private:
typedef typename TargetSpace::ctype RT;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
/**
* \param A The argument of the projection
* \param polar The image of the projection, i.e., the polar factor of A
*/
static std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> derivativeOfProjection(const FieldMatrix<field_type,3,3>& A,
FieldMatrix<field_type,3,3>& polar)
{
std::array<std::array<FieldMatrix<field_type,3,3>, 3>, 3> result;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
result[i][j][k][l] = (i==k) and (j==l);
polar = A;
// Use Heron's method
const size_t maxIterations = 100;
double tol = 0.001;
for (size_t iteration=0; iteration<maxIterations; iteration++)
{
auto oldPolar = polar;
auto polarInvert = polar;
polarInvert.invert();
for (size_t i=0; i<polar.N(); i++)
for (size_t j=0; j<polar.M(); j++)
polar[i][j] = 0.5 * (polar[i][j] + polarInvert[j][i]);
// Alternative name to align the code better with a description in a math text
const auto& dQT = result;
// Multiply from the right with Q^{-T}
decltype(result) tmp2;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
tmp2[i][j][k][l] = 0.0;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
for (int m=0; m<3; m++)
for (int o=0; o<3; o++)
tmp2[i][j][k][l] += polarInvert[m][i] * dQT[o][m][k][l] * polarInvert[j][o];
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
result[i][j][k][l] = 0.5 * (result[i][j][k][l] - tmp2[i][j][k][l]);
oldPolar -= polar;
if (oldPolar.frobenius_norm() < tol) {
break;
}
}
return result;
}
public:
/** \brief The type used for derivatives */
typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
/** \brief Constructor
* \param localFiniteElement A Lagrangian finite element that provides the interpolation points
* \param coefficients Values of the function at the Lagrange points
*/
LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& coefficients)
: localFiniteElement_(localFiniteElement),
coefficients_(coefficients)
{
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
/** \brief Rebind the FEFunction to another TargetSpace */
template<class U>
struct rebind
{
using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
};
/** \brief The number of Lagrange points */
unsigned int size() const
{
return localFiniteElement_.size();
}
/** \brief The type of the reference element */
Dune::GeometryType type() const
{
return localFiniteElement_.type();
}
/** \brief The scalar finite element used as the interpolation weights
*
* \note This method was added for InterpolationDerivatives, which needs it
* to construct a copy of a LocalGeodesicFEFunction with ADOL-C's adouble
* number type. This is not optimal, because the localFiniteElement
* really is an implementation detail of LocalGeodesicFEFunction and
* should not be needed just to copy an entire object. Other non-Euclidean
* interpolation rules may not have such a finite element at all.
* Therefore, this method may disappear again eventually.
*/
const LocalFiniteElement& localFiniteElement() const
{
return localFiniteElement_;
}
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
Rotation<field_type,3> result;
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
// Interpolate in R^{3x3}
FieldMatrix<field_type,3,3> interpolatedMatrix(0);
for (size_t i=0; i<coefficients_.size(); i++)
{
FieldMatrix<field_type,3,3> coefficientAsMatrix;
coefficients_[i].matrix(coefficientAsMatrix);
interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
}
// Project back onto SO(3)
result.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));
return result;
}
/** \brief Evaluate the derivative of the function */
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
// the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
// Actually compute the derivative
return evaluateDerivative(local,q);
}
/** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
* \param local Local coordinates in the reference element where to evaluate the derivative
* \param q Value of the local function at 'local'. If you provide something wrong here the result will be wrong, too!
*/
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const
{
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
std::vector<Dune::FieldMatrix<ctype,1,dim> > wDer;
localFiniteElement_.localBasis().evaluateJacobian(local,wDer);
// Compute matrix representations for all coefficients (we only have them in quaternion representation)
std::vector<Dune::FieldMatrix<field_type,3,3> > coefficientsAsMatrix(coefficients_.size());
for (size_t i=0; i<coefficients_.size(); i++)
coefficients_[i].matrix(coefficientsAsMatrix[i]);
// Interpolate in R^{3x3}
FieldMatrix<field_type,3,3> interpolatedMatrix(0);
for (size_t i=0; i<coefficients_.size(); i++)
interpolatedMatrix.axpy(w[i][0], coefficientsAsMatrix[i]);
Tensor3<RT,dim,3,3> derivative(0);
for (size_t dir=0; dir<dim; dir++)
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++)
for (size_t k=0; k<coefficients_.size(); k++)
derivative[dir][i][j] += wDer[k][0][dir] * coefficientsAsMatrix[k][i][j];
FieldMatrix<field_type,3,3> polarFactor;
auto derivativeOfProjection = this->derivativeOfProjection(interpolatedMatrix,polarFactor);
Tensor3<field_type,dim,3,3> intermediateResult(0);
for (size_t dir=0; dir<dim; dir++)
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++)
for (size_t k=0; k<3; k++)
for (size_t l=0; l<3; l++)
intermediateResult[dir][i][j] += derivativeOfProjection[i][j][k][l]*derivative[dir][k][l];
// One more application of the chain rule: we need to go from orthogonal matrices to quaternions
Tensor3<field_type,4,3,3> derivativeOfMatrixToQuaternion = Rotation<field_type,3>::derivativeOfMatrixToQuaternion(polarFactor);
DerivativeType result(0);
for (size_t dir0=0; dir0<4; dir0++)
for (size_t dir1=0; dir1<dim; dir1++)
for (size_t i=0; i<3; i++)
for (size_t j=0; j<3; j++)
result[dir0][dir1] += derivativeOfMatrixToQuaternion[dir0][i][j] * intermediateResult[dir1][i][j];
return result;
}
/** \brief Evaluate the value and the derivative of the interpolation function
*
* \return A std::pair containing the value and the first derivative of the interpolation function.
* If the interpolation is conforming then the first member of the pair will be a TargetSpace.
* Otherwise it will be a RealTuple.
*/
auto evaluateValueAndDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
// Construct the type of the result -- it depends on whether the interpolation
// is conforming or not.
using Value = std::conditional_t<conforming,TargetSpace,RealTuple<RT,embeddedDim> >;
std::pair<Value,DerivativeType> result;
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
// Interpolate in R^{3x3}
FieldMatrix<field_type,3,3> interpolatedMatrix(0);
for (size_t i=0; i<coefficients_.size(); i++)
{
FieldMatrix<field_type,3,3> coefficientAsMatrix;
coefficients_[i].matrix(coefficientAsMatrix);
interpolatedMatrix.axpy(w[i][0], coefficientAsMatrix);
}
// Project back onto SO(3)
static_assert(conforming, "Nonconforming interpolation into SO(3) is not implemented!");
result.first.set(Dune::GFE::PolarDecomposition()(interpolatedMatrix));
result.second = evaluateDerivative(local, result.first);
return result;
}
/** \brief Get the i'th base coefficient. */
const TargetSpace& coefficient(int i) const
{
return coefficients_[i];
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
* \todo We really only need the local basis
*/
const LocalFiniteElement& localFiniteElement_;
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
};
/** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for R^3 x SO(3)
*
* \tparam dim Dimension of the reference element
* \tparam ctype Type used for coordinates on the reference element
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
*/
template <int dim, class ctype, class LocalFiniteElement, class field_type>
class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
{
public:
using TargetSpace = ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >;
private:
typedef typename TargetSpace::ctype RT;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
public:
/** \brief The type used for derivatives */
typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
/** \brief Constructor
* \param localFiniteElement A Lagrangian finite element that provides the interpolation points
* \param coefficients Values of the function at the Lagrange points
*/
LocalProjectedFEFunction(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& coefficients)
: localFiniteElement_(localFiniteElement),
coefficients_(coefficients),
translationCoefficients_(coefficients.size())
{
assert(localFiniteElement.localBasis().size() == coefficients.size());
using namespace Dune::Indices;
for (size_t i=0; i<coefficients.size(); i++)
translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
for (size_t i=0; i<coefficients.size(); i++)
orientationCoefficients[i] = coefficients[i][_1];
orientationFunction_ = std::make_unique<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
}
/** \brief Rebind the FEFunction to another TargetSpace */
template<class U>
struct rebind
{
using other = LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,U>;
};
/** \brief The number of Lagrange points */
unsigned int size() const
{
return localFiniteElement_.size();
}
/** \brief The type of the reference element */
Dune::GeometryType type() const
{
return localFiniteElement_.type();
}
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
using namespace Dune::Indices;
TargetSpace result;
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
result[_0] = FieldVector<field_type,3>(0.0);
for (size_t i=0; i<w.size(); i++)
result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
result[_1] = orientationFunction_->evaluate(local);
return result;
}
/** \brief Evaluate the derivative of the function */
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
// the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
// Actually compute the derivative
return evaluateDerivative(local,q);
}
/** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
* \param local Local coordinates in the reference element where to evaluate the derivative
* \param q Value of the local function at 'local'. If you provide something wrong here the result will be wrong, too!
*/
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const
{
using namespace Dune::Indices;
DerivativeType result(0);
// get translation part
std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
for (size_t i=0; i<translationCoefficients_.size(); i++)
for (int j=0; j<3; j++)
result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
// get orientation part
Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_->evaluateDerivative(local,q[_1]);
for (int i=0; i<4; i++)
for (int j=0; j<dim; j++)
result[3+i][j] = qResult[i][j];
return result;
}
/** \brief Get the i'th base coefficient. */
const TargetSpace& coefficient(int i) const
{
return coefficients_[i];
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
* \todo We really only need the local basis
*/
const LocalFiniteElement& localFiniteElement_;
// The coefficients of this interpolation rule
std::vector<TargetSpace> coefficients_;
// The coefficients again. Yes, we store it twice: To evaluate the interpolation rule efficiently
// we need access to the coefficients for the various factor spaces separately.
std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
std::unique_ptr<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type, 3> > > orientationFunction_;
};
}
}
#endif
dune/gfe/localquickanddirtyfefunction.hh deleted 100644 → 0
View file @ 5db924c4
#ifndef DUNE_GFE_LOCALQUICKANDDIRTYFEFUNCTION_HH
#define DUNE_GFE_LOCALQUICKANDDIRTYFEFUNCTION_HH
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/gfe/linearalgebra.hh>
#include <dune/gfe/spaces/rotation.hh>
namespace Dune {
namespace GFE {
/** \brief Interpolate on a manifold, as fast as we can
*
* This class implements interpolation of values on a manifold in a 'quick-and-dirty' way.
* No particular interpolation rule is used consistently for all manifolds. Rather, it is
* used whatever is quickest and 'reasonable' for any given space. The reason to have this
* is to provide initial iterates for the iterative solvers used to solve the local
* geodesic-FE minimization problems.
*
* \note In particular, we currently do not implement the derivative of the interpolation,
* because it is not needed for producing initial iterates!
*
* \tparam dim Dimension of the reference element
* \tparam ctype Type used for coordinates on the reference element
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
* \tparam TargetSpace The manifold that the function takes its values in
*/
template <int dim, class ctype, class LocalFiniteElement, class TS>
class LocalQuickAndDirtyFEFunction
{
public:
using TargetSpace=TS;
private:
typedef typename TargetSpace::ctype RT;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
static const int spaceDim = TargetSpace::TangentVector::dimension;
public:
/** \brief The type used for derivatives */
typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
/** \brief Constructor
* \param localFiniteElement A Lagrangian finite element that provides the interpolation points
* \param coefficients Values of the function at the Lagrange points
*/
LocalQuickAndDirtyFEFunction(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& coefficients)
: localFiniteElement_(localFiniteElement),
coefficients_(coefficients)
{
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
/** \brief The number of Lagrange points */
unsigned int size() const
{
return localFiniteElement_.localBasis().size();
}
/** \brief The type of the reference element */
Dune::GeometryType type() const
{
return localFiniteElement_.type();
}
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Get the i'th base coefficient. */
TargetSpace coefficient(int i) const
{
return coefficients_[i];
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
* \todo We really only need the local basis
*/
const LocalFiniteElement& localFiniteElement_;
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
};
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
TargetSpace LocalQuickAndDirtyFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
// Special-casing for the Rotations. This is quite a challenge:
//
// Projection-based interpolation in the
// unit quaternions is not a good idea, here, because you may end up on the
// wrong sheet of the covering of SO(3) by the unit quaternions. The minimization
// solver for the weighted distance functional may then actually converge
// to the wrong local minimizer -- the interpolating functions wraps "the wrong way"
// around SO(3).
//
// Projection-based interpolation in SO(3) is not a good idea, because it is about as
// expensive as the geodesic interpolation itself. That's not good if all we want
// is an initial iterate.
//
// Averaging in R^{3x3} followed by QR decomposition is not a good idea either:
// The averaging can give you singular matrices quite quickly (try two matrices that
// differ by a rotation of pi/2). Then, the QR factorization of the identity may
// not be the identity (but differ in sign)! I tried that with the implementation
// from Numerical Recipes.
//
// What do we do instead? We start from the coefficient that produces the lowest
// value of the objective functional.
if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
{
AverageDistanceAssembler averageDistance(coefficients_,w);
auto minDistance = averageDistance.value(coefficients_[0]);
std::size_t minCoefficient = 0;
for (std::size_t i=1; i<coefficients_.size(); i++)
{
auto distance = averageDistance.value(coefficients_[i]);
if (distance < minDistance)
{
minDistance = distance;
minCoefficient = i;
}
}
return coefficients_[minCoefficient];
}
else
{
typename TargetSpace::CoordinateType c(0);
for (size_t i=0; i<coefficients_.size(); i++)
c.axpy(w[i][0], coefficients_[i].globalCoordinates());
return TargetSpace(c);
}
}
} // namespace GFE
} // namespace Dune
#endif
dune/gfe/localtangentfefunction.hh deleted 100644 → 0
View file @ 5db924c4
#ifndef DUNE_GFE_LOCALTANGENTFEFUNCTION_HH
#define DUNE_GFE_LOCALTANGENTFEFUNCTION_HH
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
namespace Dune {
namespace GFE {
/** \brief Interpolate on a tangent space
*
* \tparam dim Dimension of the reference element
* \tparam ctype Type used for coordinates on the reference element
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
* \tparam TargetSpace The manifold that the function takes its values in
*/
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
class LocalTangentFEFunction
{
typedef typename TargetSpace::ctype RT;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
static const int spaceDim = TargetSpace::TangentVector::dimension;
public:
/** \brief The type used for derivatives */
typedef Dune::FieldMatrix<RT, embeddedDim, dim> DerivativeType;
/** \brief Constructor
* \param localFiniteElement A Lagrangian finite element that provides the interpolation points
* \param coefficients Values of the function at the Lagrange points
*/
LocalTangentFEFunction(const LocalFiniteElement& localFiniteElement,
const std::vector<TargetSpace>& coefficients)
: localFiniteElement_(localFiniteElement),
coefficients_(coefficients)
{
assert(localFiniteElement_.localBasis().size() == coefficients_.size());
}
/** \brief The number of Lagrange points */
unsigned int size() const
{
return localFiniteElement_.localBasis().size();
}
/** \brief The type of the reference element */
Dune::GeometryType type() const
{
return localFiniteElement_.type();
}
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Evaluate the derivative of the function */
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const;
/** \brief Evaluate the derivative of the function, if you happen to know the function value (much faster!)
* \param local Local coordinates in the reference element where to evaluate the derivative
* \param q Value of the local gfe function at 'local'. If you provide something wrong here the result will be wrong, too!
*/
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const;
/** \brief Get the i'th base coefficient. */
TargetSpace coefficient(int i) const
{
return coefficients_[i];
}
private:
/** \brief The scalar local finite element, which provides the weighting factors
* \todo We really only need the local basis
*/
const LocalFiniteElement& localFiniteElement_;
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
};
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
TargetSpace LocalTangentFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
// This is the base point of the tangent space that we are interpolating on
// HACK: This point is reasonable only for spheres
typename TargetSpace::CoordinateType baseEmb(0);
baseEmb[0] = 1.0;
TargetSpace base(baseEmb);
// Map all points onto the tangent space at 'base'
std::vector<typename TargetSpace::EmbeddedTangentVector> tangentCoefficients(coefficients_.size());
for (size_t i=0; i<coefficients_.size(); i++)
tangentCoefficients[i] = TargetSpace::log(base,coefficients_[i]);
// Interpolate on the tangent space
typename TargetSpace::EmbeddedTangentVector c(0);
for (size_t i=0; i<tangentCoefficients.size(); i++)
c.axpy(w[i][0], tangentCoefficients[i]);
return TargetSpace::exp(base,c);
}
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
typename LocalTangentFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
LocalTangentFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
// the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
// Actually compute the derivative
return evaluateDerivative(local,q);
}
template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
typename LocalTangentFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
LocalTangentFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace& q) const
{
DUNE_THROW(NotImplemented, "evaluateDerivative not implemented yet!");
}
}
}
#endif
dune/gfe/maxnormtrustregion.hh
View file @ 48b9358e
......@@ -5,95 +5,101 @@
#include <dune/solvers/common/boxconstraint.hh>
template <int blocksize, class field_type=double>
class MaxNormTrustRegion
namespace Dune::GFE
{
public:
MaxNormTrustRegion(size_t size, field_type initialRadius)
: obstacles_(size)
template <int blocksize, class field_type=double>
class MaxNormTrustRegion
{
set(initialRadius);
}
public:
MaxNormTrustRegion(size_t size, field_type initialRadius)
: obstacles_(size)
{
set(initialRadius);
}
void set(field_type radius) {
void set(field_type radius) {
radius_ = radius;
radius_ = radius;
for (size_t i=0; i<obstacles_.size(); i++) {
for (size_t i=0; i<obstacles_.size(); i++) {
for (int k=0; k<blocksize; k++) {
for (int k=0; k<blocksize; k++) {
obstacles_[i].lower(k) = -radius;
obstacles_[i].upper(k) = radius;
obstacles_[i].lower(k) = -radius;
obstacles_[i].upper(k) = radius;
}
}
}
}
/** \brief Set trust region radius with a separate scaling for each vector block component
*/
void set(field_type radius, const Dune::FieldVector<field_type, blocksize>& scaling) {
/** \brief Set trust region radius with a separate scaling for each vector block component
*/
void set(field_type radius, const Dune::FieldVector<field_type, blocksize>& scaling) {
radius_ = radius;
radius_ = radius;
for (size_t i=0; i<obstacles_.size(); i++) {
for (size_t i=0; i<obstacles_.size(); i++) {
for (int k=0; k<blocksize; k++) {
for (int k=0; k<blocksize; k++) {
obstacles_[i].lower(k) = -radius*scaling[k];
obstacles_[i].upper(k) = radius*scaling[k];
obstacles_[i].lower(k) = -radius*scaling[k];
obstacles_[i].upper(k) = radius*scaling[k];
}
}
}
}
/** \brief Return true if the given vector is not contained in the trust region */
bool violates(const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& v) const
{
assert(v.size() == obstacles_.size());
for (size_t i=0; i<v.size(); i++)
for (size_t j=0; j<blocksize; j++)
if (v[i][j] < obstacles_[i].lower(j) or v[i][j] > obstacles_[i].upper(j))
return true;
/** \brief Return true if the given vector is not contained in the trust region */
bool violates(const Dune::BlockVector<Dune::FieldVector<double,blocksize> >& v) const
{
assert(v.size() == obstacles_.size());
for (size_t i=0; i<v.size(); i++)
for (size_t j=0; j<blocksize; j++)
if (v[i][j] < obstacles_[i].lower(j) or v[i][j] > obstacles_[i].upper(j))
return true;
return false;
}
return false;
}
field_type radius() const {
return radius_;
}
field_type radius() const {
return radius_;
}
void scale(field_type factor) {
void scale(field_type factor) {
radius_ *= factor;
radius_ *= factor;
for (size_t i=0; i<obstacles_.size(); i++) {
for (size_t i=0; i<obstacles_.size(); i++) {
for (int k=0; k<blocksize; k++) {
for (int k=0; k<blocksize; k++) {
obstacles_[i].lower(k) *= factor;
obstacles_[i].upper(k) *= factor;
obstacles_[i].lower(k) *= factor;
obstacles_[i].upper(k) *= factor;
}
}
}
}
const std::vector<BoxConstraint<field_type,blocksize> >& obstacles() const {
return obstacles_;
}
const std::vector<BoxConstraint<field_type,blocksize> >& obstacles() const {
return obstacles_;
}
private:
private:
std::vector<BoxConstraint<field_type,blocksize> > obstacles_;
std::vector<BoxConstraint<field_type,blocksize> > obstacles_;
field_type radius_;
field_type radius_;
};
};
} // namespace Dune::GFE
#endif
dune/gfe/mixedriemanniantrsolver.cc
View file @ 48b9358e
......@@ -25,7 +25,7 @@ template <class GridType,
class Basis,
class Basis0, class TargetSpace0,
class Basis1, class TargetSpace1>
void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::
void Dune::GFE::MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::
setup(const GridType& grid,
const MixedGFEAssembler<Basis, TargetSpace>* assembler,
const Basis0& tmpBasis0,
......@@ -262,7 +262,7 @@ template <class GridType,
class Basis,
class Basis0, class TargetSpace0,
class Basis1, class TargetSpace1>
void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::solve()
void Dune::GFE::MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,TargetSpace1>::solve()
{
int argc = 0;
char** argv;
......
dune/gfe/mixedriemanniantrsolver.hh
View file @ 48b9358e
......@@ -22,168 +22,174 @@
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#include <dune/gfe/spaces/productmanifold.hh>
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
template <class GridType,
class Basis,
class Basis0, class TargetSpace0,
class Basis1, class TargetSpace1>
class MixedRiemannianTrustRegionSolver
: public NumProc
{
const static int blocksize0 = TargetSpace0::TangentVector::dimension;
const static int blocksize1 = TargetSpace1::TangentVector::dimension;
const static int gridDim = GridType::dimension;
using TargetSpace = Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1>;
// Centralize the field type here
typedef double field_type;
// Some types that I need
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize0, blocksize0> > MatrixType00;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize0, blocksize1> > MatrixType01;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize1, blocksize0> > MatrixType10;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize1, blocksize1> > MatrixType11;
typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixType00,MatrixType01>,
Dune::MultiTypeBlockVector<MatrixType10,MatrixType11> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize0> > CorrectionType0;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize1> > CorrectionType1;
typedef Dune::MultiTypeBlockVector<CorrectionType0, CorrectionType1> CorrectionType;
typedef Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > SolutionType;
/** \brief Records information about the last run of the RiemannianTrustRegionSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
{
std::size_t finalIteration;
field_type finalEnergy;
};
public:
MixedRiemannianTrustRegionSolver()
: NumProc(NumProc::FULL),
namespace Dune::GFE
{
h1SemiNorm0_(nullptr), h1SemiNorm1_(nullptr)
{
std::fill(std::get<0>(scaling_).begin(), std::get<0>(scaling_).end(), 1.0);
std::fill(std::get<1>(scaling_).begin(), std::get<1>(scaling_).end(), 1.0);
}
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const MixedGFEAssembler<Basis, TargetSpace>* assembler,
const Basis0& basis0,
const Basis1& basis1,
const SolutionType& x,
const Dune::BitSetVector<blocksize0>& dirichletNodes0,
const Dune::BitSetVector<blocksize1>& dirichletNodes1,
double tolerance,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
int multigridIterations,
double mgTolerance,
int mu,
int nu1,
int nu2,
int baseIterations,
double baseTolerance,
bool instrumented);
void setScaling(const Dune::FieldVector<double,blocksize0+blocksize1>& scaling)
/** \brief Riemannian trust-region solver for geodesic finite-element problems */
template <class GridType,
class Basis,
class Basis0, class TargetSpace0,
class Basis1, class TargetSpace1>
class MixedRiemannianTrustRegionSolver
: public NumProc
{
for (int i=0; i<blocksize0; i++)
std::get<0>(scaling_)[i] = scaling[i];
for (int i=0; i<blocksize1; i++)
std::get<1>(scaling_)[i] = scaling[i+blocksize0];
}
const static int blocksize0 = TargetSpace0::TangentVector::dimension;
const static int blocksize1 = TargetSpace1::TangentVector::dimension;
const static int gridDim = GridType::dimension;
using TargetSpace = Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1>;
// Centralize the field type here
typedef double field_type;
// Some types that I need
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize0, blocksize0> > MatrixType00;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize0, blocksize1> > MatrixType01;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize1, blocksize0> > MatrixType10;
typedef Dune::BCRSMatrix<Dune::FieldMatrix<field_type, blocksize1, blocksize1> > MatrixType11;
typedef Dune::MultiTypeBlockMatrix<Dune::MultiTypeBlockVector<MatrixType00,MatrixType01>,
Dune::MultiTypeBlockVector<MatrixType10,MatrixType11> > MatrixType;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize0> > CorrectionType0;
typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize1> > CorrectionType1;
typedef Dune::MultiTypeBlockVector<CorrectionType0, CorrectionType1> CorrectionType;
typedef Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > SolutionType;
/** \brief Records information about the last run of the RiemannianTrustRegionSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
{
std::size_t finalIteration;
field_type finalEnergy;
};
public:
MixedRiemannianTrustRegionSolver()
: NumProc(NumProc::FULL),
h1SemiNorm0_(nullptr), h1SemiNorm1_(nullptr)
{
std::fill(std::get<0>(scaling_).begin(), std::get<0>(scaling_).end(), 1.0);
std::fill(std::get<1>(scaling_).begin(), std::get<1>(scaling_).end(), 1.0);
}
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const MixedGFEAssembler<Basis, TargetSpace>* assembler,
const Basis0& basis0,
const Basis1& basis1,
const SolutionType& x,
const Dune::BitSetVector<blocksize0>& dirichletNodes0,
const Dune::BitSetVector<blocksize1>& dirichletNodes1,
double tolerance,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
int multigridIterations,
double mgTolerance,
int mu,
int nu1,
int nu2,
int baseIterations,
double baseTolerance,
bool instrumented);
void setScaling(const Dune::FieldVector<double,blocksize0+blocksize1>& scaling)
{
for (int i=0; i<blocksize0; i++)
std::get<0>(scaling_)[i] = scaling[i];
for (int i=0; i<blocksize1; i++)
std::get<1>(scaling_)[i] = scaling[i+blocksize0];
}
#if 0
void setIgnoreNodes(const Dune::BitSetVector<blocksize0>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
assert(loopSolver);
loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
}
void setIgnoreNodes(const Dune::BitSetVector<blocksize0>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
Dune::shared_ptr<LoopSolver<CorrectionType> > loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
assert(loopSolver);
loopSolver->iterationStep_->ignoreNodes_ = ignoreNodes_;
}
#endif
void solve();
void solve();
void setInitialIterate(const SolutionType& x)
{
x_ = x;
}
void setInitialIterate(const SolutionType& x)
{
x_ = x;
}
SolutionType getSol() const
{
return x_;
}
SolutionType getSol() const
{
return x_;
}
const Statistics& getStatistics() const {return statistics_;}
const Statistics& getStatistics() const {return statistics_;}
protected:
protected:
#if 0
std::unique_ptr<GUIndex> guIndex_;
std::unique_ptr<GUIndex> guIndex_;
#endif
/** \brief The grid */
const GridType* grid_;
/** \brief The grid */
const GridType* grid_;
/** \brief The solution vectors */
SolutionType x_;
/** \brief The solution vectors */
SolutionType x_;
/** \brief The initial trust-region radius in the maximum-norm */
double initialTrustRegionRadius_;
/** \brief The initial trust-region radius in the maximum-norm */
double initialTrustRegionRadius_;
/** \brief Trust-region norm scaling */
std::tuple<Dune::FieldVector<double,blocksize0>, Dune::FieldVector<double,blocksize1> > scaling_;
/** \brief Trust-region norm scaling */
std::tuple<Dune::FieldVector<double,blocksize0>, Dune::FieldVector<double,blocksize1> > scaling_;
/** \brief Maximum number of trust-region steps */
int maxTrustRegionSteps_;
/** \brief Maximum number of trust-region steps */
int maxTrustRegionSteps_;
/** \brief Maximum number of multigrid iterations */
int innerIterations_;
/** \brief Maximum number of multigrid iterations */
int innerIterations_;
/** \brief Error tolerance of the multigrid QP solver */
double innerTolerance_;
/** \brief Error tolerance of the multigrid QP solver */
double innerTolerance_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
const MixedGFEAssembler<Basis, TargetSpace>* assembler_;
/** \brief The assembler for the material law */
const MixedGFEAssembler<Basis, TargetSpace>* assembler_;
/** \brief TEMPORARY: The two separate matrices */
std::unique_ptr<Basis0> basis0_;
std::unique_ptr<Basis1> basis1_;
/** \brief TEMPORARY: The two separate matrices */
std::unique_ptr<Basis0> basis0_;
std::unique_ptr<Basis1> basis1_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;
MonotoneMGStep<MatrixType00, CorrectionType0>* mmgStep0;
MonotoneMGStep<MatrixType11, CorrectionType1>* mmgStep1;
MonotoneMGStep<MatrixType00, CorrectionType0>* mmgStep0;
MonotoneMGStep<MatrixType11, CorrectionType1>* mmgStep1;
double tolerance_;
double tolerance_;
/** \brief Contains 'true' everywhere -- the trust-region is bounded */
Dune::BitSetVector<blocksize0> hasObstacle0_;
Dune::BitSetVector<blocksize1> hasObstacle1_;
/** \brief Contains 'true' everywhere -- the trust-region is bounded */
Dune::BitSetVector<blocksize0> hasObstacle0_;
Dune::BitSetVector<blocksize1> hasObstacle1_;
/** \brief The Dirichlet nodes */
const Dune::BitSetVector<blocksize0>* ignoreNodes0_;
const Dune::BitSetVector<blocksize1>* ignoreNodes1_;
/** \brief The Dirichlet nodes */
const Dune::BitSetVector<blocksize0>* ignoreNodes0_;
const Dune::BitSetVector<blocksize1>* ignoreNodes1_;
/** \brief The norm used to measure multigrid convergence */
H1SemiNorm<CorrectionType0>* h1SemiNorm0_;
H1SemiNorm<CorrectionType1>* h1SemiNorm1_;
/** \brief The norm used to measure multigrid convergence */
H1SemiNorm<CorrectionType0>* h1SemiNorm0_;
H1SemiNorm<CorrectionType1>* h1SemiNorm1_;
/** \brief Store information about solver runs for testing */
Statistics statistics_;
};
/** \brief Store information about solver runs for testing */
Statistics statistics_;
};
} // namespace Dune::GFE
#include "mixedriemanniantrsolver.cc"
......
dune/gfe/neumannenergy.hh
View file @ 48b9358e
......@@ -5,12 +5,17 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/elasticity/assemblers/localenergy.hh>
#include <dune/gfe/assemblers/localenergy.hh>
#include <dune/gfe/spaces/productmanifold.hh>
namespace Dune::GFE
{
namespace Dune::GFE {
/** \brief Integrate a density over a part of the domain boundary
*
* This is typically used to implement Neumann boundary conditions. Hence the name.
* This is typically used to implement Neumann boundary conditions for Cosserat materials.
* Essentially it works only for that: It is implicitly assumed that the target space
* is a product, and that the first factor is a RealTuple.
*/
template<class Basis, class ... TargetSpaces>
class NeumannEnergy
......@@ -24,6 +29,10 @@ namespace Dune::GFE {
using RT = typename Dune::GFE::LocalEnergy<Basis,TargetSpace>::RT;
constexpr static int dim = GridView::dimension;
constexpr static int dimworld = GridView::dimensionworld;
// TODO: Remove the hard-coded first factor space!
using WorldVector = typename std::tuple_element_t<0, std::tuple<TargetSpaces...> >::EmbeddedTangentVector;
public:
......@@ -31,7 +40,7 @@ namespace Dune::GFE {
* \param parameters The material parameters
*/
NeumannEnergy(const std::shared_ptr<BoundaryPatch<GridView> >& neumannBoundary,
std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<DT,dim>)> neumannFunction)
std::function<WorldVector(Dune::FieldVector<DT,dimworld>)> neumannFunction)
: neumannBoundary_(neumannBoundary),
neumannFunction_(neumannFunction)
{}
......@@ -77,9 +86,9 @@ namespace Dune::GFE {
std::vector<Dune::FieldVector<DT,1> > shapeFunctionValues;
localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
Dune::FieldVector<RT,dim> value(0);
WorldVector value(0);
for (size_t i=0; i<localFiniteElement.size(); i++)
for (int j=0; j<dim; j++)
for (int j=0; j<WorldVector::dimension; j++)
value[j] += shapeFunctionValues[i] * localSolution[i][j];
// Value of the Neumann data at the current position
......@@ -101,7 +110,7 @@ namespace Dune::GFE {
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_;
std::function<WorldVector(Dune::FieldVector<DT,dimworld>)> neumannFunction_;
};
} // namespace Dune::GFE
......
dune/gfe/parallel/globalmapper.hh
View file @ 48b9358e
......@@ -7,7 +7,8 @@
#include <dune/parmg/parallel/datahandle.hh>
#endif
namespace Dune {
namespace Dune::GFE
{
template <class Basis>
class GlobalMapper
......@@ -63,5 +64,6 @@ namespace Dune {
std::size_t size_;
};
}
} // namespace Dune::GFE
#endif // DUNE_GFE_PARALLEL_GLOBALMAPPER_HH
dune/gfe/parallel/globalp1mapper.hh
View file @ 48b9358e
......@@ -6,7 +6,8 @@
#include <dune/gfe/parallel/globalmapper.hh>
namespace Dune {
namespace Dune::GFE
{
template <class Basis>
class GlobalP1Mapper : public GlobalMapper<Basis>
......@@ -79,5 +80,6 @@ namespace Dune {
#endif
};
}
} // namespace Dune::GFE
#endif /* DUNE_GFE_PARALLEL_GLOBALP1MAPPER_HH */
dune/gfe/parallel/globalp2mapper.hh
View file @ 48b9358e
......@@ -7,10 +7,11 @@
#include <dune/gfe/parallel/globalmapper.hh>
namespace Dune {
namespace Dune::GFE
{
template <class Basis>
class GlobalP2Mapper : public Dune::GlobalMapper<Basis>
class GlobalP2Mapper : public GlobalMapper<Basis>
{
using GridView = typename Basis::GridView;
using P2BasisMapper = Functions::LagrangeBasis<GridView,2>;
......@@ -119,5 +120,7 @@ namespace Dune {
P2BasisMapper p2Mapper_;
};
}
} // namespace Dune::GFE
#endif // DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
dune/gfe/parallel/mapperfactory.hh
View file @ 48b9358e
......@@ -19,7 +19,7 @@ namespace Dune::GFE
template <typename GridView>
struct MapperFactory<Functions::LagrangeBasis<GridView,1> >
{
typedef Dune::GlobalP1Mapper<Functions::LagrangeBasis<GridView,1> > GlobalMapper;
typedef GlobalP1Mapper<Functions::LagrangeBasis<GridView,1> > GlobalMapper;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
......@@ -30,7 +30,7 @@ namespace Dune::GFE
template <typename GridView>
struct MapperFactory<Dune::Functions::LagrangeBasis<GridView,2> >
{
typedef Dune::GlobalP2Mapper<Functions::LagrangeBasis<GridView,2> > GlobalMapper;
typedef GlobalP2Mapper<Functions::LagrangeBasis<GridView,2> > GlobalMapper;
typedef P2BasisMapper<GridView> LocalMapper;
static LocalMapper createLocalMapper(const GridView& gridView)
{
......
dune/gfe/parallel/matrixcommunicator.hh
View file @ 48b9358e
......@@ -9,169 +9,174 @@
#include <dune/gfe/parallel/mpifunctions.hh>
template<typename RowGlobalMapper, typename GridView1, typename GridView2, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColumnGlobalMapper=RowGlobalMapper>
class MatrixCommunicator {
namespace Dune::GFE
{
struct TransferMatrixTuple {
typedef typename MatrixType::block_type EntryType;
template<typename RowGlobalMapper, typename GridView1, typename GridView2, typename MatrixType, typename LocalMapper1, typename LocalMapper2, typename ColumnGlobalMapper=RowGlobalMapper>
class MatrixCommunicator {
size_t row, col;
EntryType entry;
struct TransferMatrixTuple {
typedef typename MatrixType::block_type EntryType;
TransferMatrixTuple() {}
TransferMatrixTuple(const size_t& r, const size_t& c, const EntryType& e) : row(r), col(c), entry(e) {}
};
size_t row, col;
EntryType entry;
void transferMatrix(const MatrixType& localMatrix) {
// Create vector for transfer data
std::vector<TransferMatrixTuple> localMatrixEntries;
TransferMatrixTuple() {}
TransferMatrixTuple(const size_t& r, const size_t& c, const EntryType& e) : row(r), col(c), entry(e) {}
};
// Convert local matrix to serializable array
typedef typename MatrixType::row_type::ConstIterator ColumnIterator;
void transferMatrix(const MatrixType& localMatrix) {
// Create vector for transfer data
std::vector<TransferMatrixTuple> localMatrixEntries;
for (typename MatrixType::ConstIterator rIt = localMatrix.begin(); rIt != localMatrix.end(); ++rIt)
for (ColumnIterator cIt = rIt->begin(); cIt != rIt->end(); ++cIt) {
const int i = rIt.index();
const int j = cIt.index();
// Convert local matrix to serializable array
typedef typename MatrixType::row_type::ConstIterator ColumnIterator;
localMatrixEntries.push_back(TransferMatrixTuple(localToGlobal1_[i], localToGlobal2_[j], *cIt));
}
for (typename MatrixType::ConstIterator rIt = localMatrix.begin(); rIt != localMatrix.end(); ++rIt)
for (ColumnIterator cIt = rIt->begin(); cIt != rIt->end(); ++cIt) {
const int i = rIt.index();
const int j = cIt.index();
// Get number of matrix entries on each process
std::vector<int> localMatrixEntriesSizes(MPIFunctions::shareSizes(communicator_, localMatrixEntries.size()));
localMatrixEntries.push_back(TransferMatrixTuple(localToGlobal1_[i], localToGlobal2_[j], *cIt));
}
// Get matrix entries from every process
globalMatrixEntries = MPIFunctions::gatherv(communicator_, localMatrixEntries, localMatrixEntriesSizes, root_rank);
}
// Get number of matrix entries on each process
std::vector<int> localMatrixEntriesSizes(MPIFunctions::shareSizes(communicator_, localMatrixEntries.size()));
public:
MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const GridView1& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
: rowGlobalMapper_(rowGlobalMapper),
columnGlobalMapper_(rowGlobalMapper),
localMapper1_(localMapper1),
localMapper2_(localMapper2),
communicator_(gridView.comm()),
root_rank(root)
{
setLocalToGlobalRows(gridView);
setLocalToGlobalColumns(gridView);
}
// Get matrix entries from every process
globalMatrixEntries = MPIFunctions::gatherv(communicator_, localMatrixEntries, localMatrixEntriesSizes, root_rank);
}
MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const ColumnGlobalMapper& columnGlobalMapper,
const GridView1& gridView1, const GridView2& gridView2,
const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
: rowGlobalMapper_(rowGlobalMapper),
columnGlobalMapper_(columnGlobalMapper),
localMapper1_(localMapper1),
localMapper2_(localMapper2),
communicator_(gridView1.comm()),
root_rank(root)
{
setLocalToGlobalRows(gridView1);
setLocalToGlobalColumns(gridView2);
}
public:
MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const GridView1& gridView, const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
: rowGlobalMapper_(rowGlobalMapper),
columnGlobalMapper_(rowGlobalMapper),
localMapper1_(localMapper1),
localMapper2_(localMapper2),
communicator_(gridView.comm()),
root_rank(root)
{
setLocalToGlobalRows(gridView);
setLocalToGlobalColumns(gridView);
}
MatrixType reduceAdd(const MatrixType& local)
{
transferMatrix(local);
MatrixCommunicator(const RowGlobalMapper& rowGlobalMapper, const ColumnGlobalMapper& columnGlobalMapper,
const GridView1& gridView1, const GridView2& gridView2,
const LocalMapper1& localMapper1, const LocalMapper2& localMapper2, const int& root)
: rowGlobalMapper_(rowGlobalMapper),
columnGlobalMapper_(columnGlobalMapper),
localMapper1_(localMapper1),
localMapper2_(localMapper2),
communicator_(gridView1.comm()),
root_rank(root)
{
setLocalToGlobalRows(gridView1);
setLocalToGlobalColumns(gridView2);
}
MatrixType globalMatrix;
MatrixType reduceAdd(const MatrixType& local)
{
transferMatrix(local);
// Create occupation pattern in matrix
Dune::MatrixIndexSet occupationPattern;
MatrixType globalMatrix;
occupationPattern.resize(rowGlobalMapper_.size(), columnGlobalMapper_.size());
// Create occupation pattern in matrix
Dune::MatrixIndexSet occupationPattern;
for (size_t k = 0; k < globalMatrixEntries.size(); ++k)
occupationPattern.add(globalMatrixEntries[k].row, globalMatrixEntries[k].col);
occupationPattern.resize(rowGlobalMapper_.size(), columnGlobalMapper_.size());
occupationPattern.exportIdx(globalMatrix);
for (size_t k = 0; k < globalMatrixEntries.size(); ++k)
occupationPattern.add(globalMatrixEntries[k].row, globalMatrixEntries[k].col);
// Initialize matrix to zero
globalMatrix = 0;
occupationPattern.exportIdx(globalMatrix);
// Move entries to matrix
for(size_t k = 0; k < globalMatrixEntries.size(); ++k)
globalMatrix[globalMatrixEntries[k].row][globalMatrixEntries[k].col] += globalMatrixEntries[k].entry;
// Initialize matrix to zero
globalMatrix = 0;
return globalMatrix;
}
// Move entries to matrix
for(size_t k = 0; k < globalMatrixEntries.size(); ++k)
globalMatrix[globalMatrixEntries[k].row][globalMatrixEntries[k].col] += globalMatrixEntries[k].entry;
MatrixType reduceCopy(const MatrixType& local)
{
transferMatrix(local);
return globalMatrix;
}
MatrixType globalMatrix;
MatrixType reduceCopy(const MatrixType& local)
{
transferMatrix(local);
// Create occupation pattern in matrix
Dune::MatrixIndexSet occupationPattern;
MatrixType globalMatrix;
occupationPattern.resize(rowGlobalMapper_.size(), columnGlobalMapper_.size());
// Create occupation pattern in matrix
Dune::MatrixIndexSet occupationPattern;
for (size_t k = 0; k < globalMatrixEntries.size(); ++k)
occupationPattern.add(globalMatrixEntries[k].row, globalMatrixEntries[k].col);
occupationPattern.resize(rowGlobalMapper_.size(), columnGlobalMapper_.size());
occupationPattern.exportIdx(globalMatrix);
for (size_t k = 0; k < globalMatrixEntries.size(); ++k)
occupationPattern.add(globalMatrixEntries[k].row, globalMatrixEntries[k].col);
// Move entries to matrix
for(size_t k = 0; k < globalMatrixEntries.size(); ++k)
globalMatrix[globalMatrixEntries[k].row][globalMatrixEntries[k].col] = globalMatrixEntries[k].entry;
occupationPattern.exportIdx(globalMatrix);
// Move entries to matrix
for(size_t k = 0; k < globalMatrixEntries.size(); ++k)
globalMatrix[globalMatrixEntries[k].row][globalMatrixEntries[k].col] = globalMatrixEntries[k].entry;
return globalMatrix;
}
private:
return globalMatrix;
}
void setLocalToGlobalRows(const GridView1& gridView)
{
localToGlobal1_.resize(localMapper1_.size());
private:
for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
for (int codim = 0; codim <= GridView1::dimension; codim++)
for (size_t i=0; i<it->subEntities(codim); i++)
{
typename LocalMapper1::Index localIdx;
typename RowGlobalMapper::Index globalIdx;
if (localMapper1_.contains(*it,i,codim,localIdx)
and rowGlobalMapper_.contains(*it,i,codim,globalIdx))
localToGlobal1_[localIdx] = globalIdx;
}
void setLocalToGlobalRows(const GridView1& gridView)
{
localToGlobal1_.resize(localMapper1_.size());
}
void setLocalToGlobalColumns(const GridView2& gridView)
{
localToGlobal2_.resize(localMapper2_.size());
for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
for (int codim = 0; codim <= GridView2::dimension; codim++)
for (size_t i=0; i<it->subEntities(codim); i++)
{
typename LocalMapper2::Index localIdx;
typename ColumnGlobalMapper::Index globalIdx;
if (localMapper2_.contains(*it,i,codim,localIdx)
and columnGlobalMapper_.contains(*it,i,codim,globalIdx))
localToGlobal2_[localIdx] = globalIdx;
}
for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
for (int codim = 0; codim <= GridView1::dimension; codim++)
for (size_t i=0; i<it->subEntities(codim); i++)
{
typename LocalMapper1::Index localIdx;
typename RowGlobalMapper::Index globalIdx;
if (localMapper1_.contains(*it,i,codim,localIdx)
and rowGlobalMapper_.contains(*it,i,codim,globalIdx))
localToGlobal1_[localIdx] = globalIdx;
}
}
}
void setLocalToGlobalColumns(const GridView2& gridView)
{
localToGlobal2_.resize(localMapper2_.size());
// Mappers for the global numbering
const RowGlobalMapper& rowGlobalMapper_;
const ColumnGlobalMapper& columnGlobalMapper_;
for (auto it = gridView.template begin<0>(); it != gridView.template end<0>(); ++it)
for (int codim = 0; codim <= GridView2::dimension; codim++)
for (size_t i=0; i<it->subEntities(codim); i++)
{
typename LocalMapper2::Index localIdx;
typename ColumnGlobalMapper::Index globalIdx;
if (localMapper2_.contains(*it,i,codim,localIdx)
and columnGlobalMapper_.contains(*it,i,codim,globalIdx))
localToGlobal2_[localIdx] = globalIdx;
}
// Mappers for the local numbering
const LocalMapper1& localMapper1_;
const LocalMapper2& localMapper2_;
}
const typename GridView1::Communication& communicator_;
int root_rank;
// Mappers for the global numbering
const RowGlobalMapper& rowGlobalMapper_;
const ColumnGlobalMapper& columnGlobalMapper_;
std::vector<typename RowGlobalMapper::Index> localToGlobal1_;
std::vector<typename ColumnGlobalMapper::Index> localToGlobal2_;
// Mappers for the local numbering
const LocalMapper1& localMapper1_;
const LocalMapper2& localMapper2_;
const typename GridView1::Communication& communicator_;
int root_rank;
std::vector<typename RowGlobalMapper::Index> localToGlobal1_;
std::vector<typename ColumnGlobalMapper::Index> localToGlobal2_;
std::vector<TransferMatrixTuple> globalMatrixEntries;
};
std::vector<TransferMatrixTuple> globalMatrixEntries;
};
} // namespace Dune::GFE
#endif
dune/gfe/parallel/mpifunctions.hh
View file @ 48b9358e
......@@ -4,67 +4,72 @@
#include <numeric>
#include <vector>
struct MPIFunctions {
static std::vector<int> offsetsFromSizes(const std::vector<int>& sizes) {
std::vector<int> offsets(sizes.size());
for (size_t k = 0; k < offsets.size(); ++k)
offsets[k] = std::accumulate(sizes.begin(), sizes.begin() + k, 0);
namespace Dune::GFE
{
return offsets;
}
struct MPIFunctions {
static std::vector<int> offsetsFromSizes(const std::vector<int>& sizes) {
std::vector<int> offsets(sizes.size());
template<typename Communicator>
static std::vector<int> shareSizes(const Communicator& communicator, const int& shareRef) {
std::vector<int> sizesVector(communicator.size());
for (size_t k = 0; k < offsets.size(); ++k)
offsets[k] = std::accumulate(sizes.begin(), sizes.begin() + k, 0);
int share = shareRef;
communicator.template allgather<int>(&share, 1, sizesVector.data());
return offsets;
}
template<typename Communicator>
static std::vector<int> shareSizes(const Communicator& communicator, const int& shareRef) {
std::vector<int> sizesVector(communicator.size());
return sizesVector;
}
int share = shareRef;
communicator.template allgather<int>(&share, 1, sizesVector.data());
template<typename Communicator, typename T>
static void scatterv(const Communicator& communicator, std::vector<T>& localVec, std::vector<T>& globalVec, std::vector<int>& sizes, int root_rank) {
int mysize = localVec.size();
std::vector<int> offsets(offsetsFromSizes(sizes));
return sizesVector;
}
communicator.template scatterv(globalVec.data(), sizes.data(), offsets.data(), localVec.data(), mysize, root_rank);
}
template<typename Communicator, typename T>
static void scatterv(const Communicator& communicator, std::vector<T>& localVec, std::vector<T>& globalVec, std::vector<int>& sizes, int root_rank) {
int mysize = localVec.size();
template<typename Communicator, typename T>
static std::vector<T> gatherv(const Communicator& communicator, std::vector<T>& localVec, std::vector<int>& sizes, int root_rank) {
int mysize = localVec.size();
std::vector<int> offsets(offsetsFromSizes(sizes));
std::vector<T> globalVec;
communicator.scatterv(globalVec.data(), sizes.data(), offsets.data(), localVec.data(), mysize, root_rank);
}
if (communicator.rank() == root_rank)
globalVec.resize(std::accumulate(sizes.begin(), sizes.end(), 0));
template<typename Communicator, typename T>
static std::vector<T> gatherv(const Communicator& communicator, std::vector<T>& localVec, std::vector<int>& sizes, int root_rank) {
int mysize = localVec.size();
std::vector<int> offsets(offsetsFromSizes(sizes));
std::vector<T> globalVec;
communicator.template gatherv(localVec.data(), mysize, globalVec.data(), sizes.data(), offsets.data(), root_rank);
if (communicator.rank() == root_rank)
globalVec.resize(std::accumulate(sizes.begin(), sizes.end(), 0));
std::vector<int> offsets(offsetsFromSizes(sizes));
return globalVec;
}
communicator.gatherv(localVec.data(), mysize, globalVec.data(), sizes.data(), offsets.data(), root_rank);
template<typename Communicator, typename T>
static std::vector<T> allgatherv(const Communicator& communicator, std::vector<T>& localVec, std::vector<int>& sizes) {
int mysize = localVec.size();
std::vector<T> globalVec(std::accumulate(sizes.begin(), sizes.end(), 0));
return globalVec;
}
std::vector<int> offsets(offsetsFromSizes(sizes));
template<typename Communicator, typename T>
static std::vector<T> allgatherv(const Communicator& communicator, std::vector<T>& localVec, std::vector<int>& sizes) {
int mysize = localVec.size();
communicator.template allgatherv(localVec.data(), mysize, globalVec.data(), sizes.data(), offsets.data());
std::vector<T> globalVec(std::accumulate(sizes.begin(), sizes.end(), 0));
std::vector<int> offsets(offsetsFromSizes(sizes));
return globalVec;
}
};
communicator.allgatherv(localVec.data(), mysize, globalVec.data(), sizes.data(), offsets.data());
return globalVec;
}
};
} // namespace Dune::GFE
#endif
dune/gfe/parallel/p2mapper.hh
View file @ 48b9358e
......@@ -8,45 +8,51 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
/** \brief Mimic a dune-grid mapper for a P2 space, using the dune-functions dof ordering of such a space
*/
template<class GridView>
class P2BasisMapper
namespace Dune::GFE
{
typedef typename GridView::Grid::template Codim<0>::Entity Element;
public:
typedef typename Dune::Functions::LagrangeBasis<GridView,2>::MultiIndex::value_type Index;
/** \brief Mimic a dune-grid mapper for a P2 space, using the dune-functions dof ordering of such a space
*/
template<class GridView>
class P2BasisMapper
{
typedef typename GridView::Grid::template Codim<0>::Entity Element;
public:
P2BasisMapper(const GridView& gridView)
: p2Basis_(gridView)
{}
typedef typename Dune::Functions::LagrangeBasis<GridView,2>::MultiIndex::value_type Index;
std::size_t size() const
{
return p2Basis_.size();
}
P2BasisMapper(const GridView& gridView)
: p2Basis_(gridView)
{}
template <class Entity>
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
auto localView = p2Basis_.localView();
localView.bind(entity);
std::size_t size() const
{
return p2Basis_.size();
}
for (size_t i=0; i<localView.size(); i++)
template <class Entity>
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
auto localView = p2Basis_.localView();
localView.bind(entity);
for (size_t i=0; i<localView.size(); i++)
{
result = localView.index(i);
return true;
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
{
result = localView.index(i);
return true;
}
}
return false;
}
return false;
}
Dune::Functions::LagrangeBasis<GridView,2> p2Basis_;
};
Dune::Functions::LagrangeBasis<GridView,2> p2Basis_;
};
} // namespace Dune::GFE
#endif
dune/gfe/parallel/vectorcommunicator.hh
View file @ 48b9358e
......@@ -7,98 +7,103 @@
#include <dune/gfe/parallel/mpifunctions.hh>
template<typename GUIndex, typename Communicator, typename VectorType>
class VectorCommunicator {
namespace Dune::GFE
{
struct TransferVectorTuple {
typedef typename VectorType::value_type EntryType;
template<typename GUIndex, typename Communicator, typename VectorType>
class VectorCommunicator {
size_t globalIndex_;
EntryType value_;
struct TransferVectorTuple {
typedef typename VectorType::value_type EntryType;
TransferVectorTuple() {}
TransferVectorTuple(const size_t& r, const EntryType& e)
: globalIndex_(r),
value_(e) {}
};
size_t globalIndex_;
EntryType value_;
TransferVectorTuple() {}
TransferVectorTuple(const size_t& r, const EntryType& e)
: globalIndex_(r),
value_(e) {}
};
private:
void transferVector(const VectorType& localVector) {
// Create vector for transfer data
std::vector<TransferVectorTuple> localVectorEntries;
private:
void transferVector(const VectorType& localVector) {
// Create vector for transfer data
std::vector<TransferVectorTuple> localVectorEntries;
// Translate vector entries
for (size_t k=0; k<localVector.size(); k++)
localVectorEntries.push_back(TransferVectorTuple(guIndex.index(k), localVector[k]));
// Translate vector entries
for (size_t k=0; k<localVector.size(); k++)
localVectorEntries.push_back(TransferVectorTuple(guIndex.index(k), localVector[k]));
// Get number of vector entries on each process
localVectorEntriesSizes = MPIFunctions::shareSizes(communicator_, localVectorEntries.size());
// Get number of vector entries on each process
localVectorEntriesSizes = MPIFunctions::shareSizes(communicator_, localVectorEntries.size());
// Get vector entries from every process
globalVectorEntries = MPIFunctions::gatherv(communicator_, localVectorEntries, localVectorEntriesSizes, root_rank);
}
// Get vector entries from every process
globalVectorEntries = MPIFunctions::gatherv(communicator_, localVectorEntries, localVectorEntriesSizes, root_rank);
}
public:
VectorCommunicator(const GUIndex& gi,
const Communicator& communicator,
const int& root)
: guIndex(gi), communicator_(communicator), root_rank(root)
{}
public:
VectorCommunicator(const GUIndex& gi,
const Communicator& communicator,
const int& root)
: guIndex(gi), communicator_(communicator), root_rank(root)
{}
VectorType reduceAdd(const VectorType& localVector)
{
transferVector(localVector);
VectorType reduceAdd(const VectorType& localVector)
{
transferVector(localVector);
VectorType globalVector(guIndex.size());
globalVector = 0;
VectorType globalVector(guIndex.size());
globalVector = 0;
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVector[globalVectorEntries[k].globalIndex_] += globalVectorEntries[k].value_;
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVector[globalVectorEntries[k].globalIndex_] += globalVectorEntries[k].value_;
return globalVector;
}
return globalVector;
}
VectorType reduceCopy(const VectorType& localVector)
{
transferVector(localVector);
VectorType reduceCopy(const VectorType& localVector)
{
transferVector(localVector);
VectorType globalVector(guIndex.size());
VectorType globalVector(guIndex.size());
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVector[globalVectorEntries[k].globalIndex_] = globalVectorEntries[k].value_;
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVector[globalVectorEntries[k].globalIndex_] = globalVectorEntries[k].value_;
return globalVector;
}
return globalVector;
}
VectorType scatter(const VectorType& global)
{
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVectorEntries[k].value_ = global[globalVectorEntries[k].globalIndex_];
VectorType scatter(const VectorType& global)
{
for (size_t k = 0; k < globalVectorEntries.size(); ++k)
globalVectorEntries[k].value_ = global[globalVectorEntries[k].globalIndex_];
const int localSize = localVectorEntriesSizes[communicator_.rank()];
const int localSize = localVectorEntriesSizes[communicator_.rank()];
// Create vector for transfer data
std::vector<TransferVectorTuple> localVectorEntries(localSize);
// Create vector for transfer data
std::vector<TransferVectorTuple> localVectorEntries(localSize);
MPIFunctions::scatterv(communicator_, localVectorEntries, globalVectorEntries, localVectorEntriesSizes, root_rank);
MPIFunctions::scatterv(communicator_, localVectorEntries, globalVectorEntries, localVectorEntriesSizes, root_rank);
// Create vector for local solution
VectorType x(localSize);
// Create vector for local solution
VectorType x(localSize);
// And translate solution again
for (size_t k = 0; k < localVectorEntries.size(); ++k)
x[k] = localVectorEntries[k].value_;
// And translate solution again
for (size_t k = 0; k < localVectorEntries.size(); ++k)
x[k] = localVectorEntries[k].value_;
return x;
}
return x;
}
private:
const GUIndex& guIndex;
const Communicator& communicator_;
int root_rank;
private:
const GUIndex& guIndex;
const Communicator& communicator_;
int root_rank;
std::vector<int> localVectorEntriesSizes;
std::vector<TransferVectorTuple> globalVectorEntries;
};
std::vector<int> localVectorEntriesSizes;
std::vector<TransferVectorTuple> globalVectorEntries;
};
} // namespace Dune::GFE
#endif
dune/gfe/polardecomposition.hh
View file @ 48b9358e
......@@ -13,7 +13,8 @@
// Two Methods to compute the Polar Factor of a 3x3 matrix
///////////////////////////////////////////////////////////////////////////////////////////
namespace Dune::GFE {
namespace Dune::GFE
{
class PolarDecomposition
{
......@@ -321,6 +322,7 @@ namespace Dune::GFE {
return v;
}
};
}
} // namespace Dune::GFE
#endif
dune/gfe/riemannianpnsolver.cc
View file @ 48b9358e
......@@ -19,7 +19,7 @@
#include <dune/gfe/parallel/vectorcommunicator.hh>
template <class Basis, class TargetSpace, class Assembler>
void RiemannianProximalNewtonSolver<Basis, TargetSpace, Assembler>::
void Dune::GFE::RiemannianProximalNewtonSolver<Basis, TargetSpace, Assembler>::
setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
......@@ -65,7 +65,7 @@ setup(const GridType& grid,
}
template <class Basis, class TargetSpace, class Assembler>
void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
void Dune::GFE::RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::
setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
......@@ -153,8 +153,8 @@ setup(const GridType& grid,
// Write all intermediate solutions, if requested
if (instrumented_
&& dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
&& dynamic_cast<::IterativeSolver<CorrectionType>*>(innerSolver_.get()))
dynamic_cast<::IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
// ////////////////////////////////////////////////////////////
// Create Hessian matrix and its occupation structure
......@@ -168,7 +168,7 @@ setup(const GridType& grid,
template <class Basis, class TargetSpace, class Assembler>
void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
void Dune::GFE::RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
{
int rank = grid_->comm().rank();
......
dune/gfe/riemannianpnsolver.hh
View file @ 48b9358e
......@@ -20,150 +20,155 @@
#include <dune/gfe/parallel/mapperfactory.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> >
namespace Dune::GFE
{
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 Dune::GFE::MapperFactory<Basis>::GlobalMapper GlobalMapper;
typedef typename Dune::GFE::MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
/** \brief Records information about the last run of the RiemannianProximalNewtonSolver
*
* This is used primarily for unit testing.
*/
struct Statistics
/** \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> >
{
std::size_t finalIteration;
typedef typename Basis::GridView::Grid GridType;
field_type finalEnergy;
};
const static int blocksize = TargetSpace::TangentVector::dimension;
public:
const static int gridDim = GridType::dimension;
RiemannianProximalNewtonSolver()
: IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
hessianMatrix_(nullptr), h1SemiNorm_(NULL)
{
std::fill(scaling_.begin(), scaling_.end(), 1.0);
}
/** \brief Set up the solver using a choldmod or umfpack 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);
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
const Dune::ParameterTree& parameterSet);
void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
{
scaling_ = scaling;
}
// Centralize the field type here
typedef double field_type;
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
innerSolver_->ignoreNodes_ = ignoreNodes_;
}
void solve();
[[deprecated]]
void setInitialSolution(const SolutionType& x) {
x_ = x;
}
void setInitialIterate(const SolutionType& x) {
x_ = x;
}
SolutionType getSol() const {return x_;}
// 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;
const Statistics& getStatistics() const {return statistics_;}
#if HAVE_MPI
typedef typename Dune::GFE::MapperFactory<Basis>::GlobalMapper GlobalMapper;
typedef typename Dune::GFE::MapperFactory<Basis>::LocalMapper LocalMapper;
#endif
protected:
/** \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)
{
std::fill(scaling_.begin(), scaling_.end(), 1.0);
}
/** \brief Set up the solver using a choldmod or umfpack 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);
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
const Dune::ParameterTree& parameterSet);
void setScaling(const Dune::FieldVector<double,blocksize>& scaling)
{
scaling_ = scaling;
}
void setIgnoreNodes(const Dune::BitSetVector<blocksize>& ignoreNodes)
{
ignoreNodes_ = &ignoreNodes;
innerSolver_->ignoreNodes_ = ignoreNodes_;
}
void solve();
[[deprecated]]
void setInitialSolution(const SolutionType& x) {
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_;
std::unique_ptr<GlobalMapper> globalMapper_;
#endif
/** \brief The grid */
const GridType* grid_;
/** \brief The grid */
const GridType* grid_;
/** \brief The solution vector */
SolutionType x_;
/** \brief The solution vector */
SolutionType x_;
/** \brief The initial regularization parameter for the proximal newton step */
double initialRegularization_;
double tolerance_;
/** \brief The initial regularization parameter for the proximal newton step */
double initialRegularization_;
double tolerance_;
/** \brief Regularization scaling */
Dune::FieldVector<double,blocksize> scaling_;
/** \brief Regularization scaling */
Dune::FieldVector<double,blocksize> scaling_;
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
/** \brief Maximum number of proximal-newton steps */
std::size_t maxProximalNewtonSteps_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief Hessian matrix */
std::unique_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
const Assembler* assembler_;
/** \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 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 Dirichlet nodes */
const Dune::BitSetVector<blocksize>* ignoreNodes_;
/** \brief The norm used to measure convergence for statistics*/
std::shared_ptr<H1SemiNorm<CorrectionType> > h1SemiNorm_;
/** \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 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 If set to true we log convergence speed and other stuff */
bool instrumented_;
/** \brief Output path for instrumented log */
std::string instrumentedPath_ = "/tmp";
/** \brief Output path for instrumented log */
std::string instrumentedPath_ = "/tmp";
/** \brief Norm type used for stopping criterion (default is infinity norm) */
enum class ErrorNormType {infinity, H1semi} normType_ = ErrorNormType::infinity;
/** \brief Norm type used for stopping criterion (default is infinity norm) */
enum class ErrorNormType {infinity, H1semi} normType_ = ErrorNormType::infinity;
/** \brief Norm type used for regularization term (default is Euclidean) */
enum class RegularizationNormType {Euclidean, H1semi, H1, L2} regNormType_ = RegularizationNormType::Euclidean;
/** \brief Norm type used for regularization term (default is Euclidean) */
enum class RegularizationNormType {Euclidean, H1semi, H1, L2} regNormType_ = RegularizationNormType::Euclidean;
/** \brief Store information about solver runs for unit testing */
Statistics statistics_;
/** \brief Store information about solver runs for unit testing */
Statistics statistics_;
};
};
} // namespace Dune::GFE
#include "riemannianpnsolver.cc"
......
dune/gfe/riemanniantrsolver.cc
View file @ 48b9358e
......@@ -27,7 +27,7 @@
#include <dune/gfe/parallel/vectorcommunicator.hh>
template <class Basis, class TargetSpace, class Assembler>
void RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
void Dune::GFE::RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
......@@ -66,7 +66,7 @@ setup(const GridType& grid,
}
template <class Basis, class TargetSpace, class Assembler>
void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
void Dune::GFE::RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
setup(const GridType& grid,
const Assembler* assembler,
const SolutionType& x,
......@@ -121,11 +121,11 @@ setup(const GridType& grid,
// Hack: the two-norm may not scale all that well, but it is fast!
auto baseNorm = std::make_shared<TwoNorm<CorrectionType> >();
auto baseSolver = std::make_shared<::LoopSolver<CorrectionType> >(baseSolverStep,
baseIterations,
baseTolerance,
baseNorm,
Solver::QUIET);
auto baseSolver = std::make_shared<Solvers::LoopSolver<CorrectionType> >(baseSolverStep,
baseIterations,
baseTolerance,
baseNorm,
Solver::QUIET);
#endif
#if HAVE_MPI
......@@ -179,11 +179,11 @@ setup(const GridType& grid,
#endif
h1SemiNorm_ = std::make_shared<H1SemiNorm<CorrectionType> >(A);
innerSolver_ = std::make_shared<::LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
innerTolerance_,
h1SemiNorm_,
Solver::QUIET);
innerSolver_ = std::make_shared<Solvers::LoopSolver<CorrectionType> >(mmgStep,
innerIterations_,
innerTolerance_,
h1SemiNorm_,
Solver::QUIET);
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a mass matrix to create a norm that's equivalent to the L2-norm
......@@ -204,8 +204,8 @@ setup(const GridType& grid,
// Write all intermediate solutions, if requested
if (instrumented_
&& dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
&& dynamic_cast<Solvers::IterativeSolver<CorrectionType>*>(innerSolver_.get()))
dynamic_cast<Solvers::IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
// ////////////////////////////////////////////////////////////
// Create Hessian matrix and its occupation structure
......@@ -248,7 +248,7 @@ setup(const GridType& grid,
#if HAVE_MPI
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef Dune::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1> > GlobalLeafP1Mapper;
typedef GFE::GlobalP1Mapper<Dune::Functions::LagrangeBasis<typename Basis::GridView,1> > GlobalLeafP1Mapper;
GlobalLeafP1Mapper p1Index(grid_->leafGridView());
typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> LeafP1LocalMapper;
......@@ -282,7 +282,7 @@ setup(const GridType& grid,
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
typedef Dune::Functions::LagrangeBasis<typename GridType::LevelGridView, 1> FEBasis;
typedef Dune::GlobalP1Mapper<FEBasis> GlobalLevelP1Mapper;
typedef GFE::GlobalP1Mapper<FEBasis> GlobalLevelP1Mapper;
GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+1));
GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i));
......@@ -328,15 +328,15 @@ setup(const GridType& grid,
template <class Basis, class TargetSpace, class Assembler>
void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
void Dune::GFE::RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
{
int rank = grid_->comm().rank();
MonotoneMGStep<MatrixType,CorrectionType>* mgStep = nullptr; // Non-shared pointer -- the innerSolver keeps the ownership
// if the inner solver is a monotone multigrid set up a max-norm trust-region
if (dynamic_cast<LoopSolver<CorrectionType>*>(innerSolver_.get())) {
auto loopSolver = std::dynamic_pointer_cast<LoopSolver<CorrectionType> >(innerSolver_);
if (dynamic_cast<Solvers::LoopSolver<CorrectionType>*>(innerSolver_.get())) {
auto loopSolver = std::dynamic_pointer_cast<Solvers::LoopSolver<CorrectionType> >(innerSolver_);
mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&loopSolver->getIterationStep());
}
......