function solveBV
% Codes showing how to solve a boundary value problem by command lines
%
% 1. Construct geometry by a geometrical matrix
% 2. Define boundary conditions (BCs) by a BC matrix
% 3. Match coefficients between the MATLAB elliptic solver and the PDE to
% be solved
% 4. Discretize PDE and BCs in the geometrical domain
% 5. Solve the BC problem and do post-processing
%
% These code need MATLAB PDE Toolbox
%
% By Hongxue Cai (h-cai@northwestern.edu)
%
% =============== Problem to be solved ===============>
%
% An incompressible fluid is in the domain deifned by the geometrical
% matrix gGeo below. This 2D domain is a cross section of a fluid chamber
% (scala tympani)
% in the cochlea, a spiral fluid-solid interacted structure in the inner ear.
% A deformable boundary (geometrical element #2 defined by the column# 2 in gGeo)
% represents a vibrating membrane (the basilar membrane). We want to find
% the distribution of the fluid pressure when a motion of the memebrane is
% given
%
% Equation to be solved: L(P) = k^2* P (1)
% where P: fluid pressure
% L: represent the Laplace operator
% k: the wave number of the traveling wave on the memebrane
%
% Note: This equation is derived from the Navier-Stokes equation and the
% incomressibility of the fluid. In fact, the P is the dominant term of the
% WKB method applied to the pressure (see references)
%
% Bounary conditions:
% dP/dn = r*w^2*U on deformable boundary (2)
% dP/dn = 0 on rigid walls (3)
% where: r: the intensity of the fluid
% n: the normal direction on a boundary point
% w = 2*Pi*freq, with freq, the vibrating frequency
% U: normal displacement component in the direction of n
%
% MATLAB solver used: -grad(c*grad(P)) + a*P = f, with the following genral
% Neuman condtion: n.(c*grad(P)) + q*P = g
% where c, a, f, q, g are coefficients to be defined by matching the PDE
% equatioin and BCs with the MATLAB solver, and we have
% c = -1
% a = k^2
% f = 0
% q = 0
% g = 0 or -r*w^2*U
%
%References:
% 1. Cai, H., Manoussaki, D. and Chadwick, R.S., Effects of coiling on the micromechanics
% of the mammalian cochlea, J. R. Soc. Interface, 2(4):341-348, 2005.
% 2. Cai, H., Shoelson, B and Chadwick R.S., Evidence of tectorial membrane radial motion
% in a propagating mode of a complex cochlear model,
% Proc. Natl. Acad. Sci. USA, 101(16):6243-6248, 2004.
% 3. Cai, H. and Chadwick, R.S., Radial structure of traveling waves in the inner ear,
% SIAM J. Appl. Math., 63(4):1105-1120, 2003.
%
% < ======================================================
% Geometrical metrix -------------------------------------------->
BM1 = -0.0079;
BM2 = 0.0171;
R0 = 0.022;
xc = 0;
yc = 0;
yBM = -R0/2.5;
xR = -sqrt(R0*R0 - yBM*yBM);
xL = -xR;
gGeo = [
2 2 2 1
xR BM1 BM2 xR
BM1 BM2 xL xL
yBM yBM yBM yBM
yBM yBM yBM yBM
0 0 0 1
1 1 1 0
0 0 0 xc
0 0 0 yc
0 0 0 R0
];
%geometry matrix, in which, each column defines a geometrical element that
%can be either a line or a circular arc. (Note that any complex geometry can
%be composed by line and circular arts.) In each column, the number in row#1
%indicates the type of the geometrica element, 2 for line element, 1 for
%circular arts. Then, row#2 is the x-coordinate of the starting point;
%row#3, that of the ending point; Row #4 and #5 are, respectively the
%y-coordinates of the the starting and ending point. The number in the 6th
%anf 7th
%rows indicate the location of geometric domain defined, (0,1) on the left side of this
%geometric element (When you walk from the begining point to the ending point along the element,
%the domain is on your left), (1,0) on th right. The rows 8, 9, 10 are zero for line,
%and repesent the x, y-coordinates and the radius of the center point for a
%circular arc.
%
%The values that define the geometry come from an application which models
%the baislar
%membrane in the inner ear. These values should be different for different
%applications, therefore one should not pay much attention to the
%individual numbers themselves.
% plotting geometry
figure
pdegplot(gGeo) %MATLAB PDE Toolbox needed
%meshing
[p,e,t]=initmesh(gGeo); %initial mesh
[p,e,t]=refinemesh(gGeo,p,e,t); % refine mesh
[p,e,t]=refinemesh(gGeo,p,e,t);
%meshing data are stored in p, e, t
% display mesh
figure
pdemesh(p,e,t); %PDE Toolbox needed
% < -------------------------------------------
% ----------------- Define constant or known variable -------------->
freq = 100; %frequency, Hz
w = 2*pi*freq;
r = 1; %fluid intensity, g/cm^3
rw2 = r*w*w;
%
k = 6 - 0.03i; %complex wavenumber
%Note that the geometrical domain (the cross section) is in the x-y plane, and the wave
%propogates in the z- direction with a wavenumber k. The geometry of the
%cochlea chages slowly in the z-direction, and the WKB method allows us to
%solve approximately a 3D flow in a 2D doamin (See references)
%
BMa = pi*(BM1+BM2)/2/(BM1-BM2); %no initial mode change
BMb = -pi/(BM1-BM2); %no initial mode change
Am = 1e-9; %Amplitude, 1 nm. Initially, the defomable element (the basilar memebrane) is given a motion
% y = Am*cos(BMa + BMb*x), a clamped beam with maximum
% displacement at the middle. When x = BM1, BM2, y = 0 (BMa +
% BMb*x = -pi/2, pi/2)
% < --------------------------------------------------------------------
% define PDE coefficients --------------------->
c = -1;
a = k^2;
f = 0;
% < -------------------------------------------
% Define boudary condition ------------------------------->
%
% In this Ex, the geometrical element #2 (define by the column#2 in the
% matrix gGeo, subject to a non-zero Neumann BC; all other boundaries (#1,
% 3, and 4) are rigid walls, therefore a zero Neumann BC applied to.
%
% How many geometrical elements
n_Geo = size(gGeo,2);
n_mb = 2; %No. of deformable boundary in gGeo
q = 48; % q = 0
% g 4 deformable boundary
clear g
temp = [mat2str(-rw2*Am),'*cos(x*',mat2str(BMb),'+',mat2str(BMa),')'];
g(1:length(temp),1)=(temp*1)';
% to generate special codes in the BC matrix
% boundary condition matrix
bc = bcNeumann(n_mb, q, g, n_Geo);
%solve the PDE
disp('Elliptic problem (fluid in the lower chamber....>')
press = assempde(bc, p, e, t, c, a, f);
% post-processing ...........
%plotting the pressure distribution
figure
pdeplot(p, e, t, 'xydata', real(press))
% ============ --------- end of sloveBV -------------========
%
%
%sub-function
function bc = bcNeumann(n_mb, q, g, nCol)
% -------- Neumann boundary condition matrix (scalar case) -------
% n_mb: Non-zero Neumann boundary segments (moving segments) index
% nCol: Total number of boundary segments, nCol-length(n_mb)=
% zero Neumann boundary segments number
%
% ------- n.(c grad u)+ q u = g -----------------
%
Lmb = length(n_mb); % how many deformable boundaries
%
Lq = length(q);
Lg = length(g);
%
% number of rows in boundary condition matrix bc
row_b = max(Lq + Lg);
row_b = row_b + 4;
%
% set default values to 0
bc = zeros(row_b, nCol);
%
bc(1,:)=1; %dimension of the system N = 1 (scalar case)
bc(2,:)=0; %Number of Dirichlet boundary conditions M = 0
bc(3,:)=1; %Lq = 1 for zero Neumann boundary segments (q = 0)
bc(4,:)=1; %Lg = 1 for zero Neumann boundary segments (g = 0)
%
%Lq for non-zero Neumann boundary segments
for ij=1:Lmb
bc(3,n_mb(ij)) = Lq(ij);
end
%Lg for non-zero Neumann boundary segments
for ij=1:Lmb
bc(4,n_mb(ij)) = Lg(ij);
end
bc(5:6,:)=48; %q = g = 0 (48) for zero Neumann boundary segments
%q for non-zero Neumann boundary segments
nRow(1:Lmb) = 4;
for ij=1:Lmb
bc(nRow(ij)+1:nRow(ij)+Lq(ij),n_mb(ij)) = q(1:Lq(ij),ij);
end
nRow=nRow+Lq;
%g for non-zero Neumann boundary segments
for ij=1:Lmb
bc(nRow(ij) + 1:nRow(ij) + Lg(ij), n_mb(ij)) = g(1:Lg(ij),ij);
end
nRow = nRow + Lg;
row_b = length(bc);
%
return
