nonlinelast.cc

#include <config.h>

#include <dune/common/bitfield.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 "src/ogdenassembler.hh"

#include <dune/istl/io.hh>

#include <dune/ag-common/prolongboundarypatch.hh>
#include <dune/ag-common/sampleonbitfield.hh>
#include <dune/ag-common/readbitfield.hh>
#include <dune/ag-common/quadraticipopt.hh>
#include <dune/ag-common/iterativesolver.hh>
#include <dune/ag-common/energynorm.hh>
#include <dune/ag-common/mmgstep.hh>

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

//#define IPOPT_BASE

// The grid dimension
const int dim = 2;

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("nonlinelast.parset");

    // 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;
    GridType grid;

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

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

    dirichletValues[0] *= scaling;

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

    PatchProlongator<GridType>::prolong(dirichletBoundary);

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

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

    }

    sampleOnBitField(grid, dirichletValues, dirichletNodes);

    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));
    
    // Initial solution
    x = 0;
    corr = 0;
    rhs = 0;

    // Create a solver
#if defined IPOPT

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

    solver.dirichletNodes_ = &dirichletNodes[maxlevel];

#elif defined GAUSS_SEIDEL

    L2Norm<ContactVector<dim> > l2Norm;

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

    IterativeSolver<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;

    IterativeSolver<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;

    IterativeSolver<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 OgdenAssembler<GridType> Assembler;
    Assembler ogdenAssembler(grid);
    
    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*dim+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.assembleMatrix(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*dim+l])
                      rhs[k][l] = 0;

          solver.setProblem(*ogdenAssembler.matrix_, 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.leafIndexSet());
  amiramesh.write("resultGrid");

 } catch (Exception e) {

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

 }