#ifndef _RHEOLEF_PROBLEM_H
#define _RHEOLEF_PROBLEM_H
//
// This file is part of Rheolef.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
// =========================================================================
// AUTHORS: Pierre.Saramito@imag.fr
// DATE:   21 february 2020

namespace rheolef {
/**
@classfile problem linear solver

Description
===========
The `problem` class solves a given linear system
for PDEs in variational formulation.

Description
===========
The degrees of freedom are splitting between *unknown*
degrees of freedom and *blocked* one.
See also @ref form_2 and @ref space_2.
Let `a` be a bilinear @ref form_2 and `lh` be the right-hand-side,
as in the previous example.
The linear system expands as:

        [ a.uu  a.ub ] [ uh.u ]   [ lh.u ]
        [            ] [      ] = [      ]
        [ a.bu  a.bb ] [ uh.b ]   [ lh.b ]

The `uh.b` are blocked degrees of freedom:
their values are prescribed and the corresponding values
are move to the right-hand-side of the system that reduces to:

        a.uu*uh.u =  lh.u - a.ub*uh.b

This writes:

	problem p (a);
        p.solve (lh, uh);

Observe that, during the `p.solve` call,
`uh` is both an input variable, for the `uh.b`
contribution to the right-hand-side,
and an output variable, with `uh.u`.
When using an iterative resolution,
the details about its convergence, e.g. the number
of iterations and the final residue,
can be obtain via the `p.option()` member function,
see @ref solver_option_4.
Finally, the previous linear system is solved via the @ref solver_4 class:
the `problem` class is simply a convenient wrapper around the @ref solver_4 one.

Example
=======
See @ref dirichlet.cc

Customization
=============
The @ref solver_4 could be customized via the constructor optional
@ref solver_option_4 argument:

        problem p (a, sopt);

When using a direct @ref solver_4, the determinant of the 
linear system matrix is also available as `p.det()`.
When using an iterative @ref solver_4, the preconditionner
could be customized:

	p.set_preconditionner (m);

TODO
====
The `solve` method could return a boolean when success.

Implementation
==============
@showfromfile
The `problem` class is simply an alias to the `problem_basic` class

@snippet problem.h verbatim_problem
@par

The `problem_basic` class provides a generic interface:

@snippet problem.h verbatim_problem_basic
@snippet problem.h verbatim_problem_basic_cont
*/
} // namespace rheolef

#include "rheolef/form.h"

namespace rheolef {

// [verbatim_problem_basic]
template <class T, class M = rheo_default_memory_model>
class problem_basic {
public:

// typedefs:
 
  typedef typename solver_basic<T,M>::size_type         size_type;
  typedef typename solver_basic<T,M>::determinant_type  determinant_type;

// allocators:

  problem_basic ();
  problem_basic (const form_basic<T,M>& a,
                 const solver_option& sopt = solver_option());

  void update_value (const form_basic<T,M>& a);

  void set_preconditioner (const solver_basic<T,M>&);

// accessor:

  void solve       (const field_basic<T,M>& lh, field_basic<T,M>& uh) const;
  void trans_solve (const field_basic<T,M>& lh, field_basic<T,M>& uh) const;

  determinant_type det() const;
  const solver_option& option() const;
  bool initialized() const;
// [verbatim_problem_basic]

// internals:

  // used by problem_mixed:
  const solver_basic<T,M>& get_solver() const { return _sa; }

// data:
protected:
  form_basic<T,M>   _a;
  solver_basic<T,M> _sa;

// [verbatim_problem_basic_cont]
};
// [verbatim_problem_basic_cont]

// [verbatim_problem]
typedef problem_basic<Float> problem;
// [verbatim_problem]

// -----------------------------------------------------------------------------
// inlined
// -----------------------------------------------------------------------------
template<class T, class M>
inline
problem_basic<T,M>::problem_basic()
 : _a(),  
   _sa() 
{
}
template<class T, class M>
inline
problem_basic<T,M>::problem_basic (
  const form_basic<T,M>& a,
  const solver_option&   sopt)
 : _a(a),  
   _sa (a.uu(), sopt) 
{
}
template<class T, class M>
inline
void
problem_basic<T,M>::update_value (const form_basic<T,M>& a)
{
  _a  = a;
  _sa.update_value (a.uu());
}
template<class T, class M>
inline
void
problem_basic<T,M>::set_preconditioner (const solver_basic<T,M>& m)
{
  _sa.set_preconditioner (m);
}
template<class T, class M>
inline
void
problem_basic<T,M>::solve (const field_basic<T,M>& lh, field_basic<T,M>& uh) const
{
  uh.set_u() = _sa.solve (lh.u() - _a.ub()*uh.b());
}
template<class T, class M>
void
problem_basic<T,M>::trans_solve (const field_basic<T,M>& lh, field_basic<T,M>& uh) const
{
  uh.set_u() = _sa.trans_solve (lh.u() - _a.bu().trans_mult(uh.b()));
}
template<class T, class M>
inline
bool
problem_basic<T,M>::initialized() const
{
  return _sa.initialized();
}
template<class T, class M>
inline
typename problem_basic<T,M>::determinant_type
problem_basic<T,M>::det() const
{
  return _sa.det();
}
template<class T, class M>
inline
const solver_option&
problem_basic<T,M>::option() const
{
  return _sa.option();
}

} // namespace rheolef
#endif // _RHEOLEF_PROBLEM_H