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

% Geometry
g1=circ2('1','base','center','pos',{'0','0'},'rot','0');
g2=circ2('100','base','center','pos',{'0','0'},'rot','0');
g3=geomcomp({g2,g1},'ns',{'g2','g1'},'sf','g2-g1','edge','none');

% Analyzed geometry
clear s
s.objs={g3};
s.name={'CO1'};
s.tags={'g3'};

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

% Initialize mesh
fem.mesh=meshinit(fem, ...
    'hauto',5, ...
    'hmaxedg',[3,0.05,4,0.05,6,0.05,7,0.05]);

% (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.Fx = {0,'x'};
bnd.Fy = {0,'y'};
bnd.ind = [1,1,2,2,1,2,2,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
% 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
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;
postplot(fem, ...
    'tridata',{'srr','cont','internal','unit','Pa'}, ...    'tridata',{'srr','cont','internal','unit','Pa'}, ...
    'trimap','jet(1024)',...
    'tridlim',[-1,0]);
axis equal;
axis ([-1,1,-1,1]*5);
set(gca,'fontsize',18);
title('Radial stress','fontsize',24)
% set(1,'renderer','painters');
% , ...
%     'title','Surface: von Mises stress [Pa]', ...
%     'axis',[-1,1,-1.037037037037037,1.037037037037037]);
%%
figure(2); clf;
orient tall
% Plot in cross-section or along domain
subplot 211; hold on;
postcrossplot(fem,1,[0 10;0 0], ...
    'lindata','-srr*(r^2)', ...
    'linxdata','r',...
    'lincolor','r',...
    'linewidth',2);
postcrossplot(fem,1,[0 0;0 10], ...
    'lindata','-srr*(r^2)', ...
    'linxdata','r',...
    'lincolor','g',...
    'linewidth',2);
postcrossplot(fem,1,[0 10;0 10], ...
    'lindata','-srr*(r^2)', ...
    'linxdata','r',...
    'lincolor','b',...
    'linewidth',2);
legend('x-axis','y-axis','45 degrees','location','northwest')
axis([0,10,0.9,1.1]);
set(gca,'fontsize',18);
xlabel('Radial distance r','fontsize',24)
ylabel('Radial stress*r^2','fontsize',24);box on;
%
subplot 212; hold on;
postcrossplot(fem,1,[0 10;0 0], ...
    'lindata','stt*(r^2)', ...
    'linxdata','r',...
    'lincolor','r',...
    'linewidth',2);
postcrossplot(fem,1,[0 0;0 10], ...
    'lindata','stt*(r^2)', ...
    'linxdata','r',...
    'lincolor','g',...
    'linewidth',2);
postcrossplot(fem,1,[0 10;0 10], ...
    'lindata','stt*(r^2)', ...
    'linxdata','r',...
    'lincolor','b',...
    'linewidth',2);
axis([0,10,0.9,1.1]);
set(gca,'fontsize',18);
xlabel('Radial distance r','fontsize',24)
ylabel('Hoop stress*r^2','fontsize',24);

box on;
set(2,'renderer','painters')

