From 00b4f2b8368d75238fff3c316104006cb65fbdd3 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 22 May 2018 15:49:09 +0200
Subject: [PATCH] Update harmonicenergytest.cc
This makes the test build again.
It also removes the FD test for the gradient of the harmonic energy.
Since the HarmonicEnergyStiffness class does not implement the
gradient anymore, there is nothing to test here.
---
test/harmonicenergytest.cc | 92 +++++++++-----------------------------
1 file changed, 22 insertions(+), 70 deletions(-)
diff --git a/test/harmonicenergytest.cc b/test/harmonicenergytest.cc
index 560d6552..b0b94669 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);
}
--
GitLab