function out = model
clc
%
% MMB.m
%
% Model exported on Sep 2 2018, 21:52 by COMSOL 5.3.0.223.

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

model = ModelUtil.create('Model');



model.label('MMB.mph');

model.comments(['Untitled\n\n']);

%% Parameter
model.param.set('L', '2.0[m]');
model.param.set('H', '1.0[m]');

%% Create component
model.component.create('comp1', true);

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


model.component('comp1').geom('geom1').create('r1', 'Rectangle');
model.component('comp1').geom('geom1').feature('r1').set('pos', [0 0]);
model.component('comp1').geom('geom1').feature('r1').set('size', {'L' 'H'});
model.component('comp1').geom('geom1').run;

model.component('comp1').physics.create('solid', 'SolidMechanics', 'geom1');
model.component('comp1').physics('solid').create('fix1', 'Fixed', 1);
model.component('comp1').physics('solid').feature('fix1').selection.set([1]);
model.component('comp1').physics('solid').create('bndl1', 'BoundaryLoad', 1);
model.component('comp1').physics('solid').feature('bndl1').selection.set([4]);

%% Mesh
model.component('comp1').mesh.create('mesh1');
model.component('comp1').mesh('mesh1').create('fq1', 'FreeQuad');
model.component('comp1').mesh('mesh1').feature('fq1').selection.geom('geom1', 2);
model.component('comp1').mesh('mesh1').feature('fq1').selection.set([1]);

model.component('comp1').mesh('mesh1').feature('size').set('custom', 'on');
model.component('comp1').mesh('mesh1').feature('size').set('hmax', 0.1);

model.component('comp1').mesh('mesh1').run;


model.component('comp1').view('view1').axis.set('xmin', -0.02999996766448021);
model.component('comp1').view('view1').axis.set('xmax', 1.2300000190734863);
model.component('comp1').view('view1').axis.set('ymin', -0.7800000309944153);
model.component('comp1').view('view1').axis.set('ymax', 1.380000114440918);
model.component('comp1').view('view1').axis.set('abstractviewlratio', -0.04999994859099388);
model.component('comp1').view('view1').axis.set('abstractviewrratio', 0.04999994859099388);
model.component('comp1').view('view1').axis.set('abstractviewbratio', -0.6314285397529602);
model.component('comp1').view('view1').axis.set('abstractviewtratio', 0.6314285397529602);
model.component('comp1').view('view1').axis.set('abstractviewxscale', 0.0025142855010926723);
model.component('comp1').view('view1').axis.set('abstractviewyscale', 0.002514285733923316);
%% Materials
model.component('comp1').material.create('mat1', 'Common');
model.component('comp1').material('mat1').propertyGroup('def').set('youngsmodulus', '2e9');
model.component('comp1').material('mat1').propertyGroup('def').set('poissonsratio', '0.35');
model.component('comp1').material('mat1').propertyGroup('def').set('density', '7800');

%% BCs for load
model.component('comp1').physics('solid').feature('bndl1').set('FperArea', {'0'; '-10000*(Y>((0.5-0.1)*H))*(Y<((0.5+0.1)*H))'; '0'});

ModelUtil.tags


%% Study
model.study.create('std1');
model.study('std1').create('stat', 'Stationary');

%% run study
model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').attach('std1');
model.sol('sol1').create('st1', 'StudyStep');
model.sol('sol1').create('v1', 'Variables');
model.sol('sol1').create('s1', 'Stationary');
model.sol('sol1').feature('s1').create('fc1', 'FullyCoupled');
model.sol('sol1').feature('s1').feature.remove('fcDef');
model.sol('sol1').attach('std1');
model.sol('sol1').runAll;
%% Ploting
model.result.create('pg1', 'PlotGroup2D');
model.result('pg1').create('surf1', 'Surface');
model.result('pg1').feature('surf1').create('def', 'Deform');

%% Obtain matrix and vector
fu=mphxmeshinfo(model);
MA = mphmatrix(model ,'sol1', ...
    'Out', {'Kc','Lc','Null','ud','uscale'},...
    'initmethod','sol','initsol','zero');
U = mphgetu(model);

Uc = MA.Null*(MA.Kc\MA.Lc);
U0 = Uc+MA.ud;
U1 = (1+U0).*MA.uscale;
%%
HOME=pwd;
xx=strcat(HOME,'/MMB.mph')
model.save(xx)
ModelUtil.remove('MMB.mph')

out = model;