% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.2 (COMSOL 3.2.0.222, $Date: 2005/09/01 18:02:30 $)

flclear fem

%% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.2';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 222;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2005/09/01 18:02:30 $';
fem.version = vrsn;


%% Constants : Note different notation from central load problem
fem.const.E=70e9; %Young's modulus
fem.const.nu=0.3; %Poisson ratio
fem.const.H=34e3; %Thickness of elastic plate
fem.const.D=fem.const.E*fem.const.H^3/(12*(1-fem.const.nu^2)); %Rigidity
fem.const.rhop=3300; %Density of mantle
fem.const.rhov=3000; %Density of volcanic load
fem.const.rhow=1000; %Density of water
fem.const.g=9.81; %Gravity
fem.const.Hv=35e3; %Height of volcanic pile
fem.const.Wv=150e3; %Width of volcanic pile
fem.const.yv=0; %location of tip of the chain
fem.const.xv=0; %location of tip of the chain
fem.const.B=5e3; %Regional bathymetry (for Winkler force)

%% Geometry: Just a simple 1000x2000 km rectangle 
% Construct plate as rectangular 1000x2000 km solid object (type 's')
g1=rect2('1000e3','2000e3','base','corner','pos',{'0','-1000e3'},'rot','0');
s.objs={g1};
s.name={'plate'};
s.tags={'g1'};
% add point at tip of chain, mainly for mesh refinement.
p1=point2(0,fem.const.yv);
p.objs={p1};
p.name={'tip'};
p.tags={'p1'};
% assemble solid and point
fem.draw=struct('s',s,'p',p);
fem.geom=geomcsg(fem);
% visualize geometry with labels on each edge
figure(1); clf; geomplot(fem,'edgelabel','on')

%% Initialize mesh 
% triangular mesh refined alog edge number 1
fem.mesh=meshinit(fem, ...
    'hmax',20e3,...
    'hmaxvtx',[3,1e3],...
    'hmaxedg',[3,5e3]);
figure(1); clf; meshplot(fem)

%% (Default values are not included)

% Application mode 1: Mindlin plate
clear appl
appl.mode.class = 'SmeMindlin';
appl.module = 'SME';
appl.gporder = 4;
appl.cporder = 1;
appl.assignsuffix = '_smdrm';
clear bnd
% Boundary conditions: 
% For symmetry, impose no rotation
% far edges are fixed
% Top boundary is free, but could be no rotation either
bnd.constrlocaltype = {'free','fixed','rotation'};
bnd.ind = [3,2,3,1,2];
appl.bnd = bnd;
clear equ
% Enter physical properties
equ.rho = 'rhop'; %does not actually enters problem
equ.nu = 'nu';
equ.Fz = 'F';
equ.thickness = 'H';
equ.E = 'E';
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;
fem.border = 1;
fem.units = 'SI';

%% Subdomain expressions
clear equ
equ.ind = [1];
% Name three degrees of freedom of application
equ.dim = {'w','thx','thy'};
% Define various expression
equ.expr.r='sqrt((x-xv)^2+(y-yv)^2)';
equ.expr.topov='((y>=yv)*max(1-abs(x-xv)/Wv,0)+(y<yv)*max(1-r/Wv,0))*Hv';
equ.expr.Fvolc='-topov*rhov*g';
% Density contrast if always underwater...
% equ.expr.drho='rhop-rhow';
% Modified for water depth of B
equ.expr.drho='rhop-(rhow*((topov+w-B)<0))';
equ.expr.Fwinkler='-drho*g*w';
equ.expr.F='Fwinkler+Fvolc';
% Note, this is the COMSOL original output:
% equ.expr = {'Fwinkler','-rhop*g*w', ...    
%   'r','sqrt((x-xv)^2+(y-yv)^2)', ...
%     'topov','((y>=yv)*max(1-abs(x-xv)/Wv,0)+(y<yv)*max(1-r/Wv,0))*Hv', ...
%   'Fvolc','-topov*rhov*g', ...
%   'F','Fwinkler+Fvolc'};
fem.equ = equ;

%% Assemble problem
% Multiphysics
fem=multiphysics(fem);
% Extend mesh
fem.xmesh=meshextend(fem);

%% Solve problem
% Note that an iterative solver (nlin) is required because the density
% contrast is a non-linear function of the final topography, which is unknown.
% fem.sol=femnlin(fem, ...
%     'solcomp',{'thx','w','thy','thn_smdrm'}, ...
%     'outcomp',{'thx','w','thy','thn_smdrm'});
% if the chain was always underwater or always above, a linear solver
% is sufficient (and much faster)
fem.sol=femlin(fem, ...
    'solcomp',{'thx','w','thy','thn_smdrm'}, ...
    'outcomp',{'thx','w','thy','thn_smdrm'});
           
           % Save current fem structure for restart purposes
fem0=fem;

%% Plot solution: map view
figure(1)
postplot(fem, ...
    'tridata',{'w+topov-B','cont','internal'}, ...
    'trimap','jet(1024)', ...
    'title','Topography (m)', ...
    'refine',1);

%% Plot solution: 3D reconstruction
figure(2)
clf;
postplot(fem, ...
    'tridata',{'min(max(w,-1000),1000)','cont','internal'}, ...
    'tridlim',[-1,1]*1000,...
    'triz','w+topov', ...
    'trimap','jet(1024)', ...
    'title','Plate deflection ploted on bathymetry (m) VE: 20x1', ...
    'refine',5, ...
    'grid','on');
set(gca,'dataaspectratio',[1,1,5e-2]);
axis([0,1000,-1000,1000,0,10]*1000)

view(30,40);
HL1=lightangle(120,05);
HL2=lightangle(50,75);

box on
grid off
