--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/generate_rhs.m	Sat Jul 27 17:21:03 2013 +0200
@@ -0,0 +1,155 @@
+
+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\
+                  octave_idx_type nr = A.size (0), nc = A.size (1);\n\
+                  // nz shoud be estimated in a better way\n\
+                  double alpha = 0.005;\n\
+                  octave_idx_type nz = ceil (alpha * nr * nc);\n\
+                  SparseMatrix sm (nr, nc, nz);\n\
+\n\
+                  octave_idx_type ii = 0;\n\
+                  sm.cidx (0) = 0;\n\
+                  for (octave_idx_type j = 0; j < nc; ++j)\n\
+                   {\n\
+                     for (octave_idx_type i = 0; i < nr; ++i)\n\
+                       {\n\
+                         double tmp = A(i, j);\n\
+                         if (tmp != 0.)\n\
+                           {\n\
+                             if (ii == nz)\n\
+                               {\n\
+                                 nz = ceil ((nc * nr * ii) / (j * nr + 1));\n\
+                                 sm.change_capacity (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