Magnetostatics example in dolfinx

dolfinx
2

For projection, you can consider the approach taken in dolfiny by either copying this approach or installing this convenience library.

Also it is beneficial to have a look at:

Finally, I have made a minimal example highlighting the combination of these operations:

import dolfinx
import dolfinx.io

import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 5,5, dolfinx.cpp.mesh.CellType.quadrilateral)
num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
# Create a meshtag for the first and last local indices of the mesh
indices = np.array([0,1,2,3, num_cells_local-4, num_cells_local-3, num_cells_local-2, num_cells_local-1], dtype=np.int32)
values = np.array([1,1,1,1,2,2,2,2], dtype=np.int32)
mt = dolfinx.MeshTags(mesh, mesh.topology.dim, indices, values)


# Create a DG function for these cells
V = dolfinx.FunctionSpace(mesh, ("DG",0))
v = dolfinx.Function(V)

for i in range(num_cells_local):
    if i in mt.indices:
        local_index = np.flatnonzero(mt.indices == i)[0]
        if mt.values[local_index] == 1:
            v.vector.setValueLocal(i, 1.1)
        elif mt.values[local_index] == 2:
            v.vector.setValueLocal(i, 2.1)
    else:
        v.vector.setValueLocal(i, 0.1)


# Project DG function to CG space
W= dolfinx.FunctionSpace(mesh, ("CG",1))
u, w = ufl.TrialFunction(W), ufl.TestFunction(W)
a= ufl.inner(u, w) * ufl.dx
l = ufl.inner(v, w) * ufl.dx
uh = dolfinx.Function(W)
A = dolfinx.fem.assemble_matrix(a)
A.assemble()
b = dolfinx.fem.assemble_vector(l)

solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(A)
solver.solve(b, uh.vector)
uh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Write mesh tags and projected function to output files
xdmf_file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w")
xdmf_file.write_mesh(mesh)
mesh.topology.create_connectivity(mesh.topology.dim-1,mesh.topology.dim)
xdmf_file.write_meshtags(mt)
xdmf_file.close()
dolfinx.io.VTKFile("projectedDG.pvd").write(uh)
3 Likes