// Copyright (C) 2008 Marc Spiegelman
// Licensed under the GNU LGPL Version 2.1.
//
// based on Poisson
//
// This program illustrates the use of the DOLFIN for solving a nonlinear PDE
// by solving the linear Compaction pressure 
//
//     - div f^3 grad u(x, y) + c f u = -div f^3 \khat
//
//
// on the unit square with porosity f given by
//
//     f(x, y) = exp(
//
// and boundary conditions given by
//
//     u(x, y)     = t  for x = 1
//     du/dn(x, y) = 0  otherwise
//
// where t is pseudo time.
//
// This is equivalent to solving: 
// F(u) = (grad(v), (1-u^2)*grad(u)) - f(x,y) = 0

#include <dolfin.h>
#include "PressureEquation.h"
  
using namespace dolfin;

// porosity
class Source : public Function, public TimeDependent
{
  public:

    Source(Mesh& mesh) : Function(mesh) {}

    real eval(const real* x) const
    {
      real dx = x[0] - 0.5;
      real dy = x[1] - 0.5;
      return 1.+3.*exp(-(dx*dx + dy*dy)/0.02);
    }
};


// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary) const
  {
    return ( std::abs(x[1] - 1.0) < DOLFIN_EPS || abs(x[1]) < DOLFIN_EPS) && on_boundary;
  }
};



int main(int argc, char* argv[])
{
  dolfin_init(argc, argv);
 
  // Set up problem
  UnitSquare mesh(32, 32);


  // Create source function
  Source f(mesh);

  // Dirichlet boundary conditions
  Function u1(mesh,1.0);
  DirichletBoundary boundary; 
  DirichletBC bc(u1, mesh, boundary);

  // Solution function
  Function u;

  // Create forms and nonlinear PDE
  PressureEquationBilinearForm a(u);
  PressureEquationLinearForm L(f);
  LinearPDE pde(a, L, mesh, bc);

  // Solve
  pde.solve(u);

  // Plot solution
  plot(u);

  // Save function to file
  File f0("pressure.pvd");
  f0 << u;

  //save f to file?
  File f1("porosity.pvd");
  f1 << f;

  return 0;
}