Scaling Laplacian from 2D to 3D

mesh
1

Hello all!

I am attempting to model the Laplacian of electric potential across a solid electrolyte, with an equation that is just the Laplacian of phi = 0. I have implemented this in two dimensions (4x2 rectangle), but I am struggling to create a 3D mesh that will accept the domain of a 4x2x2 box. How can I create this 3D mesh? Also, when I would be ‘solving’ the equation, would it be solving for just the faces of the box or actually within the box itself (which is what I want)?

Here is my current code implementation of the 2D Laplacian.

from fenics import *
from mshr import *
import numpy as np

# Create mesh and define function space
domain = Rectangle(Point(0, 0), Point(4, 2))
mesh = generate_mesh(domain, 64) # The higher the value the more resolved the mesh is

V = FunctionSpace(mesh, 'P', 1) # 1st order FEM

# Define boundaries 
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_markers.set_all(0)

# Neumann conditions are already 0 if I don't include them, they should be 0 for my problem.
Voltage = 1.0  # Applied voltage

# Dirichlet conditions
phi_anode = Constant(0.0)  # Potential at the anode, which can be set to zero
phi_cathode = Constant(Voltage)  # Potential at the cathode, where Voltage is the applied voltage

# Apply Dirichlet conditions at the left and right boundaries of the domain
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
left_boundary = CompiledSubDomain('near(x[0], 0)')
right_boundary = CompiledSubDomain('near(x[0], 4)')

left_boundary.mark(boundaries, 1)
right_boundary.mark(boundaries, 2)

bc_anode = DirichletBC(V, phi_anode, boundaries, 1)
bc_cathode = DirichletBC(V, phi_cathode, boundaries, 2)

bcs = [bc_anode, bc_cathode]  # List of Dirichlet boundary conditions

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = Constant(0.0)*v*dx  # Since the right-hand side of Poisson's equation is zero

Thanks for the help!

George

2

Welcome to the FEniCS discourse forum. As a general suggestion, when you post code please make sure that it contains all required imports. Without those, it is harder for people to run your code and, in this specific case, realize that you are generating the mesh through mshr.

You can use

domain = Box(Point(0, 0, 0), Point(4, 2, 2))

to get the 3D domain. The PDE will be solved in the interior of the box.

If you are new to FEniCS, you may want to read The new DOLFINx solver is now recommended over DOLFIN.

3

My apologies, I just fixed those imports. However, I understand how to create the 3D box domain. But, will this Box() create a volume mesh which I can create a function space to solve the PDE in?

4

Yes, it will. You can export the mesh to file (either VTK or XDMF format) and open it with paraview if you want to get a closer look and convince yourself that it actually is a mesh which extends to the interior of the box.

5

Thank you for that clarification. Sorry to be a pain, but how would I save it as a PVD to view? I can’t seem to save it.

6

File("mesh.pvd") << domain

If you are new to fenics i would strongly suggest using DOLFINx

Described in

with demos at
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos.html
and The FEniCSx tutorial — FEniCSx tutorial