Some questions about meshes and boundaries for fenicsx

Hi all, first time poster here, although I have used the legacy version of fenics in the past.

I would like to ask two things about migrating my calculations to fenicsx.

Is there an equivalent way to create simple 3d meshes in fenicsx? For example (in legacy fenics) sphere:

sphere = Sphere(Point(1, 2, 3), 5)
mesh = generate_mesh(sphere, 32)

or cone

cone = Cone(Point(1, 2, 3), Point(4, 5, 6), 7)
mesh = generate_mesh(cone, 32)

Also with the command

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
points = domain.geometry.x

we can save the coordinates of the whole mesh. How can I get an array of only the points that are in the boundary? Something like

boundary_mesh = BoundaryMesh(mesh, 'exterior')
boundary_nodes = boundary_mesh.coordinates()  

for the legacy fenics.

Apologies if the questions are naive.

The mesh generator mshr that was used in legacy Dolfin has been long deprecated. To create simple meshes I would suggest using Gmsh or pygmsh.

Both of these modules uses the gmsh.model which can be directly converted into a DOLFINx mesh using
dolfinx.io.gmshio.model_to_mesh;

You could combine

facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_vertices=dolfinx.cpp.mesh.entities_to_geometry(mesh, mesh.topology.dim-1, facets)
coordinates=mesh.geometry.x[boundary_vertices.reshape(-1)]
1 Like

Thanks for your answer. I tried what you posted. Seems that everything works but I am not sure about something, that is probably trivial. Consider, for example, a simple cone:

import pygmsh
import gmsh
import numpy

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
    return out_mesh


mesh_size_1 = 0.1
geom = pygmsh.occ.Geometry()
model3D = geom.__enter__()
cone0 =  model3D.add_cone([0,0,0],[0,0,5],0,1,mesh_size=mesh_size_1)
model3D.synchronize()
geom.generate_mesh(dim=3)
gmsh.write("cone.msh")
model3D.__exit__()

import meshio
mesh_from_file = meshio.read("cone.msh")

triangle_mesh = create_mesh(mesh_from_file, "triangle")
meshio.write("cone.xdmf", triangle_mesh)

Unless I am doing something wrong the cone seems… hollow in paraview for lack of a better word. Shouldn’t the cone had points in the interior?

Im not sure what you expect.
You are only saving triangles (which meshio only picks those marked as entity groups in gmsh).
If you want to mesh the 3D mesh, I would expect you to save tetrahedral cells, i.e.

tetra_mesh = create_mesh(mesh_from_file, "tetra")
meshio.write("cone.xdmf", tetra_mesh)
2 Likes

Thanks, works fine now.

I have a hopefully final question, thankfully not very naive. In my legacy fenics code I use for a variety of reasons external functions for my boundaries that have numerical calculations on them. I did something like this

class BCfunction(UserExpression):
    def eval(self, value, x):
        value[0] = f(x)  #  some other extrernal function of x with non trivial numerics
    def value_shape(self):
        return () 


u_D = BCfunction()

def boundary(x, on_boundary):
  return on_boundary

bc = DirichletBC(V, u_D, boundary)

However I am not sure what is the better way to do this in fenicsx

You would be more explicit in DOLFINx.
You should:

  1. create a function u_bc = dolfinx.fem.Function(V).
  2. Create a python function
def bc_func(x):
    return f(x) 

f(x) is here assumed to be vectorized, i.e. it takes in a 2D array of points
x = ((x_0, x_1, ... x_N), (y_0, y_1, ..., y_N), (z_0, z_1, ...., z_N)
If your function is not vectorized, you would need to transpose the input

def bc_func(x):
    xT = x.T
    values = np.zeros(x.shape[1])
    for (i,coord) in enumerate(xT)):
         values[i] = f(coord)
    return values
  1. Interpolate: u_bc.interpolate(bc_func)
  2. Use u_bc in the dirichletbc.
2 Likes

Thank you very much, it worked perfectly.

As a last note, I guess if my function f(x) calculates the boundary conditions only at the actual boundaries of my geometry (that is why I asked how to get them in my original question), i.e. an analytical function in position space x but with parameters from other calculations, is the interpolation step optional?

In DOLFINx, you can find all cells incident to the boundary, as explained in: Interpolating results on an edge - #2 by dokken
and only interpolate onto those cells.

Remember to add u_bc.x.scatter_forward() after an interpolation on a subset of cells.