How can I output the coordinates of a specific boundary and its corresponding computation results as two lists in dolfinx?

dolfinx
1

I have read a lot of related content, but I still cannot achieve the desired functionality in dolfinx：https://fenicsproject.discourse.group/t/how-to-extract-meshs-coordinates-and-values-according-to-facet-markers/11420/2?u=french_fries
I can only compute all the coordinates：

coordinates = mesh.geometry.x

Here is my current code：

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx.io import gmshio
from ufl import (Measure, TestFunction, TrialFunction, dot, dx, grad)
from dolfinx.fem import (Function, functionspace, dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
import material_properties

T_hp = 100
g = material_properties.g * material_properties.K

# Create mesh and define function space
mesh, cell_markers, facet_markers = gmshio.read_from_msh("trapezoid.msh", MPI.COMM_WORLD, gdim=2)
VT = functionspace(mesh, ("CG", 1))
ds = Measure('ds', domain=mesh, subdomain_data=facet_markers)

# Define boundary condition
facets = facet_markers.find(562)
fdim = mesh.topology.dim - 1
dofs = locate_dofs_topological(VT, fdim, facets)
bcT = dirichletbc(T_hp, dofs, VT)

# Define variational problem
T_ = TestFunction(VT)
dT = TrialFunction(VT)
aT = dot(grad(dT), grad(T_)) * dx
LT = -g * T_ * ds(563)

# Compute solution
Delta_T = Function(VT, name="Temperature")
problem = LinearProblem(aT, LT, bcs=[bcT], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Delta_T = problem.solve()

# obtain Temperature on  boundary 561
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = facet_markers.find(561)
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 = VT.dofmap.index_map.size_local+ VT.dofmap.index_map.num_ghosts
dof_to_geometry_map = np.full(num_dofs, -1, dtype=np.int32)
dofmap = VT.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

I don’t know what to do next.and I hope to get some help to achieve the same functionality as in my previous post (old version of dolfin)：https://fenicsproject.discourse.group/t/how-to-output-the-data-results-on-the-boundary/13227/2?u=french_fries
Thank you!

2

Couldn’t you just add mesh.geometry.x[node] to extract the coordinate for that given vertex?
It is unclear to me what part of the code you are struggling with porting.
You can access the equvialent value of the function with DeltaT.x.array[dof].

3

I would like to transfer the following Dolfin code to Dolfinx (outputting the coordinates and values of a specific boundary as a list)：

# obtain Temperature on boundary 2
facets = boundary.where_equal(2)
send_coordinates = []
send_values = []
send_global_index = []

global_vertices = mesh.topology().global_indices(0)
vertex_values = Delta_T.compute_vertex_values()
coordinates = mesh.coordinates()

for facet_index in facets:
    facet = Facet(mesh, facet_index)
    vertices = facet.entities(0)
    send_global_index.extend(vertices)
unique_vertices = np.unique(send_global_index)

for vertex in unique_vertices:
    send_coordinates.append(coordinates[vertex])
    send_values.append(vertex_values[vertex])

send_coordinates = np.asarray(send_coordinates).reshape(-1, coordinates.shape[1])

facet_values = np.zeros(len(unique_vertices))
facet_coordinates = np.zeros((len(unique_vertices), coordinates.shape[1]))
c = np.vstack(send_coordinates)
v = np.hstack(send_values)
i = np.hstack(unique_vertices)

for coord, value, index in zip(c, v, i):
    local_idx = np.flatnonzero(unique_vertices == index)[0]  
    facet_values[local_idx] = value
    facet_coordinates[local_idx] = coord

for (coord, value) in zip(facet_coordinates, facet_values):
        print(coord, value)

I have made some modifications, and it seems to be running. Could you please help me review the modified code for any potential issues? Thank you.Here is the modified code：

# obtain Temperature on boundary 2
facets = facet_markers.find(2)
send_coordinates = []
send_values = []

vertex_values = Delta_T.x.array
coordinates = mesh.geometry.x       
unique_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, facets, mesh.topology.dim-1, 0)

for vertex in unique_vertices:
    send_coordinates.append(coordinates[vertex])
    send_values.append(vertex_values[vertex])

send_coordinates = np.asarray(send_coordinates).reshape(-1, coordinates.shape[1])

facet_values = np.zeros(len(unique_vertices))
facet_coordinates = np.zeros((len(unique_vertices), coordinates.shape[1]))
c = np.vstack(send_coordinates)
v = np.hstack(send_values)
i = np.hstack(unique_vertices)

for coord, value, index in zip(c, v, i):
    local_idx = np.flatnonzero(unique_vertices == index)[0] 
    facet_values[local_idx] = value
    facet_coordinates[local_idx] = coord

for (coord, value) in zip(facet_coordinates, facet_values):
        print(coord, value)
4

Asking people to review partial code snippets is asking alot.

I would suggest you make a test where you know the answer to check the validity of your code.

5

Sure, I’ve done the test and everything should be fine. Similarly, I would like to implement this section of dolfin code in dolfinx (defining non-uniform boundary conditions)：https://fenicsproject.discourse.group/t/how-to-represent-non-uniform-first-and-second-type-boundary-conditions/13654/6?u=french_fries
Below is the modified code I’m using：

import dolfinx
import numpy as np
from mpi4py import MPI
from dolfinx import io, default_scalar_type
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem
from ufl import (TestFunction, TrialFunction, dot, dx, grad)
from dolfinx.fem import (Constant, dirichletbc, functionspace, Function, locate_dofs_topological)
from pathlib import Path

mesh, cell_markers, facet_markers = gmshio.read_from_msh("test.msh", MPI.COMM_WORLD, gdim=2)
VT = functionspace(mesh, ("CG", 1))
T_data = np.array([0, 100, 50, 80, 60])
T_coordinates = np.array([[0, 0, 0], [0, 0.25, 0], [0, 0.5, 0], [0, 0.75, 0], [0, 1, 0]])
Tf = Function(VT)

facets = facet_markers.find(6)
vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, facets, mesh.topology.dim-1, 0)
fdim = mesh.topology.dim - 1
boundary_dofs = locate_dofs_topological(VT, fdim, facets)

# Get all boundary 6 coordinates
send_coordinates = []
coordinates = mesh.geometry.x
for vertex in vertices:
    send_coordinates.append(coordinates[vertex])
boundary_coords = np.asarray(send_coordinates).reshape(-1, coordinates.shape[1])

# Compute distances between boundary coordinates and q_coordinates
points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(T_coordinates, 0)
distances = np.sum(np.square(points_A - points_B), axis=2)
is_close = distances < 1e2 * 1e-14
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    Tf.vector[boundary_dofs[row]] = T_data[col]
    print(boundary_coords[row][1], T_data[col])

bcT = dirichletbc(Tf, boundary_dofs, VT)

# Define variational problem
f = Constant(mesh, default_scalar_type(100))
T_ = TestFunction(VT)
dT = TrialFunction(VT)
aT = dot(grad(dT), grad(T_)) * dx
LT = f * T_ * dx

# Compute solution
Delta_T = Function(VT)
problem = LinearProblem(aT, LT, bcs=[bcT], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Delta_T = problem.solve()

The results show that the assignment is correct, but there appears to be a type error for Tf in the dirichletbc function：

1.0 60
0.75 80
0.5 50
0.25 100
0.0 0
Traceback (most recent call last):
  File "/home/dyfluid/Desktop/dolfinx/test.py", line 41, in <module>
    bcT = dirichletbc(Tf, boundary_dofs, VT)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/bcs.py", line 189, in dirichletbc
    bc = bctype(_value, dofs, V._cpp_object)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: __init__(): incompatible function arguments. The following argument types are supported:

Could you help me fix it?If necessary, I can provide the mesh file (2D unit rectangle, 5 nodes on each edge). Thank you.

6

Simply remove the latter argument, i.e.
bcT = dirichletbc(Tf, boundary_dofs)
as Tf is in VT

7

Thank you，it works very well！