From fadc99d3739ee880f4f50be9ff63b70c73737c21 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 14 Jun 2019 16:51:21 +0200
Subject: [PATCH] Remove the program rod-eoc
It combined computing rod configurations and their discretization
errors. Nowadays the way to do that is to use rod3d.cc for the
rod configurations, and compute-disc-errors.cc for the errors.
Hence rod-eoc.cc is not needed anymore and can go.
---
src/CMakeLists.txt | 3 +-
src/rod-eoc.cc | 248 ---------------------------------------------
2 files changed, 1 insertion(+), 250 deletions(-)
delete mode 100644 src/rod-eoc.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 98aea9b5..15adcea1 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -3,8 +3,7 @@ set(programs compute-disc-error
film-on-substrate
gradient-flow
harmonicmaps
- rod3d
- rod-eoc)
+ rod3d)
# rodobstacle)
foreach(_program ${programs})
diff --git a/src/rod-eoc.cc b/src/rod-eoc.cc
deleted file mode 100644
index 5926b97b..00000000
--- a/src/rod-eoc.cc
+++ /dev/null
@@ -1,248 +0,0 @@
-#include <config.h>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#include <dune/grid/onedgrid.hh>
-
-#include <dune/istl/io.hh>
-
-#include <dune/functions/functionspacebases/lagrangebasis.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/norms/energynorm.hh>
-
-#include <dune/gfe/rigidbodymotion.hh>
-#include <dune/gfe/geodesicdifference.hh>
-#include <dune/gfe/rotation.hh>
-#include <dune/gfe/rodassembler.hh>
-#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/geodesicfefunctionadaptor.hh>
-#include <dune/gfe/rodwriter.hh>
-
-typedef Dune::OneDGrid GridType;
-
-typedef RigidBodyMotion<double,3> TargetSpace;
-typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
-
-const int blocksize = TargetSpace::TangentVector::dimension;
-
-using namespace Dune;
-using std::string;
-
-void solve (const GridType& grid,
- SolutionType& x,
- int numLevels,
- const TargetSpace& dirichletValue,
- 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 rod parameter settings
- const double A = parameters.get<double>("A");
- const double J1 = parameters.get<double>("J1");
- const double J2 = parameters.get<double>("J2");
- const double E = parameters.get<double>("E");
- const double nu = parameters.get<double>("nu");
-
- // Create a function space basis
- typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView, 1> FEBasis;
- FEBasis feBasis(grid.leafGridView());
-
- // Transitional: the same basis as a dune-fufem object
- typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
- FufemFEBasis fufemFeBasis(feBasis);
-
- // Create a local assembler
- RodLocalStiffness<OneDGrid::LeafGridView,double> localStiffness(grid.leafGridView(),
- A, J1, J2, E, nu);
-
-
-
-
- x.resize(grid.size(1));
-
- // //////////////////////////
- // Initial solution
- // //////////////////////////
-
- for (size_t i=0; i<x.size(); i++) {
- x[i].r[0] = 0;
- x[i].r[1] = 0;
- x[i].r[2] = double(i)/(x.size()-1);
- x[i].q = Rotation<double,3>::identity();
- }
-
- // set Dirichlet value
- x.back() = dirichletValue;
-
- // Both ends are Dirichlet
- BitSetVector<blocksize> dirichletNodes(grid.size(1));
- dirichletNodes.unsetAll();
-
- dirichletNodes[0] = dirichletNodes.back() = true;
-
- // ///////////////////////////////////////////
- // Create a solver for the rod problem
- // ///////////////////////////////////////////
-
- RodAssembler<FEBasis,3> rodAssembler(feBasis, &localStiffness);
-
- RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
-#if 1
- rodSolver.setup(grid,
- &rodAssembler,
- x,
- dirichletNodes,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- innerTolerance,
- 1, 3, 3,
- 100, // iterations of the base solver
- 1e-8, // base tolerance
- false); // instrumentation
-#else
- rodSolver.setupTCG(grid,
- &rodAssembler,
- x,
- dirichletNodes,
- tolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- innerTolerance,
- false); // instrumentation
-#endif
-
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
-
- rodSolver.setInitialIterate(x);
- rodSolver.solve();
-
- x = rodSolver.getSol();
-}
-
-int main (int argc, char *argv[]) try
-{
- // parse data file
- ParameterTree parameterSet;
- if (argc==2)
- ParameterTreeParser::readINITree(argv[1], parameterSet);
- else
- ParameterTreeParser::readINITree("rod-eoc.parset", parameterSet);
-
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
- const int numRodBaseElements = parameterSet.get<int>("numRodBaseElements");
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
-
- RigidBodyMotion<double,3> dirichletValue;
- dirichletValue.r[0] = parameterSet.get<double>("dirichletValueX");
- dirichletValue.r[1] = parameterSet.get<double>("dirichletValueY");
- dirichletValue.r[2] = parameterSet.get<double>("dirichletValueZ");
-
- FieldVector<double,3> axis;
- axis[0] = parameterSet.get<double>("dirichletAxisX");
- axis[1] = parameterSet.get<double>("dirichletAxisY");
- axis[2] = parameterSet.get<double>("dirichletAxisZ");
- double angle = parameterSet.get<double>("dirichletAngle");
-
- dirichletValue.q = Rotation<double,3>(axis, M_PI*angle/180);
-
- // ///////////////////////////////////////////////////////////
- // First compute the 'exact' solution on a very fine grid
- // ///////////////////////////////////////////////////////////
-
- // Create the reference grid
-
- GridType referenceGrid(numRodBaseElements, 0, 1);
-
- referenceGrid.globalRefine(numLevels-1);
-
- // Solve the rod Dirichlet problem
- SolutionType referenceSolution;
- solve(referenceGrid, referenceSolution, numLevels, dirichletValue, parameterSet);
-
-
- // //////////////////////////////////////////////////////////////////////
- // Compute mass matrix and laplace matrix to emulate L2 and H1 norms
- // //////////////////////////////////////////////////////////////////////
-
- typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
- FEBasis basis(referenceGrid.leafGridView());
- OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
-
- 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
- // ///////////////////////////////////////////////////////////
-
- for (int i=1; i<=numLevels; i++) {
-
- GridType grid(numRodBaseElements, 0, 1);
- grid.globalRefine(i-1);
-
- // compute again
- SolutionType solution;
- solve(grid, solution, i, dirichletValue, parameterSet);
-
- // Prolong solution to the very finest grid
- for (int j=i; j<numLevels; j++)
- GeodesicFEFunctionAdaptor<FEBasis,TargetSpace>::geodesicFEFunctionAdaptor(grid, solution);
-
- std::stringstream numberAsAscii;
- numberAsAscii << i;
- writeRod(solution, "rodGrid_" + numberAsAscii.str());
-
- assert(referenceSolution.size() == solution.size());
-
- BlockVector<TargetSpace::TangentVector> difference = computeGeodesicDifference(solution,referenceSolution);
-
- H1SemiNorm< BlockVector<TargetSpace::TangentVector> > h1Norm(laplace);
- H1SemiNorm< BlockVector<TargetSpace::TangentVector> > l2Norm(massMatrix);
-
- // Compute max-norm difference
- std::cout << "Level: " << i-1
- << ", max-norm error: " << difference.infinity_norm()
- << std::endl;
-
- std::cout << "Level: " << i-1
- << ", L2 error: " << l2Norm(difference)
- << std::endl;
-
- std::cout << "Level: " << i-1
- << ", H1 error: " << h1Norm(difference)
- << std::endl;
-
- }
-
-} catch (Exception& e)
-{
- std::cout << e.what() << std::endl;
-}
--
GitLab