Thanks Yves and Konstantinos! The iterative solver was used in my case. So I tried to re-compile getFEM with mumps and parallel computing enabled, and installed again (before this I installed METIS and Mumps), but for some unknown reason it made Matlab crash once the gf_model_get(md, 'solve') is executed.
So I realized I must have messed up everything. Then I deleted getFEM and downloaded the latest version of getFEM (v5.2) and re-installed it. Now it always crashes Matlab when tries to solve the model (a segmentation violation error was thrown out by Matlab). I did a "make check" and found the check_all.sh failed. Do you know what might go wrong with my installation? Below is the commands I used for installation: ./configure --enable-matlab --enable-mumps --with-pic make cp interface/src/matlab/gf_matlab.mexa64 interface/src/matlab/gf_matlab sudo make install sudo cp /usr/local/getfem_toolbox/gf_matlab /usr/local/getfem_toolbox/gf_matlab.mexa64 Thanks again for your help! On Thu, Jun 22, 2017 at 8:05 AM, Konstantinos Poulios < [email protected]> wrote: > 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> >> > > -- 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/ <http://neuralengr.com/members/yu-%28andy%29-huang>
