Select Git revision
dune-elasticity.pc.in
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
viscoelast.cc 13.61 KiB
#include <config.h>
#include <dune/common/bitsetvector.hh>
#include <dune/common/configparser.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/disc/operators/p1operator.hh>
#include <dune/ag-common/sampleonbitfield.hh>
#include <dune/ag-common/prolongboundarypatch.hh>
#include <dune/ag-common/readbitfield.hh>
#include <dune/ag-common/blockgsstep.hh>
#include <dune/ag-common/multigridstep.hh>
#include <dune/ag-common/transferoperators/compressedmultigridtransfer.hh>
#include <dune/ag-common/iterativesolver.hh>
#include <dune/ag-common/computestress.hh>
#include <dune/ag-common/neumannassembler.hh>
#include <dune/ag-common/norms/energynorm.hh>
#include <dune/ag-common/hierarchicestimator.hh>
#include <dune/disc/elasticity/linearelasticityassembler.hh>
#include "src/viscoelasticassembler.hh"
//#include <fstream>
//Using Maxwell-Voigt consititutive law for viscoelastic materials : stress = C * strain + N * strain_rate
// where C elasticitytensor
// N viscositytensor
// 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> > OperatorType;
typedef BlockVector<FieldVector<double, dim> > VectorType;
// parse data file
ConfigParser parameterSet;
parameterSet.parseFile("viscoelast.parset");
// 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;
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<BoundaryPatch<GridType> > dirichletBoundary(maxLevel+1);
dirichletBoundary[0].setup(grid, 0);
readBoundaryPatch(dirichletBoundary[0], path + dirichletFile);
std::vector<VectorType> dirichletValues(maxLevel+1);
dirichletValues[0].resize(grid.size(0, dim));
AmiraMeshReader<int>::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]);
*/
std::vector<BitSetVector<dim> > neumannNodes(maxLevel+1);
std::vector<BoundaryPatch<GridType> > neumannBoundary(maxLevel+1);
std::vector<VectorType> neumannValues(maxLevel+1);
//Read neumann boundary values
if (neumannBV) {
neumannBoundary[0].setup(grid, 0);
readBoundaryPatch(neumannBoundary[0], path + neumannFile);
neumannValues[0].resize(grid.size(0,dim));
AmiraMeshReader<int>::readFunction(neumannValues[0], 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
PatchProlongator<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
if(neumannBV){
PatchProlongator<GridType>::prolong(neumannBoundary);
for (int i=0; i<=grid.maxLevel(); i++) {
neumannNodes[i].resize(grid.size(i,dim));
for (int j=0; j<grid.size(i,dim); j++)
for (int k=0; k<dim; k++)
neumannNodes[i][j][k] = neumannBoundary[i].containsVertex(j);
}
sampleOnBitField(grid, neumannValues, neumannNodes);
}
//PatchProlongator<GridType>::prolong(dirichletBoundaryAll);
//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];
}
// Assemble neumann boundary values
if (neumannBV)
assembleAndAddNeumannTerm(neumannBoundary[grid.maxLevel()],neumannValues[grid.maxLevel()],rhs);
// Assemble elastic part of the problem
LeafP1Function<GridType,double,dim> u(grid),f(grid);
LinearElasticityLocalStiffness<GridType::LeafGridView,double> lstiff(E, nu);
LeafP1OperatorAssembler<GridType,double,dim> elasticbilinearForm(grid);
elasticbilinearForm.assemble(lstiff,u,f);
//Assemble viscous part of the problem
LinearViscosityLocalStiffness<GridType::LeafGridView,double> vstiff(nu_shear,nu_bulk);
LeafP1OperatorAssembler<GridType,double,dim> viscousbilinearForm(grid);
viscousbilinearForm.assemble(vstiff,u,f);
/*
//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 = (*viscousbilinearForm);
*elasticbilinearForm *= timestep;
lgs += *elasticbilinearForm;
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, grid.maxLevel()+1);
multigridStep.setMGType(1, nu1, nu2);
multigridStep.ignoreNodes_ = &dirichletNodes.back();
multigridStep.basesolver_ = &baseSolver;
multigridStep.presmoother_ = &presmoother;
multigridStep.postsmoother_ = &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.leafView());
amiramesh.write(resultpath + name,1);
(*viscousbilinearForm).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,dim>::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,dim>::getViscoelasticSurfaceStress_interpNormals(dirichletBoundaryAll[grid.maxLevel()],x,v,E, nu, nu_bulk, nu_shear).two_norm();
//Compute total strain
Stress<GridType,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;
}