Mercurial > fem-fenics-eugenio

changeset 42:741a7999f590

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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;
+}