%  Two-dimensional boundary element program
%  with stress boundary conditions for Matlab 4.1
%  This program is modified from Appendix B of the book
%  "Boundary element methods in Solid Mechanics"
%  by S.L. Crouch and A.M. Starfield, 1983
%  Last revised on 2/8/94
	
	clear		%zero all variables

	global PR1 PR2 CON CONS
	global Sxxs Syys Sxys Sxxn Syyn Sxyn Uxs Uxn Uys Uyn

%  Read elastic constants, remote stress field data, and boundary data from 
%  data files elastic.dat, remote.dat, and boundary.dat, respectively

load elastic.dat;			%  PR and E should be in one row;
	PR = elastic(1);		%  PR = Poisson's ratio;
	E = elastic(2);			%  E = Young's modulus;
	clear elastic

load remote.dat				%  PXX, PYY, PXY should be in one row;
	Pxx = remote(1);		%  Remote stress sigma xx;
	Pyy = remote(2);		%  Remote stress sigma yy;
	Pxy = remote(3);		%  Remote stress sigma xy;
	clear remote

load boundary.dat
%  Data in each row of this file pertains to one element
%  Now separate out data types by column to form column vectors
	XBEG = boundary(:,1);	%  element endpoint coordinate
	YBEG = boundary(:,2);	%  element endpoint coordinate
	XEND = boundary(:,3);	%  element endpoint coordinate
	YEND = boundary(:,4);	%  element endpoint coordinate
	BVs = boundary(:,5);	%  Shear traction boundary condition
	BVn = boundary(:,6);	%  Normal traction boundary condition

	NUM = length(BVs);		% NUM = number of elements
	clear boundary

%  Define new elastic constants to be used internally
	CON = 1/(4*pi*(1-PR));
	CONS = E/(1+PR);
	PR1 = 1-2*PR;
	PR2 = 2*(1-PR);

%  Define locations, lengths, orientations, and boundary conditions for boundary elements
%  This information is stored here as column vectors of size (NUMx1)
	x = (XBEG + XEND)/2;				% element midpoints
	y = (YBEG + YEND)/2;				% element midpoints
	XE = x;
	YE = y;
	XD = XEND-XBEG;
	YD = YEND-YBEG;
	A = sqrt(XD.*XD +YD.*YD)/2;			% half-length of element
%  In the lines below, B = beta = orientation of element relative to global x-axis
	SINB = YD./(2*A);					% sine beta
	COSB = XD./(2*A);					% cosine beta
	SIN2B = 2*SINB.*COSB;				% sine 2*beta
	COS2B = COSB.*COSB-SINB.*SINB;		% cosine 2*beta
	SINB2 = SINB.*SINB;					% sine squared of beta
	COSB2 = COSB.*COSB;					% cosine squared of beta

%  Adjust stress boundary conditions to account for remote stresses.
%  First resolve components of the remote stress onto the plane of each element...
	SIGs = (Pyy-Pxx).*SIN2B/2 + Pxy.*COS2B;
	SIGn = Pxx.*SINB2 - Pxy.*SIN2B + Pyy.*COSB2;
%  and then subtract the resolved remote stress from the boundary stresses.
	BVs = BVs - SIGs;
	BVn = BVn - SIGn;
%  The remote stress will be added back at the end of the solution.

%  Compute influence coefficients between boundary elements.
%  C(i,j) = effect at obs pt i due to a load at element j.  
%  The organization of the coefficients here is slightly different 
%  from Crouch & Starfield
%	First dimension the vectors storing the influence coefficients
	nobs=NUM*NUM;
	set1;			
	for i=1:NUM
		first = (i-1)*NUM+1;
		last = i*NUM;
		xe = XE(i);
		ye = YE(i);
		COSBe = COSB(i);
		SINBe = SINB(i);
		a = A(i);
		Sd = 1;
		Nd = 1;
		coeff;
		% Append sequentially produced vectors to master column vectors 
		append1;		 
	end

% Calculate trig functions
% and direction cosine matrices for stress transformations
	dircos1

%  Reshape column vectors into square matrices, then clear column vectors
	dimx = NUM;
	dimy = NUM;
	square1;

%  Convert xy stresses returned from coeff to normal and shear tractions on elements
	Cisjs = (Syys-Sxxs).*(SINN2B)/2 + Sxys.*COSS2B;
	Cisjn = (Syyn-Sxxn).*(SINN2B)/2 + Sxyn.*COSS2B;
	Cinjs = Sxxs.*SINNB2 - Sxys.*SINN2B + Syys.*COSSB2;
	Cinjn = Sxxn.*SINNB2 - Sxyn.*SINN2B + Syyn.*COSSB2;

%  Solve system of algebraic equations to get the displacement discontinuities
%  First regroup the influence coefficient and boundary condition submatrices as shown:
%  |   [Cisjs]	[Cisjn]	|	| [Ds]	|	=	| [Bs]	|
%  |   [Cinjs]	[Cinjn] |	| [Dn]	|		| [Bn]	|
	C = [Cisjs, Cisjn; Cinjs, Cinjn];		% C = influence coefficients;
%	clear Cisjs Cisjn Cinjs Cinjn;
	B = [BVs; BVn];							% B= Stress boundary conditions;
%  Solve the matrix system [C][D]=[B] to get the displacement discontinuity vector D
	D= C\B;
%  Separate the D vector into subvectors Ds and Dn
	Ds = D(1:NUM);				% Ds elements are in upper half of D
	Dn = D(NUM+1:2*NUM);		% Dn elements are in lower half of D
%	clear B C D

%  Compute displacements and stresses on boundary elements
%  The stresses obtained should match the boundary conditions

%  First find the stresses and displacements, in an x-y reference frame, by superposition
	UxNEG = Uxs*Ds + Uxn*Dn;			% x-displacement on negative side of i;
	UyNEG = Uys*Ds + Uyn*Dn;			% y-displacement on negative side of i;
	SIGxx = Pxx + Sxxs*Ds + Sxxn*Dn;	% Note that remote stress is added back in;
	SIGyy = Pyy + Syys*Ds + Syyn*Dn;	% Note that remote stress is added back in;
	SIGxy = Pxy + Sxys*Ds + Sxyn*Dn;	% Note that remote stress is added back in;

%  Now resolve DISPLACEMENTS in the xy frame into shear and normal components
%  Start by solving for displacements on the NEGATIVE sides of the elements...;
	UsNEG = UxNEG.*COSB + UyNEG.*SINB;
	UnNEG = -UxNEG.*SINB + UyNEG.*COSB;

%  and then add the displacement discontinuity at each element to get 
%  the displacement components on the positive sides of the elements
	UsPOS = UsNEG - Ds;
	UnPOS = UnNEG - Dn;

%  The last step is to resolve the xy STRESSES to normal and shear stresses on the elements.
	SIGs = (SIGyy-SIGxx).*(SIN2B/2)  +  SIGxy.*COS2B;
	SIGn = SIGxx.*SINB2  -  SIGxy.*SIN2B  +  SIGyy.*COSB2;
%  Any writing to output files should occur now....
%	profile		% To plot profiles of displacements along boundaries

%  Clear arrays in preparation for next section, saving vectors Ds and Dn 
%	clear xi yi xm ym
%	clear Uxs Uxn Uys Uyn Sxxs Sxxn Syys Syyn Sxys Sxyn;
%	clear SIGxx SIGyy SIGxy 
%	clear SIGs SIGn
%	clear UxPOS UxNEG UyPOS UyNEG UsPOS UsNEG UnPOS UnNEG
%	clear SINNB COSSB SINN2B COSS2B SINNB2 COSSB2 

%  COMPUTATION OF DISPLACEMENTS AND STRESSES AT SPECIFIED POINTS I IN BODY
check = 0;
if check ==1
%  Begin by defining the coordinates of the gridpoints
	twoddgrid

%  Calculate influence coefficients at gridpoints
	x = X(:);			% Grid nodes in column vector form
	y = Y(:);			% Grid nodes in column vector form
%	First dimension the vectors storing the influence coefficients
	nobs=NUM2*NUM;
	set1;
% Now loop through grid element by element
% The end result will be ten column vectors each having NUM2*NUM rows			
	for i=1:NUM
		first = (i-1)*NUM2+1;
		last = i*NUM2;
		xe = XE(i);
		ye = YE(i);
		COSBe = COSB(i);
		SINBe = SINB(i);
		a = A(i);
		Sd = Ds(i);
		Nd = Dn(i);
		coeff;
		% Append sequentially produced vectors to master column vectors 
		append1;		 
	end

% Reshape column vectors into square matrices the size of the
% observation grid, and in the process sum the contributions
% from each element 
	square2;

% Add back the remote field
	Sxx = Pxx + Sxx;
	Syy = Pyy + Syy;
	Sxy = Pxy + Sxy;

	
S1 = SIG1(Sxx,Syy,Sxy);
S2 = SIG2(Sxx,Syy,Sxy);
tau = taumax(Sxx,Syy,Sxy);
mean = ave(Sxx,Syy,Sxy);
theta = angp(Sxx,Syy,Sxy); 
end