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 <[email protected]> Date: Tue, Jun 20, 2017 at 3:10 AM Subject: Re: [Getfem-users] Simulating electric field distribution To: "Yu (Andy) Huang" <[email protected]>, [email protected] 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]> 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 <%28646%29%20509-8798> > Email: [email protected] > *[email protected]* <[email protected]> > 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 <(646)%20509-8798> Email: [email protected] *[email protected]* <[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 --------- -- 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] *[email protected]* <[email protected]> http://www.parralab.org/people/yu-andy-huang/ <http://neuralengr.com/members/yu-%28andy%29-huang>
