Current Examples
This demo is implemented in a single Python file. Download here: tutorial_current.ipynb
This demo illustrates how to:
Define currents in the electrodes.
from module1_mesh import*
from module2_forward import*
from module3_inverse import*
from module4_auxiliar import*
import matplotlib.pyplot as plt
Mesh
mesh_inverse, mesh_direct=MyMesh(r=1, n=8, n_vertex=121)
plt.figure(figsize=(8, 8))
plt.subplot(1,2,1)
plot(mesh_direct);
plt.subplot(1,2,2)
plot(mesh_inverse);
![../../_images/mesh8.png](../../_images/mesh8.png)
Current Examples
Method 1
"Current"
n_g=3 #Number currents
I_all=current_method(n_g, value=1, method=1) #Creating current
#Plotting
for i in range(n_g):
mesh=mesh_direct
VD=FiniteElement('CG',mesh.ufl_cell(),1)
g_u=interpolate(I_all[i], FunctionSpace(mesh,VD))
g_u=getBoundaryVertex(mesh, g_u)
bond=plot_boundary(mesh, data=g_u, name='boundary g'+str(i))
![../../_images/method1_1.png](../../_images/method1_1.png)
![../../_images/method1_2.png](../../_images/method1_2.png)
![../../_images/method1_3.png](../../_images/method1_3.png)
>>> print("Mesh Direct:")
>>> Verifyg(I_all, mesh_direct)
>>> print("\n Mesh Inverse:")
>>> Verifyg(I_all, mesh_inverse)
Mesh Direct:
Integral boundary: 2.42861286636753e-16 0
Integral boundary: -1.3357370765021415e-16 1
Integral boundary: -3.122502256758253e-16 2
Integral boundary g(0)*g(1): 0.0
Integral boundary g(0)*g(2): 0.0
Integral boundary g(1)*g(2): 0.0
Mesh Inverse:
Integral boundary: 3.469446951953614e-18 0
Integral boundary: -9.020562075079397e-17 1
Integral boundary: -3.469446951953614e-17 2
Integral boundary g(0)*g(1): 0.0
Integral boundary g(0)*g(2): 0.0
Integral boundary g(1)*g(2): 0.0
Method 2
"Current"
n_g=3 #Number currents
I_all=current_method(n_g, value=1, method=2) #Creating current
#Plotting
for i in range(n_g):
mesh=mesh_direct
VD=FiniteElement('CG',mesh.ufl_cell(),1)
g_u=interpolate(I_all[i], FunctionSpace(mesh,VD))
g_u=getBoundaryVertex(mesh, g_u)
bond=plot_boundary(mesh, data=g_u, name='boundary g'+str(i))
>>> print("Mesh Direct:")
>>> Verifyg(I_all, mesh_direct)
>>> print("\n Mesh Inverse:")
>>> Verifyg(I_all, mesh_inverse)
Mesh Direct:
Integral boundary: 8.270294171719428e-16 0
Integral boundary: 4.163336342344337e-16 1
Integral boundary: -3.469446951953614e-17 2
Integral boundary g(0)*g(1): 1.0598076236045806e-17
Integral boundary g(0)*g(2): 0.0010576671174781075
Integral boundary g(1)*g(2): 4.217546450968612e-17
Mesh Inverse:
Integral boundary: 7.4593109467002705e-16 0
Integral boundary: -2.7582103268031233e-16 1
Integral boundary: 3.122502256758253e-17 2
Integral boundary g(0)*g(1): 3.8294020732188017e-16
Integral boundary g(0)*g(2): -7.741203511546502e-17
Integral boundary g(1)*g(2): -4.85722573273506e-17
![../../_images/method2_1.png](../../_images/method2_1.png)
![../../_images/method2_2.png](../../_images/method2_2.png)
![../../_images/method2_3.png](../../_images/method2_3.png)
Example 1
myI1=Expression(" sin(x[0]*pi) ",degree=2)
g_u=interpolate(myI1, FunctionSpace(mesh,VD))
g_u2=getBoundaryVertex(mesh, g_u)
bond=plot_boundary(mesh, data=g_u2, name='boundary g')
![../../_images/example1.png](../../_images/example1.png)
>>> Verifyg([g_u], mesh_direct)
Integral boundary: -1.4991805540043313e-16 0
Example 2
myI2=Expression(" x[1]>0 ? 1 :-1 ",degree=1)
g_u=interpolate(myI2, FunctionSpace(mesh,VD))
g_u2=getBoundaryVertex(mesh, g_u)
bond=plot_boundary(mesh, data=g_u2, name='boundary g'+str(i))
![../../_images/example2_1.png](../../_images/example2_1.png)
>>> print(assemble(g_u*ds(mesh))) #Integral boundary #Like Verifyg
-1.6306400674181987e-15
Example 3
value=2
n_g=2
myI3=[Expression(f" x[1]>=0 ? {value}*sin(acos(x[0])*{i+1}) : {value}*sin((-acos(x[0]))*{i+1})",degree=1) for i in range(0,n_g)]
for i in range(n_g):
mesh=mesh_direct
VD=FiniteElement('CG',mesh.ufl_cell(),1)
g_u=interpolate(myI3[i], FunctionSpace(mesh,VD))
g_u2=getBoundaryVertex(mesh, g_u)
bond=plot_boundary(mesh, data=g_u2, name='boundary g'+str(i))
![../../_images/example3_1.png](../../_images/example3_1.png)
![../../_images/example3_2.png](../../_images/example3_2.png)