Hello everyone,
I am trying to do mechanical analysis of a plate with a hole, as shown in figure.
A parametric study for the y-displacement at the midpoint of top surface as a function of the length of the geometry is carried out. The bottom face is fixed, while the top is subjected to traction boundary condition (Tr = 1). The sides have a periodic boundary condition. Height of the block is D ( = 1) and width is W( = 1). D is the parameter to be varied for the study.
When I do the parametric study using a for loop in FEniCSx the first iteration gives correct result and the iterations thereafter give wrong results. When I do the same parametric study manually i.e. run the simulation for different parameter values sequentially, I get different values for y-displacement. For example,
D = 1.00, B[1] (manually) = 1.10308, B[1] (for-loop) = 1.10308
D = 1.10, B[1] (manually) = 1.16850, B[1] (for-loop) = 0.98714
Probably, there is some issue with my implementation of the multi-point constraint function within the for-loop. I am using FEniCSx 0.9 on Ubuntu. Can someone please help me with the issue. A MWE is attached below
import gmsh
import os
import numpy as np
import pyvista
import dolfinx
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element
from dolfinx import plot, default_scalar_type, mesh, geometry
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, functionspace, Expression, assemble_scalar, dirichletbc, form, locate_dofs_topological, locate_dofs_geometrical, set_bc, petsc, Function)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc, LinearProblem)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import (create_mesh, meshtags_from_entities, locate_entities_boundary)
from ufl import (FacetNormal, Identity, Measure, TestFunction, lhs, TrialFunction, as_vector, div, dot, ds, dx, inner, grad, nabla_grad, rhs, sym, system)
import ufl
from dolfinx import io
import dolfinx_mpc
from dolfinx_mpc import LinearProblem, MultiPointConstraint
from dolfinx_mpc.utils import create_point_to_point_constraint
from dolfinx.common import Timer, TimingType, list_timings
os.system('clear')
print(f"DOLFINx version: {dolfinx.__version__}")
gmsh.initialize()
#==================================================================
no_iter = 2
B_D = np.zeros((no_iter, 2))
for iter in range(no_iter):
W = 1.00
D = 1.00 +iter*0.1
r = 0.20
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, W, D)
circle = gmsh.model.occ.addDisk(W/2, D/3, 0, r, r)
bool1 = gmsh.model.occ.cut([(gdim, rectangle)],[(gdim, circle)])[0]
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim, [bool1[0][1]], 1*(iter+1))
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
gmsh.option.setNumber("Mesh.Algorithm", 11)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 1)
gmsh.option.setNumber("Mesh.RecombineAll", 0)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2)
gmsh.model.mesh.generate(gdim)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
V = functionspace(domain, ("Lagrange", 1, (gdim,)))
u_sol = Function(V, name="Displacement")
E = 1.0
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = E*nu / ((1+nu)*(1-2*nu))
boundaries = [(1, lambda x: np.isclose(x[1], 0)),
(2, lambda x: np.isclose(x[1], D))]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = dolfinx.mesh.locate_entities(domain, 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 = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
def on_bottom(x):
return np.isclose(x[1], 0)
bottom_dofs = locate_dofs_geometrical(V, on_bottom)
bc_bottom = dirichletbc(np.zeros((2,)), bottom_dofs, V)
bcs = [bc_bottom]
def on_side(x):
return np.isclose(x[0], 1)
def eps(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lmbda * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*eps(u)
##Periodic boundary condition =============================
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = 1 - x[0]
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
with Timer("~Elasticity: Initialize MPC"):
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, on_side, periodic_relation, bcs)
mpc.finalize()
##Variational problem =====================================
ds = ufl.Measure("ds", domain = domain, subdomain_data = facet_tag)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
Tr = Constant(domain, default_scalar_type((0, 1.00)))
a = ufl.inner(sigma(u), eps(v)) * ufl.dx
L = ufl.dot(Tr, v) * ds(2)
problem = LinearProblem(a, L, mpc, bcs= bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#===========================================================
point = np.zeros((3,1))
point[0] = W/2
point[1] = D
tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(tree, point.T)
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, point.T)
for i, pt in enumerate(point.T):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(pt)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
B = uh.eval(points_on_proc, cells)
B_D[iter] = [D, B[1]]
print(B_D)