Mercurial > fem-fenics-eugenio
changeset 68:4ef8dbafcdfc
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