# 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__()
model3D.synchronize()
geom.generate_mesh(dim=3)
gmsh.write("cone.msh")
model3D.__exit__()

import meshio

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 = 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)
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.