From f826cd5c3667846cc1f7ee2356f65afb80a48660 Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Wed, 3 Feb 2016 11:47:49 +0100
Subject: [PATCH] Clean-up the whole solver a bit
---
viscoelast.cc | 576 ++++++++++++++++++++------------------------------
1 file changed, 225 insertions(+), 351 deletions(-)
diff --git a/viscoelast.cc b/viscoelast.cc
index 1046147..9c434fa 100644
--- a/viscoelast.cc
+++ b/viscoelast.cc
@@ -16,381 +16,255 @@
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/fufem/sampleonbitfield.hh>
#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/utilities/dirichletbcassembler.hh>
#include <dune/fufem/readbitfield.hh>
-#include <dune/fufem/computestress.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>
#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
+#include <dune/fufem/utilities/dirichletbcassembler.hh>
+#include <dune/functions/functionspacebases/pq1nodalbasis.hh>
-//#include <fstream>
+/** Solver for a quasi-static (viscoelastic) Kelvin-Voigt material
+ *
+ * stress = C * strain + N * strain_rate,
+ * where
+ * C Hooke tensor
+ * N viscosity tensor
+ */
+// The grid dimension
+const int dim = 3;
+using namespace Dune;
+
+int main (int argc, char *argv[]) try {
+
+ // Some types that I need
+ typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
+ typedef BlockVector<FieldVector<double, dim> > VectorType;
+
+ // parse data file
+ ParameterTree parameterSet;
+ if (argc==2)
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+ else
+ ParameterTreeParser::readINITree("viscoelast.parset", parameterSet);
+
+ // read solver settings
+ const int minLevel = parameterSet.get<int>("minLevel");
+ const int numIt = 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 baseIt = parameterSet.get<int>("baseIt");
+ const double tolerance = parameterSet.get<double>("tolerance");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+ const double timestep = parameterSet.get<double>("timestep");
+ const int num_timesteps = parameterSet.get<int>("num_timesteps");
+ const double E = parameterSet.get<double>("E");
+ const double nu = parameterSet.get<double>("nu");
+ const double nu_shear = parameterSet.get<double>("nu_shear");
+ const double nu_bulk = parameterSet.get<double>("nu_bulk");
+ const bool neumannBV = parameterSet.get("neumannBV", int(0));
+
+ // read problem settings
+ std::string path = parameterSet.get<std::string>("path");
+ std::string resultpath = parameterSet.get<std::string>("resultpath");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+
+ // /////////////////////////////
+ // Generate the grid
+ // /////////////////////////////
+
+ typedef UGGrid<dim> GridType;
+ typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
+ typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
+
+ GridType* grid;
+ if (parameterSet.hasKey("parFile")) {
+ std::string parFile = parameterSet.get<std::string>("parFile");
+ grid = AmiraMeshReader<GridType>::read(path + gridFile, PSurfaceBoundary<dim-1>::read(path + parFile));
+ } else
+ grid = AmiraMeshReader<GridType>::read(path + gridFile);
+ grid->setRefinementType(GridType::COPY);
+ // Read Dirichlet boundary values
+ BitSetVector<dim> dirichletNodes;
+ VectorType dirichletValues;
+ DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,parameterSet,dirichletNodes,dirichletValues,path);
+ //Read neumann boundary values
+ LevelBoundaryPatch coarseNeumannBoundary;
+ VectorType coarseNeumannValues;
+ if (neumannBV) {
-//Using Maxwell-Voigt consititutive law for viscoelastic materials : stress = C * strain + N * strain_rate
-// where C elasticitytensor
-// N viscositytensor
+ std::string neumannFile = parameterSet.get<std::string>("neumannFile");
+ std::string neumannValuesFile = parameterSet.get<std::string>("neumannValuesFile");
+ readBoundaryPatch<GridType>(coarseNeumannBoundary, path + neumannFile);
+ coarseNeumannValues.resize(grid->size(0,dim));
+ AmiraMeshReader<GridType>::readFunction(coarseNeumannValues, path + neumannValuesFile);
+ }
-// The grid dimension
-const int dim = 3;
+ // refine uniformly until minlevel
+ for (int i=0; i<minLevel; i++)
+ grid->globalRefine(1);
-using namespace Dune;
+ // initial solution
+ VectorType x = dirichletValues;
+
+ //initial velocity needed for the computation of the stress(->strainrate)
+ VectorType v(grid->size(dim));
+ v = 0;
+
+ //Determine Neumann dofs
+
+ BitSetVector<dim> neumannNodes(grid->size(dim));
+ VectorType neumannValues(neumannNodes.size());
+
+ BoundaryPatch<GridType::LeafGridView> neumannBoundary;
+ if(neumannBV) {
+ BitSetVector<1> scalarNeumannNodes;
+ DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,*coarseNeumannBoundary.getVertices(),
+ coarseNeumannValues,scalarNeumannNodes,
+ neumannValues);
+ neumannBoundary.setup(grid->leafGridView(),scalarNeumannNodes);
+
+ for (size_t j=0; j<neumannNodes.size(); j++)
+ for (int k=0; k<dim; k++)
+ neumannNodes[j][k] = neumannBoundary.containsVertex(j);
+ }
+
+
+ //right-hand-side rhs = bodyforces + neumanboundary
+ VectorType rhs(grid->size(dim));
+ rhs = 0;
+
+ //right-hand-side of the implicit euler scheme N*d_k + timestep*rhs N viscosity_stiffness_matrix
+ VectorType erhs = rhs;
+
+ MatrixType elasticPart,viscousPart;
+
+ using DuneP1Basis = Dune::Functions::PQ1NodalBasis<typename GridType::LeafGridView>;
+ using P1Basis = DuneFunctionsBasis<DuneP1Basis>;
+ P1Basis p1Basis(grid->leafGridView());
+ OperatorAssembler<P1Basis,P1Basis> globalAssembler(p1Basis,p1Basis);
+ StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(E, nu);
+ globalAssembler.assemble(localAssembler, elasticPart);
+
+ ViscosityAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> visLocalAssembler(nu_shear, nu_bulk);
+ globalAssembler.assemble(visLocalAssembler, viscousPart);
+
+ if (neumannBV) {
+ BasisGridFunction<P1Basis,VectorType> neumannFunction(p1Basis, neumannValues);
+ NeumannBoundaryAssembler<GridType, Dune::FieldVector<double,dim> > neumannBoundaryAssembler(neumannFunction);
+
+ BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(p1Basis, neumannBoundary);
+ boundaryFunctionalAssembler.assemble(neumannBoundaryAssembler, rhs, true);
+ }
+
+ /////////////////////////////
+ // Create a solver
+ // ///////////////////////////
+
+ // First create a base solver
+ BlockGSStep<MatrixType, VectorType> baseSolverStep;
+
+ EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep);
+
+ ::LoopSolver<VectorType> baseSolver(&baseSolverStep,
+ baseIt,
+ baseTolerance,
+ &baseEnergyNorm,
+ Solver::QUIET);
+
+ baseSolver.verbosity_ = Solver::QUIET;
+ baseSolver.tolerance_ = baseTolerance;
+
+ // Make pre and postsmoothers
+ BlockGSStep<MatrixType, VectorType> presmoother;
+ BlockGSStep<MatrixType, VectorType> postsmoother;
+
+
+
+ //FEM leads to:
+ //Quasistatic equation of the form : N*D'(t)+K*D(t)=rhs(t) N viscous_stiffness K elastic_stiffness D displacement
+ //Using implicit euler scheme leads to: (N+timestep*K)*D_{k+1}=N*D_{k}+timestep*rhs(t_k+1)
+
+ MatrixType lgs = viscousPart;
+ elasticPart *= timestep;
+ lgs += elasticPart;
+ rhs *= timestep;
+
+ for (int j=1;j<=num_timesteps;j++) {
+
+ // time[j-1]=j*timestep;
+
+ //N*D_{0}=0 because x=0
+ erhs += rhs ;
+
+ MultigridStep<MatrixType, VectorType> multigridStep(lgs, x, erhs);
+
+ multigridStep.setMGType(mu, nu1, nu2);
+ multigridStep.ignoreNodes_ = &dirichletNodes;
+ multigridStep.basesolver_ = &baseSolver;
+ multigridStep.setSmoother(&presmoother,&postsmoother);
+
+ std::vector<CompressedMultigridTransfer<VectorType>* > transfers(grid->maxLevel());
+ for (size_t i=0; i<transfers.size(); i++) {
+ transfers[i] = new CompressedMultigridTransfer<VectorType>;
+ transfers[i]->setup(*grid, i, i+1);
+ }
+ multigridStep.setTransferOperators(transfers);
+
+ EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep);
+
+ ::LoopSolver<VectorType> solver(&multigridStep,
+ numIt,
+ tolerance,
+ &energyNorm,
+ Solver::FULL);
+
+ solver.preprocess();
+
+ // Compute solution
+ solver.solve();
+
+ x = multigridStep.getSol();
+
+ std::string name= "resultGrid";
+
+ // make string from int (always four digits, with leading zeros)
+ std::stringstream numberString;
+ numberString << std::setw(4) << std::setfill('0') << j;
+ std::string num = numberString.str();
+
+ name+=num;
+ LeafAmiraMeshWriter<GridType> amiramesh;
+ amiramesh.addLeafGrid(*grid,true);
+ amiramesh.addVertexData(x, grid->leafGridView());
+ amiramesh.write(resultpath + name,1);
+
+ viscousPart.mv(x,erhs);
+
+ v+=x;
+ v*=(1/timestep);
+ }
+} catch (Exception e) {
+ std::cout << e << std::endl;
+}
-int main (int argc, char *argv[]) try
- {
- // Some types that I need
- typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
- typedef BlockVector<FieldVector<double, dim> > VectorType;
-
- // parse data file
- ParameterTree parameterSet;
- if (argc==2)
- ParameterTreeParser::readINITree(argv[1], parameterSet);
- else
- ParameterTreeParser::readINITree("viscoelast.parset", parameterSet);
-
- // read solver settings
- const int minLevel = parameterSet.get<int>("minLevel");
- const int maxLevel = parameterSet.get<int>("maxLevel");
- const int numIt = 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 baseIt = parameterSet.get<int>("baseIt");
- const double tolerance = parameterSet.get<double>("tolerance");
- const double baseTolerance = parameterSet.get<double>("baseTolerance");
- const bool nestedIteration = parameterSet.get<int>("nestedIteration");
- const bool paramBoundaries = parameterSet.get<int>("paramBoundaries");
- const double timestep = parameterSet.get<double>("timestep");
- const int num_timesteps = parameterSet.get<int>("num_timesteps");
- const double threshold = parameterSet.get<double>("threshold");
- const double E = parameterSet.get<double>("E");
- const double nu = parameterSet.get<double>("nu");
- const double nu_shear = parameterSet.get<double>("nu_shear");
- const double nu_bulk = parameterSet.get<double>("nu_bulk");
- const bool neumannBV = parameterSet.get("neumannBV", int(0));
-
-
-
- // read problem settings
- std::string path = parameterSet.get<std::string>("path");
- std::string resultpath = parameterSet.get<std::string>("resultpath");
- std::string gridFile = parameterSet.get<std::string>("gridFile");
- std::string parFile = parameterSet.get("parFile", "");
- std::string dirichletFile = parameterSet.get<std::string>("dirichletFile");
- std::string dirichletValuesFile = parameterSet.get<std::string>("dirichletValuesFile");
- std::string neumannFile = parameterSet.get("neumannFile", "");
- std::string neumannValuesFile = parameterSet.get("neumannValuesFile", "");
-
- /*
- //Only needed for computation of surface stress on the whole boundary
- std::string dirichletFileAll = parameterSet.get<std::string>("dirichletFileAll");
- */
-
- // /////////////////////////////
- // Generate the grid
- // /////////////////////////////
-
- typedef UGGrid<dim> GridType;
- typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
- typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
-
- GridType grid;
-
- grid.setRefinementType(GridType::COPY);
-
-
- if (paramBoundaries)
- Dune::AmiraMeshReader<GridType>::read(grid, path + gridFile, path + parFile);
- else
- Dune::AmiraMeshReader<GridType>::read(grid, path + gridFile);
-
-
- // Read Dirichlet boundary values
- std::vector<BitSetVector<dim> > dirichletNodes(maxLevel+1);
-
- std::vector<LevelBoundaryPatch> dirichletBoundary(maxLevel+1);
- dirichletBoundary[0].setup(grid.levelGridView(0));
- readBoundaryPatch<GridType>(dirichletBoundary[0], path + dirichletFile);
-
- std::vector<VectorType> dirichletValues(maxLevel+1);
- dirichletValues[0].resize(grid.size(0, dim));
- AmiraMeshReader<GridType>::readFunction(dirichletValues[0], path + dirichletValuesFile);
-
- /*
- //Create boundary patch for whole boundary (only needed for computation of surface stress on the whole boundary)
-
- std::vector<BitField> dirichletNodesAll(maxLevel+1);
- readBitField(dirichletNodesAll[0], grid.size(0, dim), path + dirichletFileAll);
-
- std::vector<BoundaryPatch<GridType> > dirichletBoundaryAll(maxLevel+1);
- dirichletBoundaryAll[0].setup(grid, 0, dirichletNodesAll[0]);
- */
-
- LevelBoundaryPatch coarseNeumannBoundary;
- VectorType coarseNeumannValues;
-
- //Read neumann boundary values
- if (neumannBV) {
-
- coarseNeumannBoundary.setup(grid.levelGridView(0));
- readBoundaryPatch<GridType>(coarseNeumannBoundary, path + neumannFile);
- coarseNeumannValues.resize(grid.size(0,dim));
- AmiraMeshReader<GridType>::readFunction(coarseNeumannValues, path + neumannValuesFile);
- }
-
-
-
- // refine uniformly until minlevel
- for (int i=0; i<minLevel; i++)
- grid.globalRefine(1);
-
- // initial solution
- VectorType x(grid.size(dim));
- x = 0;
-
- //initial velocity -in this quasistatic equation the velocity is needed for computation of the stress(->strainrate)
- VectorType v(grid.size(dim));
- v = 0;
-
-
-
- // //////////////////////////////////////////////////
- // Refinement loop
- // //////////////////////////////////////////////////
-
-
-
- // Determine Dirichlet dofs
- BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);
-
- for (int i=0; i<=grid.maxLevel(); i++) {
-
- dirichletNodes[i].resize(grid.size(i,dim)*dim);
- for (int j=0; j<grid.size(i,dim); j++)
- dirichletNodes[i][j] = dirichletBoundary[i].containsVertex(j);
-
- }
-
- sampleOnBitField(grid, dirichletValues, dirichletNodes);
-
- //Determine Neumann dofs
-
- BitSetVector<dim> neumannNodes(grid.size(dim));
- VectorType neumannValues(neumannNodes.size());
- BoundaryPatch<GridType::LeafGridView> neumannBoundary(grid.leafGridView());
- if(neumannBV) {
- BoundaryPatchProlongator<GridType>::prolong(coarseNeumannBoundary,neumannBoundary);
-
- for (size_t j=0; j<neumannNodes.size(); j++)
- for (int k=0; k<dim; k++)
- neumannNodes[j][k] = neumannBoundary.containsVertex(j);
-
- sampleOnBitField(grid, coarseNeumannValues, neumannValues, neumannNodes);
- }
-
-
- //right-hand-side rhs = bodyforces + neumanboundary
- VectorType rhs(grid.size(dim));
- rhs = 0;
-
- //right-hand-side of the implicit euler scheme N*d_k + timestep*rhs N viscosity_stiffness_matrix
- VectorType erhs(grid.size(dim));
- erhs = 0;
-
- for (int i=0; i<rhs.size(); i++)
- for (int j=0; j<dim; j++) {
- if (dirichletNodes[grid.maxLevel()][i][j])
- x[i][j] = dirichletValues[grid.maxLevel()][i][j];
-
- }
-
-
- OperatorType elasticPart,viscousPart;
- typedef P1NodalBasis<GridType::LeafGridView,double> P1Basis;
- P1Basis p1Basis(grid.leafGridView());
- OperatorAssembler<P1Basis,P1Basis> globalAssembler(p1Basis,p1Basis);
- StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(E, nu);
- globalAssembler.assemble(localAssembler, elasticPart);
-
- ViscosityAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> visLocalAssembler(nu_shear, nu_bulk);
- globalAssembler.assemble(visLocalAssembler, viscousPart);
-
- if (neumannBV) {
- BasisGridFunction<P1Basis,VectorType> neumannFunction(p1Basis, neumannValues);
- NeumannBoundaryAssembler<GridType, Dune::FieldVector<double,dim> > neumannBoundaryAssembler(neumannFunction);
-
- BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(p1Basis, neumannBoundary);
- boundaryFunctionalAssembler.assemble(neumannBoundaryAssembler, rhs, true);
- }
- /*
- //test vectors with time ,total surfacestress and total strain data
- BlockVector<FieldVector<double,1> > time(num_timesteps),surfacestress(num_timesteps), strain(num_timesteps);
- time=0.0;
- surfacestress=0.0;
- strain=0.0;
- */
-
-
-
-
- /////////////////////////////
- // Create a solver
- // ///////////////////////////
-
- // First create a base solver
- BlockGSStep<OperatorType, VectorType> baseSolverStep;
-
- EnergyNorm<OperatorType, VectorType> baseEnergyNorm(baseSolverStep);
-
- ::LoopSolver<VectorType> baseSolver(&baseSolverStep,
- baseIt,
- baseTolerance,
- &baseEnergyNorm,
- Solver::QUIET);
-
- baseSolver.verbosity_ = Solver::QUIET;
- baseSolver.tolerance_ = baseTolerance;
-
- // Make pre and postsmoothers
- BlockGSStep<OperatorType, VectorType> presmoother;
- BlockGSStep<OperatorType, VectorType> postsmoother;
-
-
-
- //FEM leads to:
- //Quasistatic equation of the form : N*D'(t)+K*D(t)=rhs(t) N viscous_stiffness K elastic_stiffness D displacement
- //Using implicit euler scheme leads to: (N+timestep*K)*D_{k+1}=N*D_{k}+timestep*rhs(t_k+1)
-
- OperatorType lgs = viscousPart;
- elasticPart *= timestep;
- lgs += elasticPart;
- rhs *= timestep;
-
- for (int j=1;j<=num_timesteps;j++) {
-
- // time[j-1]=j*timestep;
-
- //N*D_{0}=0 because x=0
- erhs += rhs ;
-
- MultigridStep<OperatorType, VectorType> multigridStep(lgs, x, erhs);
-
- multigridStep.setMGType(1, nu1, nu2);
-
- multigridStep.ignoreNodes_ = &dirichletNodes.back();
-
- multigridStep.basesolver_ = &baseSolver;
- multigridStep.setSmoother(&presmoother,&postsmoother);
-
- multigridStep.mgTransfer_.resize(grid.maxLevel());
- for (int i=0; i<multigridStep.mgTransfer_.size(); i++) {
- CompressedMultigridTransfer<VectorType>* newTransferOp = new CompressedMultigridTransfer<VectorType>;
- newTransferOp->setup(grid, i, i+1);
- multigridStep.mgTransfer_[i] = newTransferOp;
- }
-
- EnergyNorm<OperatorType, VectorType> energyNorm(multigridStep);
-
- ::LoopSolver<VectorType> solver(&multigridStep,
- numIt,
- tolerance,
- &energyNorm,
- Solver::FULL);
-
-
-
- solver.preprocess();
- multigridStep.preprocess();
-
- // Compute solution
- solver.solve();
-
- x = multigridStep.getSol();
-
- std::string name= "resultGrid";
-
- // make string from int (always four digits, with leading zeros)
- std::stringstream numberString;
- numberString << std::setw(4) << std::setfill('0') << j;
- std::string num = numberString.str();
-
- name+=num;
- LeafAmiraMeshWriter<GridType> amiramesh;
- amiramesh.addLeafGrid(grid,true);
- amiramesh.addVertexData(x, grid.leafGridView());
- amiramesh.write(resultpath + name,1);
-
-
-
- viscousPart.mv(x,erhs);
-
-
- v+=x;
- v*=(1/timestep);
-
- //compute von mises stress for each element at its center of gravity
- BlockVector<FieldVector<double,1> > stress;
- std::string namestress= "resultGridstress";
- namestress+=num;
- Stress<GridType::LevelGridView>::getViscoelasticStress(grid, x, v, stress, E, nu, nu_bulk, nu_shear);
- LeafAmiraMeshWriter<GridType>::writeBlockVector(grid, stress, resultpath + namestress);
-
-
-
- /*
- //Compute the total surface stress
-
- surfacestress[j-1]=Stress<GridType::LevelGridView,dim>::getViscoelasticSurfaceStress_interpNormals(dirichletBoundaryAll[grid.maxLevel()],x,v,E, nu, nu_bulk, nu_shear).two_norm();
-
- //Compute total strain
- Stress<GridType::LevelGridView,dim>::getTotalStrain(grid,x,strain[j-1]);
- */
-
-
- v=x;
- v*=-1;
-
-
- }
-
- /*
- //save stress data in stress.dat file in format which can be easily plotted using matlab
- std::fstream file;
- file.open("stress.dat", std::ios::out);
- for (int z=0;z<num_timesteps;z++){
- file << time[z] << " " << surfacestress[z] << std::endl;
- }
- file.close();
-
- std::fstream strainfile;
- strainfile.open("strain.dat", std::ios::out);
- for (int z=0;z<num_timesteps;z++){
- strainfile << time[z] << " " << strain[z] << std::endl;
- }
- strainfile.close();
- */
-
- } catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }
-
--
GitLab