% topopt problem by Kristian

clc; clear all; close(gcf); close(gcf); close(gcf); close(gcf); close(gcf);
import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');

model.modelNode.create('mod1');

% DEFINE GEOMETRY
model.geom.create('geom1', 2);
model.geom('geom1').feature.create('r1', 'Rectangle');
model.geom('geom1').feature('r1').set('base', 'corner');
model.geom('geom1').feature('r1').set('pos', {'0' '0'});
model.geom('geom1').feature('r1').set('lx', '2');
model.geom('geom1').feature('r1').set('ly', '1');

model.geom('geom1').feature.create('r2', 'Rectangle');
model.geom('geom1').feature('r2').set('base', 'corner');
model.geom('geom1').feature('r2').set('pos', {'2' '0'});
model.geom('geom1').feature('r2').set('lx', '3');
model.geom('geom1').feature('r2').set('ly', '1');

model.geom('geom1').feature.create('r3', 'Rectangle');
model.geom('geom1').feature('r3').set('base', 'corner');
model.geom('geom1').feature('r3').set('pos', {'5' '0'});
model.geom('geom1').feature('r3').set('lx', '2');
model.geom('geom1').feature('r3').set('ly', '1');

model.geom('geom1').feature.create('pt1', 'Point');
model.geom('geom1').feature('pt1').set('p', {'3.5' '0.5'});
model.geom('geom1').run;

%%% DEFINE MESH
model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').feature.create('ftri3', 'FreeTri');
model.mesh('mesh1').feature('size').set('custom', 'on');
model.mesh('mesh1').feature('size').set('hmax', '0.1');
model.mesh('mesh1').feature('size').set('hmin', '0.01');
model.mesh('mesh1').run;

%%% DEFINE PHYSICS
model.physics.create('spf', 'LaminarFlow', 'geom1');
model.physics('spf').prop('CompressibilityProperty').set('Compressibility', 1, 'Incompressible');
model.physics('spf').prop('EquationForm').set('form', 1, 'Stationary');
model.physics('spf').prop('ShapeProperty').set('order_fluid', 1, '2');
model.physics('spf').prop('Stabilization').set('StreamlineDiffusion', 1, '0');
model.physics('spf').prop('Stabilization').set('CrosswindDiffusion', 1, '0');

model.physics('spf').feature('fp1').set('rho_mat', 1, 'userdef');
model.physics('spf').feature('fp1').set('rho', 1, '1');
model.physics('spf').feature('fp1').set('mu_mat', 1, 'userdef');
model.physics('spf').feature('fp1').set('mu', 1, '1');


model.physics('spf').feature.create('bs1', 'BoundaryStress', 1);
model.physics('spf').feature('bs1').set('BoundaryCondition', 1, 'NormalStressNormalFlow');
model.physics('spf').feature('bs1').set('f0', 1, '0');
model.physics('spf').feature('bs1').selection.set([10]);
model.physics('spf').feature.create('bs2', 'BoundaryStress', 1);
model.physics('spf').feature('bs2').set('BoundaryCondition', 1, 'NormalStressNormalFlow');
model.physics('spf').feature('bs2').selection.set([1]);
model.physics('spf').feature('bs2').set('f0', 1, '1');
model.physics('spf').feature.create('vf1', 'VolumeForce', 2);
model.physics('spf').feature('vf1').selection.set([2]);
model.physics('spf').feature('vf1').set('F', {'-alpha*u' '-alpha*v' '0'});

model.physics.create('sens', 'Sensitivity', 'geom1');
model.physics('sens').feature.create('iobj1', 'IntegralObjective', 0);
model.physics('sens').feature('iobj1').selection.set([5]);
model.physics('sens').feature('iobj1').set('objectiveExpression', 1, 'u');
model.physics('sens').feature.create('cvar1', 'ControlVariableField', 2);
model.physics('sens').feature('cvar1').selection.set([2]);
model.physics('sens').feature('cvar1').set('fieldVariableName', 1, 'gamma');
model.physics('sens').feature('cvar1').set('order', 1, '1');
model.physics('sens').feature('cvar1').set('initialValue', 1, 'gamma_t');

%%% DEFINE STUDY
model.study.create('std1');
model.study('std1').feature.create('stat', 'Stationary');

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.create('s1', 'Stationary');
model.sol('sol1').feature('s1').feature.create('fc1', 'FullyCoupled');
model.sol('sol1').feature('s1').feature.create('i1', 'Iterative');
model.sol('sol1').feature('s1').feature.create('sn1', 'Sensitivity');
model.sol('sol1').feature('s1').feature('sn1').set('sensmethod', 'adjoint');
model.sol('sol1').feature('s1').feature('sn1').set('sensfunc', 'mod1.sens.iobj1');
model.sol('sol1').attach('std1');

model.variable.create('var1');
model.variable('var1').model('mod1');
model.variable('var1').set('alpha', 'alphamin+(alphamax-alphamin)*q*(1-gamma)/(q+gamma)');
model.variable('var1').set('gamma_t', 1);
model.variable('var1').set('q', 0.03);
model.variable('var1').set('alphamax', 1e4);
model.variable('var1').set('alphamin', 0);

%%% CALCULATE SOLUTION
model.sol('sol1').runAll;

%%%% PLOT RESULTS
model.result.create('pg1', 2);
model.result('pg1').feature.create('surf1', 'Surface');
model.result('pg1').feature('surf1').set('expr', 'spf.U');
model.result.create('pg2', 2);
model.result('pg2').feature.create('surf1', 'Surface');
model.result('pg2').feature('surf1').set('expr', 'gamma');

figure(1); subplot(2,1,1); mphplot(model,'pg1','rangenum',1); 
figure(1); subplot(2,1,2); mphplot(model,'pg2','rangenum',1); 

u0 = mphgetu(model);
xmi = model.sol('sol1').xmeshInfo;
names = cell(xmi.dofs.dofNames);
gname = 'mod1.gamma';
gidx = strmatch(gname, names)-1; % zero based
igamma = find(xmi.dofs.nameInds==gidx);

% SENSITIVITY ANALYSIS
phi = mphglobal(model, 'sens.iobj1');
dPhidgamma = mpheval(model, 'fsens(gamma)', 'selection', 2);

% Update u0. This is NOT the right expression!!!
u0(igamma) = u0(igamma)+dPhidgamma.d1';
% u0(igamma) = u0(igamma)-0.5;
model.sol('sol1').setU(u0);
model.sol('sol1').createSolution;

%%% CALCULATE SOLUTION
model.sol('sol1').feature('v1').set('initmethod', 'sol');
model.sol('sol1').feature('v1').set('initsol', 'sol1');
% model.sol('sol1').feature('v1').set('notsolmethod', 'sol');
model.sol('sol1').feature('v1').set('notsol', 'sol1');
% model.sol('sol1').feature('v1').set('storesol', 'u');
model.sol('sol1').runAll;

% Compare design field before and after solution (should be identical)
gamma_old = u0(igamma);
u0 = mphgetu(model);
gamma = u0(igamma);
gamma_old(1:10)
gamma(1:10)

%%%% PLOT NEW RESULTS 
figure(2); subplot(2,1,1); mphplot(model,'pg1','rangenum',1); 
figure(2); subplot(2,1,2); mphplot(model,'pg2','rangenum',1);