% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.3 (COMSOL 3.3.0.405, $Date: 2006/08/31 18:03:47 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.3';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 405;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2006/08/31 18:03:47 $';
fem.version = vrsn;

% Geometry
g1=rect2('80','40','base','corner','pos',{'-40','-40'},'rot','0');

fem.const.TM=1370;%degC, mantle temperature

% Geometry
% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

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

% Create mapped quad mesh
fem.mesh=meshmap(fem, ...
    'edgegroups',{{[2],[4],[3],[1]}}, ...
    'hauto',1);

% (Default values are not included)

% Application mode 1: convection and conduction
clear appl
appl.mode.class = 'FlConvCond';
appl.assignsuffix = '_cc';
clear prop %static analysis
prop.analysis='static';
appl.prop = prop;
clear bnd
% Boundary type 1: temperature=0   (surface)
% Boundary type 2: convective flux (edges)
% Boundary type 3: temperature=1   (bottom)
bnd.type = {'T','qc','T'};
bnd.T0 = {0,0,1};
bnd.ind = [2,3,1,2];
appl.bnd = bnd;
clear equ
equ.k = 1;  % thermal conductivity
equ.rho = 1;% density
equ.C = 1;  % heat capacity
equ.v = 'vz';%analytical velocity
equ.u = 'vx';%analytical velocity
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Subdomain settings
clear equ
equ.ind = [1];
equ.dim = {'T'};

% Subdomain expressions
equ.expr = {'r','sqrt(x^2+y^2)', ... %radius from axis
    'th','asin(x/r)', ... %angle from vertical
    'vx','(2/pi)*(th-sin(th)*cos(th))', ...%corner flow solution
    'vz','(2/pi)*cos(th)^2'}; %corner flow solution
fem.equ = equ;


%% Multiphysics
fem=multiphysics(fem);

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

% Solve problem: problem is linear; Default solver recognizes that
fem.sol=femstatic(fem, ...
    'solcomp',{'T'}, ...
    'outcomp',{'T'});

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

%% Plot solution
figure(1); clf
postplot(fem, ...
    'contdata',{'TM*T','cont','internal'}, ... scale temperature
    'contlevels',[0:100:1300], ...
    'contmap',hot(1024),...
    'contdlim',[0,1300],...
    'contbar','on' ,...
    'flowdata',{'vx','vz'}, ...
    'flowstart',{linspace(-20,20,20),linspace(-20,-20,20)}, ...
    'flowcolor',[0,0,0] ,...
    'arrowdata',{'vx','vz'}, ...
    'arrowxspacing',15, ...
    'arrowyspacing',15, ...
    'arrowtype','arrow', ...
    'arrowstyle','proportional', ...
    'arrowcolor',[0.0,0.0,1.0]);
axis equal;
axis([-20,20,-20,0])

