1## Copyright (C) 2006-2013 Carlo de Falco, Massimiliano Culpo 2## 3## This file is part of: 4## BIM - Diffusion Advection Reaction PDE Solver 5## 6## BIM is free software; you can redistribute it and/or modify 7## it under the terms of the GNU General Public License as published by 8## the Free Software Foundation; either version 2 of the License, or 9## (at your option) any later version. 10## 11## BIM is distributed in the hope that it will be useful, 12## but WITHOUT ANY WARRANTY; without even the implied warranty of 13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14## GNU General Public License for more details. 15## 16## You should have received a copy of the GNU General Public License 17## along with BIM; If not, see <http://www.gnu.org/licenses/>. 18## 19## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> 20## author: Matteo Porro <meoo85 _AT_ users.sourceforge.net> 21 22## -*- texinfo -*- 23## 24## @deftypefn {Function File} {[@var{norm_u}]} = @ 25## bim2c_norm(@var{mesh},@var{u},@var{norm_type}) 26## 27## Compute the @var{norm_type}-norm of function @var{u} on the domain described 28## by the triangular grid @var{mesh}. 29## 30## The input function @var{u} can be either a piecewise linear conforming scalar 31## function or an elementwise constant scalar or vector function. 32## 33## The string parameter @var{norm_type} can be one among 'L2', 'H1' and 'inf'. 34## 35## Should the input function be piecewise constant, the H1 norm will not be 36## computed and the function will return an error message. 37## 38## For the numerical integration of the L2 norm the second order middle point 39## quadrature rule is used. 40## 41## @seealso{bim1c_norm, bim3c_norm} 42## 43## @end deftypefn 44 45function [norm_u] = bim2c_norm (m, u, norm_type) 46 47 ## Check input 48 if (nargin != 3) 49 error ("bim2c_norm: wrong number of input parameters."); 50 elseif (! (isstruct (m) && isfield (m,"p") 51 && isfield (m, "t") 52 && isfield (m, "e"))) 53 error ("bim2c_norm: first input is not a valid mesh structure."); 54 endif 55 56 nnodes = columns (m.p); 57 nel = columns (m.t); 58 59 if (isequal (size (u), [2, nel])) 60 u = u'; 61 endif 62 63 if ((numel (u) != nnodes) && (rows (u) != nel)) 64 error ("bim2c_norm: numel(u) != nnodes and rows(u) != nel."); 65 endif 66 67 if (! (strcmp (norm_type,'L2') 68 || strcmp (norm_type,'inf') 69 || strcmp (norm_type,'H1'))) 70 error ("bim2c_norm: invalid norm type parameter."); 71 endif 72 73 if (strcmp (norm_type,'inf')) 74 norm_u = max (abs (u(:))); 75 else 76 if (numel (u) == nnodes) 77 78 M = __mass_matrix__ (m); 79 80 if (strcmp (norm_type, 'H1')) 81 A = bim2a_laplacian (m, 1, 1); 82 M += A; 83 endif 84 85 norm_u = sqrt(u' * M * u); 86 87 else 88 89 if (strcmp (norm_type, 'H1')) 90 error (["bim2c_norm: cannot compute the H1 norm ", ... 91 "of an elementwise constant function."]); 92 endif 93 94 norm_u = m.area' * (norm (u, 2, 'rows').^2); 95 norm_u = sqrt (norm_u); 96 97 endif 98 endif 99 100endfunction 101 102function M = __mass_matrix__ (mesh) 103 104 t = mesh.t; 105 nnodes = columns (mesh.p); 106 nelem = columns (t); 107 108 ## Local contributions 109 110 Mref = 1/12 * [2 1 1; 1 2 1; 1 1 2]; 111 area = reshape (mesh.area, 1, 1, nelem); 112 113 ## Computation 114 for inode = 1:3 115 for jnode = 1:3 116 ginode(inode,jnode,:) = t(inode,:); 117 gjnode(inode,jnode,:) = t(jnode,:); 118 endfor 119 endfor 120 Mloc = area .* Mref; 121 122 ## assemble global matrix 123 M = sparse (ginode(:), gjnode(:), Mloc(:), nnodes, nnodes); 124 125endfunction 126 127%!test 128%!shared L, V, x, y, m 129%! L = rand (1); V = rand (1); x = linspace (0,L,4); y = x; 130%! m = msh2m_structured_mesh (x,y,1,1:4); 131%! m.area = msh2m_geometrical_properties (m, 'area'); 132%! m.shg = msh2m_geometrical_properties (m, 'shg'); 133%! u = V * ones (columns(m.p),1); 134%! uinf = bim2c_norm (m, u, 'inf'); 135%! uL2 = bim2c_norm (m, u, 'L2'); 136%! uH1 = bim2c_norm (m, u, 'H1'); 137%! assert ([uinf, uL2, uH1], [V, V*L, V*L], 1e-12); 138%!test 139%! u = V * (m.p(1,:) + 2*m.p(2,:))'; 140%! uinf = bim2c_norm (m, u, 'inf'); 141%! uL2 = bim2c_norm (m, u, 'L2'); 142%! uH1 = bim2c_norm (m, u, 'H1'); 143%! assert ([uinf, uL2, uH1], 144%! [3*L*V, V*L^2*sqrt(8/3), V*sqrt(8/3*L^4 + 5*L^2)], 145%! 1e-12); 146%!test 147%! u = V * ones (columns(m.t),1); 148%! uinf = bim2c_norm (m, u, 'inf'); 149%! uL2 = bim2c_norm (m, u, 'L2'); 150%! assert ([uinf, uL2], [V, V*L], 1e-12); 151%!test 152%! u = V * ones (columns(m.t),1); 153%! uvect = [u, 2*u]; 154%! uinf = bim2c_norm (m, uvect, 'inf'); 155%! uL2 = bim2c_norm (m, uvect, 'L2'); 156%! assert ([uinf, uL2], [2*V, V*L*sqrt(5)], 1e-12); 157