%% FEATool Multiphysics Version 1.11, Build 19.07.210, Model M-Script file.
%% Created on 13-Feb-2020 13:43:45 with MATLAB 9.5.0.944444 (R2018b) PCWIN64.


%% Starting new model.
fea.sdim = { 'x', 'y', 'z' };
fea.geom = struct;
fea = addphys( fea, @conductivemediadc, { 'V' } );

%% Geometry operations.
gobj = gobj_block( [0], [20], [0], [20], [0], [20], 'B1' );
fea = geom_add_gobj( fea, gobj );
gobj = gobj_ellipsoid( [0 0 0], [1], [1], [1], [1], [16], 'E1' );
fea = geom_add_gobj( fea, gobj );
fea = geom_apply_transformation( fea, 'E1', [10 10 10 ], [ 1  1  1 ], 0, 1 );

%% Grid operations.
fea.grid = gridgen( fea, 'gridgen', 'gmsh', 'hmax', 0.3, 'nsm', 25 );

%% Equation settings.
fea.phys.dc.dvar = { 'V' };
fea.phys.dc.sfun = { 'sflag1' };
fea.phys.dc.eqn.coef = { 'dts_dc', 'd_ts', 'Time scaling coefficient', { '0', '0' };
                         's_dc', 'sigma', 'Electric conductivity', { '.02', '.12' };
                         'q_dc', 'Q', 'Current source', { '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '0' } };
fea.phys.dc.eqn.seqn = 'dts_dc*V'' - s_dc*(Vx_x + Vy_y + Vz_z) = q_dc';
fea.phys.dc.eqn.vars = { 'Electric potential, V', 'V';
                         'Electric field, -grad(V)', { '-Vx', '-Vy', '-Vz' } };
fea.phys.dc.prop.isaxi = 0;
fea.phys.dc.prop.active = [ 1, 1 ];
fea.phys.dc.prop.intb = 0;
fea.phys.dc.dvar = { 'V' };
fea.phys.dc.sfun = { 'sflag1' };
fea.phys.dc.eqn.coef = { 'dts_dc', 'd_ts', 'Time scaling coefficient', { '0', '0' };
                         's_dc', 'sigma', 'Electric conductivity', { '.02', '.12' };
                         'q_dc', 'Q', 'Current source', { '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '0' } };
fea.phys.dc.eqn.seqn = 'dts_dc*V'' - s_dc*(Vx_x + Vy_y + Vz_z) = q_dc';
fea.phys.dc.eqn.vars = { 'Electric potential, V', 'V';
                         'Electric field, -grad(V)', { '-Vx', '-Vy', '-Vz' } };
fea.phys.dc.prop.isaxi = 0;
fea.phys.dc.prop.active = [ 1, 1 ];
fea.phys.dc.prop.intb = 0;
fea.phys.dc.dvar = { 'V' };
fea.phys.dc.sfun = { 'sflag1' };
fea.phys.dc.eqn.coef = { 'dts_dc', 'd_ts', 'Time scaling coefficient', { '0', '0' };
                         's_dc', 'sigma', 'Electric conductivity', { '.02', '.12' };
                         'q_dc', 'Q', 'Current source', { '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '0' } };
fea.phys.dc.eqn.seqn = 'dts_dc*V'' - s_dc*(Vx_x + Vy_y + Vz_z) = q_dc';
fea.phys.dc.eqn.vars = { 'Electric potential, V', 'V';
                         'Electric field, -grad(V)', { '-Vx', '-Vy', '-Vz' } };
fea.phys.dc.prop.isaxi = 0;
fea.phys.dc.prop.active = [ 1, 1 ];
fea.phys.dc.prop.intb = 0;

%% Boundary settings.
fea.phys.dc.bdr.sel = [ 2, 2, 2, 2, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1 ];
fea.phys.dc.bdr.coef = { 'bcp_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, 0, 0, '100', 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                         'bcf_dc', '-n.(-sigmagrad V) = J_0', 'Current flow', { 'J_0' }, { 0, 0, 0, 0, 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.dc.bdr.coefi = { 'bcic_dc', 'n.(F_1-F_2) = 0, F_i=-sigma_igrad V_i', '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 };
                          'bcij_dc', '-n.(F_1-F_2) = J_0, F_i=-sigma_igrad V_i', 'Flux discontinuity', { 'J_0' }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                          'bcir_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 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 } };
fea.phys.dc.bdr.vars = { 'Normal current density', '-s_dc*(nx*Vx+ny*Vy+nz*Vz)' };
fea.phys.dc.prop.intb = 0;
fea.phys.dc.bdr.sel = [ 2, 2, 1, 2, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1 ];
fea.phys.dc.bdr.coef = { 'bcp_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, '0', 0, '100', 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                         'bcf_dc', '-n.(-sigmagrad V) = J_0', 'Current flow', { 'J_0' }, { 0, 0, 0, 0, 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.dc.bdr.coefi = { 'bcic_dc', 'n.(F_1-F_2) = 0, F_i=-sigma_igrad V_i', '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 };
                          'bcij_dc', '-n.(F_1-F_2) = J_0, F_i=-sigma_igrad V_i', 'Flux discontinuity', { 'J_0' }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                          'bcir_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 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 } };
fea.phys.dc.bdr.vars = { 'Normal current density', '-s_dc*(nx*Vx+ny*Vy+nz*Vz)' };
fea.phys.dc.prop.intb = 0;
fea.phys.dc.bdr.sel = [ 2, 2, 1, 2, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1 ];
fea.phys.dc.bdr.coef = { 'bcp_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, [], { 0, 0, '0', 0, '100', 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                         'bcf_dc', '-n.(-sigmagrad V) = J_0', 'Current flow', { 'J_0' }, { 0, 0, 0, 0, 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.dc.bdr.coefi = { 'bcic_dc', 'n.(F_1-F_2) = 0, F_i=-sigma_igrad V_i', '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 };
                          'bcij_dc', '-n.(F_1-F_2) = J_0, F_i=-sigma_igrad V_i', 'Flux discontinuity', { 'J_0' }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, [], { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
                          'bcir_dc', 'V = V_0', 'Electric potential', { 'V_0' }, { 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 } };
fea.phys.dc.bdr.vars = { 'Normal current density', '-s_dc*(nx*Vx+ny*Vy+nz*Vz)' };
fea.phys.dc.prop.intb = 0;

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

fea.sol.u = solvestat( fea, ...
                       'iupw', 0, ...
                       'linsolv', 'mumps', ...
                       'icub', 'auto', ...
                       'nlrlx', 1, ...
                       'toldef', 1e-06, ...
                       'tolchg', 1e-06, ...
                       'reldef', 0, ...
                       'relchg', 1, ...
                       'maxnit', 20, ...
                       'init', { 'V0_dc' }, ...
                       'solcomp', [ 1, 1 ] );

%% Postprocessing.
postplot( fea, ...
          'surfexpr', 'V', ...
          'alpha', 1, ...
          'colormap', 'jet', ...
          'colorbar', 1, ...
          'solnum', 1 );

%% Solver call.

