Hi all,
I am working on a simple elasticity problem and I would like to add a new variable to that which should be equal to the displacment gradient:
here is the variational formulation:
πππ: the stress tensor, πππ : the strain tensor, and π’π : the displacement field. Ξ¨jk is the new variable that should equals the π’π,π (du/dxj). I used Lagrange multiplier to force this relation.
Since I am using itertaive solvers, I need to set up the null space for this variational form.
I think I have to provide 18 vectors that span the null space:( 6 for rigid body motions, 6 for gradient of the displacement field π’π,π,and 6 arbitrary functions for πΌππ).
I followed the post here to set up the null space and was able to do that. But I donβt know how to make sure this is the correct null space for my problem. What are the option to make sure that I have set up the correct nulll space?
I tried to test it with
assert nullspace.test(A)
but I got the assertaion error.
Here is MWE:
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx.fem import (Function, functionspace, dirichletbc, form, locate_dofs_topological, Constant)
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.mesh import locate_entities_boundary, create_box, CellType, create_unit_square
from ufl import (Identity, TestFunctions, TrialFunctions, dx, indices, as_tensor)
from basix.ufl import element, mixed_element
from dolfinx import la
from dolfinx import default_scalar_type
import dolfinx
Fl = 1e-8
Fs = 1e-5
Fp = 1e-4
# ....................
C_star = Fl ** 2 / Fs
h_star = 1 / Fl
F_star = Fs
# ........................................ Parameters
E = 1 * 152e9 * C_star
nu = 0.33
c1 = E * nu / ((1 + nu) * (1 - 2 * nu))
c2 = E / (2 * (1 + nu))
# ......................................... Reading mesh file
h = 0.76e-3 / Fl
H = h
mesh = create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],[1, 1, 1], cell_type=CellType.hexahedron)
# ......................................... Define Elements
v = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)) # Displacement Element
s = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3)) # Displacement gradient Element
lm = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3)) # Lagrange multipliers Element
# ......................................... Define mixed element
MIXEL = mixed_element([v, s, lm])
Space = functionspace(mesh, MIXEL)
# ......................................... Define Trail and Test Functions
(u, gradu, lamu) = TrialFunctions(Space)
(del_u, del_gradu, del_lamu) = TestFunctions(Space)
# ....................................... Mark boundary regions
tol = 1e-16
def Bottom(x):
return np.isclose(x[2], 0, atol=tol)
# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
clamp = Function(S0)
clamp.x.array[:] = 0
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
# .........Fixed BC
clamp_dofs = locate_dofs_topological((Space.sub(0), S0), mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(clamp, clamp_dofs, Space.sub(0))
bc = [bc0]
# ......................................... Define tensor
delta = Identity(len(u))
i, j, k = indices(3)
# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])
def StressT(u):
return as_tensor((c1 * StrainT(u)[i, i] * delta[j, k] + 2 * c2 * StrainT(u)[j, k]), [j, k])
def Constraint(u, gradu): # CONSTRAINT
return as_tensor(gradu[j, k] - u[k].dx(j), (j, k))
# ......................................... Define Variational Formulation
a = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a += Constraint(u, gradu)[j, k] * del_lamu[j, k] * dx
a += Constraint(del_u, del_gradu)[j, k] * lamu[j, k] * dx
#.................................................Assembling
bilinear_form = form(a)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
# ################################### Apply Nullspace
def build_nullspace(V: Space):
bs = Space.dofmap.index_map_bs # number of block
length0 = Space.dofmap.index_map.size_local # number of dof in the cell
basis = [la.vector(Space.dofmap.index_map, bs=bs, dtype=PETSc.ScalarType) for i in range(18)] # 18 is the number of basis we need
b = [b.array for b in basis]
# Get dof indices for each subspace (x, y and z dofs)
dof_u = [Space.sub(0).sub(i).dofmap.list.flatten() for i in
range(3)] # creates dof associated with each nodal variables
# Set the three translational rigid body modes
for i in range(3):
b[i][dof_u[i]] = 1.0
# Set the three rotational rigid body modes
S0, S0_to_Space = Space.sub(0).collapse()
x = S0.tabulate_dof_coordinates()
dofs_block = S0.dofmap.list.flatten()
x0, x1, x2 = x[dofs_block, 0], x[dofs_block, 1], x[dofs_block, 2]
b[3][dof_u[0]] = -x1
b[3][dof_u[1]] = x0
b[4][dof_u[0]] = x2
b[4][dof_u[2]] = -x0
b[5][dof_u[2]] = x1
b[5][dof_u[1]] = -x2
##################################
dof_gu = [Space.sub(1).sub(i).dofmap.list.flatten() for i in range(9)]
b[6][dof_gu[0]] = 1.0
b[7][dof_gu[4]] = 1.0
b[8][dof_gu[8]] = 1.0
b[9][dof_gu[5]] = 1.0
b[9][dof_gu[7]] = -1.0
b[10][dof_gu[2]] = 1.0
b[10][dof_gu[6]] = -1.0
b[11][dof_gu[1]] = 1.0
b[11][dof_gu[3]] = -1.0
#######################
dof_lam = [Space.sub(2).sub(i).dofmap.list.flatten() for i in
range(9)] # creates dof associated with each nodal variables
b[12][dof_lam[0]] = 1.0
b[13][dof_lam[1]] = 1.0
b[13][dof_lam[3]] = 1.0
b[14][dof_lam[2]] = 1.0
b[14][dof_lam[6]] = 1.0
b[15][dof_lam[4]] = 1.0
b[16][dof_lam[5]] = 1.0
b[16][dof_lam[7]] = 1.0
b[17][dof_lam[8]] = 1.0
_basis = [x._cpp_object for x in basis]
dolfinx.cpp.la.orthonormalize(_basis)
assert dolfinx.cpp.la.is_orthonormal(_basis)
basis_petsc = [
PETSc.Vec().createWithArray(x[: bs * length0], bsize=1, comm=Space.mesh.comm) # type: ignore
for x in b
]
return PETSc.NullSpace().create(vectors=basis_petsc)
nullspace = build_nullspace(Space)
assert nullspace.test(A)
β
β