/*===========================================================================

 Copyright (C) 2006-2020 Yves Renard, Julien Pommier.

 This file is a part of GetFEM

 GetFEM  is  free software;  you  can  redistribute  it  and/or modify it
 under  the  terms  of the  GNU  Lesser General Public License as published
 by  the  Free Software Foundation;  either version 3 of the License,  or
 (at your option) any later version along with the GCC Runtime Library
 Exception either version 3.1 or (at your option) any later version.
 This program  is  distributed  in  the  hope  that it will be useful,  but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License and GCC Runtime Library Exception for more details.
 You  should  have received a copy of the GNU Lesser General Public License
 along  with  this program;  if not, write to the Free Software Foundation,
 Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.

===========================================================================*/

/**
   @file bilaplacian.cc
   @brief Bilaplacian problem. A dummy
   bilaplacian problem is solved on a regular mesh, and is compared to
   the analytical solution.

   This program is used to check that getfem++ is working. This is also 
   a good example of use of GetFEM.

   @see laplacian.cc
*/

#include "getfem/getfem_config.h"
#include "getfem/getfem_assembling.h"
#include "getfem/getfem_export.h"
#include "getfem/getfem_regular_meshes.h"
#include "getfem/getfem_fourth_order.h"
#include "getfem/getfem_model_solvers.h"
#include "gmm/gmm.h"
#include "getfem/getfem_derivatives.h"
#include "gmm/gmm_inoutput.h"

using std::endl; using std::cout; using std::cerr;
using std::ends; using std::cin;

/* some GetFEM types that we will be using */
using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
using bgeot::base_node;  /* geometrical nodes(derived from base_small_vector)*/
using bgeot::scalar_type; /* = double */
using bgeot::size_type;   /* = unsigned long */
using bgeot::base_matrix; /* small dense matrix. */

typedef getfem::model_real_plain_vector  plain_vector;

/**************************************************************************/
/*  Exact solution.                                                       */
/**************************************************************************/

scalar_type FT = 0.0;
scalar_type D  = 1.0;
scalar_type nu = 0.0;
scalar_type pressure  = 1.0 ;

#if 1

scalar_type sol_u(const base_node &x)
{ return sin(FT*std::accumulate(x.begin(), x.end(), 0.0)); }
scalar_type sol_lapl_u(const base_node &x)
{ return -FT*FT*sol_u(x) * x.size(); }
scalar_type sol_f(const base_node &x)
{ return FT*FT*FT*FT*sol_u(x)*gmm::sqr(x.size()); }
base_small_vector sol_du(const base_node &x) {
  base_small_vector res(x.size());
  std::fill(res.begin(), res.end(),
	    FT * cos(FT*std::accumulate(x.begin(), x.end(), 0.0)));
  return res;
}
base_small_vector neumann_val(const base_node &x)
{ return -FT*FT*sol_du(x) * scalar_type(x.size()); }

base_matrix sol_hessian(const base_node &x) {
  base_matrix m(x.size(), x.size());
  std::fill(m.begin(), m.end(), -FT*FT*sol_u(x));
  return m;
}

#else

scalar_type sol_u(const base_node &x) {
 return pressure/(2. * D) * x[0] * x[0] * ( x[0] * x[0] / 12. - x[0] / 3. + 1./2. ); }
 
scalar_type sol_lapl_u(const base_node &) { 
return pressure/(2. * D) * ( x[0] * x[0] - 2. * x[0] + 1. ); }

scalar_type sol_f(const base_node &) { return pressure ; }
base_small_vector sol_du(const base_node &x) {
  base_small_vector res(x.size()); 
  res[0] = pressure/(2. * D) * x[0] * ( x[0] * x[0] / 3. - x[0] + 1. ); 
  res[1] = 0.; 
  return res;
}
base_small_vector neumann_val(const base_node &x)
{ cout << "erreur : KL = 1 doit etre selectionne \n";
  return 0 ; /*base_small_vector(x.size()); */ }

base_matrix sol_hessian(const base_node &x)
{ base_matrix m(x.size(), x.size()); 
  m(0,0) = sol_lapl_u(x);
  return m; }


// scalar_type sol_u(const base_node &x) { return x[0]*x[1]; }
// scalar_type sol_lapl_u(const base_node &) { return 0.0; }
// scalar_type sol_f(const base_node &) { return 0.0; }
// base_small_vector sol_du(const base_node &x) {
//   base_small_vector res(x.size()); res[0] = x[1]; res[1] = x[0]; 
//   return res;
// }
// base_small_vector neumann_val(const base_node &x)
// { return base_small_vector(x.size());  }
// 
// base_matrix sol_hessian(const base_node &x)
// { base_matrix m(x.size(), x.size()); m(1,0) = m(0,1) = 1.0; return m; }

#endif

base_matrix sol_mtensor(const base_node &x) { 
  base_matrix m = sol_hessian(x), mm(x.size(), x.size());
  scalar_type l = sol_lapl_u(x);
  for (size_type i = 0; i < x.size(); ++i) mm(i,i) = l * nu;
  gmm::scale(m, (1-nu));
  gmm::add(mm, m);
  gmm::scale(m, -D);
  return m;
}

base_small_vector sol_bf(const base_node &x)
{ return -D * neumann_val(x); }

/**************************************************************************/
/*  Structure for the bilaplacian problem.                                */
/**************************************************************************/

struct bilaplacian_problem {

  enum { CLAMPED_BOUNDARY_NUM = 0, SIMPLE_SUPPORT_BOUNDARY_NUM = 1,
	 FORCE_BOUNDARY_NUM = 2, MOMENTUM_BOUNDARY_NUM = 3};
  /* Boundary conditions :
    CLAMPED_BOUNDARY_NUM        --> normal derivative of the solution
    SIMPLE_SUPPORT_BOUNDARY_NUM --> value of the solution
    FORCE_BOUNDARY_NUM          --> the "big" condition
    MOMENTUM_BOUNDARY_NUM       --> flexure moment condition  */

  getfem::mesh mesh;        /* the mesh */
  getfem::mesh_im mim;      /* the integration methods.                     */
  getfem::mesh_fem mf_u;    /* main mesh_fem, for the bilaplacian solution  */
  getfem::mesh_fem mf_mult; /* mesh_fem for the Dirichlet condition.        */
  getfem::mesh_fem mf_rhs;  /* mesh_fem for the right hand side (f(x),..)   */

  scalar_type residual;     /* max residual for the iterative solvers       */
  size_type dirichlet_version;

  std::string datafilename;
  bgeot::md_param PARAM;

  bool KL;
  size_type NX ;
  size_type boundary_ref ;
  /* boundary_ref = 0 corresponds to clamped edge on the 4 edges   
     boundary_ref = 1 corresponds to :    
                                    
                                               free edge
                                              __________
                                             |          |
                                   simple    |          |   simple
                                   support   |          |   support
                                             |          |
                                             |__________|
                                    
                                               clamped
                                               edge
  */
                                 
  bool solve(plain_vector &U);
  void init(void);
  void compute_error(plain_vector &U);
  bilaplacian_problem(void) : mim(mesh),mf_u(mesh), mf_mult(mesh),
			      mf_rhs(mesh) {}
};

/* Read parameters from the .param file, build the mesh, set finite element
 * and integration methods and selects the boundaries.
 */
void bilaplacian_problem::init(void) {
  std::string MESH_FILE = PARAM.string_value("MESH_FILE");
  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
  std::string FEM_TYPE  = PARAM.string_value("FEM_TYPE","FEM name");
  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
					       "Name of integration method");
  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
  cout << "FEM_TYPE="  << FEM_TYPE << "\n";
  cout << "INTEGRATION=" << INTEGRATION << "\n";

  size_type N;
  if (!MESH_FILE.empty()) {
    mesh.read_from_file(MESH_FILE);
    MESH_TYPE = bgeot::name_of_geometric_trans
      (mesh.trans_of_convex(mesh.convex_index().first_true()));
    cout << "MESH_TYPE=" << MESH_TYPE << "\n";
    N = mesh.dim();
  } else {
    cout << "MESH_TYPE=" << MESH_TYPE << "\n";
    
    /* First step : build the mesh */
    bgeot::pgeometric_trans pgt = 
      bgeot::geometric_trans_descriptor(MESH_TYPE);
    N = pgt->dim();
    std::vector<size_type> nsubdiv(N);
    NX = PARAM.int_value("NX", "Nomber of space steps ") ;
    std::fill(nsubdiv.begin(),nsubdiv.end(), NX );
    getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
			      PARAM.int_value("MESH_NOISED") != 0);

    bgeot::base_matrix M(N,N);
    for (size_type i=0; i < N; ++i) {
      static const char *t[] = {"LX","LY","LZ"};
      M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
    }
    /* scale the unit mesh to [LX,LY,..] and incline it */
    mesh.transformation(M);
  }

  int dv = int(PARAM.int_value("DIRICHLET_VERSION", "Dirichlet version"));
  dirichlet_version = size_type(dv);
  datafilename=PARAM.string_value("ROOTFILENAME","Base name of data files.");
  residual=PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;
  FT = PARAM.real_value("FT"); if (FT == 0.0) FT = 1.0;
  pressure = PARAM.real_value("pressure") ; 
  KL = (PARAM.int_value("KL", "Kirchhoff-Love model or not") != 0);
  D = PARAM.real_value("D", "Flexion modulus");
  if (KL) nu = PARAM.real_value("NU", "Poisson ratio");
  boundary_ref = PARAM.int_value("BOUNDARY_REF");
  
  /* set the finite element on the mf_u */
  getfem::pfem pf_u = getfem::fem_descriptor(FEM_TYPE);
  getfem::pintegration_method ppi = 
    getfem::int_method_descriptor(INTEGRATION);

  mim.set_integration_method(mesh.convex_index(), ppi);
  mf_u.set_finite_element(mesh.convex_index(), pf_u);

  std::string dirichlet_fem_name = PARAM.string_value("DIRICHLET_FEM_TYPE");
  if (dirichlet_fem_name.size() == 0) {
    mf_mult.set_finite_element(mesh.convex_index(), pf_u);
    cout << "No DIRICHLET_FEM_TYPE\n";
  } else {
    cout << "DIRICHLET_FEM_TYPE="  << dirichlet_fem_name << "\n";
    mf_mult.set_finite_element(mesh.convex_index(), 
			       getfem::fem_descriptor(dirichlet_fem_name));
  }

  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
     not used in the .param file */
  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
  if (data_fem_name.size() == 0) {
    if (!pf_u->is_lagrange()) {
      GMM_ASSERT1(false, "You are using a non-lagrange FEM. "
		  << "In that case you need to set "
		  << "DATA_FEM_TYPE in the .param file");
    }
    mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
  } else {
    mf_rhs.set_finite_element(mesh.convex_index(), 
			      getfem::fem_descriptor(data_fem_name));
  }
  
  /* set boundary conditions
   * (Neuman on the upper face, Dirichlet elsewhere) */
  cout << "Selecting Neumann and Dirichlet boundaries\n";
  getfem::mesh_region border_faces;
  getfem::outer_faces_of_mesh(mesh, border_faces);
  
  if (boundary_ref == 0) { // clamped plate
     for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
        mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f());
	mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f());
     }
  }
  
  if (boundary_ref == 1) { // mixed boundary conditions
     for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
       base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());
       un /= gmm::vect_norm2(un);
       if (gmm::abs(un[N-1] - 1.0) <= 0.5) {
         mesh.region(FORCE_BOUNDARY_NUM).add(i.cv(), i.f());
         mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f());
       }
       else {
         mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f());
         if (gmm::abs(un[N-1] + 1.0) <= 0.5)
   	    mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f());
         else
   	    mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f());
       }
     }
   }
   
   if (boundary_ref == 2) { // clamped plate on the left edge, free boundary elsewhere
     for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {
        base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());
        un /= gmm::vect_norm2(un);
	if (un[0] < - 0.5 ) { 
           mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f());
 	   mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f());
	} else { 
	  mesh.region(FORCE_BOUNDARY_NUM).add(i.cv(), i.f());
	  mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f());
	}
     }
  }
}

/* compute the relative error with respect to the exact solution */
void bilaplacian_problem::compute_error(plain_vector &U) {

  GMM_ASSERT1(!mf_rhs.is_reduced(), "To be adapted");
  std::vector<scalar_type> V(mf_rhs.nb_basic_dof());
  getfem::interpolation(mf_u, mf_rhs, U, V);
  for (size_type i = 0; i < mf_rhs.nb_basic_dof(); ++i)
    V[i] -= sol_u(mf_rhs.point_of_basic_dof(i));
  cout.precision(16);
  cout  << "L2 error = " << getfem::asm_L2_norm(mim, mf_rhs, V)  << endl
        << "H1 error = " << getfem::asm_H1_norm(mim, mf_rhs, V)  << endl
        << "H2 error = " << getfem::asm_H2_norm(mim, mf_rhs, V)  << endl
        /*<< "Linfty error = " << gmm::vect_norminf(V)  << endl*/; 
  cout  << "semi-norme H1 = " << getfem::asm_H1_semi_norm(mim, mf_rhs, V)  << endl 
        << "semi-norme H2 = " << getfem::asm_H2_semi_norm(mim, mf_rhs, V)  << endl ;
       
}



/**************************************************************************/
/*  Model.                                                                */
/**************************************************************************/


bool bilaplacian_problem::solve(plain_vector &U) {
  size_type nb_dof_rhs = mf_rhs.nb_dof();
  size_type N = mesh.dim();
  
  cout << "Number of dof for u: " << mf_u.nb_dof() << endl;
  
  getfem::model model;
  
  // Main unknown of the problem.
  model.add_fem_variable("u", mf_u);
  
  // Bilaplacian brick.
  model.add_initialized_scalar_data("D", D);
  if (KL) {
    model.add_initialized_scalar_data("nu", nu);
    getfem::add_bilaplacian_brick_KL(model, mim, "u", "D", "nu");
  } else {
    getfem::add_bilaplacian_brick(model, mim, "u", "D");
  }
  
  // Volumic source term brick.
  plain_vector F(nb_dof_rhs);
  getfem::interpolation_function(mf_rhs, F, sol_f);
  model.add_initialized_fem_data("VolumicData", mf_rhs, F);
  getfem::add_source_term_brick(model, mim, "u", "VolumicData");
  
  
  // Defining the prescribed momentum.
  if (KL) {
    gmm::resize(F, nb_dof_rhs*N*N);
    getfem::interpolation_function(mf_rhs, F,
                                   sol_mtensor, MOMENTUM_BOUNDARY_NUM);
    gmm::scale(F, -1.0);
  } else {
    gmm::resize(F, nb_dof_rhs);
    getfem::interpolation_function(mf_rhs, F, sol_lapl_u,
                                   MOMENTUM_BOUNDARY_NUM);
  }
  
  // Prescribed momentum on the boundary
  model.add_initialized_fem_data("Momentum", mf_rhs, F);
  getfem::add_normal_derivative_source_term_brick
    (model, mim, "u", "Momentum", MOMENTUM_BOUNDARY_NUM);
  
  // Defining the Neumann condition right hand side.
  plain_vector H(nb_dof_rhs*N*N);
  gmm::resize(F, nb_dof_rhs*N);
  if (KL) {
    getfem::interpolation_function(mf_rhs, F,sol_bf,FORCE_BOUNDARY_NUM);
    getfem::interpolation_function(mf_rhs, H,sol_mtensor,FORCE_BOUNDARY_NUM);
    model.add_initialized_fem_data("M", mf_rhs, F);
    model.add_initialized_fem_data("H", mf_rhs, H);
    add_Kirchhoff_Love_Neumann_term_brick
      (model, mim, "u", "H", "M", FORCE_BOUNDARY_NUM);
  }
  else {
    getfem::interpolation_function(mf_rhs, F,neumann_val,FORCE_BOUNDARY_NUM);
    gmm::scale(F, -1.0);
    model.add_initialized_fem_data("Neumanndata", mf_rhs, F);
    add_normal_source_term_brick
      (model, mim, "u", "Neumanndata", FORCE_BOUNDARY_NUM);
  }
  
  // Normal derivative Dirichlet condition brick
  gmm::resize(F, nb_dof_rhs*N);
  gmm::clear(F);
  getfem::interpolation_function(mf_rhs, F, sol_du, CLAMPED_BOUNDARY_NUM);
  model.add_initialized_fem_data("DerDirdata", mf_rhs, F);
  if (dirichlet_version == 0)
    add_normal_derivative_Dirichlet_condition_with_multipliers
      (model, mim, "u", mf_u, CLAMPED_BOUNDARY_NUM, "DerDirdata", false);
  else
    add_normal_derivative_Dirichlet_condition_with_penalization
      (model, mim, "u", PARAM.real_value("EPS_DIRICHLET_PENAL"),
       CLAMPED_BOUNDARY_NUM, "DerDirdata", false);
  
  
  // Dirichlet condition brick.
  gmm::resize(F, nb_dof_rhs);
  getfem::interpolation_function
    (mf_rhs, F, sol_u, SIMPLE_SUPPORT_BOUNDARY_NUM);
  model.add_initialized_fem_data("Dirichletdata", mf_rhs, F);
  if (dirichlet_version == 0)
    add_Dirichlet_condition_with_multipliers
      (model, mim, "u", mf_u, SIMPLE_SUPPORT_BOUNDARY_NUM, "Dirichletdata");
  else
    add_Dirichlet_condition_with_penalization
      (model, mim, "u", PARAM.real_value("EPS_DIRICHLET_PENAL"),
       SIMPLE_SUPPORT_BOUNDARY_NUM, "Dirichletdata", &mf_mult);
  
  // Generic solve.
  cout << "Total number of variables : " << model.nb_dof() << endl;
  gmm::iteration iter(residual, 1, 40000);
  getfem::standard_solve(model, iter);
  gmm::resize(U, mf_u.nb_dof());
  gmm::copy(model.real_variable("u"), U);
  return (iter.converged());
 
}
  
/**************************************************************************/
/*  main program.                                                         */
/**************************************************************************/

int main(int argc, char *argv[]) {

  GETFEM_MPI_INIT(argc, argv); // For parallelized version
  GMM_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
  FE_ENABLE_EXCEPT;        // Enable floating point exception for Nan.

  try {
    bilaplacian_problem p;
    p.PARAM.read_command_line(argc, argv);
    p.init();
    plain_vector U(p.mf_u.nb_dof()), V(p.mf_u.nb_dof()) ;
    if (!p.solve(U)) GMM_ASSERT1(false, "Solve has failed");

    p.compute_error(U);

    if (p.PARAM.int_value("VTK_EXPORT") ) {
      cout << "export to " << p.datafilename + ".vtk" << "..\n";
      getfem::vtk_export exp(p.datafilename + ".vtk",
			     p.PARAM.int_value("VTK_EXPORT")==1);
      exp.exporting(p.mf_u); 
      exp.write_point_data(p.mf_u, U, "bilaplacian_displacement");
      cout << "export done, you can view the data file with (for example)\n"
	   << "mayavi2 -d " << p.datafilename
	   << ".vtk -f WarpScalar -m Surface -m Outline\n";
    }
  }
  GMM_STANDARD_CATCH_ERROR;

  GETFEM_MPI_FINALIZE;

  return 0; 
}