Thanks for you previous help on my silly questions! Now I manage to simulate the electric field on a toy sphere with two layers, each layer having a different electrical conductivity. I'm just not sure if I did it properly, because when I compare the results to those I got from Abaqus (a commercial FEM solver), I see some difference that I don't understand (see attached screenshots). I used the same mesh, with the same boundary condition and same conductivity values, but the distribution and absolute values of voltage and field are all different between getFEM and Abaqus. My major concern is the way I coded the boundary conditions. The problem I'm solving is a Laplacian equation of electric potential u, with the following BC:
I only coded explicitly (1) and (2) so not sure if it's good enough. I put part of my code below (Matlab code). Any advice on the code is much appreciated!
% ===================================================
mesh = gfMesh('import','gmsh', 'toySphere.msh');
mfu = gf_mesh_fem(mesh, 1); % scalar-field (electric potential)
mfE = gf_mesh_fem(mesh, 3); % 3d vector-field (electric field)
gf_mesh_fem_set(mfu, 'fem', gf_fem('FEM_PK(3,1)')); % P1 Lagrange
gf_mesh_fem_set(mfE, 'fem', gf_fem('FEM_PK(3,1)')); % P1 Lagrange
mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(1)')); % integration method
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mfu);
rid = gf_mesh_get(mesh,'regions');
reg1 = gf_mesh_get(mesh, 'region', rid(1));
reg2 = gf_mesh_get(mesh, 'region', rid(2));
region1 = 1;
region2 = 2;
gf_mesh_set(mesh, 'region', region1, reg1);
gf_mesh_set(mesh, 'region', region2, reg2);
% governing equation and conductivities
sigma = [0.276;0.126]; % conductivity values for the two layers in the sphere
gf_model_set(md, 'add linear generic assembly brick', mim, [num2str(sigma(1)) '*(Grad_u.Grad_Test_u)'],region1);
gf_model_set(md, 'add linear generic assembly brick', mim, [num2str(sigma(2)) '*(Grad_u.Grad_Test_u)'],region2);
fb1 = gf_mesh_get(mesh, 'outer faces with direction', [0 0 1], 0.1, cvid(indAnode));
fb2 = gf_mesh_get(mesh, 'outer faces with direction', [0 0 -1], 0.1, cvid(indCathode));
% the boundary condition is injecting 1 A/m^2 current density at the north pole of the sphere, with the south pole as ground
% here indAnode and indCathode is the index of the element corresponding to the north and south pole
anode_area = 3;
cathode_area = 4;
gf_mesh_set(mesh, 'region', anode_area, fb1);
gf_mesh_set(mesh, 'region', cathode_area, fb2);
gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mfu, cathode_area);
gf_model_set(md, 'add initialized data','Jn', ones(gf_mesh_get(mesh, 'nbpts'),1));
% here is the line that I suspect mostly. I have to pass 'Jn' a vector, otherwise it won't solve.
gf_model_set(md, 'add source term brick', mim, 'u', [num2str(sigma(1)) '*(-Grad_u.Normal)'], anode_area, 'Jn');
% Neumann BC of electric potential
% solve
gf_model_get(md, 'solve');
% extracted solution
u = gf_model_get(md, 'variable', 'u');
E = gf_model_get(md, 'interpolation', '-Grad_u', mfu); % electric field
% display
figure; gf_plot(mfu, u, 'mesh','on','cvlst', get(mesh, 'outer faces')); colormap(jet); colorbar
figure; gf_plot(mfE, E, 'mesh','on','norm','on','cvlst', get(mesh, 'outer faces')); colormap(jet); colorbar
compareAbaqusGetfem.pdf