Dear Andy, Have you checked which linear solver is used by your:
gf_model_get(md, 'solve'); ? If you have compiled GetFEM with mumps, the default behavior would be to use mumps for such large systems. Can you confirm? BR Kostas On Thu, Jun 22, 2017 at 1:31 AM, Yu (Andy) Huang <[email protected]> wrote: > 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 <%28646%29%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 <%28646%29%20509-8798> > Email: [email protected] > *[email protected]* <[email protected]> > http://www.parralab.org/people/yu-andy-huang/ > <http://neuralengr.com/members/yu-%28andy%29-huang> >
