Periodic boundary conditions Homogenization

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'

Hi,
this demo with dolfinx_mpc will appear soon. Try with something like:

def periodic_relation_left_right(x):
    out_x = np.zeros(x.shape)
    out_x[0] = l - x[0]
    out_x[1] = x[1]
    out_x[2] = x[2]
    return out_x


def periodic_relation_bottom_top(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0]
    out_x[1] = l - x[1]
    out_x[2] = x[2]
    return out_x


mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, right, periodic_relation_left_right, [])
mpc.create_periodic_constraint_geometrical(V, top, periodic_relation_bottom_top, [])
mpc.finalize()

Also you should use the LinearProblem from dolfinx_mpc

2 Likes

Thank you! It seems that it worked!

Demo is now available here: Periodic homogenization of linear elasticity — Computational Mechanics Numerical Tours with FEniCSx

1 Like

Great work! Thank ypu again for the feedback!