Mercurial > fem-fenics-eugenio
changeset 42:741a7999f590
New functions used to get the rhs and lhs of the variational problem.
* fem_rhs.cc: take as input the functionspace, a list of coefficients (optional)
a list of bc (optional) and return a sparse matrix
* fem_lhs.cc: as above burt return a vector
* Makefile.in: compile also fem_rhs.cc fem_lhs.cc
author | gedeone-octave <marco.vassallo@outlook.com> |
---|---|
date | Sat, 20 Jul 2013 23:49:23 +0200 |
parents | cedb2b5d6455 |
children | de8427de316e |
files | src/Makefile.in src/fem_lhs.cc src/fem_rhs.cc |
diffstat | 3 files changed, 255 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Makefile.in Sat Jul 20 23:46:03 2013 +0200 +++ b/src/Makefile.in Sat Jul 20 23:49:23 2013 +0200 @@ -1,10 +1,16 @@ MKOCTFILE ?= mkoctfile -OCTFILES=fem_init_env.oct fem_init_mesh.oct fem_get_mesh.oct fem_fs.oct fem_bc.oct fem_coeff.oct +OCTFILES=fem_init_env.oct\ + fem_init_mesh.oct \ + fem_get_mesh.oct \ + fem_fs.oct \ + fem_bc.oct \ + fem_coeff.oct \ + fem_lhs.oct \ + fem_rhs.oct \ LIBS += -ldolfin - all: $(OCTFILES) fem_init_env.oct: fem_init_env.o @@ -46,6 +52,18 @@ fem_coeff.o: fem_coeff.cc expression.h coefficient.h $(MKOCTFILE) $(CPPFLAGS) -c fem_coeff.cc $(LDFLAGS) -o $@ -I. +fem_lhs.oct: fem_lhs.o fem_init_env.o + $(MKOCTFILE) $(CPPFLAGS) -s fem_init_env.o fem_lhs.o -o $@ $(LDFLAGS) $(LIBS) + +fem_lhs.o: fem_lhs.cc expression.h Laplace.h + $(MKOCTFILE) $(CPPFLAGS) -c fem_lhs.cc $(LDFLAGS) -o $@ -I. + +fem_rhs.oct: fem_rhs.o fem_init_env.o + $(MKOCTFILE) $(CPPFLAGS) -s fem_init_env.o fem_rhs.o -o $@ $(LDFLAGS) $(LIBS) + +fem_rhs.o: fem_rhs.cc expression.h Laplace.h + $(MKOCTFILE) $(CPPFLAGS) -c fem_rhs.cc $(LDFLAGS) -o $@ -I. + Laplace.h: Laplace.ufl ffc -l dolfin Laplace.ufl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/fem_lhs.cc Sat Jul 20 23:49:23 2013 +0200 @@ -0,0 +1,110 @@ +/* + Copyright (C) 2013 Marco Vassallo + + 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 2 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 "Laplace.h" +#include <dolfin.h> +#include "functionspace.h" +#include "boundarycondition.h" +#include "coefficient.h" + +DEFUN_DLD (fem_lhs, args, , "-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{bc}]} = \ +fem_rhs (@var{Functional Space}, @var{Coefficient}, \ +@var{Boundary Condition}) \n\ +The input parameters are\n\ +@itemize @bullet \n\ +@item @var{Functional Space} is the fem-fenics functional space where\ +the bilinear form is defined\n\ +@item @var{Boundary Condition} contains the value of the essential bc that\ +we want to apply to the Bilinear Form\n\ +@item @var{Coefficient} is a variable of type coefficient which contains\ +the value of the coefficient for the bilinear form\n\ +The output @var{A} is a sparse matrix which represents the bilinear form\n\ +@seealso{fem_init_mesh, fem_fs}\n\ +@end deftypefn") +{ + + int nargin = args.length (); + octave_value retval; + + if (nargin < 1) + print_usage (); + else + { + if (args(0).type_id () == functionspace::static_type_id ()) + { + const functionspace & fspo + = static_cast<const functionspace&> (args(0).get_rep ()); + + if (! error_state) + { + const dolfin::FunctionSpace V = fspo.get_fsp (); + Laplace::LinearForm L (V); + std::size_t ncoef = L.num_coefficients (), nc = 0; + + for (std::size_t i = 1; i < ncoef + 1; ++i) + { + if (args(i).type_id () == coefficient::static_type_id ()) + { + const coefficient & cf + = static_cast <const coefficient&> (args(i).get_rep ()); + + std::size_t n = L.coefficient_number (cf.get_str ()); + const boost::shared_ptr<const expression> & pexp = cf.get_expr (); + L.set_coefficient (n, pexp); + ++nc; + } + } + + if (nc != ncoef) + error ("Wrong number of coefficient"); + else + { + + dolfin::Vector b; + dolfin::assemble (b, L); + + for (std::size_t i = 1 + nc; i < nargin; ++i) + { + 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) + pbc[j]->apply(b); + } + } + + 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(i) = b[i]; + + retval = myb; + } + } + } + } + return retval; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/fem_rhs.cc Sat Jul 20 23:49:23 2013 +0200 @@ -0,0 +1,125 @@ +/* + Copyright (C) 2013 Marco Vassallo + + 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 2 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 "Laplace.h" +#include <dolfin.h> +#include "functionspace.h" +#include "boundarycondition.h" +#include "coefficient.h" + +DEFUN_DLD (fem_rhs, args, , "-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{bc}]} = \ +fem_rhs (@var{Functional Space}, @var{Coefficient}, \ +@var{Boundary Condition}) \n\ +The input parameters are\n\ +@itemize @bullet \n\ +@item @var{Functional Space} is the fem-fenics functional space where\ +the bilinear form is defined\n\ +@item @var{Boundary Condition} contains the value of the essential bc that\ +we want to apply to the Bilinear Form\n\ +@item @var{Coefficient} is a variable of type coefficient which contains\ +the value of the coefficient for the bilinear form\n\ +The output @var{A} is a sparse matrix which represents the bilinear form\n\ +@seealso{fem_init_mesh, fem_fs}\n\ +@end deftypefn") +{ + + int nargin = args.length (); + octave_value retval; + + if (nargin < 1) + print_usage (); + else + { + if (args(0).type_id () == functionspace::static_type_id ()) + { + const functionspace & fspo + = static_cast<const functionspace&> (args(0).get_rep ()); + + if (! error_state) + { + const dolfin::FunctionSpace V = fspo.get_fsp (); + Laplace::BilinearForm a (V, V); + std::size_t ncoef = a.num_coefficients (), nc = 0; + + for (std::size_t i = 1; i < ncoef + 1; ++i) + { + if (args(i).type_id () == coefficient::static_type_id ()) + { + const coefficient & cf + = static_cast <const coefficient&> (args(i).get_rep ()); + + std::size_t n = a.coefficient_number (cf.get_str ()); + const boost::shared_ptr<const expression> & pexp = cf.get_expr (); + a.set_coefficient (n, pexp); + ++nc; + } + } + + if (nc != ncoef) + error ("Wrong number of coefficient"); + else + { + + dolfin::Matrix A; + dolfin::assemble (A, a); + + for (std::size_t i = 1 + nc; i < nargin; ++i) + { + 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) + pbc[j]->apply(A); + } + } + + int nr = A.size (0), nc = A.size (1); + // nz shoud be estimated in a better way + int nz = nr * nc; + SparseMatrix sm (nr, nc, nz); + + int ii = 0; + sm.cidx (0) = 0; + for (int j = 0; j < nc; ++j) + { + for (int i = 0; i < nr; i++) + { + double tmp = A(i, j); + if (tmp != 0.) + { + sm.data(ii) = tmp; + sm.ridx(ii) = i; + ii++; + } + } + sm.cidx(j+1) = ii; + } + sm.maybe_compress (); + + retval = sm; + + } + } + } + } + return retval; +}