%% Multi_SolderedBar_2D_Thermal
%
% M-code simulation model "\\Home-NAS\NAS_Data\Prabhu\Simulations\FEATool\Multi_SolderedBar_2D_Thermal.m"
% Created with <https://www.featool.com FEATool Multiphysics> version 1.17.1, build 24.10.294,
% 14-Dec-2024 17:38:20 with MATLAB 9.7.0.1737446 (R2019b) Update 9 PCWIN64.
%
% Keywords: _heat_
%FEAToolMultiphysicsApp

%% Initialization
clear fea
fea.sdim = { 'x', 'y' };
fea = addphys( fea, @heattransfer, { 'T' } );

%% Geometry definition
gobj = gobj_rectangle( 0.17886, 0.66838, 0.5, 0.57186, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( [0], [5], [-0.06], [0.06], 'R1' );
fea.geom.objects{1} = gobj;
gobj = gobj_rectangle( 0, 5, 0.2, 0.30682, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( [0], [5], [0.06], [0.075], 'R2' );
fea.geom.objects{2} = gobj;
gobj = gobj_rectangle( 0, 4.7188, -0.323, -0.10339, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( [0], [5], [-0.075], [-0.06], 'R3' );
fea.geom.objects{3} = gobj;
gobj = gobj_rectangle( 0, 5, 0.2, 0.33911, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( [0], [5], [0.075], [0.3], 'R4' );
fea.geom.objects{4} = gobj;
gobj = gobj_rectangle( 0, 5, -0.32305, -0.13115, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( [0], [5], [-0.075], [-0.3], 'R5' );
fea.geom.objects{5} = gobj;
gobj = gobj_rectangle( [-0.05], [5], [0.075], [0.3], 'R4' );
fea.geom.objects{4} = gobj;
gobj = gobj_rectangle( [-0.05], [5], [-0.3], [-0.075], 'R5' );
fea.geom.objects{5} = gobj;
gobj_tags = { 'R1', 'R2', 'R3', 'R4' };
[fea.geom, gobj_tags] = copy_geometry_object( gobj_tags, fea.geom, { [ 0  0.375]} );
gobj_tags = { 'R1', 'R2', 'R3', 'R5' };
[fea.geom, gobj_tags] = copy_geometry_object( gobj_tags, fea.geom, { [ 0  -0.375]} );

%%
clf
plotgeom( fea );

%% Mesh generation
%%
fea.grid = gridgen( fea, 'gridgen', 'default', 'hmax', 0.26, 'grading', 0.3 );
%%
fea = gridrefine( fea );
%%
fea = gridrefine( fea );

%%
clf
plotgrid( fea );

% %% Equation settings
fea.phys.ht.dvar = { 'T' };
fea.phys.ht.sfun = { 'sflag1' };
fea.phys.ht.eqn.coef = { 'rho_ht', 'rho', 'Density', { 'rho_Matl1', 'rho_Matl2', 'rho_Matl2', 'rho_Matl3', 'rho_Matl3', 'rho_Matl1', 'rho_Matl2', 'rho_Matl2', 'rho_Matl3', 'rho_Matl1', 'rho_Matl2', 'rho_Matl2', 'rho_Matl3' };
                         'cp_ht', 'Cp', 'Heat capacity', { 'cp_Matl1', 'cp_Matl2', 'cp_Matl2', 'cp_Matl3', 'cp_Matl3', 'cp_Matl1', 'cp_Matl2', 'cp_Matl2', 'cp_Matl3', 'cp_Matl1', 'cp_Matl2', 'cp_Matl2', 'cp_Matl3' };
                         'k_ht', 'k', 'Thermal conductivity', { 'k_Matl1', 'k_Matl2', 'k_Matl2', 'k_Matl3', 'k_Matl3', 'k_Matl1', 'k_Matl2', 'k_Matl2', 'k_Matl3', 'k_Matl1', 'k_Matl2', 'k_Matl2', 'k_Matl3' };
                         'u_ht', 'u', 'Convection velocity in x-direction', { '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0' };
                         'v_ht', 'v', 'Convection velocity in y-direction', { '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0' };
                         'q_ht', 'Q', 'Heat source', { '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0' };
                         'T0_ht', 'T0', 'Initial condition for T', { '25', '25', '25', '25', '25', '25', '25', '25', '25', '25', '25', '25', '25' } };
fea.phys.ht.eqn.seqn = 'rho_ht*cp_ht*T'' - k_ht*(Tx_x + Ty_y) + rho_ht*cp_ht*(u_ht*Tx_t + v_ht*Ty_t) = q_ht';
fea.phys.ht.eqn.sdiff = ' - k_ht*(Tx_x + Ty_y)';
fea.phys.ht.eqn.sconv = 'rho_ht*cp_ht*(u_ht*Tx_t + v_ht*Ty_t)';
fea.phys.ht.eqn.vars = { 'Temperature, T', 'T';
                         'Temperature gradient, T', { 'Tx', 'Ty' } };
fea.phys.ht.prop.isaxi = 0;
fea.phys.ht.prop.artstab.id = 0;
fea.phys.ht.prop.artstab.id_coef = 0.5;
fea.phys.ht.prop.artstab.sd = 0;
fea.phys.ht.prop.artstab.sd_coef = 0.25;
fea.phys.ht.prop.artstab.iupw = 0;
fea.phys.ht.prop.active = [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ];
fea.phys.ht.prop.intb = 1;
%% Constants and expressions
fea.expr = { 'rho_Matl1', '5.32e-3';
             'cp_Matl1', '0.325';
             'k_Matl1', '0.055';
             'rho_Matl2', '14.52e-3';
             'cp_Matl2', '0.150';
             'k_Matl2', '0.0573';
             'rho_Matl3', '17e-3';
             'cp_Matl3', '0.16';
             'k_Matl3', '0.150' };


%% Boundary settings
fea.phys.ht.bdr.sel = [ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 3, 3, 3, 3, 3, 3, 3, 1, 3, 3, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, -1, -2, -1, -2, -1, -1, -1, -1, -1, -1, -2, -1 ];
fea.phys.ht.bdr.coef = { 'bcr_ht', 'T = T0', 'Temperature', { 'T0' }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, '25', 0, 0, 0, 0, 0, 0, 0, '25', 0, 0, 0, '25', 0, 0, 0, 0, 0, 0, 0, 0, 0, '25', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                         'bcc_ht', 'n.(-kgrad T) = 0', 'Convective flux/outflow', [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                         'bci_ht', 'n.(-kgrad T + rhoCpuT) = 0', 'Thermal insulation/symmetry', [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 'rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T' }, { '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0' };
                         'bcf_ht', '-n.(-kgrad T + rhoCpuT) = q0 + h (Tinf -T) + Const (Tamb4-T4)', 'Heat flux', { { 'q0', 'h', 'Tinf', 'Const', 'Tamb' } }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 'rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T+(', ')+(', ')*((', ')-T)+(', ')*((', ')^4-T^4)' }, { { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' }, { '0', '0', '0', '0', '0' } } };
fea.phys.ht.bdr.coefi = { 'bcic_ht', 'n.(F1-F2) = 0, Fi=-kigrad Ti  + rhoiCp,iuiTi', 'Continuity', [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                          'bcij_ht', '-n.(F1-F2) = q0, Fi=-kigrad Ti  + rhoiCp,iuiTi', 'Flux discontinuity', { 'q0' }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 'rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T' }, { '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '0', '20', '0', '20', '0', '0', '0', '0', '0', '0', '20', '0' };
                          'bcir_ht', 'T = T0', 'Temperature', { 'T0' }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } };

fea.phys.ht.bdr.vars = { 'Normal conductive heat flux, T', '-k_ht*(nx*Tx+ny*Ty)';
                         'Normal convective heat flux, T', 'rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T';
                         'Normal total heat flux, T', '-k_ht*(nx*Tx+ny*Ty)+rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T' };
fea.phys.ht.prop.intb = 1;


%% Solver
fea = parsephys( fea );
fea = parseprob( fea );

%% Modify boundary conditions for time dependant heat sources
fea.bdr.n{1,1}{36} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(20)';
fea.bdr.n{1,1}{38} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(20)';
fea.bdr.n{1,1}{45} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(20)';

[fea.sol.u,fea.sol.t] = solvetime( fea, ...
                                   'iupw', 0, ...
                                   'linsolv', 'mumps', ...
                                   'tstep', 1e-06, ...
                                   'tmax', 0.0003, ...
                                   'tstop', 1e-20, ...
                                   'ischeme', 2, ...
                                   'imass', 4, ...
                                   'icub', 'auto', ...
                                   'nlrlx', 1, ...
                                   'toldef', 1e-06, ...
                                   'tolchg', 1e-06, ...
                                   'reldef', 0, ...
                                   'relchg', 1, ...
                                   'maxnit', 20, ...
                                   'hook', [], ...
                                   'nproc', 12, ...
                                   'init', { 'T0_ht' }, ...
                                   'solcomp', [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ] );

%% Postprocessing
first.u = fea.sol.u;
first.t = fea.sol.t;
%%
clf
postplot( fea, ...
          'surfexpr', 'T', ...
          'title', 'surface: Temperature, T, Time = 0.01001', ...
          'solnum', size(fea.sol.u,2) );

%% Solver

fea.bdr.n{1,1}{36} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(0)';
fea.bdr.n{1,1}{38} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(0)';
fea.bdr.n{1,1}{45} = '(rho_ht*cp_ht*(nx*u_ht+ny*v_ht)*T)+(0)';

[fea.sol.u,fea.sol.t] = solvetime( fea, ...
                                   'iupw', 0, ...
                                   'linsolv', 'mumps', ...
                                   'tstep', 1e-05, ...
                                   'tmax', 0.01, ...
                                   'tstop', 1e-20, ...
                                   'ischeme', 2, ...
                                   'imass', 4, ...
                                   'icub', 'auto', ...
                                   'nlrlx', 1, ...
                                   'toldef', 1e-06, ...
                                   'tolchg', 1e-06, ...
                                   'reldef', 0, ...
                                   'relchg', 1, ...
                                   'maxnit', 20, ...
                                   'hook', [], ...
                                   'nproc', 12, ...
                                   'init', 'restart' , ...
                                   'solcomp', [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ] );

%% Postprocessing

%%
clf
postplot( fea, ...
          'surfexpr', 'T', ...
          'title', 'surface: Temperature, T, Time = 0.01001', ...
          'solnum', size(fea.sol.u,2) );

%% Add the two solutions
sol = horzcat(solution.u,fea.sol.u);
t_temp = fea.sol.t+solution.t(1,size(solution.t,2));
time = horzcat(solution.t,t_temp);
fea.sol.u=sol;
fea.sol.t=time;