diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 49dd39a82be03f70bff6e3e37f3ba3c6384c8122..b0fa64926e3517dacfd7dd44b82b3c64529c40e9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,12 +1,13 @@
---
-image: duneci/dune-fufem:2.4
+image: duneci/dune:2.6
+
+before_script:
+ - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-matrix-vector.git
+ - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
+ - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-fufem.git
dune-fufem--gcc:
- script:
- - dunecontrol --current all
- - dunecontrol --current make test
+ script: duneci-standard-test
dune-fufem--clang:
- script:
- - CMAKE_FLAGS="-DCMAKE_CXX_COMPILER=/usr/bin/clang++" dunecontrol --current all
- - dunecontrol --current make test
+ script: duneci-standard-test --opts=/duneci/opts.clang
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9fbd1a322c391729e9ccc28d4d8d261d0623f68a..6b1dac1137ce1bbfa9e4f9ef1adbfe559bc432ea 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.8.6)
+cmake_minimum_required(VERSION 2.8.12)
project(dune-gfe CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
@@ -17,7 +17,7 @@ include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
-
+dune_enable_all_packages()
add_subdirectory("src")
add_subdirectory("dune")
add_subdirectory("doc")
diff --git a/dune-gfe.pc.in b/dune-gfe.pc.in
index aa613574ece11eb6345c63c669159dace96886ef..5fa55f710e3a42565ae7bb9d89f9a8059be29655 100644
--- a/dune-gfe.pc.in
+++ b/dune-gfe.pc.in
@@ -9,6 +9,7 @@ DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: Dune dune-gfe module, contains the implementation of geodesic finite elements
-Requires: ${DEPENDENCIES}
-Libs:
+URL:
+Requires: dune-grid dune-uggrid dune-istl dune-localfunctions dune-functions dune-solvers dune-fufem
+Libs: -L${libdir}
Cflags: -I${includedir}
diff --git a/dune.module b/dune.module
index 7184b149308fb2b7d98aa91e0d1a30dfe000ba6e..26bcb99d069dd57cf83d3f23f07f8e188620457c 100644
--- a/dune.module
+++ b/dune.module
@@ -1,11 +1,11 @@
-
-#dune module information file#
-##############################
+################################
+# Dune module information file #
+################################
#Name of the module
-Module:dune-gfe
+Module: dune-gfe
Version: svn
Maintainer: oliver.sander@tu-dresden.de
#depending on
-Depends: dune-grid dune-istl dune-localfunctions dune-functions dune-solvers dune-fufem
+Depends: dune-grid dune-uggrid dune-istl dune-localfunctions dune-functions dune-solvers dune-fufem
Suggests: dune-foamgrid
diff --git a/dune/gfe/CMakeLists.txt b/dune/gfe/CMakeLists.txt
index 7ffd4e92a85db2846bb920bb776fc5f24aebd128..0cb24b316392d5583f7624b81ebda0e1b2457bf7 100644
--- a/dune/gfe/CMakeLists.txt
+++ b/dune/gfe/CMakeLists.txt
@@ -1,3 +1,2 @@
#install headers
-install(FILES gfe.hh DESTINATION include/dune/gfe)
-
+install(FILES gfe.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe)
diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index b4def6ba1165b4aac40784f58edd8049e2c078f9..3c23d23285b334c81971459c86a6cbd50b3af045 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -69,21 +69,6 @@ public:
orthonormalFrame.mv(embeddedGradient,gradient);
}
- void assembleEmbeddedHessian(const TargetSpace& x,
- Dune::FieldMatrix<ctype,embeddedSize,embeddedSize>& matrix) const
- {
- matrix = 0;
- for (size_t i=0; i<coefficients_.size(); i++) {
- Dune::SymmetricMatrix<ctype,embeddedSize> tmp = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x);
- for (int j=0; j<embeddedSize; j++)
- for (int k=0; k<=j; k++) {
- matrix[j][k] += weights_[i] * tmp(j,k);
- if (j!=k)
- matrix[k][j] += weights_[i] * tmp(j,k);
- }
- }
- }
-
void assembleEmbeddedHessian(const TargetSpace& x,
Dune::SymmetricMatrix<ctype,embeddedSize>& matrix) const
{
@@ -92,6 +77,7 @@ public:
matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
}
+ // TODO Use a Symmetric matrix for the result
void assembleHessian(const TargetSpace& x,
Dune::FieldMatrix<ctype,size,size>& matrix) const
{
diff --git a/dune/gfe/averageinterface.hh b/dune/gfe/averageinterface.hh
index cdcd8e83af660086ccc9b462021df575ef7ef89d..c57c98a5ab8749cae92a86a018d8377a1eaa2696 100644
--- a/dune/gfe/averageinterface.hh
+++ b/dune/gfe/averageinterface.hh
@@ -9,6 +9,8 @@
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
+#include dune/matrix-vector/crossproduct.hh>
+
#include <dune/fufem/dgindexset.hh>
#include <dune/fufem/arithmetic.hh>
#include <dune/fufem/surfmassmatrix.hh>
@@ -500,7 +502,7 @@ void computeTotalForceAndTorque(const BoundaryPatch<GridView>& interface,
neumannFunction.evaluateLocal(*it->inside(), quadPos, value);
totalForce.axpy(quad[ip].weight() * integrationElement, value);
- totalTorque.axpy(quad[ip].weight() * integrationElement, Arithmetic::crossProduct(worldPos-center,value));
+ totalTorque.axpy(quad[ip].weight() * integrationElement, MatrixVector::crossProduct(worldPos-center,value));
}
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index 1ff24720b0da85f5b683c107a8a5519095989c9e..bcc3b95092dbabde8eb16a88d23ee2c6995eb0a0 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -45,10 +45,10 @@ public:
/** \brief Specialize for scalar bases, here we cannot call tree().child() */
template <class GridView, int order, std::size_t i>
-class LocalFiniteElementFactory<Dune::Functions::PQkNodalBasis<GridView,order>,i>
+class LocalFiniteElementFactory<Dune::Functions::LagrangeBasis<GridView,order>,i>
{
public:
- static auto get(const typename Dune::Functions::PQkNodalBasis<GridView,order>::LocalView& localView,
+ static auto get(const typename Dune::Functions::LagrangeBasis<GridView,order>::LocalView& localView,
std::integral_constant<std::size_t, i> iType)
-> decltype(localView.tree().finiteElement())
{
diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index 9abe1085e076a6370fde5187a6bdedb7a7d21e6a..97239597566cf33e209608d489acf58253925f0d 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -1,10 +1,12 @@
#ifndef COSSERAT_VTK_WRITER_HH
#define COSSERAT_VTK_WRITER_HH
+#include <dune/common/version.hh>
+
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/io/file/vtk/pvtuwriter.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
@@ -129,7 +131,7 @@ public:
// Downsample 3rd-order functions onto a P2-space. That's all VTK can visualize today.
if (order>=3)
{
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,2> P2Basis;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,2> P2Basis;
P2Basis p2Basis(gridView);
std::vector<RigidBodyMotion<double,3> > downsampledConfig;
@@ -147,8 +149,8 @@ public:
for (const auto t : gridView.indexSet().types(0))
numElements[t] = 0;
- for (auto it = gridView.template begin<0,Dune::Interior_Partition>(); it != gridView.template end<0,Dune::Interior_Partition>(); ++it)
- numElements[it.type()]++;
+ for (auto&& t : elements(gridView, Dune::Partitions::interior))
+ numElements[t.type()]++;
std::size_t totalNumElements = 0;
for (const auto nE : numElements)
@@ -169,23 +171,30 @@ public:
connectivitySize += ((vtkOrder==2) ? 8 : 4) * nE.second;
else if (nE.first.isTriangle())
connectivitySize += ((vtkOrder==2) ? 6 : 3) * nE.second;
+ else if (nE.first.isHexahedron())
+ connectivitySize += ((vtkOrder==2) ? 20 : 8) * nE.second;
else
DUNE_THROW(Dune::IOError, "Unsupported element type '" << nE.first << "' found!");
}
std::vector<int> connectivity(connectivitySize);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis.localIndexSet();
+#endif
size_t i=0;
for (const auto& element : elements(gridView, Dune::Partitions::interior))
{
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
if (element.type().isQuadrilateral())
{
if (vtkOrder==2)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(2);
connectivity[i++] = localIndexSet.index(8);
@@ -195,31 +204,127 @@ public:
connectivity[i++] = localIndexSet.index(5);
connectivity[i++] = localIndexSet.index(7);
connectivity[i++] = localIndexSet.index(3);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(2);
+ connectivity[i++] = localView.index(8);
+ connectivity[i++] = localView.index(6);
+
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(5);
+ connectivity[i++] = localView.index(7);
+ connectivity[i++] = localView.index(3);
+#endif
}
else // first order
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(1);
connectivity[i++] = localIndexSet.index(3);
connectivity[i++] = localIndexSet.index(2);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(3);
+ connectivity[i++] = localView.index(2);
+#endif
}
}
if (element.type().isTriangle())
{
if (vtkOrder==2)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ connectivity[i++] = localIndexSet.index(0);
+ connectivity[i++] = localIndexSet.index(2);
+ connectivity[i++] = localIndexSet.index(5);
+ connectivity[i++] = localIndexSet.index(1);
+ connectivity[i++] = localIndexSet.index(4);
+ connectivity[i++] = localIndexSet.index(3);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(2);
+ connectivity[i++] = localView.index(5);
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(4);
+ connectivity[i++] = localView.index(3);
+#endif
+ }
+ else // first order
+ {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ connectivity[i++] = localIndexSet.index(0);
+ connectivity[i++] = localIndexSet.index(1);
+ connectivity[i++] = localIndexSet.index(2);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(2);
+#endif
+ }
+ }
+ if (element.type().isHexahedron())
+ {
+ if (vtkOrder==2)
+ {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(2);
connectivity[i++] = localIndexSet.index(5);
connectivity[i++] = localIndexSet.index(1);
connectivity[i++] = localIndexSet.index(4);
connectivity[i++] = localIndexSet.index(3);
+#else
+ // Corner dofs
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(2);
+ connectivity[i++] = localView.index(8);
+ connectivity[i++] = localView.index(6);
+
+ connectivity[i++] = localView.index(18);
+ connectivity[i++] = localView.index(20);
+ connectivity[i++] = localView.index(26);
+ connectivity[i++] = localView.index(24);
+
+ // Edge dofs
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(5);
+ connectivity[i++] = localView.index(7);
+ connectivity[i++] = localView.index(3);
+
+ connectivity[i++] = localView.index(19);
+ connectivity[i++] = localView.index(23);
+ connectivity[i++] = localView.index(25);
+ connectivity[i++] = localView.index(21);
+
+ connectivity[i++] = localView.index(9);
+ connectivity[i++] = localView.index(11);
+ connectivity[i++] = localView.index(17);
+ connectivity[i++] = localView.index(15);
+#endif
}
else // first order
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
connectivity[i++] = localIndexSet.index(0);
connectivity[i++] = localIndexSet.index(1);
+ connectivity[i++] = localIndexSet.index(3);
connectivity[i++] = localIndexSet.index(2);
+ connectivity[i++] = localIndexSet.index(4);
+ connectivity[i++] = localIndexSet.index(5);
+ connectivity[i++] = localIndexSet.index(7);
+ connectivity[i++] = localIndexSet.index(6);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(1);
+ connectivity[i++] = localView.index(3);
+ connectivity[i++] = localView.index(2);
+ connectivity[i++] = localView.index(4);
+ connectivity[i++] = localView.index(5);
+ connectivity[i++] = localView.index(7);
+ connectivity[i++] = localView.index(6);
+#endif
}
}
}
@@ -237,6 +342,9 @@ public:
if (element.type().isTriangle())
offsetCounter += (vtkOrder==2) ? 6 : 3;
+ if (element.type().isHexahedron())
+ offsetCounter += (vtkOrder==2) ? 20 : 8;
+
offsets[i++] += offsetCounter;
}
@@ -251,6 +359,9 @@ public:
if (element.type().isTriangle())
cellTypes[i++] = (vtkOrder==2) ? 22 : 5;
+
+ if (element.type().isHexahedron())
+ cellTypes[i++] = (vtkOrder==2) ? 25 : 12;
}
vtkFile.cellTypes_ = cellTypes;
@@ -300,7 +411,7 @@ public:
std::vector<RealTuple<double,3> > displacementConfiguration = deformationConfiguration;
typedef typename GridType::LeafGridView GridView;
- typedef Dune::Functions::PQkNodalBasis<GridView,2> P2DeformationBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView,2> P2DeformationBasis;
P2DeformationBasis p2DeformationBasis(gridView);
if (order == 3)
@@ -337,7 +448,11 @@ public:
// dump point coordinates
writer.beginPoints();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
writer.addArray<float>("Coordinates", 3);
+#else
+ writer.addArray("Coordinates", 3, Dune::VTK::Precision::float32);
+#endif
writer.endPoints();
for (int i=0; i<gridView.comm().size(); i++)
@@ -352,9 +467,8 @@ public:
/////////////////////////////////////////////////////////////////////////////////
// Stupid: I can't directly get the number of Interior_Partition elements
- size_t numElements = 0;
- for (const auto& element : elements(gridView, Dune::Partitions::interior))
- numElements++;
+ size_t numElements = std::distance(gridView.template begin<0, Dune::Interior_Partition>(),
+ gridView.template end<0, Dune::Interior_Partition>());
std::ofstream outFile(fullfilename);
@@ -376,15 +490,20 @@ public:
outFile << " <Cells>" << std::endl;
outFile << " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2DeformationBasis.localIndexSet();
+#endif
for (const auto& element : elements(gridView, Dune::Partitions::interior))
{
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
outFile << " ";
if (element.type().isQuadrilateral())
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
outFile << localIndexSet.index(0) << " ";
outFile << localIndexSet.index(2) << " ";
outFile << localIndexSet.index(8) << " ";
@@ -394,6 +513,17 @@ public:
outFile << localIndexSet.index(5) << " ";
outFile << localIndexSet.index(7) << " ";
outFile << localIndexSet.index(3) << " ";
+#else
+ outFile << localView.index(0) << " ";
+ outFile << localView.index(2) << " ";
+ outFile << localView.index(8) << " ";
+ outFile << localView.index(6) << " ";
+
+ outFile << localView.index(1) << " ";
+ outFile << localView.index(5) << " ";
+ outFile << localView.index(7) << " ";
+ outFile << localView.index(3) << " ";
+#endif
outFile << std::endl;
}
}
diff --git a/dune/gfe/embeddedglobalgfefunction.hh b/dune/gfe/embeddedglobalgfefunction.hh
index 3bad3c1d3c634fe09120012245ec2c30138fe4bb..86731241e589908cbcbe83106dfdad849463a2d3 100644
--- a/dune/gfe/embeddedglobalgfefunction.hh
+++ b/dune/gfe/embeddedglobalgfefunction.hh
@@ -68,17 +68,25 @@ public:
typename TargetSpace::CoordinateType operator()(const Element& element, const Dune::FieldVector<ctype,gridDim>& local) const
{
auto localView = basis_.localView();
- auto localIndexSet = basis_.localIndexSet();
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ auto localIndexSet = basis_.localIndexSet();
localIndexSet.bind(localView);
auto numOfBaseFct = localIndexSet.size();
+#else
+ auto numOfBaseFct = localView.size();
+#endif
// Extract local coefficients
std::vector<TargetSpace> localCoeff(numOfBaseFct);
for (size_t i=0; i<numOfBaseFct; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localCoeff[i] = coefficients_[localIndexSet.index(i)];
+#else
+ localCoeff[i] = coefficients_[localView.index(i)];
+#endif
// create local gfe function
LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
@@ -96,17 +104,26 @@ public:
Dune::FieldMatrix<ctype, embeddedDim, gridDim> derivative(const Element& element, const Dune::FieldVector<ctype,gridDim>& local) const
{
auto localView = basis_.localView();
- auto localIndexSet = basis_.localIndexSet();
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ auto localIndexSet = basis_.localIndexSet();
localIndexSet.bind(localView);
int numOfBaseFct = localIndexSet.size();
+#else
+ auto numOfBaseFct = localView.size();
+#endif
+
// Extract local coefficients
std::vector<TargetSpace> localCoeff(numOfBaseFct);
- for (int i=0; i<numOfBaseFct; i++)
+ for (decltype(numOfBaseFct) i=0; i<numOfBaseFct; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localCoeff[i] = coefficients_[localIndexSet.index(i)];
+#else
+ localCoeff[i] = coefficients_[localView.index(i)];
+#endif
// create local gfe function
LocalInterpolationRule localInterpolationRule(localView.tree().finiteElement(),localCoeff);
diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index ffea522fe78a18a854e46244d0d35db8be3e5d85..0e737066c64874735400fdf95ac192b77b93046f 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -76,7 +76,9 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
@@ -85,7 +87,9 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
// Bind the local FE basis view to the current element
localView.bind(*it);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
const auto& lfe = localView.tree().finiteElement();
@@ -93,8 +97,13 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
for (size_t j=0; j<lfe.size(); j++) {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto iIdx = localIndexSet.index(i)[0];
auto jIdx = localIndexSet.index(j)[0];
+#else
+ auto iIdx = localView.index(i);
+ auto jIdx = localView.index(j);
+#endif
nb.add(iIdx, jIdx);
@@ -128,7 +137,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endit = basis_.gridView().template end<0,Dune::Interior_Partition> ();
@@ -136,7 +147,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
for( ; it != endit; ++it ) {
localView.bind(*it);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
const int numOfBaseFct = localView.tree().size();
@@ -144,7 +157,11 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
std::vector<TargetSpace> localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
+#else
+ localSolution[i] = sol[localView.index(i)];
+#endif
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
@@ -154,11 +171,19 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto row = localIndexSet.index(i)[0];
+#else
+ auto row = localView.index(i);
+#endif
for (int j=0; j<numOfBaseFct; j++ ) {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto col = localIndexSet.index(j)[0];
+#else
+ auto col = localView.index(j);
+#endif
hessian[row][col] += localStiffness_->A_[i][j];
}
@@ -166,7 +191,11 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
gradient[localIndexSet.index(i)[0]] += localGradient[i];
+#else
+ gradient[localView.index(i)] += localGradient[i];
+#endif
}
@@ -185,7 +214,9 @@ assembleGradient(const std::vector<TargetSpace>& sol,
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
@@ -194,7 +225,9 @@ assembleGradient(const std::vector<TargetSpace>& sol,
for (; it!=endIt; ++it) {
localView.bind(*it);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
// A 1d grid has two vertices
const auto nDofs = localView.tree().size();
@@ -203,7 +236,11 @@ assembleGradient(const std::vector<TargetSpace>& sol,
std::vector<TargetSpace> localSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
+#else
+ localSolution[i] = sol[localView.index(i)];
+#endif
// Assemble local gradient
std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
@@ -212,8 +249,11 @@ assembleGradient(const std::vector<TargetSpace>& sol,
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
grad[localIndexSet.index(i)[0]] += localGradient[i];
-
+#else
+ grad[localView.index(i)[0]] += localGradient[i];
+#endif
}
}
@@ -230,7 +270,9 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
ElementIterator it = basis_.gridView().template begin<0,Dune::Interior_Partition>();
ElementIterator endIt = basis_.gridView().template end<0,Dune::Interior_Partition>();
@@ -239,7 +281,9 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
for (; it!=endIt; ++it) {
localView.bind(*it);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
// Number of degrees of freedom on this element
size_t nDofs = localView.tree().size();
@@ -247,7 +291,11 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
std::vector<TargetSpace> localSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
+#else
+ localSolution[i] = sol[localView.index(i)[0]];
+#endif
energy += localStiffness_->energy(localView, localSolution);
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 645b519a5ad3ce17fc4cea5648d9aff54627cfc8..253bd102fda9c3fd1f60311e1e84379a306d24ad 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -49,25 +49,25 @@ energy(const typename Basis::LocalView& localView,
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = quad[pt].position();
- const double integrationElement = element.geometry().integrationElement(quadPos);
+ const auto integrationElement = element.geometry().integrationElement(quadPos);
const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
- double weight = quad[pt].weight() * integrationElement;
+ auto weight = quad[pt].weight() * integrationElement;
// The derivative of the local function defined on the reference element
auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos);
- // The derivative of the function defined on the actual element
- typename LocalInterpolationRule::DerivativeType derivative(0);
-
- for (size_t comp=0; comp<referenceDerivative.N(); comp++)
- jacobianInverseTransposed.umv(referenceDerivative[comp], derivative[comp]);
-
- // Add the local energy density
- // The Frobenius norm is the correct norm here if the metric of TargetSpace is the identity.
- // (And if the metric of the domain space is the identity, which it always is here.)
- energy += weight * derivative.frobenius_norm2();
+ // Compute the Frobenius norm squared of the derivative. This is the correct term
+ // if both domain and target space use the metric inherited from an embedding.
+ for (size_t i=0; i<jacobianInverseTransposed.N(); i++)
+ for (int j=0; j<TargetSpace::embeddedDim; j++)
+ {
+ RT entry = 0;
+ for (size_t k=0; k<jacobianInverseTransposed.M(); k++)
+ entry += jacobianInverseTransposed[i][k] * referenceDerivative[j][k];
+ energy += weight * entry * entry;
+ }
}
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 1552f595fefad0f6710ed2a108530fe023af1011..3502f371dc54ec4756c2a7701604a6eda35d4153 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -111,7 +111,7 @@ public:
}
private:
- static Dune::FieldMatrix<RT,embeddedDim,embeddedDim> pseudoInverse(const Dune::FieldMatrix<RT,embeddedDim,embeddedDim>& dFdq,
+ static Dune::SymmetricMatrix<RT,embeddedDim> pseudoInverse(const Dune::SymmetricMatrix<RT,embeddedDim>& dFdq,
const TargetSpace& q)
{
const int shortDim = TargetSpace::TangentVector::dimension;
@@ -120,25 +120,25 @@ private:
Dune::FieldMatrix<RT,shortDim,embeddedDim> O = q.orthonormalFrame();
// compute A = O dFDq O^T
+ // TODO A really is symmetric, too
Dune::FieldMatrix<RT,shortDim,shortDim> A;
for (int i=0; i<shortDim; i++)
for (int j=0; j<shortDim; j++) {
A[i][j] = 0;
for (int k=0; k<embeddedDim; k++)
for (int l=0; l<embeddedDim; l++)
- A[i][j] += O[i][k]*dFdq[k][l]*O[j][l];
+ A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) *O[j][l];
}
A.invert();
- Dune::FieldMatrix<RT,embeddedDim,embeddedDim> result;
+ Dune::SymmetricMatrix<RT,embeddedDim> result;
+ result = 0.0;
for (int i=0; i<embeddedDim; i++)
- for (int j=0; j<embeddedDim; j++) {
- result[i][j] = 0;
+ for (int j=0; j<=i; j++)
for (int k=0; k<shortDim; k++)
for (int l=0; l<shortDim; l++)
- result[i][j] += O[k][i]*A[k][l]*O[l][j];
- }
+ result(i,j) += O[k][i]*A[k][l]*O[l][j];
return result;
}
@@ -195,10 +195,7 @@ evaluate(const Dune::FieldVector<ctype, dim>& local) const
solver.setup(&assembler,
initialIterate,
1e-14, // tolerance
- 100, // maxTrustRegionSteps
- 2, // initial trust region radius
- 100, // inner iterations
- 1e-14 // inner tolerance
+ 100 // maxNewtonSteps
);
solver.solve();
@@ -313,7 +310,7 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
- Dune::FieldMatrix<RT,embeddedDim,embeddedDim> dFdq(0);
+ Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
assembler.assembleEmbeddedHessian(q,dFdq);
const int shortDim = TargetSpace::TangentVector::dimension;
@@ -328,7 +325,7 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
A[i][j] = 0;
for (int k=0; k<embeddedDim; k++)
for (int l=0; l<embeddedDim; l++)
- A[i][j] += O[i][k]*dFdq[k][l]*O[j][l];
+ A[i][j] += O[i][k] * ((k>=l) ? dFdq(k,l) : dFdq(l,k)) * O[j][l];
}
//
@@ -421,7 +418,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
/** \todo Use a symmetric matrix here */
- Dune::FieldMatrix<RT,embeddedDim,embeddedDim> dFdq(0);
+ Dune::SymmetricMatrix<RT,embeddedDim> dFdq;
assembler.assembleEmbeddedHessian(q,dFdq);
@@ -436,7 +433,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
// dFDq is not invertible, if the target space is embedded into a higher-dimensional
// Euclidean space. Therefore we use its pseudo inverse. I don't think that is the
// best way, though.
- Dune::FieldMatrix<RT,embeddedDim,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
+ Dune::SymmetricMatrix<RT,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq,q);
//
Tensor3<RT,embeddedDim,embeddedDim,embeddedDim> dvDqF
@@ -488,7 +485,8 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
for (int j=0; j<embeddedDim; j++)
for (int k=0; k<dim; k++)
for (int l=0; l<embeddedDim; l++)
- result[i][j][k] += dFdqPseudoInv[j][l] * foo[i][l][k];
+ // TODO Smarter implementation of the product with a symmetric matrix
+ result[i][j][k] += ((j>=l) ? dFdqPseudoInv(j,l) : dFdqPseudoInv(l,j)) * foo[i][l][k];
}
@@ -623,7 +621,7 @@ public:
/** \brief Evaluate the derivative of the function */
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
- Dune::FieldMatrix<ctype, embeddedDim, dim> result(0);
+ DerivativeType result(0);
// get translation part
std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
@@ -634,7 +632,7 @@ public:
result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
// get orientation part
- Dune::FieldMatrix<ctype,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local);
+ Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local);
for (int i=0; i<4; i++)
for (int j=0; j<dim; j++)
result[3+i][j] = qResult[i][j];
diff --git a/dune/gfe/localgfetestfunctionbasis.hh b/dune/gfe/localgfetestfunctionbasis.hh
index ed1af0a3fd5045e7447325a8824ec258b9e88cf8..60062767fbc8bfd12da4f85a5caaab81d0d58486 100644
--- a/dune/gfe/localgfetestfunctionbasis.hh
+++ b/dune/gfe/localgfetestfunctionbasis.hh
@@ -109,7 +109,7 @@ 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>,1> Traits;
+ std::array<Dune::FieldMatrix<ctype, embeddedDim, dim>,spaceDim>> Traits;
/** \brief Constructor
*/
diff --git a/dune/gfe/localquickanddirtyfefunction.hh b/dune/gfe/localquickanddirtyfefunction.hh
index 277d98a33c2c4ddc33b2aace83b4dee0a3f37edc..5f46c8ad63304417666059814386400281aaf847 100644
--- a/dune/gfe/localquickanddirtyfefunction.hh
+++ b/dune/gfe/localquickanddirtyfefunction.hh
@@ -100,11 +100,54 @@ namespace Dune {
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());
+ // 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 recipies.
+ //
+ // 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);
+ return TargetSpace(c);
+ }
}
} // namespace GFE
diff --git a/dune/gfe/mixedgfeassembler.hh b/dune/gfe/mixedgfeassembler.hh
index 52d8dcf96386a32e8b4cfbe99ccc42f11b83c3c9..941b035f465d76b29e9fb030997a7c45b122dae2 100644
--- a/dune/gfe/mixedgfeassembler.hh
+++ b/dune/gfe/mixedgfeassembler.hh
@@ -95,13 +95,16 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
// Loop over grid elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
// Add element stiffness matrix onto the global stiffness matrix
@@ -114,6 +117,18 @@ getMatrixPattern(Dune::MatrixIndexSet& nb00,
{
// The global index of the j-th local degree of freedom of the element 'e'
auto col = localIndexSet.index(j);
+#else
+ // Add element stiffness matrix onto the global stiffness matrix
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ // The global index of the i-th local degree of freedom of the element 'e'
+ auto row = localView.index(i);
+
+ for (size_t j=0; j<localView.size(); j++ )
+ {
+ // The global index of the j-th local degree of freedom of the element 'e'
+ auto col = localView.index(j);
+#endif
if (row[0]==0 and col[0]==0)
nb00.add(row[1],col[1]);
@@ -167,13 +182,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
using namespace Dune::TypeTree::Indices;
const int nDofs0 = localView.tree().child(_0).finiteElement().size();
@@ -185,10 +204,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
for (int i=0; i<nDofs0+nDofs1; i++)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
+#else
+ if (localView.index(i)[0] == 0)
+ localConfiguration0[i] = configuration0[localView.index(i)[1]];
+ else
+ localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
+#endif
}
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
@@ -202,11 +228,19 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0+nDofs1; i++)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto row = localIndexSet.index(i);
+#else
+ auto row = localView.index(i);
+#endif
for (int j=0; j<nDofs0+nDofs1; j++ )
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto col = localIndexSet.index(j);
+#else
+ auto col = localView.index(j);
+#endif
if (row[0]==0 and col[0]==0)
hessian[_0][_0][row[1]][col[1]] += localStiffness_->A00_[i][j];
@@ -222,10 +256,17 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
}
// Add local gradient to global gradient
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
gradient0[localIndexSet.index(i)[1]] += localGradient0[i];
else
gradient1[localIndexSet.index(i)[1]] += localGradient1[i-nDofs0];
+#else
+ if (localView.index(i)[0] == 0)
+ gradient0[localView.index(i)[1]] += localGradient0[i];
+ else
+ gradient1[localView.index(i)[1]] += localGradient1[i-nDofs0];
+#endif
}
}
@@ -287,14 +328,18 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
// A view on the FE basis on a single element
auto localView = basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = basis_.localIndexSet();
+#endif
// Loop over all elements
for (const auto& element : elements(basis_.gridView(), Dune::Partitions::interior))
{
// Bind the local FE basis view to the current element
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
// Number of degrees of freedom on this element
using namespace Dune::TypeTree::Indices;
@@ -306,10 +351,17 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
for (int i=0; i<nDofs0+nDofs1; i++)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
if (localIndexSet.index(i)[0] == 0)
localConfiguration0[i] = configuration0[localIndexSet.index(i)[1]];
else
localConfiguration1[i-nDofs0] = configuration1[localIndexSet.index(i)[1]];
+#else
+ if (localView.index(i)[0] == 0)
+ localConfiguration0[i] = configuration0[localView.index(i)[1]];
+ else
+ localConfiguration1[i-nDofs0] = configuration1[localView.index(i)[1]];
+#endif
}
energy += localStiffness_->energy(localView,
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 5282aaa6d4127534f13d93d093af6951ce8fba30..dff40c95860ada21929500f252d4af44dd57de5d 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -6,7 +6,7 @@
#include <dune/istl/io.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
@@ -81,12 +81,12 @@ setup(const GridType& grid,
// ////////////////////////////////
#ifdef HAVE_IPOPT
// First create an IPOpt base solver
- auto baseSolver0 = new QuadraticIPOptSolver<MatrixType00,CorrectionType0>;
- baseSolver0->verbosity_ = NumProc::QUIET;
- baseSolver0->tolerance_ = baseTolerance;
- auto baseSolver1 = new QuadraticIPOptSolver<MatrixType11,CorrectionType1>;
- baseSolver1->verbosity_ = NumProc::QUIET;
- baseSolver1->tolerance_ = baseTolerance;
+ auto baseSolver0 = std::make_shared<QuadraticIPOptSolver<MatrixType00,CorrectionType0>>();
+ baseSolver0->setVerbosity(NumProc::QUIET);
+ baseSolver0->setTolerance(baseTolerance);
+ auto baseSolver1 = std::make_shared<QuadraticIPOptSolver<MatrixType11,CorrectionType1>>();
+ baseSolver1->setVerbosity(NumProc::QUIET);
+ baseSolver1->setTolerance(baseTolerance);
#else
// First create a Gauss-seidel base solver
TrustRegionGSStep<MatrixType00, CorrectionType0>* baseSolverStep0 = new TrustRegionGSStep<MatrixType00, CorrectionType0>;
@@ -122,29 +122,29 @@ setup(const GridType& grid,
// Make pre and postsmoothers
- TrustRegionGSStep<MatrixType00, CorrectionType0>* presmoother0 = new TrustRegionGSStep<MatrixType00, CorrectionType0>;
- TrustRegionGSStep<MatrixType00, CorrectionType0>* postsmoother0 = new TrustRegionGSStep<MatrixType00, CorrectionType0>;
+ auto presmoother0 = std::make_shared<TrustRegionGSStep<MatrixType00, CorrectionType0> >();
+ auto postsmoother0 = std::make_shared<TrustRegionGSStep<MatrixType00, CorrectionType0> >();
mmgStep0 = new MonotoneMGStep<MatrixType00, CorrectionType0>;
mmgStep0->setMGType(mu, nu1, nu2);
mmgStep0->ignoreNodes_ = globalDirichletNodes0;
- mmgStep0->basesolver_ = baseSolver0;
+ mmgStep0->setBaseSolver(baseSolver0);
mmgStep0->setSmoother(presmoother0, postsmoother0);
- mmgStep0->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType0>();
- mmgStep0->verbosity_ = Solver::QUIET;
+ mmgStep0->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<CorrectionType0>>());
+ mmgStep0->setVerbosity(Solver::QUIET);
- TrustRegionGSStep<MatrixType11, CorrectionType1>* presmoother1 = new TrustRegionGSStep<MatrixType11, CorrectionType1>;
- TrustRegionGSStep<MatrixType11, CorrectionType1>* postsmoother1 = new TrustRegionGSStep<MatrixType11, CorrectionType1>;
+ auto presmoother1 = std::make_shared<TrustRegionGSStep<MatrixType11, CorrectionType1> >();
+ auto postsmoother1 = std::make_shared<TrustRegionGSStep<MatrixType11, CorrectionType1> >();
mmgStep1 = new MonotoneMGStep<MatrixType11, CorrectionType1>;
mmgStep1->setMGType(mu, nu1, nu2);
mmgStep1->ignoreNodes_ = globalDirichletNodes1;
- mmgStep1->basesolver_ = baseSolver1;
+ mmgStep1->setBaseSolver(baseSolver1);
mmgStep1->setSmoother(presmoother1, postsmoother1);
- mmgStep1->obstacleRestrictor_= new MandelObstacleRestrictor<CorrectionType1>();
- mmgStep1->verbosity_ = Solver::QUIET;
+ mmgStep1->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<CorrectionType1>>());
+ mmgStep1->setVerbosity(Solver::QUIET);
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
@@ -193,12 +193,6 @@ setup(const GridType& grid,
// Create the transfer operators
// ////////////////////////////////////
- for (size_t k=0; k<mmgStep0->mgTransfer_.size(); k++)
- delete(mmgStep0->mgTransfer_[k]);
-
- for (size_t k=0; k<mmgStep1->mgTransfer_.size(); k++)
- delete(mmgStep1->mgTransfer_[k]);
-
mmgStep0->mgTransfer_.resize(numLevels-1);
mmgStep1->mgTransfer_.resize(numLevels-1);
@@ -206,20 +200,20 @@ setup(const GridType& grid,
{
if (numLevels>1) {
typedef typename TruncatedCompressedMGTransfer<CorrectionType0>::TransferOperatorType TransferOperatorType;
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> >,
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> >,
FufemBasis0>(pkToP1TransferMatrix,p1Basis,*basis0_);
- mmgStep0->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType0>;
+ mmgStep0->mgTransfer_.back() = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType0>>();
Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
- dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType0>*>(mmgStep0->mgTransfer_.back())->setMatrix(topTransferOperator);
+ std::dynamic_pointer_cast<TruncatedCompressedMGTransfer<CorrectionType0>>(mmgStep0->mgTransfer_.back())->setMatrix(topTransferOperator);
for (size_t i=0; i<mmgStep0->mgTransfer_.size()-1; i++){
// Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
+ auto newTransferOp0 = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType0>>();
newTransferOp0->setup(*grid_,i+1,i+2);
mmgStep0->mgTransfer_[i] = newTransferOp0;
@@ -234,7 +228,7 @@ setup(const GridType& grid,
for (size_t i=0; i<mmgStep0->mgTransfer_.size(); i++){
// Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType0>* newTransferOp0 = new TruncatedCompressedMGTransfer<CorrectionType0>;
+ auto newTransferOp0 = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType0>>();
newTransferOp0->setup(*grid_,i,i+1);
mmgStep0->mgTransfer_[i] = newTransferOp0;
@@ -246,20 +240,20 @@ setup(const GridType& grid,
{
if (numLevels>1) {
typedef typename TruncatedCompressedMGTransfer<CorrectionType1>::TransferOperatorType TransferOperatorType;
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> >,
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> >,
FufemBasis1>(pkToP1TransferMatrix,p1Basis,*basis1_);
- mmgStep1->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType1>;
- Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
- dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType1>*>(mmgStep1->mgTransfer_.back())->setMatrix(topTransferOperator);
+ mmgStep1->mgTransfer_.back() = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType1>>();
+ std::shared_ptr<TransferOperatorType> topTransferOperator = std::make_shared<TransferOperatorType>(pkToP1TransferMatrix);
+ std::dynamic_pointer_cast<TruncatedCompressedMGTransfer<CorrectionType1>>(mmgStep1->mgTransfer_.back())->setMatrix(topTransferOperator);
for (size_t i=0; i<mmgStep1->mgTransfer_.size()-1; i++){
// Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType1>* newTransferOp1 = new TruncatedCompressedMGTransfer<CorrectionType1>;
+ auto newTransferOp1 = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType1>>();
newTransferOp1->setup(*grid_,i+1,i+2);
mmgStep1->mgTransfer_[i] = newTransferOp1;
}
@@ -273,7 +267,7 @@ setup(const GridType& grid,
for (size_t i=0; i<mmgStep1->mgTransfer_.size(); i++){
// Construct the local multigrid transfer matrix
- TruncatedCompressedMGTransfer<CorrectionType1>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType1>;
+ auto newTransferOp = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType1>>();
newTransferOp->setup(*grid_,i,i+1);
mmgStep1->mgTransfer_[i] = newTransferOp;
}
@@ -291,8 +285,8 @@ setup(const GridType& grid,
hasObstacle0_.resize(assembler->basis_.size({0}), true);
hasObstacle1_.resize(assembler->basis_.size({1}), true);
#endif
- mmgStep0->hasObstacle_ = &hasObstacle0_;
- mmgStep1->hasObstacle_ = &hasObstacle1_;
+ mmgStep0->setHasObstacles(hasObstacle0_);
+ mmgStep1->setHasObstacles(hasObstacle1_);
}
}
@@ -403,13 +397,13 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
mmgStep0->setProblem(stiffnessMatrix[_0][_0], corr_global[_0], residual[_0]);
trustRegionObstacles0 = trustRegion0.obstacles();
- mmgStep0->obstacles_ = &trustRegionObstacles0;
+ mmgStep0->setObstacles(trustRegionObstacles0);
mmgStep0->preprocess();
mmgStep1->setProblem(stiffnessMatrix[_1][_1], corr_global[_1], residual[_1]);
trustRegionObstacles1 = trustRegion1.obstacles();
- mmgStep1->obstacles_ = &trustRegionObstacles1;
+ mmgStep1->setObstacles(trustRegionObstacles1);
mmgStep1->preprocess();
diff --git a/dune/gfe/mixedriemanniantrsolver.hh b/dune/gfe/mixedriemanniantrsolver.hh
index d291e0e8a7b40c6bd83759b85ed9261169120e6a..197388783636b2427c11172d5cfb546e5f36d8c3 100644
--- a/dune/gfe/mixedriemanniantrsolver.hh
+++ b/dune/gfe/mixedriemanniantrsolver.hh
@@ -11,7 +11,7 @@
#include <dune/grid/utility/globalindexset.hh>
-#include <dune/functions/common/tuplevector.hh>
+#include <dune/common/tuplevector.hh>
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
@@ -47,7 +47,7 @@ class MixedRiemannianTrustRegionSolver
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::Functions::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > SolutionType;
+ typedef Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > SolutionType;
public:
diff --git a/dune/gfe/nonplanarcosseratshellenergy.hh b/dune/gfe/nonplanarcosseratshellenergy.hh
index b995ccaa9c677e98a15d9cb6c704ecd50bff6e17..7378f35614ae49253c39431f16f711830aa95718 100644
--- a/dune/gfe/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/nonplanarcosseratshellenergy.hh
@@ -5,6 +5,8 @@
#include <dune/common/parametertree.hh>
#include <dune/geometry/quadraturerules.hh>
+#include <dune/matrix-vector/crossproduct.hh>
+
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/boundarypatch.hh>
@@ -277,7 +279,7 @@ energy(const typename Basis::LocalView& localView,
const auto& localFiniteElement = localView.tree().finiteElement();
////////////////////////////////////////////////////////////////////////////////////
- // Construct a linear (i.e., non-constant!) normal field on the surface
+ // Construct a linear (i.e., non-constant!) normal field on each element
////////////////////////////////////////////////////////////////////////////////////
auto gridView = localView.globalBasis().gridView();
@@ -312,8 +314,6 @@ energy(const typename Basis::LocalView& localView,
const DT integrationElement = geometry.integrationElement(quadPos);
- DT weight = quad[pt].weight() * integrationElement;
-
// The value of the local function
RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
@@ -351,7 +351,7 @@ energy(const typename Basis::LocalView& localView,
aCovariant[i][j] = 0.0;
}
- aCovariant[2] = Arithmetic::crossProduct(aCovariant[0], aCovariant[1]);
+ aCovariant[2] = Dune::MatrixVector::crossProduct(aCovariant[0], aCovariant[1]);
aCovariant[2] /= aCovariant[2].two_norm();
auto aContravariant = aCovariant;
@@ -428,15 +428,18 @@ energy(const typename Basis::LocalView& localView,
//////////////////////////////////////////////////////////
// Add the membrane energy density
- energy += weight * (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_m(Ee)
- + weight * (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_m(Ee*b + c*Ke)
- + weight * Dune::Power<3>::eval(thickness_) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
- + weight * Dune::Power<5>::eval(thickness_) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
+ auto energyDensity = (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_m(Ee)
+ + (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_m(Ee*b + c*Ke)
+ + Dune::Power<3>::eval(thickness_) / 6.0 * W_mixt(Ee, c*Ke*b - 2*H*c*Ke)
+ + Dune::Power<5>::eval(thickness_) / 80.0 * W_mp( (Ee*b + c*Ke)*b);
// Add the bending energy density
- energy += weight * (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_curv(Ke)
- + weight * (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_curv(Ke*b)
- + weight * Dune::Power<5>::eval(thickness_) / 80.0 * W_curv(Ke*b*b);
+ energyDensity += (thickness_ - K*Dune::Power<3>::eval(thickness_) / 12.0) * W_curv(Ke)
+ + (Dune::Power<3>::eval(thickness_) / 12.0 - K * Dune::Power<5>::eval(thickness_) / 80.0)*W_curv(Ke*b)
+ + Dune::Power<5>::eval(thickness_) / 80.0 * W_curv(Ke*b*b);
+
+ // Add energy density
+ energy += quad[pt].weight() * integrationElement * energyDensity;
///////////////////////////////////////////////////////////
// Volume load contribution
diff --git a/dune/gfe/parallel/globalp1mapper.hh b/dune/gfe/parallel/globalp1mapper.hh
index 8e1371798ce300c90082f66614fd4ac557ffa86e..e7da2815b0f19eb91479592c20c2663d525585ff 100644
--- a/dune/gfe/parallel/globalp1mapper.hh
+++ b/dune/gfe/parallel/globalp1mapper.hh
@@ -14,11 +14,6 @@
// Include Dune header files
#include <dune/common/version.hh>
-/** include parallel capability */
-#if HAVE_MPI
- #include <dune/common/parallel/mpihelper.hh>
-#endif
-
namespace Dune {
template <class GridView>
diff --git a/dune/gfe/parallel/globalp2mapper.hh b/dune/gfe/parallel/globalp2mapper.hh
index a2f636baeda94c00fd4be074f3e7ac5787a7c12c..1bb6a3ede3732e9e59d3c1e50b1cc0b586ebe1f1 100644
--- a/dune/gfe/parallel/globalp2mapper.hh
+++ b/dune/gfe/parallel/globalp2mapper.hh
@@ -14,12 +14,7 @@
// Include Dune header files
#include <dune/common/version.hh>
-/** include parallel capability */
-#if HAVE_MPI
- #include <dune/common/parallel/mpihelper.hh>
-#endif
-
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
namespace Dune {
@@ -38,6 +33,9 @@ namespace Dune {
{
static_assert(GridView::dimension==2, "Only implemented for two-dimensional grids");
+ if (gridView.size(GeometryTypes::triangle)>1)
+ DUNE_THROW(NotImplemented, "GlobalP2Mapper only works for quad grids!");
+
GlobalIndexSet<GridView> globalVertexIndex(gridView,2);
GlobalIndexSet<GridView> globalEdgeIndex(gridView,1);
GlobalIndexSet<GridView> globalElementIndex(gridView,0);
@@ -46,21 +44,33 @@ namespace Dune {
size_ = globalVertexIndex.size(2) + globalEdgeIndex.size(1) + globalElementIndex.size(0);
auto localView = p2Mapper_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = p2Mapper_.localIndexSet();
+#endif
// Determine
for (const auto& element : elements(gridView))
{
localView.bind(element);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
// Loop over all local degrees of freedom
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
+#else
+ for (size_t i=0; i<localView.size(); i++)
+#endif
{
int codim = localView.tree().finiteElement().localCoefficients().localKey(i).codim();
int entity = localView.tree().finiteElement().localCoefficients().localKey(i).subEntity();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
int localIndex = localIndexSet.index(i);
+#else
+ auto localIndex = localView.index(i);
+#endif
int globalIndex;
switch (codim)
{
@@ -105,19 +115,29 @@ namespace Dune {
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
auto localView = p2Mapper_.localView();
- auto localIndexSet = p2Mapper_.localIndexSet();
localView.bind(entity);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ auto localIndexSet = p2Mapper_.localIndexSet();
localIndexSet.bind(localView);
+#endif
Index localIndex;
bool dofFound = false;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
+#else
+ for (size_t i=0; i<localView.size(); i++)
+#endif
{
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
{
dofFound = true;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndex = localIndexSet.index(i);
+#else
+ localIndex = localView.index(i);
+#endif
break;
}
}
@@ -134,14 +154,13 @@ namespace Dune {
return size_;
}
- Functions::PQkNodalBasis<GridView,2> p2Mapper_;
+ Functions::LagrangeBasis<GridView,2> p2Mapper_;
IndexMap localGlobalMap_;
- size_t nOwnedLocalEntity_;
size_t size_;
};
}
-#endif /* GLOBALUNIQUEINDEX_HH_ */
+#endif // DUNE_GFE_PARALLEL_GLOBALP2MAPPER_HH
diff --git a/dune/gfe/parallel/p2mapper.hh b/dune/gfe/parallel/p2mapper.hh
index d6eed91ebe341a57b976ac1e2707cb2f09932e4e..3f9384d3e290c001006b964790c47ebd4ef303c3 100644
--- a/dune/gfe/parallel/p2mapper.hh
+++ b/dune/gfe/parallel/p2mapper.hh
@@ -6,7 +6,7 @@
#include <dune/geometry/type.hh>
#include <dune/common/typetraits.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#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
*/
@@ -16,7 +16,7 @@ class P2BasisMapper
typedef typename GridView::Grid::template Codim<0>::Entity Element;
public:
- typedef typename Dune::Functions::PQkNodalBasis<GridView,2>::MultiIndex::value_type Index;
+ typedef typename Dune::Functions::LagrangeBasis<GridView,2>::MultiIndex::value_type Index;
P2BasisMapper(const GridView& gridView)
: p2Basis_(gridView)
@@ -31,18 +31,26 @@ public:
bool contains(const Entity& entity, uint subEntity, uint codim, Index& result) const
{
auto localView = p2Basis_.localView();
- auto localIndexSet = p2Basis_.localIndexSet();
localView.bind(entity);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ auto localIndexSet = p2Basis_.localIndexSet();
localIndexSet.bind(localView);
+#endif
- Index localIndex;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
for (size_t i=0; i<localIndexSet.size(); i++)
+#else
+ for (size_t i=0; i<localView.size(); i++)
+#endif
{
if (localView.tree().finiteElement().localCoefficients().localKey(i).subEntity() == subEntity
and localView.tree().finiteElement().localCoefficients().localKey(i).codim() == codim)
{
- localIndex = i;
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
result = localIndexSet.index(i)[0];
+#else
+ result = localView.index(i);
+#endif
return true;
}
}
@@ -50,7 +58,7 @@ public:
return false;
}
- Dune::Functions::PQkNodalBasis<GridView,2> p2Basis_;
+ Dune::Functions::LagrangeBasis<GridView,2> p2Basis_;
};
#endif
diff --git a/dune/gfe/periodic1dpq1nodalbasis.hh b/dune/gfe/periodic1dpq1nodalbasis.hh
index df31bc5a67c23c10d7fe27cab0d3aecdbef85844..78852d98420986db2c232baa51986842ebd47e41 100644
--- a/dune/gfe/periodic1dpq1nodalbasis.hh
+++ b/dune/gfe/periodic1dpq1nodalbasis.hh
@@ -7,11 +7,9 @@
#include <dune/localfunctions/lagrange/pqkfactory.hh>
-#include <dune/typetree/leafnode.hh>
-
#include <dune/functions/functionspacebases/nodes.hh>
-#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
namespace Dune {
@@ -20,48 +18,46 @@ namespace Functions {
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
-// PQ1NodeFactory
+// PQ1PreBasis
// PQ1NodeIndexSet
// PQ1Node
//
-// The factory allows to create the others and is the owner of possible shared
+// The pre-basis allows to create the others and is the owner of possible shared
// state. These three components do _not_ depend on the global basis or index
// set and can be used without a global basis.
// *****************************************************************************
-template<typename GV, typename ST, typename TP>
+template<typename GV>
class Periodic1DPQ1Node;
-template<typename GV, class MI, class TP, class ST>
+template<typename GV, class MI>
class Periodic1DPQ1NodeIndexSet;
-template<typename GV, class MI, class ST>
-class Periodic1DPQ1NodeFactory;
-
-template<typename GV, class MI, class ST>
-class Periodic1DPQ1NodeFactory
+template<typename GV, class MI>
+class Periodic1DPQ1PreBasis
{
static const int dim = GV::dimension;
public:
- /** \brief The grid view that the FE space is defined on */
+ //! The grid view that the FE basis is defined on
using GridView = GV;
- using size_type = ST;
- template<class TP>
- using Node = Periodic1DPQ1Node<GV, size_type, TP>;
+ //! Type used for indices and size information
+ using size_type = std::size_t;
- template<class TP>
- using IndexSet = Periodic1DPQ1NodeIndexSet<GV, MI, TP, ST>;
+ using Node = Periodic1DPQ1Node<GV>;
+
+ using IndexSet = Periodic1DPQ1NodeIndexSet<GV, MI>;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
- using SizePrefix = Dune::ReservedVector<size_type, 2>;
+ //! Type used for prefixes handed to the size() method
+ using SizePrefix = Dune::ReservedVector<size_type, 1>;
- /** \brief Constructor for a given grid view object */
- Periodic1DPQ1NodeFactory(const GridView& gv) :
+ //! Constructor for a given grid view object
+ Periodic1DPQ1PreBasis(const GridView& gv) :
gridView_(gv)
{}
@@ -75,16 +71,20 @@ public:
return gridView_;
}
- template<class TP>
- Node<TP> node(const TP& tp) const
+ //! Update the stored grid view, to be called if the grid has changed
+ void update (const GridView& gv)
+ {
+ gridView_ = gv;
+ }
+
+ Node makeNode() const
{
- return Node<TP>{tp};
+ return Node{};
}
- template<class TP>
- IndexSet<TP> indexSet() const
+ IndexSet makeIndexSet() const
{
- return IndexSet<TP>{*this};
+ return IndexSet{*this};
}
size_type size() const
@@ -102,7 +102,7 @@ public:
DUNE_THROW(RangeError, "Method size() can only be called for prefixes of length up to one");
}
- /** \todo This method has been added to the interface without prior discussion. */
+ //! Get the total dimension of the space spanned by this basis
size_type dimension() const
{
return size()-1;
@@ -119,25 +119,22 @@ public:
-template<typename GV, typename ST, typename TP>
+template<typename GV>
class Periodic1DPQ1Node :
- public LeafBasisNode<ST, TP>
+ public LeafBasisNode
{
static const int dim = GV::dimension;
static const int maxSize = StaticPower<2,GV::dimension>::power;
- using Base = LeafBasisNode<ST,TP>;
using FiniteElementCache = typename Dune::PQkLocalFiniteElementCache<typename GV::ctype, double, dim, 1>;
public:
- using size_type = ST;
- using TreePath = TP;
+ using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = typename FiniteElementCache::FiniteElementType;
- Periodic1DPQ1Node(const TreePath& treePath) :
- Base(treePath),
+ Periodic1DPQ1Node() :
finiteElement_(nullptr),
element_(nullptr)
{}
@@ -174,24 +171,25 @@ protected:
-template<typename GV, class MI, class TP, class ST>
+template<typename GV, class MI>
class Periodic1DPQ1NodeIndexSet
{
enum {dim = GV::dimension};
public:
- using size_type = ST;
+ using size_type = std::size_t;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
- using NodeFactory = Periodic1DPQ1NodeFactory<GV, MI, ST>;
+ using PreBasis = Periodic1DPQ1PreBasis<GV, MI>;
- using Node = typename NodeFactory::template Node<TP>;
+ using Node = Periodic1DPQ1Node<GV>;
- Periodic1DPQ1NodeIndexSet(const NodeFactory& nodeFactory) :
- nodeFactory_(&nodeFactory)
+ Periodic1DPQ1NodeIndexSet(const PreBasis& preBasis) :
+ preBasis_(&preBasis),
+ node_(nullptr)
{}
/** \brief Bind the view to a grid element
@@ -215,6 +213,7 @@ public:
*/
size_type size() const
{
+ assert(node_ != nullptr);
return node_->finiteElement().size();
}
@@ -222,7 +221,7 @@ public:
MultiIndex index(size_type i) const
{
Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
- const auto& gridIndexSet = nodeFactory_->gridView().indexSet();
+ const auto& gridIndexSet = preBasis_->gridView().indexSet();
const auto& element = node_->element();
//return {{ gridIndexSet.subIndex(element,localKey.subEntity(),dim) }};
@@ -237,18 +236,18 @@ public:
}
protected:
- const NodeFactory* nodeFactory_;
+ const PreBasis* preBasis_;
const Node* node_;
};
/** \brief Nodal basis of a scalar first-order Lagrangian finite element space
+ * on a one-dimensional domain with periodic boundary conditions
*
* \tparam GV The GridView that the space is defined on
- * \tparam ST The type used for local indices; global indices are FlatMultiIndex<ST>
*/
-template<typename GV, class ST = std::size_t>
-using Periodic1DPQ1NodalBasis = DefaultGlobalBasis<Periodic1DPQ1NodeFactory<GV, FlatMultiIndex<ST>, ST> >;
+template<typename GV>
+using Periodic1DPQ1NodalBasis = DefaultGlobalBasis<Periodic1DPQ1PreBasis<GV, FlatMultiIndex<std::size_t> > >;
} // end namespace Functions
} // end namespace Dune
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 1f53e565ad5abfb5560b098ae973c9003d37804b..0c3958e0a2176c154dde43b6e5708e7c2bffc606 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -5,7 +5,6 @@
#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
@@ -66,7 +65,7 @@ setup(const GridType& grid,
//////////////////////////////////////////////////////////////////
#if HAVE_MPI
- globalMapper_ = std::unique_ptr<GlobalMapper>(new GlobalMapper(grid_->leafGridView()));
+ globalMapper_ = std::make_unique<GlobalMapper>(grid_->leafGridView());
#endif
// ////////////////////////////////
@@ -130,7 +129,7 @@ setup(const GridType& grid,
operatorAssembler.assemble(laplaceStiffness, localA);
#if HAVE_MPI
- LocalMapper localMapper(grid_->leafGridView());
+ LocalMapper localMapper = MapperFactory<typename Basis::GridView,Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
@@ -195,7 +194,7 @@ setup(const GridType& grid,
// On the lower grid levels a hierarchy of P1/Q1 spaces is used again.
////////////////////////////////////////////////////////////////////////
- bool isP1Basis = std::is_same<Basis,Dune::Functions::PQkNodalBasis<typename Basis::GridView,1> >::value;
+ bool isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView,1> >::value;
if (isP1Basis)
mmgStep->mgTransfer_.resize(numLevels-1);
@@ -206,11 +205,11 @@ setup(const GridType& grid,
if (not isP1Basis)
{
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> > p1Basis(grid_->leafGridView());
TransferOperatorType pkToP1TransferMatrix;
assembleBasisInterpolationMatrix<TransferOperatorType,
- DuneFunctionsBasis<Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> >,
+ DuneFunctionsBasis<Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> >,
FufemBasis>(pkToP1TransferMatrix,p1Basis,basis);
#if HAVE_MPI
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
@@ -218,8 +217,8 @@ setup(const GridType& grid,
typedef Dune::GlobalP1Mapper<typename GridType::LeafGridView> GlobalLeafP1Mapper;
GlobalLeafP1Mapper p1Index(grid_->leafGridView());
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView, Dune::MCMGVertexLayout> LeafP1LocalMapper;
- LeafP1LocalMapper leafP1LocalMapper(grid_->leafGridView());
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> LeafP1LocalMapper;
+ LeafP1LocalMapper leafP1LocalMapper(grid_->leafGridView(), Dune::mcmgVertexLayout());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
@@ -251,9 +250,9 @@ setup(const GridType& grid,
GlobalLevelP1Mapper fineGUIndex(grid_->levelGridView(i+1));
GlobalLevelP1Mapper coarseGUIndex(grid_->levelGridView(i));
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LevelGridView, Dune::MCMGVertexLayout> LevelLocalMapper;
- LevelLocalMapper fineLevelLocalMapper(grid_->levelGridView(i+1));
- LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i));
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LevelGridView> LevelLocalMapper;
+ LevelLocalMapper fineLevelLocalMapper(grid_->levelGridView(i+1), Dune::mcmgVertexLayout());
+ LevelLocalMapper coarseLevelLocalMapper(grid_->levelGridView(i), Dune::mcmgVertexLayout());
#endif
typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
#if HAVE_MPI
@@ -296,12 +295,12 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
{
int rank = grid_->comm().rank();
- std::shared_ptr<MonotoneMGStep<MatrixType,CorrectionType> > mgStep;
+ 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_);
- mgStep = std::dynamic_pointer_cast<MonotoneMGStep<MatrixType,CorrectionType> >(loopSolver->iterationStep_);
+ mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&loopSolver->getIterationStep());
}
@@ -348,7 +347,7 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace>::solve()
grid_->leafGridView().comm(),
0);
- LocalMapper localMapper(grid_->leafGridView());
+ LocalMapper localMapper = MapperFactory<typename Basis::GridView,Basis>::createLocalMapper(grid_->leafGridView());
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
typename GridType::LeafGridView,
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 36e61e7077577646e1d7bfd7202d78e8b86048e7..add3d56d61c31853ab262531d23a93892bd08b1a 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -8,6 +8,8 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/norms/h1seminorm.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -26,12 +28,16 @@ template <typename GridView, typename Basis>
struct MapperFactory
{};
-/** \brief Specialization for PQ1NodalBasis */
+/** \brief Specialization for LagrangeBasis<1> */
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::PQkNodalBasis<GridView,1> >
+struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,1> >
{
typedef Dune::GlobalP1Mapper<GridView> GlobalMapper;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> LocalMapper;
+ static LocalMapper createLocalMapper(const GridView& gridView)
+ {
+ return LocalMapper(gridView, Dune::mcmgVertexLayout());
+ }
};
// This case is not going to actually work, but I need the specialization to make
@@ -44,15 +50,19 @@ struct MapperFactory<GridView, Dune::Functions::Periodic1DPQ1NodalBasis<GridView
};
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::PQkNodalBasis<GridView,2> >
+struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,2> >
{
typedef Dune::GlobalP2Mapper<GridView> GlobalMapper;
typedef P2BasisMapper<GridView> LocalMapper;
+ static LocalMapper createLocalMapper(const GridView& gridView)
+ {
+ return LocalMapper(gridView);
+ }
};
-/** \brief Specialization for PQ3NodalBasis */
+/** \brief Specialization for LagrangeBasis<3> */
template <typename GridView>
-struct MapperFactory<GridView, Dune::Functions::PQkNodalBasis<GridView,3> >
+struct MapperFactory<GridView, Dune::Functions::LagrangeBasis<GridView,3> >
{
// Error: we don't currently have a global P3 mapper
};
diff --git a/dune/gfe/rodassembler.cc b/dune/gfe/rodassembler.cc
index 558d284f8d4427d17cadfd8403eac1205ca8d492..9fef383e8321c874f2dde938c47329462e534de7 100644
--- a/dune/gfe/rodassembler.cc
+++ b/dune/gfe/rodassembler.cc
@@ -26,7 +26,9 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
// A view on the FE basis on a single element
auto localView = this->basis_.localView();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
auto localIndexSet = this->basis_.localIndexSet();
+#endif
ElementIterator it = this->basis_.gridView().template begin<0>();
ElementIterator endIt = this->basis_.gridView().template end<0>();
@@ -35,7 +37,9 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
for (; it!=endIt; ++it) {
localView.bind(*it);
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localIndexSet.bind(localView);
+#endif
// A 1d grid has two vertices
static const int nDofs = 2;
@@ -44,7 +48,11 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
std::vector<RigidBodyMotion<double,3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
localSolution[i] = sol[localIndexSet.index(i)[0]];
+#else
+ localSolution[i] = sol[localView.index(i)];
+#endif
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
@@ -55,7 +63,11 @@ assembleGradient(const std::vector<RigidBodyMotion<double,3> >& sol,
// Add to global gradient
for (int i=0; i<nDofs; i++)
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
grad[localIndexSet.index(i)[0]] += localGradient[i];
+#else
+ grad[localView.index(i)] += localGradient[i];
+#endif
}
diff --git a/dune/gfe/rodfactory.hh b/dune/gfe/rodfactory.hh
index edfd47705fd4132f10c7fb19867fce4fdce975e8..94e06e451f012a8fe396f6124877f578a59d0624 100644
--- a/dune/gfe/rodfactory.hh
+++ b/dune/gfe/rodfactory.hh
@@ -5,6 +5,8 @@
#include <dune/common/fvector.hh>
#include <dune/fufem/arithmetic.hh>
+#include <dune/matrix-vector/crossproduct.hh>
+
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
@@ -39,7 +41,7 @@ template <int dim>
Dune::FieldVector<double,3> zAxis(0);
zAxis[2] = 1;
- Dune::FieldVector<double,3> axis = Arithmetic::crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
+ Dune::FieldVector<double,3> axis = MatrixVector::crossProduct(Dune::FieldVector<double,3>(end-beginning), zAxis);
if (axis.two_norm() != 0)
axis /= -axis.two_norm();
diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index c215606277aa95e56c815b1be69c44286a90d70e..8d268fbd4ec251a107dbf9dd4a52f367bf464643 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -7,17 +7,17 @@
#include <dune/istl/matrix.hh>
#include <dune/geometry/quadraturerules.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include "localgeodesicfestiffness.hh"
#include "rigidbodymotion.hh"
template<class GridView, class RT>
class RodLocalStiffness
- : public LocalGeodesicFEStiffness<Dune::Functions::PQkNodalBasis<GridView,1>, RigidBodyMotion<RT,3> >
+ : public LocalGeodesicFEStiffness<Dune::Functions::LagrangeBasis<GridView,1>, RigidBodyMotion<RT,3> >
{
typedef RigidBodyMotion<RT,3> TargetSpace;
- typedef Dune::Functions::PQkNodalBasis<GridView,1> Basis;
+ typedef Dune::Functions::LagrangeBasis<GridView,1> Basis;
// grid types
typedef typename GridView::Grid::ctype DT;
@@ -63,13 +63,10 @@ public:
//! Constructor
RodLocalStiffness (const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
- : gridView_(gridView)
- {
- for (int i=0; i<3; i++) {
- K_[i] = K[i];
- A_[i] = A[i];
- }
- }
+ : K_(K),
+ A_(A),
+ gridView_(gridView)
+ {}
/** \brief Constructor setting shape constants and material parameters
\param A The rod section area
@@ -465,7 +462,7 @@ getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
// multiply with jacobian inverse
Dune::FieldVector<double,1> tmp(0);
inv.umv(shapeGrad[dof][0], tmp);
- shapeGrad[dof] = tmp;
+ shapeGrad[dof][0] = tmp;
}
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 1d9bb430fc1dfdc7df2662a53c1453134059bdf8..aa023750e5188f60422f3d30f516b97b443e5cec 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -892,6 +892,27 @@ public:
}
+ /** \brief Derivative of the map from unit quaternions to orthogonal matrices
+ */
+ static Tensor3<T,3,3,4> derivativeOfQuaternionToMatrix(const Dune::FieldVector<T,4>& p)
+ {
+ Tensor3<T,3,3,4> result;
+
+ result[0][0] = { 2*p[0], -2*p[1], -2*p[2], 2*p[3]};
+ result[0][1] = { 2*p[1], 2*p[0], -2*p[3], -2*p[2]};
+ result[0][2] = { 2*p[2], 2*p[3], 2*p[0], 2*p[1]};
+
+ result[1][0] = { 2*p[1], 2*p[0], 2*p[3], 2*p[2]};
+ result[1][1] = {-2*p[0], 2*p[1], -2*p[2], 2*p[3]};
+ result[1][2] = {-2*p[3], 2*p[2], 2*p[1], -2*p[0]};
+
+ result[2][0] = { 2*p[2], -2*p[3], 2*p[0], -2*p[1]};
+ result[2][1] = { 2*p[3], 2*p[2], 2*p[1], 2*p[0]};
+ result[2][2] = {-2*p[0], -2*p[1], 2*p[2], 2*p[3]};
+
+ return result;
+ }
+
/** \brief Set rotation from orthogonal matrix
We tacitly assume that the matrix really is orthogonal */
@@ -1127,6 +1148,30 @@ public:
return result;
}
+ /** \brief Project a vector in R^4 onto the unit quaternions
+ *
+ * \warning This is NOT the standard projection from R^{3 \times 3} onto SO(3)!
+ */
+ static Rotation<T,3> projectOnto(const CoordinateType& p)
+ {
+ Rotation<T,3> result(p);
+ result /= result.two_norm();
+ return result;
+ }
+
+ /** \brief Derivative of the projection of a vector in R^4 onto the unit quaternions
+ *
+ * \warning This is NOT the standard projection from R^{3 \times 3} onto SO(3)!
+ */
+ static auto derivativeOfProjection(const Dune::FieldVector<T,4>& p)
+ {
+ Dune::FieldMatrix<T,4,4> result;
+ for (int i=0; i<4; i++)
+ for (int j=0; j<4; j++)
+ result[i][j] = ( (i==j) - p[i]*p[j] / p.two_norm2() ) / p.two_norm();
+ return result;
+ }
+
/** \brief The global coordinates, if you really want them */
const CoordinateType& globalCoordinates() const {
return *this;
diff --git a/dune/gfe/symmetricmatrix.hh b/dune/gfe/symmetricmatrix.hh
index d6d32d7cddd4cd1ed791ea7eb0b37f0173694aec..a8d847bb45e34d10d7a466544112f0b4346a8f78 100644
--- a/dune/gfe/symmetricmatrix.hh
+++ b/dune/gfe/symmetricmatrix.hh
@@ -23,10 +23,7 @@ public:
enum {blocklevel = 0};
- /** \brief Default constructor
- *
- * Tensor is initialized containing zeros if no argument is given.
- * \param eye if true tensor is initialized as identity
+ /** \brief Default constructor, creates uninitialized matrix
*/
SymmetricMatrix()
{}
@@ -47,7 +44,7 @@ public:
* \param i line index
* \param j column index
* \note You need to know what you are doing: You can only access the lower
- * left triangular submatrix using this methods. It requires i<=j.
+ * left triangular submatrix using this methods. It requires i>=j.
*/
T& operator() (int i, int j)
{
@@ -59,7 +56,7 @@ public:
* \param i line index
* \param j column index
* \note You need to know what you are doing: You can only access the lower
- * left triangular submatrix using this methods. It requires i<=j.
+ * left triangular submatrix using this methods. It requires i>=j.
*/
const T& operator() (int i, int j) const
{
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index 2e11e9ea3c6a5c490c7e717302dd0513763af07f..b00b59e1b3d4e79c3a1c56786ae26ea0b7a9e883 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -13,39 +13,14 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::
setup(const AverageDistanceAssembler<TargetSpace>* assembler,
const TargetSpace& x,
double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int innerIterations,
- double innerTolerance)
+ int maxNewtonSteps)
{
assembler_ = assembler;
x_ = x;
tolerance_ = tolerance;
- maxTrustRegionSteps_ = maxTrustRegionSteps;
- initialTrustRegionRadius_ = initialTrustRegionRadius;
- innerIterations_ = innerIterations;
- innerTolerance_ = innerTolerance;
+ maxNewtonSteps_ = maxNewtonSteps;
this->verbosity_ = NumProc::QUIET;
minNumberOfIterations_ = 1;
-
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- // ////////////////////////////////
- // Create a projected gauss-seidel solver
- // ////////////////////////////////
-
- // First create a Gauss-seidel base solver
- innerSolverStep_ = std::auto_ptr<TrustRegionGSStep<MatrixType, CorrectionType> >(new TrustRegionGSStep<MatrixType, CorrectionType>);
-
- energyNorm_ = std::auto_ptr<EnergyNorm<MatrixType, CorrectionType> >(new EnergyNorm<MatrixType, CorrectionType>(*innerSolverStep_.get()));
-
- innerSolver_ = std::auto_ptr< ::LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(innerSolverStep_.get(),
- innerIterations,
- innerTolerance,
- energyNorm_.get(),
- Solver::QUIET));
-
- innerSolver_->useRelativeError_ = false;
-#endif
}
@@ -54,68 +29,34 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
{
assert(minNumberOfIterations_ > 0);
- MaxNormTrustRegion<blocksize,field_type> trustRegion(1, // we have only one block
- initialTrustRegionRadius_);
-
- field_type energy = assembler_->value(x_);
-
// /////////////////////////////////////////////////////
- // Trust-Region Solver
+ // Newton Solver
// /////////////////////////////////////////////////////
- for (size_t i=0; i<maxTrustRegionSteps_; i++) {
+ for (size_t i=0; i<maxNewtonSteps_; i++) {
if (this->verbosity_ == Solver::FULL) {
std::cout << "----------------------------------------------------" << std::endl;
- std::cout << " Trust-Region Step Number: " << i
- << ", radius: " << trustRegion.radius()
- << ", energy: " << energy << std::endl;
+ std::cout << " Newton Step Number: " << i
+ << ", energy: " << assembler_->value(x_) << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
}
- CorrectionType rhs(1); // length is 1 _block_
- CorrectionType corr(1); // length is 1 _block_
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- corr = 0;
-#endif
+ CorrectionType corr;
- MatrixType hesseMatrix(1,1);
+ MatrixType hesseMatrix;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- assembler_->assembleGradient(x_, rhs[0]);
- assembler_->assembleHessian(x_, hesseMatrix[0][0]);
-#else
- /** \todo Fix this sense copying */
- typename TargetSpace::EmbeddedTangentVector foo;
- assembler_->assembleEmbeddedGradient(x_, foo);
- rhs[0] = foo;
- assembler_->assembleEmbeddedHessian(x_, hesseMatrix[0][0]);
-#endif
+ typename TargetSpace::EmbeddedTangentVector rhs;
+ assembler_->assembleEmbeddedGradient(x_, rhs);
+ assembler_->assembleEmbeddedHessian(x_, hesseMatrix);
// The right hand side is the _negative_ gradient
rhs *= -1;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- dynamic_cast<LinearIterationStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->setProblem(hesseMatrix, corr, rhs);
-
- dynamic_cast<TrustRegionGSStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->obstacles_ = &trustRegion.obstacles();
-
- innerSolver_->preprocess();
-#endif
// /////////////////////////////
// Solve !
// /////////////////////////////
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- innerSolver_->solve();
-#else
Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
- GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
-#endif
-
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- corr = innerSolver_->iterationStep_->getSol();
-#endif
-
- //std::cout << "Corr: " << corr << std::endl;
+ GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix, corr, rhs, basis);
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
@@ -125,84 +66,31 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ std::cout << i+1 << " Newton steps were taken." << std::endl;
break;
}
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
-
+#if 0
TargetSpace newIterate = x_;
- newIterate = TargetSpace::exp(newIterate, corr[0]);
-
- field_type oldEnergy = energy;
- field_type energy = assembler_->value(newIterate);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- CorrectionType tmp(corr.size());
- tmp = 0;
- hesseMatrix.umv(corr, tmp);
- field_type modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-#else
- field_type modelDecrease = (rhs*corr) - 0.5 * hesseMatrix[0][0].energyScalarProduct(corr[0],corr[0]);
-#endif
+ newIterate = TargetSpace::exp(newIterate, corr);
- if (this->verbosity_ == NumProc::FULL) {
- std::cout << "Absolute model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "Relative model decrease: " << modelDecrease / energy
- << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
- }
-
- assert(modelDecrease >= -1e-15);
+ if (this->verbosity_ == NumProc::FULL)
+ {
+ field_type oldEnergy = assembler_->value(x_);
+ field_type energy = assembler_->value(newIterate);
- if (energy >= oldEnergy) {
- if (this->verbosity_ == NumProc::FULL)
+ if (energy >= oldEnergy)
printf("Richtung ist keine Abstiegsrichtung!\n");
}
- if (energy >= oldEnergy &&
- i>minNumberOfIterations_-1 &&
- (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Suspecting rounding problems" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
-
- x_ = newIterate;
- break;
- }
-
- // //////////////////////////////////////////////
- // Check for acceptance of the step
- // //////////////////////////////////////////////
- if ( (oldEnergy-energy) / modelDecrease > 0.9) {
- // very successful iteration
-
- x_ = newIterate;
- trustRegion.scale(2);
-
- } else if ( (oldEnergy-energy) / modelDecrease > 0.01
- || std::abs(oldEnergy-energy) < 1e-12) {
- // successful iteration
- x_ = newIterate;
-
- } else {
- // unsuccessful iteration
- trustRegion.scale(0.5);
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Unsuccessful iteration!" << std::endl;
- }
-
- // Write current energy
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "--- Current energy: " << energy << " ---" << std::endl;
-
+ // Actually take the step
+ x_ = newIterate;
+#else
+ x_ = TargetSpace::exp(x_, corr);
+#endif
}
}
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 10e2d87f60f9b6e29bc7dca3b8d05b72156ed133..568d8314ac7912b66d48d124c9aafca0bf9f2689 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -1,17 +1,7 @@
#ifndef TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
#define TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
-
-//#define USE_GAUSS_SEIDEL_SOLVER
-
-#include <dune/istl/matrix.hh>
-
-#include <dune/solvers/common/boxconstraint.hh>
-#ifdef USE_GAUSS_SEIDEL_SOLVER
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
-#endif
-#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/common/numproc.hh>
#include <dune/gfe/symmetricmatrix.hh>
@@ -28,31 +18,16 @@ class TargetSpaceRiemannianTRSolver
// Centralize the field type here
typedef typename TargetSpace::ctype field_type;
- // Some types that I need
- // The types have the dynamic outer type because the dune-solvers solvers expect
- // this sort of type.
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- typedef Dune::Matrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
- typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
-#else
- typedef Dune::Matrix<Dune::SymmetricMatrix<field_type, embeddedBlocksize> > MatrixType;
- typedef Dune::BlockVector<Dune::FieldVector<field_type, embeddedBlocksize> > CorrectionType;
-#endif
+ typedef Dune::SymmetricMatrix<field_type, embeddedBlocksize> MatrixType;
+ typedef Dune::FieldVector<field_type, embeddedBlocksize> CorrectionType;
public:
- TargetSpaceRiemannianTRSolver()
- // : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL)
- {}
-
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const AverageDistanceAssembler<TargetSpace>* assembler,
const TargetSpace& x,
double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int innerIterations,
- double innerTolerance);
+ int maxNewtonSteps);
void solve();
@@ -70,32 +45,13 @@ protected:
/** \brief Tolerance of the solver */
double tolerance_;
- /** \brief The initial trust-region radius in the maximum-norm */
- double initialTrustRegionRadius_;
-
- /** \brief Maximum number of trust-region steps */
- size_t maxTrustRegionSteps_;
-
- /** \brief Maximum number of multigrid iterations */
- int innerIterations_;
-
- /** \brief Error tolerance of the multigrid QP solver */
- double innerTolerance_;
+ /** \brief Maximum number of Newton steps */
+ size_t maxNewtonSteps_;
/** \brief The assembler for the average-distance functional */
const AverageDistanceAssembler<TargetSpace>* assembler_;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- /** \brief The solver for the quadratic inner problems */
- std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;
-
- /** \brief The iteration step for the quadratic inner problems */
- std::auto_ptr<TrustRegionGSStep<MatrixType, CorrectionType> > innerSolverStep_;
-
- /** \brief Norm for the quadratic inner problems */
- std::auto_ptr<EnergyNorm<MatrixType, CorrectionType> > energyNorm_;
-#endif
- /** \brief Specify a minimal number of iterations the trust-region solver has to do
+ /** \brief Specify a minimal number of iterations the Newton solver has to do
*
* This is needed when working with automatic differentiation. While a very low
* number of iterations may be enough to precisely compute the value of a
diff --git a/dune/gfe/tensorssd.hh b/dune/gfe/tensorssd.hh
index 43680fe5e7283f00ba434c1a59de2b9cda2893e4..ac83b94ce4d57e7721b1e6bed113f62101043eca 100644
--- a/dune/gfe/tensorssd.hh
+++ b/dune/gfe/tensorssd.hh
@@ -8,7 +8,8 @@
#include <array>
#include <dune/common/fmatrix.hh>
-
+#include <dune/istl/matrix.hh>
+
/** \brief A third-rank tensor with two static (SS) and one dynamic (D) dimension
*
* \tparam T Type of the entries
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index a1145af47f0c2f1990cd7031e6bbaedeeb8ac1cc..734188488d5accf9b290a535ec518148507d1350 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -413,10 +413,13 @@ public:
static DerivativeOfProjection derivativeOfProjection(const Dune::FieldVector<T,N>& p)
{
+ auto normSquared = p.two_norm2();
+ auto norm = std::sqrt(normSquared);
+
Dune::FieldMatrix<T,N,N> result;
for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
- result[i][j] = ( (i==j) - p[i]*p[j] / p.two_norm2() ) / p.two_norm();
+ result[i][j] = ( (i==j) - p[i]*p[j] / normSquared ) / norm;
return result;
}
diff --git a/dune/gfe/vertexnormals.hh b/dune/gfe/vertexnormals.hh
index fb899e0ee3af64a1e4ca608727db615999d0c113..ac067464b5bc2a3c36d4397dff7330b582f4d6eb 100644
--- a/dune/gfe/vertexnormals.hh
+++ b/dune/gfe/vertexnormals.hh
@@ -9,6 +9,8 @@
#include <dune/geometry/referenceelements.hh>
+#include <dune/matrix-vector/crossproduct.hh>
+
#include <dune/gfe/unitvector.hh>
/** \brief Compute averaged vertex normals for a 2d-in-3d grid
@@ -25,11 +27,11 @@ std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const
for (const auto& element : elements(gridView))
{
- for (int i=0; i<element.subEntities(2); i++)
+ for (std::size_t i=0; i<element.subEntities(2); i++)
{
auto cornerPos = Dune::ReferenceElements<double,2>::general(element.type()).position(i,2);
auto tangent = element.geometry().jacobianTransposed(cornerPos);
- auto cornerNormal = Arithmetic::crossProduct(tangent[0], tangent[1]);
+ auto cornerNormal = Dune::MatrixVector::crossProduct(tangent[0], tangent[1]);
cornerNormal /= cornerNormal.two_norm();
unscaledNormals[indexSet.subIndex(element,i,2)] += cornerNormal;
@@ -50,4 +52,4 @@ std::vector<UnitVector<typename GridView::ctype,3> > computeVertexNormals(const
}
-#endif
\ No newline at end of file
+#endif
diff --git a/dune/gfe/vtkfile.hh b/dune/gfe/vtkfile.hh
index db460f03996d8b96b01e1324a14d9c3b3a99348e..cf14e8ef653547c5f879bcf459f138c3cf5707e6 100644
--- a/dune/gfe/vtkfile.hh
+++ b/dune/gfe/vtkfile.hh
@@ -12,6 +12,7 @@
#endif
#include <dune/common/fvector.hh>
+#include <dune/common/version.hh>
// For parallel infrastructure stuff:
#include <dune/grid/io/file/vtk.hh>
@@ -55,19 +56,34 @@ namespace Dune {
writer.beginMain();
writer.beginPointData();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
writer.addArray<float>("director0", 3);
writer.addArray<float>("director1", 3);
writer.addArray<float>("director2", 3);
writer.addArray<float>("zCoord", 1);
+#else
+ writer.addArray("director0", 3, VTK::Precision::float32);
+ writer.addArray("director1", 3, VTK::Precision::float32);
+ writer.addArray("director2", 3, VTK::Precision::float32);
+ writer.addArray("zCoord", 1, VTK::Precision::float32);
+#endif
writer.endPointData();
writer.beginCellData();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
writer.addArray<float>("mycelldata", 1);
+#else
+ writer.addArray("mycelldata", 1, VTK::Precision::float32);
+#endif
writer.endCellData();
// dump point coordinates
writer.beginPoints();
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
writer.addArray<float>("Coordinates", 3);
+#else
+ writer.addArray("Coordinates", 3, VTK::Precision::float32);
+#endif
writer.endPoints();
for (int i=0; i<mpiHelper.size(); i++)
@@ -88,7 +104,11 @@ namespace Dune {
// Write vertex coordinates
outFile << " <Points>" << std::endl;
{ // extra parenthesis to control destruction of the pointsWriter object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<float> pointsWriter(outFile, "Coordinates", 3, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter pointsWriter(outFile, "Coordinates", 3, Dune::Indent(4), VTK::Precision::float32);
+#endif
for (size_t i=0; i<points_.size(); i++)
for (int j=0; j<3; j++)
pointsWriter.write(points_[i][j]);
@@ -98,19 +118,31 @@ namespace Dune {
// Write elements
outFile << " <Cells>" << std::endl;
{ // extra parenthesis to control destruction of the cellConnectivityWriter object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<int> cellConnectivityWriter(outFile, "connectivity", 1, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter cellConnectivityWriter(outFile, "connectivity", 1, Dune::Indent(4), VTK::Precision::int32);
+#endif
for (size_t i=0; i<cellConnectivity_.size(); i++)
cellConnectivityWriter.write(cellConnectivity_[i]);
}
{ // extra parenthesis to control destruction of the writer object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<int> cellOffsetsWriter(outFile, "offsets", 1, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter cellOffsetsWriter(outFile, "offsets", 1, Dune::Indent(4), VTK::Precision::int32);
+#endif
for (size_t i=0; i<cellOffsets_.size(); i++)
cellOffsetsWriter.write(cellOffsets_[i]);
}
{ // extra parenthesis to control destruction of the writer object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<unsigned int> cellTypesWriter(outFile, "types", 1, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter cellTypesWriter(outFile, "types", 1, Dune::Indent(4), VTK::Precision::uint32);
+#endif
for (size_t i=0; i<cellTypes_.size(); i++)
cellTypesWriter.write(cellTypes_[i]);
}
@@ -124,7 +156,11 @@ namespace Dune {
// Z coordinate for better visualization of wrinkles
{ // extra parenthesis to control destruction of the writer object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<float> zCoordWriter(outFile, "zCoord", 1, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter zCoordWriter(outFile, "zCoord", 1, Dune::Indent(4), VTK::Precision::float32);
+#endif
for (size_t i=0; i<zCoord_.size(); i++)
zCoordWriter.write(zCoord_[i]);
}
@@ -132,7 +168,11 @@ namespace Dune {
// The three director fields
for (size_t i=0; i<3; i++)
{
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<float> directorWriter(outFile, "director" + std::to_string(i), 3, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter directorWriter(outFile, "director" + std::to_string(i), 3, Dune::Indent(4), VTK::Precision::float32);
+#endif
for (size_t j=0; j<directors_[i].size(); j++)
for (int k=0; k<3; k++)
directorWriter.write(directors_[i][j][k]);
@@ -148,7 +188,11 @@ namespace Dune {
{
outFile << " <CellData>" << std::endl;
{ // extra parenthesis to control destruction of the writer object
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
Dune::VTK::AsciiDataArrayWriter<float> cellDataWriter(outFile, "mycelldata", 1, Dune::Indent(4));
+#else
+ Dune::VTK::AsciiDataArrayWriter cellDataWriter(outFile, "mycelldata", 1, Dune::Indent(4), VTK::Precision::float32);
+#endif
for (size_t i=0; i<cellData_.size(); i++)
cellDataWriter.write(cellData_[i]);
}
diff --git a/dune/gfe/vtkreader.hh b/dune/gfe/vtkreader.hh
index 566deead39a87ad62ed886dadd5c153921f5ab04..1279a655e0233f591cd881ff1aaf0e04bb903447 100644
--- a/dune/gfe/vtkreader.hh
+++ b/dune/gfe/vtkreader.hh
@@ -35,12 +35,9 @@ public:
factory.insertVertex(pos);
}
- Dune::GeometryType triangle;
- triangle.makeTriangle();
-
for (size_t i=0; i<vtkFile.cellConnectivity_.size(); i+=3)
{
- factory.insertElement(triangle, {vtkFile.cellConnectivity_[i],
+ factory.insertElement(Dune::GeometryTypes::triangle, {vtkFile.cellConnectivity_[i],
vtkFile.cellConnectivity_[i+1],
vtkFile.cellConnectivity_[i+2]});
@@ -51,4 +48,4 @@ public:
};
-#endif
\ No newline at end of file
+#endif
diff --git a/problems/cantilever-dirichlet-values.py b/problems/cantilever-dirichlet-values.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d1d60a5c0b33499677be2da2c10d2e6ea91b17c
--- /dev/null
+++ b/problems/cantilever-dirichlet-values.py
@@ -0,0 +1,17 @@
+import math
+
+class DirichletValues:
+ def __init__(self, homotopyParameter):
+ self.homotopyParameter = homotopyParameter
+
+ def deformation(self, x):
+ # Dirichlet b.c. simply clamp the shell in the reference configuration
+ out = [x[0], x[1], 0]
+
+ return out
+
+
+ def orientation(self, x):
+ rotation = [[1,0,0], [0, 1, 0], [0, 0, 1]]
+ return rotation
+
diff --git a/problems/cosserat-continuum-cantilever.parset b/problems/cosserat-continuum-cantilever.parset
new file mode 100644
index 0000000000000000000000000000000000000000..4244121458e5add3ca49043149afe4c7340eb6f8
--- /dev/null
+++ b/problems/cosserat-continuum-cantilever.parset
@@ -0,0 +1,110 @@
+#############################################
+# Grid parameters
+#############################################
+
+structuredGrid = true
+
+# bounding box
+lower = 0 0
+upper = 100 10
+
+elements = 10 1
+
+# Number of grid levels
+numLevels = 1
+
+#############################################
+# Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-3
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 200
+
+trustRegionScaling = 1 1 1 0.01 0.01 0.01
+
+# Initial trust-region radius
+initialTrustRegionRadius = 3.125
+
+# Number of multigrid iterations per trust-region step
+numIt = 400
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 1
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-5
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+# Measure convergence
+instrumented = 0
+
+############################
+# Material parameters
+############################
+
+
+## For the Wriggers L-shape example
+[materialParameters]
+
+# shell thickness
+thickness = 0.6
+
+# Lame parameters
+# corresponds to E = 71240 N/mm^2, nu=0.31
+# However, we use units N/m^2
+mu = 2.7191e+4
+lambda = 4.4364e+4
+
+# Cosserat couple modulus
+mu_c = 0
+
+# Length scale parameter
+L_c = 0.6
+
+# Curvature exponent
+q = 2
+
+# Shear correction factor
+kappa = 1
+
+[]
+
+#############################################
+# Boundary values
+#############################################
+
+problem = cantilever
+
+### Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+### Python predicate specifying all Neumann grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 99.99"
+
+### Neumann values
+neumannValues = 0 0 3
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], 0]"
+
+#startFromFile = yes
+#initialIterateFilename = initial_iterate.vtu
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d0e141f0357ba9068d2f0b7a1c355d71c58392d7..45a09e9929fa2e21f5abe2511d677218954de0e0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,16 +2,13 @@ set(programs compute-disc-error
cosserat-continuum
gradient-flow
harmonicmaps
- mixed-cosserat-continuum
rod-eoc)
# rodobstacle)
foreach(_program ${programs})
add_executable(${_program} ${_program}.cc)
- add_dune_adolc_flags(${_program})
- add_dune_ipopt_flags(${_program})
+ target_link_dune_default_libraries(${_program})
add_dune_pythonlibs_flags(${_program})
- add_dune_ug_flags(${_program})
- add_dune_mpi_flags(${_program})
+ target_link_libraries(${_program} tinyxml2)
# target_compile_options(${_program} PRIVATE "-fpermissive")
endforeach()
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index 14b9551f3514ea19cff4a9ca4b97eb3c71b04dd4..97cabf4093320723277103a98e7a3bd3fc21e26c 100644
--- a/src/compute-disc-error.cc
+++ b/src/compute-disc-error.cc
@@ -6,15 +6,23 @@
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/uggrid.hh>
+#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#if HAVE_DUNE_FOAMGRID
+#include <dune/foamgrid/foamgrid.hh>
+#else
+#include <dune/grid/onedgrid.hh>
+#endif
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/matrix-vector/genericvectortools.hh>
#include <dune/fufem/discretizationerror.hh>
#include <dune/fufem/dunepython.hh>
+#include <dune/fufem/makesphere.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/unitvector.hh>
@@ -29,6 +37,116 @@ const int dimworld = 2;
using namespace Dune;
+/** \brief Check whether given local coordinates are contained in the reference element
+ \todo This method exists in the Dune grid interface! But we need the eps.
+*/
+static bool checkInside(const Dune::GeometryType& type,
+ const Dune::FieldVector<double, dim> &loc,
+ double eps)
+{
+ switch (type.dim())
+ {
+ case 0: // vertex
+ return false;
+
+ case 1: // line
+ return -eps <= loc[0] && loc[0] <= 1+eps;
+
+ case 2:
+
+ if (type.isSimplex()) {
+ return -eps <= loc[0] && -eps <= loc[1] && (loc[0]+loc[1])<=1+eps;
+ } else if (type.isCube()) {
+ return -eps <= loc[0] && loc[0] <= 1+eps
+ && -eps <= loc[1] && loc[1] <= 1+eps;
+ } else
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+
+ case 3:
+ if (type.isSimplex()) {
+ return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+ && (loc[0]+loc[1]+loc[2]) <= 1+eps;
+ } else if (type.isPyramid()) {
+ return -eps <= loc[0] && -eps <= loc[1] && -eps <= loc[2]
+ && (loc[0]+loc[2]) <= 1+eps
+ && (loc[1]+loc[2]) <= 1+eps;
+ } else if (type.isPrism()) {
+ return -eps <= loc[0] && -eps <= loc[1]
+ && (loc[0]+loc[1])<= 1+eps
+ && -eps <= loc[2] && loc[2] <= 1+eps;
+ } else if (type.isCube()) {
+ return -eps <= loc[0] && loc[0] <= 1+eps
+ && -eps <= loc[1] && loc[1] <= 1+eps
+ && -eps <= loc[2] && loc[2] <= 1+eps;
+ }else
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+
+ default:
+ DUNE_THROW(Dune::GridError, "checkInside(): ERROR: Unknown type " << type << " found!");
+ }
+
+}
+
+/** \param element An element of the source grid
+ */
+template <class GridType>
+auto findSupportingElement(const GridType& sourceGrid,
+ const GridType& targetGrid,
+ typename GridType::template Codim<0>::Entity element,
+ FieldVector<double,dim> pos)
+{
+ while (element.level() != 0)
+ {
+ pos = element.geometryInFather().global(pos);
+ element = element.father();
+
+ assert(checkInside(element.type(), pos, 1e-7));
+ }
+
+ //////////////////////////////////////////////////////////////////////
+ // Find the corresponding coarse grid element on the adaptive grid.
+ // This is a linear algorithm, but we expect the coarse grid to be small.
+ //////////////////////////////////////////////////////////////////////
+ LevelMultipleCodimMultipleGeomTypeMapper<GridType> sourceP0Mapper (sourceGrid, 0, mcmgElementLayout());
+ LevelMultipleCodimMultipleGeomTypeMapper<GridType> targetP0Mapper(targetGrid, 0, mcmgElementLayout());
+ const auto coarseIndex = sourceP0Mapper.index(element);
+
+ auto targetLevelView = targetGrid.levelGridView(0);
+
+ for (auto&& targetElement : elements(targetLevelView))
+ if (targetP0Mapper.index(targetElement) == coarseIndex)
+ {
+ element = std::move(targetElement);
+ break;
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ // Find a corresponding point (not necessarily vertex) on the leaf level
+ // of the adaptive grid.
+ //////////////////////////////////////////////////////////////////////////
+
+ while (!element.isLeaf())
+ {
+ FieldVector<double,dim> childPos;
+ assert(checkInside(element.type(), pos, 1e-7));
+
+ auto hIt = element.hbegin(element.level()+1);
+ auto hEndIt = element.hend(element.level()+1);
+
+ for (; hIt!=hEndIt; ++hIt)
+ {
+ childPos = hIt->geometryInFather().local(pos);
+ if (checkInside(hIt->type(), childPos, 1e-7))
+ break;
+ }
+
+ element = *hIt;
+ pos = childPos;
+ }
+
+ return std::tie(element, pos);
+}
+
template <class GridView, int order, class TargetSpace>
void measureDiscreteEOC(const GridView gridView,
const GridView referenceGridView,
@@ -40,7 +158,7 @@ void measureDiscreteEOC(const GridView gridView,
// Construct the scalar function space bases corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
- typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
FEBasis referenceFEBasis(referenceGridView);
@@ -121,23 +239,54 @@ void measureDiscreteEOC(const GridView gridView,
auto element = hierarchicSearch.findEntity(globalPos);
auto localPos = element.geometry().local(globalPos);
- auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
+ auto refValue = referenceSolution(rElement, qp.position());
+ auto numValue = numericalSolution(element, localPos);
+ auto diff = refValue - numValue;
assert(diff.size()==7);
+ // Compute error of the displacement degrees of freedom
for (int i=0; i<3; i++)
deformationL2ErrorSquared += integrationElement * qp.weight() * diff[i] * diff[i];
- for (int i=3; i<7; i++)
- orientationL2ErrorSquared += integrationElement * qp.weight() * diff[i] * diff[i];
-
auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
for (int i=0; i<3; i++)
deformationH1ErrorSquared += integrationElement * qp.weight() * derDiff[i].two_norm2();
- for (int i=3; i<7; i++)
- orientationH1ErrorSquared += integrationElement * qp.weight() * derDiff[i].two_norm2();
+ // Compute error of the orientation degrees of freedom
+ // We need to transform from quaternion coordinates to matrix coordinates first.
+ FieldMatrix<double,3,3> referenceValue, numericalValue;
+ Rotation<double,3> referenceRotation(std::array<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
+ referenceRotation.matrix(referenceValue);
+
+ Rotation<double,3> numericalRotation(std::array<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
+ numericalRotation.matrix(numericalValue);
+ auto orientationDiff = referenceValue - numericalValue;
+
+ orientationL2ErrorSquared += integrationElement * qp.weight() * orientationDiff.frobenius_norm2();
+
+ auto referenceDerQuat = referenceSolution.derivative(rElement, qp.position());
+ auto numericalDerQuat = numericalSolution.derivative(element, localPos);
+
+ // Transform to matrix coordinates
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{refValue[3], refValue[4], refValue[5], refValue[6]});
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(FieldVector<double,4>{numValue[3], numValue[4], numValue[5], numValue[6]});
+
+ Tensor3<double,3,3,dim> refDerivative(0);
+ Tensor3<double,3,3,dim> numDerivative(0);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<dim; k++)
+ for (int l=0; l<4; l++)
+ {
+ refDerivative[i][j][k] += derivativeQuaternionToMatrixRef[i][j][l] * referenceDerQuat[3+l][k];
+ numDerivative[i][j][k] += derivativeQuaternionToMatrixNum[i][j][l] * numericalDerQuat[3+l][k];
+ }
+
+ auto derOrientationDiff = refDerivative - numDerivative; // compute the difference
+ orientationH1ErrorSquared += derOrientationDiff.frobenius_norm2() * qp.weight() * integrationElement;
}
}
@@ -152,6 +301,69 @@ void measureDiscreteEOC(const GridView gridView,
<< "H^1 error orientation: " << std::sqrt(orientationH1ErrorSquared)
<< std::endl;
}
+ else if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
+ {
+ double l2ErrorSquared = 0;
+ double h1ErrorSquared = 0;
+
+ for (const auto& rElement : elements(referenceGridView))
+ {
+ const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
+
+ for (const auto& qp : quadRule)
+ {
+ auto integrationElement = rElement.geometry().integrationElement(qp.position());
+
+ auto globalPos = rElement.geometry().global(qp.position());
+
+ auto element = hierarchicSearch.findEntity(globalPos);
+ auto localPos = element.geometry().local(globalPos);
+
+ FieldMatrix<double,3,3> referenceValue, numericalValue;
+ auto refValue = referenceSolution(rElement, qp.position());
+ Rotation<double,3> referenceRotation(refValue);
+ referenceRotation.matrix(referenceValue);
+
+ auto numValue = numericalSolution(element, localPos);
+ Rotation<double,3> numericalRotation(numValue);
+ numericalRotation.matrix(numericalValue);
+
+ auto diff = referenceValue - numericalValue;
+
+ l2ErrorSquared += integrationElement * qp.weight() * diff.frobenius_norm2();
+
+ auto referenceDerQuat = referenceSolution.derivative(rElement, qp.position());
+ auto numericalDerQuat = numericalSolution.derivative(element, localPos);
+
+ // Transform to matrix coordinates
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
+
+ Tensor3<double,3,3,dim> refDerivative(0);
+ Tensor3<double,3,3,dim> numDerivative(0);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<dim; k++)
+ for (int l=0; l<4; l++)
+ {
+ refDerivative[i][j][k] += derivativeQuaternionToMatrixRef[i][j][l] * referenceDerQuat[l][k];
+ numDerivative[i][j][k] += derivativeQuaternionToMatrixNum[i][j][l] * numericalDerQuat[l][k];
+ }
+
+ auto derDiff = refDerivative - numDerivative; // compute the difference
+ h1ErrorSquared += derDiff.frobenius_norm2() * qp.weight() * integrationElement;
+
+ }
+ }
+
+ std::cout << "levels: " << gridView.grid().maxLevel()+1
+ << " "
+ << "L^2 error: " << std::sqrt(l2ErrorSquared)
+ << " "
+ << "h^1 error: " << std::sqrt(h1ErrorSquared)
+ << std::endl;
+ }
else
{
double l2ErrorSquared = 0;
@@ -165,10 +377,13 @@ void measureDiscreteEOC(const GridView gridView,
{
auto integrationElement = rElement.geometry().integrationElement(qp.position());
- auto globalPos = rElement.geometry().global(qp.position());
-
- auto element = hierarchicSearch.findEntity(globalPos);
- auto localPos = element.geometry().local(globalPos);
+ // Given a point with local coordinates qp.position() on element rElement of the reference grid
+ // find the element and local coordinates on the grid of the numerical simulation
+ auto supportingElement = findSupportingElement(referenceGridView.grid(),
+ gridView.grid(),
+ rElement, qp.position());
+ auto element = std::get<0>(supportingElement);
+ auto localPos = std::get<1>(supportingElement);
auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
@@ -200,7 +415,7 @@ void measureAnalyticalEOC(const GridView gridView,
// Construct the scalar function space bases corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
- typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
//////////////////////////////////////////////////////////////////////////////////
@@ -229,7 +444,7 @@ void measureAnalyticalEOC(const GridView gridView,
/////////////////////////////////////////////////////////////////
// Read reference solution and its derivative into a PythonFunction
- typedef VirtualDifferentiableFunction<FieldVector<double, dim>, typename TargetSpace::CoordinateType> FBase;
+ typedef VirtualDifferentiableFunction<FieldVector<double, dimworld>, typename TargetSpace::CoordinateType> FBase;
Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
@@ -250,22 +465,116 @@ void measureAnalyticalEOC(const GridView gridView,
// QuadratureRule for the integral of the L^2 error
QuadratureRuleKey quadKey(dim,6);
- // Compute the embedded L^2 error
- double l2Error = DiscretizationError<GridView>::computeL2Error(numericalSolution.get(),
- referenceSolution.get(),
- quadKey);
+ // Compute errors in the L2 norm and the h1 seminorm.
+ // SO(3)-valued maps need special treatment, because they are stored as quaternions,
+ // but the errors need to be computed in matrix space.
+ if constexpr (std::is_same<TargetSpace,Rotation<double,3> >::value)
+ {
+ constexpr int blocksize = TargetSpace::CoordinateType::dimension;
- // Compute the embedded H^1 error
- double h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
- numericalSolution.get(),
- referenceSolution.get(),
- quadKey);
+ // The error to be computed
+ double l2ErrorSquared = 0;
+ double h1ErrorSquared = 0;
- std::cout << "elements: " << gridView.size(0)
- << " "
- << "L^2 error: " << l2Error
- << " ";
- std::cout << "h^1 error: " << std::sqrt(h1Error) << std::endl;
+ for (auto&& element : elements(gridView))
+ {
+ // Get quadrature formula
+ quadKey.setGeometryType(element.type());
+ const auto& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
+
+ for (auto quadPoint : quad)
+ {
+ auto quadPos = quadPoint.position();
+ const auto integrationElement = element.geometry().integrationElement(quadPos);
+ const auto weight = quadPoint.weight();
+
+ // Evaluate function a
+ FieldVector<double,blocksize> numValue;
+ numericalSolution.get()->evaluateLocal(element, quadPos,numValue);
+
+ // Evaluate function b. If it is a grid function use that to speed up the evaluation
+ FieldVector<double,blocksize> refValue;
+ if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >>(referenceSolution))
+ std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >>(referenceSolution)->evaluateLocal(element,
+ quadPos,
+ refValue
+ );
+ else
+ referenceSolution->evaluate(element.geometry().global(quadPos), refValue);
+
+ // Get error in matrix space
+ Rotation<double,3> numRotation(numValue);
+ FieldMatrix<double,3,3> numValueMatrix;
+ numRotation.matrix(numValueMatrix);
+
+ Rotation<double,3> refRotation(refValue);
+ FieldMatrix<double,3,3> refValueMatrix;
+ refRotation.matrix(refValueMatrix);
+
+ // Evaluate derivatives in quaternion space
+ FieldMatrix<double,blocksize,dimworld> num_di;
+ FieldMatrix<double,blocksize,dimworld> ref_di;
+
+ if (dynamic_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >*>(numericalSolution.get()))
+ dynamic_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >*>(numericalSolution.get())->evaluateDerivativeLocal(element,
+ quadPos,
+ num_di);
+ else
+ numericalSolution->evaluateDerivative(element.geometry().global(quadPos), num_di);
+
+ if (std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >>(referenceSolution))
+ std::dynamic_pointer_cast<const VirtualGridViewFunction<GridView,FieldVector<double,blocksize> >>(referenceSolution)->evaluateDerivativeLocal(element,
+ quadPos,
+ ref_di);
+ else
+ referenceSolution->evaluateDerivative(element.geometry().global(quadPos), ref_di);
+
+ // Transform into matrix space
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixNum = Rotation<double,3>::derivativeOfQuaternionToMatrix(numValue);
+ Tensor3<double,3,3,4> derivativeQuaternionToMatrixRef = Rotation<double,3>::derivativeOfQuaternionToMatrix(refValue);
+
+ Tensor3<double,3,3,dim> numDerivative(0);
+ Tensor3<double,3,3,dim> refDerivative(0);
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<dim; k++)
+ for (int l=0; l<blocksize; l++)
+ {
+ numDerivative[i][j][k] += derivativeQuaternionToMatrixNum[i][j][l] * num_di[l][k];
+ refDerivative[i][j][k] += derivativeQuaternionToMatrixRef[i][j][l] * ref_di[l][k];
+ }
+
+ // integrate error
+ l2ErrorSquared += (numValueMatrix - refValueMatrix).frobenius_norm2() * weight * integrationElement;
+ auto diff = numDerivative - refDerivative;
+ h1ErrorSquared += diff.frobenius_norm2() * weight * integrationElement;
+ }
+ }
+
+ std::cout << "elements: " << gridView.size(0)
+ << " "
+ << "L^2 error: " << std::sqrt(l2ErrorSquared)
+ << " ";
+ std::cout << "h^1 error: " << std::sqrt(h1ErrorSquared) << std::endl;
+ }
+ else
+ {
+ auto l2Error = DiscretizationError<GridView>::computeL2Error(numericalSolution.get(),
+ referenceSolution.get(),
+ quadKey);
+
+ auto h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
+ numericalSolution.get(),
+ referenceSolution.get(),
+ quadKey);
+
+ std::cout << "elements: " << gridView.size(0)
+ << " "
+ << "L^2 error: " << l2Error
+ << " ";
+ std::cout << "h^1 error: " << std::sqrt(h1Error) << std::endl;
+ }
}
template <class GridType, class TargetSpace>
@@ -277,6 +586,8 @@ void measureEOC(const std::shared_ptr<GridType> grid,
if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
{
+ referenceGrid->globalRefine(parameterSet.get<int>("numReferenceLevels")-1);
+
switch (order)
{
case 1:
@@ -294,6 +605,7 @@ void measureEOC(const std::shared_ptr<GridType> grid,
default:
DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
}
+ return; // Success
}
if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
@@ -315,7 +627,10 @@ void measureEOC(const std::shared_ptr<GridType> grid,
default:
DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
}
+ return; // Success
}
+
+ DUNE_THROW(NotImplemented, "Unknown discretization error mode encountered!");
}
int main (int argc, char *argv[]) try
@@ -347,7 +662,12 @@ int main (int argc, char *argv[]) try
/////////////////////////////////////////
// Create the grids
/////////////////////////////////////////
- typedef UGGrid<dim> GridType;
+#if HAVE_DUNE_FOAMGRID
+ typedef std::conditional<dim==1 or dim!=dimworld,FoamGrid<dim,dimworld>,UGGrid<dim> >::type GridType;
+#else
+ static_assert(dim==dimworld, "You need to have dune-foamgrid installed for dim != dimworld!");
+ typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
+#endif
const int numLevels = parameterSet.get<int>("numLevels");
@@ -356,29 +676,41 @@ int main (int argc, char *argv[]) try
FieldVector<double,dimworld> lower(0), upper(1);
std::string structuredGridType = parameterSet["structuredGrid"];
- if (structuredGridType != "false" )
+ if (structuredGridType == "sphere")
+ {
+ if constexpr (dimworld==3)
+ {
+ grid = makeSphereOnOctahedron<GridType>({0,0,0},1);
+ referenceGrid = makeSphereOnOctahedron<GridType>({0,0,0},1);
+ } else
+ DUNE_THROW(Exception, "Dimworld must be 3 to create a sphere grid!");
+ }
+ else if (structuredGridType == "simplex" || structuredGridType == "cube")
{
lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
auto elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
if (structuredGridType == "simplex")
- grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements);
- else if (structuredGridType == "cube")
- grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
- else
+ {
+ grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements);
+ referenceGrid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements);
+ } else if (structuredGridType == "cube")
+ {
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ referenceGrid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ } else
DUNE_THROW(Exception, "Unknown structured grid type '" << structuredGridType << "' found!");
}
else
{
std::string path = parameterSet.get<std::string>("path");
std::string gridFile = parameterSet.get<std::string>("gridFile");
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
referenceGrid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
}
grid->globalRefine(numLevels-1);
- referenceGrid->globalRefine(parameterSet.get<int>("numReferenceLevels")-1);
// Do the actual measurement
const int targetDim = parameterSet.get<int>("targetDim");
@@ -486,7 +818,7 @@ int main (int argc, char *argv[]) try
return 0;
}
-catch (Exception e)
+catch (Exception& e)
{
std::cout << e << std::endl;
return 1;
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 95c87570a740a9e0c545281a52180cd5d9d5391c..9029c25248891d0c0e8614d1eca0ad72de6c7911 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -1,6 +1,7 @@
#include <config.h>
#include <fenv.h>
+#include <array>
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
@@ -10,9 +11,11 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+#include <dune/common/typetraits.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
+#include <dune/common/tuplevector.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
@@ -23,7 +26,9 @@
#include <dune/foamgrid/foamgrid.hh>
#endif
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
@@ -36,6 +41,7 @@
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/cosseratenergystiffness.hh>
#include <dune/gfe/nonplanarcosseratshellenergy.hh>
#include <dune/gfe/cosseratvtkwriter.hh>
@@ -45,19 +51,30 @@
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/vertexnormals.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
+#include <dune/gfe/mixedgfeassembler.hh>
+#include <dune/gfe/mixedriemanniantrsolver.hh>
// grid dimension
const int dim = 2;
const int dimworld = 2;
-// Order of the approximation space
-const int order = 2;
+//#define MIXED_SPACE
+
+// Order of the approximation space for the displacement
+const int displacementOrder = 2;
+
+// Order of the approximation space for the microrotations
+const int rotationOrder = 2;
+
+#ifndef MIXED_SPACE
+static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
// Image space of the geodesic fe functions
typedef RigidBodyMotion<double,3> TargetSpace;
// Tangent vector of the image space
const int blocksize = TargetSpace::TangentVector::dimension;
+#endif
using namespace Dune;
@@ -116,7 +133,13 @@ int main (int argc, char *argv[]) try
<< std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems/')"
<< std::endl;
+ using namespace TypeTree::Indices;
+#ifdef MIXED_SPACE
+ using SolutionType = TupleVector<std::vector<RealTuple<double,3> >,
+ std::vector<Rotation<double,3> > >;
+#else
typedef std::vector<TargetSpace> SolutionType;
+#endif
// parse data file
ParameterTree parameterSet;
@@ -157,12 +180,13 @@ int main (int argc, char *argv[]) try
FieldVector<double,dimworld> lower(0), upper(1);
- if (parameterSet.get<bool>("structuredGrid")) {
+ std::string structuredGridType = parameterSet["structuredGrid"];
+ if (structuredGridType == "cube") {
lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
- array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
+ auto elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
} else {
@@ -191,11 +215,35 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView, order> FEBasis;
+ using namespace Dune::Functions::BasisFactory;
+#ifdef MIXED_SPACE
+ auto compositeBasis = makeBasis(
+ gridView,
+ composite(
+ lagrange<displacementOrder>(),
+ lagrange<rotationOrder>()
+ )
+ );
+
+ typedef Dune::Functions::LagrangeBasis<GridView,displacementOrder> DeformationFEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView,rotationOrder> OrientationFEBasis;
+
+ DeformationFEBasis deformationFEBasis(gridView);
+ OrientationFEBasis orientationFEBasis(gridView);
+
+ // Construct fufem-style function space bases to ease the transition to dune-functions
+ typedef DuneFunctionsBasis<DeformationFEBasis> FufemDeformationFEBasis;
+ FufemDeformationFEBasis fufemDeformationFEBasis(deformationFEBasis);
+
+ typedef DuneFunctionsBasis<OrientationFEBasis> FufemOrientationFEBasis;
+ FufemOrientationFEBasis fufemOrientationFEBasis(orientationFEBasis);
+#else
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, displacementOrder> FEBasis;
FEBasis feBasis(gridView);
typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
FufemFEBasis fufemFeBasis(feBasis);
+#endif
// /////////////////////////////////////////
// Read Dirichlet values
@@ -232,23 +280,71 @@ int main (int argc, char *argv[]) try
if (mpiHelper.rank()==0)
std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+#ifdef MIXED_SPACE
+ BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemDeformationFEBasis,deformationDirichletNodes);
- BitSetVector<1> dirichletNodes(feBasis.indexSet().size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemFeBasis,dirichletNodes);
+ BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary,fufemDeformationFEBasis,neumannNodes);
+
+ BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
+ for (size_t i=0; i<deformationFEBasis.size(); i++)
+ if (deformationDirichletNodes[i][0])
+ for (int j=0; j<3; j++)
+ deformationDirichletDofs[i][j] = true;
+
+ BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemOrientationFEBasis,orientationDirichletNodes);
+
+ BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
+ for (size_t i=0; i<orientationFEBasis.size(); i++)
+ if (orientationDirichletNodes[i][0])
+ for (int j=0; j<3; j++)
+ orientationDirichletDofs[i][j] = true;
+#else
+ BitSetVector<1> dirichletNodes(feBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
- BitSetVector<1> neumannNodes(feBasis.indexSet().size(), false);
- constructBoundaryDofs(neumannBoundary,fufemFeBasis,neumannNodes);
+ BitSetVector<1> neumannNodes(feBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
- BitSetVector<blocksize> dirichletDofs(feBasis.indexSet().size(), false);
- for (size_t i=0; i<feBasis.indexSet().size(); i++)
+ BitSetVector<blocksize> dirichletDofs(feBasis.size(), false);
+ for (size_t i=0; i<feBasis.size(); i++)
if (dirichletNodes[i][0])
for (int j=0; j<5; j++)
dirichletDofs[i][j] = true;
+#endif
// //////////////////////////
// Initial iterate
// //////////////////////////
+#ifdef MIXED_SPACE
+ SolutionType x;
+
+ x[_0].resize(deformationFEBasis.size());
+
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
+ std::vector<FieldVector<double,3> > v;
+ ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
+
+ for (size_t i=0; i<x[_0].size(); i++)
+ x[_0][i] = v[i];
+
+ x[_1].resize(orientationFEBasis.size());
+#if 0
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
+
+ std::vector<FieldVector<double,3> > v;
+ Functions::interpolate(feBasis, v, pythonInitialDeformation);
+
+ for (size_t i=0; i<x.size(); i++)
+ xDisp[i] = v[i];
+#endif
+#else
SolutionType x(feBasis.size());
if (parameterSet.hasKey("startFromFile"))
@@ -256,7 +352,7 @@ int main (int argc, char *argv[]) try
std::shared_ptr<GridType> initialIterateGrid;
if (parameterSet.get<bool>("structuredGrid"))
{
- std::array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
+ std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
initialIterateGrid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
}
else
@@ -270,13 +366,19 @@ int main (int argc, char *argv[]) try
SolutionType initialIterate;
GFE::CosseratVTKReader::read(initialIterate, parameterSet.get<std::string>("initialIterateFilename"));
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView, 2> InitialBasis;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, 2> InitialBasis;
InitialBasis initialBasis(initialIterateGrid->leafGridView());
- GFE::EmbeddedGlobalGFEFunction<InitialBasis,TargetSpace> initialFunction(initialBasis,initialIterate);
+#ifdef PROJECTED_INTERPOLATION
+ using LocalInterpolationRule = LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+#else
+ using LocalInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+#endif
+ GFE::EmbeddedGlobalGFEFunction<InitialBasis,LocalInterpolationRule,TargetSpace> initialFunction(initialBasis,initialIterate);
std::vector<FieldVector<double,7> > v;
Dune::Functions::interpolate(feBasis,v,initialFunction);
+ DUNE_THROW(NotImplemented, "Replace scalar basis by power basis!");
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(v[i]);
@@ -291,13 +393,20 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); i++)
x[i].r = v[i];
}
+#endif
////////////////////////////////////////////////////////
// Main homotopy loop
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
+#ifdef MIXED_SPACE
+ CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
+ orientationFEBasis,x[_1],
+ resultPath + "mixed-cosserat_homotopy_0");
+#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_0");
+#endif
for (int i=0; i<numHomotopySteps; i++) {
@@ -316,17 +425,32 @@ int main (int argc, char *argv[]) try
neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
homotopyParameter);
- shared_ptr<VolumeLoad> volumeLoad;
- if (parameterSet.hasKey("volumeLoad"))
- volumeLoad = make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
+ shared_ptr<VolumeLoad> volumeLoad;
+ if (parameterSet.hasKey("volumeLoad"))
+ volumeLoad = make_shared<VolumeLoad>(parameterSet.get<FieldVector<double,3> >("volumeLoad"),
homotopyParameter);
- if (mpiHelper.rank() == 0) {
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
- }
+ if (mpiHelper.rank() == 0) {
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
+ }
// Assembler using ADOL-C
+#ifdef MIXED_SPACE
+ CosseratEnergyLocalStiffness<decltype(compositeBasis),
+ 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
+ &neumannBoundary,
+ neumannFunction,
+ nullptr);
+
+ MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
+ RealTuple<double,3>,
+ Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
+
+ MixedGFEAssembler<decltype(compositeBasis),
+ RealTuple<double,3>,
+ Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
+#else
using LocalEnergyBase = LocalGeodesicFEStiffness<FEBasis,RigidBodyMotion<adouble,3> >;
std::shared_ptr<LocalEnergyBase> cosseratEnergyADOLCLocalStiffness;
@@ -352,16 +476,33 @@ int main (int argc, char *argv[]) try
TargetSpace> localGFEADOLCStiffness(cosseratEnergyADOLCLocalStiffness.get());
GeodesicFEAssembler<FEBasis,TargetSpace> assembler(gridView, &localGFEADOLCStiffness);
+#endif
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
+#ifdef MIXED_SPACE
+ MixedRiemannianTrustRegionSolver<GridType,
+ decltype(compositeBasis),
+ DeformationFEBasis, RealTuple<double,3>,
+ OrientationFEBasis, Rotation<double,3> > solver;
+#else
RiemannianTrustRegionSolver<FEBasis,TargetSpace> solver;
+#endif
solver.setup(*grid,
&assembler,
+#ifdef MIXED_SPACE
+ deformationFEBasis,
+ orientationFEBasis,
+#endif
x,
+#ifdef MIXED_SPACE
+ deformationDirichletDofs,
+ orientationDirichletDofs,
+#else
dirichletDofs,
+#endif
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
@@ -390,9 +531,21 @@ int main (int argc, char *argv[]) try
PythonFunction<FieldVector<double,dimworld>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemFeBasis, ddV, deformationDirichletValues, dirichletDofs);
-
std::vector<FieldMatrix<double,3,3> > dOV;
+
+#ifdef MIXED_SPACE
+ ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+ ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+
+ for (size_t j=0; j<x[_0].size(); j++)
+ if (deformationDirichletNodes[j][0])
+ x[_0][j] = ddV[j];
+
+ for (size_t j=0; j<x[_1].size(); j++)
+ if (orientationDirichletNodes[j][0])
+ x[_1][j].set(dOV[j]);
+#else
+ ::Functions::interpolate(fufemFeBasis, ddV, deformationDirichletValues, dirichletDofs);
::Functions::interpolate(fufemFeBasis, dOV, orientationDirichletValues, dirichletDofs);
for (size_t j=0; j<x.size(); j++)
@@ -401,6 +554,7 @@ int main (int argc, char *argv[]) try
x[j].r = ddV[j];
x[j].q.set(dOV[j]);
}
+#endif
// /////////////////////////////////////////////////////
// Solve!
@@ -414,7 +568,13 @@ int main (int argc, char *argv[]) try
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
+#ifdef MIXED_SPACE
+ CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
+ orientationFEBasis,x[_1],
+ resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
+#else
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "cosserat_homotopy_" + iAsAscii.str());
+#endif
}
@@ -422,6 +582,7 @@ int main (int argc, char *argv[]) try
// Output result
// //////////////////////////////
+#ifndef MIXED_SPACE
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
// This data may be used by other applications measuring the discretization error
BlockVector<TargetSpace::CoordinateType> xEmbedded(x.size());
@@ -429,26 +590,31 @@ int main (int argc, char *argv[]) try
xEmbedded[i] = x[i].globalCoordinates();
std::ofstream outFile("cosserat-continuum-result-" + std::to_string(numLevels) + ".data", std::ios_base::binary);
- GenericVector::writeBinary(outFile, xEmbedded);
+ MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
+#endif
// finally: compute the average deformation of the Neumann boundary
// That is what we need for the locking tests
FieldVector<double,3> averageDef(0);
+#ifdef MIXED_SPACE
+ for (size_t i=0; i<x[_0].size(); i++)
+ if (neumannNodes[i][0])
+ averageDef += x[_0][i].globalCoordinates();
+#else
for (size_t i=0; i<x.size(); i++)
if (neumannNodes[i][0])
averageDef += x[i].r;
+#endif
averageDef /= neumannNodes.count();
- if (mpiHelper.rank()==0)
+ if (mpiHelper.rank()==0 and parameterSet.hasKey("neumannValues"))
{
std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
<< ", average deflection: " << averageDef << std::endl;
}
- // //////////////////////////////
- } catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }
+} catch (Exception& e)
+{
+ std::cout << e.what() << std::endl;
+}
diff --git a/src/gradient-flow.cc b/src/gradient-flow.cc
index 3d49ced53f294554a93e8f5782f2ec2a72c491e2..d74c21a2b6237013323989dfff01bf0b8e0fde5f 100644
--- a/src/gradient-flow.cc
+++ b/src/gradient-flow.cc
@@ -2,13 +2,14 @@
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
+#include <array>
+
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
-#include <array>
-
+#include <dune/common/typetraits.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
@@ -22,7 +23,7 @@
#include <dune/grid/io/file/vtk.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/basisinterpolator.hh>
@@ -141,7 +142,7 @@ int main (int argc, char *argv[]) try
// Construct the scalar function space basis corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView, order> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, order> FEBasis;
//typedef Dune::Functions::Periodic1DPQ1NodalBasis<typename GridType::LeafGridView> FEBasis;
FEBasis feBasis(grid->leafGridView());
@@ -201,7 +202,8 @@ int main (int argc, char *argv[]) try
auto l2DistanceSquaredEnergy = std::make_shared<L2DistanceSquaredEnergy<FEBasis, ATargetSpace> >();
std::vector<std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > > addends(2);
- addends[0] = std::make_shared<HarmonicEnergyLocalStiffness<FEBasis, ATargetSpace> >();
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ addends[0] = std::make_shared<HarmonicEnergyLocalStiffness<FEBasis, GeodesicInterpolationRule, ATargetSpace> >();
addends[1] = l2DistanceSquaredEnergy;
std::vector<double> weights = {1.0, 1.0/(2*timeStepSize)};
@@ -244,13 +246,13 @@ int main (int argc, char *argv[]) try
TypeTree::hybridTreePath(),
xEmbedded);
- SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),0);
+ SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),refinementLevels(0));
vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::scalar, xEmbedded[0].size()));
vtkWriter.write("gradientflow_result_0");
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
std::ofstream outFile("gradientflow_result_0.data", std::ios_base::binary);
- GenericVector::writeBinary(outFile, xEmbedded);
+ MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
///////////////////////////////////////////////////////
@@ -284,21 +286,21 @@ int main (int argc, char *argv[]) try
TypeTree::hybridTreePath(),
xEmbedded);
- SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),0);
+ SubsamplingVTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView(),refinementLevels(0));
vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::scalar, xEmbedded[0].size()));
vtkWriter.write("gradientflow_result_" + std::to_string(i+1));
// Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
std::ofstream outFile("gradientflow_result_" + std::to_string(i+1) + ".data", std::ios_base::binary);
- GenericVector::writeBinary(outFile, xEmbedded);
+ MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
}
return 0;
}
-catch (Exception e)
+catch (Exception& e)
{
- std::cout << e << std::endl;
+ std::cout << e.what() << std::endl;
return 1;
}
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 4717e6c181663cc0ffe9a3a2bbe25c3b305e0de5..00479169a3d0f26ac02e1d265e6d6c0222053401 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -1,6 +1,7 @@
#include <config.h>
#include <fenv.h>
+#include <array>
// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
@@ -8,30 +9,27 @@
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/typetraits.hh>
-namespace Dune {
- template <>
- struct IsNumber<adouble>
- {
- constexpr static bool value = true;
- };
-}
-
-#include <array>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/uggrid.hh>
-#include <dune/grid/onedgrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/io/file/vtk.hh>
+#if HAVE_DUNE_FOAMGRID
+#include <dune/foamgrid/foamgrid.hh>
+#else
+#include <dune/grid/onedgrid.hh>
+#endif
+
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/bsplinebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/fufem/boundarypatch.hh>
@@ -45,6 +43,7 @@ namespace Dune {
#include <dune/gfe/rotation.hh>
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/realtuple.hh>
+#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
@@ -56,10 +55,12 @@ namespace Dune {
// grid dimension
const int dim = 2;
+const int dimworld = 2;
// Image space of the geodesic fe functions
// typedef Rotation<double,2> TargetSpace;
// typedef Rotation<double,3> TargetSpace;
+// typedef RigidBodyMotion<double,3> TargetSpace;
// typedef UnitVector<double,2> TargetSpace;
typedef UnitVector<double,3> TargetSpace;
// typedef UnitVector<double,4> TargetSpace;
@@ -74,6 +75,51 @@ const int order = 1;
using namespace Dune;
+template <typename Writer, typename Basis, typename SolutionType>
+void fillVTKWriter(Writer& vtkWriter, const Basis& feBasis, const SolutionType& x, std::string filename)
+{
+ typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
+ EmbeddedVectorType xEmbedded(x.size());
+ for (size_t i=0; i<x.size(); i++)
+ xEmbedded[i] = x[i].globalCoordinates();
+
+ if constexpr (std::is_same<TargetSpace, Rotation<double,3> >::value)
+ {
+ std::array<BlockVector<FieldVector<double,3> >,3> director;
+ for (int i=0; i<3; i++)
+ director[i].resize(x.size());
+
+ for (size_t i=0; i<x.size(); i++)
+ {
+ FieldMatrix<double,3,3> m;
+ x[i].matrix(m);
+ director[0][i] = {m[0][0], m[1][0], m[2][0]};
+ director[1][i] = {m[0][1], m[1][1], m[2][1]};
+ director[2][i] = {m[0][2], m[1][2], m[2][2]};
+ }
+
+ auto dFunction0 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(feBasis,director[0]);
+ auto dFunction1 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(feBasis,director[1]);
+ auto dFunction2 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(feBasis,director[2]);
+
+ vtkWriter.addVertexData(dFunction0, VTK::FieldInfo("director0", VTK::FieldInfo::Type::vector, 3));
+ vtkWriter.addVertexData(dFunction1, VTK::FieldInfo("director1", VTK::FieldInfo::Type::vector, 3));
+ vtkWriter.addVertexData(dFunction2, VTK::FieldInfo("director2", VTK::FieldInfo::Type::vector, 3));
+
+ // Needs to be in this scope; otherwise the stack-allocated dFunction?-objects will get
+ // destructed before 'write' is called.
+ vtkWriter.write(filename);
+ }
+ else
+ {
+ auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,TypeTree::hybridTreePath(),xEmbedded);
+
+ vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::vector, xEmbedded[0].size()));
+
+ vtkWriter.write(filename);
+ }
+}
+
int main (int argc, char *argv[])
{
@@ -124,17 +170,22 @@ int main (int argc, char *argv[])
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
+#if HAVE_DUNE_FOAMGRID
+ typedef std::conditional<dim==1 or dim!=dimworld,FoamGrid<dim,dimworld>,UGGrid<dim> >::type GridType;
+#else
+ static_assert(dim==dimworld, "You need to have dune-foamgrid installed for dim != dimworld!");
typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
+#endif
shared_ptr<GridType> grid;
- FieldVector<double,dim> lower(0), upper(1);
+ FieldVector<double,dimworld> lower(0), upper(1);
std::array<unsigned int,dim> elements;
std::string structuredGridType = parameterSet["structuredGrid"];
if (structuredGridType != "false" ) {
- lower = parameterSet.get<FieldVector<double,dim> >("lower");
- upper = parameterSet.get<FieldVector<double,dim> >("upper");
+ lower = parameterSet.get<FieldVector<double,dimworld> >("lower");
+ upper = parameterSet.get<FieldVector<double,dimworld> >("upper");
elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
if (structuredGridType == "simplex")
@@ -161,7 +212,7 @@ int main (int argc, char *argv[])
using GridView = GridType::LeafGridView;
GridView gridView = grid->leafGridView();
#ifdef LAGRANGE
- typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
#else
typedef Dune::Functions::BSplineBasis<GridView> FEBasis;
@@ -180,7 +231,7 @@ int main (int argc, char *argv[])
// Make Python function that computes which vertices are on the Dirichlet boundary,
// based on the vertex positions.
std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+ PythonFunction<FieldVector<double,dimworld>, bool> pythonDirichletVertices(Python::evaluate(lambda));
for (auto&& vertex : vertices(gridView))
{
@@ -205,13 +256,22 @@ int main (int argc, char *argv[])
// Read initial iterate into a PythonFunction
Python::Module module = Python::import(parameterSet.get<std::string>("initialIterate"));
- auto pythonInitialIterate = Python::makeFunction<TargetSpace::CoordinateType(const FieldVector<double,dim>&)>(module.get("f"));
+ auto pythonInitialIterate = Python::makeFunction<TargetSpace::CoordinateType(const FieldVector<double,dimworld>&)>(module.get("f"));
std::vector<TargetSpace::CoordinateType> v;
+ using namespace Functions::BasisFactory;
+
+ auto powerBasis = makeBasis(
+ gridView,
+ power<TargetSpace::CoordinateType::dimension>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
#ifdef LAGRANGE
- Dune::Functions::interpolate(feBasis, v, pythonInitialIterate);
+ Dune::Functions::interpolate(powerBasis, v, pythonInitialIterate);
#else
- Dune::Functions::interpolate(feBasis, v, pythonInitialIterate, lower, upper, elements, order);
+ Dune::Functions::interpolate(powerBasis, v, pythonInitialIterate, lower, upper, elements, order);
#endif
for (size_t i=0; i<x.size(); i++)
@@ -225,9 +285,8 @@ int main (int argc, char *argv[])
// ////////////////////////////////////////////////////////////
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
- typedef LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
- //typedef GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace> LocalInterpolationRule;
- std::cout << "Using local interpolation: " << className<LocalInterpolationRule>() << std::endl;
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
+ using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<dim, double, FEBasis::LocalView::Tree::FiniteElement, ATargetSpace>;
// Assembler using ADOL-C
std::shared_ptr<LocalGeodesicFEStiffness<FEBasis,ATargetSpace> > localEnergy;
@@ -235,13 +294,26 @@ int main (int argc, char *argv[])
std::string energy = parameterSet.get<std::string>("energy");
if (energy == "harmonic")
{
-
- localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, LocalInterpolationRule, ATargetSpace>);
+ if (parameterSet["interpolationMethod"] == "geodesic")
+ localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, GeodesicInterpolationRule, ATargetSpace>);
+ else if (parameterSet["interpolationMethod"] == "projected")
+ localEnergy.reset(new HarmonicEnergyLocalStiffness<FEBasis, ProjectedInterpolationRule, ATargetSpace>);
+ else
+ DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
} else if (energy == "chiral_skyrmion")
{
-
- localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, LocalInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
+// // Doesn't work: we are not inside of a template
+// if constexpr (std::is_same<TargetSpace, UnitVector<double,3> >::value)
+// {
+ if (parameterSet["interpolationMethod"] == "geodesic")
+ localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, GeodesicInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
+ else if (parameterSet["interpolationMethod"] == "projected")
+ localEnergy.reset(new GFE::ChiralSkyrmionEnergy<FEBasis, ProjectedInterpolationRule, adouble>(parameterSet.sub("energyParameters")));
+ else
+ DUNE_THROW(Exception, "Unknown interpolation method " << parameterSet["interpolationMethod"] << " requested!");
+// } else
+// DUNE_THROW(Exception, "Build program with TargetSpace = UnitVector<3> for the ChiralSkyrmion energy!");
} else
DUNE_THROW(Exception, "Unknown energy type '" << energy << "'");
@@ -273,9 +345,6 @@ int main (int argc, char *argv[])
// Solve!
// /////////////////////////////////////////////////////
- std::cout << "Energy: " << assembler.computeEnergy(x) << std::endl;
- //exit(0);
-
solver.setInitialIterate(x);
solver.solve();
@@ -285,20 +354,16 @@ int main (int argc, char *argv[])
// Output result
// //////////////////////////////
+ SubsamplingVTKWriter<GridView> vtkWriter(gridView,Dune::refinementLevels(order-1));
+ std::string baseName = "harmonicmaps-result-" + std::to_string(order) + "-" + std::to_string(numLevels);
+ fillVTKWriter(vtkWriter, feBasis, x, resultPath + baseName);
+
+ // Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
EmbeddedVectorType xEmbedded(x.size());
for (size_t i=0; i<x.size(); i++)
- xEmbedded[i] = x[i].globalCoordinates();
-
- auto xFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<TargetSpace::CoordinateType>(feBasis,TypeTree::hybridTreePath(),xEmbedded);
+ xEmbedded[i] = x[i].globalCoordinates();
- std::string baseName = "harmonicmaps-result-" + std::to_string(order) + "-" + std::to_string(numLevels);
-
- SubsamplingVTKWriter<GridView> vtkWriter(gridView,order-1);
- vtkWriter.addVertexData(xFunction, VTK::FieldInfo("orientation", VTK::FieldInfo::Type::vector, xEmbedded[0].size()));
- vtkWriter.write(resultPath + baseName);
-
- // Write the corresponding coefficient vector: verbatim in binary, to be completely lossless
std::ofstream outFile(baseName + ".data", std::ios_base::binary);
MatrixVector::Generic::writeBinary(outFile, xEmbedded);
outFile.close();
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
deleted file mode 100644
index 5702b0326df654ce4859cf3271e5aac3d75aab39..0000000000000000000000000000000000000000
--- a/src/mixed-cosserat-continuum.cc
+++ /dev/null
@@ -1,409 +0,0 @@
-#include <config.h>
-
-#include <fenv.h>
-
-// Includes for the ADOL-C automatic differentiation library
-// Need to come before (almost) all others.
-#include <adolc/adouble.h>
-#include <adolc/drivers/drivers.h> // use of "Easy to Use" drivers
-#include <adolc/taping.h>
-
-#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
-
-#include <array>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/onedgrid.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-
-#include <dune/grid/io/file/gmshreader.hh>
-
-#include <dune/functions/common/tuplevector.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
-#include <dune/functions/functionspacebases/compositebasis.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
-#include <dune/fufem/dunepython.hh>
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/norms/energynorm.hh>
-
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
-#include <dune/gfe/cosseratenergystiffness.hh>
-#include <dune/gfe/cosseratvtkwriter.hh>
-#include <dune/gfe/mixedgfeassembler.hh>
-#include <dune/gfe/mixedriemanniantrsolver.hh>
-
-// grid dimension
-const int dim = 2;
-
-using namespace Dune;
-
-/** \brief A constant vector-valued function, for simple Neumann boundary values */
-struct NeumannFunction
- : public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
-{
- NeumannFunction(const FieldVector<double,3> values,
- double homotopyParameter)
- : values_(values),
- homotopyParameter_(homotopyParameter)
- {}
-
- void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
- out = 0;
- out.axpy(homotopyParameter_, values_);
- }
-
- FieldVector<double,3> values_;
- double homotopyParameter_;
-};
-
-
-int main (int argc, char *argv[]) try
-{
- // initialize MPI, finalize is done automatically on exit
- Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
-
- // Start Python interpreter
- Python::start();
- Python::Reference main = Python::import("__main__");
- Python::run("import math");
-
- //feenableexcept(FE_INVALID);
- Python::runStream()
- << std::endl << "import sys"
- << std::endl << "sys.path.append('/home/sander/dune/dune-gfe/problems/')"
- << std::endl;
-
- using namespace Dune::TypeTree::Indices;
- typedef Dune::Functions::TupleVector<std::vector<RealTuple<double,3> >,
- std::vector<Rotation<double,3> > > SolutionType;
-
- // parse data file
- ParameterTree parameterSet;
- if (argc < 2)
- DUNE_THROW(Exception, "Usage: ./mixed-cosserat-continuum <parameter file>");
-
- ParameterTreeParser::readINITree(argv[1], parameterSet);
-
- ParameterTreeParser::readOptions(argc, argv, parameterSet);
-
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
- int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
- const double tolerance = parameterSet.get<double>("tolerance");
- const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
- const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
- const int multigridIterations = parameterSet.get<int>("numIt");
- const int nu1 = parameterSet.get<int>("nu1");
- const int nu2 = parameterSet.get<int>("nu2");
- const int mu = parameterSet.get<int>("mu");
- const int baseIterations = parameterSet.get<int>("baseIt");
- const double mgTolerance = parameterSet.get<double>("mgTolerance");
- const double baseTolerance = parameterSet.get<double>("baseTolerance");
- const bool instrumented = parameterSet.get<bool>("instrumented");
- std::string resultPath = parameterSet.get("resultPath", "");
-
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
- typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
-
- shared_ptr<GridType> grid;
-
- FieldVector<double,dim> lower(0), upper(1);
-
- if (parameterSet.get<bool>("structuredGrid")) {
-
- lower = parameterSet.get<FieldVector<double,dim> >("lower");
- upper = parameterSet.get<FieldVector<double,dim> >("upper");
-
- auto elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
- grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
-
- } else {
- std::string path = parameterSet.get<std::string>("path");
- std::string gridFile = parameterSet.get<std::string>("gridFile");
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
- }
-
- grid->globalRefine(numLevels-1);
-
- grid->loadBalance();
-
- if (mpiHelper.rank()==0)
- std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
-
- typedef GridType::LeafGridView GridView;
- GridView gridView = grid->leafGridView();
-
- using namespace Dune::Functions::BasisBuilder;
-
- auto compositeBasis = makeBasis(
- gridView,
- composite(
- pq<2>(),
- pq<1>()
- )
- );
-
- typedef Dune::Functions::PQkNodalBasis<GridView,2> DeformationFEBasis;
- typedef Dune::Functions::PQkNodalBasis<GridView,1> OrientationFEBasis;
-
- DeformationFEBasis deformationFEBasis(gridView);
- OrientationFEBasis orientationFEBasis(gridView);
-
- // Construct fufem-style function space bases to ease the transition to dune-functions
- typedef DuneFunctionsBasis<DeformationFEBasis> FufemDeformationFEBasis;
- FufemDeformationFEBasis fufemDeformationFEBasis(deformationFEBasis);
-
- typedef DuneFunctionsBasis<OrientationFEBasis> FufemOrientationFEBasis;
- FufemOrientationFEBasis fufemOrientationFEBasis(orientationFEBasis);
-
-
-
- std::cout << "Deformation: " << deformationFEBasis.size() << ", orientation: " << orientationFEBasis.size() << std::endl;
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
-
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
-
- GridType::Codim<dim>::LeafIterator vIt = gridView.begin<dim>();
- GridType::Codim<dim>::LeafIterator vEndIt = gridView.end<dim>();
-
- const GridView::IndexSet& indexSet = gridView.indexSet();
-
- // Make Python function that computes which vertices are on the Dirichlet boundary,
- // based on the vertex positions.
- std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
-
- // Same for the Neumann boundary
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
-
- for (; vIt!=vEndIt; ++vIt) {
-
- bool isDirichlet;
- pythonDirichletVertices.evaluate(vIt->geometry().corner(0), isDirichlet);
- dirichletVertices[indexSet.index(*vIt)] = isDirichlet;
-
- bool isNeumann;
- pythonNeumannVertices.evaluate(vIt->geometry().corner(0), isNeumann);
- neumannVertices[indexSet.index(*vIt)] = isNeumann;
-
- }
-
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
-
- BitSetVector<1> neumannNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(neumannBoundary,fufemDeformationFEBasis,neumannNodes);
-
-
- if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
-
-
- BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemDeformationFEBasis,deformationDirichletNodes);
-
- BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
- for (size_t i=0; i<deformationFEBasis.size(); i++)
- if (deformationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- deformationDirichletDofs[i][j] = true;
-
- BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,fufemOrientationFEBasis,orientationDirichletNodes);
-
- BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
- for (size_t i=0; i<orientationFEBasis.size(); i++)
- if (orientationDirichletNodes[i][0])
- for (int j=0; j<3; j++)
- orientationDirichletDofs[i][j] = true;
-
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
-
- SolutionType x;
-
- x[_0].resize(deformationFEBasis.size());
-
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
- std::vector<FieldVector<double,3> > v;
- ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
-
- for (size_t i=0; i<x[_0].size(); i++)
- x[_0][i] = v[i];
-
- x[_1].resize(orientationFEBasis.size());
-#if 0
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
-
- std::vector<FieldVector<double,3> > v;
- Functions::interpolate(feBasis, v, pythonInitialDeformation);
-
- for (size_t i=0; i<x.size(); i++)
- xDisp[i] = v[i];
-#endif
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
-
- // Output initial iterate (of homotopy loop)
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
- orientationFEBasis,x[_1],
- resultPath + "mixed-cosserat_homotopy_0");
-
- for (int i=0; i<numHomotopySteps; i++) {
-
- double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
- if (mpiHelper.rank()==0)
- std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
-
-
- // ////////////////////////////////////////////////////////////
- // Create an assembler for the energy functional
- // ////////////////////////////////////////////////////////////
-
- const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
-
- shared_ptr<NeumannFunction> neumannFunction;
- if (parameterSet.hasKey("neumannValues"))
- neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
- homotopyParameter);
-
-
- if (mpiHelper.rank() == 0) {
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
- }
-
- // Assembler using ADOL-C
- CosseratEnergyLocalStiffness<decltype(compositeBasis),
- 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
- &neumannBoundary,
- neumannFunction,
- nullptr);
-
- MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
-
- MixedGFEAssembler<decltype(compositeBasis),
- RealTuple<double,3>,
- Rotation<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
-
- // /////////////////////////////////////////////////
- // Create a Riemannian trust-region solver
- // /////////////////////////////////////////////////
-
- MixedRiemannianTrustRegionSolver<GridType,
- decltype(compositeBasis),
- DeformationFEBasis, RealTuple<double,3>,
- OrientationFEBasis, Rotation<double,3> > solver;
- solver.setup(*grid,
- &assembler,
- deformationFEBasis,
- orientationFEBasis,
- x,
- deformationDirichletDofs,
- orientationDirichletDofs,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- instrumented);
-
- solver.setScaling(parameterSet.get<FieldVector<double,6> >("trustRegionScaling"));
-
- ////////////////////////////////////////////////////////
- // Set Dirichlet values
- ////////////////////////////////////////////////////////
-
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
-
- Python::Callable C = dirichletValuesClass.get("DirichletValues");
-
- // Call a constructor.
- Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
-
- // Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > deformationDirichletValues(dirichletValuesPythonObject.get("deformation"));
- PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(dirichletValuesPythonObject.get("orientation"));
-
- std::vector<FieldVector<double,3> > ddV;
- ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
-
- std::vector<FieldMatrix<double,3,3> > dOV;
- ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
-
- for (size_t j=0; j<x[_0].size(); j++)
- if (deformationDirichletNodes[j][0])
- x[_0][j] = ddV[j];
-
- for (size_t j=0; j<x[_1].size(); j++)
- if (orientationDirichletNodes[j][0])
- x[_1][j].set(dOV[j]);
-
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- solver.setInitialIterate(x);
- solver.solve();
-
- x = solver.getSol();
-
- // Output result of each homotopy step
- std::stringstream iAsAscii;
- iAsAscii << i+1;
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,x[_0],
- orientationFEBasis,x[_1],
- resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
-
- }
-
- // //////////////////////////////
- // Output result
- // //////////////////////////////
-
- // finally: compute the average deformation of the Neumann boundary
- // That is what we need for the locking tests
- FieldVector<double,3> averageDef(0);
- for (size_t i=0; i<x[_0].size(); i++)
- if (neumannNodes[i][0])
- averageDef += x[_0][i].globalCoordinates();
- averageDef /= neumannNodes.count();
-
- if (mpiHelper.rank()==0)
- {
- std::cout << "Neumann values = " << parameterSet.get<FieldVector<double, 3> >("neumannValues") << " "
- << ", average deflection: " << averageDef << std::endl;
- }
-
- } catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 4035fd75736fbd73201cdafbfd0b88c74c463443..b1a268f13fe44ede63c4e355b5e1640b3e61dfce 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -18,12 +18,7 @@ set(TESTS
targetspacetest
)
-add_directory_test_target(_test_target)
-add_dependencies(${_test_target} ${TESTS})
-
foreach(_test ${TESTS})
- add_executable(${_test} ${_test}.cc)
- target_link_libraries(${_test} ${DUNE_LIBS})
- add_dune_ug_flags(${_test})
- add_test(${_test} ${_test})
+ dune_add_test(SOURCES ${_test}.cc)
+ target_compile_options(${_test} PRIVATE -g)
endforeach()
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 603db78ffafb3bd102314070f213876e44419ac7..523d7406fba0780826d01b43adff77092d9216a0 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -25,6 +25,17 @@ typedef double FDType;
#include <adolc/taping.h>
#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/typetraits.hh>
+// This specialization is needed (in particular) to make copying of FieldMatrices work.
+namespace Dune {
+ template <>
+ struct IsNumber<adouble>
+ {
+ constexpr static bool value = true;
+ };
+}
+
#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
@@ -33,7 +44,7 @@ typedef double FDType;
#include <dune/istl/io.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
@@ -429,7 +440,7 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
- typedef Functions::PQkNodalBasis<GridView,2> FEBasis;
+ typedef Functions::LagrangeBasis<GridView,1> FEBasis;
FEBasis feBasis(gridView);
// /////////////////////////////////////////
@@ -577,7 +588,7 @@ int main (int argc, char *argv[]) try
}
// //////////////////////////////
- } catch (Exception e) {
+ } catch (Exception& e) {
std::cout << e << std::endl;
diff --git a/test/averagedistanceassemblertest.cc b/test/averagedistanceassemblertest.cc
index 0f49c4a40b4e511cd225eeec82fe5490f692a5ae..5aa37cc1d06755541bbf6d4f8f3ce2db1901c458 100644
--- a/test/averagedistanceassemblertest.cc
+++ b/test/averagedistanceassemblertest.cc
@@ -66,7 +66,7 @@ void testWeightSet(const std::vector<TargetSpace>& corners,
int quadOrder = 3;
const Dune::QuadratureRule<double, dim>& quad
- = Dune::QuadratureRules<double, dim>::rule(GeometryType(GeometryType::simplex,dim), quadOrder);
+ = Dune::QuadratureRules<double, dim>::rule(GeometryTypes::simplex(dim), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 6881cb5d5f0649c95a7961fe3f6f158f2b03dd4c..d74cc9a4f372fc22a5f9f2494a1ba08ed443b4f3 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -6,7 +6,7 @@
#include <dune/geometry/type.hh>
#include <dune/geometry/quadraturerules.hh>
-#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/fufem/functions/constantfunction.hh>
@@ -119,7 +119,7 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
Tensor3<double,3,3,domainDim> DR;
- typedef Dune::Functions::PQkNodalBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
// evaluate fd approximation of derivative
@@ -157,7 +157,7 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1";
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> FEBasis;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> FEBasis;
FEBasis feBasis(grid->leafGridView());
CosseratEnergyLocalStiffness<FEBasis,3> assembler(materialParameters,
@@ -349,7 +349,7 @@ int main(int argc, char** argv)
// Create a local assembler object
////////////////////////////////////////////////////////////////////////////
- typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> Basis;
+ typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> Basis;
Basis basis(grid->leafGridView());
std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
diff --git a/test/frameinvariancetest.cc b/test/frameinvariancetest.cc
index 95180a69423e3c017acd89001dd6c8c8500f0c7f..24b4151c157ea72326730a4bc91ccc78a8d3adf8 100644
--- a/test/frameinvariancetest.cc
+++ b/test/frameinvariancetest.cc
@@ -23,8 +23,6 @@ using std::string;
int main (int argc, char *argv[]) try
{
// Some types that I need
- typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > MatrixType;
- typedef BlockVector<FieldVector<double, blocksize> > CorrectionType;
typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
// Problem settings
@@ -35,8 +33,10 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////
typedef OneDGrid GridType;
GridType grid(numRodBaseElements, 0, 1);
+ using GridView = GridType::LeafGridView;
+ GridView gridView = grid.leafGridView();
- SolutionType x(grid.size(1));
+ SolutionType x(gridView.size(1));
// //////////////////////////
// Initial solution
@@ -57,10 +57,7 @@ int main (int argc, char *argv[]) try
// Create a second, rotated copy of the configuration
// /////////////////////////////////////////////////////////////////////
- FieldVector<double,3> displacement;
- displacement[0] = 0;
- displacement[1] = 0;
- displacement[2] = 0;
+ FieldVector<double,3> displacement {0, 0, 0};
FieldVector<double,3> axis(0); axis[0]=1;
Rotation<double,3> rotation(axis,M_PI/2);
@@ -88,20 +85,20 @@ int main (int argc, char *argv[]) try
writeRod(x,"rod");
writeRod(rotatedX, "rotated");
- RodLocalStiffness<GridType::LeafGridView,double> assembler(grid.leafGridView(),
+ RodLocalStiffness<GridView,double> assembler(gridView,
1,1,1,1e6,0.3);
for (int i=1; i<2; i++) {
double p = double(i)/2;
- assembler.getStrain(x,*grid.lbegin<0>(0), p);
- assembler.getStrain(rotatedX,*grid.lbegin<0>(0), p);
+ assembler.getStrain(x,*gridView.begin<0>(), p);
+ assembler.getStrain(rotatedX,*gridView.begin<0>(), p);
}
- } catch (Exception e) {
+ } catch (Exception& e) {
- std::cout << e << std::endl;
+ std::cout << e.what() << std::endl;
}
diff --git a/test/globalgfetestfunctionbasistest.cc b/test/globalgfetestfunctionbasistest.cc
index beb1d7197b7c233dd603ebd2ff31243f93f39015..16e8cc5399c75cbbdcc4601652760d5a8a79cc0b 100644
--- a/test/globalgfetestfunctionbasistest.cc
+++ b/test/globalgfetestfunctionbasistest.cc
@@ -41,19 +41,17 @@ void test()
typedef GlobalGFETestFunctionBasis<P1Basis,TargetSpace> GlobalBasis;
GlobalBasis basis(p1Basis,testPoints);
- typedef typename OneDGrid::Codim<0>::LeafIterator ElementIterator;
- ElementIterator eIt = grid.leafbegin<0>();
- ElementIterator eEndIt = grid.leafend<0>();
-
- for (; eIt != eEndIt; ++eIt) {
-
- const typename GlobalBasis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
+ for (const auto element : elements(grid.leafGridView()))
+ {
+ const typename GlobalBasis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(element);
FieldVector<double,1> stupidTestPoint(0);
std::vector<std::array<typename TargetSpace::EmbeddedTangentVector, TargetSpace::TangentVector::dimension> > values;
lfe.localBasis().evaluateFunction(stupidTestPoint, values);
for(size_t i=0;i<values.size();i++) {
- std::cout<<values[i]<<std::endl;
+ for (auto v : values[i])
+ std::cout << v << " ";
+ std::cout << std::endl;
std::cout<<lfe.localCoefficients().localKey(i)<<std::endl;
}
//int i = basis.index(*eIt,1);
diff --git a/test/harmonicenergytest.cc b/test/harmonicenergytest.cc
index 560d6552f121aeadbe64c2d88cfcbc907e8b241f..b0b9466927bf87194a8c862c19577e0f41d7d0e2 100644
--- a/test/harmonicenergytest.cc
+++ b/test/harmonicenergytest.cc
@@ -2,10 +2,11 @@
#include <dune/grid/uggrid.hh>
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/harmonicenergystiffness.hh>
+#include <dune/gfe/localgeodesicfefunction.hh>
#include "multiindex.hh"
#include "valuefactory.hh"
@@ -16,18 +17,20 @@ typedef UnitVector<double,3> TargetSpace;
using namespace Dune;
-template <class GridType>
-void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients) {
+template <class Basis>
+void testEnergy(const Basis& basis, const std::vector<TargetSpace>& coefficients)
+{
+ using GridView = typename Basis::GridView;
+ GridView gridView = basis.gridView();
- PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
- typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
-
- //LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
+ using GeodesicInterpolationRule = LocalGeodesicFEFunction<dim, double, typename Basis::LocalView::Tree::FiniteElement, TargetSpace>;
-
- HarmonicEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,TargetSpace> assembler;
+ HarmonicEnergyLocalStiffness<Basis,GeodesicInterpolationRule,TargetSpace> assembler;
std::vector<TargetSpace> rotatedCoefficients(coefficients.size());
+ auto localView = basis.localView();
+ localView.bind(*gridView.template begin<0>());
+
for (int i=0; i<10; i++) {
Rotation<double,3> rotation(FieldVector<double,3>(1), double(i));
@@ -40,68 +43,13 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
rotatedCoefficients[j] = tmp;
}
- std::cout << "energy: " << assembler.energy(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
+ std::cout << "energy: " << assembler.energy(localView,
rotatedCoefficients) << std::endl;
- std::vector<typename TargetSpace::EmbeddedTangentVector> rotatedGradient;
- assembler.assembleEmbeddedGradient(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
- rotatedCoefficients,
- rotatedGradient);
-
- for (size_t j=0; j<coefficients.size(); j++) {
- FieldVector<double,3> tmp;
- matrix.mtv(rotatedGradient[j], tmp);
- std::cout << "gradient: " << tmp << std::endl;
- }
-
}
}
-template <class GridType>
-void testGradientOfEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients)
-{
- PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
- typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
-
- //LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
- HarmonicEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,TargetSpace> assembler;
-
- std::vector<typename TargetSpace::EmbeddedTangentVector> gradient;
- assembler.assembleEmbeddedGradient(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
- coefficients,
- gradient);
-
- ////////////////////////////////////////////////////////////////////////////////////////
- // Assemble finite difference approximation of the gradient
- ////////////////////////////////////////////////////////////////////////////////////////
-
- std::vector<typename TargetSpace::EmbeddedTangentVector> fdGradient;
- assembler.assembleEmbeddedFDGradient(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
- coefficients,
- fdGradient);
-
- double diff = 0;
- for (size_t i=0; i<fdGradient.size(); i++)
- diff = std::max(diff, (gradient[i] - fdGradient[i]).infinity_norm());
-
- if (diff > 1e-6) {
-
- std::cout << "gradient: " << std::endl;
- for (size_t i=0; i<gradient.size(); i++)
- std::cout << gradient[i] << std::endl;
-
- std::cout << "fd gradient: " << std::endl;
- for (size_t i=0; i<fdGradient.size(); i++)
- std::cout << fdGradient[i] << std::endl;
-
- }
-
-}
template <int domainDim>
void testUnitVector3d()
@@ -126,10 +74,15 @@ void testUnitVector3d()
std::vector<unsigned int> v(domainDim+1);
for (int i=0; i<domainDim+1; i++)
v[i] = i;
- factory.insertElement(GeometryType(GeometryType::simplex,dim), v);
+ factory.insertElement(GeometryTypes::simplex(dim), v);
+
+ auto grid = std::unique_ptr<const GridType>(factory.createGrid());
+ using GridView = typename GridType::LeafGridView;
+ auto gridView = grid->leafGridView();
- const GridType* grid = factory.createGrid();
-
+ // A function space basis
+ using Basis = Functions::LagrangeBasis<GridView,1>;
+ Basis basis(gridView);
// //////////////////////////////////////////////////////////
// Test whether the energy is invariant under isometries
@@ -149,8 +102,7 @@ void testUnitVector3d()
for (int j=0; j<dim+1; j++)
coefficients[j] = testPoints[index[j]];
- testEnergy<GridType>(grid, coefficients);
- testGradientOfEnergy(grid, coefficients);
+ testEnergy<Basis>(basis, coefficients);
}
diff --git a/test/interillustration.cc b/test/interillustration.cc
index 89a0164526dc87b4d337e6c94c422b06cc767034..fdddd9e052c01e3fe861d8e25014dbbfbc19333c 100644
--- a/test/interillustration.cc
+++ b/test/interillustration.cc
@@ -7,7 +7,7 @@
#include <dune/localfunctions/lagrange/p1.hh>
#include <dune/localfunctions/lagrange/p2.hh>
-#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discretescalarglobalbasisfunction.hh>
#include <dune/gfe/unitvector.hh>
@@ -143,7 +143,7 @@ void interpolate(const LocalFEFunctionType& localGeodesicFEFunction,
// stupid, can't instantiate deformedGrid with a const grid
DeformedGridType deformedGrid(const_cast<GridType&>(*grid), deformationFunction);
- typedef Functions::PQ1NodalBasis<typename DeformedGridType::LeafGridView > FEBasis;
+ typedef Functions::LagrangeBasis<typename DeformedGridType::LeafGridView, 1> FEBasis;
FEBasis feBasis(deformedGrid.leafGridView());
Functions::DiscreteScalarGlobalBasisFunction<decltype(feBasis),decltype(variation0)> variation0Function(feBasis,variation0);
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index aea15a5580e3e36ef2a85cf7c732dc9a6fd4671b..37f874f4a1e3aef70e84610db620b1cd68cfc7b8 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -342,39 +342,34 @@ int main()
std::cout << std::setw(15) << std::setprecision(12);
- GeometryType element;
-
////////////////////////////////////////////////////////////////
// Test functions on 1d elements
////////////////////////////////////////////////////////////////
- element.makeSimplex(1);
- test<RealTuple<double,1>,1>(element);
- test<UnitVector<double,2>,1>(element);
- test<UnitVector<double,3>,1>(element);
- test<Rotation<double,3>,1>(element);
- test<RigidBodyMotion<double,3>,1>(element);
+ test<RealTuple<double,1>,1>(GeometryTypes::simplex(1));
+ test<UnitVector<double,2>,1>(GeometryTypes::simplex(1));
+ test<UnitVector<double,3>,1>(GeometryTypes::simplex(1));
+ test<Rotation<double,3>,1>(GeometryTypes::simplex(1));
+ test<RigidBodyMotion<double,3>,1>(GeometryTypes::simplex(1));
////////////////////////////////////////////////////////////////
// Test functions on 2d simplex elements
////////////////////////////////////////////////////////////////
- element.makeSimplex(2);
- test<RealTuple<double,1>,2>(element);
- test<UnitVector<double,2>,2>(element);
- test<UnitVector<double,3>,2>(element);
- test<Rotation<double,3>,2>(element);
- test<RigidBodyMotion<double,3>,2>(element);
+ test<RealTuple<double,1>,2>(GeometryTypes::simplex(2));
+ test<UnitVector<double,2>,2>(GeometryTypes::simplex(2));
+ test<UnitVector<double,3>,2>(GeometryTypes::simplex(2));
+ test<Rotation<double,3>,2>(GeometryTypes::simplex(2));
+ test<RigidBodyMotion<double,3>,2>(GeometryTypes::simplex(2));
////////////////////////////////////////////////////////////////
// Test functions on 2d quadrilateral elements
////////////////////////////////////////////////////////////////
- element.makeCube(2);
- test<RealTuple<double,1>,2>(element);
- test<UnitVector<double,2>,2>(element);
- test<UnitVector<double,3>,2>(element);
- test<Rotation<double,3>,2>(element);
- test<RigidBodyMotion<double,3>,2>(element);
+ test<RealTuple<double,1>,2>(GeometryTypes::cube(2));
+ test<UnitVector<double,2>,2>(GeometryTypes::cube(2));
+ test<UnitVector<double,3>,2>(GeometryTypes::cube(2));
+ test<Rotation<double,3>,2>(GeometryTypes::cube(2));
+ test<RigidBodyMotion<double,3>,2>(GeometryTypes::cube(2));
}
diff --git a/test/localgfetestfunctiontest.cc b/test/localgfetestfunctiontest.cc
index dee684c72e298c63352e8df82b833b65b00a4489..4b9769c95f7348e66363ed16f9eaa88d8ec5ae38 100644
--- a/test/localgfetestfunctiontest.cc
+++ b/test/localgfetestfunctiontest.cc
@@ -40,15 +40,12 @@ void test()
PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
- GeometryType simplex;
- simplex.makeSimplex(domainDim);
-
for (int i=0; i<numIndices; i++, ++index) {
for (int j=0; j<domainDim+1; j++)
coefficients[j] = testPoints[index[j]];
- LocalGfeTestFunctionFiniteElement<LocalFiniteElement,TargetSpace> testFunctionSet(feCache.get(simplex),coefficients);
+ LocalGfeTestFunctionFiniteElement<LocalFiniteElement,TargetSpace> testFunctionSet(feCache.get(GeometryTypes::simplex(domainDim)),coefficients);
FieldVector<double,domainDim> stupidTestPoint(0);
@@ -83,7 +80,7 @@ int main() try
test<Rotation<double,3>, 1>();
test<Rotation<double,3>, 2>();
-} catch (Exception e) {
+} catch (Exception& e) {
std::cout << e << std::endl;
}
diff --git a/test/rodassemblertest.cc b/test/rodassemblertest.cc
index 6542a07417e99dd8c3611b6b5516f9c954f707ac..61cb237383027791a494967816d9d4a23ba556be 100644
--- a/test/rodassemblertest.cc
+++ b/test/rodassemblertest.cc
@@ -1,5 +1,7 @@
#include <config.h>
+#include <array>
+
#include <dune/grid/onedgrid.hh>
#include <dune/istl/io.hh>
@@ -28,7 +30,7 @@ void infinitesimalVariation(RigidBodyMotion<double,3>& c, double eps, int i)
template <class GridType>
void strainFD(const std::vector<RigidBodyMotion<double,3> >& x,
double pos,
- Dune::array<Dune::FieldMatrix<double,2,6>, 6>& fdStrainDerivatives,
+ std::array<FieldMatrix<double,2,6>, 6>& fdStrainDerivatives,
const RodAssembler<GridType,3>& assembler)
{
assert(x.size()==2);
@@ -73,8 +75,8 @@ void strainFD(const std::vector<RigidBodyMotion<double,3> >& x,
template <class GridType>
void strain2ndOrderFD(const std::vector<RigidBodyMotion<double,3> >& x,
double pos,
- Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer,
- Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer,
+ std::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer,
+ std::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer,
const RodAssembler<GridType,3>& assembler)
{
assert(x.size()==2);
@@ -184,8 +186,8 @@ void strain2ndOrderFD2(const std::vector<RigidBodyMotion<double,3> >& x,
double pos,
Dune::FieldVector<double,1> shapeGrad[2],
Dune::FieldVector<double,1> shapeFunction[2],
- Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer,
- Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer,
+ std::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer,
+ std::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer,
const RodAssembler<GridType,3>& assembler)
{
assert(x.size()==2);
@@ -214,10 +216,10 @@ void strain2ndOrderFD2(const std::vector<RigidBodyMotion<double,3> >& x,
infinitesimalVariation(forwardSolution[k], eps, l+3);
infinitesimalVariation(backwardSolution[k], -eps, l+3);
- Dune::array<Dune::FieldMatrix<double,2,6>, 6> forwardStrainDer;
+ std::array<FieldMatrix<double,2,6>, 6> forwardStrainDer;
assembler.strainDerivative(forwardSolution, pos, shapeGrad, shapeFunction, forwardStrainDer);
- Dune::array<Dune::FieldMatrix<double,2,6>, 6> backwardStrainDer;
+ std::array<FieldMatrix<double,2,6>, 6> backwardStrainDer;
assembler.strainDerivative(backwardSolution, pos, shapeGrad, shapeFunction, backwardStrainDer);
for (int i=0; i<2; i++) {
@@ -571,7 +573,7 @@ int main (int argc, char *argv[]) try
hessianFDCheck(x, hessianMatrix, rodAssembler);
// //////////////////////////////
- } catch (Exception e) {
+ } catch (Exception& e) {
std::cout << e << std::endl;
diff --git a/test/rotationtest.cc b/test/rotationtest.cc
index bf349d2f9c21bc4f762d92fc62df652cb28c2493..ac4f2f6a32aeceb42bd568c64157ad318e595458 100644
--- a/test/rotationtest.cc
+++ b/test/rotationtest.cc
@@ -395,7 +395,7 @@ int main (int argc, char *argv[]) try
return not passed;
- } catch (Exception e) {
+ } catch (Exception& e) {
std::cout << e << std::endl;
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index 820f3e81c8199abf86be9a924fe6ecbac1bb5af0..cadfe4818e3f3c0c8bfce2ff69c7b0e22acaa9cb 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -384,7 +384,7 @@ int main() try
//
// test<HyperbolicHalfspacePoint<double,2> >();
-} catch (Exception e) {
+} catch (Exception& e) {
std::cout << e << std::endl;