Mercurial > fem-fenics-eugenio
changeset 91:51bfdf30b822
New DLD functions for assembling a matrix/vector from a bilinear/linear form
* assemble.cc: assemble a form, applying the bc specified. Simmetry not
mantained.
* assemble_system.cc: assemble both the bilinear and the linear form at once,
simmetry is mantained (if any).
| author | gedeone-octave <marcovass89@hotmail.it> |
|---|---|
| date | Wed, 07 Aug 2013 18:20:59 +0200 |
| parents | 7cd7bd1fc2b5 |
| children | 2ecab16e88e1 |
| files | src/assemble.cc src/assemble_system.cc |
| diffstat | 2 files changed, 248 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/assemble.cc Wed Aug 07 18:20:59 2013 +0200 @@ -0,0 +1,135 @@ +#include "form.h" +#include "boundarycondition.h" + +DEFUN_DLD (assemble, args, , "A = assemble (FORM, BC)") +{ + + int nargin = args.length (); + octave_value retval; + + if (nargin < 1) + print_usage (); + else + { + if (! form_type_loaded) + { + form::register_type (); + form_type_loaded = true; + mlock (); + } + if (args(0).type_id () == form::static_type_id ()) + { + const form & frm + = static_cast<const form&> (args(0).get_rep ()); + + if (! error_state) + { + const dolfin::Form & a = frm.get_form (); + a.check (); + + if (a.rank () == 2) + { + dolfin::Matrix A; + dolfin::assemble (A, a); + if (! boundarycondition_type_loaded) + { + boundarycondition::register_type (); + boundarycondition_type_loaded = true; + mlock (); + } + + for (std::size_t i = 1; 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); + } + else + error ("assemble: unknown argument type"); + } + + std::size_t nr = A.size (0), nc = A.size (1); + // nz shoud be estimated in a better way + double alpha = 0.0005; + octave_idx_type nz = alpha * nr * nc; + SparseMatrix sm (nr, nc, nz); + + std::size_t ii = 0; + sm.cidx (0) = 0; + for (std::size_t j = 0; j < nc; ++j) + { + for (std::size_t i = 0; i < nr; ++i) + { + double tmp = A(i, j); + if (tmp != 0.) + { + if (ii == nz) + { + nz = 1.2 * ((nc * ii) / (j + 1)); + sm.change_capacity (nz); + } + sm.data(ii) = tmp; + sm.ridx(ii) = i; + ++ii; + } + } + sm.cidx(j+1) = ii; + } + sm.maybe_compress (); + retval = sm; + } + else if (a.rank () == 1) + { + dolfin::Vector A; + dolfin::assemble (A, a); + + if (! boundarycondition_type_loaded) + { + boundarycondition::register_type (); + boundarycondition_type_loaded = true; + mlock (); + } + + for (std::size_t i = 1; 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); + } + else + 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(i) = A[i]; + + retval = myb; + } + + else + error ("assemble: unknown size"); + } + } + } + return retval; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/assemble_system.cc Wed Aug 07 18:20:59 2013 +0200 @@ -0,0 +1,113 @@ +#include "form.h" +#include "boundarycondition.h" + +DEFUN_DLD (assemble_system, args, , "A = assemble (FORM, BC)") +{ + + int nargin = args.length (); + octave_value_list retval; + + if (nargin < 2) + print_usage (); + else + { + if (! form_type_loaded) + { + form::register_type (); + form_type_loaded = true; + mlock (); + } + if (! boundarycondition_type_loaded) + { + boundarycondition::register_type (); + boundarycondition_type_loaded = true; + mlock (); + } + + if (args(0).type_id () == form::static_type_id () + && args(1).type_id () == form::static_type_id ()) + { + const form & frm1 + = static_cast<const form&> (args(0).get_rep ()); + + const form & frm2 + = static_cast<const form&> (args(1).get_rep ()); + + std::vector<const dolfin::DirichletBC*> bcs (0); + for (std::size_t i = 2; 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) + bcs.push_back (& *(pbc[j])); + } + } + + + if (! error_state) + { + const dolfin::Form & a = frm1.get_form (); + a.check (); + const dolfin::Form & L = frm2.get_form (); + L.check (); + + if (a.rank () == 2 && L.rank () == 1) + { + dolfin::Matrix A; + dolfin::Vector b; + dolfin::assemble_system (A, b, a, L, bcs); + + std::size_t nr = A.size (0), nc = A.size (1); + // nz shoud be estimated in a better way + double alpha = 0.0005; + octave_idx_type nz = alpha * nr * nc; + SparseMatrix sm (nr, nc, nz); + + std::size_t ii = 0; + sm.cidx (0) = 0; + for (std::size_t j = 0; j < nc; ++j) + { + for (std::size_t i = 0; i < nr; ++i) + { + double tmp = A(i, j); + if (tmp != 0.) + { + if (ii == nz) + { + nz = 1.2 * ((nc * ii) / (j + 1)); + sm.change_capacity (nz); + } + sm.data(ii) = tmp; + sm.ridx(ii) = i; + ++ii; + } + } + sm.cidx(j+1) = ii; + } + sm.maybe_compress (); + retval (0) = sm; + + 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(1) = myb; + } + + else + error ("assemble: unknown size"); + } + } + } + return retval; +}