function [ fea, out ] = helmholtz2DFEA( varargin )
% M-code simulation model "C:\Users\chern\OneDrive - CYCLOMED TECHNOLOGIES\Documentos - Cyclomed\General\Proyectos\EST01\3 Estudio EM\1 MATLAB\helmholtz2DFEA.m"
% Created with <https://www.featool.com FEATool Multiphysics> version 1.17.3, build 25.02.45,
% 26-May-2025 13:43:50 with MATLAB 9.14.0.2206163 (R2023a) PCWIN64.
%
% Keywords: _custom, pde_

%% Initialization
clear fea
cOptDef = {...
            'Je', 1e8;
            'air', 2;
            'superRin',0.5;
            'superRout', 0.55;
            'superHdown', 0.2;
            'superHup', 0.4;
            'ironPole', 0.4;
            'ironIntHeight', 1;
            'ironExtHeight', 1.5;
            'ironIntFace', 0.8;
            'ironExtFace', 1.4;
            'fid', 1};

% cOptDef = {...
%             'Je', varargin{1};
%             'air', varargin{2};
%             'superRin',varargin{3};
%             'superRout', varargin{4};
%             'superHdown', varargin{5};
%             'superHup', varargin{6};
%             'ironPole', varargin{7};
%             'ironIntHeight', varargin{8};
%             'ironExtHeight', varargin{9};
%             'ironIntFace', varargin{10};
%             'ironExtFace', varargin{11};
%             'fid', varargin{12}};
[got,opt] = parseopt(cOptDef,varargin{:});

fea.sdim = { 'r', 'z' };
fea = addphys( fea, {@customeqn, 1}, { 'A' } );

%% Geometry definition
gobj = gobj_rectangle( 0, opt.air, 0, opt.air, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( opt.superRin, opt.superRout, opt.superHdown, opt.superHup, 'R2');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( 0, opt.ironPole, opt.superHdown, opt.ironIntHeight, 'R3');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( 0, opt.ironExtFace, opt.ironIntHeight, opt.ironExtHeight, 'R4');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( opt.ironIntFace, opt.ironExtFace, 0, opt.ironIntHeight, 'R5');
fea = geom_add_gobj( fea, gobj );
fea.geom = geom_apply_formula( fea.geom, 'R3+R4+R5' );

%% Mesh generation

% grid_temp = gridgen(fea, 'hmax', 0.5);  % Coarse temporary grid
% n_bdr = max(grid_temp.b(3,:));

% Define different mesh sizes for subdomains
hmax = [0.1, 0.1, 0.4, 0.02];  

% % Define boundary mesh sizes (zeros for default, specific values for targeted boundaries)
% hmaxb = zeros(1, n_bdr);
% hmaxb([2]) = 0.1;  % Finer mesh on boundaries 1 and 3

% Generate grid with specified mesh sizes
fea.grid = gridgen(fea, 'hmax', hmax, 'gridgen','gridgen2d','grading', 0.3);
fea = gridrefine( fea );

% plotgrid(fea)

%% Equation settings
fea.phys.ce.dvar = { 'A' };
fea.phys.ce.sfun = { 'sflag1' };
fea.phys.ce.eqn.seqn = '1/(mu0*(mu_rel*mu_rel_iron(r,A,Ar,Az)+1))*(Ar_r+Az_z-A_t/r^2+Ar_t/r) = J';
fea.phys.ce.eqn.coef = { 'A0_ce', 'A0', 'Initial condition for A', { '0', 0, 0 } };
fea.phys.ce.eqn.vars = { 'Dependent variable, A', 'A';
                         'Gradient of A', { 'Ar', 'Az' } };
fea.phys.ce.prop.isaxi = 1;
fea.phys.ce.prop.active = [ 1, 1, 1, 1 ];
fea.phys.ce.prop.intb = 0;

% plotsubd(fea)
%% Constants and expressions

% % Option 0: Constant permeability
fea.expr = { 'Je', opt.Je;
             'ironIntFace', opt.ironIntFace;
             'ironIntHeight', opt.ironIntHeight ;
             'mu_iron', 1;
             'mu0', 'pi*4e-7';
             'Br', '-Az';
             'Bz', '(r>0)*(A/r + Ar)+(r<=0)*2*Ar';
             'B', 'sqrt(Br^2+Bz^2)';
             'J', '0 0 0 Je';
             'mu_rel', '0 1 0 0';
             'Bzz', '2*Ar' };

% Option 1: External function
% fea.expr = { 'Je', opt.Je;
%              'mu_iron', 1;
%              'mu0', 'pi*4e-7';
%              'Br', '-Az';
%              'Bz', '(r>0)*(A/r + Ar)+(r<=0)*2*Ar';
%              'B', 'sqrt(Br^2+Bz^2)';
%              'J', '0 0 Je';
%              'mu_rel', '1 mu_rel_iron(r,A,Ar,Az) 1';
%              'Bzz', '2*Ar' };

% Option 2: Analytical function
% fea.expr = { 'Je', opt.Je;
%              'mu_iron', 1;
%              'mu0', 'pi*4e-7';
%              'Br', '-Az';
%              'Bz', '(r>0)*(A/r + Ar)+(r<=0)*2*Ar';
%              'B', 'sqrt(Br^2+Bz^2)';
%              'J', '0 0 Je';
%              'mu_rel', '1 2890*exp(-((((A./r + Ar).^2+(-Az).^2).^(1/2)-0.7)/0.66)^2)+1 1';
%              'Bzz', '2*Ar' };

%% Boundary settings
fea.phys.ce.bdr.sel = [ 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 ];
fea.phys.ce.bdr.coef = { 'bcnd_ce', 'Dirichlet/Neumann boundary conditions', 'Dirichlet/Neumann boundary conditions', { 'Dirichlet, r_u', 'Neumann, g_u' }, { 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, { { 'r*(', ')' };
                         [] }, { '0', '0', '0', 0, 0, '0', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 } };
fea.phys.ce.bdr.vars = [];
fea.phys.ce.bdr.coefi = { 'bcndi_ce', 'Dirichlet/Neumann boundary conditions', 'Dirichlet/Neumann boundary conditions', { 'Dirichlet, r_u', 'Neumann, g_u' }, { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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.ce.prop.intb = 0;

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

fea.sol.u = solvestat( fea, ...
                       'iupw', 0, ...
                       'linsolv', 'bigcstab', ...
                       'icub', 'auto', ...
                       'nlrlx', 0.1, ...
                       'toldef', 1e-06, ...
                       'tolchg', 1e-06, ...
                       'reldef', 0, ...
                       'relchg', 1, ...
                       'maxnit', 100, ...
                       'hook', [], ...
                       'nproc', 1, ...
                       'init', { 'A0_ce' }, ...
                       'solcomp', [ 1, 1, 1, 1 ] );

% fea = fenics( fea, 'fid', fid, ...
%               'nlrlx', 0.5, ...           % Reduced relaxation for stability
%               'maxnit', 50, ...           % Increased iterations
%               'toldef', 1e-6, ...         % Solver tolerance
%               'solcomp', [1, 1, 1] );

%% Postprocessing
% Plot magnetic field
clf
postplot( fea, ...
          'surfexpr', 'B', ...
          'title', 'surface: B', ...
          'solnum', 1 );

% Point evaluation "Bzz"
p_eval = [ 0; 0];
out.Bcenter = evalexpr( 'Bzz', p_eval, fea, 1 );

% Maximum magnetic field SC
[min_val, out.Bsuper, min_coord, max_coord] = minmaxsubd('B', fea, 4);
[min_val, out.Biron1] = minmaxsubd('(r>ironIntFace)*(r<2)*B', fea);
[min_val, out.Biron2] = minmaxsubd('(z>ironIntHeight)*(z<2)*B', fea);

out.Biron = max(out.Biron1, out.Biron2);
out.Biron = out.Biron1;
end