%mesh resolution
resolution=10000;
riftres=500;

md=model;

%refine mesh for rift
expcoarsen('Rifts2.exp','Rifts.exp',riftres);

%md=mesh(md,'DomainOutline.exp',50000);
md=mesh(md,'DomainOutline.exp','Rifts2.exp',resolution);
md=meshprocessrifts(md);
md.numberofelements
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');

md.rifts.fill=MelangeEnum();
md.rifts.fraction=0;
md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 80 ';
md.cluster='astrid';
md.np=8;
md.penalty_lock=2;
md.rifts.fractionincrement=.01;
md=solve(md,'analysis_type','diagnostic');md=tres(md,'diagnostic');
plotmodel(md,'data','riftrelvel','data','riftvel','data','vel');