| [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hstcrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w); |
function [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hstcrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w);
persistent deltax deltay delxsq delysq i id2 id3 id4 ierr1 j mp mperod np nperod s st2 twdelx twdely twdysq ;
if isempty(deltax), deltax=0; end;
if isempty(deltay), deltay=0; end;
if isempty(delxsq), delxsq=0; end;
if isempty(delysq), delysq=0; end;
if isempty(s), s=0; end;
if isempty(st2), st2=0; end;
if isempty(twdelx), twdelx=0; end;
if isempty(twdely), twdely=0; end;
if isempty(twdysq), twdysq=0; end;
%***BEGIN PROLOGUE HSTCRT
%***PURPOSE Solve the standard five-point finite difference
% approximation on a staggered grid to the Helmholtz equation
% in Cartesian coordinates.
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B1A1A
%***TYPE SINGLE PRECISION (HSTCRT-S)
%***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% HSTCRT solves the standard five-point finite difference
% approximation on a staggered grid to the Helmholtz equation in
% Cartesian coordinates
%
% (d/dX,dU/dX) + (d/dY,dU/dY) + LAMBDA*U = F(X,Y)
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
% * * * * * * * * Parameter Description * * * * * * * * * *
%
% * * * * * * On Input * * * * * *
%
% A,B
% The range of X, i.e. A .LE. X .LE. B. A must be less than B.
%
% M
% The number of grid points in the interval (A,B). The grid points
% in the X-direction are given by X(I) = A + (I-0.5)dX for
% I=1,2,...,M where dX =(B-A)/M. M must be greater than 2.
%
% MBDCND
% Indicates the type of boundary conditions at X = A and X = B.
%
% = 0 If the solution is periodic in X,
% U(M+I,J) = U(I,J).
%
% = 1 If the solution is specified at X = A and X = B.
%
% = 2 If the solution is specified at X = A and the derivative
% of the solution with respect to X is specified at X = B.
%
% = 3 If the derivative of the solution with respect to X is
% specified at X = A and X = B.
%
% = 4 If the derivative of the solution with respect to X is
% specified at X = A and the solution is specified at X = B.
%
% BDA
% A one-dimensional array of length N that specifies the boundary
% values (if any) of the solution at X = A. When MBDCND = 1 or 2,
%
% BDA(J) = U(A,Y(J)) , J=1,2,...,N.
%
% When MBDCND = 3 or 4,
%
% BDA(J) = (d/dX)U(A,Y(J)) , J=1,2,...,N.
%
% BDB
% A one-dimensional array of length N that specifies the boundary
% values of the solution at X = B. When MBDCND = 1 or 4
%
% BDB(J) = U(B,Y(J)) , J=1,2,...,N.
%
% When MBDCND = 2 or 3
%
% BDB(J) = (d/dX)U(B,Y(J)) , J=1,2,...,N.
%
% C,D
% The range of Y, i.e. C .LE. Y .LE. D. C must be less
% than D.
%
% N
% The number of unknowns in the interval (C,D). The unknowns in
% the Y-direction are given by Y(J) = C + (J-0.5)DY,
% J=1,2,...,N, where DY = (D-C)/N. N must be greater than 2.
%
% NBDCND
% Indicates the type of boundary conditions at Y = C
% and Y = D.
%
% = 0 If the solution is periodic in Y, i.e.
% U(I,J) = U(I,N+J).
%
% = 1 If the solution is specified at Y = C and Y = D.
%
% = 2 If the solution is specified at Y = C and the derivative
% of the solution with respect to Y is specified at Y = D.
%
% = 3 If the derivative of the solution with respect to Y is
% specified at Y = C and Y = D.
%
% = 4 If the derivative of the solution with respect to Y is
% specified at Y = C and the solution is specified at Y = D.
%
% BDC
% A one dimensional array of length M that specifies the boundary
% values of the solution at Y = C. When NBDCND = 1 or 2,
%
% BDC(I) = U(X(I),C) , I=1,2,...,M.
%
% When NBDCND = 3 or 4,
%
% BDC(I) = (d/dY)U(X(I),C), I=1,2,...,M.
%
% When NBDCND = 0, BDC is a dummy variable.
%
% BDD
% A one-dimensional array of length M that specifies the boundary
% values of the solution at Y = D. When NBDCND = 1 or 4,
%
% BDD(I) = U(X(I),D) , I=1,2,...,M.
%
% When NBDCND = 2 or 3,
%
% BDD(I) = (d/dY)U(X(I),D) , I=1,2,...,M.
%
% When NBDCND = 0, BDD is a dummy variable.
%
% ELMBDA
% The constant LAMBDA in the Helmholtz equation. If LAMBDA is
% greater than 0, a solution may not exist. However, HSTCRT will
% attempt to find a solution.
%
% F
% A two-dimensional array that specifies the values of the right
% side of the Helmholtz equation. For I=1,2,...,M and J=1,2,...,N
%
% F(I,J) = F(X(I),Y(J)) .
%
% F must be dimensioned at least M X N.
%
% IDIMF
% The row (or first) dimension of the array F as it appears in the
% program calling HSTCRT. This parameter is used to specify the
% variable dimension of F. IDIMF must be at least M.
%
% W
% A one-dimensional array that must be provided by the user for
% work space. W may require up to 13M + 4N + M*INT(log2(N))
% locations. The actual number of locations used is computed by
% HSTCRT and is returned in the location W(1).
%
%
% * * * * * * On Output * * * * * *
%
% F
% contains the solution U(I,J) of the finite difference
% approximation for the grid point (X(I),Y(J)) for
% I=1,2,...,M, J=1,2,...,N.
%
% PERTRB
% If a combination of periodic or derivative boundary conditions is
% specified for a Poisson equation (LAMBDA = 0), a solution may not
% exist. PERTRB is a constant, calculated and subtracted from F,
% which ensures that a solution exists. HSTCRT then computes this
% solution, which is a least squares solution to the original
% approximation. This solution plus any constant is also a
% solution; hence, the solution is not unique. The value of PERTRB
% should be small compared to the right side F. Otherwise, a
% solution is obtained to an essentially different problem. This
% comparison should always be made to insure that a meaningful
% solution has been obtained.
%
% IERROR
% An error flag that indicates invalid input parameters.
% Except for numbers 0 and 6, a solution is not attempted.
%
% = 0 No error
%
% = 1 A .GE. B
%
% = 2 MBDCND .LT. 0 or MBDCND .GT. 4
%
% = 3 C .GE. D
%
% = 4 N .LE. 2
%
% = 5 NBDCND .LT. 0 or NBDCND .GT. 4
%
% = 6 LAMBDA .GT. 0
%
% = 7 IDIMF .LT. M
%
% = 8 M .LE. 2
%
% Since this is the only means of indicating a possibly
% incorrect call to HSTCRT, the user should test IERROR after
% the call.
%
% W
% W(1) contains the required length of W.
%
% *Long Description:
%
% * * * * * * * program Specifications * * * * * * * * * * * *
%
% Dimension of BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
if isempty(i), i=0; end;
if isempty(id2), id2=0; end;
if isempty(id3), id3=0; end;
if isempty(id4), id4=0; end;
if isempty(ierr1), ierr1=0; end;
if isempty(j), j=0; end;
if isempty(mp), mp=0; end;
if isempty(mperod), mperod=0; end;
if isempty(np), np=0; end;
if isempty(nperod), nperod=0; end;
% Arguments W(See argument list)
%
% Latest June 1, 1977
% Revision
%
% Subprograms HSTCRT,POISTG,POSTG2,GENBUN,POISD2,POISN2,POISP2,
% Required COSGEN,MERGE,TRIX,TRI3,PIMACH
%
% Special NONE
% Conditions
%
% Common NONE
% Blocks
%
% I/O NONE
%
% Precision Single
%
% Specialist Roland Sweet
%
% Language FORTRAN
%
% History Written by Roland Sweet at NCAR in January , 1977
%
% Algorithm This subroutine defines the finite-difference
% equations, incorporates boundary data, adjusts the
% right side when the system is singular and calls
% either POISTG or GENBUN which solves the linear
% system of equations.
%
% Space 8131(decimal) = 17703(octal) locations on the
% Required NCAR Control Data 7600
%
% Timing and The execution time T on the NCAR Control Data
% Accuracy 7600 for subroutine HSTCRT is roughly proportional
% to M*N*log2(N). Some typical values are listed in
% the table below.
% The solution process employed results in a loss
% of no more than FOUR significant digits for N and M
% as large as 64. More detailed information about
% accuracy can be found in the documentation for
% subroutine POISTG which is the routine that
% actually solves the finite difference equations.
%
%
% M(=N) MBDCND NBDCND T(MSECS)
% ----- ------ ------ --------
%
% 32 1-4 1-4 56
% 64 1-4 1-4 230
%
% Portability American National Standards Institute Fortran.
% The machine dependent constant PI is defined in
% function PIMACH.
%
% Required COS
% Resident
% Routines
%
% Reference Schumann, U. and R. Sweet,'A Direct Method For
% The Solution Of Poisson's Equation With Neumann
% Boundary Conditions On A Staggered Grid Of
% Arbitrary Size,' J. COMP. PHYS. 20(1976),
% PP. 171-182.
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES U. Schumann and R. Sweet, A direct method for the
% solution of Poisson's equation with Neumann boundary
% conditions on a staggered grid of arbitrary size,
% Journal of Computational Physics 20, (1976),
% pp. 171-182.
%***ROUTINES CALLED GENBUN, POISTG
%***REVISION HISTORY (YYMMDD)
% 801001 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE HSTCRT
%
%
f_shape=size(f);f=reshape([f(:).',zeros(1,ceil(numel(f)./prod([idimf])).*prod([idimf])-numel(f))],idimf,[]);
bda_shape=size(bda);bda=reshape(bda,1,[]);
bdb_shape=size(bdb);bdb=reshape(bdb,1,[]);
bdc_shape=size(bdc);bdc=reshape(bdc,1,[]);
bdd_shape=size(bdd);bdd=reshape(bdd,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
%***FIRST EXECUTABLE STATEMENT HSTCRT
ierror = 0;
if( a>=b )
ierror = 1;
end;
if( mbdcnd<0 || mbdcnd>4 )
ierror = 2;
end;
if( c>=d )
ierror = 3;
end;
if( n<=2 )
ierror = 4;
end;
if( nbdcnd<0 || nbdcnd>4 )
ierror = 5;
end;
if( idimf<m )
ierror = 7;
end;
if( m<=2 )
ierror = 8;
end;
if( ierror~=0 )
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
return;
end;
nperod = fix(nbdcnd);
mperod = 0;
if( mbdcnd>0 )
mperod = 1;
end;
deltax =(b-a)./m;
twdelx = 1../deltax;
delxsq = 2../deltax.^2;
deltay =(d-c)./n;
twdely = 1../deltay;
delysq = deltay.^2;
twdysq = 2../delysq;
np = fix(nbdcnd + 1);
mp = fix(mbdcnd + 1);
%
% DEFINE THE A,B,C COEFFICIENTS IN W-ARRAY.
%
id2 = fix(m);
id3 = fix(id2 + m);
id4 = fix(id3 + m);
s =(deltay./deltax).^2;
st2 = 2..*s;
for i = 1 : m;
w(i) = s;
j = fix(id2 + i);
w(j) = -st2 + elmbda.*delysq;
j = fix(id3 + i);
w(j) = s;
end; i = fix(m+1);
%
% ENTER BOUNDARY DATA FOR X-BOUNDARIES.
%
if( mp~=1 )
if( mp==4 || mp==5 )
for j = 1 : n;
f(1,j) = f(1,j) + bda(j).*twdelx;
end; j = fix(n+1);
w(id2+1) = w(id2+1) + w(1);
else;
for j = 1 : n;
f(1,j) = f(1,j) - bda(j).*delxsq;
end; j = fix(n+1);
w(id2+1) = w(id2+1) - w(1);
end;
if( mp~=1 )
if( mp==3 || mp==4 )
for j = 1 : n;
f(m,j) = f(m,j) - bdb(j).*twdelx;
end; j = fix(n+1);
w(id3) = w(id3) + w(1);
else;
for j = 1 : n;
f(m,j) = f(m,j) - bdb(j).*delxsq;
end; j = fix(n+1);
w(id3) = w(id3) - w(1);
end;
end;
end;
%
% ENTER BOUNDARY DATA FOR Y-BOUNDARIES.
%
if( np~=1 )
if( np==4 || np==5 )
for i = 1 : m;
f(i,1) = f(i,1) + bdc(i).*twdely;
end; i = fix(m+1);
else;
for i = 1 : m;
f(i,1) = f(i,1) - bdc(i).*twdysq;
end; i = fix(m+1);
end;
if( np~=1 )
if( np==3 || np==4 )
for i = 1 : m;
f(i,n) = f(i,n) - bdd(i).*twdely;
end; i = fix(m+1);
else;
for i = 1 : m;
f(i,n) = f(i,n) - bdd(i).*twdysq;
end; i = fix(m+1);
end;
end;
end;
for i = 1 : m;
for j = 1 : n;
f(i,j) = f(i,j).*delysq;
end; j = fix(n+1);
end; i = fix(m+1);
if( mperod~=0 )
w(1) = 0.;
w(id4) = 0.;
end;
pertrb = 0.;
if( elmbda>=0 )
if( elmbda~=0 )
ierror = 6;
elseif( mp~=2 && mp~=3 && mp~=5 ) ;
if( np~=2 && np~=3 && np~=5 )
%
% FOR SINGULAR PROBLEMS MUST ADJUST DATA TO INSURE THAT A SOLUTION
% WILL EXIST.
%
s = 0.;
for j = 1 : n;
for i = 1 : m;
s = s + f(i,j);
end; i = fix(m+1);
end; j = fix(n+1);
pertrb = s./(m.*n);
for j = 1 : n;
for i = 1 : m;
f(i,j) = f(i,j) - pertrb;
end; i = fix(m+1);
end; j = fix(n+1);
pertrb = pertrb./delysq;
end;
end;
end;
%
% SOLVE THE EQUATION.
%
if( nperod==0 )
[nperod,n,mperod,m,dumvar5,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=genbun(nperod,n,mperod,m,w(sub2ind(size(w),max(1,1)):end),w(sub2ind(size(w),max(id2+1,1)):end),w(sub2ind(size(w),max(id3+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(id4+1,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(1,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(id2+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(id3+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(id4+1,1)):end))~=(dumvar11)); w(1-1+dumvar5i)=dumvar5(dumvar5i); w(id2+1-1+dumvar6i)=dumvar6(dumvar6i); w(id3+1-1+dumvar7i)=dumvar7(dumvar7i); w(id4+1-1+dumvar11i)=dumvar11(dumvar11i);
else;
[nperod,n,mperod,m,dumvar5,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=poistg(nperod,n,mperod,m,w(sub2ind(size(w),max(1,1)):end),w(sub2ind(size(w),max(id2+1,1)):end),w(sub2ind(size(w),max(id3+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(id4+1,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(1,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(id2+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(id3+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(id4+1,1)):end))~=(dumvar11)); w(1-1+dumvar5i)=dumvar5(dumvar5i); w(id2+1-1+dumvar6i)=dumvar6(dumvar6i); w(id3+1-1+dumvar7i)=dumvar7(dumvar7i); w(id4+1-1+dumvar11i)=dumvar11(dumvar11i);
end;
w(1) = w(id4+1) + 3.*m;
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end %subroutine hstcrt
%DECK HSTCS1
|
|
