Mercurial > fem-fenics-eugenio
changeset 116:77eefe47f7ab
Maint: improve the txt message
author | gedeone-octave <marcovass89@hotmail.it> |
---|---|
date | Tue, 27 Aug 2013 22:32:57 +0200 |
parents | 41614aa82191 |
children | 70322135b25f |
files | src/assemble.cc src/assemble_system.cc src/boundarycondition.h |
diffstat | 3 files changed, 225 insertions(+), 193 deletions(-) [+] |
line wrap: on
line diff
--- a/src/assemble.cc Tue Aug 27 21:45:15 2013 +0200 +++ b/src/assemble.cc Tue Aug 27 22:32:57 2013 +0200 @@ -1,28 +1,70 @@ +/* + 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 "form.h" #include "boundarycondition.h" -DEFUN_DLD (assemble, args, nargout, "A = assemble (FORM, BC)") +DEFUN_DLD (assemble, args, nargout, "-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{A}], [@var{x}(Optional)]} = \ +assemble (@var{form a}, @var{DirichletBC}(Optional), @var{...}) \n\ +The input arguments are\n\ +@itemize @bullet\n\ +@item @var{form a} which is the form to assemble. It can be a form of rank 2\n\ +(bilinear), a form of rank 1 (linear) or a form of rank 0 (functional).\n\ +@item @var{DirichletBC} represents the optional BC that you wish to apply to\n\ +the system. If more than one BC has to be applied, just list them.\n\ +@enditemize \n\ +The output @var{A} is a discretized representation of the @var{form a}:\n\ +@itemize @bullet\n\ +@item @var{A} is a sparse Matrix if @var{form a} is a bilinear form\n\ +@item @var{A} is a Vector if @var{form a} is a linear form\n\ +@item @var{A} is a Double if @var{form a} is a functional\n\ +@enditemize \n\ +@If you need to apply boundary condition to a vector for a nonlinear problem \n\ +@then you should provide as 2nd argument the vector and you will receive it back\n\ +@as the second output argument. For an example of this situation, you can look\n\ +@the example HyperElasticity.m\n\ +@end deftypefn") { - int nargin = args.length (); octave_value_list retval; + if (! boundarycondition_type_loaded) + { + boundarycondition::register_type (); + boundarycondition_type_loaded = true; + mlock (); + } + + if (! form_type_loaded) + { + form::register_type (); + form_type_loaded = true; + mlock (); + } + if (nargout == 1) { 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 ()); + const form & frm = static_cast<const form&> (args(0).get_rep ()); if (! error_state) { @@ -31,14 +73,10 @@ if (a.rank () == 2) { + std::cout << "Assembling Matrix from the bilinear form..." + << std::endl; 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) { @@ -57,7 +95,6 @@ error ("assemble: unknown argument type"); } - std::size_t nr = A.size (0), nc = A.size (1); std::vector<double> data_tmp; std::vector<std::size_t> cidx_tmp; @@ -87,7 +124,6 @@ ii = nz; } - ridx.resize (ii); cidx.resize (ii); data.resize (ii); @@ -95,18 +131,14 @@ SparseMatrix sm (data, ridx, cidx, nr, nc); retval(0) = sm; } + else if (a.rank () == 1) { + std::cout << "Assembling Vector from the linear form..." + << std::endl; 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 ()) @@ -135,34 +167,33 @@ retval(0) = myb; } + else if (a.rank () == 0) { + std::cout << "Assembling double from the functional form..." + << std::endl; double b = dolfin::assemble (a); retval(0) = octave_value (b); } else - error ("assemble: unknown size"); + error ("assemble: unknown form size"); } } } } else if (nargout == 2) { + std::cout << "Assemble: apply boundary condition to a vector for a nonlinear problem..." + << std::endl; + if (nargin < 2) 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 ()); + const form & frm = static_cast<const form &> (args(0).get_rep ()); const Array<double> myx = args(1).array_value (); if (! error_state) @@ -172,6 +203,8 @@ if (a.rank () == 1) { + std::cout << "Assembling Vector from the linear form..." + << std::endl; dolfin::Vector A; dolfin::assemble (A, a); @@ -179,13 +212,6 @@ for (std::size_t i = 0; i < myx.length (); ++i) x.setitem (i, myx.xelem (i)); - if (! boundarycondition_type_loaded) - { - boundarycondition::register_type (); - boundarycondition_type_loaded = true; - mlock (); - } - for (std::size_t i = 2; i < nargin; ++i) { if (args(i).type_id () == boundarycondition::static_type_id ()) @@ -200,7 +226,7 @@ pbc[j]->apply(A, x); } else - error ("assemble NL vector: unknown argument type"); + error ("assemble: unknown argument type"); } dim_vector dims; @@ -220,7 +246,7 @@ } else - error ("assemble NL vector: unknown size"); + error ("assemble: unknown size"); } } }
--- a/src/assemble_system.cc Tue Aug 27 21:45:15 2013 +0200 +++ b/src/assemble_system.cc Tue Aug 27 22:32:57 2013 +0200 @@ -1,11 +1,59 @@ +/* + 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 "form.h" #include "boundarycondition.h" -DEFUN_DLD (assemble_system, args, nargout, "A = assemble (FORM, BC)") +DEFUN_DLD (assemble_system, args, nargout, "-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{A}], [@var{b}], [@var{x}(Optional)]} = \ +assemble_system (@var{form a}, @var{form L}, @var{DirichletBC}(Optional), @var{...}) \n\ +The input arguments are\n\ +@itemize @bullet\n\ +@item @var{form a} the bilinear form to assemble.\n\ +@item @var{form a} the linear form to assemble.\n\ +@item @var{DirichletBC} represents the optional BC that you wish to apply to\n\ +the system. If more than one BC has to be applied, just list them.\n\ +@enditemize \n\ +The output @var{A} is a discretized representation of the system:\n\ +@itemize @bullet\n\ +@item @var{A} is the sparse Matrix corresponding to the @var{form a}\n\ +@item @var{A} is the Vector corresponding to the @var{form L}\n\ +@enditemize \n\ +@If you need to apply boundary condition to a system for a nonlinear problem \n\ +@then you should provide as 3rd argument the vector and you will receive it back\n\ +@as the third output argument. For an example of this situation, you can look\n\ +@the example HyperElasticity.m\n\ +@end deftypefn") { - int nargin = args.length (); octave_value_list retval; + 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 (nargout == 2) { @@ -13,20 +61,12 @@ print_usage (); else { - if (! form_type_loaded) - { - form::register_type (); - form_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 ()); + const form & frm1 = static_cast<const form&> (args(0).get_rep ()); + const form & frm2 = static_cast<const form&> (args(1).get_rep ()); if (! error_state) { @@ -42,13 +82,6 @@ dolfin::Vector B; dolfin::assemble (B, b); - if (! boundarycondition_type_loaded) - { - boundarycondition::register_type (); - boundarycondition_type_loaded = true; - mlock (); - } - for (std::size_t i = 2; i < nargin; ++i) { if (args(i).type_id () == boundarycondition::static_type_id ()) @@ -63,9 +96,110 @@ pbc[j]->apply(A, B); } else - error ("assemble: unknown argument type"); + error ("assemble_system: unknown argument type"); } + std::size_t nr = A.size (0), nc = A.size (1); + std::vector<double> data_tmp; + std::vector<std::size_t> cidx_tmp; + + octave_idx_type dims = A.size (0), nz = 0, ii = 0; + ColumnVector ridx (dims), cidx (dims), data (dims); + + for (std::size_t i = 0; i < nr; ++i) + { + A.getrow (i, cidx_tmp, data_tmp); + nz += cidx_tmp.size (); + if (dims < nz) + { + dims = 1.2 * ((nr * nz) / (i + 1));; + ridx.resize (dims); + cidx.resize (dims); + data.resize (dims); + } + + for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j) + { + ridx.xelem (ii + j) = i + 1; + cidx.xelem (ii + j) = cidx_tmp [j] + 1; + data.xelem (ii + j) = data_tmp [j]; + } + ii = nz; + } + + ridx.resize (ii); + cidx.resize (ii); + data.resize (ii); + + SparseMatrix sm (data, ridx, cidx, nr, nc); + retval(0) = sm; + + dim_vector dim; + dim.resize (2); + dim(0) = B.size (); + dim(1) = 1; + Array<double> myb (dim); + + for (std::size_t i = 0; i < B.size (); ++i) + myb.xelem (i) = B[i]; + + retval(1) = myb; + } + } + else + error ("assemble_system: unknown size"); + } + } + } + else if (nargout == 3) + { + std::cout << "Assemble_system: apply boundary condition to a vector for a nonlinear problem..." + << std::endl; + if (nargin < 3) + print_usage (); + else + { + 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 ()); + const Array<double> myx = args(2).array_value (); + + if (! error_state) + { + const dolfin::Form & a = frm1.get_form (); + const dolfin::Form & b = frm2.get_form (); + a.check (); + b.check (); + + if (a.rank () == 2 && b.rank () == 1) + { + dolfin::Matrix A; + dolfin::assemble (A, a); + dolfin::Vector B; + dolfin::assemble (B, b); + dolfin::Vector x (myx.length ()); + + for (std::size_t i = 0; i < myx.length (); ++i) + x.setitem (i, myx.xelem (i)); + + for (std::size_t i = 3; 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, B, x); + } + else + error ("assemble_system: unknown argument type"); + } std::size_t nr = A.size (0), nc = A.size (1); std::vector<double> data_tmp; @@ -96,7 +230,6 @@ ii = nz; } - ridx.resize (ii); cidx.resize (ii); data.resize (ii); @@ -104,129 +237,6 @@ SparseMatrix sm (data, ridx, cidx, nr, nc); retval(0) = sm; - - dim_vector dim; - dim.resize (2); - dim(0) = B.size (); - dim(1) = 1; - Array<double> myb (dim); - - for (std::size_t i = 0; i < B.size (); ++i) - myb.xelem (i) = B[i]; - - retval(1) = myb; - } - } - else - error ("assemble: unknown size"); - } - } - } - else if (nargout == 3) - { - if (nargin < 3) - print_usage (); - else - { - if (! form_type_loaded) - { - form::register_type (); - form_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 ()); - - const Array<double> myx = args(2).array_value (); - - if (! error_state) - { - const dolfin::Form & a = frm1.get_form (); - const dolfin::Form & b = frm2.get_form (); - a.check (); - b.check (); - - if (a.rank () == 2 && b.rank () == 1) - { - dolfin::Matrix A; - dolfin::assemble (A, a); - dolfin::Vector B; - dolfin::assemble (B, b); - - dolfin::Vector x (myx.length ()); - - for (std::size_t i = 0; i < myx.length (); ++i) - x.setitem (i, myx.xelem (i)); - - if (! boundarycondition_type_loaded) - { - boundarycondition::register_type (); - boundarycondition_type_loaded = true; - mlock (); - } - - for (std::size_t i = 3; 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, B, x); - } - else - error ("assemble: unknown argument type"); - } - - - std::size_t nr = A.size (0), nc = A.size (1); - std::vector<double> data_tmp; - std::vector<std::size_t> cidx_tmp; - - octave_idx_type dims = A.size (0), nz = 0, ii = 0; - ColumnVector ridx (dims), cidx (dims), data (dims); - - for (std::size_t i = 0; i < nr; ++i) - { - A.getrow (i, cidx_tmp, data_tmp); - nz += cidx_tmp.size (); - if (dims < nz) - { - dims = 1.2 * ((nr * nz) / (i + 1));; - ridx.resize (dims); - cidx.resize (dims); - data.resize (dims); - } - - for (octave_idx_type j = 0; j < cidx_tmp.size (); ++j) - { - ridx.xelem (ii + j) = i + 1; - cidx.xelem (ii + j) = cidx_tmp [j] + 1; - data.xelem (ii + j) = data_tmp [j]; - } - - ii = nz; - } - - - ridx.resize (ii); - cidx.resize (ii); - data.resize (ii); - - SparseMatrix sm (data, ridx, cidx, nr, nc); - retval(0) = sm; - - dim_vector dim; dim.resize (2); dim(0) = B.size (); @@ -244,7 +254,7 @@ } } else - error ("assemble: unknown size"); + error ("assemble_system: unknown size"); } } }
--- a/src/boundarycondition.h Tue Aug 27 21:45:15 2013 +0200 +++ b/src/boundarycondition.h Tue Aug 27 22:32:57 2013 +0200 @@ -36,17 +36,13 @@ os << "Boundary condition : " << bcu[i]->str (true) << std::endl; } - ~boundarycondition (void) {} - bool is_defined (void) const { return true; } - const std::vector< boost::shared_ptr <const dolfin::DirichletBC> > & get_bc (void) const { return bcu; } - void add_bc (const boost::shared_ptr <const dolfin::FunctionSpace> & V, boost::shared_ptr<const dolfin::GenericFunction> f, std::size_t n) {