No BSD License  

Different material properties in different sub-domains

by

 

21 Feb 2006 (Updated )

Codes showing how to use different material properties in different sub-domains ...

subDemo
function subDemo
%  ========   Codes showing how to use different material properties in different sub-domains
%           in mechanical modeling
%  ========   This is also a good example to show how to use the boundary
%           condition matrix ("bcMatrix" in a recent submission) in solving
%           boundary value problems
%
% These codes need MATLAB PDE Toolbox
%
%       by Hongxue Cai (h-cai@northwestern.edu)
%
% ===============  Problem to be solved   ===============>
%
% The organ of corti (OC) vibrates and interacts with cochlear fluid. It is composed of different cells, 
% e.g. inner and  outer hair cells (IHCs and OHCs), Pillar cells, ... 
% Different cells have different material properties (Young's modulus). 
%
% Equation to be solved:     - div (X) + r d2U/dt2  = 0             
% where     X:  stress tensor
%           U:  displacement vector having u, and v as its x and y
%               componebts
%           t:  time
%           r:  density
% Expressing X by u and v, we can arrive at a PDE system of order 2 for u
% and v.
%
% Bounary conditions:
%       1. Xn = P n + T s = [g1 g2]'       
% for deformable boundaries in contact with fluid.
% where:    Xn:     normal stress vector
%           n:      the normal direction on a boundary point
%           s:      the tangential direction on a boundary point
%           P:      fluid pressure
%           T:      tangetial stress (see sub-function g1g2 below)
%
%        2. u = v = 0, on rigid walls
%
%        3. u = 0, v = Y_BM for the boundaries in contact with the BM
%        (basilar memebrane)
% where: Y_BM is the deflection of the BM (known)
%
% We want to solve the displacement field in the OC when the surounding
% fluid distribution is known.
%
% MATLAB solver used: elliptic system solver
% 
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%												%
%  FEM: Geomtry -> mesh -> material properties -> boundary condition ->  %
%       solving  -> post processing 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 
%
clc
disp('Plane strain for the organ of Corti (OC)')
%
% % -------     Geometry    --------->
gOC = geoMatrix;        %see sub-function geoMatrix below
%
% % Display geometry
figure
pdegplot(gOC)
% % <-------     Geometry    --------
%
%
% % 			--------  meshing ------------->
disp('Meshing ...')
[p,e,t] = initmesh(gOC);
[p,e,t] = refinemesh(gOC,p,e,t);
%
Ntri    = size(t,2);	% number of triangles in OC mesh
%
% %       --- Display Mesh --->
figure
pdemesh(p,e,t);
% 			<--------  meshing -------------
%
% Known parameters 
freq    =  1000;        %Hz, vibratory frequency
omg     = 2*pi*freq;    % radian frequency
k       = 39.5-5.5i;    % wavenumber
Am      = 1e-9;         % Initial Y_BM = Am*cos(a+b*x), (-pi/2 ~ pi/2)
%
%	--------- Meterial properties  -------------->
nu    = 0.49;         % Poisson's ratio of OC
%
rou_OC  = 1;            % g/cm3     density of OC
rou     = 1;			% g/cm3     density of cochlear fluid
v_fluid = 0.01;         % cm2/s		fluid viscosity

cai     = rou*omg*omg;
tao     = sqrt(i*omg*v_fluid)/i/omg;
%
E_OC    = 2e3;          % Young's modulus of OC: dyn/cm2
E_ISC   = E_OC;         % Young's modulus for the inner sulcus cell
E_HS    = 1e4;          % Young's modulus of Hensen cell 2e3

E_OHC   = 6e4;          % Young's modulus of OHCs 
E_IHC   = 4e4;          % Young's modulus of IHCs
E_De    = 1e5;          % Young's modulus of Deiter's celle
E_RL    = 3e5;          % Young's modulus of RL (reticulat lamina)
E_Pi    = 4e5;          % Young's modulus of pillar cells
%
% To handle damping, we model the OC as a "Voigt solid." i.e. we use a
% complex Young's modulus such that E  =  E + i E'
coeE = 0.1;
E_ISC     = E_ISC + i*coeE*E_ISC;
E_OC      = E_OC + i*coeE*E_OC;
E_OC      = E_OC + i*coeE*E_OC;
E_Pi      = E_Pi + i*coeE*E_Pi;
E_OHC     = E_OHC + i*coeE*E_OHC;
E_IHC     = E_IHC + i*coeE*E_IHC;
E_De      = E_De + i*coeE*E_De;
E_RL      = E_RL + i*coeE*E_RL;
E_HS      = E_HS + i*coeE*E_HS;

% Discretizing material properties in OC subdomains ...
% e.g. subdomain no. 2 is Pillar cells, so E(pdesdt(t,2)) = E_Pi 
E = zeros(1,Ntri);
E(pdesdt(t,1))  = E_OC ;

E(pdesdt(t,3))  = E_OC ;    E(pdesdt(t,10)) = E_OC ;
E(pdesdt(t,11)) = E_OC ;    E(pdesdt(t,25)) = E_OC ;

E(pdesdt(t,12)) = E_HS ;
E(pdesdt(t,2))  = E_Pi ;
E(pdesdt(t,4))  = E_OHC ;   E(pdesdt(t,5)) = E_OHC ;    E(pdesdt(t,6)) = E_OHC ;
E(pdesdt(t,7))  = E_De ;    E(pdesdt(t,8)) = E_De ;     E(pdesdt(t,9)) = E_De ;
E(pdesdt(t,13)) = E_IHC ;
E(pdesdt(t,15)) = E_RL ;	E(pdesdt(t,17)) = E_RL ;	E(pdesdt(t,19)) = E_RL ;
E(pdesdt(t,21)) = E_RL ;	E(pdesdt(t,23)) = E_RL ;    E(pdesdt(t,24)) = E_RL ;
E(pdesdt(t,14)) = E_RL ;
E(pdesdt(t,16)) = E_RL ;
E(pdesdt(t,18)) = E_RL ;
E(pdesdt(t,20)) = E_RL ;
E(pdesdt(t,22)) = E_RL ;

E(pdesdt(t,26)) = E_ISC ;E(pdesdt(t,27)) = E_ISC ;E(pdesdt(t,28)) = E_ISC ;
E(pdesdt(t,29)) = E_ISC ;E(pdesdt(t,30)) = E_ISC ;
%
% % display E values (different in different sub-domains) -- log plot
figure
pdeplot(p, e, t, 'xydata', log10(E));
%
%PDE coefficient for MATLAB elliptic system solver: -div(c u) + au  =  f
% %f
xForce  = 0.0;          % x body force
yForce  = 0.0;          % y body force
f       = [
    xForce
    yForce
]; 
%
% %a
% 4 axial coupling
B       = E/2/(1+nu);
ak      = B*k*k;
a       = zeros(4,Ntri);
a(1,:)  = ak-rou_OC*omg*omg;	%+i*C_OC*omg;
a(4,:)  = ak-rou_OC*omg*omg;
%
% %c
c       = zeros(10,Ntri);
c(1,:)  = E/(1+nu)+E/(1+nu)*nu/(1-2*nu);
c(10,:) = E/(1+nu)+E/(1+nu)*nu/(1-2*nu);
c(3,:)  = E/(2*(1+nu));
c(5,:)  = E/(2*(1+nu));
c(8,:)  = E/(2*(1+nu));
c(6,:)  = E/(1+nu)*nu/(1-2*nu);
%
%
%   ================= Boundary conditions   =======================>
%
%  -------  No-zero Neumann segments  ------->
in_N = [3 50:59 69 31 79 80 83 86 89 45 70 49 73];
%
%coefficients for pressure p = polyval(c_p, x);
c_p = [
 -2.64972744951575 - 0.57062166587713i
 -0.15416754647399 - 0.03318840490841i
 -0.00358712479416 - 0.00077243509093i
 -0.00004171517817 - 0.00000897860636i
 -0.00000024276194 - 0.00000005228013i
 -0.00000000056410 - 0.00000000012134i
 ]*1e7;
%
% coefficients for pressure gradient in x direction
c_dpdx = [
 -5.72231461839305 + 0.38821140575060i
 -0.33304900905839 + 0.02177104182281i
 -0.00775245003761 + 0.00048778370550i
 -0.00009020197081 + 0.00000545957122i
 -0.00000052471448 + 0.00000003050022i
 -0.00000000122046 + 0.00000000006815i
 ]*1e9;
%
% coefficients for pressure gradient in y direction
c_dpdy = [
  3.33679411558593 + 1.75748861745406i
  0.19235170038884 + 0.10160526109744i
  0.00443394045499 + 0.00234934946743i
  0.00005106948251 + 0.00002714246053i
  0.00000029411536 + 0.00000015681108i
  0.00000000067665 + 0.00000000036194i
]*1e9;

[g1OC, g2OC] = g1g2(in_N, c_p, c_dpdx, c_dpdy, tao);

%   <---- gap segment needs special treatment

%  ----	Zero Dirichlet segments: rigd walls -------->
in_D0   = [1 8 30 76 77 81 84 87 90];

% -----  Non-zero Dirichlet segments (BM in gOC)-------->
% for geometrical elements in contact with the basilar membrane (BM)
in_D    = [11 75 25:29];
BM1     = 0.0125/5.4;
BM2     = 0.0125 + (BM1 + 0.0125);
BMa     = pi*(BM1 + BM2)/2/(BM1 - BM2);         %no initial mode change
BMb     = -pi/(BM1 - BM2);     
temp    = [mat2str(Am),'*cos(x*',mat2str(BMb),'+',mat2str(BMa),')'];

clear r2OC Lr2OC r1OC Lr1OC
for ij = 1:length(in_D)
    r2OC(1:length(temp),ij) = (temp)';
    r1OC(1,ij) = '0';
end
% <-----  Non-zero Dirichlet segments (BM in gOC)--------

clear qN gN hD rD
for mj  =  1:4
    for mk  =  1:length(in_N)
        qN{mj, mk}  =  '0';
    end
    for mk  =  1:length(in_D)
        hD{mj, mk}  =  '0';
    end
end
for mk  =  1:length(in_D)
    hD{1, mk}  =  '1';      %h11  =  1
    hD{4, mk}  =  '1';      %h22  =  1
    rD{1, mk}  =  r1OC(:, mk);
    rD{2, mk}  =  r2OC(:, mk);
end
for mk  =  1:length(in_N)
    gN{1, mk}  =  g1OC(:, mk);
    gN{2, mk}  =  g2OC(:, mk);
end

% OC subdomain segments:
in_nb = [2 4:7 9 10 12:24 32:44 46:48 60:68 71 72 74 78 82 85 88];

%No zero Neumman bc
in_N0 = 0;

%No mixed bc
in_M = 0;

%boundary condition matrix
bcm  =  bcMatrix(2, in_N, qN, gN, in_N0, in_D, hD, rD, in_D0, in_M, 0, 0, 0, 0, in_nb);

uv  = assempde(bcm, p, e, t, c, a, f); 	%uv = (u,v)T
u   = uv(1:length(uv)/2);
v   = uv(length(uv)/2+1:length(uv));

% %post processing   ------------>
%
% %plotting displacement
figure
pdeplot(p, e, t, 'xydata', real(v))
title('Displacement in the y dirextion (v)')
xlabel('x (cm)')
ylabel('y (cm)')
%
% %plotting stress
syy = pdesmech(p,t,c,uv,'tensor', 'syy');
syyNode=pdeprtni(p, t, syy);
figure
pdeplot(p, e, t, 'xydata', real(syyNode))
title('Stress in the y dirextion (Xyy)')
xlabel('x (cm)')
ylabel('y (cm)')


%	--------- end of main program  ---------------------------


%sub-functions  ==============>
function bcm  =  bcMatrix(N, in_N, qN, gN, in_N0, in_D, hD, rD, in_D0...
		,in_M, qM, gM, hM, rM, in_nb)
%	 ======   Providing boundary conditions matrix bcm   ========= 
%   =======  (system case, e.g. for 2D plane-strain problem, 2 PDEs for 2 dependent variables)  ===== 
%  
% A crutial step for solving boundary value problems (BVP) is to define the boundary conditions (BCs). 
% In a recent submission,
% "Solve boundary value probelm by command lines," an simple example is
% given for the scalar case (1 PDE for 1 dependent variable)
% subject to Neumann BC only (sub-function "bcNeumann" in that submission). 
% I now try to show how to construct the BC
% matrix for system case (Dimession N, i.e. N PDEs for N dependent variables). 
% Dirichlet, Neumann, and mixed BCs are considered.
% A more complex example showing how to use this function will be provided
% soon
%
%   ===== 	By Hongxue Cai (h-cai@northwestern.edu)   ======= 
%
%
% ---- Input:  ---------->
%       N:      dimnession of the PDE system (e.g. N = 1 for scalar case, N = 2
%               for 2D plane-strain problem)
%       in_N:   index of non-zero Neumann boundary segments (e.g. if the # 2,
%               6, and 7 geometrical elements are subject to non-zero Neumann BC,
%               then in_N  =  [2 6 7], or [2 7 6] ...;
%               in_N = 0 if no non-zero Neumman bc
%       qN:     q  for non-zero Neumann boundary segments, cell array of
%               size N^2 x length(in_N) (Note: cell array is convinient as
%               q are text exprssions of different length)
%               e.g. for 2D plane-strain problem, N = 2, qN{1, 1:length(in_N)} = q11, 
%               and if in_N = [2 6 7], then qN{1, 1} is q11 for boundary #2, 
%               qN{1, 2} is q11 for boundary #6; qN{1, 3} is q11 for boundary #7;
%               Similarly, 
%               qN{2, 1:length(in_N)}= q21; qN{3, 1:length(in_N)} = q12; and qN{4, 1:length(in_N)} = q22.
%               qN = 0 if no non-zero Neumman bc
%       gN:     g  for non-zero Neumann boundary segments, cell array of
%               size N x length(in_N)
%               gN = 0 if no non-zero Neumman bc
%
%       in_N0:  index of zero Neumann boundary segments
%               (g1 = g2 = ... = gN = 0)
%
%       in_D:   index of non-zero Dirichlet boundary segments
%       hD:     h  for non-zero Dirichlet boundary segments, cell array
%       rD:     r  for non-zero Dirichlet boundary segments, cell array
%
%       in_D0:  index of zero Dirichlet boundary conditions (r1 = r2  = ... = rN = 0)
%
%       in_M:	index of mixed boundary segments
%       qM:     q  for mixed boundary segments, cell array
%       gM:     g  for mixed boundary segments, cell array
%       hM:     h  for mixed boundary segments, cell array
%       rM:     r  for mixed boundary segments, cell array
%Note: in the case of mixed boundary condition, 1 Dirichlet bc is assumed
%(M =1)
%
%       in_nb:  index for non-boundary segements, thoses inside the
%               geometrical domain who seperate the sub-domains of the geometry
%       (an example showing how to use this function and how to use
%       different material properties in different sub-domains will be
%       given soon

% if non-cell input for q, g, h, r, convert to cell 
if isnumeric(qN)
    qN = num2cell(qN);
end
if isnumeric(gN)
    gN = num2cell(gN);
end
if isnumeric(hD)
    hD = num2cell(hD);
end
if isnumeric(rD)
    rD = num2cell(rD);
end
if isnumeric(qM)
    qM = num2cell(qM);
end   
if isnumeric(gM)
    gM = num2cell(gM);
end
if isnumeric(hM)
    hM = num2cell(hM);
end
if isnumeric(rM)
    rM = num2cell(rM);
end   

% length of q and g for non-zero Neumann boundaries
for k = 1:length(in_N)
    for j = 1:N^2
        LqN(j, k) = length(qN{j, k});
    end
    for j = 1:N
        LgN(j, k) = length(gN{j, k});
    end
    Lsumqg(k) = sum(LqN(:, k)) + sum(LgN(:, k));
end
LN = max(Lsumqg);       %maximum rows of text expression

% length of h and r for non-zero Dirichlet boundaries
for k = 1:length(in_D)
    for j = 1:N^2
        LhD(j, k) = length(hD{j, k});
    end
    for j = 1:N
        LrD(j, k) = length(rD{j, k});
    end
    Lsumhr(k) = sum(LhD(:, k)) + sum(LrD(:, k));
end
LD = max(Lsumhr);       %maximum rows of text expression


% length of q, g, h, and r for mixed boundaries
if in_M ~= 0
    for k = 1:length(in_M)
        for j = 1:N^2
            LqM(j, k) = length(qM{j, k});
%             LhM(j, k) = length(hM{j, k});
        end
        for j = 1:N    
            LgM(j, k) = length(gM{j, k});
            LhM(j, k) = length(hM{j, k});
        end
        LrM(1, k) = length(rM{1, k});       %for mixed bc, 1 Dirichlet bc assumed
        Lsumqghr(k) = sum(LqM(:, k)) + sum(LgM(:, k)) + sum(LhM(:, k)) + sum(LrM(:, k));
    end
    LM = max(Lsumqghr);       %maximum rows of text expression
else
    LM = 0;
end

L_N = length(in_N);
L_N0 = length(in_N0);
L_D = length(in_D);
L_D0 = length(in_D0);
L_M = length(in_M);
Lnb = length(in_nb);

%number of column of the bc matrix
coLnbcm = L_N + L_N0 + L_D + L_D0 + L_M + Lnb;

% row # before text expression begins
row0_N = 3 + N^2 + N + - 1;         % non-zero Neumann boundaries
M = N;
row0_D = 3 + M - 1;                 % non-zero Dirichlet boundaries   
M = 1;                              % 1 Dirichlet bc is assumed
row0_M = 3 + N^2 + N + M*N + M - 1; % mixed boundaries

% number of rows of the bc matrix
row_bcm = max([LN+row0_N, LD+row0_D, LM+row0_M]);

bcm = zeros(row_bcm, coLnbcm);

bcm(1,:) = N;	% row#1 --> Dimension of the system


%	------------------ Non-zero Neumman boundaries   ---------------------->
if in_N ~= 0
    for ij = 1:L_N
        bcm(2, in_N(ij)) = 0;          % M = 0 -> non Dirichlet bc       
        bcm(3:3+N^2-1,in_N(ij)) = LqN(:, ij);       %Length for q11, q21, ..., qNN
        bcm(3+N^2:3+N^2+N-1,in_N(ij)) = LgN(ij);	%Length for g1, g2, ..., gN
        %Note: as M = 0, no need for Lh, Lr, h, r
        %
        %q for non-zero Nmann boundary segments
        nRowN = 3 + N^2 + N;
        for j = 1:N^2
            bcm(nRowN:nRowN+LqN(j, ij)-1, in_N(ij)) = qN{j,ij};
            nRowN = nRowN + LqN(j, ij);
        end

        %g for non-zero Nmann boundary segments
        for j = 1:N
            bcm(nRowN:nRowN+LgN(j, ij)-1, in_N(ij)) = gN{j,ij};   
            nRowN = nRowN + LgN(j, ij);
        end
    end
end

%	<------------------ end of Non-zero Neumman boundaries  --------
%
%
%	------------------ Zero Neumman boundaries   ---------------------->
if in_N0 ~= 0
    bcm(2, in_N0) = 0;          % M = 0 -> non Dirichlet bc
    bcm(3:3+N^2-1,in_N0) = 1;       %Length for q11, q21, ..., qNN
    bcm(3+N^2:3+N^2+N-1,in_N0) = 1;	%Length for g1, g2, ..., gN
    %Note: as M = 0, no need for Lh, Lr, h, r
    %
    %q
    nRowN = 3 + N^2 + N;
    bcm(nRowN:nRowN+N^2-1, in_N0) = '0';
    nRowN = nRowN + N^2;

    %g for non-zero Nmann boundary segments
    bcm(nRowN:nRowN+N-1, in_N0) = '0';
    nRowN = nRowN + N;
end

%	<------------------ Zero Neumman boundaries  --------
%
%
%
%	-------------------     Non-zero Dirichlet boundaries  --------------->
if in_D ~= 0
    M = N;
    for ij = 1:L_D
        bcm(2, in_D(ij)) = M;                                          %M = N -> Dirichlet bc only
        bcm(3:3+N^2-1, in_D(ij)) = 1;                           %Length for q11, q21, ..., qNN
        bcm(3+N^2:3+N^2+N-1, in_D(ij)) = 1;                     %Length for g1, g2, ..., gN
        bcm(3+N^2+N:3+N^2+N+M*N-1, in_D(ij)) = LhD(:,ij);       %Length for h11, h21, ..., hNN
        bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_D(ij)) = LrD(:,ij); %Length for r1, r2, ..., rN
        %
        %q = 0
        nRowD = 3+N^2+N+M*N+M;
        for j = 1:N^2
            bcm(nRowD:nRowD+N^2-1,in_D(ij)) = '0'; 
        end
        nRowD = nRowD + N^2;
        
        %g = 0
        for j = 1:N
            bcm(nRowD:nRowD+N-1,in_D(ij)) = '0'; 
        end
        nRowD = nRowD + N;
        
        %h
        for j = 1:M^2
            bcm(nRowD:nRowD+LhD(j, ij)-1, in_D(ij)) = hD{j,ij};
            nRowD = nRowD + LhD(j, ij);
        end

        %r
        for j = 1:M
            bcm(nRowD:nRowD+LrD(j, ij)-1, in_D(ij)) = rD{j,ij};
            nRowD = nRowD + LrD(j, ij);
        end
    end
end
%	<-------------------  end of  Non-zero Dirichlet boundaries
%
%
%	-------  Zero Dirichlet boundaries   (rigid walls) ----->
if in_D0 ~= 0
    M = N;
    bcm(2, in_D0) = M;                                      %M = N -> Dirichlet bc only
    bcm(3:3+N^2-1, in_D0) = 1;                          %Length for q11, q21, ..., qNN
    bcm(3+N^2:3+N^2+N-1, in_D0) = 1;                    %Length for g1, g2, ..., gN
    bcm(3+N^2+N:3+N^2+N+M*N-1, in_D0) = 1;              %Length for h11, h21, ..., hNN
    bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_D0) = 1;        %Length for r1, r2, ..., rN
    %
    %qij = 0
    nRowD = 3+N^2+N+M*N+M;
    for j = 1:N^2
        bcm(nRowD:nRowD+N^2-1,in_D0) = '0';
    end
    nRowD = nRowD + N^2;

    %g = 0
    for j = 1:N
        bcm(nRowD:nRowD+N-1,in_D0) = '0';
    end
    nRowD = nRowD + N;

    %hii = 1, hij = 0
    bcm(nRowD:nRowD+M^2-1, in_D0) = '0';
    bcm(nRowD:M+1:nRowD+M^2-1, in_D0) = '1';
    nRowD = nRowD + M^2;

    %r = 0
    bcm(nRowD:nRowD+M-1, in_D0) = '0';
       
end
%	<-------  end of Zero-Dirichlet boundaries   (rigid walls) -----
%
%
%	------------------     Mixed boundaries  -------------------- ----->
if in_M ~= 0
    M =1;
    for ij = 1:L_M
        bcm(2, in_M(ij)) = M;                                   %M = 1 -> 1 Dirichlet bc
        bcm(3:3+N^2-1, in_M(ij)) = LqM(:,ij);;                  %Length for q11, q21, ..., qNN
        bcm(3+N^2:3+N^2+N-1, in_M(ij)) = LgM(:,ij);;            %Length for g1, g2, ..., gN
        bcm(3+N^2+N:3+N^2+N+M*N-1, in_M(ij)) = LhM(:,ij);       %Length for h11, h12, ..., h1N
        bcm(3+N^2+N+M*N:3+N^2+N+M*N+M-1, in_M(ij)) = LrM(:,ij); %Length for r1, r2, ..., rN
        %
        %q
        nRowM = 3+N^2+N+M*N+M;
        for j = 1:N^2
            bcm(nRowM:nRowM+LqM(j, ij)-1, in_M(ij)) = qM{j,ij};
            nRowM = nRowM + LqM(j, ij);
        end
        
        %g 
        for j = 1:N
            bcm(nRowM:nRowM+LgM(j, ij)-1, in_M(ij)) = gM{j,ij};
            nRowM = nRowM + LgM(j, ij); 
        end
        
        %h
        for j = 1:N     %1 Dirichlet BC assumed
            bcm(nRowM:nRowM+LhM(j, ij)-1, in_M(ij)) = hM{j,ij};
            nRowM = nRowM + LhM(j, ij);
        end

        %r
        j = 1;       %1, not N, 1 Dirichlet BC assumed
        bcm(nRowM:nRowM+LrM(j, ij)-1, in_M(ij)) = rM{j,ij};
        nRowM = nRowM + LrM(j, ij);
    end
end
%	<------------------     end of mixed boundaries  --------
%
%
%	----------     non-boundary segments, sub-domain segments   ----->
if in_nb ~= 0
    for ij = 1:Lnb
        bcm(1,in_nb(ij)) = 0;       %system number N =  0 for border segments
        bcm(2,in_nb(ij)) = N;       %Dirichlet boundary number M =  N
    end
end
%	<----------   end of non-boundary segments, sub-domain segments  ----
return

function gOC = geoMatrix
%   Provide geometry matrix for the organ of Corti (OC)
yBM    = -0.0088;       %y coordinate of basilar membrane (BM)

R(1)   = 0.0037;
xc(1)  = -0.0154;
yc(1)  = -1.8462e-004 - sqrt(R(1)*R(1) - (xc(1) + 0.0176)*(xc(1) + 0.0176));

x(1)   = -0.0121;
y(1)   = yc(1) - sqrt(R(1)*R(1) - (xc(1) - x(1))*(xc(1) - x(1)));
y(2)   = yBM + 2*(y(1)-yBM)/3;
x(3)   = -0.0110;
y(3)   = -2.2000e-004;

x(4)   = 0.0044;

tem1 = 0.0066 - 0.0022*sqrt(3);
tem2 = tem1 + sqrt(0.0044*0.0044-0.0033*0.0033);
tem3 = tem2 + 0.0055*sqrt(3);
tem4 = tem3 - 0.011;
tem5 = tem4-sqrt(0.0011*0.0011);
tem6 = tem5-sqrt(0.0011*0.0011);
y(4)   = y(3) + tem6*(x(4) - x(3))/0.0187;
calpha = atan((y(4) - y(3))/(x(4) - x(3)));

xc(2)  = x(4);
yc(2)  = 0.5*((x(4) - 0.0125)*(x(4) - 0.0125) + y(2)*y(2) - y(4)*y(4))/(y(2) - y(4));
R(2)   = y(4) - yc(2);

%  4 Gap_TM=Gap_RL
x(5)   = x(3);
m2     = (y(4)-y(3))/(x(4)-x(3));
b2     = 0.0020-m2*0.0044;
y(5)   = m2*x(5)+b2;
x(6)   = xc(1);
y(6)   = yc(1)-sqrt(R(1)*R(1)-(xc(1)-x(6))*(xc(1)-x(6)));

%  ----- 4 sterocillia   -------->
Xsc(1)  = 0.5*(x(3) + x(4)) + 8.8000e-004;
Ysc(1)  = y(3)+(Xsc(1)-x(3))*tan(calpha);
temp    = 0.0022;
Xsc(2)  = Xsc(1) + temp;
Ysc(2)  = y(3) + (Xsc(2) - x(3))*tan(calpha);
Xsc(3)  = Xsc(1) + 2*temp;
Ysc(3)  = y(3)+(Xsc(3)-x(3))*tan(calpha);

dOHC    = 1e-3;         %diameter of outer hair cell (OHC)
lOHC    = 90e-4;        %length of outer hair cell (OHC)
cgamma   = pi/7;
cbeta    = cgamma - calpha;
RL_OHC  = dOHC/cos(cbeta);

Xsc(7)  = Xsc(1)+RL_OHC*cos(calpha);
Ysc(7)  = y(3)+(Xsc(7)-x(3))*tan(calpha);
Xsc(8)  = Xsc(2)+RL_OHC*cos(calpha);
Ysc(8)  = y(3)+(Xsc(8)-x(3))*tan(calpha);
Xsc(9)  = Xsc(3)+RL_OHC*cos(calpha);
Ysc(9)  = y(3)+(Xsc(9)-x(3))*tan(calpha);

% 4 Pillar cells
xPi(1)  = -0.0125 + 0.025/5.4;
xPi(2)  = 0;
xPi(3)  = -0.0125 + 0.025/3;
xPi(4)  = -0.00968;
xPi(7)  = 0.00176;
xPi(5)  = -0.00484;
yPi(5)  = y(3)+(xPi(5)-x(3))*tan(calpha);
xPi(6)  =  -0.00418;
yPi(6)  = y(3)+(xPi(6)-x(3))*tan(calpha);
yPi(3)  = -0.00308;

% 4 BM
BM1     = xPi(1);
BM2     = 0.0125 + (BM1 + 0.0125);

% apex of OHCs
xOHC(1) = Xsc(1)+lOHC*sin(cgamma);
yOHC(1) = Ysc(1)-lOHC*cos(cgamma);
xOHC(2) = Xsc(7)+lOHC*sin(cgamma);
yOHC(2) = Ysc(7)-lOHC*cos(cgamma);
xOHC(3) = Xsc(2)+lOHC*sin(cgamma);
yOHC(3) = Ysc(2)-lOHC*cos(cgamma);
xOHC(4) = Xsc(8)+lOHC*sin(cgamma);
yOHC(4) = Ysc(8)-lOHC*cos(cgamma);
xOHC(5) = Xsc(3)+lOHC*sin(cgamma);
yOHC(5) = Ysc(3)-lOHC*cos(cgamma);
xOHC(6) = Xsc(9)+lOHC*sin(cgamma);
yOHC(6) = Ysc(9)-lOHC*cos(cgamma);

x(7)   = 0.0121;
y(7)   = yc(2)-sqrt(R(2)*R(2)-(xc(2)-x(7))*(xc(2)-x(7)));

% 4 inner hair cell (IHC) & Reticular lamina (RL)  ------------->
xIHC(1) = x(3)+RL_OHC*cos(calpha);
yIHC(1) = y(3)+(xIHC(1)-x(3))*tan(calpha);
xIHC(3) = xIHC(1)+0.8*RL_OHC*cos(calpha);
yIHC(3) = y(3)+(xIHC(3)-x(3))*tan(calpha);

xRL(11) = x(1)+4.5*(x(3)-x(1))/5;
yRL(11) = (xRL(11)-x(1))*(y(3)-y(1))/(x(3)-x(1))+y(1);
xRL(1)  = xRL(11)+RL_OHC*cos(calpha);
yRL(1)  = yRL(11)+(xRL(1)-xRL(11))*tan(calpha);
Lt1     = sqrt((xIHC(1)-xRL(1))*(xIHC(1)-xRL(1))+(yIHC(1)-yRL(1))*(yIHC(1)-yRL(1)));
L_IHC   = 0.0044;       

xIHC(2) = (xRL(1)-xIHC(1))*L_IHC/Lt1+xIHC(1);
yIHC(2) = (yRL(1)-yIHC(1))*L_IHC/Lt1+yIHC(1);

xIHC(4)     = xIHC(2)+0.8*RL_OHC*cos(calpha);
yIHC(4)     = yIHC(2)+RL_OHC*tan(calpha)-0.012*0.022;

[xRL(12),yRL(12)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xIHC(3),yIHC(3),xIHC(4),yIHC(4));

[xRL(2),yRL(2)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xPi(4),yBM,xPi(5),yPi(5));
[xRL(3),yRL(3)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xPi(7),yBM,xPi(6),yPi(6));
[xRL(4),yRL(4)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(1),yOHC(1),Xsc(1),Ysc(1));
[xRL(5),yRL(5)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(2),yOHC(2),Xsc(7),Ysc(7));
[xRL(6),yRL(6)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(3),yOHC(3),Xsc(2),Ysc(2));
[xRL(7),yRL(7)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(4),yOHC(4),Xsc(8),Ysc(8));
[xRL(8),yRL(8)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(5),yOHC(5),Xsc(3),Ysc(3));
[xRL(9),yRL(9)] = cintersection(xRL(11),yRL(11),xRL(1),yRL(1),xOHC(6),yOHC(6),Xsc(9),Ysc(9));

m9      = (yRL(1)-yRL(11))/(xRL(1)-xRL(11));
b9      = yRL(11)-m9*xRL(11);
a9      = 1+m9*m9;
bb9     = 2*(m9*(b9-yc(2))-xc(2));
c9      = xc(2)*xc(2)+(b9-yc(2))*(b9-yc(2))-R(2)*R(2);
x10a    = (-bb9+sqrt(bb9*bb9-4*a9*c9))/2/a9;
x10b    = (-bb9-sqrt(bb9*bb9-4*a9*c9))/2/a9;
xRL(10) = max(x10a,x10b);
yRL(10) = m9*xRL(10)+b9;

y(8)   = -3.6e-3;
x(8)   = xc(1)-sqrt(R(1)*R(1)-(y(8)-yc(1))*(y(8)-yc(1)));
x(9)   = -19.5e-3; 
y(9)   = -4.8e-3;

R(3)   = 0.0045;
[xc(3),yc(3)] = xcyc2PR(x(9), y(9), x(6), yBM, R(3),1);

x(10)   = -0.0192;
y(10)   = yc(3)-sqrt(R(3)*R(3)-(x(10)-xc(3))*(x(10)-xc(3)));
x(11)   = -0.0187;
y(11)   = yc(1)-sqrt(R(1)*R(1)-(x(11)-xc(1))*(x(11)-xc(1)));
x(12)   = -0.0185;
y(12)   = yc(3)-sqrt(R(3)*R(3)-(x(12)-xc(3))*(x(12)-xc(3)));
x(13)   = -0.0179;
y(13)   = yc(1)-sqrt(R(1)*R(1)-(x(13)-xc(1))*(x(13)-xc(1)));
x(14)   = -0.0171;
y(14)   = yc(3)-sqrt(R(3)*R(3)-(x(14)-xc(3))*(x(14)-xc(3)));
x(15)   = -0.0167;
y(15)   = yc(1)-sqrt(R(1)*R(1)-(x(15)-xc(1))*(x(15)-xc(1)));
x(16)   = -0.0135; 
y(16)   = yBM;
x(17)   = -0.0137;
y(17)   = yc(1)-sqrt(R(1)*R(1)-(x(17)-xc(1))*(x(17)-xc(1)));

% 4 rotation
xc(4) = 0;
yc(4) = 0;
rangle  = pi/8;
RR      = sqrt((x(3)-xc(2)).*(x(3)-xc(2))+(y(3)-yc(2)).*(y(3)-yc(2)));
aa      = atan((yc(2)-y(3))./(xc(2)-x(3)));
xc(2) = x(3)+RR.*cos(rangle+aa);
yc(2) = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-x(4)).*(x(3)-x(4))+(y(3)-y(4)).*(y(3)-y(4)));
aa      = atan((y(4)-y(3))./(x(4)-x(3)));
x(4)   = x(3)+RR.*cos(rangle+aa);
y(4)   = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-x(5)).*(x(3)-x(5))+(y(3)-y(5)).*(y(3)-y(5)));
aa      = pi/2;
x(5)   = x(3)+RR.*cos(rangle+aa);
y(5)   = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-x(7)).*(x(3)-x(7))+(y(3)-y(7)).*(y(3)-y(7)));
aa      = atan((y(7)-y(3))./(x(7)-x(3)));
x(7)   = x(3)+RR.*cos(rangle+aa);
y(7)   = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-Xsc).*(x(3)-Xsc)+(y(3)-Ysc).*(y(3)-Ysc));
aa      = atan((Ysc-y(3))./(Xsc-x(3)));
Xsc     = x(3)+RR.*cos(rangle+aa);
Ysc     = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xRL(1:10)).*(x(3)-xRL(1:10))+(y(3)-yRL(1:10)).*(y(3)-yRL(1:10)));
aa      = atan((yRL(1:10)-y(3))./(xRL(1:10)-x(3)));
xRL(1:10) = x(3)+RR.*cos(rangle+aa);
yRL(1:10) = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xRL(12)).*(x(3)-xRL(12))+(y(3)-yRL(12)).*(y(3)-yRL(12)));
aa      = atan((yRL(12)-y(3))./(xRL(12)-x(3)));
xRL(12) = x(3)+RR.*cos(rangle+aa);
yRL(12) = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xIHC(1:2:3)).*(x(3)-xIHC(1:2:3))+(y(3)-yIHC(1:2:3)).*(y(3)-yIHC(1:2:3)));
aa      = atan((yIHC(1:2:3)-y(3))./(xIHC(1:2:3)-x(3)));
xIHC(1:2:3) = x(3)+RR.*cos(rangle+aa);
yIHC(1:2:3) = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xOHC).*(x(3)-xOHC)+(y(3)-yOHC).*(y(3)-yOHC));
aa      = atan((yOHC-y(3))./(xOHC-x(3)));
xOHC    = x(3)+RR.*cos(rangle+aa);
yOHC    = y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xPi(5:6)).*(x(3)-xPi(5:6))+(y(3)-yPi(5:6)).*(y(3)-yPi(5:6)));
aa      = atan((yPi(5:6)-y(3))./(xPi(5:6)-x(3)));
xPi(5:6)= x(3)+RR.*cos(rangle+aa);
yPi(5:6)= y(3)+RR.*sin(rangle+aa);
RR      = sqrt((x(3)-xPi(3)).*(x(3)-xPi(3))+(y(3)-yPi(3)).*(y(3)-yPi(3)));
aa      = atan((yPi(3)-y(3))./(xPi(3)-x(3)));

xPi(3)  = x(3) + RR.*cos(rangle + aa) - 0.001;
yPi(3)  = y(3) + RR.*sin(rangle + aa);

xc(4)  = -0.01012;
tem1    = (xIHC(2)-xc(4))*(xIHC(2)-xc(4))+yIHC(2)*yIHC(2);
tem2    = (xIHC(4)-xc(4))*(xIHC(4)-xc(4))+yIHC(4)*yIHC(4);
yc(4)  = (tem1-tem2)*0.5/(yIHC(2)-yIHC(4));
R(4)   = sqrt((xIHC(2)-xc(4))*(xIHC(2)-xc(4))+(yIHC(2)-yc(4))*(yIHC(2)-yc(4)));

xc(5)  = -0.25*0.022;
tem1    = (xIHC(2)-xc(5))*(xIHC(2)-xc(5))+yIHC(2)*yIHC(2);
tem2    = (xRL(1)-xc(5))*(xRL(1)-xc(5))+yRL(1)*yRL(1);
yc(5)  = (tem1-tem2)*0.5/(yIHC(2)-yRL(1));
R(5)   = sqrt((xIHC(2)-xc(5))*(xIHC(2)-xc(5))+(yIHC(2)-yc(5))*(yIHC(2)-yc(5)));

xc(6)   = -0.6*0.022;
tem1    = (xRL(12)-xc(6))*(xRL(12)-xc(6))+yRL(12)*yRL(12);
tem2    = (xIHC(4)-xc(6))*(xIHC(4)-xc(6))+yIHC(4)*yIHC(4);
yc(6)   = (tem1-tem2)*0.5/(yRL(12)-yIHC(4));
R(6)    = sqrt((xRL(12)-xc(6))*(xRL(12)-xc(6))+(yRL(12)-yc(6))*(yRL(12)-yc(6)));

% % 4 Deiter's cells  ----------->
Dx1     = 0.005*0.022;
Dy1     = 0.05*0.022;
R_D(1)  = 0.05*0.022;
xD(1)   = xOHC(1)-Dx1;	xD(3)=xOHC(3)-Dx1;	xD(5)=xOHC(5)-Dx1;
yD(1)   = yOHC(1)-Dy1;	yD(3)=yOHC(3)-Dy1;	yD(5)=yOHC(5)-Dy1;
AA1     = 0.5*sqrt(Dx1^2+Dy1^2);
bb1     = acos(AA1/R_D(1))-atan(Dx1/Dy1);
xDc(1)  = xOHC(1)+R_D(1)*sin(bb1);	xDc(3)=xOHC(3)+R_D(1)*sin(bb1);	xDc(5)=xOHC(5)+R_D(1)*sin(bb1);
yDc(1)  = yOHC(1)-R_D(1)*cos(bb1);	yDc(3)=yOHC(3)-R_D(1)*cos(bb1);	yDc(5)=yOHC(5)-R_D(1)*cos(bb1);
xD(2)   = 0.18*0.022;	xD(4)=0.27*0.022;	xD(6)=0.36*0.022;

R_D(2)  = 0.1*0.022;
R_D(3)  = R_D(2);
AA2     = 0.5*sqrt((xOHC(2)-xD(3))^2+(yOHC(2)-yD(3))^2);
bb2     = acos(AA2/R_D(2))-atan((xD(3)-xOHC(2))/(yOHC(2)-yD(3)));
xDc(2)  = xOHC(2)-R_D(2)*sin(bb2);	xDc(4)=xOHC(4)-R_D(2)*sin(bb2);
yDc(2)  = yOHC(2)-R_D(2)*cos(bb2);	yDc(4)=yOHC(4)-R_D(2)*cos(bb2);

xD(7)   = xD(5)+(xD(5)-xD(3));
yD(7)   = yD(5)+(yD(5)-yD(3));
xD(8)   = 0.45*0.022;
xDc(6)  = xOHC(6)-R_D(2)*sin(bb2);
yDc(6)  = yOHC(6)-R_D(2)*cos(bb2);

xc(7)   = 0.7*(x(7)+BM2);
tem1    = (x(7)-xc(7))*(x(7)-xc(7))+y(7)*y(7);
tem2    = (BM2-xc(7))*(BM2-xc(7))+y(2)*y(2);
yc(7)  = (tem1-tem2)*0.5/(y(7)-y(2));
R(7)    = sqrt((x(7)-xc(7))*(x(7)-xc(7))+(y(7)-yc(7))*(y(7)-yc(7)));

gOC     = [
    2    	2   	2	    2	    2	    2	    2	    2	    2	    2	2
    x(16)  x(6)    x(1)	xRL(12)	xPi(4)	xRL(2)	xRL(3)	xPi(4)	xPi(1)	xPi(3)	xPi(2)
    xPi(4) x(6)    xRL(11)	xRL(2)	xRL(2)	xRL(3)	xPi(7)	xPi(1)	xPi(3)	xPi(2)	xPi(7)
    yBM 	yBM  	y(1)	yRL(12)	yBM     yRL(2)	yRL(3)	yBM     yBM     yPi(3)	yBM
    yBM 	y(6) 	yRL(11)	yRL(2)	yRL(2)	yRL(3)	yBM     yBM     yPi(3)	yBM     yBM
    1    	29    	0	    15	    1	    16	    3	    2	    2	    2	    2
    0    	30    	1	    1	    2	    2	    2	    0	    25	    25	    0
    0    	0    	0	    0	    0	    0	    0	    0	    0	    0	    0
    0    	0    	0	    0	    0	    0	    0	    0	    0	    0	    0
    0    	0    	0	    0	    0	    0	    0	    0	    0	    0	    0
];

R_OHCtip    = 0.001*3;
[xcOHCtip,ycOHCtip] = xcyc2PR(xOHC(1),yOHC(1),xOHC(2),yOHC(2),R_OHCtip,1);
gOC(:,12:15)=[
    2	    2	    2	    1
    xRL(4)	xRL(5)	xRL(4)	xOHC(1)
    xRL(5)	xOHC(2)	xOHC(1)	xOHC(2)
    yRL(4)	yRL(5)	yRL(4)	yOHC(1)
    yRL(5)	yOHC(2)	yOHC(1)	yOHC(2)
    18	    10	    4	    4
    4	    4  	    3	    7
    0	    0	    0	    xcOHCtip
    0	    0	    0	    ycOHCtip
    0	    0	    0	    R_OHCtip
];

[xcOHCtip,ycOHCtip]=xcyc2PR(xOHC(3),yOHC(3),xOHC(4),yOHC(4),R_OHCtip,1);
gOC(:,16:19)=[
    2	    2	    2	    1
    xRL(6)	xRL(7)	xRL(6)	xOHC(3)
    xRL(7)	xOHC(4)	xOHC(3)	xOHC(4)
    yRL(6)	yRL(7)	yRL(6)	yOHC(3)
    yRL(7)	yOHC(4)	yOHC(3)	yOHC(4)
    20	    11	    5	    5
    5	    5	    10	    8
    0	    0	    0	    xcOHCtip
    0	    0	    0	    ycOHCtip
    0	    0	    0	    R_OHCtip
];

[xcOHCtip,ycOHCtip]=xcyc2PR(xOHC(5),yOHC(5),xOHC(6),yOHC(6),R_OHCtip,1);
gOC(:,20:23)=[
    2	    2	    2	    1
    xRL(8)	xRL(9)	xRL(8)	xOHC(5)
    xRL(9)	xOHC(6)	xOHC(5)	xOHC(6)
    yRL(8)	yRL(9)	yRL(8)	yOHC(5)
    yRL(9)	yOHC(6)	yOHC(5)	yOHC(6)
    22	    12	    6	    6
    6	    6	    11	    9
    0	    0	    0	    xcOHCtip
    0	    0	    0	    ycOHCtip
    0	    0	    0	    R_OHCtip
];

gOC(:,24:34)=[
    2	    2	    2	    2	    2	    2	    2	    1	    2	    2	    2
    xRL(3)	xPi(7)	xD(2)	xD(4)	xD(6)	xD(8)	BM2	    x(7)	xRL(5)	xRL(7)	xRL(9)
    xRL(4)	xD(2)	xD(4)	xD(6)	xD(8)	BM2	    BM2	    xRL(10)	xRL(6)	xRL(8)	xRL(10)
    yRL(3)	yBM	yBM	yBM	yBM	yBM	yBM	y(7)	yRL(5)	yRL(7)	yRL(9)
    yRL(4)	yBM	yBM	yBM	yBM	yBM	y(2)	yRL(10)	yRL(6)	yRL(8)	yRL(10)
    17	    3	    7	    8	    9	    12	    12	    12	    19	    21	    23
    3	    0	    0	    0	    0	    0	    0	    0	    10	    11	    12
    0	    0	    0	    0	    0	    0	    0	    xc(2)	0	    0	    0
    0	    0	    0	    0	    0	    0	    0	    yc(2)	0	    0	    0
    0	    0	    0	    0	    0	    0	    0	    R(2)	0	    0	    0
    
];

gOC(:,35:44)=[
    1	    2	    1	    2	    1	    2	    1	    1	    1	    2
    xOHC(1)	xD(1)	xOHC(3)	xD(3)	xOHC(5)	xD(5)	xD(3)	xD(5)	xD(7)	xD(7)	
    xD(1)	xD(2)	xD(3)	xD(4)	xD(5)	xD(6)	xOHC(2)	xOHC(4)	xOHC(6)	xD(8)
    yOHC(1)	yD(1)	yOHC(3)	yD(3)	yOHC(5)	yD(5)	yD(3)	yD(5)	yD(7)	yD(7)
    yD(1)	yBM	yD(3)	yBM	yD(5)	yBM	yOHC(2)	yOHC(4)	yOHC(6)	yBM
    7	    7	    8	    8	    9	    9	    7	    8	    9	    12
    3	    3	    10	    7	    11	    8	    10	    11	    12	    9
    xDc(1)	0	    xDc(3)	0	    xDc(5)	0	    xDc(2)	xDc(4)	xDc(6)	0
    yDc(1)	0	    yDc(3)	0	    yDc(5)	0	    yDc(2)	yDc(4)	yDc(6)	0
    R_D(1)	0	    R_D(1)	0	    R_D(1)	0	    R_D(2)	R_D(2)	R_D(3)	0
];

gOC(:,45:48)=[
    1	    1	    1	    2
    x(17)	xRL(1)	xIHC(2)	xRL(11)
    x(1)	xIHC(2)	xIHC(4)	xRL(1)
    y(17)	yRL(1)	yIHC(2)	yRL(11)
    y(1)	yIHC(2)	yIHC(4)	yRL(1)
    0	    13	    13	    14
    1	    1	    1	    1
    xc(1)	xc(5)	xc(4)	0
    yc(1)	yc(5)	yc(4)	0
    R(1)	R(5)	R(4)	0
];

gOC(:,49:58)=[
    2	    2	    2	    2	    2	    2	    2	    2	    2	    2
    x(3)	xIHC(3)	xPi(5)	xPi(6)	Xsc(1)	Xsc(7)	Xsc(2)	Xsc(8)	Xsc(3)	Xsc(9)
    xIHC(1)	xPi(5)	xPi(6)	Xsc(1)	Xsc(7)	Xsc(2)	Xsc(8)	Xsc(3)	Xsc(9)	x(4)
    y(3)	yIHC(3)	yPi(5)	yPi(6)	Ysc(1)	Ysc(7)	Ysc(2)	Ysc(8)	Ysc(3)	Ysc(9)
    yIHC(1)	yPi(5)	yPi(6)	Ysc(1)	Ysc(7)	Ysc(2)	Ysc(8)	Ysc(3)	Ysc(9)	y(4)
    0	    0	    0	    0	    0	    0	    0	    0	    0	    0
    14	    15	    16  	17	    18	    19	    20  	21  	22	    23
    0	    0	    0	    0	    0	    0	    0   	0   	0	    0
    0	    0	    0   	0   	0	    0	    0	    0   	0	    0	
    0   	0	    0   	0   	0   	0	    0   	0   	0	    0
];
gOC(:,59:69)=[
    2	    2	    2	    2	    2	    2	    2	    2	    2	    2	    1
    xRL(11)	xRL(1)	xRL(2)	xRL(3)	xRL(4)	xRL(5)	xRL(6)	xRL(7)	xRL(8)	xRL(9)	xRL(10)
    x(3)	xIHC(1)	xPi(5)	xPi(6)	Xsc(1)	Xsc(7)	Xsc(2)	Xsc(8)	Xsc(3)	Xsc(9)	x(4)
    yRL(11)	yRL(1)	yRL(2)	yRL(3)	yRL(4)	yRL(5)	yRL(6)	yRL(7)	yRL(8)	yRL(9)	yRL(10)
    y(3)	yIHC(1)	yPi(5)	yPi(6)	Ysc(1)	Ysc(7)	Ysc(2)	Ysc(8)	Ysc(3)	Ysc(9)	y(4)
    0	    14	    15	    16	    17	    18	    19	    20	    21	    22	    23
    14	    24	    16	    17	    18	    19	    20	    21	    22	    23	    0
    0	    0	    0	    0	    0	    0   	0	    0	    0	    0	    xc(2)
    0	    0	    0	    0   	0	    0   	0	    0	    0	    0	    yc(2)
    0	    0	    0	    0	    0	    0	    0	    0   	0   	0	    R(2)
];
gOC(:,70:75)=[
    1	    2	    1	    2	    2	    2
    x(7)	xRL(12)	xIHC(4)	xIHC(1)	xRL(1)	xPi(1)
    BM2	    xIHC(3)	xRL(12)	xIHC(3)	xRL(12)	xPi(2)
    y(7)	yRL(12)	yIHC(4)	yIHC(1)	yRL(1)	yBM
    y(2)	yIHC(3)	yRL(12)	yIHC(3)	yRL(12)	yBM
    0	    24	    13	    0	    24	    25
    12	    15	    1	    24	    13	    0
    xc(7)	0	    xc(6)	0	    0	    0
    yc(7)	0	    yc(6)	0	    0	    0
    R(7)	0	    R(6)	0	    0	    0
];

gOC(:,76:90)=[
    1	    2       2       1       1       1       2       1       1       2       1       1       2       1       2
    x(9)	x(8)   x(10)   x(8)     x(11)   x(10)   x(12)   x(13)   x(12)   x(14)   x(15)   x(14)   x(16)   x(6)    x(6)
    x(10)   x(9)   x(11)   x(11)    x(13)   x(12)   x(13)   x(15)   x(14)   x(15)   x(6)    x(6)    x(17)   x(17)   x(16)
    y(9)    y(8)   y(10)   y(8)     y(11)   y(10)   y(12)   y(13)   y(12)   y(14)   y(15)   y(14)   yBM     y(6)    yBM
    y(10)	y(9)   y(11)   y(11)    y(13)   y(12)   y(13)   y(15)   y(14)   y(15)   y(6)    yBM     y(17)   y(17)   yBM
    26	    26      26      0       0       27      27      0       28      28      0       29      30      0       30
    0	    0       27      26      27      0       28      28      0       29      29      0       1       30      0
    xc(3)   0       0       xc(1)   xc(1)   xc(3)   0       xc(1)   xc(3)   0       xc(1)   xc(3)   0       xc(1)   0
    yc(3)   0       0       yc(1)   yc(1)   yc(3)   0       yc(1)   yc(3)   0       yc(1)   yc(3)   0       yc(1)   0
    R(3)    0       0       R(1)    R(1)    R(3)    0       R(1)    R(3)    0       R(1)    R(3)    0       R(1)    0
];
%
%Note:  The numbers in rows # 6 & 7 indicate the sub-domain numbers on the
%       left and right of the corresponding element, respectively. 0 means no material (outside
%       the geometrical domain)
return
%
%
function [g1N, g2N] = g1g2(in_N, c_p, c_dpdx, c_dpdy, tao)
%	======  Providing g1, g2 for non-zero Nmann boundary condition of OC  ===== 
%	======  Note: for simplicity, the tangential stress duo to wall movement is ignored  =======
L_N = length(in_N);
for ij=1:L_N
	p   = ['(polyval(', mat2str(c_p(:)), ',x))'];                % pressure
	px  = ['(polyval(', mat2str(c_dpdx(:)), ',x))'];             % dpdx, pressure gradient in x dicrection
	py  = ['(polyval(', mat2str(c_dpdy(:)), ',x))'];             % dpdy
    ps  = ['(-', px, '.*ny+', py, '.*nx)'];                      % dpds
    tao1 = ['(', mat2str(tao), ')*', ps];                        % tao1: tangential stress

    %g1
	temp_g1 = ['-', p, '.*nx-', tao1,'.*ny'];
	g1N(1:length(temp_g1),ij) = (temp_g1)';

    %g2
	temp_g2 = ['-', p, '.*ny+', tao1, '.*nx'];
	g2N(1:length(temp_g2),ij) = (temp_g2)';
end
return

function [x,y] = cintersection(x1,y1,x2,y2,x3,y3,x4,y4)

%	====== intersection of two lines =====>

% 	line 1 passes (x1,y1) and (x2,y2)
% 	line 2 passes (x3,y3) and (x4,y4)

m1  = (y2-y1)/(x2-x1);
b1  = y1-m1*x1;
if x3 ~= x4
	m2  =(y4-y3)/(x4-x3);
	b2  = y3-m2*x3;

	x   =(b1-b2)/(m2-m1);
	y   =(m2*b1-m1*b2)/(m2-m1);
else
	x   = x3;
	y   = m1*(x-x1)+y1;
end

function [xc, yc] = xcyc2PR(x1, y1, x2, y2, R, id)
% 4 finding the center point of a cricle passing 2 points (x1,y1), (x2,y2) with a
% radius R
% id = 1 ---> up or right;  id = 2 -----> lower or left

if x1 ~= x2
    alpha   = (y2-y1)/(x1-x2);
    beta    = 0.5*(x1*x1-x2*x2+y1*y1-y2*y2)/(x1-x2);
    a       = alpha*alpha+1;
    b       = 2*(-y1+alpha*(beta-x1));
    c       = (beta-x1)*(beta-x1)+y1*y1-R*R;
    if id == 1
        yc  = 0.5*(-b+sqrt(b*b-4*a*c))/a;
    else
        yc  = 0.5*(-b-sqrt(b*b-4*a*c))/a;
    end
    xc      = alpha*yc+beta;
else
    yc      = 0.5*(y1+y2);
    h       = 0.5*abs(y2-y1);
    if id == 1
        xc  = x1+sqrt(R*R-h*h);
    else
        xc  = x1-sqrt(R*R-h*h);
    end
end

Contact us