--- a/src/assemble.cc
+++ b/src/assemble.cc
@@ -14,16 +14,19 @@
  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, "-*- texinfo -*-\n\
+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{form a} which is the form to assemble. \n\
+It can be a form of rank 2 (bilinear), a form of rank 1 (linear) or a form\n\
+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\
 @end itemize \n\
@@ -36,7 +39,7 @@
 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\
+at the example HyperElasticity.m\n\
 @seealso{BilinearForm, LinearForm, ResidualForm, JacobianForm}\n\
 @end deftypefn")
 {
@@ -80,14 +83,18 @@
 
                       for (std::size_t i = 1; i < nargin; ++i)
                         {
-                          if (args(i).type_id () == boundarycondition::static_type_id ())
+                          if (args(i).type_id () == 
+                              boundarycondition::static_type_id ())
                             {
                               const boundarycondition & bc
-                                = static_cast<const boundarycondition&> (args(i).get_rep ());
-
-                              const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > & pbc
-                                    = bc.get_bc ();
-
+                                = static_cast<const boundarycondition&> 
+                                (args(i).get_rep ());
+
+                              const 
+                                std::vector<boost::shared_ptr
+                                            <const dolfin::DirichletBC> > 
+                                & pbc = bc.get_bc ();
+                              
                               for (std::size_t j = 0; j < pbc.size (); ++j)
                                 pbc[j]->apply(A);
                             }
@@ -95,8 +102,12 @@
                             error ("assemble: unknown argument type");
                         }
 
-//                    It provides an upper boung for the nnz elements
-                      boost::tuples::tuple<const std::size_t*, const std::size_t*, const double*, int> aa = A.data ();
+                      // Get capacity of the dolfin sparse matrix
+                      boost::tuples::tuple<const std::size_t*, 
+                                           const std::size_t*, 
+                                           const double*, int> 
+                        aa = A.data ();
+
                       int nnz = aa.get<3> ();
                       std::size_t nr = A.size (0), nc = A.size (1);
                       std::vector<double> data_tmp;
@@ -104,7 +115,9 @@
 
                       dim_vector dims (nnz, 1);
                       octave_idx_type nz = 0, ii = 0;
-                      Array<octave_idx_type> ridx (dims, 0), cidx (dims, 0);
+                      Array<octave_idx_type> 
+                        ridx (dims, 0), 
+                        cidx (dims, 0);
                       Array<double> data (dims, 0);
 
                       octave_idx_type* orow = ridx.fortran_vec ();
@@ -116,7 +129,8 @@
                          A.getrow (i, cidx_tmp, data_tmp);
                          nz += cidx_tmp.size ();
 
-                         for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j)
+                         for (octave_idx_type j = 0; 
+                              j < cidx_tmp.size (); ++j)
                            {
                              orow [ii + j] = i;
                              oc [ii + j] = cidx_tmp [j];
@@ -142,13 +156,16 @@
 
                       for (std::size_t i = 1; i < nargin; ++i)
                         {
-                          if (args(i).type_id () == boundarycondition::static_type_id ())
+                          if (args(i).type_id () == 
+                              boundarycondition::static_type_id ())
                             {
                               const boundarycondition & bc
-                                = static_cast<const boundarycondition&> (args(i).get_rep ());
-
-                              const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > & pbc
-                                    = bc.get_bc ();
+                                = static_cast<const boundarycondition&> 
+                                (args(i).get_rep ());
+
+                              const std::vector<boost::shared_ptr 
+                                                <const dolfin::DirichletBC> >
+                                & pbc = bc.get_bc ();
 
                               for (std::size_t j = 0; j < pbc.size (); ++j)
                                 pbc[j]->apply(A);
@@ -189,7 +206,8 @@
         {
           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)
@@ -208,15 +226,19 @@
 
                       for (std::size_t i = 2; i < nargin; ++i)
                         {
-                          if (args(i).type_id () == boundarycondition::static_type_id ())
+                          if (args(i).type_id () == 
+                              boundarycondition::static_type_id ())
                             {
                               const boundarycondition & bc
-                                = static_cast<const boundarycondition&> (args(i).get_rep ());
-
-                              const std::vector<boost::shared_ptr <const dolfin::DirichletBC> > & pbc
-                                    = bc.get_bc ();
-
-                              for (std::size_t j = 0; j < pbc.size (); ++j)
+                                = static_cast<const boundarycondition&> 
+                                (args(i).get_rep ());
+
+                              const std::vector<boost::shared_ptr 
+                                                <const dolfin::DirichletBC> > 
+                                & pbc = bc.get_bc ();
+
+                              for (std::size_t j = 0; 
+                                   j < pbc.size (); ++j)
                                 pbc[j]->apply(A, x);
                             }
                           else