diff --git a/examples_MATLAB/helmholtz2D_wifi.m b/examples_MATLAB/helmholtz2D_wifi.m new file mode 100644 index 0000000..21aaa0c --- /dev/null +++ b/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) diff --git a/mole_MATLAB/NeumannSidedBC.m b/mole_MATLAB/NeumannSidedBC.m new file mode 100644 index 0000000..b5bebdf --- /dev/null +++ b/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 diff --git a/mole_MATLAB/NeumannSidedBC2D.m b/mole_MATLAB/NeumannSidedBC2D.m new file mode 100644 index 0000000..d72f5ff --- /dev/null +++ b/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 diff --git a/mole_MATLAB/NeumannSidedBC3D.m b/mole_MATLAB/NeumannSidedBC3D.m new file mode 100644 index 0000000..fc99986 --- /dev/null +++ b/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 diff --git a/mole_MATLAB/robinBC2D.m b/mole_MATLAB/robinBC2D.m index 86b9420..e948110 100644 --- a/mole_MATLAB/robinBC2D.m +++ b/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;