function output = generate_rhs (ufl_name)

STRING ="\n\
#include ""@@UFL_NAME@@.h""\n\
#include ""functionspace.h""\n\
#include ""boundarycondition.h""\n\
#include ""coefficient.h""\n\
#include ""function.h""\n\
\n\
DEFUN_DLD (fem_rhs_@@UFL_NAME@@, args, , ""A = fem_rhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF, BC)"")\n\
{\n\
\n\
  int nargin = args.length ();\n\
  octave_value retval;\n\
\n\
  if (nargin < 1)\n\
    print_usage ();\n\
  else\n\
    {\n\
      if (! functionspace_type_loaded)\n\
        {\n\
          functionspace::register_type ();\n\
          functionspace_type_loaded = true;\n\
          mlock ();\n\
        }\n\
      if (args(0).type_id () == functionspace::static_type_id ())\n\
        {\n\
          const functionspace & fspo\n\
            = static_cast<const functionspace&> (args(0).get_rep ());\n\
\n\
          if (! error_state)\n\
            {\n\
              const dolfin::FunctionSpace V = fspo.get_fsp ();\n\
              @@UFL_NAME@@::BilinearForm a (V, V);\n\
              std::size_t ncoef = a.num_coefficients (), nc = 0;\n\
\n\
              if (! coefficient_type_loaded)\n\
                {\n\
                  coefficient::register_type ();\n\
                  coefficient_type_loaded = true;\n\
                  mlock ();\n\
                }\n\
\n\
              if (! function_type_loaded)\n\
                {\n\
                  function::register_type ();\n\
                  function_type_loaded = true;\n\
                  mlock ();\n\
                }\n\
\n\
              for (std::size_t i = 1; i < nargin; ++i)\n\
                {\n\
                  if (args(i).type_id () == coefficient::static_type_id ())\n\
                    {\n\
                      const coefficient & cf\n\
                        = static_cast <const coefficient&> (args(i).get_rep ());\n\
\n\
                      std::size_t n = a.coefficient_number (cf.get_str ());\n\
                      const boost::shared_ptr<const expression> & pexp = cf.get_expr ();\n\
                      a.set_coefficient (n, pexp);\n\
                      ++nc;\n\
                    }\n\
\n\
                  if (args(i).type_id () == function::static_type_id ())\n\
                    {\n\
                      const function & fun\n\
                        = static_cast <const function&> (args(i).get_rep ());\n\
\n\
                      std::size_t n = a.coefficient_number (fun.get_str ());\n\
                      const boost::shared_ptr<const dolfin::Function> & pfun = fun.get_pfun ();\n\
                      a.set_coefficient (n, pfun);\n\
                      ++nc;\n\
                    }\n\
                 }\n\
\n\
              if (nc != ncoef)\n\
                error (""Wrong number of coefficient"");\n\
              else\n\
                {\n\
\n\
                  dolfin::Matrix A;\n\
                  dolfin::assemble (A, a);\n\
\n\
                  if (! boundarycondition_type_loaded)\n\
                    {\n\
                      boundarycondition::register_type ();\n\
                      boundarycondition_type_loaded = true;\n\
                      mlock ();\n\
                    }\n\
\n\
                  for (std::size_t i = 1; i < nargin; ++i)\n\
                    {\n\
                      if (args(i).type_id () == boundarycondition::static_type_id ())\n\
                        {\n\
                          const boundarycondition & bc\n\
                            = static_cast<const boundarycondition&> (args(i).get_rep ());\n\
\n\
                          const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > & pbc\n\
                                = bc.get_bc ();\n\
                              for (std::size_t j = 0; j < pbc.size (); ++j)\n\
                                pbc[j]->apply(A);\n\
                        }\n\
                    }\n\
\n\
                  std::size_t nr = A.size (0), nc = A.size (1);\n\
                  // nz shoud be estimated in a better way\n\
                  double alpha = 0.005;\n\
                  std::size_t nz = alpha * nr * nc;\n\
                  SparseMatrix sm (nr, nc, ceil (nz));\n\
\n\
                  std::size_t ii = 0;\n\
                  sm.cidx (0) = 0;\n\
                  for (std::size_t j = 0; j < nc; ++j)\n\
                   {\n\
                     for (std::size_t i = 0; i < nr; ++i)\n\
                       {\n\
                         double tmp = A(i, j);\n\
                         if (tmp != 0.)\n\
                           {\n\
                             if (ii == nz)\n\
                               {\n\
                                 nz = 1.1 * ((nc * ii) / (j + 1));\n\
                                 sm.change_capacity (ceil (nz));\n\
                               }\n\
                             sm.data(ii) = tmp;\n\
                             sm.ridx(ii) = i;\n\
                             ++ii;\n\
                           }\n\
                       }\n\
                     sm.cidx(j+1) = ii;\n\
                  }\n\
                  sm.maybe_compress ();\n\
\n\
                  retval = sm;\n\
\n\
                }\n\
            }\n\
        }\n\
    }\n\
  return retval;\n\
}";

STRING =  strrep (STRING, "@@UFL_NAME@@", ufl_name);

fid = fopen (sprintf ("fem_rhs_%s.cc", ufl_name), 'w');
if (fid >= 0)
  fputs (fid, STRING);
  output = fclose (fid);
else
  error ("cannot open file");
  output = 1;
endif

endfunction