--- a/src/assemble.cc	Tue Aug 27 21:45:15 2013 +0200
+++ b/src/assemble.cc	Tue Aug 27 22:32:57 2013 +0200
@@ -1,28 +1,70 @@
+/*
+ Copyright (C) 2013 Marco Vassallo
+
+ This program 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.
+
+ 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 General Public License for more
+ details.
+
+ You should have received a copy of the GNU General Public License along with
+ this program; if not, see <http://www.gnu.org/licenses/>.
+*/
 #include "form.h"
 #include "boundarycondition.h"
 
-DEFUN_DLD (assemble, args, nargout, "A = assemble (FORM, BC)")
+DEFUN_DLD (assemble, args, nargout, "-*- texinfo -*-\n\
+@deftypefn {Function File} {[@var{A}], [@var{x}(Optional)]} = \
+assemble (@var{form a}, @var{DirichletBC}(Optional), @var{...}) \n\
+The input arguments are\n\
+@itemize @bullet\n\
+@item @var{form a} which is the form to assemble. It can be a form of rank 2\n\
+(bilinear), a form of rank 1 (linear) or a form of rank 0 (functional).\n\
+@item @var{DirichletBC} represents the optional BC that you wish to apply to\n\
+the system. If more than one BC has to be applied, just list them.\n\
+@enditemize \n\
+The output @var{A} is a discretized representation of the @var{form a}:\n\
+@itemize @bullet\n\
+@item @var{A} is a sparse Matrix if @var{form a} is a bilinear form\n\
+@item @var{A} is a Vector if @var{form a} is a linear form\n\
+@item @var{A} is a Double if @var{form a} is a functional\n\
+@enditemize \n\
+@If you need to apply boundary condition to a vector for a nonlinear problem \n\
+@then you should provide as 2nd argument the vector and you will receive it back\n\
+@as the second output argument. For an example of this situation, you can look\n\
+@the example HyperElasticity.m\n\
+@end deftypefn")
 {
-
   int nargin = args.length ();
   octave_value_list retval;
 
+  if (! boundarycondition_type_loaded)
+    {
+      boundarycondition::register_type ();
+      boundarycondition_type_loaded = true;
+      mlock ();
+    }
+
+  if (! form_type_loaded)
+    {
+      form::register_type ();
+      form_type_loaded = true;
+      mlock ();
+    }
+
   if (nargout == 1)
     {
       if (nargin < 1)
         print_usage ();
       else
         {
-          if (! form_type_loaded)
-            {
-              form::register_type ();
-              form_type_loaded = true;
-              mlock ();
-            }
           if (args(0).type_id () == form::static_type_id ())
             {
-              const form & frm
-                = static_cast<const form&> (args(0).get_rep ());
+              const form & frm = static_cast<const form&> (args(0).get_rep ());
 
               if (! error_state)
                 {
@@ -31,14 +73,10 @@
 
                   if (a.rank () == 2)
                     {
+                      std::cout << "Assembling Matrix from the bilinear form..."
+                        << std::endl;
                       dolfin::Matrix A;
                       dolfin::assemble (A, a);
-                      if (! boundarycondition_type_loaded)
-                        {
-                          boundarycondition::register_type ();
-                          boundarycondition_type_loaded = true;
-                          mlock ();
-                        }
 
                       for (std::size_t i = 1; i < nargin; ++i)
                         {
@@ -57,7 +95,6 @@
                             error ("assemble: unknown argument type");
                         }
 
-
                       std::size_t nr = A.size (0), nc = A.size (1);
                       std::vector<double> data_tmp;
                       std::vector<std::size_t> cidx_tmp;
@@ -87,7 +124,6 @@
                          ii = nz;
                        }
 
-
                       ridx.resize (ii);
                       cidx.resize (ii);
                       data.resize (ii);
@@ -95,18 +131,14 @@
                       SparseMatrix sm (data, ridx, cidx, nr, nc);
                       retval(0) = sm;
                     }
+
                   else if (a.rank () == 1)
                     {
+                      std::cout << "Assembling Vector from the linear form..."
+                        << std::endl;
                       dolfin::Vector A;
                       dolfin::assemble (A, a);
 
-                      if (! boundarycondition_type_loaded)
-                        {
-                          boundarycondition::register_type ();
-                          boundarycondition_type_loaded = true;
-                          mlock ();
-                        }
-
                       for (std::size_t i = 1; i < nargin; ++i)
                         {
                           if (args(i).type_id () == boundarycondition::static_type_id ())
@@ -135,34 +167,33 @@
 
                       retval(0) = myb;
                     }
+
                   else if (a.rank () == 0)
                     {
+                      std::cout << "Assembling double from the functional form..."
+                        << std::endl;
                       double b = dolfin::assemble (a);
                       retval(0) = octave_value (b);
                     }
 
                   else
-                    error ("assemble: unknown size");
+                    error ("assemble: unknown form size");
                 }
             }
         }
     }
   else if (nargout == 2)
     {
+      std::cout << "Assemble: apply boundary condition to a vector for a nonlinear problem..."
+        << std::endl;
+
       if (nargin < 2)
         print_usage ();
       else
         {
-          if (! form_type_loaded)
-            {
-              form::register_type ();
-              form_type_loaded = true;
-              mlock ();
-            }
           if (args(0).type_id () == form::static_type_id ())
             {
-              const form & frm
-                = static_cast<const form&> (args(0).get_rep ());
+              const form & frm = static_cast<const form &> (args(0).get_rep ());
               const Array<double> myx = args(1).array_value ();
 
               if (! error_state)
@@ -172,6 +203,8 @@
 
                   if (a.rank () == 1)
                     {
+                      std::cout << "Assembling Vector from the linear form..."
+                        << std::endl;
                       dolfin::Vector A;
                       dolfin::assemble (A, a);
 
@@ -179,13 +212,6 @@
                       for (std::size_t i = 0; i < myx.length (); ++i)
                         x.setitem (i, myx.xelem (i));
 
-                      if (! boundarycondition_type_loaded)
-                        {
-                          boundarycondition::register_type ();
-                          boundarycondition_type_loaded = true;
-                          mlock ();
-                        }
-
                       for (std::size_t i = 2; i < nargin; ++i)
                         {
                           if (args(i).type_id () == boundarycondition::static_type_id ())
@@ -200,7 +226,7 @@
                                 pbc[j]->apply(A, x);
                             }
                           else
-                            error ("assemble NL vector: unknown argument type");
+                            error ("assemble: unknown argument type");
                         }
 
                       dim_vector dims;
@@ -220,7 +246,7 @@
                     }
 
                   else
-                    error ("assemble NL vector: unknown size");
+                    error ("assemble: unknown size");
                 }
             }
         }