### Fenics loses geometry data

338
views
0
13 months ago by
Hi all, I've been having a weird problem that I can't seem to fix. I've been using gmsh to create geometry and meshes, then running them through Fenics to solve poisson's equation in 3D, and outputting the results as a .pvd to view in ParaView. The solver runs fine and I get understandable results. However, one of the geometries (the bottom rectangular prism labeled "Substrate", id number 2) I defined in gmsh doesn't make it through to the .pvd. It just doesn't display, even though the boundary condition applied to it works fine. I've attached pictures, and the fenics code. Any ideas why this might be? I'm thinking it might have something to with the fact that I specified a couple surfaces as members of multiple Physical Groups in gmsh, but I'm not sure why that would cause problems...let me know if any other info is needed!

Fenics code:
from fenics import *
from mshr import *
%matplotlib notebook

mesh = Mesh("/home/fenics/shared/electrodes_fem/gmsh/electrodes_3d_4.xml")
domains = MeshFunction('size_t', mesh, "/home/fenics/shared/electrodes_fem/gmsh/electrodes_3d_4_physical_region.xml")
boundaries = MeshFunction('size_t', mesh, "/home/fenics/shared/electrodes_fem/gmsh/electrodes_3d_4_facet_region.xml")

V = FunctionSpace(mesh, 'Lagrange', 2)

bc_border = DirichletBC(V, Constant(0), boundaries, 1)
bc_e1 = DirichletBC(V, Constant(1000), boundaries, 3)
bc_e2 = DirichletBC(V, Constant(-1000), boundaries, 4)
bc_substrate = DirichletBC(V, Constant(0), boundaries, 2)

bcs = [bc_border, bc_e1, bc_e2, bc_substrate]

dx = Measure('dx', domain=mesh, subdomain_data=domains)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Save solution in VTK format
file = File("../shared/electrodes_fem_3d/electrode_field_fem_3d_4.pvd")
file << u

# Plot solution and mesh
plot(u)​

Gmsh screenshot:

ParaView screenshot:
Community: FEniCS Project