Mercurial > fem-fenics-eugenio

changeset 68:4ef8dbafcdfc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New octave script for creation of function on the fly
author gedeone-octave <marco.vassallo@outlook.com>
date Sat, 27 Jul 2013 17:21:03 +0200
parents 5a781e79bee1
children 6dbe66dcf1c5
files inst/generate_fs.m inst/generate_lhs.m inst/generate_makefile.m inst/generate_rhs.m
diffstat 4 files changed, 404 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/generate_fs.m	Sat Jul 27 17:21:03 2013 +0200
@@ -0,0 +1,56 @@
+
+function output = generate_fs (ufl_name)
+
+STRING ="\n\
+#include ""functionspace.h""\n\
+#include ""mesh.h""\n\
+#include ""@@UFL_NAME@@.h""\n\
+\n\
+DEFUN_DLD (fem_fs_@@UFL_NAME@@, args, , ""initialize a fs from a mesh declared with fem_init_mesh"")\n\
+{\n\
+  int nargin = args.length ();\n\
+  octave_value retval;\n\
+\n\
+  if (nargin < 1 || nargin > 1)\n\
+    print_usage ();\n\
+  else\n\
+    {\n\
+\n\
+      if (! mesh_type_loaded)\n\
+        {\n\
+          mesh::register_type ();\n\
+          mesh_type_loaded = true;\n\
+          mlock ();\n\
+        }\n\
+\n\
+      if (args(0).type_id () == mesh::static_type_id ())\n\
+        {\n\
+          const mesh & msho = static_cast<const mesh&> (args(0).get_rep ());\n\
+          const dolfin::Mesh & mshd = msho.get_msh ();\n\
+          boost::shared_ptr <const dolfin::FunctionSpace> g (new @@UFL_NAME@@::FunctionSpace (mshd));\n\
+\n\
+          if (! functionspace_type_loaded)\n\
+            {\n\
+              functionspace::register_type ();\n\
+              functionspace_type_loaded = true;\n\
+              mlock ();\n\
+            }\n\
+\n\
+          retval = new functionspace(g);\n\
+        }\n\
+    }\n\
+  return retval;\n\
+}";
+
+STRING =  strrep (STRING, "@@UFL_NAME@@", ufl_name);
+
+fid = fopen (sprintf ("fem_fs_%s.cc", ufl_name), 'w');
+if (fid >= 0)
+  fputs (fid, STRING);
+  output = fclose (fid);
+else
+  error ("cannot open file");
+  output = 1;
+endif
+
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/generate_lhs.m	Sat Jul 27 17:21:03 2013 +0200
@@ -0,0 +1,133 @@
+
+function output = generate_lhs (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_lhs_@@UFL_NAME@@, args, , "" b = fem_lhs_@@UFL_NAME@@ (FUNCTIONAL SPACE, COEFF, BC)"")\n\
+{\n\
+  int nargin = args.length ();\n\
+  octave_value retval;\n\
+\n\
+  if (nargin < 1)\n\
+    print_usage ();\n\
+  else\n\
+    {\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@@::LinearForm L (V);\n\
+              std::size_t ncoef = L.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 = L.coefficient_number (cf.get_str ());\n\
+                      const boost::shared_ptr<const expression> & pexp = cf.get_expr ();\n\
+                      L.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 = L.coefficient_number (fun.get_str ());\n\
+                      const boost::shared_ptr<const dolfin::Function> & pfun = fun.get_pfun ();\n\
+                      L.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::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\
+                }\n\
+            }\n\
+        }\n\
+    }\n\
+  return retval;\n\
+}";
+
+STRING =  strrep (STRING, "@@UFL_NAME@@", ufl_name);
+
+fid = fopen (sprintf ("fem_lhs_%s.cc", ufl_name), 'w');
+if (fid >= 0)
+  fputs (fid, STRING);
+  output = fclose (fid);
+else
+  error ("cannot open file");
+  output = 1;
+endif
+
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/generate_makefile.m	Sat Jul 27 17:21:03 2013 +0200
@@ -0,0 +1,60 @@
+
+function output = generate_makefile (ufl_name, path)
+
+STRING ="\n\
+DIR = @@PATH@@\n\
+CPPFLAGS=\n\
+LIBS= -ldolfin\n\
+MKOCTFILE = mkoctfile\n\
+FFC = ffc\n\
+\n\
+OCTFILES = fem_fs_@@UFL_NAME@@.oct fem_rhs_@@UFL_NAME@@.oct fem_lhs_@@UFL_NAME@@.oct\n\
+\n\
+all : $(OCTFILES)\n\
+fs : fem_fs_@@UFL_NAME@@.oct\n\
+rhs : fem_rhs_@@UFL_NAME@@.oct\n\
+lhs : fem_lhs_@@UFL_NAME@@.oct\n\
+\n\
+fem_fs_@@UFL_NAME@@.o: fem_fs_@@UFL_NAME@@.cc\n\
+	$(MKOCTFILE) -I$(DIR) -I. $(CPPFLAGS) $(LDFLAGS) $< -c $@\n\
+\n\
+fem_fs_@@UFL_NAME@@.oct: @@UFL_NAME@@.h fem_fs_@@UFL_NAME@@.o\n\
+	$(MKOCTFILE) -s fem_fs_@@UFL_NAME@@.o $(DIR)fem_init_env.o $(LDFLAGS) $(LIBS) -o $@\n\
+\n\
+fem_rhs_@@UFL_NAME@@.o: fem_rhs_@@UFL_NAME@@.cc\n\
+	$(MKOCTFILE) -I$(DIR) -I. $(LDFLAGS) $(CPPFLAGS) $< -c $@\n\
+\n\
+fem_rhs_@@UFL_NAME@@.oct: @@UFL_NAME@@.h fem_rhs_@@UFL_NAME@@.o\n\
+	$(MKOCTFILE) -s fem_rhs_@@UFL_NAME@@.o $(DIR)fem_init_env.o $(LDFLAGS) $(LIBS) -o $@\n\
+\n\
+fem_lhs_@@UFL_NAME@@.o: fem_lhs_@@UFL_NAME@@.cc\n\
+	$(MKOCTFILE) -I$(DIR) -I. $(LDFLAGS) $(CPPFLAGS) $< -c $@\n\
+\n\
+fem_lhs_@@UFL_NAME@@.oct: @@UFL_NAME@@.h fem_lhs_@@UFL_NAME@@.o\n\
+	$(MKOCTFILE) -s fem_lhs_@@UFL_NAME@@.o $(DIR)fem_init_env.o $(LDFLAGS) $(LIBS) -o $@\n\
+\n\
+@@UFL_NAME@@.h: @@UFL_NAME@@.ufl\n\
+	$(FFC) -l dolfin @@UFL_NAME@@.ufl\n\
+\n\
+.PHONY: clean\n\
+\n\
+clean:\n\
+	 rm -f fem_fs_@@UFL_NAME@@.o fem_fs_@@UFL_NAME@@.cc @@UFL_NAME@@.h\n\
+	 rm -f fem_rhs_@@UFL_NAME@@.o fem_rhs_@@UFL_NAME@@.cc\n\
+	 rm -f fem_lhs_@@UFL_NAME@@.o fem_lhs_@@UFL_NAME@@.cc\n\
+	 rm -f Makefile_@@UFL_NAME@@\n\
+";
+
+STRING =  strrep (STRING, "@@UFL_NAME@@", ufl_name);
+STRING =  strrep (STRING, "@@PATH@@", path);
+
+fid = fopen (sprintf ("Makefile_%s", ufl_name), 'w');
+if (fid >= 0)
+  fputs (fid, STRING);
+  output = fclose (fid);
+else
+  error ("cannot open file");
+  output = 1;
+endif
+
+endfunction
--- /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