Added Neumann BC operator with ability to chose which side it is impl…

…emented on

examples_MATLAB/helmholtz2D_wifi.m

@@ -0,0 +1,78 @@
+% Solves the 2D Helmholtz equation
+% lap(E) + k^2/n^2 E = 0
+% k = angular wave number
+% n = refractive index
+% This equation models the wifi signal propagation in a room
+
+clc
+close all
+
+addpath('../mole_MATLAB')
+
+% define the walls
+wall = @(X,Y) (X>=10.0 & X<=39.0 & Y>=20.0 & Y<=21.0) | (X>=30.0 & X<=31.0 & Y>=1.0 & Y<=16.0) | ...
+   (X<=0.5 | X>=39.5 | Y<=0.5 | Y>=39.5);
+
+% wall reflection coefficient
+aa = 0.7
+% wall absorption coefficient
+bb = 0.9*0.5
+% angular wave number
+wn = 6;
+
+%hotspot position
+hsx = 2.0
+hsy = 10.0
+hsr = 1.0
+
+% Spatial discretization
+k = 2;   % Order of accuracy
+m = 500;  % Number of cells along the x-axis
+n = 500;   % Number of cells along the y-axis
+a = 0;   % West
+b = 40;   % East
+c = 0;   % South
+d = 40;   % North
+dx = (b-a)/m;  % Step length along the x-axis
+dy = (d-c)/n;  % Step length along the y-axis
+
+% 2D Staggered grid
+xgrid = [a a+dx/2 : dx : b-dx/2 b];
+ygrid = [c c+dy/2 : dy : d-dy/2 d];
+
+[X, Y] = meshgrid(xgrid, ygrid);
+
+% complex coefficient c=k^2/n^2
+c = (wn^2./(1+(aa+i*bb)*wall(X,Y)).^2)';
+c = reshape(c, (m+2)*(n+2), 1);
+
+%to detect the hotspot's position
+HS = ((X-hsx).^2 + (Y-hsy).^2 < hsr^2)';
+HS = reshape(HS, (m+2)*(n+2), 1);
+ind = find(HS>0);
+freenodes = setdiff(1:(m+2)*(n+2),ind);
+
+% Mimetic operator (Laplacian)
+L = lap2D(k, m, dx, n, dy);
+id_op = diag(sparse(c));
+L = L + id_op;
+L = L + robinBC2D(k, m, dx, n, dy, 0, 1); % Neumann BC
+
+% RHS
+RHS = zeros(m+2,n+2);
+% aux = -2*(X.^2+Y.^2-X-Y) + (X.^2-X).*(Y.^2-Y);
+% RHS(2:end-1,2:end-1) = aux(2:end-1,2:end-1)';
+
+RHS = reshape(RHS, (m+2)*(n+2), 1);
+RHS = RHS - L*HS;
+
+SOL = zeros((m+2)*(n+2),1);
+SOL(freenodes) = L(freenodes,freenodes)\RHS(freenodes);
+SOL(ind) = ones(length(ind),1);
+
+% graph the logarithm of the modulus of SOL
+logSOL = log(abs(SOL));
+surf(X, Y, reshape(logSOL, m+2, n+2)');
+shading interp
+colorbar
+view(2)

mole_MATLAB/NeumannSidedBC.m

@@ -0,0 +1,33 @@
+function BC = NeumannSidedBC(k, m, dx, left, left_coef, right, right_coef)
+% Returns a m+2 by m+2 one-dimensional mimetic boundary operator that 
+% imposes a boundary condition of Neumann  type
+%
+% Parameters:
+%                k : Order of accuracy
+%                m : Number of cells
+%               dx : Step size
+%              left: T/F (binary) build left BC
+%        left_coef : Value for left BC
+%             right: T/F (binary) build right BC
+%       right_coef : Value for right BC
+
+    % Dirichlet could be made in a similar fashion
+    %A = sparse(m+2, m+2);
+
+    B = sparse(m+2, m+1);
+
+    if ( left )
+        B(1, 1) = -left_coef;
+    end
+
+    if ( right )
+        B(end, end) = right_coef;
+    end
+
+    G = grad(k, m, dx);
+
+    BC = B*G;
+
+    % Or if making Robin/Dirichlet
+    %BC = A + B*G;
+end

mole_MATLAB/NeumannSidedBC2D.m

@@ -0,0 +1,34 @@
+function BC = NeumannSidedBC2D(k, m, dx, n, dy, left, left_coef, right, right_coef, front, front_coef, back, back_coef)
+% Returns a two-dimensional mimetic boundary operator that 
+% imposes a boundary condition of Robin's type
+%
+% Parameters:
+%                k : Order of accuracy
+%                m : Number of cells along x-axis
+%               dx : Step size along x-axis
+%                n : Number of cells along y-axis
+%               dy : Step size along y-axis
+%              left: T/F (binary) build left BC
+%        left_coef : Value for left BC
+%            right : T/F (binary) build right BC
+%       right_coef : Value for right BC
+%            front : T/F (binary) build front BC
+%       front_coef : Value for front(bottom) BC
+%              back: T/F (binary) build back BC
+%        back_coef : Value for back(top) BC
+
+    % 1-D boundary operator
+    Bm = NeumannSidedBC(k, m, dx, left, left_coef, right, right_coef);
+    Bn = NeumannSidedBC(k, n, dy, front, front_coef, back, back_coef);
+
+    Im = speye(m+2);
+    In = speye(n+2);
+
+    In(1, 1) = 0;
+    In(end, end) = 0;
+
+    BC1 = kron(In, Bm);
+    BC2 = kron(Bn, Im);
+
+    BC = BC1 + BC2;
+end

mole_MATLAB/NeumannSidedBC3D.m

@@ -0,0 +1,53 @@
+function BC = NeumannSidedBC3D(k, m, dx, left, left_coef, ...
+                                        right, right_coef, ...
+                                        front, front_coef, ...
+                                         back, back_coef, ...
+                                          top, top_coef, ...
+                                       bottom, bottom_coef)
+% Returns a three-dimensional mimetic boundary operator that 
+% imposes a boundary condition of Robin's type
+%
+% Parameters:
+%                k : Order of accuracy
+%                m : Number of cells along x-axis
+%               dx : Step size along x-axis
+%                n : Number of cells along y-axis
+%               dy : Step size along y-axis
+%                o : Number of cells along z-axis
+%               dz : Step size along z-axis
+%              left: T/F (binary) build left BC
+%        left_coef : Value for left BC
+%             right: T/F (binary) build right BC
+%       right_coef : Value for right BC
+%             front: T/F (binary) build front BC
+%       front_coef : Value for front BC
+%             back : T/F (binary) build back BC
+%        back_coef : Value for back BC
+%              top : T/F (binary) build top BC
+%         top_coef : Value for top BC
+%           bottom : T/F (binary) build bottom BC
+%      bottom_coef : Value for bottom BC
+
+
+    % 1-D boundary operator
+    Bm = NeumannSidedBC(k, m, dx, left, left_coef, right, right_coef);
+    Bn = NeumannSidedBC(k, n, dy, front, front_coef, back, back_coef);
+    Bo = NeumannSidedBC(k, o, dz, top, top_coef, bottom, bottom_coef);
+
+    Im = speye(m+2);
+    In = speye(n+2);
+    Io = speye(o+2);
+
+    Io(1, 1) = 0;
+    Io(end, end) = 0;
+
+    In2 = In;
+    In2(1, 1) = 0;
+    In2(end, end) = 0;
+
+    BC1 = kron(kron(Io, In2), Bm);
+    BC2 = kron(kron(Io, Bn), Im);
+    BC3 = kron(kron(Bo, In), Im);
+
+    BC = BC1 + BC2 + BC3;
+end

mole_MATLAB/robinBC2D.m

@@ -16,8 +16,8 @@
     Bn = robinBC(k, n, dy, a, b);
 
     Im = speye(m+2);
-
     In = speye(n+2);
+
     In(1, 1) = 0;
     In(end, end) = 0;