3D Dirichlet boundary conditions on subdomains

81
views
0
3 months ago by
Hey,
anybody have 3D equivalent snippet to this?
# Define boundaries
boundary_markers = MeshFunction('size_t', mesh, 2, mesh.domains())
boundaries = MeshFunction('size_t', mesh, 1, mesh.domains())

# Use the cell domains associated with each facet to set the boundary
facets_domains = [(f, [boundary_markers[c] for c in cells(f)]) for f in facets(mesh)]
for f, d in facets_domains:
boundaries[f] = 0
boundaries[f] = 1 if 1 in d else boundaries[f]
boundaries[f] = 2 if 2 in d else boundaries[f]

# Boundary on central shape
bc1 = DirichletBC(V, Constant(0.0), boundaries, 1)
bc2 = DirichletBC(V, Constant(100.0), boundaries, 2)

# Outer Boundary
bc_domain = DirichletBC(V, Constant(0.0), 'on_boundary')​

I'm trying to switch from 2D (works marvelously) where I used mshr (3D mshr doesn't support subdomains).

After fighting with dolphin-convert that consistently fails on
File "/usr/lib/python2.7/dist-packages/dolfin_utils/meshconvert/meshconvert.py", line 288, in gmsh2xml
elem_type = int(element[1])
ValueError: invalid literal for int() with base 10: ''​

I've switched to pygmsh + meshio and I'm able to import model and domains with
mesh = Mesh("shared/mesh/mesh.xml")
boundary_markers = MeshFunction("size_t", mesh, "shared/mesh/mesh_gmsh:physical.xml")
boundaries = MeshFunction("size_t", mesh, "shared/mesh/mesh_gmsh:geometrical.xml")​

Also, I might be doing things wrong (and I expect there should be easier way to do this), if so - any advice is welcome.

Thanks
Community: FEniCS Project