Skip to content
Snippets Groups Projects
Forked from agnumpde / dune-elasticity
405 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
nonlinelast.cc 9.46 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/elasticity/assemblers/ogdenassembler.hh>

#include <dune/istl/io.hh>

#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/fufem/sampleonbitfield.hh>
#include <dune/fufem/readbitfield.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>

#include <dune/solvers/solvers/quadraticipopt.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>

// Choose a solver
#define IPOPT
//#define GAUSS_SEIDEL
//#define MULTIGRID

//#define IPOPT_BASE

// 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
    ParameterTree parameterSet;
    if (argc==2)
        ParameterTreeParser::readINITree(argv[1], parameterSet);
    else
        ParameterTreeParser::readINITree("nonlinelast.parset", parameterSet);

    // read solver settings
    const int numLevels        = parameterSet.get<int>("numLevels");
    const int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
    const int numNewtonSteps   = parameterSet.get<int>("numNewtonSteps");
    const int numIt            = parameterSet.get("numIt", int(0));
    const int nu1              = parameterSet.get("nu1", int(0));
    const int nu2              = parameterSet.get("nu2", int(0));
    const int mu               = parameterSet.get("mu", int(0));
    const int baseIt           = parameterSet.get("baseIt", int(0));
    const double tolerance     = parameterSet.get("tolerance", double(0));
    const double baseTolerance = parameterSet.get("baseTolerance", double(0));
    const double scaling       = parameterSet.get<double>("scaling");

    // read problem settings
    std::string gridFile       = parameterSet.get<std::string>("gridFile");
    std::string dnFile         = parameterSet.get<std::string>("dnFile");
    std::string dvFile         = parameterSet.get<std::string>("dvFile");

    // /////////////////////////////
    //   Generate the grid
    // /////////////////////////////
    typedef UGGrid<dim> GridType;
    typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;

    GridType grid;

    Dune::AmiraMeshReader<GridType>::read(grid, gridFile);
    
    // Read Dirichlet boundary values
    std::vector<LevelBoundaryPatch>  dirichletBoundary(numLevels);
    dirichletBoundary[0].setup(grid.levelView(0));
    readBoundaryPatch<GridType>(dirichletBoundary[0], dnFile);

    std::vector<VectorType> dirichletValues(numLevels);
    dirichletValues[0].resize(grid.size(0, dim));
    AmiraMeshReader<GridType>::readFunction(dirichletValues[0], dvFile);

    dirichletValues[0] *= scaling;

    for (int i=0; i<numLevels-1; i++)
        grid.globalRefine(1);

    BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary);

    std::vector<BitSetVector<dim> >  dirichletNodes(numLevels);
    for (int i=0; i<numLevels; i++) {
        dirichletNodes[i].resize(grid.size(i, dim));
        for (int j=0; j<grid.size(i, dim); j++) {

            for (int k=0; k<dim; k++)
                dirichletNodes[i][j][k] = dirichletBoundary[i].containsVertex(j);
        }

    }

    // hack Dirichlet data
    dirichletNodes[numLevels-1][2*61  ] = false;
    dirichletNodes[numLevels-1][2*61+1] = false;
    dirichletNodes[numLevels-1][2*62  ] = false;
    dirichletNodes[numLevels-1][2*62+1] = false;
    dirichletNodes[numLevels-1][2*56  ] = false;
    dirichletNodes[numLevels-1][2*56+1] = false;
    dirichletNodes[numLevels-1][2*58  ] = false;
    dirichletNodes[numLevels-1][2*58+1] = false;
    dirichletNodes[numLevels-1][2*79  ] = false;
    dirichletNodes[numLevels-1][2*79+1] = false;
    dirichletNodes[numLevels-1][2*74  ] = false;
    dirichletNodes[numLevels-1][2*74+1] = false;
    dirichletNodes[numLevels-1][2*76  ] = false;
    dirichletNodes[numLevels-1][2*76+1] = false;
    dirichletNodes[numLevels-1][2*73  ] = false;
    dirichletNodes[numLevels-1][2*73+1] = false;

    dirichletNodes[numLevels-1][2*4  ] = false;
    dirichletNodes[numLevels-1][2*4+1] = false;
    dirichletNodes[numLevels-1][2*1  ] = false;
    dirichletNodes[numLevels-1][2*1+1] = false;
    dirichletNodes[numLevels-1][2*11  ] = false;
    dirichletNodes[numLevels-1][2*11+1] = false;
    dirichletNodes[numLevels-1][2*9  ] = false;
    dirichletNodes[numLevels-1][2*9+1] = false;
    dirichletNodes[numLevels-1][2*27  ] = false;
    dirichletNodes[numLevels-1][2*27+1] = false;
    dirichletNodes[numLevels-1][2*25  ] = false;
    dirichletNodes[numLevels-1][2*25+1] = false;
    dirichletNodes[numLevels-1][2*33  ] = false;
    dirichletNodes[numLevels-1][2*33+1] = false;
    

    sampleOnBitField(grid, dirichletValues, dirichletNodes);
    dirichletValues[numLevels-1][0][1] = 0.2;
    dirichletValues[numLevels-1][31][1] = 0.2;

    int maxlevel = grid.maxLevel();

  // Determine Dirichlet dofs
    VectorType rhs(grid.size(grid.maxLevel(), dim));
    VectorType x(grid.size(grid.maxLevel(), dim));
    VectorType corr(grid.size(grid.maxLevel(), dim));
    OperatorType problemMatrix;
    
    // Initial solution
    x = 0;
    corr = 0;
    rhs = 0;

    // Create a solver
#if defined IPOPT

    QuadraticIPOptSolver<OperatorType,VectorType> solver;
    solver.verbosity_ = NumProc::QUIET;

    solver.ignoreNodes_ = &dirichletNodes[maxlevel];

#elif defined GAUSS_SEIDEL

    L2Norm<ContactVector<dim> > l2Norm;

    typedef ProjectedBlockGSStep<OperatorType, VectorType> SmootherType;
    SmootherType projectedBlockGSStep(*bilinearForm, x, rhs);
    projectedBlockGSStep.dirichletNodes_ = &dirichletNodes[maxlevel];

    ::LoopSolver<OperatorType, ContactVector<dim> > solver;
    solver.iterationStep = &projectedBlockGSStep;
    solver.numIt = numIt;
    solver.verbosity_ = Solver::FULL;
    solver.errorNorm_ = &l2Norm;
    solver.tolerance_ = tolerance;
    
#elif defined MULTIGRID

    L2Norm<ContactVector<dim> > l2Norm;

    // First create a base solver
#ifdef IPOPT_BASE

    LinearIPOptSolver baseSolver;

#else // Gauss-Seidel is the base solver

    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > baseSolverStep;

    ::LoopSolver<OperatorType, ContactVector<dim> > baseSolver;
    baseSolver.iterationStep = &baseSolverStep;
    baseSolver.numIt = baseIt;
    baseSolver.verbosity_ = Solver::QUIET;
    baseSolver.errorNorm_ = &l2Norm;
    baseSolver.tolerance_ = baseTolerance;
#endif

    // Make pre and postsmoothers
    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > presmoother;
    ProjectedBlockGSStep<OperatorType, ContactVector<dim> > postsmoother;
    

    ContactMMGStep<OperatorType, ContactVector<dim>, FuncSpaceType > contactMMGStep(*bilinearForm, x, rhs, numLevels);

    contactMMGStep.setMGType(1, nu1, nu2);
    contactMMGStep.dirichletNodes_ = &dirichletNodes;
    contactMMGStep.basesolver_     = &baseSolver;
    contactMMGStep.presmoother_    = &presmoother;
    contactMMGStep.postsmoother_   = &postsmoother;    
    contactMMGStep.funcSpace_      = funcSpace;

    ::LoopSolver<OperatorType, ContactVector<dim> > solver;
    solver.iterationStep = &contactMMGStep;
    solver.numIt = numIt;
    solver.verbosity_ = Solver::FULL;
    solver.errorNorm_ = &l2Norm;
    solver.tolerance_ = tolerance;

#else
    #warning You have to specify a solver!
#endif
    typedef P1NodalBasis<GridType::LeafGridView,double> P1Basis;
    typedef OgdenAssembler<P1Basis,P1Basis> Assembler;
    
    P1Basis p1Basis(grid.leafView());
    Assembler ogdenAssembler(p1Basis,p1Basis);
    
    for (int i=0; i<numHomotopySteps; i++) {

        // scale Dirichlet values
        VectorType scaledDirichlet = dirichletValues[maxlevel];
        scaledDirichlet *= (double(i)+1)/numHomotopySteps;

        // Set new Dirichlet values in solution
        for (int j=0; j<x.size(); j++)
            for (int k=0; k<dim; k++)
                if (dirichletNodes[maxlevel][j][k])
                    x[j][k] = scaledDirichlet[j][k];



        for (int j=0; j<numNewtonSteps; j++) {

        	std::cout << "iteration: " << j << std::endl;

        	// Assemble the Jacobi matrix at the current solution
            rhs = 0;
            ogdenAssembler.assembleProblem(problemMatrix,x, rhs);

            corr = 0;
          
            // Set new Dirichlet values in the right hand side
            for (int k=0; k<rhs.size(); k++)
            	for (int l=0; l<dim; l++)
            		if (dirichletNodes[maxlevel][k][l])
            			rhs[k][l] = 0;

            solver.setProblem(problemMatrix, corr, rhs);

            solver.preprocess();
#ifdef MULTIGRID
          contactMMGStep.preprocess();
#endif

          // Solve correction problem
          solver.solve();

          // Add damped correction to the current solution
          x += corr;

          std::cout << "Correction: " << corr.infinity_norm() << std::endl;

      }

  }

#ifdef MULTIGRID
  x = contactMMGStep.getSol();
#endif

  //contactAssembler.postprocess(x);

  //std::cout << x << std::endl;

  // Output result
  LeafAmiraMeshWriter<GridType> amiramesh;
  amiramesh.addLeafGrid(grid,true);
  amiramesh.addVertexData(x, grid.leafView());
  amiramesh.write("resultGrid",true);

 } catch (Exception e) {

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

 }