Mercurial > fem-fenics-eugenio

changeset 93:9e7035e0494b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Return an object of type form. * generate_rhs: create a bilinear form, setting the coefficients specified as input. The object returned is a fem-fenics::form. * generate_lhs: create a linear form, setting the coefficients specified as input. The object returned is a fem-fenics::form.
author gedeone-octave <marcovass89@hotmail.it>
date Wed, 07 Aug 2013 18:25:27 +0200
parents 2ecab16e88e1
children 4f688918f76a
files inst/generate_lhs.m inst/generate_rhs.m
diffstat 2 files changed, 22 insertions(+), 96 deletions(-) [+]
line wrap: on
line diff
--- a/inst/generate_lhs.m	Wed Aug 07 18:21:45 2013 +0200
+++ b/inst/generate_lhs.m	Wed Aug 07 18:25:27 2013 +0200
@@ -4,11 +4,11 @@
 STRING ="\n\
 #include ""@@UFL_NAME@@.h""\n\
 #include ""functionspace.h""\n\
-#include ""boundarycondition.h""\n\
+#include ""form.h""\n\
 #include ""coefficient.h""\n\
 #include ""function.h""\n\
 \n\
-DEFUN_DLD (@@UFL_NAME@@_LinearForm, args, , "" b = fem_lhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF, BC)"")\n\
+DEFUN_DLD (@@UFL_NAME@@_LinearForm, args, , "" b = fem_lhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF)"")\n\
 {\n\
   int nargin = args.length ();\n\
   octave_value retval;\n\
@@ -24,6 +24,13 @@
           functionspace_type_loaded = true;\n\
           mlock ();\n\
         }\n\
+      if (! form_type_loaded)\n\
+        {\n\
+          form::register_type ();\n\
+          form_type_loaded = true;\n\
+          mlock ();\n\
+        }\n\
+\n\
       if (args(0).type_id () == functionspace::static_type_id ())\n\
         {\n\
           const functionspace & fspo\n\
@@ -31,7 +38,7 @@
 \n\
           if (! error_state)\n\
             {\n\
-              const dolfin::FunctionSpace V = fspo.get_fsp ();\n\
+              const dolfin::FunctionSpace & V = fspo.get_fsp ();\n\
               @@UFL_NAME@@::LinearForm L (V);\n\
               std::size_t ncoef = L.num_coefficients (), nc = 0;\n\
 \n\
@@ -78,40 +85,7 @@
                 error (""Wrong number of coefficient"");\n\
               else\n\
                 {\n\
-\n\
-                  dolfin::Vector b;\n\
-                  dolfin::assemble (b, L);\n\
-\n\
-                  if (! boundarycondition_type_loaded)\n\
-                    {\n\
-                      boundarycondition::register_type ();\n\
-                      boundarycondition_type_loaded = true;\n\
-                      mlock ();\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(b);\n\
-                        }\n\
-                    }\n\
-\n\
-                  dim_vector dims;\n\
-                  dims.resize (2);\n\
-                  dims(0) = b.size ();\n\
-                  dims(1) = 1;\n\
-                  Array<double> myb (dims);\n\
-\n\
-                  for (std::size_t i = 0; i < b.size (); ++i)\n\
-                    myb(i) = b[i];\n\
-\n\
-                  retval = myb;\n\
+                  retval = new form (L);\n\
                 }\n\
             }\n\
         }\n\
--- a/inst/generate_rhs.m	Wed Aug 07 18:21:45 2013 +0200
+++ b/inst/generate_rhs.m	Wed Aug 07 18:25:27 2013 +0200
@@ -4,11 +4,11 @@
 STRING ="\n\
 #include ""@@UFL_NAME@@.h""\n\
 #include ""functionspace.h""\n\
-#include ""boundarycondition.h""\n\
+#include ""form.h""\n\
 #include ""coefficient.h""\n\
 #include ""function.h""\n\
 \n\
-DEFUN_DLD (@@UFL_NAME@@_BilinearForm, args, , ""A = fem_rhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF, BC)"")\n\
+DEFUN_DLD (@@UFL_NAME@@_BilinearForm, args, , ""A = fem_rhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF)"")\n\
 {\n\
 \n\
   int nargin = args.length ();\n\
@@ -24,6 +24,13 @@
           functionspace_type_loaded = true;\n\
           mlock ();\n\
         }\n\
+      if (! form_type_loaded)\n\
+        {\n\
+          form::register_type ();\n\
+          form_type_loaded = true;\n\
+          mlock ();\n\
+        }\n\
+\n\
       if (args(0).type_id () == functionspace::static_type_id ())\n\
         {\n\
           const functionspace & fspo\n\
@@ -31,7 +38,7 @@
 \n\
           if (! error_state)\n\
             {\n\
-              const dolfin::FunctionSpace V = fspo.get_fsp ();\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\
@@ -78,62 +85,7 @@
                 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.0005;\n\
-                  octave_idx_type nz = alpha * nr * nc;\n\
-                  SparseMatrix sm (nr, nc, 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.2 * ((nc * ii) / (j + 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\
+                  retval = new form (a);\n\
                 }\n\
             }\n\
         }\n\