% defnosclstat.m
%
% Model exported on Dec 12 2013, 21:04 by COMSOL 4.3.2.152.

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');

model.modelPath('C:\harsh\AOcomsol\coarse');

model.modelNode.create('mod1');

model.geom.create('geom1', 2);

model.mesh.create('mesh1', 'geom1');

model.physics.create('plate', 'Plate', 'geom1');

model.study.create('std1');
model.study('std1').feature.create('stat', 'Stationary');
model.study('std1').feature('stat').activate('plate', true);

model.geom('geom1').feature.create('c1', 'Circle');
model.geom('geom1').feature('c1').set('r', '2.636/2');
model.geom('geom1').run('c1');
model.geom('geom1').feature.create('c2', 'Circle');
model.geom('geom1').feature('c2').set('r', '0.697/2');
model.geom('geom1').run('c2');
model.geom('geom1').feature.create('co1', 'Compose');
model.geom('geom1').feature('co1').set('formula', 'c1-c2');
model.geom('geom1').run;

model.material.create('mat1');
model.material('mat1').propertyGroup('def').set('youngsmodulus', {'9.1e10'});
model.material('mat1').propertyGroup('def').set('poissonsratio', {'0.243'});
model.material('mat1').propertyGroup('def').set('density', {'2530'});

model.physics('plate').prop('d').set('d', 1, '1e-3[m]');
model.physics('plate').feature.create('fix1', 'Fixed', 1);
model.physics('plate').feature('fix1').selection.set([3 4 6 7]);
model.physics('plate').feature.create('el1', 'EdgeLoad', 1);
model.physics('plate').feature('el1').selection.set([1 2 5 8]);
model.physics('plate').feature('el1').set('FeperLength', {'0' '0' '1'});

model.mesh('mesh1').autoMeshSize(7);

model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'stat');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature('v1').set('control', 'stat');
model.sol('sol1').feature.create('s1', 'Stationary');
model.sol('sol1').feature('s1').feature.create('fc1', 'FullyCoupled');
model.sol('sol1').feature('s1').feature('fc1').set('termonres', 'auto');
model.sol('sol1').feature('s1').feature('fc1').set('reserrfact', 1000);
model.sol('sol1').feature('s1').feature('fc1').set('termonres', 'auto');
model.sol('sol1').feature('s1').feature('fc1').set('reserrfact', 1000);
model.sol('sol1').feature('s1').feature.remove('fcDef');
model.sol('sol1').attach('std1');
model.sol('sol1').feature('v1').set('scalemethod', 'none');

model.result.create('pg1', 2);
model.result('pg1').set('data', 'dset1');
model.result('pg1').name('Stress Top (plate)');
model.result('pg1').feature.create('surf1', 'Surface');
model.result('pg1').feature('surf1').set('expr', {'plate.mises'});
model.result('pg1').feature('surf1').feature.create('hght', 'Height');
model.result('pg1').feature('surf1').feature('hght').set('expr', {'w'});
model.result('pg1').feature('surf1').set('const', {'plate.z' '1' 'Height of evaluation in shell, [-1,1]'; 'plate.refpntx' '0' 'Reference point for moment computation, x coordinate'; 'plate.refpnty' '0' 'Reference point for moment computation, y coordinate'; 'plate.refpntz' '0' 'Reference point for moment computation, z coordinate'});

model.sol('sol1').runAll;

model.result('pg1').run;
model.result('pg1').feature.create('surf2', 'Surface');
model.result('pg1').run;

%% Part I edit 


dat1=mpheval(model,{'w'});
mphplot(dat1,'mesh','off','rangenum',1);



U=mphgetu(model); % This is true solution 

% Solution by Matrix Extraction. Followed steps given in  "Livelink for MATLAB" user guide version 4.3b  example given on page 141-142 Chapter 4 
infoC=mphsolinfo(model,'solname','sol1');
MC=mphmatrix(model,'sol1','out',{'Ec','Kc','Dc','Lc','Null','ud','uscale'},'initmethod','sol','initsol','sol1');
UC=MC.Null*(MC.Kc\MC.Lc);
U0=UC+MC.ud;
U1=(1+U0).*MC.uscale;
infoB=mphxmeshinfo(model,'solname','sol1');