/usr/share/octave/packages/bim-1.1.5/bim3a_osc_laplacian.m is in octave-bim 1.1.5-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 | ## Copyright (C) 2012 Carlo de Falco
##
## This file is part of:
## BIM - Diffusion Advection Reaction PDE Solver
##
## BIM 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.
##
## BIM 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 BIM; If not, see <http://www.gnu.org/licenses/>.
##
## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
## -*- texinfo -*-
## @deftypefn {Function File} @
## {@var{A}} = bim3a_osc_laplacian (@var{mesh}, @var{epsilon})
##
## Build the osc finite element stiffness matrix for a diffusion
## problem.
##
## For details on the Orthogonal Subdomain Collocation (OSC) method
## see: M.Putti and C.Cordes, SIAM J.SCI.COMPUT. Vol.19(4), pp.1154-1168, 1998.
##
## The equation taken into account is:
##
## - div (@var{epsilon} grad (u)) = f
##
## where @var{epsilon} is an element-wise constant scalar function.
##
## @seealso{bim3a_rhs, bim3a_reaction, bim2a_laplacian, bim3a_laplacian}
## @end deftypefn
function M = bim3a_osc_laplacian (msh, epsilon)
## Check input
if (nargin != 2)
print_usage ();
elseif (! (isstruct (msh)
&& isfield (msh, "p")
&& isfield (msh, "t")
&& isfield (msh, "e")))
error (["bim3a_osc_laplacian: first input ", ...
"is not a valid msh structure"]);
endif
nnodes = columns (msh.p);
nelem = columns (msh.t);
## Turn scalar input to a vector of appropriate size
if (isscalar (epsilon))
epsilon = epsilon * ones (nelem, 1);
endif
if (! isvector (epsilon))
error ("bim3a_osc_laplacian: coefficient is not a vector");
elseif (numel (epsilon) != nelem)
error (["bim3a_osc_laplacian: length of epsilon is ", ...
"not equal to the number of mesh elements"]);
endif
## Avoid warnings for broadcasting
warning ("off", "Octave:broadcast", "local")
## Local contributions
Lloc = __osc_local_laplacian__ (msh.p, msh.t, msh.shg,
epsilon, msh.area, nnodes,
nelem);
## Assembly
for inode = 1:4
for jnode = 1:4
ginode(inode, jnode,:) = msh.t(inode, :);
gjnode(inode, jnode,:) = msh.t(jnode, :);
endfor
endfor
M = sparse (ginode(:), gjnode(:), Lloc(:), nnodes, nnodes);
endfunction
%!shared msh, epsilon, M, nnodes, nelem, x, y, z
%!test
%! msh = bim3c_mesh_properties (msh3m_structured_mesh (0:5, 0:5, 0:5, 1, 1:6));
%! x = msh.p (1, :).';
%! y = msh.p (2, :).';
%! z = msh.p (3, :).';
%! u = ones (size (x));
%! M = bim3a_osc_laplacian (msh, 1);
%! assert (M * u, zeros (size (u)), eps * 100)
%!test
%! u = x;
%! bnd = bim3c_unknowns_on_faces (msh, [1, 2]);
%! int = setdiff (1:columns (msh.p), bnd);
%! assert (M(int, int) * u(int), -M(int, bnd) * u(bnd), 100 * eps)
%!test
%! u = y;
%! bnd = bim3c_unknowns_on_faces (msh, [3, 4]);
%! int = setdiff (1:columns (msh.p), bnd);
%! assert (M(int, int) * u(int), -M(int, bnd) * u(bnd), 100 * eps)
%!test
%! u = z;
%! bnd = bim3c_unknowns_on_faces (msh, [5, 6]);
%! int = setdiff (1:columns (msh.p), bnd);
%! assert (M(int, int) * u(int), -M(int, bnd) * u(bnd), 100 * eps)
%!test
%! u = z;
%! bnd = bim3c_unknowns_on_faces (msh, [5, 6]);
%! int = setdiff (1:columns (msh.p), bnd);
%! M = bim3a_osc_laplacian (msh, pi);
%! assert (M(int, int) * u(int), -M(int, bnd) * u(bnd), 100 * eps)
|