Some questions about meshes and boundaries for fenicsx

1

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.

2

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)]
2 Likes
3

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?

4

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
5

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

6

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
7

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?

8

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.

1 Like