finite-strain-elasticity.cc

#include <config.h>
#include <functional>

// Includes for the ADOL-C automatic differentiation library
// Need to come before (almost) all others.
#include <adolc/adouble.h>
#include <adolc/drivers/drivers.h>    // use of "Easy to Use" drivers
#include <adolc/taping.h>

#include <dune/fufem/utilities/adolcnamespaceinjections.hh>

#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>

#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>

#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/io/file/vtk.hh>

#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/dunepython.hh>

#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>

#include <dune/elasticity/common/trustregionsolver.hh>
#include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh>
#include <dune/elasticity/materials/localintegralenergy.hh>
#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
#include <dune/elasticity/materials/henckydensity.hh>
#include <dune/elasticity/materials/exphenckydensity.hh>
#include <dune/elasticity/materials/mooneyrivlindensity.hh>
#include <dune/elasticity/materials/neohookedensity.hh>
#include <dune/elasticity/materials/neumannenergy.hh>
#include <dune/elasticity/materials/sumenergy.hh>

// grid dimension
const int dim = 3;

const int order = 1;

using namespace Dune;

int main (int argc, char *argv[]) try
{
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);

  if (mpiHelper.rank()==0)
    std::cout << "ORDER = " << order << std::endl;
  // Start Python interpreter
  Python::start();
  Python::Reference main = Python::import("__main__");
  Python::run("import math");

  //feenableexcept(FE_INVALID);
  Python::runStream()
        << std::endl << "import sys"
        << std::endl << "import os"
        << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
        << std::endl;

  typedef BlockVector<FieldVector<double,dim> > SolutionType;

  // parse data file
  ParameterTree parameterSet;

  ParameterTreeParser::readINITree(argv[1], parameterSet);

  ParameterTreeParser::readOptions(argc, argv, parameterSet);

  // read solver settings
  const int numLevels                   = parameterSet.get<int>("numLevels");
  int numHomotopySteps                  = parameterSet.get<int>("numHomotopySteps");
  const double tolerance                = parameterSet.get<double>("tolerance");
  const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
  const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
  const int multigridIterations         = 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 baseIterations              = parameterSet.get<int>("baseIt");
  const double mgTolerance              = parameterSet.get<double>("mgTolerance");
  const double baseTolerance            = parameterSet.get<double>("baseTolerance");
  const double damping                  = parameterSet.get<double>("damping");
  std::string resultPath                = parameterSet.get("resultPath", "");

  // ///////////////////////////////////////
  //    Create the grid
  // ///////////////////////////////////////
  typedef UGGrid<dim> GridType;

  std::shared_ptr<GridType> grid;

  FieldVector<double,dim> lower(0), upper(1);

  if (parameterSet.get<bool>("structuredGrid")) {

    lower = parameterSet.get<FieldVector<double,dim> >("lower");
    upper = parameterSet.get<FieldVector<double,dim> >("upper");

    std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);

  } else {
    std::string path                = parameterSet.get<std::string>("path");
    std::string gridFile            = parameterSet.get<std::string>("gridFile");
    grid = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
  }

  grid->globalRefine(numLevels-1);

  grid->loadBalance();

  if (mpiHelper.rank()==0)
    std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
  
#if HAVE_DUNE_PARMG
  using GridView = GridType::LevelGridView;
  GridView gridView = grid->levelGridView(grid->maxLevel());
#else
  using GridView = GridType::LeafGridView;
  GridView gridView = grid->leafGridView();
#endif


  ////////////////////////////////////////////
  // Basis
  ////////////////////////////////////////////

  using namespace Dune::Functions::BasisFactory;
  auto basis = makeBasis(
    gridView,
    power<dim>(
      lagrange<order>(),
      blockedInterleaved()
  ));

  using PowerBasis = decltype(basis);
  using LocalView = typename PowerBasis::LocalView;

  // /////////////////////////////////////////
  //   Read Dirichlet values
  // /////////////////////////////////////////

  BitSetVector<dim> dirichletVertices(gridView.size(dim), false);
  BitSetVector<dim> neumannVertices(gridView.size(dim), false);

  const GridView::IndexSet& indexSet = gridView.indexSet();

  // Make Python function that computes which vertices are on the Dirichlet boundary,
  // based on the vertex positions.
  std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));

  // Same for the Neumann boundary
  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));

  for (auto&& v : vertices(gridView))
  {
    bool isDirichlet;
    pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
    dirichletVertices[indexSet.index(v)] = isDirichlet;

    bool isNeumann;
    pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
    neumannVertices[indexSet.index(v)] = isNeumann;
  }

  BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
  auto neumannBoundary = std::make_shared<BoundaryPatch<GridView>>(gridView, neumannVertices);

  std::cout << "On process " << mpiHelper.rank() << ": Neumann boundary has " << neumannBoundary->numFaces() << " faces\n";


  BitSetVector<dim> dirichletNodes(basis.size(), false);
  constructBoundaryDofs(dirichletBoundary,basis,dirichletNodes);

  BitSetVector<dim> neumannNodes(basis.size(), false);
  constructBoundaryDofs(*neumannBoundary,basis,neumannNodes);

  BitSetVector<dim> dirichletDofs(basis.size(), false);
  for (size_t i=0; i<basis.size(); i++)
    if (dirichletNodes[i][0])
      for (int j=0; j<dim; j++)
        dirichletDofs[i][j] = true;

  // //////////////////////////
  //   Initial iterate
  // //////////////////////////

  SolutionType x(basis.size());

  lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
  PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));

  Dune::Functions::interpolate(basis, x, pythonInitialDeformation);

  ////////////////////////////////////////////////////////
  //   Main homotopy loop
  ////////////////////////////////////////////////////////

  // Output initial iterate (of homotopy loop)
  SolutionType identity;
  Dune::Functions::interpolate(basis, identity, [](FieldVector<double,dim> x){ return x; });

  SolutionType displacement = x;
  displacement -= identity;

  auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis, displacement);
  auto localDisplacementFunction = localFunction(displacementFunction);

  //  We need to subsample, because VTK cannot natively display real second-order functions
  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
  vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
  vtkWriter.write(resultPath + "finite-strain_homotopy_0");

  Dune::Timer homotopyTimer;
  for (int i=0; i<numHomotopySteps; i++)
  {
    double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
    if (mpiHelper.rank()==0)
      std::cout << "Homotopy step: " << i << ",    parameter: " << homotopyParameter << std::endl;


    // ////////////////////////////////////////////////////////////
    //   Create an assembler for the energy functional
    // ////////////////////////////////////////////////////////////

    const ParameterTree& materialParameters = parameterSet.sub("materialParameters");

    FieldVector<double,dim> neumannValues(0.);
    if (parameterSet.hasKey("neumannValues"))
    {
      neumannValues = parameterSet.get<FieldVector<double,dim> >("neumannValues");
    }

    // create simple constant valued Neumann function
    auto neumannFunction = [&]( FieldVector<double,dim> )
    {
      auto nv = neumannValues;
      nv *= (-homotopyParameter);
      return nv;
    };

    if (mpiHelper.rank()==0)
    {
      std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,dim> >("neumannValues") << std::endl;
    }

    if (mpiHelper.rank() == 0)
    {
      std::cout << "Material parameters:" << std::endl;
      materialParameters.report();
    }

    // Assembler using ADOL-C
    if (mpiHelper.rank()==0)
    std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;

    std::shared_ptr<Elasticity::LocalDensity<dim,adouble>> elasticDensity;

    if (parameterSet.get<std::string>("energy") == "stvenantkirchhoff")
      elasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,adouble>>(materialParameters);
    if (parameterSet.get<std::string>("energy") == "neohooke")
      elasticDensity = std::make_shared<Elasticity::NeoHookeDensity<dim,adouble>>(materialParameters);
    if (parameterSet.get<std::string>("energy") == "hencky")
      elasticDensity = std::make_shared<Elasticity::HenckyDensity<dim,adouble>>(materialParameters);
    if (parameterSet.get<std::string>("energy") == "exphencky")
      elasticDensity = std::make_shared<Elasticity::ExpHenckyDensity<dim,adouble>>(materialParameters);
    if (parameterSet.get<std::string>("energy") == "mooneyrivlin")
      elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,adouble>>(materialParameters);

    if(!elasticDensity)
      DUNE_THROW(Exception, "Error: Selected energy not available!");

    auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);

    auto neumannEnergy = std::make_shared<Elasticity::NeumannEnergy<LocalView,adouble>>(neumannBoundary,neumannFunction);

    auto totalEnergy = std::make_shared<Elasticity::SumEnergy<LocalView,adouble>>(elasticEnergy, neumannEnergy);

    auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(totalEnergy);

    Elasticity::FEAssembler<PowerBasis,SolutionType> assembler(basis, localADOLCStiffness);

    // /////////////////////////////////////////////////
    //   Create a Riemannian trust-region solver
    // /////////////////////////////////////////////////

    TrustRegionSolver<PowerBasis,SolutionType> solver;
    solver.setup(*grid,
                 &assembler,
                 x,
                 dirichletDofs,
                 tolerance,
                 maxTrustRegionSteps,
                 initialTrustRegionRadius,
                 multigridIterations,
                 mgTolerance,
                 mu, nu1, nu2,
                 baseIterations,
                 baseTolerance,
                 damping
                );

    ////////////////////////////////////////////////////////
    //   Set Dirichlet values
    ////////////////////////////////////////////////////////

    // The python class that implements the Dirichlet values
    Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");

    // The member method 'DirichletValues' of that class
    Python::Callable C = dirichletValuesClass.get("DirichletValues");

    // Call a constructor.
    Python::Reference dirichletValuesPythonObject = C(homotopyParameter);

    // Extract object member functions as Dune functions
    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> >   dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));

    Dune::Functions::interpolate(basis, x, dirichletValues, dirichletDofs);

    // /////////////////////////////////////////////////////
    //   Solve!
    // /////////////////////////////////////////////////////

    solver.setInitialIterate(x);
    solver.solve();

    x = solver.getSol();

    /////////////////////////////////
    //   Output result
    /////////////////////////////////

    // Compute the displacement
    auto displacement = x;
    displacement -= identity;

    auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis, displacement);

    //  We need to subsample, because VTK cannot natively display real second-order functions
    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(2));
    vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
    vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1));
  }

  if (mpiHelper.rank()==0)
    std::cout << "Complete duration: " << homotopyTimer.elapsed() << " sec." << std::endl;

} catch (Exception& e) {
    std::cout << e.what() << std::endl;
}