When reading the codes in Fenicsx tutorial about heat equation, I have a question about the following codes:
boundary_facets=mesh.locate_entities_boundary(
domain,fdim,lambda x:np.full(x.shape[1],True,dtype=bool))
What does ‘x.shape[1]’ mean? When implementing the function ‘locate_entities_boundary’ , where does ‘x’ come from?Is it from ‘tabulate_dof_coordinates()’?Could you mind explaining what the ‘lambda’ function does in ‘locate_entities_boundary’ ?Thanks in advance!
The third input argument to locate_entities_boundary is a python function, that takes in coordinates
x, of shape (3, num_points), such that x[0] are all points in the x-coordinate, x[1] all points in y-coordinate etc.
locate_entities_boundary works on the computational domain (domain) in your case. It will find all entities of dimension fdim (facets in your case), and check if the vertices satisfy this criterion.
A lambda function is just a compact way of packing
def check_coordinates(x):
return np.full(x.shape[1], True, dtype=bool))
boundary_facets=mesh.locate_entities_boundary(
domain,fdim,check_coordinates)
In this case, since you want all facets on the boundary, we say that as long as the facet is on the boundary, it’s vertex satisfy the criterion (i.e. we return an array full of True, for all points).
You could also use
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
Thank you very much.I still have some questions :‘locate_entities_boundary’ will only look for facets which are already on boundary ,or it will look for all facets including ones inside domain?And ‘x’ is stored in (3,num_points) form, where num_points for all mesh points or dof points?Who supplies ‘x’ for the funtion locate_entities_boundary? Thanks in advance!
There is another function for looking at all facets (including interior ones, namely dolfinx.mesh.locate_entities.
locate_entities_boundary uses exterior_facet_indices to get the facets. Then for each of the facets; it finds the vertices, and checks if all vertices satisfies the condition in the Python function.
Sorry for another question: this np.full() function return an narray full of true but locate_entities_boundary need a marker function which returns a True or False? Do I misunderstand anything?
In general it needs a function that returns an array of booleans, whose shape is the same as the second axis of the input.
See for instance: Hyperelasticity — FEniCSx tutorial
on how to mark parts of a boundary using numpy.isclose
Oh,I just look for the document which says that in locate_entities_boundary, the marker function returns an array of boolean value of size num_points, which is True if it is on the boundary.So np.full() returns what it wants.
Hello dokken, I have another question about it.Since x[0] and x[1] are just different in x or y coordinates but with the same number of points, I replace x.shape[1] with x.shape[0] but this time the code throws an error : Length of array of markers is wrong.So what is the reason here?
The shape of x is (3, number of points). Thus x.shape[0] will return 3. not len(x[0]),.
len(x[0])=len(x[1])=number of points.
I see.Thank you very much.
Dear dokken,
much later I’m coming back to this. I’m also a newbie and am trying to identify the boundary entities of an imported .msh - file. I browsed some of the older forum posts, where I often stumbled into the exterior_facet_indices() - method (and that it only works on certain dolfinx versions), as you also suggested here.
This is my code snippet to read in the mesh, which failed:
comm = MPI.COMM_WORLD
dolfinx_mesh, _, _ = read_from_msh('my_mesh.msh', comm, gdim=2)
# _ _ are the cell and face tags, which are currently unused
tdim = dolfinx_mesh.topology.dim
dolfinx_mesh.topology.create_connectivity(tdim - 1, tdim)
exterior_facets_ids = dolfinx_mesh.exterior_facet_indices(dolfinx_mesh.topology)
I get the usual TraceBack though on dolfinx version 0.7.3:
AttributeError: 'Mesh' object has no attribute 'exterior_facet_indices'
So instead I would like to use the locate_entities_boundary() - method potentially. The issue is I don’t know how to define the marker function, which I need to pass into it…
How could I either:
Sincerely
Joshua
This should be
facets_ids = dolfinx.mesh.exterior_facet_indices(dolfinx_mesh.topology)
from the dolfinx.mesh module.
Thank you, that was truly an easy fix…
As a follow up - I’m now trying to plot the boundaries of my mesh differently from the rest of the domain with pyvista.
In the docs under Visualisation with PyVista there’s is the example Mesh tags and using subplots, in which cell tags are being created with the meshtags function.
From the Visualization Demo:
cell_tags = meshtags(msh, msh.topology.dim, np.arange(num_cells), in_circle(midpoints))
# Create VTK mesh
cells, types, x = plot.vtk_mesh(msh)
grid = pyvista.UnstructuredGrid(cells, types, x)
# Attach the cells tag data to the pyvista grid
grid.cell_data["Marker"] = cell_tags.values
grid.set_active_scalars("Marker")
Is there a similar way/ function to add tags (boolean values, as I understand) to the boundary facets, of which I now have the indices? So that I can then plot them seperately in pyvista?
Or am I on the wrong path. If so, what would be the easiest way to plot the boundaries in pyvista after loading the mesh in with read_from_mesh()?
Consider the following:
# Plot a meshtag object of facets using pyvista
# Author: Jørgen S. Dokken
from mpi4py import MPI
import dolfinx
import numpy as np
def left(x):
return np.isclose(x[0], 0.0)
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
# Give each exterior facet a unique value
values = np.arange(len(facets), dtype=np.int32)
# Set all values on left boundary facet to -1
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
for i, facet in enumerate(facets):
if facet in left_facets:
values[i] = -1
mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)
mesh_data = dolfinx.plot.vtk_mesh(mesh, mt.dim, mt.indices)
# Plot meshtags
import pyvista
grid = pyvista.UnstructuredGrid(*mesh_data)
# Attach the cells tag data to the pyvita grid
grid.cell_data["Marker"] = mt.values
grid.set_active_scalars("Marker")
# Create a plotter with two subplots, and add mesh tag plot to the
# first sub-window
pyvista.start_xvfb()
plotter = pyvista.Plotter(shape=(1, 1), off_screen=True)
plotter.subplot(0, 0)
plotter.add_text(
"Mesh with markers", font_size=14, color="black", position="upper_edge"
)
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
plotter.screenshot("meshtag.png")
Here I use off_screen=True and screenshot as I wanted to generate the following png. You can of course change these
