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

[Yves Renard]

Dear Andy,

Unfortunately, the parallel version using mpi is not compatible with Matlab 
(may be the openmp one is, but I never tried). You should install the 
sequential version of Mumps (without Metis) or use the python interface which 
is compatible with the parallel version of Mumps (but in that case, mpirun  as 
to be used to run the python script, see the example in the documentation).

Best regards,

Yves.

----- Original Message -----
From: "Yu (Andy) Huang" <address@hidden>
To: address@hidden
Sent: Friday, July 7, 2017 11:51:36 PM
Subject: Re: [Getfem-users] Fwd: Simulating electric field distribution

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,
>>> <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: address@hidden
>>>       address@hidden <address@hidden>
>>> 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: address@hidden
>>       address@hidden <address@hidden>
>> http://www.parralab.org/people/yu-andy-huang/
>>
>>
>> --
>>
>>   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,
>> <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: address@hidden
>>       address@hidden <address@hidden>
>> 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: address@hidden
      address@hidden <address@hidden>
http://www.parralab.org/people/yu-andy-huang/
<http://neuralengr.com/members/yu-%28andy%29-huang>