diff --git a/Makefile.am b/Makefile.am
index c066ea53e9cb146267afb75e6b84320e8c499e5e..b023fe84d786c666492be9ca224632e4f25b6c5c 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -15,7 +15,7 @@ ADOLC_LIBS = -ladolc
noinst_PROGRAMS = finite-strain-elasticity \
mixed-cosserat-continuum \
- rodobstacle harmonicmaps-eoc rod-eoc
+ rodobstacle rod-eoc
finite_strain_elasticity_SOURCES = finite-strain-elasticity.cc
finite_strain_elasticity_CXXFLAGS = $(UG_CPPFLAGS) $(IPOPT_CPPFLAGS) \
@@ -37,11 +37,6 @@ rodobstacle_SOURCES = rodobstacle.cc
rod_eoc_SOURCES = rod-eoc.cc
-harmonicmaps_eoc_SOURCES = harmonicmaps-eoc.cc
-harmonicmaps_eoc_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
-harmonicmaps_eoc_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
- $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
-
# don't follow the full GNU-standard
# we need automake 1.5
AUTOMAKE_OPTIONS = foreign 1.5
diff --git a/harmonicmaps-eoc.cc b/harmonicmaps-eoc.cc
deleted file mode 100644
index 9f3f3d448ac76b7692d8163f538ce421a3e3c62e..0000000000000000000000000000000000000000
--- a/harmonicmaps-eoc.cc
+++ /dev/null
@@ -1,382 +0,0 @@
-#include <config.h>
-
-//#define HARMONIC_ENERGY_FD_GRADIENT
-//#define HARMONIC_ENERGY_FD_INNER_GRADIENT
-//#define SECOND_ORDER
-//#define THIRD_ORDER
-const int order = 1;
-
-#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/vtk/vtkwriter.hh>
-#include <dune/grid/io/file/amirameshreader.hh>
-#include <dune/grid/io/file/gmshreader.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
-#include <dune/fufem/functiontools/boundarydofs.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functions/vtkbasisgridfunction.hh>
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/norms/h1seminorm.hh>
-
-#include <dune/gfe/unitvector.hh>
-#include <dune/gfe/harmonicenergystiffness.hh>
-#include <dune/gfe/geodesicfeassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/geodesicfefunctionadaptor.hh>
-#include <dune/gfe/gfedifferencenormsquared.hh>
-
-// grid dimension
-const int dim = 2;
-
-typedef UnitVector<double,3> TargetSpace;
-typedef std::vector<TargetSpace> SolutionType;
-
-const int blocksize = TargetSpace::TangentVector::dimension;
-
-using namespace Dune;
-using std::string;
-
-struct DirichletFunction
- : public Dune::VirtualFunction<FieldVector<double,dim>, TargetSpace::CoordinateType >
-{
- void evaluate(const FieldVector<double, dim>& x, TargetSpace::CoordinateType& out) const {
-
-#if 0
- FieldVector<double,3> axis;
- axis[0] = x[0]; axis[1] = x[1]; axis[2] = 1;
- Rotation<double,3> rotation(axis, x.two_norm()*M_PI*3);
-
- FieldMatrix<double,3,3> rMat;
- rotation.matrix(rMat);
- out = rMat[2];
-#endif
-#if 1 // This data seems to have the necessary smoothness to have optimal
- // convergence order even for 3rd-order elements
- double localX = (x[0] + 2)/4;
- double localY = (x[1] + 1)/2;
- double angleX = 0.5 * M_PI * sin(M_PI*localX);
- double angleY = 0.5 * M_PI * sin(M_PI*localY);
-
- // Rotation matrix around the x axis
- FieldMatrix<double,TargetSpace::CoordinateType::dimension,TargetSpace::CoordinateType::dimension> rotationX(0);
- rotationX[0][0] = 1;
- rotationX[1][1] = cos(angleY);
- rotationX[1][2] = -sin(angleY);
- rotationX[2][1] = sin(angleY);
- rotationX[2][2] = cos(angleY);
-
- // Rotation matrix around the y axis
- FieldMatrix<double,TargetSpace::CoordinateType::dimension,TargetSpace::CoordinateType::dimension> rotationY(0);
- rotationY[0][0] = cos(angleX);
- rotationY[0][2] = -sin(angleX);
- rotationY[1][1] = 1;
- rotationY[2][0] = sin(angleX);
- rotationY[2][2] = cos(angleX);
-
- //angle *= -4*x[1]*(x[1]-1);
- TargetSpace::CoordinateType unitVector(0); unitVector[2] = 1;
-
- TargetSpace::CoordinateType tmp;
- rotationX.mv(unitVector,tmp);
- rotationY.mv(tmp,out);
-#endif
- }
-};
-
-
-template <class GridType>
-void solve (const shared_ptr<GridType>& grid,
- SolutionType& x,
- int numLevels,
- const ParameterTree& parameters)
-{
- // read solver setting
- const double innerTolerance = parameters.get<double>("innerTolerance");
- const double tolerance = parameters.get<double>("tolerance");
- const int maxTrustRegionSteps = parameters.get<int>("maxTrustRegionSteps");
- const double initialTrustRegionRadius = parameters.get<double>("initialTrustRegionRadius");
- const int multigridIterations = parameters.get<int>("numIt");
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
-
- BitSetVector<1> allNodes(grid->size(dim));
- allNodes.setAll();
- BoundaryPatch<typename GridType::LeafGridView> dirichletBoundary(grid->leafGridView(), allNodes);
-
-#if defined THIRD_ORDER
- typedef P3NodalBasis<typename GridType::LeafGridView,double> FEBasis;
-#elif defined SECOND_ORDER
- typedef P2NodalBasis<typename GridType::LeafGridView,double> FEBasis;
-#else
- typedef P1NodalBasis<typename GridType::LeafGridView,double> FEBasis;
-#endif
- FEBasis feBasis(grid->leafGridView());
-
- BitSetVector<blocksize> dirichletNodes;
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
-
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
-
- x.resize(feBasis.size());
-
- BlockVector<TargetSpace::CoordinateType> dirichletFunctionValues;
- DirichletFunction dirichletFunction;
- Functions::interpolate(feBasis, dirichletFunctionValues, dirichletFunction);
-
- TargetSpace::CoordinateType innerValue(0);
- innerValue[0] = 1;
- innerValue[1] = 0;
-
- for (size_t i=0; i<x.size(); i++)
- x[i] = (dirichletNodes[i][0]) ? dirichletFunctionValues[i] : innerValue;
-
- // ////////////////////////////////////////////////////////////
- // Create an assembler for the Harmonic Energy Functional
- // ////////////////////////////////////////////////////////////
-
- HarmonicEnergyLocalStiffness<typename GridType::LeafGridView,typename FEBasis::LocalFiniteElement, TargetSpace> harmonicEnergyLocalStiffness;
-
- GeodesicFEAssembler<FEBasis,TargetSpace> assembler(grid->leafGridView(),
- &harmonicEnergyLocalStiffness);
-
- // ///////////////////////////////////////////
- // Create a solver for the rod problem
- // ///////////////////////////////////////////
-
- RiemannianTrustRegionSolver<GridType,TargetSpace> solver;
-
- solver.setup(*grid,
- &assembler,
- x,
- dirichletNodes,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- innerTolerance,
- 1, 3, 3,
- 100, // iterations of the base solver
- 1e-8, // base tolerance
- false); // instrumentation
-
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- solver.setInitialIterate(x);
- solver.solve();
-
- x = solver.getSol();
-}
-
-int main (int argc, char *argv[]) try
-{
- // parse data file
- ParameterTree parameterSet;
- if (argc==2)
- ParameterTreeParser::readINITree(argv[1], parameterSet);
- else
- ParameterTreeParser::readINITree("harmonicmaps-eoc.parset", parameterSet);
-
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
-
- // grid file
- std::string gridFileName = parameterSet.get<std::string>("gridFile");
-
- // only if a structured grid is used
- const int numBaseElements = parameterSet.get<int>("numBaseElements");
- FieldVector<double,dim> lowerLeft = parameterSet.get<FieldVector<double,dim> >("lowerLeft");
- FieldVector<double,dim> upperRight = parameterSet.get<FieldVector<double,dim> >("upperRight");
-
- // ///////////////////////////////////////////////////////////
- // First compute the 'exact' solution on a very fine grid
- // ///////////////////////////////////////////////////////////
-
- typedef std::conditional<dim==1,OneDGrid,UGGrid<dim> >::type GridType;
-
- // Create the reference grid
- shared_ptr<GridType> referenceGrid;
-
- if (parameterSet.get<std::string>("gridType")=="structured") {
- array<unsigned int,dim> elements;
- elements.fill(numBaseElements);
- referenceGrid = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft,
- upperRight,
- elements);
- } else {
-
- if (gridFileName.rfind(".msh")!=std::string::npos)
- referenceGrid = shared_ptr<GridType>(GmshReader<GridType>::read(gridFileName));
- else
- referenceGrid = shared_ptr<GridType>(AmiraMeshReader<GridType>::read(gridFileName));
- }
-
- referenceGrid->globalRefine(numLevels-1);
-
- // Solve the rod Dirichlet problem
- SolutionType referenceSolution;
- solve(referenceGrid, referenceSolution, numLevels, parameterSet);
-
- BlockVector<TargetSpace::CoordinateType> xEmbedded(referenceSolution.size());
- for (size_t j=0; j<referenceSolution.size(); j++)
- xEmbedded[j] = referenceSolution[j].globalCoordinates();
-
- // Set up a basis of the underlying fe space
-#ifdef THIRD_ORDER
- typedef P3NodalBasis<GridType::LeafGridView,double> FEBasis;
-#elif defined SECOND_ORDER
- typedef P2NodalBasis<GridType::LeafGridView,double> FEBasis;
-#else
- typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
-#endif
- FEBasis referenceBasis(referenceGrid->leafGridView());
-
-#if !defined THIRD_ORDER && ! defined SECOND_ORDER
- VTKWriter<GridType::LeafGridView> vtkWriter(referenceGrid->leafGridView());
-
- shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> > > vtkVectorField
- = Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> > >
- (new VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> >(referenceBasis, xEmbedded, "unit vectors"));
-
- vtkWriter.addVertexData(vtkVectorField);
- vtkWriter.write("reference_result");
-#endif
-
- // //////////////////////////////////////////////////////////////////////
- // Compute mass matrix and laplace matrix to emulate L2 and H1 norms
- // //////////////////////////////////////////////////////////////////////
-
- OperatorAssembler<FEBasis,FEBasis> operatorAssembler(referenceBasis, referenceBasis);
-
- LaplaceAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> laplaceLocalAssembler;
- MassAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> massMatrixLocalAssembler;
-
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
- ScalarMatrixType laplace, massMatrix;
-
- operatorAssembler.assemble(laplaceLocalAssembler, laplace);
- operatorAssembler.assemble(massMatrixLocalAssembler, massMatrix);
-
- // ///////////////////////////////////////////////////////////
- // Compute on all coarser levels, and compare
- // ///////////////////////////////////////////////////////////
-
- std::ofstream logFile("harmonicmaps-eoc.results");
- logFile << "# mesh size, max-norm, L2-norm, h1-seminorm" << std::endl;
-
- for (int i=1; i<numLevels; i++) {
-
- shared_ptr<GridType> grid;
- if (parameterSet.get<std::string>("gridType")=="structured") {
- array<unsigned int,dim> elements;
- elements.fill(numBaseElements);
- grid = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft,
- upperRight,
- elements);
- } else {
- if (gridFileName.rfind(".msh")!=std::string::npos)
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(gridFileName));
- else
- grid = shared_ptr<GridType>(AmiraMeshReader<GridType>::read(gridFileName));
- }
-
- grid->globalRefine(i-1);
-
- // compute again
- SolutionType solution;
- solve(grid, solution, i, parameterSet);
-
- // write solution
- std::stringstream numberAsAscii;
- numberAsAscii << i;
-
- BlockVector<TargetSpace::CoordinateType> xEmbedded(solution.size());
- for (size_t j=0; j<solution.size(); j++)
- xEmbedded[j] = solution[j].globalCoordinates();
-
- FEBasis feBasis(grid->leafGridView());
- VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
-
- shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> > > vtkVectorField
- = Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> > >
- (new VTKBasisGridFunction<FEBasis,BlockVector<TargetSpace::CoordinateType> >(feBasis, xEmbedded, "unit vectors"));
-
- vtkWriter.addVertexData(vtkVectorField);
- vtkWriter.write("harmonic_result_" + numberAsAscii.str());
-
- double newL2 = GFEDifferenceNormSquared<FEBasis,TargetSpace>::compute(*referenceGrid, referenceSolution,
- *grid, solution,
- &massMatrixLocalAssembler);
- newL2 = std::sqrt(newL2);
-
- double newH1 = GFEDifferenceNormSquared<FEBasis,TargetSpace>::compute(*referenceGrid, referenceSolution,
- *grid, solution,
- &laplaceLocalAssembler);
- newH1 = std::sqrt(newH1);
-
-
- // Prolong solution to the very finest grid
- for (int j=i; j<numLevels; j++) {
- FEBasis basis(grid->leafGridView());
-#if defined THIRD_ORDER || defined SECOND_ORDER
- GeodesicFEFunctionAdaptor<FEBasis,TargetSpace>::higherOrderGFEFunctionAdaptor<order>(basis, *grid, solution);
-#else
- GeodesicFEFunctionAdaptor<FEBasis,TargetSpace>::geodesicFEFunctionAdaptor(*grid, solution);
-#endif
- }
-
- // Interpret TargetSpace as isometrically embedded into an R^m, because this is
- // how the corresponding Sobolev spaces are defined.
-
- BlockVector<TargetSpace::CoordinateType> difference(referenceSolution.size());
-
- for (size_t j=0; j<referenceSolution.size(); j++)
- difference[j] = solution[j].globalCoordinates() - referenceSolution[j].globalCoordinates();
-
- H1SemiNorm< BlockVector<TargetSpace::CoordinateType> > h1Norm(laplace);
- H1SemiNorm< BlockVector<TargetSpace::CoordinateType> > l2Norm(massMatrix);
-
- // Compute max-norm difference
- std::cout << "Level: " << i-1 << " vertices: " << xEmbedded.size() << std::endl;
- std::cout << "h: " << std::pow(0.5, i-1) << std::endl;
- std::cout << "Level: " << i-1
- << ", max-norm error: " << difference.infinity_norm()
- << std::endl;
-
- std::cout << "Level: " << i-1
- << ", L2 error: " << l2Norm(difference) << ", new: " << newL2
- << std::endl;
-
- std::cout << "Level: " << i-1
- << ", H1 error: " << h1Norm(difference) << ", new: " << newH1
- << std::endl;
-
- logFile << std::pow(0.5, i-1) << " " << difference.infinity_norm()
- << " " << l2Norm(difference)
- << " " << h1Norm(difference)
- << std::endl;
-
- }
-
- } catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }