pntest.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include<iostream>
// dune includes
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/common/timer.hh>
#if HAVE_DUNE_ALUGRID
#include<dune/alugrid/grid.hh>
#include<dune/alugrid/dgf.hh>
#include<dune/grid/io/file/dgfparser/dgfparser.hh>
#endif
// pdelab includes
#include<dune/pdelab/finiteelementmap/pkfem.hh>
#include<dune/pdelab/finiteelementmap/qkfem.hh>
// phasefield includes
#include<dune/phasefield/porenetwork/pn_1pflow.hh>

//===============================================================
// Main program
//===============================================================


int main(int argc, char** argv)
{
  try{
    // Maybe initialize Mpi
    Dune::MPIHelper&
      helper = Dune::MPIHelper::instance(argc, argv);
    if(Dune::MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout << "Parallel code run on "
                << helper.size() << " process(es)" << std::endl;

    // open ini file
    Dune::ParameterTree ptree;
    Dune::ParameterTreeParser ptreeparser;
    ptreeparser.readINITree("pn1p.ini",ptree);
    ptreeparser.readOptions(argc,argv,ptree);

    typedef double RF;
    typedef Params_fs_r< RF,
                          //DoubleWell_limited< RF, DoubleWell_poly<RF>, Limiter_FakeDiverging<RF> >,
                          DoubleWell_limited< RF, DoubleWell_8th_order_poly<RF>, Limiter_FakeDiverging<RF> >,
                          Mobility_Const<RF>,
                          Reaction_Const<RF>,
                          VelDissipation_Quadratic<RF>
                          //VelDissipation_Quadratic_Shifted<RF>
                          > Parameters;
    Parameters param(ptree);
    std::string GridFilename = ptree.get("domain.filename","square.msh");
    std::string OutputFilename = ptree.get("output.filename","output");
    Pn_1PFlow<Parameters> pn_1pFlow(GridFilename,OutputFilename,param);

    std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
    pn_1pFlow.writeVTK(0.0);
    pn_1pFlow.grow(0.2);
    pn_1pFlow.writeVTK(1.0);
    std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;

  }
  catch (Dune::Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
    return 1;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
    return 1;
  }
}