-
Jonathan Youett authoredJonathan Youett authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
viscoelast.cc 9.87 KiB
#include <config.h>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/amirameshwriter.hh>
#include <dune/grid/io/file/amirameshreader.hh>
#include <dune/istl/io.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.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/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/lagrangebasis.hh>
/** 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;
std::unique_ptr<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) {
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);
}
// refine uniformly until minlevel
for (int i=0; i<minLevel; i++)
grid->globalRefine(1);
// 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::LagrangeBasis<typename GridType::LeafGridView, 1>;
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
auto baseSolverStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(
Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep);
::LoopSolver<VectorType> baseSolver(baseSolverStep,
baseIt,
baseTolerance,
baseEnergyNorm,
Solver::QUIET);
// Make pre and postsmoothers
auto presmoother = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(
Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
auto postsmoother = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(
Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
//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.setIgnore(dirichletNodes);
multigridStep.setBaseSolver(baseSolver);
multigridStep.setSmoother(&presmoother,&postsmoother);
std::vector<std::shared_ptr<CompressedMultigridTransfer<VectorType> > > transfers(grid->maxLevel());
for (size_t i=0; i<transfers.size(); i++) {
transfers[i] = std::make_shared<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;
}