// Copyright (C) 2008-2016 Johan Jansson and Niclas Jansson as main authors.

// Licensed under the GNU LGPL Version 2.1.
//
// This program solves the incompressible Navier-Stokes equations using
// least-squares-Galerkin stabilized FEM with a Schur-preconditioned fixed-point
// iteration between the momentunm and pressure equations and a do-nothing adaptive
// method.

#include "PrimalSolver.h"

using namespace dolfin;

constexpr real rho = 1050;
constexpr real mu = 2.65e-3;
constexpr real nu = mu/rho;

constexpr real Uin = 0.5;

constexpr real tFinal = 2;

constexpr real cfl_target = 20;

int main(int argc, char* argv[])
{
  dolfin_set("output destination","silent");
  if(dolfin::MPI::processNumber() == 0)
      dolfin_set("output destination","terminal");

  if(!ParameterSystem::parameters.defined("alpha"))
      dolfin_add("alpha", 0.0);
  if(!ParameterSystem::parameters.defined("adapt_percent"))
      dolfin_add("adapt_percent", 10.0);
  if(!ParameterSystem::parameters.defined("discrete_tolerance"))
      dolfin_add("discrete_tolerance", 1e-2);

  // Create mesh
  Mesh mesh("Domain.xml");

  dolfin_set("Mesh read in serial", true);

  // Create boundary conditions
  Mesh wallsMesh("Walls.xml");
  Mesh inletMesh("Inlet.xml");
  Mesh outletMesh("Outlet.xml");
  Mesh cylinderMesh("Cylinder.xml");

  dolfin_set("Mesh read in serial", false);

  SubDomain walls(wallsMesh);
  SubDomain inlet(inletMesh);
  SubDomain outlet(outletMesh);
  SubDomain cylinder(cylinderMesh);

  Function u_i0(mesh,Array<real>(0,0,0));
  Function inflow(mesh,Array<real>(Uin,0,0));
  Function noSlipVelocity(mesh,Array<real>(0,0,0));
  Function outflow(mesh,0);
  Function f(mesh,Array<real>(0,0,0));

  DirichletBCList dbcs_m;
  dbcs_m.push_back(std::make_pair(&inlet,&inflow));
  dbcs_m.push_back(std::make_pair(&cylinder,&noSlipVelocity));
  dbcs_m.push_back(std::make_pair(&walls,&noSlipVelocity));

  DirichletBCList dbcs_c;
  dbcs_c.push_back(std::make_pair(&outlet,&outflow));

  DirichletBCList dbcs_dm;

  OutputQuantityList outputQuantities;
  typedef PrimalSolver<NoSlip,
                       GlobalTimeStep,
                       UnitVelocityStabilization> NSESolver;

  NSESolver solver(
      mesh,
      dbcs_m,
      dbcs_c,
      dbcs_dm,
      NULL,
      &u_i0,
      outputQuantities,
      NULL,
      NULL,
      NULL,
      &f,
      cfl_target,
      nu
      );
  solver.setT(tFinal);
  solver.setNoSamples(100);
  solver.run();

  return 0;
}
