Hello everyone,
I’m a beginner in the use of Fenics Project and I want to solve an electrostatic problem in a dielectric.
Presently, I consider a two-dimensional system with \Omega = [0.1] \times [0.1].
The two functions I want to determine are :
- The bound charge density \rho(\vec{r}), a scalar field.
- The polarization \vec{m}(\vec{r}), a vector field.
These two functions obey the following system of partial differential equations:
The boundary conditions are as follows:
- The system is periodic according to y: \rho(x,y)=\rho(x,1-y) and \vec{m}(x,y) = \vec{m}(x,1-y).
- Dirichlet condition on the left edge: \rho(0, y) = \rho_0 and \vec{m}(0,y) = m_0 \vec{u}_x (m_0 and \rho_0 will be dependent on y then).
- Dirichlet condition on the right edge: \rho(0, y) = 0 and \vec{m}(0,y) = \vec{0}.
Fishing here and there for information from the tutorial and the forum, I was able to write the following piece of code:
from fenics import *
# Define parameters of PDEs
K = 0.014285714285714285
K_l = -0.0020138584850228334
beta = 1.1188102694571296e-06
# Create mesh
mesh = UnitSquareMesh(32, 32)
# Define function spaces and mixed (product) space
# From https://fenicsproject.org/olddocs/dolfin/latest/python/demos/mixed-poisson/demo_mixed-poisson.py.html
P_m = FiniteElement('RT', triangle, 1)
P_rho = FiniteElement('DP', triangle, 0)
V = FunctionSpace(mesh, P_m * P_rho)
# Define trial and test functions
m, rho = TrialFunctions(V)
m_test, rho_test = TestFunctions(V)
# Define bilinear form
## (K+1) m part
F = (K+1)*rho*rho_test*dx
## K_l * Laplacian part
F+= -K_l*(dot(grad(rho), grad(rho_test))*dx)
## beta * Bilaplacian part
## From : https://fenicsproject.org/olddocs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html
## Not so clear why the variationnal formulation of the bilaplacan is so complicated ...
### Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
### Penalty parameter
alpha = Constant(8.0)
F+= beta*(
inner(div(grad(rho)), div(grad(rho_test)))*dx \
- inner(avg(div(grad(rho))), jump(grad(rho_test), n))*dS \
- inner(jump(grad(rho), n), avg(div(grad(rho_test))))*dS \
+ alpha/h_avg*inner(jump(grad(rho),n), jump(grad(rho_test),n))*dS)
## Second PDE part
F+= (dot(m, m_test) + dot(grad(rho), m_test))*dx
# Boundary conditions
# bcs = ...
# Solve
w = Function(V)
solve(lhs(F) == rhs(F), w)
# Plot result
m, rho = w.split()
plot(m)
plot(rho)
But as you can see it is still missing the part on boundaries conditions (periodicity + Dirichlet). I searched on the internet but I couldn’t understand how to implement them with a mix space (vectorial * scalar).
If someone could help me, I would be really grateful.
G.
