How to apply external volume source

General
1

Hello,
I have a question about external volume source.

The weak formulation of my problem looks like this
image
where the highlighted term represents the volume source, which I have available in the form of raw data like this

# x y z  source_x source_y source_z
0.0525 0.0375 0 -0.955421 223.385 0
0 0.0375 0 -1541.5 1497.81 0
0.0525 0 0 0 0 0
0 0 0 -1177.21 531.919 0
0.0525 0.075 0 1.06877 297.886 0
0 0.075 0 -494.927 2098.87 0
.
.
.

This data was obtained on the same mesh as is used in fenics.

Here is part of my code implementing this

# mesh
msh, gmsh_cell_tags, gmsh_facet_tags = io.gmshio.read_from_msh("mesh.msh", MPI.COMM_WORLD, gdim=2)

# finite element space
FE_source = basix.ufl.element("Lagrange", msh.basix_cell(), degree=1, shape=(msh.geometry.dim,))
V_source = fem.functionspace(msh, FE_source)
source = fem.Function(V_source)

FE = basix.ufl.element("Lagrange", msh.basix_cell(), degree=FE_SPACE_DEGREE)
V = fem.functionspace(msh, FE)
(u, v) = ufl.TrialFunction(V), ufl.TestFunction(V) 
uh = fem.Function(V)

# constants
c = PETSc.ScalarType(343) 
dt = 1e-4 
t = fem.Constant(msh, 0.0)
T = 0.1

# central scheme
u_nn = fem.Function(V)
u_n = fem.Function(V)

dudtdt = (u - 2 * u_n + u_nn) / dt**2
dudt = (u - u_nn) /(2*dt)
du = (u + u_nn) / 2

# variational problem
F = 1/c**2 * ufl.inner(dudtdt, v) * ufl.dx \
  + 1/c * ufl.inner(dudt, v) * ds(1) \
  + ufl.inner(ufl.grad(du), ufl.grad(v)) * ufl.dx \
  + ufl.dot(source, ufl.grad(v)) * ufl.dx 
    
a, L = ufl.system(F)

# applying volume source
ids = msh.geometry.input_global_indices
dim = msh.topology.dim

data = np.loadtxt("source.raw")
source.x.array[:] = data[ids, 3:3+dim].flatten()

# solving problem
petsc_options = {"ksp_type": "preonly",
                 "pc_type": "lu", 
                 "pc_factor_mat_solver_type": "mumps"}
problem = dolfinx.fem.petsc.LinearProblem(a, L, u=uh, bcs=[], petsc_options=petsc_options)

vtx_u = io.VTXWriter(msh.comm, "results.bp", u_n, engine="BP4") 

while t.value < T:
    problem.solve()
    
    u_nn.x.array[:] = u_n.x.array
    u_n.x.array[:] = uh.x.array
      
    vtx_u.write(t.value)
      
    t.value += dt

(I have left out a few lines of codes for the sake of clarity and I am aware that the volume source is the same for all the time steps)

I visually checked if the mapping of the volume source is done right in the domain using pyvista and all seems good. But I still wonder if this approach is correct, because I am getting unexpected results and I want to be sure this is implemented correctly, before I go to look for problem elsewhere.

2

Please search the forum for similar questions, for instance Fenicsx: use an array of data for boundary conditions (with a couple of links) or Mapping 2D numpy array into dolfinx function

I wouldn’t trust

# applying volume source
ids = msh.geometry.input_global_indices
dim = msh.topology.dim

data = np.loadtxt("source.raw")
source.x.array[:] = data[ids, 3:3+dim].flatten()

as there is nothing that guarantees you that the order of the points in the file is the same at the locations of the DOFs.

3

Okay, I see your point, but I checked the source visually like this

with io.XDMFFile(msh.comm, "source_check.xdmf", "w") as xdmf:
   xdmf.write_mesh(msh)
   xdmf.write_function(source)

and it looked as I would expect. Similarly, I checked the coordinates of few random vertices with the following code

msh.geometry.x[5101] # = [4.8, 0.1875, 0.0]
data[ids[5101], :] # = [4.8 0.1875 0 0.00147893 0.0485567 0]

and it matched everytime. So is there any other way I could check for certain?

In any case, my main question wasn’t about the mapping, so let’s assume the source mapping is correct, does the rest of the implementation look good?"

4

See for instance: Wrong order of dofs after creating mesh from lists of points and cells - #2 by dokken
combined with
Application of point forces, mapping vertex indices to corresponding DOFs - #8 by Henrik_Nicolay_Finsb
Which would Give your the full: dof->vertex->input geometry node mapping

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?