Skip to content
Snippets Groups Projects
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nonlinelast.cc 14.53 KiB
#include <config.h>

#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/bitsetvector.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/fufem/functionspacebases/p2nodalbasis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
#include <dune/fufem/estimators/hierarchicalestimator.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/estimators/fractionalmarking.hh>
#include <dune/fufem/estimators/refinementindicator.hh>
#include <dune/fufem/utilities/dirichletbcassembler.hh>
#include <dune/fufem/utilities/gridconstruction.hh>

#ifdef HAVE_IPOPT
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif

#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/solvers/trustregionsolver.hh>

#include <dune/elasticity/common/nonlinearelasticityproblem.hh>
#define LAURSEN
#include <dune/elasticity/materials/neohookeanmaterial.hh>
#include <dune/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>

// The grid dimension
const int dim = 3;
typedef double field_type;
using namespace Dune;

int main (int argc, char *argv[]) try
{
    Dune::MPIHelper::instance(argc, argv);
    // 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("nonlinelast.parset", parameterSet);

    // read solver settings
    const int minLevel         = parameterSet.get<int>("minLevel");
    const int maxLevel         = parameterSet.get<int>("maxLevel");
    const bool useH1SemiNorm   = parameterSet.get<bool>("useH1SemiNorm");

    // read problem settings
    std::string path           = parameterSet.get<std::string>("path");
    std::string resultPath     = parameterSet.get<std::string>("resultPath");

    // /////////////////////////////
    //   Generate the grid
    // /////////////////////////////
    typedef UGGrid<dim> GridType;
    std::unique_ptr<GridType> grid;

    const auto& gridConfig = parameterSet.sub(parameterSet.get<std::string>("gridName"));

    if (gridConfig.get<bool>("createGrid", false)) {
      grid = GridConstruction<GridType,dim>::createGrid(gridConfig);
    } else if (gridConfig.hasKey("amira_parFile")) {
#if HAVE_AMIRAMESH
        auto gridFile = gridConfig.get<std::string>("gridFile");
        auto parFile = gridConfig.get<std::string>("amira_parFile");
        grid = AmiraMeshReader<GridType>::read(path + gridFile, PSurfaceBoundary<dim-1>::read(path + parFile));
#endif
    } else {
#if HAVE_AMIRAMESH
      grid = AmiraMeshReader<GridType>::read(path + gridConfig.get<std::string>("gridFile"));
#endif
    }

    // Read coarse Dirichlet boundary condition
    BitSetVector<dim> coarseDirichletNodes(grid->size(dim),false);
    VectorType coarseDirichletValues(grid->size(dim));
    coarseDirichletValues = 0;

    for (int i = 0; i < gridConfig.get<int>("nDirichletConditions"); ++i) {
        const auto& dirConfig = gridConfig.sub(formatString("dirichlet%d",i));
        BitSetVector<dim> dirichletNodes;
        VectorType dirichletValues;
        DirichletBCAssembler<GridType>::assembleDirichletBC(*grid, dirConfig,
                                    dirichletNodes, dirichletValues, path);
        auto scaling = dirConfig.get<double>("dvScaling", 1);

        // combine with previous ones
        coarseDirichletValues.axpy(scaling, dirichletValues);
        for (size_t j = 0; j < dirichletNodes.size(); ++j)
            for (int k = 0; k < dim; ++k)
                if (dirichletNodes[j][k])
                    coarseDirichletNodes[j][k] = true;
    }

    //  Uniform refine grid
    grid->setRefinementType(GridType::COPY);
    for (int i=0; i<minLevel; i++)
        grid->globalRefine(1);

    // Initial solution
    VectorType x(grid->size(dim));
    x = 0;

    // /////////////////////////////////////
    //setup the monotone multigrid step
    // /////////////////////////////////////

    MonotoneMGStep<MatrixType,VectorType> mmgStep;
    mmgStep.setMGType(parameterSet.get("mu",1),
            parameterSet.get("preSmoothingSteps",3),
            parameterSet.get("postSmoothingSteps",3));

    TrustRegionGSStep<MatrixType, VectorType> presmoother,postsmoother;
    mmgStep.setSmoother(&presmoother, &postsmoother);

    mmgStep.setObstacleRestrictor(MandelObstacleRestrictor<VectorType>{});
    mmgStep.setVerbosity(parameterSet.get<NumProc::VerbosityMode>("verbosity"));

    // Create a base solver
#ifdef HAVE_IPOPT
    QuadraticIPOptSolver<MatrixType,VectorType> baseSolver(parameterSet.get<field_type>("baseTolerance"),
                                                           100, NumProc::QUIET);
#endif
    mmgStep.setBaseSolver(baseSolver);

    EnergyNorm<MatrixType, VectorType> h1SemiNorm;
    EnergyNorm<MatrixType, VectorType  > energyNorm(mmgStep);

    ::LoopSolver<VectorType> solver(mmgStep,
            parameterSet.get<int>("mgIterations"),
            parameterSet.get<field_type>("mgTolerance"),
            (useH1SemiNorm) ? h1SemiNorm : energyNorm,
            Solver::FULL);


    // /////////////////////////////////
    // Setup the nonlinear material
    // ////////////////////////////////

    typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
    P1Basis p1Basis(grid->leafGridView());

#if HAVE_ADOLC
    MooneyRivlinMaterial<P1Basis> localEnergy(p1Basis,
                          parameterSet.get<double>("E"),
                          parameterSet.get<double>("nu"));
    using MaterialType = AdolcMaterial<P1Basis>;
    MaterialType material(p1Basis, localEnergy, parameterSet.get<bool>("vectorMode"));
#else
    using MaterialType = NeoHookeanMaterial<P1Basis>;
    MaterialType material(p1Basis,
                          parameterSet.get<double>("E"),
                          parameterSet.get<double>("nu"));
#endif

    // ///////////////////////////////////////////////////
    //   Do a homotopy of the Dirichlet boundary data
    // ///////////////////////////////////////////////////

    field_type loadFactor    = 0;
    field_type loadIncrement = parameterSet.get<double>("loadIncrement");

    for (int level=minLevel; level<=maxLevel; level++) {


        VectorType extForces(grid->size(dim));
        extForces = 0;

        // Dirichlet values on the fine grid
        Dune::BitSetVector<dim> dirichletNodes;
        VectorType dirichletValues;

        DirichletBCAssembler<GridType>::assembleDirichletBC(*grid,
                                coarseDirichletNodes,coarseDirichletValues,
                                dirichletNodes, dirichletValues);

        mmgStep.setIgnore(dirichletNodes);

        //////////////////////////////////
        //  Create the transfer operator
        //////////////////////////////////

        std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<VectorType> > > mgTransfers(grid->maxLevel());
        for (size_t i=0; i<mgTransfers.size(); i++) {

            mgTransfers[i] = std::make_shared<TruncatedCompressedMGTransfer<VectorType> >();
            mgTransfers[i]->setup(*grid,i,i+1);
        }

        mmgStep.setTransferOperators(mgTransfers);

        // //////////////////////////////////////////////////////////
        //   Create obstacles
        // //////////////////////////////////////////////////////////

        mmgStep.setHasObstacles(BitSetVector<dim>(dirichletNodes.size(), true));

        LaplaceAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement,MatrixType::block_type> laplace;
        OperatorAssembler<P1Basis, P1Basis> operatorAssembler(p1Basis, p1Basis);

        MatrixType laplaceMatrix;
        operatorAssembler.assemble(laplace, laplaceMatrix);
        h1SemiNorm.setMatrix(&laplaceMatrix);


        //////////////////////////////////////////////////////////////
        //  Create the nonlinear elasticity problem
        /////////////////////////////////////////////////////////////

        // Make static contact problem
        typedef NonlinearElasticityProblem<VectorType,MatrixType, P1Basis> ElasticityProblem;
        ElasticityProblem elasticProblem(material, extForces);

        // Create trust-region step
        TrustRegionSolver<ElasticityProblem, VectorType, MatrixType> trustRegionSolver;
        //TODO Allow a different norm here
        const ParameterTree& trConfig = parameterSet.sub("trConfig");
        trustRegionSolver.setup(trConfig, h1SemiNorm, elasticProblem, solver);


        typedef BasisGridFunction<P1Basis,VectorType> BasisGridFunction;
        VectorType deformedX;

        do {

            deformedX = x;

            if (loadFactor < 1) {

                // ////////////////////////////////////////////////////
                //   Compute new feasible load increment
                // ////////////////////////////////////////////////////


                loadIncrement *= 2;  // double it once, cause we're going to half it at least once, too!
                bool isNan;
                do {

                    loadIncrement /= 2;

                    // Set new Dirichlet values in solution
                    for (size_t k=0; k<dirichletNodes.size(); k++)
                        for (int j=0; j<dim; j++)
                            if (dirichletNodes[k][j])
                                deformedX[k][j] = dirichletValues[k][j] * (loadFactor + loadIncrement);

                    isNan = false;

                    auto disp = std::make_shared<BasisGridFunction>(p1Basis, deformedX);
                    field_type eng = material.energy(disp);
                    std::cout<<"Energy "<<eng<<std::endl;
                    isNan |= std::isnan(eng);
                } while (isNan);

                loadFactor += loadIncrement;

                std::cout << "####################################################" << std::endl;
                std::cout << "New load factor: " << loadFactor
                    << "    new load increment: " << loadIncrement << std::endl;
                std::cout << "####################################################" << std::endl;

            }
            // add Dirichlet values
            trustRegionSolver.setInitialIterate(deformedX);
            trustRegionSolver.solve();
            x = trustRegionSolver.getSol();

            LeafAmiraMeshWriter<GridType> amiramesh2;
            amiramesh2.addLeafGrid(*grid,true);
            amiramesh2.addVertexData(x, grid->leafGridView(),true);
            std::string name = "loadingStep" + std::to_string(loadFactor);
            amiramesh2.write(resultPath  + name, 1);

        } while (loadFactor < 1);


        // /////////////////////////////////////////////////////////////////////
        //   Refinement Loop
        // /////////////////////////////////////////////////////////////////////

        if (level==maxLevel)
            break;

        // ////////////////////////////////////////////////////////////////////////////
        //    Refine locally and transfer the current solution to the new leaf level
        // ////////////////////////////////////////////////////////////////////////////

        std::cout<<"Estimating error on  refinement level "<<level<<std::endl;

        // Create the materials
        typedef P2NodalBasis<GridType::LeafGridView> P2Basis;
        P2Basis p2Basis(grid->leafGridView());
#if HAVE_ADOLC
        typedef MooneyRivlinMaterial<P2Basis> MaterialType2;
        MaterialType2 p2localEnergy(p2Basis, parameterSet.get<field_type>("E"),
                                   parameterSet.get<field_type>("nu"));
        AdolcMaterial<P2Basis> p2Material(p2Basis, p2localEnergy, false);
#else
        using MaterialType2 = NeoHookeanMaterial<P2Basis>;
        MaterialType material(p2Basis,
                          parameterSet.get<double>("E"),
                          parameterSet.get<double>("nu"));
#endif

        // P2 Forces
        VectorType p2ExtForces(p2Basis.size());
        p2ExtForces = 0;

        typedef NonlinearElasticityProblem<VectorType,MatrixType, P2Basis> ElasticityProblemP2;
        ElasticityProblemP2 elasticProblemP2(p2Material, p2ExtForces);


        RefinementIndicator<GridType> refinementIndicator(*grid);

        BitSetVector<1> scalarDirichletField;
        DirichletBCAssembler<GridType>::inflateBitField(dirichletNodes, scalarDirichletField);
        BoundaryPatch<GridType::LeafGridView> dirichletBoundary(grid->leafGridView(),scalarDirichletField);

        HierarchicalEstimator<P2Basis,dim> estimator(*grid);
        estimator.estimate(x, dirichletBoundary, refinementIndicator,elasticProblemP2);


        // ////////////////////////////////////////////////////
        //   Refine grids
        // ////////////////////////////////////////////////////

        std::cout<<"Mark elements"<<std::endl;
        FractionalMarkingStrategy<GridType>::mark(refinementIndicator, *grid,
                                    parameterSet.get<field_type>("refinementFraction"));

        GridFunctionAdaptor<P1Basis> adaptor(p1Basis,true,true);
        grid->preAdapt();
        grid->adapt();
        grid->postAdapt();

        p1Basis.update();
        adaptor.adapt(x);


        std::cout << "########################################################" << std::endl;
        std::cout << "  Grid refined," <<" max level: " << grid->maxLevel()
                << "   vertices: " << grid->size(dim)
                << "   elements: " << grid->size(0)<<std::endl;
        std::cout << "####################################################" << std::endl;

        Dune::LeafAmiraMeshWriter<GridType> amiramesh2;
        amiramesh2.addLeafGrid(*grid,true);
        amiramesh2.addVertexData(x,grid->leafGridView(),true);
        std::ostringstream name;
        name << "nonlinelastRefined" << std::setw(6) << std::setfill('0') <<level ;
        amiramesh2.write(resultPath +name.str(),1);


    }

    LeafAmiraMeshWriter<GridType> amiramesh;
    amiramesh.addLeafGrid(*grid,true);
    amiramesh.addVertexData(x, grid->leafGridView(),true);
    amiramesh.write(resultPath +"nonlineast.result",1);

} catch (Exception e) {

    std::cout << e << std::endl;

}