function out = model
%
% levelset_wang.m
%
% Model exported on Jan 16 2019, 23:33 by COMSOL 5.4.0.246.


DomainWidth=2;
DomainHight=1;
ENPC=40;
ENPR=80;
nelx = ENPR;
nely = ENPC;  
FEAInterval=1;
PlotInterval=1;

% Step 1: Data initialization
EW = DomainWidth / ENPR;       %  The width of each finite element.
EH = DomainHight / ENPC;        %   The hight of each finite element.
M = [ ENPC + 1 , ENPR + 1 ];
[ x ,y ] = meshgrid( EW * [ -0.5 : ENPR + 0.5 ] , EH * [ -0.5 : ENPC + 0.5 ]);
[ FENd.x, FENd.y, FirstNdPCol ] = MakeNodes(ENPR,ENPC,EW,EH);
Ele.NdID = MakeElements( ENPR, ENPC, FirstNdPCol );
LSgrid.x = x(:); LSgrid.y = y(:);                           % The coordinates of each Level Set grid

 for i = 1 : length(Ele.NdID(:,1))
    Ele.LSgridID(i) = find((LSgrid.x - FENd.x(Ele.NdID(i,1)) - EW/2).^2 +...
       (LSgrid.y - FENd.y(Ele.NdID(i,1)) - EH/2).^2 <= 100*eps);
 end;

 
cx = DomainWidth / 200 * [ 33.33  100  166.67  0   66.67  133.33  200  33.33  100  166.67  0   66.67  133.33  200  33.33  100  166.67];
cy = DomainHight / 100 * [   0     0     0     25   25      25     25    50    50    50    75    75     75     75    100   100   100]; 
% % 
for i = 1 : length( cx )
       tmpPhi( :, i ) = - sqrt ( ( LSgrid . x - cx ( i ) ) .^2 + ( LSgrid . y  - cy ( i ) ) .^2 ) + DomainHight/8;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      
end;
 
LSgrid.Phi =  -(max(tmpPhi.')).'; 

LSgrid.Phi((LSgrid.x - min(LSgrid.x))  .* (LSgrid.x - max(LSgrid.x))  .* (LSgrid.y - max(LSgrid.y)) .* (LSgrid.y - min(LSgrid.y)) <= 100*eps) = -1e-6;
FENd.Phi = griddata( LSgrid.x, LSgrid.y, LSgrid.Phi, FENd.x, FENd.y, 'cubic'); 


% Generate phi
phi3 = [FENd.x, FENd.y, FENd.Phi]';
% phi3 = FENd.Phi';


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

model = ModelUtil.create('Model');

model.modelPath('C:\Users\sayinde\Desktop\wang_comsol_matlab');

model.label('levelset_wang.mph');

model.comments(['Levelset opt\n\n']);

model.param.set('volfrac', '0.5');
model.param.set('ht', '1.0');
model.param.set('alpha', '0.3');
model.param.set('d', '0.01');
model.param.set('dx', '2/3');
model.param.set('dy', '1/2');
model.param.set('x0', 'dx/2');
model.param.set('y0', 'dy/2');
model.param.set('r', '1/8');
model.param.set('gamma', '0.1');
model.param.set('sigma', '5');
model.param.set('beta', '1');
model.param.set('kappa', '0.00125');
model.param.set('x00', 'x0');
model.param.set('x01', '1');
model.param.set('x02', '2-x0');
model.param.set('y00', '0');
model.param.set('y01', '0');
model.param.set('y02', '0');
model.param.set('x10', '0');
model.param.set('x11', '2*x0');
model.param.set('x12', '4*x0');
model.param.set('x13', '2');
model.param.set('y10', '0.25');
model.param.set('y11', '0.25');
model.param.set('y12', '0.25');
model.param.set('y13', '0.25');
model.param.set('x20', 'x0');
model.param.set('x21', '1');
model.param.set('x22', '2-x0');
model.param.set('y20', '0.5');
model.param.set('y21', '0.5');
model.param.set('y22', '0.5');
model.param.set('x30', '0');
model.param.set('x31', '2*x0');
model.param.set('x32', '4*x0');
model.param.set('x33', '2');
model.param.set('y30', '0.75');
model.param.set('y31', '0.75');
model.param.set('y32', '0.75');
model.param.set('y33', '0.75');
model.param.set('x40', 'x0');
model.param.set('x41', '1');
model.param.set('x42', '2-x0');
model.param.set('y40', '1');
model.param.set('y41', '1');
model.param.set('y42', '1');
model.param.set('rho_void', '1e-3');
model.param.set('rho_mat', '1');
model.param.set('delta', '0.05');

model.component.create('mod1', false);

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

model.component('mod1').defineLocalCoord(false);

model.result.table.create('evl2', 'Table');
model.result.table.create('tbl1', 'Table');

model.component('mod1').mesh.create('mesh1');

model.component('mod1').geom('geom1').repairTolType('relative');
model.component('mod1').geom('geom1').create('R1', 'Rectangle');
model.component('mod1').geom('geom1').feature('R1').set('pos', '0.0,0.0');
model.component('mod1').geom('geom1').feature('R1').set('size', [2 1]);
model.component('mod1').geom('geom1').create('pt1', 'Point');
model.component('mod1').geom('geom1').feature('pt1').set('p', [2 0.5]);
model.component('mod1').geom('geom1').create('pt2', 'Point');
model.component('mod1').geom('geom1').feature('pt2').set('p', [0 0.5]);
model.component('mod1').geom('geom1').feature('fin').set('repairtoltype', 'relative');
model.component('mod1').geom('geom1').run;

model.variable.create('var1');
model.variable('var1').set('rho', '(rho_mat*(mod1.phi>=delta)+rho_void*(mod1.phi<-delta)+(0.75 * (rho_mat - rho_void ) * ( (mod1.phi / delta) - (mod1.phi / delta)^3 / 3.0 ) + 0.5 * ( rho_mat + rho_void ))*(-delta<=mod1.phi&&mod1.phi<delta))');
model.variable('var1').set('df', '-2*d(mod1.ps.Ws,rho)*d(rho,mod1.phi)');
model.variable('var1').set('dist01', 'sqrt((mod(x+3*x0,dx)-x0)^2+(mod(y+2*y0,dy)-y0)^2)-r');
model.variable('var1').set('dist02', 'sqrt((mod(x+2*x0,dx)-x0)^2+(mod(y+3*y0,dy)-y0)^2)-r');
model.variable('var1').set('dis00', '-sqrt((x-x00)^2+(y-y00)^2)+r');
model.variable('var1').set('dis01', '-sqrt((x-x01)^2+(y-y01)^2)+r');
model.variable('var1').set('dis02', '-sqrt((x-x02)^2+(y-y02)^2)+r');
model.variable('var1').set('dis10', '-sqrt((x-x10)^2+(y-y10)^2)+r');
model.variable('var1').set('dis11', '-sqrt((x-x11)^2+(y-y11)^2)+r');
model.variable('var1').set('dis12', '-sqrt((x-x12)^2+(y-y12)^2)+r');
model.variable('var1').set('dis13', '-sqrt((x-x13)^2+(y-y13)^2)+r');
model.variable('var1').set('dis20', '-sqrt((x-x20)^2+(y-y20)^2)+r');
model.variable('var1').set('dis21', '-sqrt((x-x21)^2+(y-y21)^2)+r');
model.variable('var1').set('dis22', '-sqrt((x-x22)^2+(y-y22)^2)+r');
model.variable('var1').set('dis30', '-sqrt((x-x30)^2+(y-y30)^2)+r');
model.variable('var1').set('dis31', '-sqrt((x-x31)^2+(y-y31)^2)+r');
model.variable('var1').set('dis32', '-sqrt((x-x32)^2+(y-y32)^2)+r');
model.variable('var1').set('dis33', '-sqrt((x-x33)^2+(y-y33)^2)+r');
model.variable('var1').set('dis40', '-sqrt((x-x40)^2+(y-y40)^2)+r');
model.variable('var1').set('dis41', '-sqrt((x-x41)^2+(y-y41)^2)+r');
model.variable('var1').set('dis42', '-sqrt((x-x42)^2+(y-y42)^2)+r');
model.variable('var1').set('phi0', '-max(max(max(max(max(max(max(max(max(max(max(max(max(max(max(max(dis00,dis01),dis02),dis10),dis11),dis12),dis13),dis20),dis21),dis22),dis30),dis31),dis32),dis33),dis40),dis41),dis42)');
% model.variable('var1').set('phi0', {num2str(phi3)});
model.variable('var1').set('df_drho', 'd(mod1.ps.Ws,rho)');
model.variable('var1').set('drho_dphi', 'd(rho,mod1.phi)');
model.component('mod1').variable.create('cplvar1');
model.component('mod1').variable('cplvar1').set('intcpl_source_volume', 'rho');
model.component('mod1').variable('cplvar1').set('intcpl_source_voldom', '1');
model.component('mod1').variable('cplvar1').set('intcpl_source_f', '2*ps.Ws');
model.component('mod1').variable('cplvar1').selection.geom('geom1', 2);
model.component('mod1').variable('cplvar1').selection.set([1]);
model.variable.create('cplvar2');
model.variable('cplvar2').set('volume', 'mod1.intop1(intcpl_source_volume)');
model.variable('cplvar2').set('voldom', 'mod1.intop1(intcpl_source_voldom)');
model.variable('cplvar2').set('f', 'mod1.intop1(intcpl_source_f)');

model.component('mod1').cpl.create('intop1', 'Integration');
model.component('mod1').cpl('intop1').selection.set([1]);

model.component('mod1').physics.create('ps', 'SolidMechanics', 'geom1');
model.component('mod1').physics('ps').identifier('ps');
model.component('mod1').physics('ps').field('displacement').component({'u' 'v' 'w_ps'});
model.component('mod1').physics('ps').create('pl1', 'PointLoad', 0);
model.component('mod1').physics('ps').feature('pl1').selection.set([5]);
model.component('mod1').physics('ps').create('disp1', 'Displacement1', 1);
model.component('mod1').physics('ps').feature('disp1').selection.set([1 3]);
model.component('mod1').physics.create('w', 'WeakFormPDE', 'geom1');
model.component('mod1').physics('w').field('dimensionless').field('phi');
model.component('mod1').physics('w').field('dimensionless').component({'phi'});

model.component('mod1').mesh('mesh1').create('map1', 'Map');
model.component('mod1').mesh('mesh1').feature('map1').create('dis1', 'Distribution');
model.component('mod1').mesh('mesh1').feature('map1').create('dis3', 'Distribution');
model.component('mod1').mesh('mesh1').feature('map1').feature('dis1').selection.set([1 3 5 6]);
model.component('mod1').mesh('mesh1').feature('map1').feature('dis3').selection.set([2 4]);

model.result.table('evl2').label('Evaluation 2D');
model.result.table('evl2').comments('Interactive 2D values');
model.result.table('tbl1').comments('Global Evaluation 1 (solid.Ws_ps)');

model.capeopen.label('Thermodynamics Package');

model.component('mod1').view('view1').axis.set('xmin', -0.364271342754364);
model.component('mod1').view('view1').axis.set('xmax', 2.412327527999878);
model.component('mod1').view('view1').axis.set('ymin', -0.7242448329925537);
model.component('mod1').view('view1').axis.set('ymax', 1.582151174545288);

model.component('mod1').cpl('intop1').set('intorder', 2);

model.component('mod1').physics('ps').prop('ShapeProperty').set('order_displacement', 1);
model.component('mod1').physics('ps').prop('EquationForm').set('form', 'Stationary');
model.component('mod1').physics('ps').prop('d').set('d', 0.01);
model.component('mod1').physics('ps').prop('cref').set('cref', 'ps.cp');
model.component('mod1').physics('ps').feature('lemm1').set('E', '1*rho');
model.component('mod1').physics('ps').feature('lemm1').set('nu', 0.3);
model.component('mod1').physics('ps').feature('lemm1').set('rho', 1);
model.component('mod1').physics('ps').feature('lemm1').set('minput_frequency', 'freq');
model.component('mod1').physics('ps').feature('pl1').set('Fp', [0; -1; 0]);
model.component('mod1').physics('ps').feature('disp1').set('Direction', [1; 1; 0]);
model.component('mod1').physics('w').prop('ShapeProperty').set('order', 1);
model.component('mod1').physics('w').feature('wfeq1').set('weak', 0);
model.component('mod1').physics('w').feature('init1').set('phi', 'phi0');

model.component('mod1').mesh('mesh1').feature('size').set('custom', 'on');
model.component('mod1').mesh('mesh1').feature('size').set('hmax', 0.05);
model.component('mod1').mesh('mesh1').feature('map1').feature('dis1').set('numelem', 20);
model.component('mod1').mesh('mesh1').feature('map1').feature('dis3').set('numelem', 80);
model.component('mod1').mesh('mesh1').run;

model.component('mod1').physics('ps').feature('lemm1').set('E_mat', 'userdef');
model.component('mod1').physics('ps').feature('lemm1').set('nu_mat', 'userdef');
model.component('mod1').physics('ps').feature('lemm1').set('rho_mat', 'userdef');
model.component('mod1').physics('ps').feature('lemm1').set('minput_strainreferencetemperature_src', 'userdef');

model.study.create('std1');
model.study('std1').create('stat', 'Stationary');
model.study('std1').feature('stat').set('activate', {'ps' 'on' 'w' 'off' 'frame:spatial1' 'on'});

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.result.create('pg1', 'PlotGroup2D');
model.result('pg1').create('surf1', 'Surface');
model.result('pg1').feature('surf1').create('def', 'Deform');
model.result.export.create('data1', 'Data');

model.sol('sol1').attach('std1');
model.sol('sol1').feature('s1').feature('aDef').set('cachepattern', true);
model.sol('sol1').runAll;

model.result('pg1').label('Stress (ps)');
model.result('pg1').feature('surf1').set('expr', 'phi0');
model.result('pg1').feature('surf1').set('descr', '');
model.result('pg1').feature('surf1').set('const', {'ps.refpntx' '0' 'Reference point for moment computation, x coordinate'; 'ps.refpnty' '0' 'Reference point for moment computation, y coordinate'; 'ps.refpntz' '0' 'Reference point for moment computation, z coordinate'});
model.result('pg1').feature('surf1').set('colortable', 'RainbowLight');
model.result('pg1').feature('surf1').set('resolution', 'normal');
model.result('pg1').feature('surf1').feature('def').set('scale', 1.750814484738484E-5);
model.result('pg1').feature('surf1').feature('def').set('scaleactive', false);
model.result.export('data1').set('expr', {'ps.Ws' 'phi'});
model.result.export('data1').set('unit', {'J/m^3' '1'});
model.result.export('data1').set('descr', {'Elastic strain energy density' 'Dependent variable phi'});
model.result.export('data1').set('const', {'ps.refpntx' '0' 'Reference point for moment computation, x coordinate'; 'ps.refpnty' '0' 'Reference point for moment computation, y coordinate'; 'ps.refpntz' '0' 'Reference point for moment computation, z coordinate'});
model.result.export('data1').set('filename', 'C:\Users\sayinde\Desktop\wang_comsol_matlab\eval.txt');
model.result.export('data1').set('header', false);
model.result.export('data1').set('fullprec', false);
model.result.export('data1').set('sort', false);
model.result.export('data1').run;


model.result.create('pg10', 'PlotGroup2D');
model.result('pg10').run;
model.result('pg10').create('con4', 'Contour');
model.result('pg10').feature('con4').set('expr', 'phi');
model.result('pg10').feature('con4').set('descr', 'phi');
% model.result('pg2').feature('con1').set('descr', 'Principal stress direction 1, Gauss-point evaluation, x component');
model.result('pg10').feature('con4').label('phi');
model.result('pg10').feature('con4').set('contourtype', 'filled');
figure(45)
mphplot(model,'pg10','rangenum',1)
model.result('pg10').run;

% model.result.export.create('data22', 'Data');
% model.result.export('data22').set('expr', {'phi0'});
% model.result.export('data22').set('descr', {'descr', 'phi0'});
% model.result.export('data22').set('filename', 'phi0.txt');
% model.result.export('data22').run;
% get(U);
model.sol('sol1').getU
U = mphgetu(model)
assignin('base', 'U', U);  

% out = model;