Confused about how to extract solution fields in MixedFunctionSpace

dolfinx
1

Hi everyone,
I’m relatively new to Dolfinx and am currently working on my master’s thesis, where I’m attempting to solve a Cosserat elasticity problem using a mixed finite element formulation.

My goal is to extend an existing 2D example into a 3D setting, where the geometry consists of a beam fixed at its base and subjected to a constant compressive force along the z-axis.

The constitutive relations for the stress tensor \bar{\sigma} and the couple stress tensor \bar{m}^R are given by:

\begin{align*} \bar{\boldsymbol{\sigma}} &= \frac{\partial \psi_e}{\partial \bar{\boldsymbol{\varepsilon}}} = \frac{\partial}{\partial \bar{\boldsymbol{\varepsilon}}} \left( \psi_e^B + \psi_e^C \right) = \underbrace{\lambda \left( \mathrm{tr} \, \bar{\boldsymbol{\varepsilon}}^{\text{sym}} \right) \mathbf{I} + (2\mu + \kappa) \bar{\boldsymbol{\varepsilon}}^{\text{sym}}}_{:= \bar{\boldsymbol{\sigma}}^B} + \underbrace{\kappa \bar{\boldsymbol{\varepsilon}}^{\text{skew}}}_{:= \bar{\boldsymbol{\sigma}}^C}, \end{align*}
\begin{align*} \bar{\boldsymbol{m}}^R &= \frac{\partial \psi_e}{\partial \bar{\boldsymbol{\phi}}} = \frac{\partial \psi_e^R}{\partial \bar{\boldsymbol{\phi}}} = \gamma \bar{\boldsymbol{\phi}}. \end{align*}

The weak formulation of the problem is:

\begin{align*} \int_{B} \boldsymbol{\bar{\varepsilon}}(\boldsymbol{\eta}, \boldsymbol{\zeta}) : \boldsymbol{\bar{\sigma}}(\boldsymbol{u}, \theta_3) \, dV - \int_{\partial B_{t_{\sigma}}} \boldsymbol{\eta} \cdot \hat{\boldsymbol{t}}_{\sigma} \, dA &= 0, \end{align*}
\begin{align*} \int_{B} \boldsymbol{\bar{\phi}}(\boldsymbol{\zeta}) \cdot {\boldsymbol{\bar{m}}}^{R}(\theta_3) \, dV - \int_{B} \boldsymbol{\zeta} : \mathbb{E} : \boldsymbol{\bar{\sigma}}^{C}(\boldsymbol{u}, \theta_3) \, dV - \int_{\partial B_{t_{m}}} \boldsymbol{\zeta} : \hat{\boldsymbol{t}}_{m} \, dA &= 0. \end{align*}

However, I’m a bit confused about mixed spaces and how to extract the solution after solving the problem.

Although the code works correctly for the deformation part of the problem, I’m still not obtaining meaningful results for the microrotations. When I extract the microrotation field using: theta_values = uh.sub(1).collapse().x.array
I’m able to save and plot the solution, but all the values are zero. On the other hand, if I use:

u, theta = uh.split()
theta_values = theta.x.array

I do obtain non-zero values for the microrotation field, but I’m unable to save or visualize them properly.

I’m not sure whether this issue stems from the boundary conditions I’ve imposed or from a potential issue in the formulation of the mixed function space.
Any help or suggestions would be greatly appreciated.
Thanks in advance for your time.

Here’s my MWE running on Ubuntu 22.04 and dolfinx 0.9

from dolfinx import mesh, fem,plot, io
from dolfinx.fem.petsc import LinearProblem
from ufl import (Identity, Measure, PermutationSymbol,
                nabla_grad, tr, TrialFunctions, TestFunctions,  
                dot, inner, sym, dx)
from basix.ufl import element, mixed_element
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import pyvista as pv
import numpy as np
# units are in N, tonn, mm, s, mJ
Nx = 5 # Ny = Nx = Nz/3, number of cells in x direction
a_x = 1.5e-4  # cell costant for x,y
a_z = 2/3 * a_x  # cell costant for z

G= 0.1  # Shear modulus [MPa]
nu= -0.125
rho= 1
l= 4*a_x  # Characteristic length [mm]
N= 0.9     # Coupling number
# Computing LAME costants
lame_lambda = G*((2*nu)/(1-2*nu))
lame_mu =  G*((1-2*N**2)/(1-N**2))
# Computing COSSERAT constants
gamma = 4*G*l**2
kappa = G*((2*N**2)/(1-N**2))

### Mesh and and Mixed Function Space --------------------
msh  = mesh.create_box(
    comm=MPI.COMM_WORLD, 
    points=([np.array([0, 0,0]),np.array([Nx*a_x, Nx*a_x, H:=2*Nx*a_x])]),
    n= (5,5,15),
    cell_type= mesh.CellType.hexahedron         
    )

dim = msh.geometry.dim
V_el = element("Lagrange", msh.topology.cell_name(), 2, shape=(dim,))
W_el = element("Lagrange", msh.topology.cell_name(), 1, shape=(dim,))
VW_el = mixed_element([V_el, W_el])
V = fem.functionspace(msh, VW_el)

V0, V0_dofs = V.sub(0).collapse()
V1, V1_dofs = V.sub(1).collapse()

(u_trial, theta_trial) = TrialFunctions(V) 
(v_test, zeta_test) = TestFunctions(V)
### Boundary Conditions for displacement and rotation for z = 0---------------------------
facets_dim = msh.topology.dim - 1

boundary_facets = mesh.locate_entities_boundary(
    msh = msh,
    dim= facets_dim, 
    marker=lambda x: np.isclose(x[2], 0.0) 
    )

dofs = fem.locate_dofs_topological(V=(V.sub(0), V0),entity_dim=facets_dim,entities=boundary_facets)

fixed_displ_expr= lambda x: np.zeros((3, x.shape[0]), dtype=ScalarType)
fixed_displ = fem.Function(V0)
fixed_displ.interpolate(fixed_displ_expr)

fixed_rot_expr= lambda x: np.zeros((3, x.shape[1]), dtype=ScalarType)
fixed_rot = fem.Function(V1)
fixed_rot.interpolate(fixed_rot_expr)

# bc for displacement and micro-rotation for z= 0 (x[2])
bc_displ = fem.dirichletbc(value=fixed_displ,dofs=dofs,V=V.sub(0))
bc_rot = fem.dirichletbc(value=fixed_rot,dofs=dofs,V=V.sub(1))

bcs= [bc_displ, bc_rot]

### Micropolar strain:

E_3 = PermutationSymbol(len(u_trial))  # Levi-Civita tensor for n = 3

ϵ = lambda u, theta : sym(nabla_grad(u)) - dot(E_3,theta)  # micropolar strain
ϵ_sym = lambda u : sym(nabla_grad(u))
ϵ_skew = lambda u, theta : sym(nabla_grad(u)) - dot(E_3, theta) 

ϕ = lambda theta : nabla_grad(theta)  # micro-curvature
### Stresses
def σ_B(u):
    """Evaluating Boltzmann part of stress tensor sigma"""
    e_sym = ϵ_sym(u)
    stress = lame_lambda*tr(e_sym)*Identity(len(u)) + \
                (2*lame_mu + kappa)*e_sym
    return stress
σ_C   = lambda u, theta: kappa*ϵ_skew(u, theta)  # micro-continuum coupling
σ   = lambda u, theta: σ_B(u) + σ_C(u, theta)  # force stress

m_R = lambda θ : gamma*ϕ(θ)  # couple stress

# Defining weak form -----------------------------------
forcing = fem.Constant(msh,ScalarType((0.0,0.0, -10)))
traction = fem.Constant(msh, ScalarType(( 0.0, 0.0, 0.0)))

ds = Measure("ds", domain=msh)

### left and right hand side of the weak form formulation--------------
weak_form_lhs = (
    inner(σ(u_trial, theta_trial), ϵ(v_test, zeta_test)) * dx +
    inner(ϕ(zeta_test),m_R(theta_trial)) * dx -
    inner(dot(zeta_test, E_3), σ_C(u_trial, theta_trial)) * dx
    )
weak_form_rhs = (dot(forcing,v_test) * dx +dot(traction, v_test) *ds)

### Compute solution --------------------------------------------------- 
solution = LinearProblem(
                        a=weak_form_lhs, 
                        L = weak_form_rhs,
                        bcs= bcs,
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
                        )

uh = solution.solve()
u, theta = uh.sub(0).collapse(), uh.sub(1).collapse()

### Displaying the micro-rotation using pv------------------------------------- 
p = pv.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(V1)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

grid["theta"] = theta.x.array.reshape(geometry.shape[0], 3)
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("theta", factor=10)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pv.OFF_SCREEN:
    p.show()

### saving solution to a .pvd file
with io.VTKFile(MPI.COMM_WORLD, "solution/micro_rotation.pvd", "w") as vtk:
    vtk.write_function(theta, 0.0)  # 0.0 = time step
2

Use u = uh.sub(0).collapse(), theta=uh.sub(1).collapse() to get proper functions in the relevant sub-space (which can then be visualized).