Find the coordinate of DOFs on a second-order mesh

mint_Yzr · October 14, 2024, 8:28am

How can I correctly obtain the coordinates associated with degrees of freedom (DOFs) when using a second-order mesh in Dolfinx?

I couldn’t find support for periodic boundary conditions in Dolfinx, so I attempted to manually modify the PETSc matrix assembled by Dolfinx. This requires knowing the coordinates of the DOFs to pair the DOFs on the left boundary with those on the right boundary. However, this approach has become quite messy.

I have created an example that demonstrates that using V.tabulate_dof_coordinates() does not seem to correctly retrieve the coordinates of the DOFs when using a second-order mesh. In this example, I attempt to get the DOFs on the left boundary, and print the DOFs I’ve got using a function. The resulting image is strange.

import gmsh
import numpy as np
from dolfinx.io import gmshio, XDMFFile
from mpi4py import MPI
import dolfinx
import basix

affix = "./result_test_move/"
r = 0.125
L = 1
H = L
def create_split_mesh(c_x, c_y):
    gmsh.initialize()
    gdim = 2
    mesh_comm = MPI.COMM_SELF
    model_rank = 0
    if mesh_comm.rank == model_rank:
        rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, L, tag=1)
        obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

        p1 = gmsh.model.occ.add_point(L / 2, 0, 0)
        p2 = gmsh.model.occ.add_point(L / 2, L, 0)
        l1 = gmsh.model.occ.add_line(p1, p2)

        gmsh.model.occ.synchronize()


    fluid_marker = 1
    if mesh_comm.rank == model_rank:
        volumes = gmsh.model.getEntities(dim=gdim)
        gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
        gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

    inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
    inflow, outflow, walls, obstacle = [], [], [], []
    if mesh_comm.rank == model_rank:
        boundaries = gmsh.model.getBoundary(volumes, oriented=False)
        for boundary in boundaries:
            center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
            if np.allclose(center_of_mass, [0, H / 2, 0]):
                inflow.append(boundary[1])
            elif np.allclose(center_of_mass, [L, H / 2, 0]):
                outflow.append(boundary[1])
            elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
                walls.append(boundary[1])
            else:
                obstacle.append(boundary[1])
        gmsh.model.addPhysicalGroup(1, walls, wall_marker)
        gmsh.model.setPhysicalName(1, wall_marker, "Walls")
        gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
        gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
        gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
        gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
        gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
        gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

    if mesh_comm.rank == model_rank:
        gmsh.model.mesh.embed(1, obstacle, 2, rectangle)
        # gmsh.model.mesh.embed(1, [l1], 2, rectangle)
        gmsh.model.occ.fragment([(2, rectangle)], [(1, l1)])

        gmsh.option.setNumber("Mesh.MeshSizeMin", r / 4)
        gmsh.option.setNumber("Mesh.MeshSizeMax", r / 3)

        gmsh.option.setNumber("Mesh.Algorithm", 8)
        gmsh.model.mesh.generate(gdim)
        gmsh.model.mesh.setOrder(2)
    
    msh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

    gmsh.finalize()
    return msh

msh0 = create_split_mesh(0.25, 0.4)

V_element = basix.ufl.element("CG", msh0.basix_cell(), 2, shape=(2,))
V = dolfinx.fem.functionspace(msh0, V_element)

dof_coord = V.tabulate_dof_coordinates()[:,:2]
def pbs_target(x):
    return x[:,0] < 0 + 1e-6
target = pbs_target(dof_coord).nonzero()[0].astype(np.int32)

f_in = dolfinx.fem.Function(V)
f_in.x.array[target] = 10000
with XDMFFile(MPI.COMM_WORLD, f"{affix}solution.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh0)
    xdmf.write_function(f_in)

dokken · October 14, 2024, 1:45pm

Periodic boundary conditions is supported using the extension DOLFINx_MPC:

See for instance:

github.com

jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic3d_topological.py

# This demo program solves Poisson's equation
#
#     - div grad u(x, y) = f(x, y)
#
# on the unit square with homogeneous Dirichlet boundary conditions
# at y = 0, 1 and periodic boundary conditions at x = 0, 1.
#
# Original implementation in DOLFIN by Kristian B. Oelgaard and Anders Logg
# This implementation can be found at:
# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/periodic/demo_periodic.py
#
# Copyright (C) Jørgen S. Dokken 2020-2022.
#
# This file is part of DOLFINX_MPCX.
#
# SPDX-License-Identifier:    MIT
from __future__ import annotations

from pathlib import Path
from typing import Dict, Union

This file has been truncated. show original

github.com

jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_geometrical.py

# This demo program solves Poisson's equation
#
#     - div grad u(x, y) = f(x, y)
#
# on the unit square with homogeneous Dirichlet boundary conditions
# at y = 0, 1 and periodic boundary conditions at x = 0, 1.
#
# Copyright (C) Jørgen S. Dokken 2020-2022.
#
# This file is part of DOLFINX_MPCX.
#
# SPDX-License-Identifier:    MIT
from __future__ import annotations

from pathlib import Path
from typing import Union

from mpi4py import MPI
from petsc4py import PETSc

This file has been truncated. show original

github.com

jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py

# Copyright (C) 2021-2022 fmonteghetti, Connor D. Pierce, and Jorgen S. Dokken
#
# This file is part of DOLFINX_MPC
#
# SPDX-License-Identifier:    MIT
# This demo solves the Laplace eigenvalue problem
#
#     - div grad u(x, y) = lambda*u(x, y)
#
# on the unit square with two sets of boundary conditions:
#
#   u(x=0) = u(x=1) and u(y=0) = u(y=1) = 0,
#
# or
#
#   u(x=0) = u(x=1) and u(y=0) = u(y=1).
#
# The weak form reads
#
#       (grad(u),grad(v)) = lambda * (u,v),

This file has been truncated. show original

For your specific problem, you haven’t mapped the degrees of freedom by their block size.
If you look at the shape of dof_coord,

you will realize that it is half the size of u.x.array.
This is because we have that dof (2i, 2i+1) share a dof coordinate in DOLFINx,
i.e.


def pbs_target2(x):
    return x[0] < 0 + 1e-6


dofs = dolfinx.fem.locate_dofs_geometrical(V, pbs_target2)
bs = V.dofmap.bs
unrolled_dofs = np.empty(len(dofs)*bs,dtype=np.int32)
for (i, dof) in enumerate(dofs):
    for b in range(bs):
        unrolled_dofs[i*bs+b] = dof*bs+b


f_in = dolfinx.fem.Function(V)
f_in.x.array[unrolled_dofs] = 10000
with XDMFFile(MPI.COMM_WORLD, f"{affix}solution.xdmf", "w") as xdmf:
    xdmf.write_mesh(msh0)
    xdmf.write_function(f_in)

will give you the correct result

mint_Yzr · October 14, 2024, 1:53pm

Wow, that helps a lot. Thank you for your prompt, informative, and precise reply.

Topic		Replies	Views
Unable to properly locate dofs dolfinx	3	515	March 30, 2021
Getting coordinates of facets in dolfinx dolfinx	2	2433	February 24, 2021
How can I output the coordinates of a specific boundary and its corresponding computation results as two lists in dolfinx? dolfinx	6	250	May 29, 2024
SpatialCoordinate of submesh dolfinx mesh	2	23	November 27, 2024
Non-matching coordinate ordering for CG1 functions defined on submeshes dolfinx mesh , dolfinx	2	22	November 22, 2024

Find the coordinate of DOFs on a second-order mesh

Related topics