| [intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,grhs,usol,idmn,w,pertrb,ierror]=sepeli(intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,grhs,usol,idmn,w,pertrb,ierror); |
function [intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,grhs,usol,idmn,w,pertrb,ierror]=sepeli(intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,grhs,usol,idmn,w,pertrb,ierror);
persistent i1 i10 i11 i12 i13 i2 i3 i4 i5 i6 i7 i8 i9 k l length linput ll logb2n loutpt ;
if isempty(i1), i1=0; end;
if isempty(i10), i10=0; end;
if isempty(i11), i11=0; end;
if isempty(i12), i12=0; end;
if isempty(i13), i13=0; end;
if isempty(i2), i2=0; end;
if isempty(i3), i3=0; end;
if isempty(i4), i4=0; end;
if isempty(i5), i5=0; end;
if isempty(i6), i6=0; end;
if isempty(i7), i7=0; end;
if isempty(i8), i8=0; end;
if isempty(i9), i9=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(length), length=0; end;
%***BEGIN PROLOGUE SEPELI
%***PURPOSE Discretize and solve a second and, optionally, a fourth
% order finite difference approximation on a uniform grid to
% the general separable elliptic partial differential
% equation on a rectangle with any combination of periodic or
% mixed boundary conditions.
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B1A2
%***TYPE SINGLE PRECISION (SEPELI-S)
%***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, SEPARABLE
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% Dimension of BDA(N+1), BDB(N+1), BDC(M+1), BDD(M+1),
% Arguments USOL(IDMN,N+1), GRHS(IDMN,N+1),
% W (see argument list)
%
% Latest Revision March 1977
%
% Purpose SEPELI solves for either the second-order
% finite difference approximation or a
% fourth-order approximation to a separable
% elliptic equation.
%
% 2 2
% AF(X)*d U/dX + BF(X)*dU/dX + CF(X)*U +
% 2 2
% DF(Y)*d U/dY + EF(Y)*dU/dY + FF(Y)*U
%
% = G(X,Y)
%
% on a rectangle (X greater than or equal to A
% and less than or equal to B; Y greater than
% or equal to C and less than or equal to D).
% Any combination of periodic or mixed boundary
% conditions is allowed.
%
% Purpose The possible boundary conditions are:
% in the X-direction:
% (0) Periodic, U(X+B-A,Y)=U(X,Y) for all Y,X
% (1) U(A,Y), U(B,Y) are specified for all Y
% (2) U(A,Y), dU(B,Y)/dX+BETA*U(B,Y) are
% specified for all Y
% (3) dU(A,Y)/dX+ALPHA*U(A,Y),dU(B,Y)/dX+
% BETA*U(B,Y) are specified for all Y
% (4) dU(A,Y)/dX+ALPHA*U(A,Y),U(B,Y) are
% specified for all Y
%
% in the Y-direction:
% (0) Periodic, U(X,Y+D-C)=U(X,Y) for all X,Y
% (1) U(X,C),U(X,D) are specified for all X
% (2) U(X,C),dU(X,D)/dY+XNU*U(X,D) are specified
% for all X
% (3) dU(X,C)/dY+GAMA*U(X,C),dU(X,D)/dY+
% XNU*U(X,D) are specified for all X
% (4) dU(X,C)/dY+GAMA*U(X,C),U(X,D) are
% specified for all X
%
% Arguments
%
% On Input INTL
% = 0 On initial entry to SEPELI or if any of
% the arguments C, D, N, NBDCND, COFY are
% changed from a previous call
% = 1 If C, D, N, NBDCND, COFY are unchanged
% from the previous call.
%
% IORDER
% = 2 If a second-order approximation is sought
% = 4 If a fourth-order approximation is sought
%
% A,B
% The range of the X-independent variable;
% i.e., X is greater than or equal to A and
% less than or equal to B. A must be less than
% B.
%
% M
% The number of panels into which the interval
% [A,B] is subdivided. Hence, there will be
% M+1 grid points in the X-direction given by
% XI=A+(I-1)*DLX for I=1,2,...,M+1 where
% DLX=(B-A)/M is the panel width. M must be
% less than IDMN and greater than 5.
%
% MBDCND
% Indicates the type of boundary condition at
% X=A and X=B
% = 0 If the solution is periodic in X; i.e.,
% U(X+B-A,Y)=U(X,Y) for all Y,X
% = 1 If the solution is specified at X=A and
% X=B; i.e., U(A,Y) and U(B,Y) are
% specified for all Y
% = 2 If the solution is specified at X=A and
% the boundary condition is mixed at X=B;
% i.e., U(A,Y) and dU(B,Y)/dX+BETA*U(B,Y)
% are specified for all Y
% = 3 If the boundary conditions at X=A and X=B
% are mixed; i.e., dU(A,Y)/dX+ALPHA*U(A,Y)
% and dU(B,Y)/dX+BETA*U(B,Y) are specified
% for all Y
% = 4 If the boundary condition at X=A is mixed
% and the solution is specified at X=B;
% i.e., dU(A,Y)/dX+ALPHA*U(A,Y) and U(B,Y)
% are specified for all Y
%
% BDA
% A one-dimensional array of length N+1 that
% specifies the values of dU(A,Y)/dX+
% ALPHA*U(A,Y) at X=A, when MBDCND=3 or 4.
% BDA(J) = dU(A,YJ)/dX+ALPHA*U(A,YJ);
% J=1,2,...,N+1
% when MBDCND has any other value, BDA is a
% dummy parameter.
%
% On Input ALPHA
% The scalar multiplying the solution in case
% of a mixed boundary condition at X=A (see
% argument BDA). If MBDCND = 3,4 then ALPHA is
% a dummy parameter.
%
% BDB
% A one-dimensional array of length N+1 that
% specifies the values of dU(B,Y)/dX+
% BETA*U(B,Y) at X=B. When MBDCND=2 or 3
% BDB(J) = dU(B,YJ)/dX+BETA*U(B,YJ);
% J=1,2,...,N+1
% When MBDCND has any other value, BDB is a
% dummy parameter.
%
% BETA
% The scalar multiplying the solution in case
% of a mixed boundary condition at X=B (see
% argument BDB). If MBDCND=2,3 then BETA is a
% dummy parameter.
%
% C,D
% The range of the Y-independent variable;
% i.e., Y is greater than or equal to C and
% less than or equal to D. C must be less than
% D.
%
% N
% The number of panels into which the interval
% [C,D] is subdivided. Hence, there will be
% N+1 grid points in the Y-direction given by
% YJ=C+(J-1)*DLY for J=1,2,...,N+1 where
% DLY=(D-C)/N is the panel width. In addition,
% N must be greater than 4.
%
% NBDCND
% Indicates the types of boundary conditions at
% Y=C and Y=D
% = 0 If the solution is periodic in Y; i.e.,
% U(X,Y+D-C)=U(X,Y) for all X,Y
% = 1 If the solution is specified at Y=C and
% Y = D, i.e., U(X,C) and U(X,D) are
% specified for all X
% = 2 If the solution is specified at Y=C and
% the boundary condition is mixed at Y=D;
% i.e., U(X,C) and dU(X,D)/dY+XNU*U(X,D)
% are specified for all X
% = 3 If the boundary conditions are mixed at
% Y=C and Y=D; i.e., dU(X,D)/dY+GAMA*U(X,C)
% and dU(X,D)/dY+XNU*U(X,D) are specified
% for all X
% = 4 If the boundary condition is mixed at Y=C
% and the solution is specified at Y=D;
% i.e. dU(X,C)/dY+GAMA*U(X,C) and U(X,D)
% are specified for all X
%
% BDC
% A one-dimensional array of length M+1 that
% specifies the value of dU(X,C)/dY+GAMA*U(X,C)
% at Y=C. When NBDCND=3 or 4
% BDC(I) = dU(XI,C)/dY + GAMA*U(XI,C);
% I=1,2,...,M+1.
% When NBDCND has any other value, BDC is a
% dummy parameter.
%
% GAMA
% The scalar multiplying the solution in case
% of a mixed boundary condition at Y=C (see
% argument BDC). If NBDCND=3,4 then GAMA is a
% dummy parameter.
%
% BDD
% A one-dimensional array of length M+1 that
% specifies the value of dU(X,D)/dY +
% XNU*U(X,D) at Y=C. When NBDCND=2 or 3
% BDD(I) = dU(XI,D)/dY + XNU*U(XI,D);
% I=1,2,...,M+1.
% When NBDCND has any other value, BDD is a
% dummy parameter.
%
% XNU
% The scalar multiplying the solution in case
% of a mixed boundary condition at Y=D (see
% argument BDD). If NBDCND=2 or 3 then XNU is
% a dummy parameter.
%
% COFX
% A user-supplied subprogram with
% parameters X, AFUN, BFUN, CFUN which
% returns the values of the X-dependent
% coefficients AF(X), BF(X), CF(X) in
% the elliptic equation at X.
%
% COFY
% A user-supplied subprogram with
% parameters Y, DFUN, EFUN, FFUN which
% returns the values of the Y-dependent
% coefficients DF(Y), EF(Y), FF(Y) in
% the elliptic equation at Y.
%
% NOTE: COFX and COFY must be declared external
% in the calling routine. The values returned in
% AFUN and DFUN must satisfy AFUN*DFUN greater
% than 0 for A less than X less than B,
% C less than Y less than D (see IERROR=10).
% The coefficients provided may lead to a matrix
% equation which is not diagonally dominant in
% which case solution may fail (see IERROR=4).
%
% GRHS
% A two-dimensional array that specifies the
% values of the right-hand side of the elliptic
% equation; i.e., GRHS(I,J)=G(XI,YI), for
% I=2,...,M; J=2,...,N. At the boundaries,
% GRHS is defined by
%
% MBDCND GRHS(1,J) GRHS(M+1,J)
% ------ --------- -----------
% 0 G(A,YJ) G(B,YJ)
% 1 * *
% 2 * G(B,YJ) J=1,2,...,N+1
% 3 G(A,YJ) G(B,YJ)
% 4 G(A,YJ) *
%
% NBDCND GRHS(I,1) GRHS(I,N+1)
% ------ --------- -----------
% 0 G(XI,C) G(XI,D)
% 1 * *
% 2 * G(XI,D) I=1,2,...,M+1
% 3 G(XI,C) G(XI,D)
% 4 G(XI,C) *
%
% where * means these quantities are not used.
% GRHS should be dimensioned IDMN by at least
% N+1 in the calling routine.
%
% USOL
% A two-dimensional array that specifies the
% values of the solution along the boundaries.
% At the boundaries, USOL is defined by
%
% MBDCND USOL(1,J) USOL(M+1,J)
% ------ --------- -----------
% 0 * *
% 1 U(A,YJ) U(B,YJ)
% 2 U(A,YJ) * J=1,2,...,N+1
% 3 * *
% 4 * U(B,YJ)
%
% NBDCND USOL(I,1) USOL(I,N+1)
% ------ --------- -----------
% 0 * *
% 1 U(XI,C) U(XI,D)
% 2 U(XI,C) * I=1,2,...,M+1
% 3 * *
% 4 * U(XI,D)
%
% where * means the quantities are not used in
% the solution.
%
% If IORDER=2, the user may equivalence GRHS
% and USOL to savemlv space. Note that in this
% case the tables specifying the boundaries of
% the GRHS and USOL arrays determine the
% boundaries uniquely except at the corners.
% If the tables call for both G(X,Y) and
% U(X,Y) at a corner then the solution must be
% chosen. For example, if MBDCND=2 and
% NBDCND=4, then U(A,C), U(A,D), U(B,D) must be
% chosen at the corners in addition to G(B,C).
%
% If IORDER=4, then the two arrays, USOL and
% GRHS, must be distinct.
%
% USOL should be dimensioned IDMN by at least
% N+1 in the calling routine.
%
% IDMN
% The row (or first) dimension of the arrays
% GRHS and USOL as it appears in the program
% calling SEPELI. This parameter is used to
% specify the variable dimension of GRHS and
% USOL. IDMN must be at least 7 and greater
% than or equal to M+1.
%
% W
% A one-dimensional array that must be provided
% by the user for work space. Let
% K=INT(log2(N+1))+1 and set L=2**(K+1).
% then (K-2)*L+K+10*N+12*M+27 will suffice
% as a length of W. THE actual length of W in
% the calling routine must be set in W(1) (see
% IERROR=11).
%
% On Output USOL
% contains the approximate solution to the
% elliptic equation. USOL(I,J) is the
% approximation to U(XI,YJ) for I=1,2...,M+1
% and J=1,2,...,N+1. The approximation has
% error O(DLX**2+DLY**2) if called with
% IORDER=2 and O(DLX**4+DLY**4) if called with
% IORDER=4.
%
% W
% contains intermediate values that must not be
% destroyed if SEPELI is called again with
% INTL=1. In addition W(1) contains the exact
% minimal length (in floating point) required
% for the work space (see IERROR=11).
%
% PERTRB
% If a combination of periodic or derivative
% boundary conditions (i.e., ALPHA=BETA=0 if
% MBDCND=3; GAMA=XNU=0 if NBDCND=3) is
% specified and if the coefficients of U(X,Y)
% in the separable elliptic equation are zero
% (i.e., CF(X)=0 for X greater than or equal to
% A and less than or equal to B; FF(Y)=0 for
% Y greater than or equal to C and less than
% or equal to D) then a solution may not exist.
% PERTRB is a constant calculated and
% subtracted from the right-hand side of the
% matrix equations generated by SEPELI which
% insures that a solution exists. SEPELI then
% computes this solution which is a weighted
% minimal least squares solution to the
% original problem.
%
% IERROR
% An error flag that indicates invalid input
% parameters or failure to find a solution
% = 0 No error
% = 1 If A greater than B or C greater than D
% = 2 If MBDCND less than 0 or MBDCND greater
% than 4
% = 3 If NBDCND less than 0 or NBDCND greater
% than 4
% = 4 If attempt to find a solution fails.
% (the linear system generated is not
% diagonally dominant.)
% = 5 If IDMN is too small (see discussion of
% IDMN)
% = 6 If M is too small or too large (see
% discussion of M)
% = 7 If N is too small (see discussion of N)
% = 8 If IORDER is not 2 or 4
% = 9 If INTL is not 0 or 1
% = 10 If AFUN*DFUN less than or equal to 0 for
% some interior mesh point (XI,YJ)
% = 11 If the work space length input in W(1)
% is less than the exact minimal work
% space length required output in W(1).
%
% NOTE (concerning IERROR=4): for the
% coefficients input through COFX, COFY, the
% discretization may lead to a block
% tridiagonal linear system which is not
% diagonally dominant (for example, this
% happens if CFUN=0 and BFUN/(2.*DLX) greater
% than AFUN/DLX**2). In this case solution may
% fail. This cannot happen in the limit as
% DLX, DLY approach zero. Hence, the condition
% may be remedied by taking larger values for M
% or N.
%
% Entry Points SEPELI, SPELIP, CHKPRM, CHKSNG, ORTHOG, MINSOL,
% TRISP, DEFER, DX, DY, BLKTRI, BLKTR1, INDXB,
% INDXA, INDXC, PROD, PRODP, CPROD, CPRODP,
% PPADD, PSGF, BSRH, PPSGF, PPSPF, COMPB,
% TRUN1, STOR1, TQLRAT
%
% Special Conditions NONE
%
% Common Blocks SPLP, CBLKT
%
% I/O NONE
%
% Precision Single
%
% Specialist John C. Adams, NCAR, Boulder, Colorado 80307
%
% Language FORTRAN
%
% History Developed at NCAR during 1975-76.
if isempty(linput), linput=0; end;
if isempty(ll), ll=0; end;
if isempty(logb2n), logb2n=0; end;
if isempty(loutpt), loutpt=0; end;
%
% Algorithm SEPELI automatically discretizes the separable
% elliptic equation which is then solved by a
% generalized cyclic reduction algorithm in the
% subroutine, BLKTRI. The fourth-order solution
% is obtained using 'Deferred Corrections' which
% is described and referenced in sections,
% references and method.
%
% Space Required 14654 (octal) = 6572 (decimal)
%
% Accuracy and Timing The following computational results were
% obtained by solving the sample problem at the
% end of this write-up on the Control Data 7600.
% The op count is proportional to M*N*log2(N).
% In contrast to the other routines in this
% chapter, accuracy is tested by computing and
% tabulating second- and fourth-order
% discretization errors. Below is a table
% containing computational results. The times
% given do not include initialization (i.e.,
% times are for INTL=1). Note that the
% fourth-order accuracy is not realized until the
% mesh is sufficiently refined.
%
% Second-order Fourth-order Second-order Fourth-order
% M N Execution Time Execution Time Error Error
% (M SEC) (M SEC)
% 6 6 6 14 6.8E-1 1.2E0
% 14 14 23 58 1.4E-1 1.8E-1
% 30 30 100 247 3.2E-2 9.7E-3
% 62 62 445 1,091 7.5E-3 3.0E-4
% 126 126 2,002 4,772 1.8E-3 3.5E-6
%
% Portability There are no machine-dependent constants.
%
% Required Resident SQRT, ABS, LOG
% Routines
%
% References Keller, H.B., 'Numerical Methods for Two-point
% Boundary-value Problems', Blaisdel (1968),
% Waltham, Mass.
%
% Swarztrauber, P., and R. Sweet (1975):
% 'Efficient FORTRAN Subprograms for The
% Solution of Elliptic Partial Differential
% Equations'. NCAR Technical Note
% NCAR-TN/IA-109, pp. 135-137.
%
%***REFERENCES H. B. Keller, Numerical Methods for Two-point
% Boundary-value Problems, Blaisdel, Waltham, Mass.,
% 1968.
% P. N. Swarztrauber and R. Sweet, Efficient Fortran
% subprograms for the solution of elliptic equations,
% NCAR TN/IA-109, July 1975, 138 pp.
%***ROUTINES CALLED CHKPRM, SPELIP
%***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 SEPELI
%
grhs_shape=size(grhs);grhs=reshape([grhs(:).',zeros(1,ceil(numel(grhs)./prod([idmn])).*prod([idmn])-numel(grhs))],idmn,[]);
usol_shape=size(usol);usol=reshape([usol(:).',zeros(1,ceil(numel(usol)./prod([idmn])).*prod([idmn])-numel(usol))],idmn,[]);
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 SEPELI
[intl,iorder,a,b,m,mbdcnd,c,d,n,nbdcnd,cofx,cofy,idmn,ierror]=chkprm(intl,iorder,a,b,m,mbdcnd,c,d,n,nbdcnd,cofx,cofy,idmn,ierror);
if( ierror~=0 )
grhs_shape=zeros(grhs_shape);grhs_shape(:)=grhs(1:numel(grhs_shape));grhs=grhs_shape;
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_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;
%
% COMPUTE MINIMUM WORK SPACE AND CHECK WORK SPACE LENGTH INPUT
%
l = fix(n + 1);
if( nbdcnd==0 )
l = fix(n);
end;
logb2n = fix(fix(log(l+0.5)./log(2.0)) + 1);
ll = fix(2.^(logb2n+1));
k = fix(m + 1);
l = fix(n + 1);
length =fix((logb2n-2).*ll + logb2n + max(2.*l,6.*k) + 5);
if( nbdcnd==0 )
length = fix(length + 2.*l);
end;
ierror = 11;
linput = fix(fix(w(1)+0.5));
loutpt = fix(length + 6.*(k+l) + 1);
w(1) = loutpt;
if( loutpt>linput )
grhs_shape=zeros(grhs_shape);grhs_shape(:)=grhs(1:numel(grhs_shape));grhs=grhs_shape;
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_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;
ierror = 0;
%
% SET WORK SPACE INDICES
%
i1 = fix(length + 2);
i2 = fix(i1 + l);
i3 = fix(i2 + l);
i4 = fix(i3 + l);
i5 = fix(i4 + l);
i6 = fix(i5 + l);
i7 = fix(i6 + l);
i8 = fix(i7 + k);
i9 = fix(i8 + k);
i10 = fix(i9 + k);
i11 = fix(i10 + k);
i12 = fix(i11 + k);
i13 = 2;
[intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,dumvar21,dumvar22,dumvar23,dumvar24,dumvar25,dumvar26,dumvar27,dumvar28,dumvar29,dumvar30,dumvar31,dumvar32,grhs,usol,idmn,dumvar36,pertrb,ierror]=spelip(intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,w(sub2ind(size(w),max(i1,1)):end),w(sub2ind(size(w),max(i2,1)):end),w(sub2ind(size(w),max(i3,1)):end),w(sub2ind(size(w),max(i4,1)):end),w(sub2ind(size(w),max(i5,1)):end),w(sub2ind(size(w),max(i6,1)):end),w(sub2ind(size(w),max(i7,1)):end),w(sub2ind(size(w),max(i8,1)):end),w(sub2ind(size(w),max(i9,1)):end),w(sub2ind(size(w),max(i10,1)):end),w(sub2ind(size(w),max(i11,1)):end),w(sub2ind(size(w),max(i12,1)):end),grhs,usol,idmn,w(sub2ind(size(w),max(i13,1)):end),pertrb,ierror); dumvar21i=find((w(sub2ind(size(w),max(i1,1)):end))~=(dumvar21));dumvar22i=find((w(sub2ind(size(w),max(i2,1)):end))~=(dumvar22));dumvar23i=find((w(sub2ind(size(w),max(i3,1)):end))~=(dumvar23));dumvar24i=find((w(sub2ind(size(w),max(i4,1)):end))~=(dumvar24));dumvar25i=find((w(sub2ind(size(w),max(i5,1)):end))~=(dumvar25));dumvar26i=find((w(sub2ind(size(w),max(i6,1)):end))~=(dumvar26));dumvar27i=find((w(sub2ind(size(w),max(i7,1)):end))~=(dumvar27));dumvar28i=find((w(sub2ind(size(w),max(i8,1)):end))~=(dumvar28));dumvar29i=find((w(sub2ind(size(w),max(i9,1)):end))~=(dumvar29));dumvar30i=find((w(sub2ind(size(w),max(i10,1)):end))~=(dumvar30));dumvar31i=find((w(sub2ind(size(w),max(i11,1)):end))~=(dumvar31));dumvar32i=find((w(sub2ind(size(w),max(i12,1)):end))~=(dumvar32));dumvar36i=find((w(sub2ind(size(w),max(i13,1)):end))~=(dumvar36)); w(i1-1+dumvar21i)=dumvar21(dumvar21i); w(i2-1+dumvar22i)=dumvar22(dumvar22i); w(i3-1+dumvar23i)=dumvar23(dumvar23i); w(i4-1+dumvar24i)=dumvar24(dumvar24i); w(i5-1+dumvar25i)=dumvar25(dumvar25i); w(i6-1+dumvar26i)=dumvar26(dumvar26i); w(i7-1+dumvar27i)=dumvar27(dumvar27i); w(i8-1+dumvar28i)=dumvar28(dumvar28i); w(i9-1+dumvar29i)=dumvar29(dumvar29i); w(i10-1+dumvar30i)=dumvar30(dumvar30i); w(i11-1+dumvar31i)=dumvar31(dumvar31i); w(i12-1+dumvar32i)=dumvar32(dumvar32i); w(i13-1+dumvar36i)=dumvar36(dumvar36i);
grhs_shape=zeros(grhs_shape);grhs_shape(:)=grhs(1:numel(grhs_shape));grhs=grhs_shape;
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_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
%DECK SEPX4
|
|
