How to export entity markers and nodal sources using io.XDMFFile

mesh
, ,
1

Hi,

Is there a way of writing the source vector (node loads in the code below) so they can be exported as an xdmf file and visualized in ParaView?

How can bc facets (like boundary_facets in the code below) be marked and exported with an xdmf file?

Thank you,
Alex

from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type, io
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
import numpy as np


mu = 1
lambda_ = 1.25


domain = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))


# Extract the coordinates of the nodes associated with the facets
node_coords = domain.geometry.x
print(str(node_coords) + '\n')



def clamped_boundary(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)


print(type(boundary_facets))


u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx

L = ufl.dot(T, v) * ds


a_compiled = fem.form(a)
L_compiled = fem.form(L)


# Assemble system, applying boundary conditions
A = assemble_matrix(a_compiled, bcs=[bc])
A.assemble()

b = assemble_vector(L_compiled)
apply_lifting(b, [a_compiled], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])

global_node = 3   # Chosen based on printed node coordinates
b.setValue(3*global_node,   1.0)
b.setValue(3*global_node+1, 1.0)
b.setValue(3*global_node+2, 1.0)


# Create solution function
uh = fem.Function(V)


# Solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)


# Compute solution
solver.solve(b, uh.x.petsc_vec)


with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
2

dolfinx comes with a set of demos, please have a look into them.
Search for the keyword above in all demos, and you’ll find the answer to both questions I read, i.e. how to mark facets and how to export them to file.

I know that my answer may be blunt, but with these demos you are given a set of tools to use as a starting point. In my opinion it is much more instructive for you to be informed that there is a place where you’ll be able to find answers for such common questions, rather than being given the answer to the specific common question you have asked this time.

I invite everyone else not to reply with the code.
Once @alex24 finds the answer in the demo the community will be grateful if he posted the updated code.