Hi everyone!
I want to extract the coordinates of the boundary points and the value of a function on them in FEniCSx.
I need the values on all boundary.
Does anyone know how I can do it?
Thank you so much.
Hi everyone!
I want to extract the coordinates of the boundary points and the value of a function on them in FEniCSx.
I need the values on all boundary.
Does anyone know how I can do it?
Thank you so much.
Are you here talking about the mesh vertices, or the degrees of freedom?
What function space are you considering?
v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
f = Function(V)
f is a function null in the domain but with values on the boundary. My goal is to deform the boundary of the mesh according to this function.
Consider the following:
from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 16, 16)
v_cg1 = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(mesh, v_cg1)
f = dolfinx.fem.Function(V)
f.interpolate(lambda x: (0.1*x[1]*abs(x[0]), 0.05*np.sin(2*x[0])))
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, boundary_facets, mesh.topology.dim-1, 0)
vertex_to_geometry = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, 0, boundary_vertices, False)
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
mesh.topology.create_connectivity(0, mesh.topology.dim)
v_to_c = mesh.topology.connectivity(0, mesh.topology.dim)
num_dofs = V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts
dof_to_geometry_map = np.full(num_dofs, -1, dtype=np.int32)
dofmap = V.dofmap
layout = dofmap.dof_layout
for i, (vertex, node) in enumerate(zip(boundary_vertices, vertex_to_geometry)):
cell = v_to_c.links(vertex)[0]
cell_dofs = dofmap.cell_dofs(cell)
cvs = c_to_v.links(cell)
local_index = np.flatnonzero(cvs == vertex)[0]
dof = cell_dofs[layout.entity_dofs(0, local_index)]
dof_to_geometry_map[dof] = node
# Create boundary perturbation vector
bs = V.dofmap.bs
perturbation_data = np.zeros((len(boundary_vertices), 3), dtype=np.float64)
geom = np.zeros(len(boundary_vertices), dtype=np.int32)
c = 0
for i, node in enumerate(dof_to_geometry_map):
if node != -1:
perturbation_data[c, :bs] = f.x.array[bs*i:bs*(i+1)]
geom[c] = node
c+=1
mesh.geometry.x[geom]+= perturbation_data
with dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
I tried it, but I get the next error:
TypeError Traceback (most recent call last)
Cell In[19], line 4
2 mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
3 boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
----> 4 boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, boundary_facets, mesh.topology.dim-1, 0)
5 vertex_to_geometry = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, 0, boundary_vertices, False)
8 c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
TypeError: compute_incident_entities(): incompatible function arguments. The following argument types are supported:
1. (mesh: dolfinx.cpp.mesh.Mesh, entities: numpy.ndarray[numpy.int32], d0: int, d1: int) -> numpy.ndarray[numpy.int32]
Invoked with: <dolfinx.cpp.mesh.Topology object at 0x7f81e628dc60>, array([ 1, 3, 8, ..., 1384673, 1384676, 1384678],
dtype=int32), 1, 0
if I replace mesh.topology by mesh, I get the next error:
AttributeError Traceback (most recent call last)
Cell In[20], line 5
3 boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
4 boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh, boundary_facets, mesh.topology.dim-1, 0)
----> 5 vertex_to_geometry = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, 0, boundary_vertices, False)
8 c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
9 mesh.topology.create_connectivity(0, mesh.topology.dim)
AttributeError: 'Mesh' object has no attribute '_cpp_object'
What version of dolfinx are you running?
I’m ruuning 0.5.2 version.
See: dolfinx/test_refinement.py at v0.5.2 · FEniCS/dolfinx · GitHub for usage (Im not at my computer, so I cannot run and adapt my code above).
Thank you @dokken . I can modify the code with this!
Hello probably I was having a similar situation. Can I have the coordinates of the boundary points and the function value immediately from your code, say coordinate_list = [...]
and vector_value_list = [...]
such that f(evaluated at coordinate_list[i]) = vector_value_list[i]
? If not, maybe I will post another question for this situation.