Mercurial > fem-fenics-eugenio
view src/assemble.cc @ 268:61830a4f9ab9
Improve formatting
author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
---|---|
date | Thu, 14 Aug 2014 12:26:55 +0200 |
parents | 5e9b5bbdc56b |
children | f4d6ae912a08 |
line wrap: on
line source
/* Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net> 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 3 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" #include "femfenics_factory.h" #include "dolfin_compat.h" DEFUN_DLD (assemble, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Function File} {[@var{A}], [@var{x}(Optional)]} = \ assemble (@var{form_a}, @var{DirichletBC}) \n\ Construct the discretization of a Form and apply essential BC.\n\ The input arguments are\n\ @itemize @bullet\n\ @item @var{form_a} which is the form to assemble.\n\ It can be a form of rank 2 (BilinearForm or JacobianForm), \ a form of rank 1 (LinearForm or ResidualForm) or a form \ of rank 0 (Functional).\n\ @item @var{DirichletBC} represents the optional BC applied to \ the system. \n\ @end itemize \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\ @end itemize \n\ If boundary condition has to be applied to a vector for a nonlinear problem \ then it should be provide as 2nd argument and it will be given back \ as the second output argument. For an example of this situation, please refer \ to the HyperElasticity example. \n\ @seealso{BilinearForm, LinearForm, ResidualForm, JacobianForm, Functional}\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 (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) { femfenics_factory factory; dolfin::Matrix A; dolfin::assemble (A, a); 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<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"); } } retval(0) = factory.matrix (A); } else if (a.rank () == 1) { femfenics_factory factory; dolfin::Vector A; dolfin::assemble (A, a); 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<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"); } } retval(0) = factory.vector (A); } else if (a.rank () == 0) { double b = dolfin::assemble (a); retval(0) = octave_value (b); } else { error ("assemble: unknown form size"); } } } } } else if (nargout == 2) { if (nargin < 2) { print_usage (); } else { if (args(0).type_id () == form::static_type_id ()) { const form & frm = static_cast<const form &> (args(0).get_rep ()); const Array<double> myx = args(1).array_value (); if (! error_state) { const dolfin::Form & a = frm.get_form (); a.check (); if (a.rank () == 1) { femfenics_factory factory; dolfin::Vector A; dolfin::assemble (A, a); #ifdef LATEST_DOLFIN dolfin::Vector x (MPI_COMM_WORLD, myx.length ()); #else dolfin::Vector x (myx.length ()); #endif for (std::size_t i = 0; i < myx.length (); ++i) { x.setitem (i, myx.xelem (i)); } 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<SHARED_PTR <const dolfin::DirichletBC> > & pbc = bc.get_bc (); for (std::size_t j = 0; j < pbc.size (); ++j) { pbc[j]->apply (A, x); } } else { error ("assemble: unknown argument type"); } } retval(0) = factory.vector (A); retval(1) = factory.vector (x); } else { error ("assemble: unknown size"); } } } } } return retval; }