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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%