Hey! I am Itryng to implement periodic boundary conditions to commpute the mechanical properties of an elastically heterogeneous material. I am trying to undestant how to do this using dolphinx_mpc. I wrote a code with Dirichlet boundary conditions which works fine. I have a hard time modifying it to implement periodic boundary conditions in the 2d square in both directions. This is my try with the following error:
import numpy as np
import ufl
import gmsh
import pyvista
from mpi4py import MPI
from dolfinx import default_scalar_type, fem, io, mesh
from dolfinx.fem import Function
import dolfinx.fem.petsc
from petsc4py import PETSc
from dolfinx.io import gmshio
from dolfinx.plot import vtk_mesh
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
import dolfinx_mpc
rank = MPI.COMM_WORLD.rank
# Definition of parameters
E1 = 1400 # Young's modulus of the nanoparticle in Mpa
nu1 = 0.3 # Poisson's ratio of the nanoparticle
E2 = 10 # Young's modulus of the matrix in Mpa
nu2 = 0.1 # Poisson's ratio of the matrix
l = 10 # Length of the simulation box
R = 2 # Radius of the particle
gdim = 2 # Geometric dimension of the mesh
volume = l**gdim
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gmsh.initialize()
if mesh_comm.rank == model_rank:
# Define geometry for iron cylinder
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, l, l, tag=1)
circle = gmsh.model.occ.addDisk(l/2, l/2, 0, R, R, tag=2)
whole_domain = gmsh.model.occ.fragment([(gdim, rectangle)], [(gdim, circle)])
gmsh.model.occ.synchronize()
# Create physical markers for the different phases.
# We use the following markers:
# - Bulk: 0
# - Particle: 1
bulk_marker = 0
particle_marker = 1
particle_surfaces = whole_domain[0][0]
bulk_surfaces = whole_domain[0][1]
print("Whole domain: ", whole_domain)
print("Bulk: ", bulk_surfaces)
print("Particle: ", particle_surfaces)
# Add marker for the vacuum
gmsh.model.addPhysicalGroup(2, [bulk_surfaces[1]], tag=bulk_marker)
gmsh.model.addPhysicalGroup(2, [particle_surfaces[1]], tag=particle_marker)
distance_field = gmsh.model.mesh.field.add("Distance")
edges = gmsh.model.getBoundary([particle_surfaces], oriented=False)
print("Edges: ", edges)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", R/6)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", R/4)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", R/4)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", R/2)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
# Generate mesh
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")
gmsh.write("RVE_2D.msh")
gmsh.finalize()
# Load mesh using original code's setup
domain, ct, ft = gmshio.read_from_msh("RVE_2D.msh", mesh_comm, model_rank, gdim=2)
# Define material properties as DG functions
Q = fem.functionspace(domain, ("DG", 0))
E = Function(Q)
nu = Function(Q)
bulk_cells = ct.find(bulk_marker)
E.x.array[bulk_cells] = np.full_like(bulk_cells, E2, dtype=default_scalar_type)
nu.x.array[bulk_cells] = np.full_like(bulk_cells, nu2, dtype=default_scalar_type)
particle_cells = ct.find(particle_marker)
E.x.array[particle_cells] = np.full_like(particle_cells, E1, dtype=default_scalar_type)
nu.x.array[particle_cells] = np.full_like(particle_cells, nu1, dtype=default_scalar_type)
# Define Lame parameters
mu = E / 2 / (1 + nu)
lamda = nu * E / (1 - 2 * nu) / (1 + nu)
# Strain and stress functions
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u, Eps):
return lamda * ufl.tr(epsilon(u) + Eps) * ufl.Identity(gdim) + 2*mu*(epsilon(u) + Eps)
# Macro-strain
def macro_strain(i):
ϵ = np.zeros((3,), dtype = np.float64)
ϵ[i] = 1
ϵ[2:] /= 2
return ufl.as_tensor([[ϵ[0], ϵ[2]],
[ϵ[2], ϵ[1]]])
# Voigt notation conversion
def stress2Voigt(s):
return ufl.as_vector([s[0,0],
s[1,1],
s[0,1]])
def strain2Voigt(s):
return ufl.as_vector([s[0,0],
s[1,1],
2*s[0,1]])
# Set up function space and periodic constraints
V = fem.functionspace(domain, ("Lagrange", 2, (gdim, )))
u_sol = fem.Function(V, name="Displacement")
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Definition of boundary conditions.
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], l)
def bottom(x):
return np.isclose(x[1], 0)
def top(x):
return np.isclose(x[1], l)
left_dofs = fem.locate_dofs_geometrical(V, left)
right_dofs = fem.locate_dofs_geometrical(V, right)
bottom_dofs = fem.locate_dofs_geometrical(V, bottom)
top_dofs = fem.locate_dofs_geometrical(V, top)
# Apply periodic boundary conditions (minimal intrusion)
tol = 1e-8
def periodic_boundary(x):
return np.isclose(x[0], l, atol=tol) | np.isclose(x[1], l, atol=tol)
def periodic_relation(x):
out_x = np.zeros(x.shape)
out_x[0] = l - x[0]
out_x[1] = l - x[1]
return out_x
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, [], default_scalar_type(1))
mpc.finalize()
# Define elasticity problem
Eps = ufl.as_tensor(((0, 0), (0, 0)))
f = fem.Constant(domain, np.zeros((2,)))
dx = ufl.Measure("dx", domain=domain)
a = ufl.inner(sigma(u, Eps), epsilon(v)) * dx
L = ufl.inner(f, v) * dx
# Homogenization calculation
Chom = np.zeros((3, 3))
for (j, case) in enumerate(["Exx", "Eyy", "Exy"]):
Eps = macro_strain(j)
problem = LinearProblem(a, L, u=u_sol, mpc=mpc, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
Sigma = np.zeros((3,))
for k in range(3):
Sigma[k] = fem.assemble_scalar(fem.form(stress2Voigt(sigma(u_sol, Eps))[k] * dx)) / volume
Chom[j, :] = Sigma
# Calculate effective properties
lmbda_hom, mu_hom = Chom[0, 1], Chom[2, 2]
E_hom = mu_hom * (3 * lmbda_hom + 2 * mu_hom) / (lmbda_hom + mu_hom)
nu_hom = lmbda_hom / (lmbda_hom + mu_hom) / 2
print("Effective Young's modulus:", E_hom)
print("Effective Poisson ratio:", nu_hom)
Traceback (most recent call last):
File "/home/totis7/homogenization_2D_periodic.py", line 185, in <module>
problem = LinearProblem(a, L, u=u_sol, mpc=mpc, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
File "/usr/lib/python3/dist-packages/dolfinx_mpc/problem.py", line 99, in __init__
raise ValueError(
ValueError: ('The input function has to be in the function space in the multi-point constraint', 'i.e. u = dolfinx.fem.Function(mpc.function_space)')
Exception ignored in: <function LinearProblem.__del__ at 0x7f16ab804160>
Traceback (most recent call last):
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 831, in __del__
self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'