Dear Andy,
I just add some commentaries on your code.
Yves.
Le 19/06/2017 à 21:27, Yu (Andy) Huang a écrit :
> Dear getFEM users,
>
> 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:
>
> 1) injecting 1 A/m^2 current density at the north pole of the sphere:
> -*n*.*J *= 1
> 2) ground at the south pole: V = 0
> 3) insulation at all outer boundary: *n*.*J* = 0
> 4) continuity for inner boundary: *n*.(*J1* - *J2*) = 0
>
> 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
The integration method is used for both the volumic term and the
boundary conditions, so it has to be of order two
-> mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(2)'));
>
> 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);
Did you check that your regions were ok (for instance with gf_plot_mesh ) ?
>
> % 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);
Did you check that your boundaries were ok (may be this is your first
graph ) ?
>
> 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');
The source term isthe right hand side of -J.n = f, so in your case just
1. So you do not need any "gf_model_set(md, 'add initialized data','Jn',
ones(gf_mesh_get(mesh, 'nbpts'),1));" and you can simply add
gf_model_set(md, 'add source term brick', mim, 'u', '1', anode_area);
> % 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
> %==============================================================================
>
>
>
>
> On Thu, Jun 8, 2017 at 10:55 PM, Yu (Andy) Huang
> <[email protected] <mailto:[email protected]>> wrote:
>
> Dear getFEM users,
>
> I'm entirely new to getFEM, and I'm trying to simulate the
> electric field distribution in the human brain when direct
> electric current is applied on the scalp surface. I know it's just
> a Laplacian equation of the electric potential, and I managed to
> simulate the voltage distribution on a toy (a cube).
>
> Now my question is: how do I simulate the electric field? should I
> add another variable of electric field? or can I just get the
> field from the voltage solution? I tried both but without any
> luck. I added the electric field as a new variable but did not
> figure out how to properly add boundary condition using
> gf_model_set(). If calculating field from voltage, I didn't find
> out which function to use to establish a relation between the
> field variable and voltage variable.
>
> Any suggestion is appreciated! The examples in the documentation
> are generally mechanical problems, and there are very limited
> online resources, so I really get stuck here.
>
> Thanks a lot!
>
> --
> Yu (Andy) Huang, Ph.D.
> Postdoc fellow at Dept. of Biomedical Engineering, City College of
> New York
> Center for Discovery and Innovation, Rm. 3.320,
>
> <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/>
> 85 St Nicholas Terrace, New York, NY 10027
>
> <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/>
> Tel: 1-646-509-8798 <tel:%28646%29%20509-8798>
> Email: [email protected] <mailto:[email protected]>
> [email protected]_
> <mailto:[email protected]>__
> http://www.parralab.org/people/yu-andy-huang/
> <http://www.parralab.org/people/yu-andy-huang/>
>
>
>
>
> --
> Yu (Andy) Huang, Ph.D.
> Postdoc fellow at Dept. of Biomedical Engineering, City College of New
> York
> Center for Discovery and Innovation, Rm. 3.320,
> <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/>
> 85 St Nicholas Terrace, New York, NY 10027
> <http://neuralengr.com/directions-to-cdi-center-for-discovery-and-innovation/>
> Tel: 1-646-509-8798
> Email: [email protected] <mailto:[email protected]>
> [email protected]_ <mailto:[email protected]>__
> http://www.parralab.org/people/yu-andy-huang/
--
Yves Renard ([email protected]) tel : (33) 04.72.43.87.08
Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29
20, rue Albert Einstein
69621 Villeurbanne Cedex, FRANCE
http://math.univ-lyon1.fr/~renard
---------
