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
cOptDef = {...
            'Je', 1e8;
            'fid', 1};
[got,opt] = parseopt(cOptDef,varargin{:});
fid = opt.fid;

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

%% Geometry definition
gobj = gobj_rectangle( 0, 2, 0, 2, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( 0.5, 0.55, 0.2, 0.4, 'R1');
fea = geom_add_gobj( fea, gobj );
gobj = gobj_rectangle( 0, 0.4, 0.2, 1, 'R1');
fea = geom_add_gobj( fea, gobj );

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

%% Equation settings
fea.phys.ce.dvar = { 'A' };
fea.phys.ce.sfun = { 'sflag1' };
fea.phys.ce.eqn.seqn = '(1/(mu0*mu_rel))*(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 ];
fea.phys.ce.prop.intb = 0;

%% Constants and expressions

% % Option 0: Constant permeability
% 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 1 1';
%              '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_2a', '2890*exp(-((sqrt((A/r+Ar)^2+Az^2)-0.7)/0.66)^2)+1';
             'mu_rel', '1 mu_rel_2(r,z) 1';
             'Bzz', '2*Ar' };

%% Boundary settings
fea.phys.ce.bdr.sel = [ 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_A', 'Neumann, g_A' }, { 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, { { 'r*(', ')' };
                         [] }, { '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_A', 'Neumann, g_A' }, { 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;

fea = parsephys( fea );
fea = parseprob( fea );

fea.sol.u = evalinit( fea );

% fea.vars = struct('name', 'mu_rel_2', ...
%                   'type', 'var', ...
%                   'sfun', 'sflag1', ...
%                   'data', 2891 * ones(size(fea.grid.p,2),1));

%% Solver
it = 0;
global fea0
fea0 = fea;
% nlrlx = 0.7;
while true
  clear mu_rel_2
  it = it + 1

  fea.sol.u = solvestat( fea );
  % fea.sol.u = nlrlx * fea.sol.u + (1 - nlrlx) * fea0.sol.u;

  ix = fea.sol.u ~= 0;
  chg_u = norm((fea.sol.u(ix) - fea0.sol.u(ix)) ./ fea.sol.u(ix))

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

  clf
  title(sprintf('it = %d', it))
  subplot(2,2,1)
  postplot(fea,'surfexpr','A','title','A')
  subplot(2,2,2)
  postplot(fea,'surfexpr','mu_rel','title','mu_rel','evalstyle','exact')
  subplot(2,2,3)
  postplot(fea,'surfexpr','Ar','title','Ar','evalstyle','exact')
  subplot(2,2,4)
  postplot(fea,'surfexpr','Az','title','Az','evalstyle','exact')
  drawnow

  if( chg_u < 1e-3 )
    break
  end

  % r = fea.grid.p(1,:)'; r(r <= 0) = sqrt(eps);
  % z = fea.grid.p(2,:)';
  % p = [r,z]';
  % A = evalexpr('A', p, fea);
  % Ar = evalexpr('Ar', p, fea);
  % Az = evalexpr('Az', p, fea);
  % mu_rel_2 = 2890*exp(-((sqrt((A./r+Ar).^2+Az.^2)-0.7)./0.66).^2)+1;
  % [min(mu_rel_2), mean(mu_rel_2), max(mu_rel_2)]
  % fea.vars.data = mu_rel_2;

  fea0 = fea;
end
