How to look at slab data?

[prusso630]

I have a mesh script where I entered the dimension for a concrete slab and three steel bars in inches:

import gmsh 
from mpi4py import MPI
gmsh.initialize()
proc = MPI.COMM_WORLD.rank 

concrete_marker = 11
rebar_marker = 12
if proc == 0:
    gmsh.model.add("concreteslab")

    concrete = gmsh.model.occ.add_box(-12, -3, -72, 24, 6, 144)
    bary = 2


    gmsh.option.setNumber('Mesh.MeshSizeFromCurvature', 5)
    cylinder1 = gmsh.model.occ.addCylinder(-10, -bary, 72,  0, 0, -144,   .5)  
    cylinder2 = gmsh.model.occ.addCylinder(0, -bary, 72,  0, 0, -144,   .5)   
    cylinder3 = gmsh.model.occ.addCylinder(10, -bary, 72,  0, 0, -144,   .5)  
    gmsh.model.occ.synchronize()


    volTags = gmsh.model.getEntities(3)
    gmsh.model.occ.fragment(volTags, volTags)


    gmsh.model.addPhysicalGroup(3, [concrete], concrete_marker)
    gmsh.model.addPhysicalGroup(3, [cylinder1, cylinder2, cylinder3], rebar_marker)

    gmsh.model.mesh.generate(3)
    gmsh.write("mesh3D.msh")
gmsh.finalize()

I modified demo_elasticity so that the mesh of the slab (mesh3D.msh) can pass through the python system:

from contextlib import ExitStack
import numpy as np
import ufl
from dolfinx import la
from dolfinx.fem import (Expression, Function, FunctionSpace,
                         VectorFunctionSpace, dirichletbc, form,
                         locate_dofs_topological, locate_dofs_geometrical)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_box,
                          locate_entities_boundary)
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
dtype = PETSc.ScalarType
from dolfinx.io import gmshio
import meshio
from petsc4py import PETSc
import pyvista 
from dolfinx import plot 

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

def build_nullspace(V):
    index_map = V.dofmap.index_map
    bs = V.dofmap.index_map_bs
    ns = [la.create_petsc_vector(index_map, bs) for i in range(6)]
    with ExitStack() as stack:
        vec_local = [stack.enter_context(x.localForm()) for x in ns]
        basis = [np.asarray(x) for x in vec_local]


        dofs = [V.sub(i).dofmap.list.array for i in range(3)]

 
        for i in range(3):
            basis[i][dofs[i]] = 1.0

 
        x = V.tabulate_dof_coordinates()
        dofs_block = V.dofmap.list.array
        x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
        basis[3][dofs[0]] = -x1
        basis[3][dofs[1]] = x0
        basis[4][dofs[0]] = x2
        basis[4][dofs[2]] = -x0
        basis[5][dofs[2]] = x1
        basis[5][dofs[1]] = -x2


    la.orthonormalize(ns)
    assert la.is_orthonormal(ns)

    return PETSc.NullSpace().create(vectors=ns)

slabdomaingmshio, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3) 

slabmeshio = meshio.read("mesh3D.msh")


tetra_mesh = create_mesh(slabmeshio, "tetra", prune_z=False)

meshio.write("mesh.xdmf", tetra_mesh)

with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    slabdomaingmshio = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(slabdomaingmshio, name="Grid")

Q = FunctionSpace(slabdomaingmshio, ("DG", 0))

Emod = Function(Q)
concrete = ct.find(11)
Emod.x.array[concrete] = np.full_like(concrete, 4.3e6, dtype=PETSc.ScalarType)
steel = ct.find(12)
Emod.x.array[steel]  = np.full_like(steel, 30e6, dtype=PETSc.ScalarType)

VMod = Function(Q)
VMod.x.array[concrete] = np.full_like(concrete, 0.15, dtype=PETSc.ScalarType) 
VMod.x.array[steel]  = np.full_like(steel, 0.28, dtype=PETSc.ScalarType)

ω, ρ = 5000, 33.3
x = ufl.SpatialCoordinate(slabdomaingmshio)
f = ufl.as_vector((ρ * ω**2 * x[0], ρ * ω**2 * x[1], 0.0))


μ = Emod / (2.0 * (1.0 + VMod))
λ = Emod * VMod / ((1.0 + VMod) * (1.0 - 2.0 * VMod))

def σ(v):
    """Return an expression for the stress σ given a displacement field"""
    return 2.0 * μ * ufl.sym(grad(v)) + λ * ufl.tr(ufl.sym(grad(v))) * ufl.Identity(len(v))



V = VectorFunctionSpace(slabdomaingmshio, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = form(inner(σ(u), grad(v)) * dx)
L = form(inner(f, v) * dx)





dofs_L = locate_dofs_geometrical(V, lambda x: np.isclose(x[2], -72))
u_L = Function(Q)
bc_L = dirichletbc(PETSc.ScalarType(0), dofs_L,Q)

dofs_R = locate_dofs_geometrical(V, lambda x: np.isclose(x[2], 72))
u_R = Function(Q)
bc_R = dirichletbc(PETSc.ScalarType(0), dofs_R,Q)
bcs = [bc_L,bc_R]



A = assemble_matrix(a, bcs=bcs)
A.assemble()

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


null_space = build_nullspace(V)
A.setNearNullSpace(null_space)

opts = PETSc.Options()
opts["ksp_type"] = "cg"
opts["ksp_rtol"] = 1.0e-10
opts["pc_type"] = "gamg"
opts["mg_levels_ksp_type"] = "chebyshev"
opts["mg_levels_pc_type"] = "jacobi"
opts["mg_levels_esteig_ksp_type"] = "cg"
opts["mg_levels_ksp_chebyshev_esteig_steps"] = 20

solver = PETSc.KSP().create(slabdomaingmshio.comm)
solver.setFromOptions()


solver.setOperators(A)

uh = Function(V)


solver.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
solver.solve(b, uh.vector)
solver.view()


uh.x.scatter_forward()



sigma_dev = σ(uh) - (1 / 3) * ufl.tr(σ(uh)) * ufl.Identity(len(uh))
sigma_vm = ufl.sqrt((3 / 2) * inner(sigma_dev, sigma_dev))

W = FunctionSpace(slabdomaingmshio, ("Discontinuous Lagrange", 0))
sigma_vm_expr = Expression(sigma_vm, W.element.interpolation_points())
sigma_vm_h = Function(W)
sigma_vm_h.interpolate(sigma_vm_expr)

with XDMFFile(slabdomaingmshio.comm, "out_elasticity/displacements.xdmf", "w") as file:
    file.write_mesh(slabdomaingmshio)
    file.write_function(uh)


with XDMFFile(slabdomaingmshio.comm, "out_elasticity/von_mises_stress.xdmf", "w") as file:
    file.write_mesh(slabdomaingmshio)
    file.write_function(sigma_vm_h)

unorm = uh.x.norm()
if slabdomaingmshio.comm.rank == 0:
    print("Solution vector norm:", unorm)


#pyvista.start_xvfb()
#pyvista.OFF_SCREEN = True
# Create plotter and pyvista grid
#p = pyvista.Plotter()
#topology, cell_types, geometry = plot.create_vtk_mesh(V)
#grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

# Attach vector values to grid and warp grid by vector
#grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
#actor_0 = p.add_mesh(grid, style="wireframe", color="k")
#warped = grid.warp_by_vector("u", factor=1.5)
#actor_1 = p.add_mesh(warped, show_edges=True)
#p.show_axes()
#if not pyvista.OFF_SCREEN:
#   p.show()
#else:
#   figure_as_array = p.screenshot("deflection.png")

#print('done')

This leads to a number of questions I have so far:

• I would like to enter a distributed load in pounds per foot. What are the units of w in this case? How can I enter a distributed load in pounds per foot?
• I would like to look at stress from the exported von_misses_stress.xdmf in PSI. I am attaching a screenshot of what I see in ParaView at the end of this post. So far units for stress (sigma) look nothing in the PSI range. How can I get an export of von_misses_stress.xdmf in PSI?
• In terms of displacements.xdmf there is one being generated. How is it possible in Paraview to get annotations of either current x,y,z of the displacement or an image map showing the largest displacement that occurs from Paraview?
• In terms of the commented out PyVista code at then end of the script, I tried running it in conda,spack and also container environments. deflection.png is the type of image I want for a deflection map. Things seem to lock up at:

p.show() <<-- seems to lock up the system
figure_as_array = p.screenshot("deflection.png") <<--also seems to lock up the system

[prusso630]

# Pounds per square inch of a 24"x144" slab, weight of concrete per square inch
# There may be a better value for density.
ω, ρ = 0.57, 0.087

I see that there is a way to slice in paraview that works well. Also there is a way to plot wireframe of displacements that works well.

[prusso630]

There is also a related post about the type of function space to use for multiple materials:

How to set up a function space consisting of multiple domains from an imported gmsh mesh for an elasticity analysis

You are assuming a 1 to 1 relationship between the dofs of a DG-1 space and the cell indices. This is not correct, as there is only a 1-1 relationship between a DG-0 space and the cell indices. So i would suggest using a DG-0 space. See also part two of the tutorial: Defining subdomains for different materials — FEniCSx tutorial which creates materials based on cell data from the mesh markers.