function IPMC_general_model(t_n,t_E,E_N,E_E,nu_N,nu_E,beta_N,wd,L,clamped,study,outputfile)
%t_N - thickness of Nafion membrane [m]
%t_E - thickness of electrodes [m]
%E_N - Young's modulus of Nafion membrane [Pa]
%E_E - Young's modulus of Nafion membrane [Pa]
%nu_N - Poisson ration of Nafion membrane [-]
%nu_E - Young's modulus of electrodes [-]
%beta_N - stress-concentration coupling coefficient [m^3 Pa/mol]
%wd - depth of the IPMC cantilever
%L - length of the free part of the IPMC cantilever
%clamped - 0: free tip; 1: fixed tip

%study - '1D3D': PNP solution from 1D extrapolated to 3D SP/'2D3D': PNP
%solution from 2D extrapolated to 3D SP;'2D2D': PNP solution from 2D
%inserted to 2D SP;'1D2D': PNP solution from 1D extrapolated to 2D SP

%outputfile - string of result outputfile name


%Study - 1D3D,1D2D,2D3D,2D
ModelUtil.showProgress(true)

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

model = ModelUtil.create('Model');

model.modelPath(tempdir);

model.name('IPMC_general_model.mph');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Global parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


model.param.set('D', '1e-10', 'diffusivity coeff [m2/s]');
model.param.set('Na', '6.022e23', 'mol-1');
model.param.set('e', '1.602e-19', 'element. charge [C]');
model.param.set('F', 'Na*e', 'Faraday const. [C mol-1]');
%model.param.set('epsilon', '25e-3', 'el. permitivity [F/m]');
model.param.set('epsilon', '2.8e-3', 'el. permitivity [F/m]');
%model.param.set('c0', '1200', 'ion concentration [mol/m3]');
model.param.set('c0', '1000', 'ion concentration [mol/m3]');
model.param.set('R', '8.31', '{J/Kg/mol}');
model.param.set('T', '293', '{k}');
model.param.set('um', 'D/R/T', 'mobility');
model.param.set('E_N', num2str(E_N), 'young modulus membrane [Pa]');
model.param.set('nu_N', num2str(nu_N), 'poisson membrane [-]');
model.param.set('E_E', num2str(E_E), 'young modulus electrodes [Pa]');
model.param.set('nu_E', num2str(nu_E), 'poisson electrodes [-]');
model.param.set('beta_N', num2str(beta_N));
model.param.set('G', 'E_N/(2*(1+nu_N))');
model.param.set('G2', 'E_E/(2*(1+nu_E))');
model.param.set('K', 'E_N/(3*(1-2*nu_N))');
model.param.set('K2', 'E_E/(3*(1-2*nu_E))');
model.param.set('KK', '1');
model.param.set('wd', num2str(wd));
model.param.set('L', num2str(L));
model.param.set('t_E', num2str(t_E));
model.param.set('t_N', num2str(t_N));
model.param.set('fix', '5e-3');
model.param.set('silomer', '0.001');
% model.param.set('Lambda_axi', 'E_N*nu_N/((1+nu_N)*(1-2*nu_N))');
% model.param.set('Lambda2_axi', 'E_E*nu_E/((1+nu_E)*(1-2*nu_E))');
% model.param.set('mu_axi', 'E_N/(2*(1+nu_N))');
% model.param.set('mu2_axi', 'E_E/(2*(1+nu_E))');
model.param.set('Var1D3D', '1');
model.param.set('Var1D2D', '1');

model.modelNode.create('mod1');
model.modelNode.create('mod2');
model.modelNode.create('mod3');
model.modelNode('mod1').name('2D model');
model.modelNode('mod2').name('3D Model');
model.modelNode('mod3').name('1D Model');

%Geometry model 1 - 2D
model.geom.create('geom1', 2);
model.geom('geom1').model('mod1');
model.geom('geom1').feature.create('r1', 'Rectangle');
model.geom('geom1').feature.create('r2', 'Rectangle');
model.geom('geom1').feature.create('r3', 'Rectangle');
model.geom('geom1').feature.create('b1', 'BezierPolygon');
model.geom('geom1').feature('r1').set('size', {'L' 't_E'});
model.geom('geom1').feature('r2').set('size', {'L' 't_N'});
model.geom('geom1').feature('r2').set('pos', {'0' 't_E'});
model.geom('geom1').feature('r3').set('size', {'L' 't_E'});
model.geom('geom1').feature('r3').set('pos', {'0' 't_N+t_E'});
model.geom('geom1').feature('b1').set('p', {'fix' 'fix'; '0' 't_N+2*t_E'});
model.geom('geom1').feature('b1').set('degree', {'1'});
model.geom('geom1').feature('b1').set('w', {'1' '1'});
model.geom('geom1').run;



%Geometry model 2 - 3D
model.geom.create('geom2', 3);
model.geom('geom2').model('mod2');
model.geom('geom2').feature.create('wp1', 'WorkPlane');
model.geom('geom2').feature('wp1').geom.feature.create('r1', 'Rectangle');
model.geom('geom2').feature('wp1').geom.feature.create('r2', 'Rectangle');
model.geom('geom2').feature('wp1').geom.feature.create('r3', 'Rectangle');
model.geom('geom2').feature('wp1').geom.feature.create('b1', 'BezierPolygon');
model.geom('geom2').feature('wp1').geom.feature.create('uni1', 'Union');
model.geom('geom2').feature.create('ext1', 'Extrude');
model.geom('geom2').feature('wp1').geom.feature('r1').set('size', {'L' 't_E'});
model.geom('geom2').feature('wp1').geom.feature('r2').set('size', {'L' 't_N'});
model.geom('geom2').feature('wp1').geom.feature('r2').set('pos', {'0' 't_E'});
model.geom('geom2').feature('wp1').geom.feature('r3').set('size', {'L' 't_E'});
model.geom('geom2').feature('wp1').geom.feature('r3').set('pos', {'0' 't_E+t_N'});
model.geom('geom2').feature('wp1').geom.feature('b1').set('degree', {'1'});
model.geom('geom2').feature('wp1').geom.feature('b1').set('p', {'fix' 'fix'; '0' '2*t_E+t_N'});
model.geom('geom2').feature('wp1').geom.feature('b1').set('w', {'1' '1'});
model.geom('geom2').feature('wp1').geom.feature('uni1').selection('input').set({'b1' 'r1' 'r2' 'r3'});
model.geom('geom2').feature('ext1').setIndex('distance', 'wd/2', 0);
model.geom('geom2').feature('ext1').selection('input').set({'wp1'});
model.geom('geom2').run;
model.geom('geom2').feature.create('mir1', 'Mirror');
model.geom('geom2').feature('mir1').set('keep', 'on');
model.geom('geom2').feature('mir1').selection('input').set({'ext1'});
model.geom('geom2').feature('mir1').setIndex('axis', '-1', 2);
model.geom('geom2').run('mir1');
model.geom('geom2').run;
model.geom('geom2').run('fin');


%Geometry model 3 - 1D
model.geom.create('geom3', 1);
model.geom('geom3').feature.create('i1', 'Interval');
model.geom('geom3').feature.create('i2', 'Interval');
model.geom('geom3').feature.create('i3', 'Interval');
model.geom('geom3').feature('i1').set('p2', 't_E');
model.geom('geom3').feature('i2').set('p1', 't_E');
model.geom('geom3').feature('i2').set('p2', 't_E+t_N');
model.geom('geom3').feature('i3').set('p1', 't_E+t_N');
model.geom('geom3').feature('i3').set('p2', 't_E+t_N+t_E');
model.geom('geom3').run;


% model.selection.create('sel1', 'Explicit');
% model.selection('sel1').model('mod1');
% model.selection('sel1').geom('geom1', 1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Global Variables definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.variable.create('var1');
%Electrical excitation definition
model.variable('var1').set('V_fce', '3*sin(0.2*pi*t)');
%1D to 3D Extrusion of the solution of concentration computed in 1D model 3
model.variable('var1').set('c_mod1D3D', 'mod3.linext2(mod3.c2)');
%2D to 3D Extrusion of the solution of concentration computed in 2D model 1
model.variable('var1').set('c_mod2D3D', 'mod1.linext1(mod1.c)');
%1D to 2D Extrusion of the solution of concentration computed in 1D model 3
model.variable('var1').set('c_mod1D2D', 'mod3.linext5(mod3.c2)');
%1D to 3D Extrusion Electrostatic force computed in 1D model 3
model.variable('var1').set('Felstat1D3D', 'mod3.linext3(mod3.Felstat1D)');
%2D to 3D Extrusion Electrostatic force computed in 2D model 1
model.variable('var1').set('Felstat2D3D', 'mod1.linext4(mod1.Felstat)');
%1D to 2D Extrusion Electrostatic force computed in 1D model 3
model.variable('var1').set('Felstat1D2D', 'mod3.linext6(mod3.Felstat1D)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2D Model 1 Variables definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.variable.create('var5');
model.variable('var5').model('mod1');
model.variable('var5').set('Sx', '4/3*G*ux-2/3*G*vy-p');
model.variable('var5').set('Sy', '4/3*G*vy-2/3*G*ux-p');
model.variable('var5').set('Sz', '-2/3*G*ux-2/3*G*vy-p');
model.variable('var5').set('Sxy', 'G*(uy+vx)');
model.variable('var5').set('p', '-K*(ux+vy-3*beta_N/K*(c*(1-Var1D2D)+c_mod1D2D*Var1D2D-c0))');
model.variable('var5').set('Felstat', '-KK*(F/2)*(abs(c-c0)-(c-c0))*phiy*wd');
model.variable('var5').selection.geom('geom1', 2);
model.variable('var5').selection.set([2 5]);
model.variable.create('var6');
model.variable('var6').model('mod1');
model.variable('var6').set('Sx', '4/3*G2*ux-2/3*G2*vy-p2');
model.variable('var6').set('Sy', '4/3*G2*vy-2/3*G2*ux-p2');
model.variable('var6').set('Sz', '-2/3*G2*ux-2/3*G2*vy-p2');
model.variable('var6').set('Sxy', 'G2*(uy+vx)');
model.variable('var6').set('p2', '-K2*(ux+vy)');
model.variable('var6').selection.geom('geom1', 2);
model.variable('var6').selection.set([1 3 4 6]);
model.variable.create('var7');
model.variable('var7').model('mod1');
model.variable('var7').set('F_pom', 'intop1(Felstat)');
model.variable('var7').name('Variables 7a');

model.variable.create('var8');
model.variable('var8').model('mod1');
model.variable('var8').set('Felstat', '-F_pom/(t_E)*(V_fce<0)');
model.variable('var8').selection.geom('geom1', 2);
model.variable('var8').selection.set([1 4]);
model.variable.create('var9');
model.variable('var9').model('mod1');
model.variable('var9').set('Felstat', '-F_pom/(t_E)*(V_fce>0)');
model.variable('var9').selection.geom('geom1', 2);
model.variable('var9').selection.set([3 6]);
model.variable.create('var12');
model.variable('var12').model('mod1');
model.variable('var12').set('c_mod1', 'linext1(c)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3D Model 2 Variables definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.variable.create('var10');
model.variable('var10').model('mod2');
% model.variable('var10').set('Sxx', 'Lambda_axi*(u3x+v3y+w3z)+2*mu_axi*u3x-p');
% model.variable('var10').set('Syy', 'Lambda_axi*(u3x+v3y+w3z)+2*mu_axi*v3y-p');
% model.variable('var10').set('Szz', 'Lambda_axi*(u3x+v3y+w3z)+2*mu_axi*w3z-p');
model.variable('var10').set('Sxx', '4/3*G*u3x-2/3*G*(v3y+w3z)-p');
model.variable('var10').set('Syy', '4/3*G*v3y-2/3*G*(u3x+w3z)-p');
model.variable('var10').set('Szz', '4/3*G*w3z-2/3*G*(u3x+v3y)-p');
model.variable('var10').set('Sxy', '2*G*(u3y+v3x)*1/2');
model.variable('var10').set('Syz', '2*G*(v3z+w3y)*1/2');
model.variable('var10').set('Szx', '2*G*(w3x+u3z)*1/2');
model.variable('var10').set('p', '-K*(u3x+v3y+w3z-3*beta_N/K*(c_mod1D3D*Var1D3D+c_mod2D3D*(1-Var1D3D)-c0))');
model.variable('var10').selection.geom('geom2', 3);
%model.variable('var10').selection.set([2 5]);
model.variable('var10').selection.set([3 4 9 10]);
model.variable('var10').name('Variables 10a');

model.variable.create('var13');
model.variable('var13').model('mod2');
% model.variable('var13').set('Sxx', 'Lambda2_axi*(u3x+v3y+w3z)+2*mu2_axi*u3x');
% model.variable('var13').set('Syy', 'Lambda2_axi*(u3x+v3y+w3z)+2*mu2_axi*v3y');
% model.variable('var13').set('Szz', 'Lambda2_axi*(u3x+v3y+w3z)+2*mu2_axi*w3z');
% model.variable('var13').set('Sxy', '2*mu2_axi*(u3y+v3x)*1/2');
% model.variable('var13').set('Syz', '2*mu2_axi*(v3z+w3y)*1/2');
% model.variable('var13').set('Szx', '2*mu2_axi*(w3x+u3z)*1/2');
model.variable('var13').set('Sxx', '4/3*G2*u3x-2/3*G2*(v3y+w3z)-p2');
model.variable('var13').set('Syy', '4/3*G2*v3y-2/3*G2*(u3x+w3z)-p2');
model.variable('var13').set('Szz', '4/3*G2*w3z-2/3*G2*(u3x+v3y)-p2');
model.variable('var13').set('Sxy', '2*G2*(u3y+v3x)*1/2');
model.variable('var13').set('Syz', '2*G2*(v3z+w3y)*1/2');
model.variable('var13').set('Szx', '2*G2*(w3x+u3z)*1/2');
model.variable('var13').set('p2', '-K2*(u3x+v3y+w3z)');

model.variable('var13').selection.geom('geom2', 3);
%model.variable('var13').selection.set([1 3 4 6]);
model.variable('var13').selection.set([1 2 5 6 7 8 11 12]);

model.variable.create('var18');
model.variable('var18').model('mod2');
model.variable('var18').selection.geom('geom2', 1);
model.variable('var18').selection.set([75]);
model.variable('var18').selection.global;
model.variable('var18').set('F_r', 'intop3(vb_lm3D)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1D Model 1 Variables definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.variable.create('var14');
model.variable('var14').model('mod3');
model.variable('var14').set('F_pom1D', 'intop2(Felstat1D)');
model.variable('var14').name('Variables 14a');
model.variable.create('var15');
model.variable('var15').model('mod3');
model.variable('var15').set('Felstat1D', '-F_pom1D/(t_E)*(V_fce<0)');
model.variable('var15').selection.geom('geom3', 1);
model.variable('var15').selection.set([1]);
model.variable.create('var16');
model.variable('var16').model('mod3');
model.variable('var16').set('Felstat1D', '-F_pom1D/(t_E)*(V_fce>0)');
model.variable('var16').selection.geom('geom3', 1);
model.variable('var16').selection.set([3]);
model.variable.create('var17');
model.variable('var17').model('mod3');
model.variable('var17').set('Felstat1D', '-KK*(F/2)*(abs(c2-c0)-(c2-c0))*phi2x*wd');
model.variable('var17').selection.geom('geom3', 1);
model.variable('var17').selection.set([2]);

% model.view('view4').tag('view6');
% model.view.create('view4', 3);
% model.view.create('view5', 2);
% model.view.create('view7', 3);
% model.view.create('view9', 3);

% model.material.create('mat1');
% model.material('mat1').model('mod2');
% model.material('mat1').propertyGroup.create('Enu', 'Young''s modulus and Poisson''s ratio');
% model.material('mat1').propertyGroup.create('Murnaghan', 'Murnaghan');
% model.material('mat1').propertyGroup.create('Lame', ['Lam' native2unicode(hex2dec('00e9'), 'Cp1252') ' constants']);
% model.material('mat1').propertyGroup.create('NeoHookean', 'Neo-Hookean');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Extrusion operators 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.cpl.create('intop1', 'Integration', 'geom1');
model.cpl('intop1').selection.geom('geom1', 1);
model.cpl('intop1').selection.set([3]);

model.cpl.create('intop3', 'Integration', 'geom2');
model.cpl('intop3').selection.geom('geom2', 1);
model.cpl('intop3').selection.set([75]);


model.cpl.create('linext1', 'LinearExtrusion', 'geom1');
model.cpl('linext1').selection.set([1 2 3 4 5 6]);
model.cpl.create('linext2', 'LinearExtrusion', 'geom3');
model.cpl('linext2').selection.set([2]);
model.cpl.create('intop2', 'Integration', 'geom3');
model.cpl('intop2').selection.set([2]);
model.cpl.create('linext3', 'LinearExtrusion', 'geom3');
model.cpl('linext3').selection.all;
model.cpl.create('linext4', 'LinearExtrusion', 'geom1');
model.cpl('linext4').selection.all;
model.cpl.create('linext5', 'LinearExtrusion', 'geom3');
model.cpl('linext5').selection.set([2]);
model.cpl.create('linext6', 'LinearExtrusion', 'geom3');
model.cpl('linext6').selection.set([1 2 3]);



model.cpl('linext1').set('usesrcmap', 'on');
model.cpl('linext1').set('dstgeom', 'geom2');
model.cpl('linext1').selection('dstvertex1').set([6]);
model.cpl('linext1').selection('dstvertex2').set([30]);
model.cpl('linext1').selection('dstvertex3').set([9]);
model.cpl('linext1').selection('srcvertex1').geom('geom1', 0);
model.cpl('linext1').selection('srcvertex1').set([2]);
model.cpl('linext1').selection('srcvertex2').geom('geom1', 0);
model.cpl('linext1').selection('srcvertex2').set([10]);
model.cpl('linext1').selection('srcvertex3').geom('geom1', 0);
model.cpl('linext1').selection('srcvertex3').set([3]);
model.cpl('linext1').selection('srcvertex4').geom('geom1', 0);
model.cpl('linext2').set('dstgeom', 'geom2');
model.cpl('linext2').selection('dstvertex1').set([18]);
model.cpl('linext2').selection('dstvertex2').set([21]);
model.cpl('linext2').selection('srcvertex1').geom('geom3', 0);
model.cpl('linext2').selection('srcvertex1').set([2]);
model.cpl('linext2').selection('srcvertex2').geom('geom3', 0);
model.cpl('linext2').selection('srcvertex2').set([3]);
model.cpl('linext2').selection('srcvertex3').geom('geom3', 0);
model.cpl('linext2').selection('srcvertex4').geom('geom3', 0);
model.cpl('linext3').set('dstgeom', 'geom2');
model.cpl('linext3').selection('dstvertex1').set([15]);
model.cpl('linext3').selection('dstvertex2').set([24]);
model.cpl('linext3').selection('srcvertex1').geom('geom3', 0);
model.cpl('linext3').selection('srcvertex1').set([1]);
model.cpl('linext3').selection('srcvertex2').geom('geom3', 0);
model.cpl('linext3').selection('srcvertex2').set([4]);
model.cpl('linext3').selection('srcvertex3').geom('geom3', 0);
model.cpl('linext3').selection('srcvertex4').geom('geom3', 0);
model.cpl('linext4').set('dstgeom', 'geom2');
model.cpl('linext4').selection('dstvertex1').set([3]);
model.cpl('linext4').selection('dstvertex2').set([27]);
model.cpl('linext4').selection('dstvertex3').set([12]);
% model.cpl('linext4').selection('dstvertex1').set([2]);
% model.cpl('linext4').selection('dstvertex2').set([18]);
% model.cpl('linext4').selection('dstvertex3').set([8]);
model.cpl('linext4').selection('srcvertex1').geom('geom1', 0);
model.cpl('linext4').selection('srcvertex1').set([1]);
model.cpl('linext4').selection('srcvertex2').geom('geom1', 0);
model.cpl('linext4').selection('srcvertex2').set([9]);
model.cpl('linext4').selection('srcvertex3').geom('geom1', 0);
model.cpl('linext4').selection('srcvertex3').set([4]);
model.cpl('linext4').selection('srcvertex4').geom('geom1', 0);
model.cpl('linext5').set('dstgeom', 'geom1');
model.cpl('linext5').selection('srcvertex1').geom('geom3', 0);
model.cpl('linext5').selection('srcvertex1').set([2]);
model.cpl('linext5').selection('srcvertex2').geom('geom3', 0);
model.cpl('linext5').selection('srcvertex2').set([3]);
model.cpl('linext5').selection('srcvertex3').geom('geom3', 0);
model.cpl('linext5').selection('srcvertex4').geom('geom3', 0);
model.cpl('linext5').selection('dstvertex1').set([6]);
model.cpl('linext5').selection('dstvertex2').set([7]);
model.cpl('linext6').set('dstgeom', 'geom1');
model.cpl('linext6').selection('srcvertex1').geom('geom3', 0);
model.cpl('linext6').selection('srcvertex1').set([1]);
model.cpl('linext6').selection('srcvertex2').geom('geom3', 0);
model.cpl('linext6').selection('srcvertex2').set([4]);
model.cpl('linext6').selection('srcvertex3').geom('geom3', 0);
model.cpl('linext6').selection('srcvertex4').geom('geom3', 0);
model.cpl('linext6').selection('dstvertex1').set([5]);
model.cpl('linext6').selection('dstvertex2').set([8]);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2D model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Nernst-Planck PDE
model.physics.create('g', 'GeneralFormPDE', 'geom1');
model.physics('g').field('dimensionless').field('c');
model.physics('g').field('dimensionless').component({'c'});
model.physics('g').selection.set([2 5]);

model.physics('g').name('Nernst-Planck 2D');
model.physics('g').prop('ShapeProperty').set('boundaryFlux', '0');
model.physics('g').feature('gfeq1').set('f', '0');
model.physics('g').feature('gfeq1').set('Ga', {'-D*cx-um*F*c*phix' '-D*cy-um*F*c*phiy'});
model.physics('g').feature('gfeq1').name('Nernst-Planck PDE 2D');
model.physics('g').feature('init1').set('c', 'c0');

%Poisson Electrostatics PDE
model.physics.create('g2', 'GeneralFormPDE', 'geom1');
model.physics('g2').identifier('g2');
model.physics('g2').field('dimensionless').field('phi');
model.physics('g2').field('dimensionless').component({'phi'});
model.physics('g2').feature.create('gfeq2', 'GeneralFormPDE', 2);
model.physics('g2').feature('gfeq2').selection.set([1 3 4 6]);
model.physics('g2').feature.create('cons1', 'Constraint', 1);
model.physics('g2').feature('cons1').selection.set([2 9]);
model.physics('g2').feature.create('cons2', 'Constraint', 1);
model.physics('g2').feature('cons2').selection.set([7 14]);

model.physics('g2').name('Poisson electrostatics 2D');
model.physics('g2').prop('ShapeProperty').set('boundaryFlux', '0');
model.physics('g2').feature('gfeq1').set('f', 'F/epsilon*(c-c0)');
model.physics('g2').feature('gfeq1').set('Ga', {'-phix' '-phiy'});
model.physics('g2').feature('gfeq1').set('da', '0');
model.physics('g2').feature('gfeq1').name('Poisson electrostatics PDE 2D');
model.physics('g2').feature('gfeq2').set('f', '0');
model.physics('g2').feature('gfeq2').set('da', '0');
model.physics('g2').feature('cons1').set('R', '-phi');
model.physics('g2').feature('cons2').set('R', '-phi+V_fce');

%Structural mechanics PDE
model.physics.create('g3', 'GeneralFormPDE', 'geom1');
model.physics('g3').field('dimensionless').component({'u' 'v'});
model.physics('g3').feature.create('dir1', 'DirichletBoundary', 1);
model.physics('g3').feature('dir1').selection.set([2 7]);

model.physics('g3').name('Structural mechanics 2D');
model.physics('g3').feature('gfeq1').set('f', {'0'; 'Felstat*(1-Var1D2D)+Felstat1D2D*Var1D2D'});
model.physics('g3').feature('gfeq1').set('Ga', {'Sx' 'Sxy'; 'Sxy' 'Sy'});
model.physics('g3').feature('gfeq1').set('da', {'0'; '0'; '0'; '0'});
model.physics('g3').feature('gfeq1').name('Structural mechanics PDE 2D');

if clamped
    model.physics('g3').feature.create('wc1', 'WeakConstraint', 0);
    model.physics('g3').feature('wc1').selection.set([12]);
    model.physics('g3').feature('wc1').set('constraintExpression', 1, 'v');
    model.physics('g3').feature('wc1').set('fieldVariableName', 1, 'vb_lm2D');
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3D model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Structural mechanics PDE
model.physics.create('g4', 'GeneralFormPDE', 'geom2');
model.physics('g4').field('dimensionless').field('u3');
model.physics('g4').field('dimensionless').component({'u3' 'v3' 'w3'});
model.physics('g4').feature.create('dir1', 'DirichletBoundary', 2);
%model.physics('g4').feature('dir1').selection.set([2 13]);
model.physics('g4').feature('dir1').selection.set([2 5 22 23]);

model.physics('g4').feature.create('dir2', 'DirichletBoundary', 2);
model.physics('g4').feature('dir2').selection.set([27 28 29]);

model.physics('g4').name('Structural mechanics 3D');
model.physics('g4').feature('gfeq1').set('f', {'0'; 'Felstat1D3D*Var1D3D+Felstat2D3D*(1-Var1D3D)'; '0'});
model.physics('g4').feature('gfeq1').set('Ga', {'Sxx' 'Sxy' 'Szx'; 'Sxy' 'Syy' 'Syz'; 'Szx' 'Syz' 'Szz'});
model.physics('g4').feature('gfeq1').set('da', {'0'; '0'; '0'; '0'; '0'; '0'; '0'; '0'; '0'});
model.physics('g4').feature('gfeq1').name('Structural mechanics PDE 3D');
model.physics('g4').feature('dir2').selection.active(false);
model.physics('g4').feature('dir2').set('r', {'0'; '1e-8*t'; '0'});
model.physics('g4').feature('dir2').set('useDirichletCondition', {'0'; '1'; '0'});
model.physics('g4').feature('dir2').active(false);

if clamped==1
    model.physics('g4').feature.create('wc1', 'WeakConstraint', 1);
    model.physics('g4').feature('wc1').selection.set([75]);
    model.physics('g4').feature('wc1').set('fieldVariableName', 1, 'vb_lm3D');
    model.physics('g4').feature('wc1').set('constraintExpression', 1, 'v3');
elseif clamped == 2
    model.physics('g4').feature.create('wc2', 'WeakConstraint', 0);
    model.physics('g4').feature('wc2').set('constraintExpression', 1, 'v3');
    model.physics('g4').feature('wc2').set('fieldVariableName', 1, 'vb_lm3Dp');
    model.physics('g4').feature('wc2').selection.set([35]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1D model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Poisson Electrostatics PDE
model.physics.create('g5', 'GeneralFormPDE', 'geom3');
model.physics('g5').field('dimensionless').field('phi2');
model.physics('g5').field('dimensionless').component({'phi2'});
model.physics('g5').feature.create('gfeq2', 'GeneralFormPDE', 1);
model.physics('g5').feature('gfeq2').selection.set([1 3]);
model.physics('g5').feature.create('cons1', 'Constraint', 0);
model.physics('g5').feature('cons1').selection.set([1]);
model.physics('g5').feature.create('cons2', 'Constraint', 0);
model.physics('g5').feature('cons2').selection.set([4]);

model.physics('g5').name('Poisson Electrostatics 1D');
model.physics('g5').feature('gfeq1').set('f', 'F/epsilon*(c2-c0)');
model.physics('g5').feature('gfeq1').set('Ga', '-phi2x');
model.physics('g5').feature('gfeq1').set('da', '0');
model.physics('g5').feature('gfeq1').name('Poisson electrostatics PDE 1D');
model.physics('g5').feature('gfeq2').set('f', '0');
model.physics('g5').feature('gfeq2').set('da', '0');
model.physics('g5').feature('cons1').set('R', '-phi2');
model.physics('g5').feature('cons2').set('R', '-phi2+V_fce');

%Nernst-Planck PDE
model.physics.create('g7', 'GeneralFormPDE', 'geom3');
model.physics('g7').identifier('g7');
model.physics('g7').field('dimensionless').field('c2');
model.physics('g7').field('dimensionless').component({'c2'});
model.physics('g7').selection.set([2]);

model.physics('g7').name('Nernst-Planck PDE 1D');
model.physics('g7').feature('gfeq1').set('f', '0');
model.physics('g7').feature('gfeq1').set('Ga', '-D*c2x-um*F*c2*phi2x');
model.physics('g7').feature('gfeq1').name('Nernst-Planck 1D');
model.physics('g7').feature('init1').set('c2', 'c0');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Mesh generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2D model
model.mesh.create('mesh2', 'geom1');
model.mesh('mesh2').feature.create('map1', 'Map');
model.mesh('mesh2').feature('map1').feature.create('dis1', 'Distribution');
model.mesh('mesh2').feature('map1').feature('dis1').selection.set([2 4 6 7]);
model.mesh('mesh2').feature('map1').feature.create('dis2', 'Distribution');
model.mesh('mesh2').feature('map1').feature('dis2').selection.set([9 11 13 14]);
model.mesh('mesh2').feature('map1').feature.create('dis3', 'Distribution');
model.mesh('mesh2').feature('map1').feature('dis3').selection.set([3 10 16]);
model.mesh('mesh2').feature('map1').feature.create('dis4', 'Distribution');
model.mesh('mesh2').feature('map1').feature('dis4').selection.set([1 5 8 12 15 17]);
model.mesh('mesh2').feature('map1').feature.create('dis5', 'Distribution');

model.mesh('mesh2').name('Mesh 2D');
model.mesh('mesh2').feature('map1').feature('dis1').set('numelem', '10');
model.mesh('mesh2').feature('map1').feature('dis2').set('numelem', num2str(L/10e-3*50));
model.mesh('mesh2').feature('map1').feature('dis3').set('elemratio', '30');
model.mesh('mesh2').feature('map1').feature('dis3').set('elemcount', '70');
model.mesh('mesh2').feature('map1').feature('dis3').set('symmetric', true);
model.mesh('mesh2').feature('map1').feature('dis3').set('type', 'predefined');
model.mesh('mesh2').feature('map1').feature('dis5').set('numelem', '20');
model.mesh('mesh2').run;

%3D model


model.mesh.create('mesh3', 'geom2');
model.mesh('mesh3').feature.create('map1', 'Map');
model.mesh('mesh3').feature('map1').selection.set([3 10 17 26 33 40]);
model.mesh('mesh3').feature('map1').feature.create('dis2', 'Distribution');
model.mesh('mesh3').feature('map1').feature('dis2').selection.set([3 11 19 26]);
model.mesh('mesh3').feature('map1').feature.create('dis3', 'Distribution');
model.mesh('mesh3').feature('map1').feature('dis3').selection.set([32 40 48 55]);
model.mesh('mesh3').feature('map1').feature.create('dis4', 'Distribution');
model.mesh('mesh3').feature('map1').feature('dis4').selection.set([10 39 65]);
model.mesh('mesh3').feature('map1').feature.create('dis5', 'Distribution');
model.mesh('mesh3').feature('map1').feature('dis5').selection.set([2 18 31 47 60 70]);

model.mesh('mesh3').feature.create('swe1', 'Sweep');
model.mesh('mesh3').feature('swe1').feature.create('dis1', 'Distribution');

model.mesh('mesh3').name('Mesh 3D');
model.mesh('mesh3').feature('map1').feature('dis2').set('numelem', '10');
model.mesh('mesh3').feature('map1').feature('dis3').set('numelem', num2str(L/10e-3*50));
model.mesh('mesh3').feature('map1').feature('dis4').set('elemratio', '30');
model.mesh('mesh3').feature('map1').feature('dis4').set('elemcount', '70');
model.mesh('mesh3').feature('map1').feature('dis4').set('symmetric', true);
model.mesh('mesh3').feature('map1').feature('dis4').set('type', 'predefined');
model.mesh('mesh3').feature('swe1').selection('sourceface').set([3 10 17 26 33 40]);
model.mesh('mesh3').feature('swe1').feature('dis1').set('numelem', '3');
model.mesh('mesh3').run('swe1');

model.mesh('mesh3').run;


% model.mesh.create('mesh3', 'geom2');
% model.mesh('mesh3').feature.create('map1', 'Map');
% model.mesh('mesh3').feature('map1').selection.set([4 8 12 17 21 25]);
% model.mesh('mesh3').feature('map1').feature.create('dis2', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis2').selection.set([5 10 15 18]);
% model.mesh('mesh3').feature('map1').feature.create('dis3', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis3').selection.set([23 28 33 36]);
% model.mesh('mesh3').feature('map1').feature.create('dis4', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis4').selection.set([9 22 27 32 42]);
% model.mesh('mesh3').feature('map1').feature.create('dis5', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis5').selection.set([4 14 22 32 39 45]);
% model.mesh('mesh3').feature('map1').feature.create('dis6', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis6').selection.set([15]);
% model.mesh('mesh3').feature('map1').feature.create('dis7', 'Distribution');
% model.mesh('mesh3').feature('map1').feature('dis7').selection.set([33]);
% model.mesh('mesh3').feature.create('swe1', 'Sweep');
% model.mesh('mesh3').feature('swe1').feature.create('dis1', 'Distribution');
% 
% model.mesh('mesh3').name('Mesh 3D');
% model.mesh('mesh3').feature('map1').feature('dis2').set('numelem', '10');
% model.mesh('mesh3').feature('map1').feature('dis3').set('numelem', '40');
% model.mesh('mesh3').feature('map1').feature('dis4').set('elemratio', '30');
% model.mesh('mesh3').feature('map1').feature('dis4').set('elemcount', '70');
% model.mesh('mesh3').feature('map1').feature('dis4').set('symmetric', true);
% model.mesh('mesh3').feature('map1').feature('dis4').set('type', 'predefined');
% model.mesh('mesh3').feature('map1').feature('dis6').active(false);
% model.mesh('mesh3').feature('map1').feature('dis6').set('numelem', '10');
% model.mesh('mesh3').feature('map1').feature('dis7').active(false);
% model.mesh('mesh3').feature('map1').feature('dis7').set('numelem', '50');
% model.mesh('mesh3').run;


%1D model
model.mesh.create('mesh4', 'geom3');
model.mesh('mesh4').feature.create('edg1', 'Edge');
model.mesh('mesh4').feature('edg1').selection.geom('geom3', 1);
model.mesh('mesh4').feature('edg1').selection.set([2]);
model.mesh('mesh4').feature('edg1').feature.create('dis1', 'Distribution');
model.mesh('mesh4').feature('edg1').feature('dis1').selection.set([1 2 3]);
model.mesh('mesh4').feature.create('edg2', 'Edge');
model.mesh('mesh4').feature('edg2').feature.create('dis1', 'Distribution');
model.mesh('mesh4').feature('edg2').feature('dis1').selection.set([1 3]);

model.mesh('mesh4').name('Mesh 1D');
model.mesh('mesh4').feature('edg1').feature('dis1').set('elemratio', '30');
model.mesh('mesh4').feature('edg1').feature('dis1').set('elemcount', '70');
model.mesh('mesh4').feature('edg1').feature('dis1').set('symmetric', true);
model.mesh('mesh4').feature('edg1').feature('dis1').set('type', 'predefined');
model.mesh('mesh4').feature('edg2').feature('dis1').set('numelem', '20');
model.mesh('mesh4').run;



switch study

    case '1D3D'
        
        model.study.create('std10');
        model.study('std10').feature.create('time', 'Transient');
        
        model.sol.create('sol10');
        model.sol('sol10').study('std10');
        model.sol('sol10').attach('std10');
        model.sol('sol10').feature.create('st1', 'StudyStep');
        model.sol('sol10').feature.create('v1', 'Variables');
        model.sol('sol10').feature.create('t1', 'Time');
        model.sol('sol10').feature('t1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol10').feature('t1').feature.create('se1', 'Segregated');
        model.sol('sol10').feature('t1').feature.remove('fcDef');
        
        model.sol('sol10').attach('std10');
        model.sol('sol10').name('Solver 1D 2D');
        model.sol('sol10').feature('st1').name('Compile Equations: Time Dependent');
        model.sol('sol10').feature('st1').set('studystep', 'time');
        model.sol('sol10').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_c').set('solvefor', false);
        if clamped == 1
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        
        model.sol('sol10').feature('t1').set('tlist', 'range(0,0.1,5)');
        model.sol('sol10').feature('t1').feature('se1').feature('ssDef').set('segvar', {'mod3_c2' 'mod3_phi2'});
        
        
        
        model.study.create('std11');
        model.study('std11').feature.create('param', 'Parametric');
        model.study('std11').feature.create('stat', 'Stationary');
        model.sol.create('sol11');
        model.sol('sol11').study('std11');
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature.create('st1', 'StudyStep');
        model.sol('sol11').feature.create('v1', 'Variables');
        model.sol('sol11').feature.create('s1', 'Stationary');
        model.sol('sol11').feature('s1').feature.create('p1', 'Parametric');
        model.sol('sol11').feature('s1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol11').feature('s1').feature.create('p2', 'Parametric');
        model.sol('sol11').feature('s1').feature.remove('fcDef');
        
        model.study('std11').feature('param').set('plistarr', {'1'});
        model.study('std11').feature('param').set('pname', {'Var1D3D'});
        model.study('std11').feature('param').set('plot', 'on');
        model.study('std11').feature('stat').set('notsolmethod', 'sol');
        model.study('std11').feature('stat').set('usesol', 'on');
        model.study('std11').feature('stat').set('activate', {'g' 'off' 'g2' 'off' 'g3' 'off' 'g4' 'on' 'g5' 'off'  ...
            'g7' 'off'});
        model.study('std11').feature('stat').set('notstudy', 'std10');
        
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature('st1').name('Compile Equations: Stationary');
        model.sol('sol11').feature('st1').set('study', 'std11');
        model.sol('sol11').feature('v1').set('notsolmethod', 'sol');
        model.sol('sol11').feature('v1').set('notsol', 'sol10');
        model.sol('sol11').feature('v1').set('control', 'stat');
        model.study('std11').feature('stat').set('notsolnum', 'all');
        model.sol('sol11').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_c').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_c2').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_phi2').set('solvefor', false);
        if clamped ==1 
            model.sol('sol11').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol11').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
%            model.sol('sol11').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);          
        end
        % if clamped
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v','mod1_vb_lm2D'});
        % else
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v'});
        % end
        %
        
        model.sol('sol11').feature('s1').set('control', 'stat');
        model.sol('sol11').feature('s1').feature('p1').set('plistarr', {'1'});
        model.sol('sol11').feature('s1').feature('p1').set('pname', {'Var1D3D'});
        model.sol('sol11').feature('s1').feature('p1').set('control', 'param');
        model.sol('sol11').feature('s1').feature('p1').set('plot', 'on');
        model.sol('sol11').feature('s1').feature('p2').set('plistarr', {'range(0,0.1,5)'});
        model.sol('sol11').feature('s1').feature('p2').set('pname', {'t'});
        model.sol('sol11').feature('s1').feature('p2').set('plot', 'on');        


    case '2D3D'
        %2D-3D Study
        model.param.set('Var1D3D', '0');
        
        model.study.create('std10');
        model.study('std10').feature.create('time', 'Transient');
        
        model.sol.create('sol10');
        model.sol('sol10').study('std10');
        model.sol('sol10').attach('std10');
        model.sol('sol10').feature.create('st1', 'StudyStep');
        model.sol('sol10').feature.create('v1', 'Variables');
        model.sol('sol10').feature.create('t1', 'Time');
        model.sol('sol10').feature('t1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol10').feature('t1').feature.create('se1', 'Segregated');
        model.sol('sol10').feature('t1').feature.remove('fcDef');
        
        model.sol('sol10').attach('std10');
        model.sol('sol10').name('Solver 2D 3D');
        model.sol('sol10').feature('st1').name('Compile Equations: Time Dependent');
        model.sol('sol10').feature('st1').set('studystep', 'time');
        model.sol('sol10').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_phi2').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_c2').set('solvefor', false);
%         if clamped
%             model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
%             model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
%         end
        if clamped == 1
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        
        
        model.sol('sol10').feature('t1').set('tlist', 'range(0,0.1,5)');
        model.sol('sol10').feature('t1').feature('se1').feature('ssDef').set('segvar', {'mod1_c' 'mod1_phi'});
        
        
        
        model.study.create('std11');
        model.study('std11').feature.create('param', 'Parametric');
        model.study('std11').feature.create('stat', 'Stationary');
        model.sol.create('sol11');
        model.sol('sol11').study('std11');
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature.create('st1', 'StudyStep');
        model.sol('sol11').feature.create('v1', 'Variables');
        model.sol('sol11').feature.create('s1', 'Stationary');
        model.sol('sol11').feature('s1').feature.create('p1', 'Parametric');
        model.sol('sol11').feature('s1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol11').feature('s1').feature.create('p2', 'Parametric');
        model.sol('sol11').feature('s1').feature.remove('fcDef');
        
        model.study('std11').feature('param').set('plistarr', {'0'});
        model.study('std11').feature('param').set('pname', {'Var1D3D'});
        model.study('std11').feature('param').set('plot', 'on');
        model.study('std11').feature('stat').set('notsolmethod', 'sol');
        model.study('std11').feature('stat').set('usesol', 'on');
        model.study('std11').feature('stat').set('activate', {'g' 'off' 'g2' 'off' 'g3' 'off' 'g4' 'on' 'g5' 'off'  ...
            'g7' 'off'});
        model.study('std11').feature('stat').set('notstudy', 'std10');
        
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature('st1').name('Compile Equations: Stationary');
        model.sol('sol11').feature('st1').set('study', 'std11');
        model.sol('sol11').feature('v1').set('notsolmethod', 'sol');
        model.sol('sol11').feature('v1').set('notsol', 'sol10');
        model.sol('sol11').feature('v1').set('control', 'stat');
        model.study('std11').feature('stat').set('notsolnum', 'all');
        model.sol('sol11').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_c').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_c2').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_phi2').set('solvefor', false);
        if clamped ==1 
            model.sol('sol11').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol11').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        elseif clamped == 2
            model.sol('sol11').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            %model.sol('sol11').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);          
        end
        % if clamped
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v','mod1_vb_lm2D'});
        % else
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v'});
        % end
        %
        
        model.sol('sol11').feature('s1').set('control', 'stat');
        model.sol('sol11').feature('s1').feature('p1').set('plistarr', {'0'});
        model.sol('sol11').feature('s1').feature('p1').set('pname', {'Var1D3D'});
        model.sol('sol11').feature('s1').feature('p1').set('control', 'param');
        model.sol('sol11').feature('s1').feature('p1').set('plot', 'on');
        model.sol('sol11').feature('s1').feature('p2').set('plistarr', {'range(0,0.1,5)'});
        model.sol('sol11').feature('s1').feature('p2').set('pname', {'t'});
        model.sol('sol11').feature('s1').feature('p2').set('plot', 'on');
 

    case '2D2D'
        % %2D-2D Study
        
        
        model.param.set('Var1D2D', '0');
        
        
        model.study.create('std10');
        model.study('std10').feature.create('time', 'Transient');
        
        model.sol.create('sol10');
        model.sol('sol10').study('std10');
        model.sol('sol10').attach('std10');
        model.sol('sol10').feature.create('st1', 'StudyStep');
        model.sol('sol10').feature.create('v1', 'Variables');
        model.sol('sol10').feature.create('t1', 'Time');
        model.sol('sol10').feature('t1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol10').feature('t1').feature.create('se1', 'Segregated');
        model.sol('sol10').feature('t1').feature.remove('fcDef');
        
        model.sol('sol10').attach('std10');
        model.sol('sol10').name('Solver 1D 2D');
        model.sol('sol10').feature('st1').name('Compile Equations: Time Dependent');
        model.sol('sol10').feature('st1').set('studystep', 'time');
        model.sol('sol10').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_phi2').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_c2').set('solvefor', false);
        if clamped == 1
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        
        model.sol('sol10').feature('t1').set('tlist', 'range(0,0.1,5)');
        model.sol('sol10').feature('t1').feature('se1').feature('ssDef').set('segvar', {'mod1_c' 'mod1_phi'});
        
        
        
        model.study.create('std11');
        model.study('std11').feature.create('param', 'Parametric');
        model.study('std11').feature.create('stat', 'Stationary');
        model.sol.create('sol11');
        model.sol('sol11').study('std11');
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature.create('st1', 'StudyStep');
        model.sol('sol11').feature.create('v1', 'Variables');
        model.sol('sol11').feature.create('s1', 'Stationary');
        model.sol('sol11').feature('s1').feature.create('p1', 'Parametric');
        model.sol('sol11').feature('s1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol11').feature('s1').feature.create('p2', 'Parametric');
        model.sol('sol11').feature('s1').feature.remove('fcDef');
        
        model.study('std11').feature('param').set('plistarr', {'0'});
        model.study('std11').feature('param').set('pname', {'Var1D2D'});
        model.study('std11').feature('param').set('plot', 'on');
        model.study('std11').feature('stat').set('notsolmethod', 'sol');
        model.study('std11').feature('stat').set('usesol', 'on');
        model.study('std11').feature('stat').set('activate', {'g' 'off' 'g2' 'off' 'g3' 'on' 'g4' 'off' 'g5' 'off'  ...
            'g7' 'off'});
        model.study('std11').feature('stat').set('notstudy', 'std10');
        
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature('st1').name('Compile Equations: Stationary');
        model.sol('sol11').feature('st1').set('study', 'std11');
        model.sol('sol11').feature('v1').set('notsolmethod', 'sol');
        model.sol('sol11').feature('v1').set('notsol', 'sol10');
        model.sol('sol11').feature('v1').set('control', 'stat');
        model.study('std11').feature('stat').set('notsolnum', 'all');
        model.sol('sol11').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_c').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_c2').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_phi2').set('solvefor', false);
        if clamped ~=0
            model.sol('sol11').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol11').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        % if clamped
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v','mod1_vb_lm2D'});
        % else
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v'});
        % end
        %
        
        model.sol('sol11').feature('s1').set('control', 'stat');
        model.sol('sol11').feature('s1').feature('p1').set('plistarr', {'0'});
        model.sol('sol11').feature('s1').feature('p1').set('pname', {'Var1D2D'});
        model.sol('sol11').feature('s1').feature('p1').set('control', 'param');
        model.sol('sol11').feature('s1').feature('p1').set('plot', 'on');
        model.sol('sol11').feature('s1').feature('p2').set('plistarr', {'range(0,0.1,5)'});
        model.sol('sol11').feature('s1').feature('p2').set('pname', {'t'});
        model.sol('sol11').feature('s1').feature('p2').set('plot', 'on');

    case '2D'
        %2D Study
         model.param.set('Var1D2D', '0');
       model.sol.create('sol10');
        model.study.create('std10');
        model.study('std10').feature.create('time', 'Transient');
        
        model.sol('sol10').study('std10');
        model.sol('sol10').attach('std10');
        model.sol('sol10').feature.create('st1', 'StudyStep');
        model.sol('sol10').feature.create('v1', 'Variables');
        model.sol('sol10').feature.create('t1', 'Time');
        model.sol('sol10').feature('t1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol10').feature('t1').feature.remove('fcDef');
        
        model.sol('sol10').attach('std10');
        model.sol('sol10').name('Solver 2D');
        model.sol('sol10').feature('st1').name('Compile Equations: Time Dependent');
        model.sol('sol10').feature('st1').set('study', 'std10');
        model.sol('sol10').feature('st1').set('studystep', 'time');
        model.sol('sol10').feature('v1').set('control', 'time');
        model.sol('sol10').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_c2').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod3_phi2').set('solvefor', false);
        if clamped == 1
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            %model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        
        model.sol('sol10').feature('t1').set('tlist', 'range(0,0.1,5)');
    case '1D2D'
        %1D-2D Study
        
        model.study.create('std10');
        model.study('std10').feature.create('time', 'Transient');
        
        model.sol.create('sol10');
        model.sol('sol10').study('std10');
        model.sol('sol10').attach('std10');
        model.sol('sol10').feature.create('st1', 'StudyStep');
        model.sol('sol10').feature.create('v1', 'Variables');
        model.sol('sol10').feature.create('t1', 'Time');
        model.sol('sol10').feature('t1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol10').feature('t1').feature.create('se1', 'Segregated');
        model.sol('sol10').feature('t1').feature.remove('fcDef');
        
        model.sol('sol10').attach('std10');
        model.sol('sol10').name('Solver 1D 2D');
        model.sol('sol10').feature('st1').name('Compile Equations: Time Dependent');
        model.sol('sol10').feature('st1').set('studystep', 'time');
        model.sol('sol10').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_u').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_v').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol10').feature('v1').feature('mod1_c').set('solvefor', false);
        if clamped == 1
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
        elseif clamped == 2
            model.sol('sol10').feature('v1').feature('mod1_vb_lm2D').set('solvefor', false);
            model.sol('sol10').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);          
        end
        
        model.sol('sol10').feature('t1').set('tlist', 'range(0,0.1,5)');
        model.sol('sol10').feature('t1').feature('se1').feature('ssDef').set('segvar', {'mod3_c2' 'mod3_phi2'});
        
        
        model.study.create('std11');
        model.study('std11').feature.create('param', 'Parametric');
        model.study('std11').feature.create('stat', 'Stationary');
        model.sol.create('sol11');
        model.sol('sol11').study('std11');
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature.create('st1', 'StudyStep');
        model.sol('sol11').feature.create('v1', 'Variables');
        model.sol('sol11').feature.create('s1', 'Stationary');
        model.sol('sol11').feature('s1').feature.create('p1', 'Parametric');
        model.sol('sol11').feature('s1').feature.create('fc1', 'FullyCoupled');
        model.sol('sol11').feature('s1').feature.create('p2', 'Parametric');
        model.sol('sol11').feature('s1').feature.remove('fcDef');
        
        model.study('std11').feature('param').set('plistarr', {'1'});
        model.study('std11').feature('param').set('pname', {'Var1D2D'});
        model.study('std11').feature('param').set('plot', 'on');
        model.study('std11').feature('stat').set('notsolmethod', 'sol');
        model.study('std11').feature('stat').set('usesol', 'on');
        model.study('std11').feature('stat').set('activate', {'g' 'off' 'g2' 'off' 'g3' 'on' 'g4' 'off' 'g5' 'off'  ...
            'g7' 'off'});
        model.study('std11').feature('stat').set('notstudy', 'std10');
        
        model.sol('sol11').attach('std11');
        model.sol('sol11').feature('st1').name('Compile Equations: Stationary');
        model.sol('sol11').feature('st1').set('study', 'std11');
        model.sol('sol11').feature('v1').set('notsolmethod', 'sol');
        model.sol('sol11').feature('v1').set('notsol', 'sol10');
        model.sol('sol11').feature('v1').set('control', 'stat');
        model.study('std11').feature('stat').set('notsolnum', 'all');
        model.sol('sol11').feature('v1').feature('mod2_u3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod2_v3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod2_w3').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_phi').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod1_c').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_c2').set('solvefor', false);
        model.sol('sol11').feature('v1').feature('mod3_phi2').set('solvefor', false);
        if clamped ~=0
            model.sol('sol11').feature('v1').feature('mod2_vb_lm3D').set('solvefor', false);
            model.sol('sol11').feature('v1').feature('mod2_vb_lm3Dp').set('solvefor', false);
        end
        % if clamped
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v','mod1_vb_lm2D'});
        % else
        %     model.sol('sol11').feature('s1').feature('se1').feature('ss1').set('segvar', {'mod1_u' 'mod1_v'});
        % end
        %
        
        model.sol('sol11').feature('s1').set('control', 'stat');
        model.sol('sol11').feature('s1').feature('p1').set('plistarr', {'1'});
        model.sol('sol11').feature('s1').feature('p1').set('pname', {'Var1D2D'});
        model.sol('sol11').feature('s1').feature('p1').set('control', 'param');
        model.sol('sol11').feature('s1').feature('p1').set('plot', 'on');
        model.sol('sol11').feature('s1').feature('p2').set('plistarr', {'range(0,0.1,5)'});
        model.sol('sol11').feature('s1').feature('p2').set('pname', {'t'});
        model.sol('sol11').feature('s1').feature('p2').set('plot', 'on');

end

if strcmp(study,'2D')
    tic;
    model.sol('sol10').runAll;
    solTime=toc;    
else
    tic;
    model.sol('sol10').runAll;
    model.sol('sol11').feature('s1').set('stol', '1e-3');
    model.sol('sol11').feature('s1').set('nonlin', 'on');
    model.sol('sol11').runAll;
    solTime=toc;
end    

save(strrep(outputfile,'.mph','.txt'),'solTime','-ascii')
mphsave(model, outputfile);
























