How to set the complex coefficient of the differential equation?


164
views
1
5 months ago by
I am modeling electric stimulation of biological tissues. I need to solve the following complex equation:
\begin{equation}
- \nabla ( [ \sigma + j \omega \varepsilon ] \nabla V) = 0
\end{equation}
I solve the real part of the equation:
\begin{equation}
- \nabla ( [ \sigma ] \nabla V) = 0
\end{equation}
using the following code:
from fenics import *

# Mesh
mesh = UnitCubeMesh(5, 5, 5)
V = FunctionSpace(mesh, 'P', 1)

# Boundary conditions
u_Left = Constant(1)
def boundary_Left(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], 0, tol)
bc_L = DirichletBC(V, u_Left, boundary_Left)

u_Right = Constant(0)
def boundary_Right(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], 1, tol)
bc_R = DirichletBC(V, u_Right, boundary_Right)

# Problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
coefficient = Constant(1) # Coefficient, which should be complex
a = coefficient * dot(grad(u), grad(v)) * dx
L = f * v * dx
u = Function(V)

# Solve
solve(a == L, u, [bc_L, bc_R])​
As I understand it, I can not just make the coefficient complex because in that case I get ComplexWarning during the simulation.
How can I specify a complex coefficient for a differential equation?
Community: FEniCS Project
I have tried to get an answer from the post https://www.allanswered.com/post/vnebp/#mboab.
Solution code:
from fenics import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(3, 3)

Vr = FiniteElement('P', mesh.ufl_cell(), 1)
Vi = FiniteElement('P', mesh.ufl_cell(), 1)
Vcomplex = MixedElement([Vr, Vi])
V = FunctionSpace(mesh, Vcomplex)

# Defining boundary conditions
u_L = Constant((1, 0))
def boundary_L(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], 0, tol)
bc_L = DirichletBC(V, u_L, boundary_L)

u_R = Constant((0, 0))
def boundary_R(x, on_boundary):
    tol = 1e-14
    return on_boundary and near(x[0], 1, tol)
bc_R = DirichletBC(V, u_R, boundary_R)

# Real and Imaginary parts of the trial and test functions
ur, ui = TrialFunctions(V)
vr, vi = TestFunctions(V)

# Phisics problem
conductivity = Constant(1.0)
permittivity = Constant(1.0)
f = Constant(0)
a = - conductivity * inner(grad(ur), grad(vr))*dx - permittivity * inner(grad(ui), grad(vi))*dx
L = f*(vr+vi)*dx

# Solution
u = Function(V)
solve(a == L, u, [bc_L, bc_R])
ur, ui = u.split()

# Plotting
plot(ur)
plt.show()
plot(ui)
plt.show()​

In the code, I define a complex FunctionSpace and solve equation in it. The code works, but it seems that in this case imaginary part always equals to zero.
In theory, the should be bouth real and imaginary parts.

It seems that mistake is in definition of a and L. How to make variables a and L correctly?

written 3 months ago by Mikhail  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »