% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.4';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 248;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2007/10/10 16:07:51 $';
fem.version = vrsn;

%% Constants : Note different notation from central load problem
fem.const.L=50; %half length of crack
fem.const.K1=sqrt(fem.const.L*pi); %K1


% Geometry
g1=rect2(2*fem.const.L,fem.const.L,'base','corner','pos',{-fem.const.L,'0'},'rot','0');
parr={point2(0,0)};
g2=geomcoerce('point',parr);
g3=circ2('1','base','center','pos',{'0','0'},'rot','0');
g4=geomcomp({g1,g3},'ns',{'g1','g3'},'sf','g1+g3','edge','none');
%g5=rect2('100','50','base','corner','pos',{-L,'0'},'rot','0');
g6=geomcomp({g4,g1},'ns',{'g4','g5'},'sf','g4*g5','edge','none');

% Analyzed geometry
clear p s
p.objs={g2};
p.name={'PT1'};
p.tags={'g2'};

s.objs={g6};
s.name={'CO2'};
s.tags={'g6'};

fem.draw=struct('p',p,'s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
    'hauto',4, ...
    'hmaxvtx',[4,0.01]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'FlPlaneStrain';
appl.gporder = 4;
appl.cporder = 2;
appl.assignsuffix = '_pn';
clear bnd
bnd.Hy = {0,0,0,0,1};
bnd.Hx = {0,1,0,0,0};
bnd.Fx = {0,0,'K1*fxx/sqrt(2*pi*r)','K1*fxy/sqrt(2*pi*r)',0};
bnd.Fy = {0,0,'K1*fxy/sqrt(2*pi*r)','K1*fyy/sqrt(2*pi*r)',0};
bnd.ind = [2,1,4,1,5,5,3,1,1];
appl.bnd = bnd;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% Define various expression
equ.expr.r='sqrt(x^2+y^2)'; %radial coordinate
equ.expr.th='acos(x/r)'; %radial coordinate
% Normal stresses;
% with normal ni, Stress is sij*ni; normal stress is sij*ni*nj
% radial direction: n=[x,y]/r
% tangential direction: n=[-y,x]/r
equ.expr.srr='(sx_pn*x*x+2*sxy_pn*x*y+sy_pn*y*y)/(r^2)'; %radial coordinate
equ.expr.stt='(sx_pn*y*y+sy_pn*x*x-2*sxy_pn*x*y)/(r^2)'; %tangential coordinate
equ.expr.srt='((sy_pn-sx_pn)*x*y+sxy_pn*(x*x-y*y))/(r^2)'; %tangential coordinate
% Stress solution angular dependence
equ.expr.fxx='cos(th/2)*(1-sin(th/2)*sin(th/2))'; %radial coordinate
equ.expr.fyy='cos(th/2)*(1+sin(th/2)*sin(th/2))'; %radial coordinate
equ.expr.fxy='sin(th/2)*(cos(th/2))^2'; %coordinate

fem.equ = equ;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femstatic(fem, ...
    'solcomp',{'u','v'}, ...
    'outcomp',{'u','v'});

% Save current fem structure for restart purposes
fem0=fem;

%% Plot solution
figure(1); clf; subplot(211)
postplot(fem, ...
    'tridata',{'stt','cont','internal','unit','Pa'}, ...
    'trimap','jet(1024)',...
    'tridlim',[-1,1]*10);
hold on;
postplot(fem, ...
    'contdata',{'stt*sqrt(2*pi*r)/K1','cont','internal','unit','Pa'}, ...
    'contmap',[0,0,0],...
    'contlevels',[0:0.1:1.5],...
    'contbar','off',...
    'contdlim',[-1,1]*10);
axis equal;
box on;
axis ([-1,1,0,0.75]*10)
set(gca,'fontsize',12)
title('color: srr; contour: srr*sqrt(2*pi*r)','fontsize',18);

% Plot in cross-section or along domain
subplot 212
postcrossplot(fem,1,[8,9], ...
    'lindata','srr*sqrt(2*pi)/K1', ...
    'cont','internal', ...
    'linxdata','acos(x)');
hold on;
postcrossplot(fem,1,[8,9], ...
    'lindata','stt*sqrt(2*pi)/K1', ...
    'cont','internal', ...
    'linxdata','acos(x)',...
    'lincolor','r');
postcrossplot(fem,1,[8,9], ...       'lindata','((sy_pn-sx_pn)*x*y+sxy_pn*(x*x-y*y))*sqrt(2*pi)/K1', ...
    'lindata','srt*sqrt(2*pi)/K1', ...
    'cont','internal', ...
    'linxdata','acos(x)',...
    'lincolor','g');
legend('srr','stt','srt','Location','southeast')
fplot(@(x)cos(x/2)*(1+(sin(x/2))^2),[0,pi],'b--')
fplot(@(x)(cos(x/2))^2,[0,pi],'r--')
fplot(@(x)sin(x/2)*cos(x/2),[0,pi],'g--');
xlabel('Angle (radians)','fontsize',18);
ylabel('stress/sqrt(50)','fontsize',18);
box on;
axis([0,pi,-1,1.5])

set(gca,'fontsize',12)