Question about selecting boundaries for boundary conditions

dolfinx
,
1

I’m quite new to fenics, this being the first project aside from playing with the tutorials / examples provided here: The FEniCSx tutorial — FEniCSx tutorial

I’m running it in a docker container using “dolfinx/lab:stable” and using the complex number compatible kernel DOLFINx.

I’m trying to simulate the Takagi-Taupin equations on a simple 2D rectangle. As in step 1 of the procedure described in this paper X-ray dynamical diffraction from single crystals with arbitrary shape and strain field: A universal approach to modeling.

The code is running without errors but the results make me suspect I’ve not implemented the boundary condition the way I expected. This is further supported by the last figure where I was trying to plot where in the mesh it thinks the boundary Gamma is not showing the top and left side of the rectangle but a bottom right corner.
If someone could tell me where I made a mistake (in the BC or elsewhere) or just tel me if I’m properly plotting the selected boundary that would be great.

import numpy as np
from mpi4py import MPI #message passage interface (How to comunicate between processes mostly relevant for parralelisiation)
import dolfinx
import ufl

from dolfinx.fem import petsc
from petsc4py import PETSc

sample_width = 300e-6 #in m
sample_height = 100e-6 #in m

res_x = 12*5
res_y = 4*5

Amplitude = 1. #Amplitude of incomming wave
lam = 25e-6 #Wavelength in m

Angle_in = 17.7 #Incident angle w.r.t. surface of incomming plane wave in degrees
Angle_in_rad = np.deg2rad(Angle_in)
s0 = np.array([np.cos(Angle_in_rad),-np.sin(Angle_in_rad),0.]) # Direction of incomming wave vector
s0 = s0 / np.linalg.norm(s0) # Normalisation
s0 = 2*np.pi*s0 / lam
s_0 = ufl.as_vector(s0[:2])
chi0 = 11.7
chih = 11.
chibh = 11.

#Generating the mesh
domain = dolfinx.mesh.create_rectangle(comm = MPI.COMM_WORLD, # message passage interface
                                       points = [[0.,0.],[sample_width,sample_height]], #corners of rectangle lowerleft - upper right
                                       n = (res_x,res_y), #horizontal and vertical resolution
                                       dtype = np.float64 #data type for mesh geometry (how accuratly the positions can be saved
                                      )

V = dolfinx.fem.functionspace(mesh = domain,
                              element = ("Lagrange",1) #use standard scale lagrange elements of continuous polynomials of order 2
                             )

#Functions defining different boundary segments
def boundary_Gamma(x):
    bound = np.logical_or(np.isclose(x[0,:], 0.), np.isclose(x[1,:], sample_height),dtype=bool)
    return bound

def boundary_notGamma(x):
    bound = np.invert(boundary_Gamma(x))
    return bound

# creating different integration measures for different boundary segments
boundaries = [(1, boundary_Gamma),
              (2, boundary_notGamma)]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
    facets = dolfinx.mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.classes.Measure("ds", domain=domain, subdomain_data=facet_tag) #create custom measure to integrate over different boundary regions seperatly

D0 = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def Da0(x):
    '''Represents the incomming wave, returns amplitude of that wave at position'''
    # print(type(x))
    # global Y
    # Y = x
    return Amplitude*np.exp(1j*np.dot(s0,x))

chi_0 = dolfinx.fem.Constant(domain,chi0)
chi_h = dolfinx.fem.Constant(domain,chih)
chi_bh = dolfinx.fem.Constant(domain,chibh)
Da_0 = dolfinx.fem.Function(V,
                            name = 'Da0',
                            dtype=complex)
Da_0.interpolate(Da0)
Dh_start = dolfinx.fem.Constant(domain, PETSc.ScalarType(0 - 0j))
const1 = dolfinx.fem.Constant(domain, PETSc.ScalarType(-1j * (np.pi / lam) * chi0))
const2 = dolfinx.fem.Constant(domain, PETSc.ScalarType(1j * (np.pi/lam) * chibh * Dh_start))


# # later to use when D_h not identically zero
# D_h = dolfinx.fem.Function(V,
#                            name = 'Dh',
#                            dtype=complex)
# D_h.x.array[:] = 0.0

#Define dirichlet BC where on the segment defined by boundary_Gamma D0 should equal the incomming wave
bc = dolfinx.fem.dirichletbc(Da_0, dolfinx.fem.locate_dofs_geometrical(V, boundary_Gamma))

#Defining lhs as a and rhs as b
a = (- ufl.inner(s_0*D0, ufl.grad(v)) +  ufl.inner(const1 *D0, v)) * ufl.dx + ufl.inner(D0,v) * ds(2) # 
l = ufl.inner(const2,v) * ufl.dx - ufl.inner(Da_0, v) * ds(1)

problem = petsc.LinearProblem(a,l,bcs=[bc])

D0_solve = problem.solve()

import pyvista
from dolfinx import plot

pyvista.start_xvfb()
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim, tdim)
topology, cell_types, geometry = plot.vtk_mesh(domain,tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("fundamentals_mesh.png")

u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

# u_grid.point_data["D_h"] = D_h.x.array.imag
# u_grid.set_active_scalars("D_h")

# u_grid.point_data["Da_0"] = Da_0.x.array.real
# u_grid.set_active_scalars("Da_0")

# u_grid.point_data["Da_0"] = Da_0.x.array.imag
# u_grid.set_active_scalars("Da_0")

# u_grid.point_data["D0"] = D0_solve.x.array.real
# u_grid.set_active_scalars("D0")

u_grid.point_data["D0"] = D0_solve.x.array.real
u_grid.set_active_scalars("D0")

# u_grid.point_data["D0"] = D0_solve.x.array.real**2 + D0_solve.x.array.imag**2
# u_grid.set_active_scalars("D0")

u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

# cells, types, x = plot.vtk_mesh(V, entities=facet_tag.find(1))

subplotter = pyvista.Plotter(shape=(1, 2))
subplotter.subplot(0, 0)
subplotter.add_text("Mesh with markers", font_size=14, color="black", position="upper_edge")
subplotter.add_mesh(grid, show_edges=True, show_scalar_bar=False)
subplotter.view_xy()

dof_gamma = dolfinx.mesh.locate_entities_boundary(domain,fdim,boundary_Gamma)
# dof_gamma = dolfinx.fem.locate_dofs_geometrical(V, boundary_Gamma)
cells, types, x = plot.vtk_mesh(V, entities =  dof_gamma)

sub_grid = pyvista.UnstructuredGrid(cells, types, x)
subplotter.subplot(0, 1)
subplotter.add_text("Subset of mesh", font_size=14, color="black", position="upper_edge")
subplotter.add_mesh(sub_grid, show_edges=True, edge_color="black")
subplotter.view_xy()

if pyvista.OFF_SCREEN:
    figsize = 100
    subplotter.screenshot(
        "2D_markers.png", transparent_background=transparent, window_size=[2 * figsize, figsize]
    )
else:
    subplotter.show()

p.s. this is the first time ever I post on one of these forums so If I didn’t properly paste the code please tell me how to do it.

2

I’ve discovered that there isn’t a mistake in the code it is working properly. The mistake was in my understanding of what the results of the simulation should have been.

I still don’t understand why I can’t just show the boundary the way I try using pyvista. As shown here

2D_markers
2D_markers2000×1000 87 KB

So if someone could explain that I would be grateful. :slight_smile:

3

What you are trying to do is to only plot the degrees of freedom associated with a set of facets. This is not possible with the call that you have tried above.
You could try:

facets_gamma = dolfinx.mesh.locate_entities_boundary(domain,fdim,boundary_Gamma)
cells_gamma = dolfinx.mesh.compute_incident_entities(msh.topology, facets_gamma, fdim, fdim + 1)
cells, types, x = plot.vtk_mesh(V, entities =  cells_gamma)

This finds all cells associated with the facets, and plots the solution on these

For further questions, I would strongly recommend you to make a smaller example. As this is just about post-processing, you don’t really need to solve a PDE for it.

4

This doesn’t change the resulting figure. It still only plots a mesh with the cells in the bottom right corner not the cells along the left and top side.

*I’ve made the edit to your code snippet replacing the msh.topology with domain.topology as I assume that is what you intended (my code has no object named msh)

5

Please note that the amount of code you have presented makes it hard to be very concrete. I’ve narrowed it down to there being a bug in vtk_mesh when the input is a function space, and not a mesh.
Smaller example illustrating this:

import numpy as np
from mpi4py import MPI
import dolfinx

sample_width = 300e-6  # in m
sample_height = 100e-6  # in m

res_x = 12 * 5
res_y = 4 * 5

# Generating the mesh
domain = dolfinx.mesh.create_rectangle(
    comm=MPI.COMM_WORLD,  # message passage interface
    points=[
        [0.0, 0.0],
        [sample_width, sample_height],
    ],  # corners of rectangle lowerleft - upper right
    n=(res_x, res_y),  # horizontal and vertical resolution
    dtype=np.float64,  # data type for mesh geometry (how accuratly the positions can be saved
)

V = dolfinx.fem.functionspace(
    mesh=domain,
    element=(
        "Lagrange",
        1,
    ),
)


def boundary_Gamma(x):
    bound = np.logical_or(
        np.isclose(x[0, :], 0.0), np.isclose(x[1, :], sample_height), dtype=bool
    )
    return bound


def boundary_notGamma(x):
    bound = np.invert(boundary_Gamma(x))
    return bound


# creating different integration measures for different boundary segments
boundaries = [(1, boundary_Gamma), (2, boundary_notGamma)]

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for marker, locator in boundaries:
    facets = dolfinx.mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(
    domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
)

import pyvista
from dolfinx import plot

if pyvista.OFF_SCREEN:
    pyvista.start_xvfb(0.1)
tdim = domain.topology.dim
domain.topology.create_connectivity(tdim, tdim)
topology, cell_types, geometry = plot.vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("fundamentals_mesh.png")

u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)


u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

# cells, types, x = plot.vtk_mesh(V, entities=facet_tag.find(1))

subplotter = pyvista.Plotter(shape=(1, 2))
subplotter.subplot(0, 0)
subplotter.add_text(
    "Mesh with markers", font_size=14, color="black", position="upper_edge"
)
subplotter.add_mesh(grid, show_edges=True, show_scalar_bar=False)
subplotter.view_xy()


facets_gamma = dolfinx.mesh.locate_entities_boundary(domain, fdim, boundary_Gamma)
cells_gamma = dolfinx.mesh.compute_incident_entities(
    domain.topology, facets_gamma, fdim, fdim + 1
)
cells, types, x = plot.vtk_mesh(V, entities=cells_gamma)


sub_grid = pyvista.UnstructuredGrid(cells, types, x)
subplotter.subplot(0, 1)
subplotter.add_text(
    "Subset of mesh", font_size=14, color="black", position="upper_edge"
)
subplotter.add_mesh(sub_grid, show_edges=True, edge_color="black")
subplotter.view_xy()

if pyvista.OFF_SCREEN:
    figsize = 100
    subplotter.screenshot(
        "2D_markers.png",
        transparent_background=False,
        window_size=[2 * figsize, figsize],
    )
else:
    subplotter.show()

where if one change
cells, types, x = plot.vtk_mesh(V, entities=cells_gamma)
to
cells, types, x = plot.vtk_mesh(domain, entities=cells_gamma)
one gets

image
image921×731 39.3 KB

6

@dokken That is it, thank you. I’ll make sure to make more minimal code examples in future.

7

I’ve made an issue at: [BUG]: vtk_mesh on subset of cells based of function space · Issue #3607 · FEniCS/dolfinx · GitHub which I might have time to fix in the next few days.