How to apply external volume source

General
5

I looked into both topics, but I’m still uncertain about how to proceed. I’ll try to explain my issue using a simpler example.

Let’s consider a 2D mesh with 16 vertices and a vector-valued finite element function space of first order

from mpi4py import MPI
import dolfinx
import basix.ufl

# mesh
quadrilateral = dolfinx.mesh.CellType.quadrilateral
msh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 3, 3, cell_type=quadrilateral)

num_vertices = msh.topology.index_map(0).size_local
num_cells = msh.topology.index_map(msh.topology.dim).size_local

# function space
FE = basix.ufl.element("Lagrange", msh.basix_cell(), degree=1, shape=(msh.geometry.dim,))
V = dolfinx.fem.functionspace(msh, FE)

image
image755×603 18.4 KB

And then I have the external data for each vertex in the following format

# x y z f_x f_y f_z
0.33 0.33 0 2.2 1 0
0 0.66 0 3.1 1 0
0.66 0 0 1.3 1 0
1 0.33 0 2.4 1 0
1 1 0 4.4 1 0
0.66 0.33 0 2.3 1 0
0 0 0 1.1 1 0
0.66 0.66 0 3.3 1 0 
1 0.66 0 3.4 1 0
0 0.33 0 2.1 1 0
0.33 0 0 1.2 1 0 
0 1 0 4.1 1 0
0.66 1 0 4.3 1 0
1 0 0 1.4 1 0 
0.33 0.66 0 3.2 1 0
0.33 1 0 4.2 1 0

which I read as

import numpy as np
data = np.loadtxt("file.dat")[:, 3:]

The values for each vertex are in random order and my goal is to assign them to

source = fem.Function(V)
source.x.array[:] = data[?, :mesh.geometry.dim].flatten()

accordingly, so the coordinates of external values match the coordinates of DOFs.

So from application-of-point-forces-mapping-vertex-indices-to-corresponding-dofs I tried

def vertex_to_dof_mapping(mesh, V):
    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    c2v = msh.topology.connectivity(msh.topology.dim, 0)

    for cell in range(num_cells):
        vertices = c2v.links(cell)
        dofs = V.dofmap.cell_dofs(cell)

        for i, vertex in enumerate(vertices):
            vertex_to_dof_map[vertex] = dofs[V.dofmap.dof_layout.entity_dofs(0, i)]

    return vertex_to_dof_map

vertex_to_dof_map = vertex_to_dof_mapping(msh, V)

This i guess should answer the dof->vertex mapping, however in my case this returns an array from 0 to 15, without any changes in the ordering?

And in case of vertex->input geometry mapping, thats where I just cycle through the external file and find the corresponding indices, where the coordinates match with mesh.geometry.x, right?