# HG changeset patch
# User Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
# Date 1405161013 -7200
# Node ID 8ca45824938edf41e64b99f23ca3a2bb45a11170
# Parent  a2991e2a085b6fb4585122fde0600aa99d97db55
Add factories hierarchy for matrices' and vectors' assembly

diff -r a2991e2a085b -r 8ca45824938e src/Makefile.in
--- a/src/Makefile.in	Sat Jun 28 16:22:13 2014 +0200
+++ b/src/Makefile.in	Sat Jul 12 12:30:13 2014 +0200
@@ -17,6 +17,7 @@
 MKOCTFILE ?= mkoctfile
 FFC ?= ffc
 
+ARFLAGS=-r -s
 CPPFLAGS='@DENSE_CPPFLAGS@ @DOLFIN_CPPFLAGS@ @EIGEN_CPPFLAGS@ -I.'
 LIBS_RAW=@EIGEN_LIBS@ @DOLFIN_LIBS@
 LIBS=$(patsubst %, "%", $(LIBS_RAW))
@@ -73,16 +74,18 @@
 Function.o: Function.cc function.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c Function.cc -o $@
 
-assemble.oct: assemble.o
-	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble.o -o $@ $(LIBS)
+assemble.oct: assemble.o libfemfenics_factories.a
+	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble.o -o $@ \
+    libfemfenics_factories.a $(LIBS)
 
-assemble.o: assemble.cc form.h boundarycondition.h
+assemble.o: assemble.cc form.h boundarycondition.h femfenics_factory.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c assemble.cc -o $@
 
-assemble_system.oct: assemble_system.o
-	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble_system.o -o $@ $(LIBS)
+assemble_system.oct: assemble_system.o libfemfenics_factories.a
+	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) assemble_system.o -o $@ \
+    libfemfenics_factories.a $(LIBS)
 
-assemble_system.o: assemble_system.cc form.h boundarycondition.h
+assemble_system.o: assemble_system.cc form.h boundarycondition.h femfenics_factory.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c assemble_system.cc -o $@
 
 save.oct: save.o mkfunction
@@ -130,8 +133,18 @@
 interpolate.o: interpolate.cc function.h
 	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c interpolate.cc -o $@
 
+femfenics_factory.o: femfenics_factory.cc femfenics_factory.h femfenics_base_factory.h\
+                     uBLAS_factory.h
+	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@
+
+uBLAS_factory.o: uBLAS_factory.cc uBLAS_factory.h femfenics_base_factory.h
+	CPPFLAGS=$(CPPFLAGS) $(MKOCTFILE) -c $< -o $@
+
+libfemfenics_factories.a: femfenics_factory.o uBLAS_factory.o
+	$(AR) $(ARFLAGS) $@ $^
+
 clean:
-	$(RM) *.o core octave-core *.oct *~ *.xml
+	$(RM) *.o core octave-core *.oct *~ *.xml *.a
 
 cleanall:
 	$(RM) *.o core octave-core *.oct *~ *.xml *.status *.log \
diff -r a2991e2a085b -r 8ca45824938e src/assemble.cc
--- a/src/assemble.cc	Sat Jun 28 16:22:13 2014 +0200
+++ b/src/assemble.cc	Sat Jul 12 12:30:13 2014 +0200
@@ -17,6 +17,7 @@
 
 #include "form.h"
 #include "boundarycondition.h"
+#include "femfenics_factory.h"
 
 DEFUN_DLD (assemble, args, nargout, 
 "-*- texinfo -*-\n\
@@ -79,7 +80,8 @@
 
                   if (a.rank () == 2)
                     {
-                      dolfin::parameters["linear_algebra_backend"] = "uBLAS";
+                      femfenics_factory factory;
+
                       dolfin::Matrix A;
                       dolfin::assemble (A, a);
 
@@ -104,55 +106,13 @@
                             error ("assemble: unknown argument type");
                         }
 
-                      // 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;
-                      std::vector<std::size_t> cidx_tmp;
-
-                      dim_vector dims (nnz, 1);
-                      octave_idx_type nz = 0, ii = 0;
-                      Array<octave_idx_type> 
-                        ridx (dims, 0), 
-                        cidx (dims, 0);
-                      Array<double> data (dims, 0);
-
-                      octave_idx_type* orow = ridx.fortran_vec ();
-                      octave_idx_type* oc = cidx.fortran_vec ();
-                      double* ov = data.fortran_vec ();
-
-                      for (std::size_t i = 0; i < nr; ++i)
-                       {
-                         A.getrow (i, cidx_tmp, data_tmp);
-                         nz += cidx_tmp.size ();
-
-                         for (octave_idx_type j = 0; 
-                              j < cidx_tmp.size (); ++j)
-                           {
-                             orow [ii + j] = i;
-                             oc [ii + j] = cidx_tmp [j];
-                             ov [ii + j] = data_tmp [j];
-                           }
-
-                         ii = nz;
-                       }
-
-                      dims(0) = ii;
-                      ridx.resize (dims);
-                      cidx.resize (dims);
-                      data.resize (dims);
-
-                      SparseMatrix sm (data, ridx, cidx, nr, nc);
-                      retval(0) = sm;
+                      retval(0) = factory.matrix (A);
                     }
 
                   else if (a.rank () == 1)
                     {
+                      femfenics_factory factory;
+
                       dolfin::Vector A;
                       dolfin::assemble (A, a);
 
@@ -176,16 +136,7 @@
                             error ("assemble: unknown argument type");
                         }
 
-                      dim_vector dims;
-                      dims.resize (2);
-                      dims(0) = A.size ();
-                      dims(1) = 1;
-                      Array<double> myb (dims);
-
-                      for (std::size_t i = 0; i < A.size (); ++i)
-                        myb.xelem (i) = A[i];
-
-                      retval(0) = myb;
+                      retval(0) = factory.vector (A);
                     }
 
                   else if (a.rank () == 0)
@@ -219,6 +170,8 @@
 
                   if (a.rank () == 1)
                     {
+                      femfenics_factory factory;
+
                       dolfin::Vector A;
                       dolfin::assemble (A, a);
 
@@ -247,20 +200,8 @@
                             error ("assemble: unknown argument type");
                         }
 
-                      dim_vector dims;
-                      dims.resize (2);
-                      dims(0) = A.size ();
-                      dims(1) = 1;
-                      Array<double> myb (dims), myc (dims);
-
-                      for (std::size_t i = 0; i < A.size (); ++i)
-                        {
-                          myb.xelem (i) = A[i];
-                          myc.xelem (i) = x[i];
-                        }
-
-                      retval(0) = myb;
-                      retval(1) = myc;
+                      retval(0) = factory.vector (A);
+                      retval(1) = factory.vector (x);
                     }
 
                   else
diff -r a2991e2a085b -r 8ca45824938e src/assemble_system.cc
--- a/src/assemble_system.cc	Sat Jun 28 16:22:13 2014 +0200
+++ b/src/assemble_system.cc	Sat Jul 12 12:30:13 2014 +0200
@@ -17,6 +17,7 @@
 
 #include "form.h"
 #include "boundarycondition.h"
+#include "femfenics_factory.h"
 
 DEFUN_DLD (assemble_system, args, nargout,
 "-*- texinfo -*-\n\
@@ -79,7 +80,8 @@
 
                   if (a.rank () == 2 && b.rank () == 1)
                     {
-                      dolfin::parameters["linear_algebra_backend"] = "uBLAS";
+                      femfenics_factory factory;
+
                       dolfin::Matrix A;
                       dolfin::assemble (A, a);
                       dolfin::Vector B;
@@ -105,62 +107,8 @@
                             error ("assemble_system: unknown argument type");
                         }
 
-                      // 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;
-                      std::vector<std::size_t> cidx_tmp;
-
-                      dim_vector dims (nnz, 1);
-                      octave_idx_type nz = 0, ii = 0;
-                      Array<octave_idx_type> 
-                        ridx (dims, 0),
-                        cidx (dims, 0);
-                      Array<double> data (dims, 0);
-
-                      octave_idx_type* orow = ridx.fortran_vec ();
-                      octave_idx_type* oc = cidx.fortran_vec ();
-                      double* ov = data.fortran_vec ();
-
-                      for (std::size_t i = 0; i < nr; ++i)
-                       {
-                         A.getrow (i, cidx_tmp, data_tmp);
-                         nz += cidx_tmp.size ();
-
-                         for (octave_idx_type j = 0;
-                              j < cidx_tmp.size (); ++j)
-                           {
-                             orow [ii + j] = i;
-                             oc [ii + j] = cidx_tmp [j];
-                             ov [ii + j] = data_tmp [j];
-                           }
-
-                         ii = nz;
-                       }
-
-                      dims(0) = ii;
-                      ridx.resize (dims);
-                      cidx.resize (dims);
-                      data.resize (dims);
-
-                      SparseMatrix sm (data, ridx, cidx, nr, nc);
-                      retval(0) = sm;
-
-                      dim_vector dim;
-                      dim.resize (2);
-                      dim(0) = B.size ();
-                      dim(1) = 1;
-                      Array<double> myb (dim);
-
-                      for (std::size_t i = 0; i < B.size (); ++i)
-                        myb.xelem (i) = B[i];
-
-                      retval(1) = myb;
+                      retval(0) = factory.matrix (A);
+                      retval(1) = factory.vector (B);
                     }
                 }
               else
@@ -192,7 +140,8 @@
 
                   if (a.rank () == 2 && b.rank () == 1)
                     {
-                      dolfin::parameters["linear_algebra_backend"] = "uBLAS";
+                      femfenics_factory factory;
+
                       dolfin::Matrix A;
                       dolfin::assemble (A, a);
                       dolfin::Vector B;
@@ -224,65 +173,9 @@
                             error ("assemble_system: unknown argument type");
                         }
 
-                      // 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;
-                      std::vector<std::size_t> cidx_tmp;
-
-                      dim_vector dims (nnz, 1);
-                      octave_idx_type nz = 0, ii = 0;
-                      Array<octave_idx_type>
-                        ridx (dims, 0),
-                        cidx (dims, 0);
-                      Array<double> data (dims, 0);
-
-                      octave_idx_type* orow = ridx.fortran_vec ();
-                      octave_idx_type* oc = cidx.fortran_vec ();
-                      double* ov = data.fortran_vec ();
-
-                      for (std::size_t i = 0; i < nr; ++i)
-                       {
-                         A.getrow (i, cidx_tmp, data_tmp);
-                         nz += cidx_tmp.size ();
-
-                         for (octave_idx_type j = 0;
-                              j < cidx_tmp.size (); ++j)
-                           {
-                             orow [ii + j] = i;
-                             oc [ii + j] = cidx_tmp [j];
-                             ov [ii + j] = data_tmp [j];
-                           }
-
-                         ii = nz;
-                       }
-
-                      dims(0) = ii;
-                      ridx.resize (dims);
-                      cidx.resize (dims);
-                      data.resize (dims);
-
-                      SparseMatrix sm (data, ridx, cidx, nr, nc);
-                      retval(0) = sm;
-
-                      dim_vector dim;
-                      dim.resize (2);
-                      dim(0) = B.size ();
-                      dim(1) = 1;
-                      Array<double> myb (dim), myc (dim);
-
-                      for (std::size_t i = 0; i < B.size (); ++i)
-                        {
-                          myb.xelem (i) = B[i];
-                          myc.xelem (i) = x[i];
-                        }
-
-                      retval(1) = myb;
-                      retval(2) = myc;
+                      retval(0) = factory.matrix (A);
+                      retval(1) = factory.vector (B);
+                      retval(2) = factory.vector (x);
                     }
                 }
               else
diff -r a2991e2a085b -r 8ca45824938e src/femfenics_base_factory.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/femfenics_base_factory.h	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,35 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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 3 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/>.
+*/
+
+#ifndef __FEMFENICS_BASE_FACTORY__
+#define __FEMFENICS_BASE_FACTORY__
+
+#include <dolfin.h>
+#include <octave/oct.h>
+
+class femfenics_base_factory
+{
+  public:
+
+  femfenics_base_factory () {}
+  virtual ~femfenics_base_factory () {}
+
+  virtual octave_value matrix (dolfin::Matrix const&) const =0;
+  virtual octave_value vector (dolfin::Vector const&) const =0;
+};
+
+#endif
\ No newline at end of file
diff -r a2991e2a085b -r 8ca45824938e src/femfenics_factory.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/femfenics_factory.cc	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,27 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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 3 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 "femfenics_factory.h"
+#include "uBLAS_factory.h"
+
+femfenics_base_factory const&
+femfenics_factory::factory (void) const
+  {
+    //FIXME: Since just one backend has its interface implemented, this method
+    //       doesn't check which one to provide, yet
+    return uBLAS_factory::instance ();
+  }
\ No newline at end of file
diff -r a2991e2a085b -r 8ca45824938e src/femfenics_factory.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/femfenics_factory.h	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,46 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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 3 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/>.
+*/
+
+#ifndef __FEMFENICS_FACTORY__
+#define __FEMFENICS_FACTORY__
+
+#include "femfenics_base_factory.h"
+
+class femfenics_factory : public femfenics_base_factory
+{
+  public:
+
+  femfenics_factory ()
+    { dolfin::parameters["linear_algebra_backend"] = linear_algebra_backend (); }
+  virtual ~femfenics_factory () {}
+
+  //FIXME: just one backend implemented
+  inline char const * linear_algebra_backend (void) const
+    { return "uBLAS"; }
+
+  virtual inline octave_value matrix (dolfin::Matrix const& A) const
+    { return factory ().matrix (A); }
+
+  virtual inline octave_value vector (dolfin::Vector const& b) const
+    { return factory ().vector (b); }
+
+  private:
+
+  femfenics_base_factory const& factory (void) const;
+};
+
+#endif
\ No newline at end of file
diff -r a2991e2a085b -r 8ca45824938e src/uBLAS_factory.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/uBLAS_factory.cc	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,98 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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 3 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 "uBLAS_factory.h"
+
+femfenics_base_factory const&
+uBLAS_factory::instance (void)
+  {
+    static uBLAS_factory const theinstance;
+    return theinstance;
+  }
+
+octave_value
+uBLAS_factory::matrix (dolfin::Matrix const& A) const
+  {
+    octave_value retval;
+
+    // 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;
+    std::vector<std::size_t> cidx_tmp;
+
+    dim_vector dims (nnz, 1);
+    octave_idx_type nz = 0, ii = 0;
+    Array<octave_idx_type> 
+      ridx (dims, 0), 
+      cidx (dims, 0);
+    Array<double> data (dims, 0);
+
+    octave_idx_type* orow = ridx.fortran_vec ();
+    octave_idx_type* oc = cidx.fortran_vec ();
+    double* ov = data.fortran_vec ();
+
+    for (std::size_t i = 0; i < nr; ++i)
+      {
+        A.getrow (i, cidx_tmp, data_tmp);
+        nz += cidx_tmp.size ();
+
+        for (octave_idx_type j = 0; 
+             j < cidx_tmp.size (); ++j)
+          {
+            orow [ii + j] = i;
+            oc [ii + j] = cidx_tmp [j];
+            ov [ii + j] = data_tmp [j];
+          }
+
+        ii = nz;
+      }
+
+    dims(0) = ii;
+    ridx.resize (dims);
+    cidx.resize (dims);
+    data.resize (dims);
+
+    SparseMatrix sm (data, ridx, cidx, nr, nc);
+    retval = sm;
+
+    return retval;
+  }
+
+octave_value
+uBLAS_factory::vector (dolfin::Vector const& b) const
+  {
+    octave_value retval;
+
+    dim_vector dims;
+    dims.resize (2);
+    dims(0) = b.size ();
+    dims(1) = 1;
+    Array<double> myb (dims);
+
+    for (std::size_t i = 0; i < b.size (); ++i)
+      myb.xelem (i) = b[i];
+
+    retval = myb;
+
+    return retval;
+  }
\ No newline at end of file
diff -r a2991e2a085b -r 8ca45824938e src/uBLAS_factory.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/uBLAS_factory.h	Sat Jul 12 12:30:13 2014 +0200
@@ -0,0 +1,41 @@
+/*
+ Copyright (C) 2014 Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
+
+ 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 3 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/>.
+*/
+
+#ifndef __FEMFENICS_UBLAS_FACTORY__
+#define __FEMFENICS_UBLAS_FACTORY__
+
+#include "femfenics_base_factory.h"
+
+class uBLAS_factory : public femfenics_base_factory
+{
+  public:
+
+  virtual ~uBLAS_factory () {}
+
+  octave_value matrix (dolfin::Matrix const&) const;
+  octave_value vector (dolfin::Vector const&) const;
+
+  static femfenics_base_factory const& instance (void);
+
+  private:
+
+  uBLAS_factory () {}
+  uBLAS_factory (uBLAS_factory const&);
+  uBLAS_factory & operator = (uBLAS_factory const&);
+};
+
+#endif
\ No newline at end of file