%% FEATool Multiphysics Version 1.13, Build 20.11.306, Model M-Script file.
%% Created on 06-Nov-2020 19:44:58 with MATLAB 9.8.0.1359463 (R2020a) Update 1 PCWIN64.


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

%% Geometry operations.
gobj = gobj_sphere( [0 0 0], [1], [1], 'S1' );
fea = geom_add_gobj( fea, gobj );
gobj = gobj_ellipsoid( [2 0 0], [1], [2], [3], [1 0 0], 'E1' );
fea = geom_add_gobj( fea, gobj );
gobj = gobj_ellipsoid( [2 0 0], [1], [2], [.5], [0 0 1], 'E1' );
fea.geom.objects{2} = gobj;
gobj = gobj_block( [-2], [4], [-3], [3], [-2], [2], 'B1' );
fea = geom_add_gobj( fea, gobj );

%% Grid operations.
fea.grid = gridgen( fea, 'gridgen', 'gmsh', 'hmax', 0.48, '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', 1, '1' };
                         's_dc', 'sigma', 'Electric conductivity', { '.12', 1, '1' };
                         'q_dc', 'Q', 'Current source', { '0', 0, '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', 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, 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', '1' };
                         's_dc', 'sigma', 'Electric conductivity', { '.12', '.12', '1' };
                         'q_dc', 'Q', 'Current source', { '0', '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '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, 1 ];
fea.phys.dc.prop.intb = 0;

%% Boundary settings.

%% 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', '0' };
                         's_dc', 'sigma', 'Electric conductivity', { '.12', '.12', '.02' };
                         'q_dc', 'Q', 'Current source', { '0', '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '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, 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', '0' };
                         's_dc', 'sigma', 'Electric conductivity', { '.12', '.12', '.02' };
                         'q_dc', 'Q', 'Current source', { '0', '0', '0' };
                         'V0_dc', 'V_0', 'Initial condition for V', { '0', '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, 1 ];
fea.phys.dc.prop.intb = 0;

%% Boundary settings.
fea.phys.dc.bdr.sel = [ -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 2, 2, 2, 2, 1, 2 ];
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, 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, '10', 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, 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, 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, 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, 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 } };
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, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 2, 2, 1, 2, 1, 2 ];
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, 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, '10', 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, 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, 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, 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, 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 } };
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, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 2, 2, 1, 2, 1, 2 ];
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, 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, '10', 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, 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, 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, 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, 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 } };
fea.phys.dc.bdr.vars = { 'Normal current density', '-s_dc*(nx*Vx+ny*Vy+nz*Vz)' };
fea.phys.dc.prop.intb = 0;

