|Subject:||Re: [Getfem-users] Fwd: Simulating electric field distribution|
|Date:||Thu, 22 Jun 2017 14:05:59 +0200|
Thanks for your generous help, Yves!! Now it worked, and it also worked on the real model!It took more than 4 hours to solve the real model (with about 1 million tetrahedra elements). Do you know any tricks to speed up the solving? I'm thinking about parallelize the solving and was reading the documentation. It says I need to compile the package again with parallel option enabled? Also do you have any empirical values for the options in calling gf_model_get(model,'solve')? e.g. the "max_iter" and "max_res"?Thanks again for your help!---------- Forwarded message ----------
From: Yves Renard <address@hidden>
Date: Tue, Jun 20, 2017 at 3:10 AM
Subject: Re: [Getfem-users] Simulating electric field distribution
To: "Yu (Andy) Huang" <address@hidden>, address@hidden
I just add some commentaries on your code.
Le 19/06/2017 à 21:27, Yu (Andy) Huang a écrit :
The integration method is used for both the volumic term and the boundary conditions, so it has to be of order twoDear 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 = 12) ground at the south pole: V = 03) insulation at all outer boundary: n.J = 04) 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
-> mim = gf_mesh_im(mesh, gf_integ('IM_TETRAHEDRON(2)'))
Did you check that your regions were ok (for instance with gf_plot_mesh ) ?
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 conductivitiessigma = [0.276;0.126]; % conductivity values for the two layers in the spheregf_model_set(md, 'add linear generic assembly brick', mim, [num2str(sigma(1)) '*(Grad_u.Grad_Test_u)'],regio
n1);gf_model_set(md, 'add linear generic assembly brick', mim, [num2str(sigma(2)) '*(Grad_u.Grad_Test_u)'],regio n2);
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 ) ?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 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');
gf_model_set(md, 'add source term brick', mim, 'u', '1', anode_area);
% Neumann BC of electric potential
% solvegf_model_get(md, 'solve');
% extracted solutionu = gf_model_get(md, 'variable', 'u');E = gf_model_get(md, 'interpolation', '-Grad_u', mfu); % electric field
% displayfigure; gf_plot(mfu, u, 'mesh','on','cvlst', get(mesh, 'outer faces')); colormap(jet); colorbarfigure; 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 <address@hidden> 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!
-- Yves Renard (address@hidden) 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/~ren
|[Prev in Thread]||Current Thread||[Next in Thread]|