getfem-users
[Top][All Lists]
Advanced
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Fwd: Simulating electric field distribution

From: Yu (Andy) Huang
Subject: Re: [Getfem-users] Fwd: Simulating electric field distribution
Date: Fri, 7 Jul 2017 17:51:36 -0400
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 <address@hidden> 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 <address@hidden> 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 <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


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 <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!

--
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,



--
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,


-- 

  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/~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,




--
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,
reply via email to
[Prev in Thread] Current Thread [Next in Thread]