ChebyshevTools: [BIG_OP,BIG_rhs]=bc_2d_new(P,Q,op_var,bc_type,bc_vals,R,Z) - File Exchange

No BSD License

Highlights from
ChebyshevTools

A=expD(M,ord,mod);
A_res=restrictA(A,m) this function extracts the submatrix A(1+m:end,1:end-m)
B=inv_deriv1(M); This function returns one matrix, for solving
B=inv_derivP(M,P) this function returns the Pth quasi inverse for Chebycheb differentiation
B=pseudo_eye(M,P) this function returns the pseudo-identity matrix I^(P), where there are
D=deriv(M); returns the first order chebyshev spectral differentiation matrix
D=deriv2(M);
P=cheb2poly(U) this function converts the spectral representation of a function in Chebyshev polynomials
S=stencil_mat(M,typ) this function returns the stencil matrix associated with the related
Y=cheb2phys2(X);
Y=cheb2phys3(X); 3d transform chebyshev to physical space
Y=phys2cheb2(X);
Y=phys2cheb3(X); 3d transform physical space to chebyshev
[BIG_OP,BIG_rhs]=bc_2d_new(P,... Generate the tau line boundary conditions and the appropriate rhs vector
[BIG_OP,BIG_rhs]=bc_2d_new_to... Generate the tau line boundary conditions and the appropriate rhs vector
[OP,rhs]=bc_1d(M,bc_type,bc_v...
f=galerkin2cheb(g,typ) this function transforms from the galerkin basis to the chebyshev basis
f=galerkin2cheb2(g,typ) this function transforms from the galerkin basis to the chebyshev basis in
f=galerkin2cheb3(g,typ) this function transforms from the galerkin basis to the chebyshev basis in
g=cheb2galerkin(f,typ) this function transforms from the chebyshev basis to the galerkin basis
g=cheb2galerkin2(f,typ) this function transforms from the galerkin basis to the chebyshev basis in
g=cheb2galerkin3(f,typ) this function transforms from the galerkin basis to the chebyshev basis in
test_problem(M,N,tc) This function is run from the command line taking argumetns
v=getdiag(A,n); extract the nth diag from matrix A
x=cheb_grid(N); this function returns a chebyshev grid on the interval [-1,1]
y=cheb2phys(F); Written by Mike Watson, November 29
y=definite_integral_cart(F,X1... Calculate the definite integral of a function defined by spectral
y=myC(p)
y=phys2cheb(f); Written by Mike Watson, November 29
deriv3.m
View all files

from ChebyshevTools by Michael Watson
Tool box for solving ODE/PDEs using spectral Chebyshev differentiation matrices.

[BIG_OP,BIG_rhs]=bc_2d_new(P,Q,op_var,bc_type,bc_vals,R,Z)

function [BIG_OP,BIG_rhs]=bc_2d_new(P,Q,op_var,bc_type,bc_vals,R,Z)
alpha=5e-1;
% Generate the tau line boundary conditions and the appropriate rhs vector


% Inputs
% P : size in the X/R direction
% Q : size in the Y/Z direction
% op_var : variable on which the operator acts
% bc_type : vector of strings containing the bc type
%   Dirichlet
%   D_R = Diriclet boundary condition on the right
%   D_L = Diriclet boundary condition on the left
%   D_P = Summation of dirichlet conditions
%   D_M = Difference of dirichlet conditions
%   Neumann
%   N_R = Neumann on the right
%   N_L = Neumann on the left
%   N_P = Summation of neumann conditions
%   N_M = Difference of neumann conditions
%   Stress Free
%   SF1_R = Stress free on the right
%   SF2_L = Stress free on the left
%   SF2_R = Stress free for azimuthal velocity on the right
%   SF2_L = Stress free for azimuthal velocity on the left
%   SF3_R = Stress free for potential functions on the right
%   SF3_L = Stress free for potential functions on the left
%   D3_R = Third Derivative on the right
%   D3_L = Third Derivative on the left
%   D3_P = Summation of Third Derivative conditions
%   D3_M = Difference of Third Derivative conditions
%   I_RRR = Integral constraint on Third Derivative
%   I_DZ = Integral constraint on First Derivative in Z

% R,Z : left/right, inner/outer, lower/upper coordinates


len=length(bc_type);
[c,len_v]=size(bc_vals);

if (len ~= len_v)
   disp(sprintf('# of bc_types and bc_vals do not match'))
   return
end

% scaling parameters
RO=R(2);
RI=R(1);
%F1=.5*(RO-RI);G1=.5*(RO+RI);

HT=Z(2);
HB=Z(1);
%F2=(HT-HB)/2;  G2=(HT+HB)/2;


% Determine "operator" vs "action"
if (strcmp(op_var,'R') | strcmp(op_var,'X'))
   M=P;
   N=Q; 
   rhs=zeros(Q,len);
   F1=.5*(RO-RI);G1=.5*(RO+RI);
   F2=(HT-HB)/2;  G2=(HT+HB)/2;
elseif (strcmp(op_var,'Z') | strcmp(op_var,'Y'))
   M=Q;
   N=P;
   rhs=zeros(N,len);
   F2=.5*(RO-RI);G2=.5*(RO+RI);
   F1=(HT-HB)/2;  G1=(HT+HB)/2;
else
   disp(sprintf('Undefined operational variable'))
end

BIG_OP=sparse(zeros(N*M,N*M));


%build the appropriate 1d case
for l=1:len
   typ=bc_type(l);
   bc_pos=M-(len-l);
   rhs_pos=l;
   
   OP=zeros(M);
   I=eye(N);

   if strcmp(typ,'D_R')
       %Dirichlet on the Right / Top
       for j=1:M
           OP(bc_pos,j)=1;
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'D_L')
       %Dirichlet on the Left / Bottom
       for j=1:M
           OP(bc_pos,j)=1*(-1)^(j+1);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'N_R')
       %Neumann on the Right / Top
       for j=2:M
           OP(bc_pos,j)=(j-1)^2;
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'N_L')
       %Neumann on the Left / Bottom
       for j=2:M
           OP(bc_pos,j)=(j-1)^2*(-1)^(j);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'D2_R')  % second derivative
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'D2_L')
       for j=1:M
           OP(bc_pos,j)=(1/3)*((j-1)^2*((j-1)^2-1))*(-1)^(j+1);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'D3_R')
       % thrid derivative at right/top
       for j=1:M
           OP(bc_pos,j)=(1/15)*(j-1)^2*((j-1)^2-1)*((j-1)^2-4);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^3;
   elseif strcmp(typ,'D3_L')
       % thrid derivative at left/bottom
       for j=1:M
           OP(bc_pos,j)=(1/15)*((j-1)^2*((j-1)^2-1)*((j-1)^2-4))*(-1)^(j);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^3;
   elseif strcmp(typ,'SF1_R')
        %f''-f'/r on the RIght / Top
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1)-(F1/RO)*((j-1)^2);  %RO
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'SF1_L')
       %f''-f'/r on the Left / Bottom
       for j=1:M
           OP(bc_pos,j)=((1/3)*(j-1)^2*((j-1)^2-1))*(-1)^(j+1)-((F1/RI)*(j-1)^2)*(-1)^(j);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'SF2_R')  
       %stress free for v at right / top : v_r - (v)/r
       for j=1:M
           OP(bc_pos,j)=(j-1)^2-(F1/RO)*(1);  %RO
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'SF2_L') 
       %stress free for v at left /  bottom : v_r - (v)/r
       for j=1:M
           OP(bc_pos,j)=((j-1)^2*(-1)^(j))-((F1/RI)*(-1)^(j+1));
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'SF3_R')  
       %stress free modified for potential function at right phi_rrr + (phi_rr)/r - (phi_r)/r^2
       for j=1:M
           OP(bc_pos,j)=(1/15)*(j-1)^2*((j-1)^2-1)*((j-1)^2-4)+(F1/RO)*(1/3)*(j-1)^2*((j-1)^2-1)-(F1^2/RO^2)*((j-1)^2);  %RO
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^3;
   elseif strcmp(typ,'SF3_L')
       %stress free modified for potential function at left phi_rrr + (phi_rr)/r - (phi_r)/r^2
       for j=1:M
           OP(bc_pos,j)=(1/15)*((j-1)^2*((j-1)^2-1)*((j-1)^2-4))*(-1)^(j)+(F1/RI)*((1/3)*(j-1)^2*((j-1)^2-1))*(-1)^(j+1)-((F1^2/RI^2)*(j-1)^2)*(-1)^(j);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^3;
   elseif strcmp(typ,'SF4_R')
       %stress free for g at right g_rr + (g_r)/r - (g)/r^2
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1)+(F1/RO)*(j-1)^2-(F1^2/RO^2)*1;
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'SF4_L')
       %stress free for g at left g_rr + (g_r)/r - (g)/r^2
       for j=1:M
           OP(bc_pos,j)=(1/3)*((j-1)^2*((j-1)^2-1))*(-1)^(j+1)+(F1/RI)*(j-1)^2*(-1)^(j)-(F1^2/RI^2)*(-1)^(j+1);
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1^2;
   elseif strcmp(typ,'C1')
       % specify the value of C1
       OP(bc_pos,:)=zeros(1,M);
       OP(bc_pos,2)=1;
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'INTC')  
       % specify the integral [f(r)] or [f(z)]
       prod=zeros(1,M);
       % cartesian component
       for j=1:2:M
           prod(j)=2*(-1/(j*(j-2))); 
       end
       OP(bc_pos,:)=prod*(F1);
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'INTV')  
       % specify the integral [r^2*f(r)]
       prod=zeros(1,M);
       % cartesian component
       for j=1:2:M
           prod(j)=2*(-1/(j*(j-2))); 
       end
       OP(bc_pos,:)=prod*(G1^2*F1);
       prod=zeros(1,M);
       % cylindrical component
       for j=1:2:M-1
          prod(j+1)=(-1)*((1/((j-2)*j))+(1/(j*(j+2)))); 
       end
       OP(bc_pos,:)=prod*(2*G1*F1^2)+OP(bc_pos,:);
       prod=zeros(1,M);
       % "r^2" component
       for j=0:2:M-1
          prod(j+1)=(-1)*((1/(2*(j-1)*(j-3)))+(1/((j-1)*(j+1)))+(1/(2*(j+1)*(j+3)))); 
       end
       OP(bc_pos,:)=prod*F1^3+OP(bc_pos,:);
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'ASYM_TB')
       % Anti-Symmetric Top/Bottom
       for j=1:2:M
           OP(bc_pos,j)=2;
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l);
   elseif strcmp(typ,'DDDINTC')
       % differentiated integral constraint for g
       % Intc for Z
       prod=zeros(1,M);
       % cartesian component
       for j=1:2:M
           prod(j)=2*(-1/(j*(j-2))); 
       end
       OP(bc_pos,:)=prod*(F1);
       % DPDDP in r
       I_R=eye(N);
       R_op=(F2*r_op(N)+G2*I_R);
       DR_op=deriv(N);
       ACT=R_op^3*DR_op^3+2*F2*R_op^2*DR_op^2-F2^2*R_op*DR_op+F2^3*I_R;
       rhs(:,rhs_pos)=F2^3*R_op^3*bc_vals(:,l);
   elseif strcmp(typ,'DR_Z_T')
       % specify DR on top
       
       %Z on top
       for j=1:M
           OP(bc_pos,j)=1;
       end
       
       % Dr
       I_R=eye(N);
       R_op=(F2*r_op(N)+G2*I_R);
       DR_op=deriv(N);
       ACT=DR_op;
       rhs(:,rhs_pos)=F2*R_op*bc_vals(:,l);
   elseif strcmp(typ,'DRP_Z_T')
       % specify DR on top
       
       %Z on top
       for j=1:M
           OP(bc_pos,j)=1;
       end
       
       % DRP = d_r(f)+f/r
       I_R=eye(N);
       R_op=(F2*r_op(N)+G2*I_R);
       DR_op=deriv(N);
       ACT=DR_op;+F2*I_R;
       rhs(:,rhs_pos)=F2*R_op*bc_vals(:,l);
   elseif strcmp(typ,'SFDEGEN_R')  
       %stress free for v at right / top : v_r - (v)/r
       for j=1:M
           OP(bc_pos,j)=(j-1)^2-(F1/RO)*(1)+1*alpha;  %RO
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'SFDEGEN_L') 
       %stress free for v at left /  bottom : v_r - (v)/r
       for j=1:M
           OP(bc_pos,j)=((j-1)^2*(-1)^(j))-((F1/RI)*(-1)^(j+1))+(-1)^(j+1)*alpha;
       end
       ACT=diag(ones(N,1));
       rhs(:,rhs_pos)=bc_vals(:,l)*F1;
   elseif strcmp(typ,'ZERO')
       % Leave a line intentially full of zeros
       % do nothing for this entry
       ACT=zeros(N);
   else
       disp(sprintf('Invalid bc_type at entry : %d ',l))
       disp(typ)
   end
   if (op_var =='R')|(op_var=='X')
      BIG_OP=BIG_OP+kron(sparse(OP),sparse(ACT));
   else
      BIG_OP=BIG_OP+kron(sparse(ACT),sparse(OP));
   end
end

% 
% figure(75)
% spy(BIG_OP)
% 
% pause

BIG_rhs=zeros(M*N,1);
for l=1:len
     temp=zeros(M,1);
     temp(M-(len-l),1)=1;
     if (op_var =='R')|(op_var=='X')
         BIG_rhs=BIG_rhs+kron(temp,rhs(:,l));
     else
         BIG_rhs=BIG_rhs+kron(rhs(:,l),temp);
     end
end

% 
% figure(76)
% spy(BIG_rhs)

Highlights from ChebyshevTools

Highlights from
ChebyshevTools