viscoelast.cc

#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>


/** 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;

    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::PQkNodalBasis<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;

}