function [C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy,...
Ux,Uy,Sxx,Syy,Sxy,S1,S2,tau,mean,theta] = twodd_func(e_dat,r_dat,b_dat,X,Y)
%  Function version of Two-dimensional boundary element program
%  All the functions twodd_func calls are internally stored subfunctions
%  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 6/29/03.  By Stephen J. Martel.
% Examples:
%load e_dat
%load r_dat
%load b_dat
%xmin = -1.9;
%xmax = +1.9;
%dx = 0.2;
%ymin = -2.9;
%ymax = +2.9;
%dy = 0.2;
%xgrid = xmin:dx:xmax;
%ygrid = ymin:dy:ymax;
%[X,Y] = meshgrid(xgrid,ygrid);
%[C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy]= twodd_func(e_dat,r_dat,b_dat)
%[C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy,...
%Ux,Uy,Sxx,Syy,Sxy,S1,S2,tau,mean,theta] = twodd_func(e_dat,r_dat,b_dat,X,Y)


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

% Elastic constants.  PR and E should be in one row;
	PR = e_dat(1);			%  PR = Poisson's ratio;
	E = e_dat(2);			%  E = Young's modulus;

% Ambient field.  PXX, PYY, PXY should be in one row;
	Pxx = r_dat(1);		%  Remote stress sigma xx;
	Pyy = r_dat(2);		%  Remote stress sigma yy;
	Pxy = r_dat(3);		%  Remote stress sigma xy;

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

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

%  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;
	[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs);			
	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;
		[Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
		coeff_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,Sd,Nd);
		% Append sequentially produced vectors to master column vectors 
		[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
		append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
		Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv);
	end

% Calculate trig functions
% and direction cosine matrices for stress transformations
	[SINNB,COSSB,SINN2B,COSS2B,SINNB2,COSSB2] = ...
	dircos1_func(SINB,COSB,SIN2B,COS2B,SINB2,COSB2,NUM);

%  Reshape column vectors into square matrices, then clear column vectors
	dimx = NUM;
	dimy = NUM;
	[Uxs,Uxn,Uys,Uyn,Sxxs,Sxxn,Syys,Syyn,Sxys,Sxyn] = ...
	square1_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,dimx,dimy);

%  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
	UxN = Uxs*Ds + Uxn*Dn;			% x-displacement on negative side of i;
	UyN = 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 (N) sides of the elements...;
	UsN = UxN.*COSB + UyN.*SINB;
	UnN = -UxN.*SINB + UyN.*COSB;

%  and then add the displacement discontinuity at each element to get 
%  the displacement components on the positive (P) sides of the elements
	UsP = UsN - Ds;
	UnP = UnN - 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;

%  COMPUTATION OF DISPLACEMENTS AND STRESSES AT SPECIFIED OBSERVATION POINTS I IN BODY
%  Proceed if output on observation grid X,Y is desired, otherwise end
if nargout > 13
	[dimx,dimy] = size(X);
	NUM2 = dimx*dimy;
%  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;
	[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs);
% 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);
		[Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
		coeff_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,Sd,Nd);
		% Append sequentially produced vectors to master column vectors 
		[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
		append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
		Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv);
	end

% Reshape column vectors into square matrices the size of the
% observation grid, and in the process sum the contributions
% from each element 
	[Ux,Uy,Sxx,Syy,Sxy] = ...
square2_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,NUM2,NUM,dimx,dimy);

% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Alphabetical listing of local subfunctions called in the primary function twodd_func 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function theta=angp(SIGxx,SIGyy,SIGxy)
% function angp.  Calculates a principal stress orientation
theta = 0.5*atan2(SIGxy,(SIGxx-SIGyy)/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv)
% Function to save influence coeffients from twodd.m in column vectors
% by appending new column vectors generated by the first set of
% calls to coeff to the previously generated column vectors
Uxsv (first:last) = Uxs;
Uysv (first:last) = Uys;
Uxnv (first:last) = Uxn;
Uynv (first:last) = Uyn;
Sxxsv (first:last) = Sxxs;
Syysv (first:last) = Syys;
Sxysv (first:last) = Sxys;
Sxxnv (first:last) = Sxxn;
Syynv (first:last) = Syyn;
Sxynv (first:last) = Sxyn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mean=ave(SIGxx,SIGyy,SIGxy)
mean = (SIGxx+SIGyy)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
coeff_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,Sd,Nd)
%  This is the MATLAB function version of SUBROUTINE COEFF from TWODD
%  It calculates the xy stress and displacement components
%  at all points (x,y) in an observation array
%  due to unit normal and unit shear displacement 
%  discontinuities at element j(xe,ye).

% Define trig functions and coordinates pertaining to the element
% and create arrays of the same size as the observation array
unit = ones(size(x));
COS2Be = (COSBe.*COSBe-SINBe.*SINBe).*unit;
SIN2Be = (2.*SINBe.*COSBe).*unit;
COSB2e = (COSBe.*COSBe).*unit;
SINB2e = (SINBe.*SINBe).*unit;
xe = xe.*unit;
ye = ye.*unit;	

% Calculate coordinates of observation points in reference frame
% centered at the element and oriented along the element.
% The transformation requires a translation and a rotation
Xb = (x-xe).*COSBe + (y-ye).*SINBe;
Yb = -(x-xe).*SINBe + (y-ye).*COSBe;
% Flag all observation points for which Yb is not 0
i1 = find(Yb);
% Flag any observation points on the element
i2 = find(Yb == 0 & abs(Xb) < a);

% Find squares of the distances from obs. pts. to element ends
R1S = (Xb-a).*(Xb-a) + Yb.*Yb;
R2S = (Xb+a).*(Xb+a) + Yb.*Yb;

%  Calculate intermediate functions F
FL1 = 0.5.*log(R1S);
FL2 = 0.5.*log(R2S);
FB2 = CON.*(FL1-FL2);
% FB3 = 0 for pts colinear with element, CON*pi for pts. on element 
% FB3 = difference of arc tangents for all other pts.
FB3 = zeros(size(x));				
FB3(i1) =  -1.*CON.*(atan((Xb(i1)+a)./Yb(i1)) - atan((Xb(i1)-a)./Yb(i1)));
FB3(i2) = CON.*pi.*ones(size(i2));
FB4 = CON.*(Yb./R1S - Yb./R2S);
FB5 = CON.*((Xb-a)./R1S - (Xb+a)./R2S);
FB6 = CON.*(((Xb-a).^2 -Yb.*Yb)./R1S.^2 - ((Xb+a).^2-Yb.*Yb)./R2S.^2);
FB7 = 2.*CON.*Yb.*((Xb-a)./R1S.^2 - (Xb+a)./R2S.^2);

%  Calculate DISPLACEMENT components due to unit SHEAR disp. discontinuity
Uxs = -PR1.*SINBe.*FB2 + PR2.*COSBe.*FB3 + Yb.*(SINBe.*FB4-COSBe.*FB5);
Uys =  PR1.*COSBe.*FB2 + PR2.*SINBe.*FB3 - Yb.*(COSBe.*FB4+SINBe.*FB5);

%  Calculate DISPLACEMENT components due to unit NORMAL disp. discontinuity
Uxn = -PR1.*COSBe.*FB2 - PR2.*SINBe.*FB3 - Yb.*(COSBe.*FB4+SINBe.*FB5);
Uyn = -PR1.*SINBe.*FB2 + PR2.*COSBe.*FB3 - Yb.*(SINBe.*FB4-COSBe.*FB5);

%  Calculate STRESS components due to unit SHEAR disp. discontinuity
Sxxs = CONS.*(2.*COSB2e.*FB4 + SIN2Be.*FB5 + Yb.*(COS2Be.*FB6-SIN2Be.*FB7));
Syys = CONS.*(2.*SINB2e.*FB4 - SIN2Be.*FB5 - Yb.*(COS2Be.*FB6-SIN2Be.*FB7));
Sxys = CONS.*(SIN2Be.*FB4 - COS2Be.*FB5 + Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));

%  Calculate STRESS components due to unit NORMAL displacement discontinuity
Sxxn = CONS.*(-FB5 + Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));
Syyn = CONS.*(-FB5 - Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));
Sxyn = CONS.*(-Yb.*(COS2Be.*FB6 - SIN2Be.*FB7));

% Multiply the components by the shear and normal displacement 
% discontintuities.  For the first pass through coeff, Sd=Nd=1. 
Uxs =  Sd*Uxs;
Uys =  Sd*Uys;
Uxn =  Nd*Uxn;
Uyn =  Nd*Uyn;
Sxxs = Sd*Sxxs;
Syys = Sd*Syys;
Sxys = Sd*Sxys;
Sxxn = Nd*Sxxn;
Syyn = Nd*Syyn;
Sxyn = Nd*Sxyn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SINNB,COSSB,SINN2B,COSS2B,SINNB2,COSSB2] = ...
dircos1_func(SINB,COSB,SIN2B,COS2B,SINB2,COSB2,NUM)
% function to calculate direction cosine matrices
% for stress transformations
%axpx = cos(x'x)			= cos(thetap - theta)
%axpy = cos(x'y) = sin(x'x)	= sin(thetap - theta)
%aypx = cos(y'x) = -sin(x'x)
%aypy = cos(y'y) = cos(x'x)
% cos(A-B) = cosA*cosB + sinA*sinB
% sin(A-B) = sinA*cosB - cosA*sinB

%axpx = COSB*(COSB') + SINB*(SINB');
%axpy = COSB*(COSB') + SINB*(SINB');
%aypx = -aypx;
%aypy = axpx;

SINNB = SINB*(ones(1,NUM));
COSSB = COSB*(ones(1,NUM));
SINN2B = SIN2B*(ones(1,NUM));
COSS2B = COS2B*(ones(1,NUM));
SINNB2 = SINB2*(ones(1,NUM));
COSSB2 = COSB2*(ones(1,NUM));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs)
% Function to dimension the vectors storing the influence coefficients
% from the first set of calls to coeff
Uxsv   = zeros(nobs,1);
Uysv   = Uxsv;
Uxnv   = Uxsv;
Uynv   = Uxsv;
Sxxsv  = Uxsv;
Syysv  = Uxsv;
Sxysv  = Uxsv;
Sxxnv  = Uxsv;
Syynv  = Uxsv;
Sxynv  = Uxsv;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S1=sig1(Sxx,Syy,Sxy)
% function S1=sig1(SIGxx,SIGyy,SIGxy)
% Calculates the most tensile 2-D principal stress magnitude
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% S1 = most tensile principal stress
% Example 
% S1=sig1(4,-4,3)

S1 = (Sxx+Syy)/2 + sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S2=sig2(Sxx,Syy,Sxy)
% function S1=sig2(Sxx,Syy,Sxy)
% Calculates the most tensile 2-D principal stress magnitude
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% S1 = most tensile principal stress
% Example 
% S1=sig2(4,-4,3)

S2 = (Sxx+Syy)/2 - sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxs,Uxn,Uys,Uyn,Sxxs,Sxxn,Syys,Syyn,Sxys,Sxyn] = ...
square1_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,dimx,dimy)
% function to turn column vectors into square matrices
% and then clear the unneeded column vectors
Uxs = reshape(Uxsv,dimx,dimy);
Uxn = reshape(Uxnv,dimx,dimy);
Uys = reshape(Uysv,dimx,dimy);
Uyn = reshape(Uynv,dimx,dimy);
Sxxs = reshape(Sxxsv,dimx,dimy);
Sxxn = reshape(Sxxnv,dimx,dimy);
Syys = reshape(Syysv,dimx,dimy);
Syyn = reshape(Syynv,dimx,dimy);
Sxys = reshape(Sxysv,dimx,dimy);
Sxyn = reshape(Sxynv,dimx,dimy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ux,Uy,Sxx,Syy,Sxy] = ...
square2_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,NUM2,NUM,dimx,dimy)
% Function square2 to help turn column vectors into square matrices
% and clears unneeded column vectors.  

% First turn the entering column vectors into matrices, 
% where the number of rows is the number of observation points (NUM2)
% and the number of columns equals the number of elements.
% The column vectors are cleared once they are done with.
Uxs = reshape(Uxsv,NUM2,NUM);
Uxn = reshape(Uxnv,NUM2,NUM);
Uys = reshape(Uysv,NUM2,NUM);
Uyn = reshape(Uynv,NUM2,NUM);
Sxxs = reshape(Sxxsv,NUM2,NUM);
Sxxn = reshape(Sxxnv,NUM2,NUM);	
Syys = reshape(Syysv,NUM2,NUM);	
Syyn = reshape(Syynv,NUM2,NUM);	
Sxys = reshape(Sxysv,NUM2,NUM);
Sxyn = reshape(Sxynv,NUM2,NUM);

% Next, sum the entries in each row to superpose
% the contributions of all the elements.
% This step will yield column vectors  
% See p. 488 of the Matlab manual for this
Uxs = sum(Uxs')';
Uxn = sum(Uxn')';
Uys = sum(Uys')';
Uyn = sum(Uyn')';
Sxxs = sum(Sxxs')';
Sxxn = sum(Sxxn')';
Syys = sum(Syys')';
Syyn = sum(Syyn')';
Sxys = sum(Sxys')';
Sxyn = sum(Sxyn')';

% Now, reshape the column vectors into matrices
% that have the dimensions of the observation grid,
% adding the contributions from the shear and
% normal displacement discontinuities 
Ux = reshape(Uxs,dimx,dimy) + reshape(Uxn,dimx,dimy);
Uy = reshape(Uys,dimx,dimy) + reshape(Uyn,dimx,dimy);
Sxx = reshape(Sxxs,dimx,dimy) + reshape(Sxxn,dimx,dimy);
Syy = reshape(Syys,dimx,dimy) + reshape(Syyn,dimx,dimy);
Sxy = reshape(Sxys,dimx,dimy) + reshape(Sxyn,dimx,dimy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tau=taumax(Sxx,Syy,Sxy)
% function tau=taumax(Sxx,Syy,Sxy)
% Calculates maximum shear stress magnitude (2-D)
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% tau = maximum shear stress
% Example 
% tau=taumax(4,-4,3)

tau = sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%