Dear all,
**I edited this post because I thought the previous version was too long and confusing.
I am currently working on a linear elasticity problem. In addition to displacement, I want to introduce a new variable, the displacement gradient, to my system, and I want to set up the null space for my formulation.
BCs:
Zero displacement at the bottom surface, constant pressure to the top surface of cube.
I followed the post here, and tried to set up the null space but I got an error.
Here is MWE:
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunctions, TrialFunctions,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split, action,
SpatialCoordinate)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import time
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
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = 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])],
[10, 10, 10], cell_type=CellType.hexahedron)
# ......................................... Define Loading
Force = 200.00000 / Fs
Area2 = h ** 2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
# ......................................... 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)
# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], H, atol=tol))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# ....................................... BCs
S0, S0_to_Space = Space.sub(0).collapse()
# Clamp BC for left side
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((la * StrainT(u)[i, i] * delta[j,k] + 2 * mu * 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
F1 = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
FC1 = Constraint(u, gradu)[j, k] * del_lamu[j, k] * dx
FC2 = Constraint(del_u, del_gradu)[j, k] * lamu[j, k] * dx
# ......................................... Preparing the form
a = F1 + FC1 + FC2
L = tr_S2[i] * del_u[i] * ds(2)
#.................................................Assembling
Result = Function(Space)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = create_vector(linear_form)
with b.localForm() as b_loc:
b_loc.set(0)
assemble_vector(b, linear_form)
apply_lifting(b, [bilinear_form], [bc])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) # if run in parrallel
set_bc(b, bc)
b.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
################################### Apply Nullspace
DG0, DG0_to_Space = Space.sub(1).collapse()
LA0, LA0_to_Space = Space.sub(2).collapse()
ns_vec = dolfinx.fem.Function(Space)
ns_vec.x.array[DG0_to_Space] = 1
ns_vec.x.array[LA0_to_Space] = 1
ns_vec2 = dolfinx.fem.Function(Space)
ns_vec2.x.array[LA0_to_Space] = 1
ns_vec.vector.normalize()
ns_vec2.vector.normalize()
nullspace = PETSc.NullSpace().create(vectors = [ns_vec.vector])
assert nullspace.test(A)
A.setNearNullSpace(nullspace)
nullspace.remove(b)