Problems with Periodic BCs in 3D

General
1

Dear all,

I try to implement a 3D problem with periodic boundary conditions in x-,y- and z-direction. The mesh consists of a box with a sphere inside and it is generated with gmsh. The physical problem is the linear thermal heat equation.

My procedure to add periodic constraint is the following:

  • At first, I try to map the slave surfaces to the opposite master surfaces without the edges where the particular slave surface is connected to the other two slave surfaces.
  • Second, I try to map the point (1,1,1) [corner] to the point (0,0,0). That means that the connecting point of all slave surfaces is mapped to the connecting point of all master surfaces.
  • Third, I try to map the edges of slave-slave intersections to the opposite master-master intersections, i.e. (x,1,1) → (x,0,0), (1,y,1) → (0,y,0), (1,1,z) → (0,0,z).

Is the procedure in principle wrong? If you do not have time to look into the attached example, is there any reference of the procedure that can be used in dolfinx_mpc to generate 3D periodic boundary coditions with the slave-master concept in the literature? Or is there any fully periodic example implemented in dolfinx that I am missing?

The following problems arise when examining the results:

  • One corner of the box is not periodic.
  • One edge is not periodic.
  • At one surface, some nodes seems to be not periodic.

I have no idea what the problem is. I would really appreciate any help! Thank you!

I tried to reduce it to a MWE. However, the implementation of periodic BCs is rather long. Here is the example:

import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh

from ufl import grad, dot, lhs, rhs, Measure, TrialFunction, TestFunction, FiniteElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

##################################################
## Mesh generation Start
##################################################
gmsh.initialize()
r = 0.1     # radius of sphere
L = 1.0     # length of box

gdim = 3         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if rank == 0:
    # Define geometry for RVE
    sphere = gmsh.model.occ.add_sphere(L/2, L/2, L/2, r)
    box = gmsh.model.occ.add_box(0.0, 0.0, 0.0, L, L, L)
    whole_domain = gmsh.model.occ.fragment([(3, box)], [(3, sphere)])
    gmsh.model.occ.synchronize()

    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [box], tag=1)         # tag for matrix
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [sphere], tag=2)      # tag for inclusion

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (1, 1), tag = 9)

    background_surfaces = []
    inclusion_surfaces = []
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        mass = gmsh.model.occ.getMass(domain[0], domain[1])
        if np.isclose(mass, 4/3 * np.pi * (r ** 3)):  # identify inclusion by its mass
            inclusion_surfaces.append(domain)
        else:
            background_surfaces.append(domain)

    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "CurvesList", [e[0] for e in edges])
    gmsh.model.mesh.field.setNumber(1, "Sampling", 300)

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.010)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.030)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.00)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.075)
    gmsh.model.mesh.field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.model.mesh.generate(gdim)

domain, ct, _ = model_to_mesh(gmsh.model, comm, rank, gdim=gdim)
gmsh.finalize()

dx = Measure('dx', domain=domain, subdomain_data=ct)
x = ufl.SpatialCoordinate(domain)
##################################################
## Mesh generation finished
##################################################

# elements and funcionspace
###########################

P1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = FunctionSpace(domain, P1)

# define function, trial function and test function
T_fluc = TrialFunction(V)
dT_fluc = TestFunction(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
kappa = Function(S)
matrix = ct.find(1)
kappa.x.array[matrix] = np.full_like(matrix, 5, dtype=PETSc.ScalarType)
inclusion = ct.find(2)
kappa.x.array[inclusion]  = np.full_like(inclusion, 1, dtype=PETSc.ScalarType)
kappa.x.scatter_forward()

# boundary conditions
#####################
bcs = []
## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - L
        return out_x

    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], L)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 3)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):

    # slave/master mapping of opposite surfaces (without slave-slave[-slave] intersections)
    def pbc_slave_to_master_map(x):
        out_x = pbc_slave_to_master_maps[i](x)
        # remove edges that are connected to another slave surface
        idx = np.logical_or(pbc_is_slave[(i + 1) % N_pbc](x),pbc_is_slave[(i + 2) % N_pbc](x))
        out_x[pbc_directions[i]][idx] = np.nan
        return out_x

    mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    def pbc_slave_to_master_map_corner(x):
        '''
        Maps the slave corner dof [intersection of slave_x, slave_y, slave_z] (1,1,1) to
        master corner dof [master_x, master_y, master_z] (0,0,0)
        '''
        out_x = x.copy()
        out_x[0] = x[0] - L
        out_x[1] = x[1] - L
        out_x[2] = x[2] - L
        idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
        out_x[0][~idx_corner] = np.nan
        out_x[1][~idx_corner] = np.nan
        out_x[2][~idx_corner] = np.nan
        return out_x

    def generate_pbc_slave_to_master_map_edges(dir_i, dir_j):
        def pbc_slave_to_master_map_edges(x):
            '''
            Maps the slave edge dofs [intersection of slave_i, slave_j] (i=1,j=1) to
            master corner dof [master_i, master_j] (i=0,j=0)
            (1) map slave_x, slave_y to master_x, master_y: i=0, j=1
            (2) map slave_x, slave_z to master_x, master_z: i=0, j=2
            (3) map slave_y, slave_z to master_y, master_z: i=1, j=2
            '''
            out_x = x.copy()
            out_x[dir_i] = x[dir_i] - L
            out_x[dir_j] = x[dir_j] - L
            idx = np.logical_and(pbc_is_slave[dir_i](x), pbc_is_slave[dir_j](x))
            out_x[dir_i][~idx] = np.nan
            out_x[dir_j][~idx] = np.nan
            idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
            out_x[dir_i][idx_corner] = np.nan
            out_x[dir_j][idx_corner] = np.nan
            print(len(out_x[dir_i][idx]))
            return out_x
        return pbc_slave_to_master_map_edges

    mapping_slave_intersections = [(0,1),(0,2),(1,2)] # pairs of slave intersections

    for ij in range(gdim):
        # mapping slave intersection node (corner) to master intersection node (opposite corner)
        mpc.create_periodic_constraint_topological(V, pbc_meshtags[0],
                                                   pbc_slave_tags[0],
                                                   pbc_slave_to_master_map_corner,
                                                   bcs)

        for inters in mapping_slave_intersections:
            # mapping slave intersections to opposite master intersections
            mpc.create_periodic_constraint_topological(V, pbc_meshtags[inters[0]],
                                                       pbc_slave_tags[inters[0]],
                                                       generate_pbc_slave_to_master_map_edges(inters[0], inters[1]),
                                                       bcs)

mpc.finalize()

# Variational problem
#####################
T_gradient = Constant(domain, PETSc.ScalarType((1,0,0)))

## weak formulation
eqn = - kappa * dot(grad(dot(T_gradient , x) + T_fluc), grad(dT_fluc)) * dx
a_form = lhs(eqn)
L_form = rhs(eqn)

## solving
func_sol = Function(V)
problem = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
func_sol = problem.solve()

# write results to file
file_results1 = XDMFFile(domain.comm, "results_3D_MWE/T_fluc.xdmf", "w")
file_results1.write_mesh(domain)
file_results1.write_function(func_sol)
file_results1.close()
2

You have the procedure correct, in principle. The following demo (a bit complicated due to its generality, I’m afraid) solves a 2D problem with fully periodic boundaries: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py

In particular, the following lines map the single dof at the upper-right corner of the domain to the lower-left corner: https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_periodic_gep.py#L294-L308.

For 3D, I think you will need to:

  • map the surfaces, excluding all dofs on the edges
  • map the edges, excluding all dofs on the vertices
  • map the vertices
3

Specifically, the mappings you’ll need are these (I think). Other groupings are possible, e.g. you could include endpoints in the edge mappings with appropriate care.

periodic bc
periodic bc2556×2112 488 KB

Dashed lines indicate that the boundary of a surface is not included in the mapping, and open circles indicate that the endpoint of an edge is not included in the mapping.

4 Likes
4

Dear Connor,
thank you very much for your effort and your explanations! That is exactly what I thought, but it is very good to know that the procedure is basically correct.
However, there must be a bug in my code so I have to go through and have a look what is going wrong.

All the best!

5

For interested users:
After some considerations of @conpierce8, the problem was related to the mesh as there were no periodicity restrictions during the mesh generation. Even though dolfinx_mpc is able to handle periodic conditions on non-periodic meshes, you need to be careful that no slave node is used for a mapping as a master node (especially in the environment of edges).

Here is the example with fully periodic boundaries on a periodic mesh:

import gmsh
import dolfinx
import dolfinx_mpc
from dolfinx.fem import FunctionSpace, Function, Constant
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io.gmshio import model_to_mesh

from ufl import grad, dot, lhs, rhs, Measure, TrialFunction, TestFunction, FiniteElement
import ufl
from petsc4py import PETSc
import numpy as np
from mpi4py import MPI

# MPI initialization
comm = MPI.COMM_WORLD
rank = MPI.COMM_WORLD.rank

##################################################
## Mesh generation Start
##################################################
gmsh.initialize()
r = 0.1     # radius of sphere
L = 1.0     # length of box

gdim = 3         # Geometric dimension of the mesh
fdim = gdim - 1  # facet dimension

if rank == 0:
    # Define geometry for RVE with periodicity of the surfaces
    sphere = gmsh.model.occ.add_sphere(L / 2, L / 2, L / 2, r)
    box = gmsh.model.occ.add_box(0.0, 0.0, 0.0, L, L, L)
    whole_domain = gmsh.model.occ.fragment([(3, box)], [(3, sphere)])
    gmsh.model.occ.synchronize()

    # Ask OpenCASCADE to compute more accurate bounding boxes of entities using
    # the STL mesh:
    gmsh.option.setNumber("Geometry.OCCBoundsUseStl", 1)

    # The periodicity transform is provided as a 4x4 affine transformation matrix given by row.
    translation_x = [1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1]
    translation_y = [1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1]
    translation_z = [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1]

    eps = 1e-3
    # We get all surfaces on the left:
    sxmin = gmsh.model.getEntitiesInBoundingBox(- eps, -eps, -eps, eps, L + eps, L + eps, 2)
    # We get all surfaces on the bottom:
    symin = gmsh.model.getEntitiesInBoundingBox(- eps, -eps, -eps, L + eps, eps, L + eps, 2)
    # We get all surfaces on the front:
    szmin = gmsh.model.getEntitiesInBoundingBox(- eps, -eps, -eps, L + eps, L + eps, eps, 2)

    for i_x in sxmin:
        # Then we get the bounding box of each left surface
        xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(i_x[0], i_x[1])
        # We translate the bounding box to the right and look for surfaces inside it:
        sxmax = gmsh.model.getEntitiesInBoundingBox(xmin - eps + L, ymin - eps,
                                                    zmin - eps, xmax + eps + L,
                                                    ymax + eps, zmax + eps, 2)
        # For all the matches, we compare the corresponding bounding boxes...
        for j_x in sxmax:
            xmin2, ymin2, zmin2, xmax2, ymax2, zmax2 = gmsh.model.getBoundingBox(j_x[0], j_x[1])
            xmin2 -= 1
            xmax2 -= 1
            # ...and if they match, we apply the periodicity constraint
            if (abs(xmin2 - xmin) < eps and abs(xmax2 - xmax) < eps
                    and abs(ymin2 - ymin) < eps and abs(ymax2 - ymax) < eps
                    and abs(zmin2 - zmin) < eps and abs(zmax2 - zmax) < eps):
                gmsh.model.mesh.setPeriodic(2, [j_x[1]], [i_x[1]], translation_x)

    for i_y in symin:
        # Then we get the bounding box of each left surface
        xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(i_y[0], i_y[1])
        # We translate the bounding box to the right and look for surfaces inside it:
        symax = gmsh.model.getEntitiesInBoundingBox(xmin - eps, ymin - eps + L,
                                                    zmin - eps, xmax + eps,
                                                    ymax + eps + L, zmax + eps, 2)
        # For all the matches, we compare the corresponding bounding boxes...
        for j_y in symax:
            xmin2, ymin2, zmin2, xmax2, ymax2, zmax2 = gmsh.model.getBoundingBox(j_y[0], j_y[1])
            xmin2 -= 1
            xmax2 -= 1
            # ...and if they match, we apply the periodicity constraint
            if (abs(xmin2 - xmin) < eps and abs(xmax2 - xmax) < eps
                    and abs(ymin2 - ymin) < eps and abs(ymax2 - ymax) < eps
                    and abs(zmin2 - zmin) < eps and abs(zmax2 - zmax) < eps):
                gmsh.model.mesh.setPeriodic(2, [j_y[1]], [i_y[1]], translation_y)


    for i_z in szmin:
        # Then we get the bounding box of each left surface
        xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(i_z[0], i_z[1])
        # We translate the bounding box to the right and look for surfaces inside it:
        szmax = gmsh.model.getEntitiesInBoundingBox(xmin - eps, ymin - eps,
                                                    zmin - eps + L, xmax + eps,
                                                    ymax + eps, zmax + eps + L, 2)
        # For all the matches, we compare the corresponding bounding boxes...
        for j_z in szmax:
            xmin2, ymin2, zmin2, xmax2, ymax2, zmax2 = gmsh.model.getBoundingBox(j_z[0], j_z[1])
            xmin2 -= 1
            xmax2 -= 1
            # ...and if they match, we apply the periodicity constraint
            if (abs(xmin2 - xmin) < eps and abs(xmax2 - xmax) < eps
                    and abs(ymin2 - ymin) < eps and abs(ymax2 - ymax) < eps
                    and abs(zmin2 - zmin) < eps and abs(zmax2 - zmax) < eps):
                gmsh.model.mesh.setPeriodic(2, [j_z[1]], [i_z[1]], translation_z)


    matrix_physical = gmsh.model.addPhysicalGroup(gdim, [box], tag=1)         # tag for matrix
    inclusion_physical = gmsh.model.addPhysicalGroup(gdim, [sphere], tag=2)      # tag for inclusion

    # Get the interface elements between rectangle and circle and tag
    interface_elements = gmsh.model.getBoundary([(gdim, inclusion_physical)])
    interface_physical = gmsh.model.addPhysicalGroup(fdim, (1, 1), tag = 9)

    background_surfaces = []
    inclusion_surfaces = []
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        mass = gmsh.model.occ.getMass(domain[0], domain[1])
        if np.isclose(mass, 4/3 * np.pi * (r ** 3)):  # identify inclusion by its mass
            inclusion_surfaces.append(domain)
        else:
            background_surfaces.append(domain)

    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(inclusion_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "CurvesList", [e[0] for e in edges])
    gmsh.model.mesh.field.setNumber(1, "Sampling", 300)

    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", 0.010)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 0.030)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 0.00)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 0.075)
    gmsh.model.mesh.field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.model.mesh.generate(gdim)

domain, ct, _ = model_to_mesh(gmsh.model, comm, rank, gdim=gdim)
gmsh.finalize()

xdmf = dolfinx.io.XDMFFile(domain.comm, "results_3D_MWE_v4/mesh.xdmf", "w")
xdmf.write_mesh(domain)
xdmf.close()
##################################################
## Mesh generation finished
##################################################

dx = Measure('dx', domain=domain, subdomain_data=ct)
x = ufl.SpatialCoordinate(domain)

# elements and funcionspace
###########################

P1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = FunctionSpace(domain, P1)

# define function, trial function and test function
T_fluc = TrialFunction(V)
dT_fluc = TestFunction(V)

# define material parameter
S = FunctionSpace(domain, ("DG", 0))
kappa = Function(S)
matrix = ct.find(1)
kappa.x.array[matrix] = np.full_like(matrix, 5, dtype=PETSc.ScalarType)
inclusion = ct.find(2)
kappa.x.array[inclusion]  = np.full_like(inclusion, 1, dtype=PETSc.ScalarType)
kappa.x.scatter_forward()

# boundary conditions
#####################
bcs = []
## periodic boundary conditions
pbc_directions = []
pbc_slave_tags = []
pbc_is_slave = []
pbc_is_master = []
pbc_meshtags = []
pbc_slave_to_master_maps = []

mpc = dolfinx_mpc.MultiPointConstraint(V)

def generate_pbc_slave_to_master_map(i):
    def pbc_slave_to_master_map(x):
        out_x = x.copy()
        out_x[i] = x[i] - L
        return out_x

    return pbc_slave_to_master_map

def generate_pbc_is_slave(i):
    return lambda x: np.isclose(x[i], L)

def generate_pbc_is_master(i):
    return lambda x: np.isclose(x[i], 0.0)

for i in range(gdim):

    pbc_directions.append(i)
    pbc_slave_tags.append(i + 3)
    pbc_is_slave.append(generate_pbc_is_slave(i))
    pbc_is_master.append(generate_pbc_is_master(i))
    pbc_slave_to_master_maps.append(generate_pbc_slave_to_master_map(i))

    facets = locate_entities_boundary(domain, fdim, pbc_is_slave[-1])
    arg_sort = np.argsort(facets)
    pbc_meshtags.append(meshtags(domain,
                                 fdim,
                                 facets[arg_sort],
                                 np.full(len(facets), pbc_slave_tags[-1], dtype=np.int32)))

N_pbc = len(pbc_directions)
for i in range(N_pbc):

    # slave/master mapping of opposite surfaces (without slave-slave[-slave] intersections)
    def pbc_slave_to_master_map(x):
        out_x = pbc_slave_to_master_maps[i](x)
        # remove edges that are connected to another slave surface
        idx = np.logical_or(pbc_is_slave[(i + 1) % N_pbc](x),pbc_is_slave[(i + 2) % N_pbc](x))
        out_x[i][idx] = np.nan
        print("Number of Nodes for Surface mapping that are excluded: ", len(out_x[i][idx]))
        print("Number of Nodes for Surface mapping: ", len(out_x[i][~idx]))
        return out_x

    mpc.create_periodic_constraint_topological(V, pbc_meshtags[i],
                                                   pbc_slave_tags[i],
                                                   pbc_slave_to_master_map,
                                                   bcs)

if len(pbc_directions) > 1:
    def generate_pbc_slave_to_master_map_edges(dir_i, dir_j):
        def pbc_slave_to_master_map_edges(x):
            '''
            Maps the slave edge dofs [intersection of slave_i, slave_j] (i=1,j=1) to
            master corner dof [master_i, master_j] (i=0,j=0)
            (1) map slave_x, slave_y to master_x, master_y: i=0, j=1
            (2) map slave_x, slave_z to master_x, master_z: i=0, j=2
            (3) map slave_y, slave_z to master_y, master_z: i=1, j=2
            '''
            out_x = x.copy()
            out_x[dir_i] = x[dir_i] - L
            out_x[dir_j] = x[dir_j] - L
            idx = np.logical_and(pbc_is_slave[dir_i](x), pbc_is_slave[dir_j](x))
            out_x[dir_i][~idx] = np.nan
            out_x[dir_j][~idx] = np.nan
            # remove corner DOFs at point (1,1,1)
            idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
            out_x[dir_i][idx_corner] = np.nan
            out_x[dir_j][idx_corner] = np.nan
            print("Number of Nodes for edge mapping: ", len(out_x[dir_i][idx]), "Minus", len(out_x[dir_i][idx_corner]), "=", len(out_x[dir_i][idx])-len(out_x[dir_i][idx_corner]))
            return out_x
        return pbc_slave_to_master_map_edges

    def pbc_slave_to_master_map_corner(x):
        '''
        Maps the slave corner dof [intersection of slave_x, slave_y, slave_z] (1,1,1) to
        master corner dof [master_x, master_y, master_z] (0,0,0)
        '''
        out_x = x.copy()
        out_x[0] = x[0] - L
        out_x[1] = x[1] - L
        out_x[2] = x[2] - L
        idx_corner = np.logical_and(pbc_is_slave[0](x), np.logical_and(pbc_is_slave[1](x), pbc_is_slave[2](x)))
        out_x[0][~idx_corner] = np.nan
        out_x[1][~idx_corner] = np.nan
        out_x[2][~idx_corner] = np.nan
        print("Number of Nodes of corner mapping:" , len(out_x[0][idx_corner]))
        return out_x

    mapping_slave_intersections = [(0,1),(0,2),(1,2)] # pairs of slave intersections

    # mapping slave intersection node (corner) to master intersection node (opposite corner)
    mpc.create_periodic_constraint_topological(V, pbc_meshtags[0],
                                               pbc_slave_tags[0],
                                               pbc_slave_to_master_map_corner,
                                               bcs)

    for inters in mapping_slave_intersections:
        # mapping slave intersections to opposite master intersections
        mpc.create_periodic_constraint_topological(V, pbc_meshtags[inters[0]],
                                                   pbc_slave_tags[inters[0]],
                                                   generate_pbc_slave_to_master_map_edges(inters[0], inters[1]),
                                                   bcs)

mpc.finalize()

# Variational problem
#####################
T_gradient = Constant(domain, PETSc.ScalarType((1,0,0)))

## weak formulation
eqn = - kappa * dot(grad(dot(T_gradient , x) + T_fluc), grad(dT_fluc)) * dx
a_form = lhs(eqn)
L_form = rhs(eqn)

## solving
func_sol = Function(V)
problem = dolfinx_mpc.LinearProblem(a_form, L_form, mpc, bcs=bcs)
func_sol = problem.solve()

# write results to file
file_results1 = XDMFFile(domain.comm, "results_3D_MWE/T_fluc.xdmf", "w")
file_results1.write_mesh(domain)
file_results1.write_function(func_sol)
file_results1.close()
6

how would you add other boundary conditions to this? I tried to remove the z-periodic boundary condition, but the code doesn’t seem to accept the Neumann or dirichlet conditions?

7

No, this example is explicit for fully periodic boundary conditions in 3D, as we create periodic constraints in each direction ( for i in range(gdim):). However, you can modify the code so that you only have periodic constraints in two directions. Please note that you must specify the boundary conditons bcs before creating periodic BCs (no more empty list). Neumann BCs are inserted via the variational formulation.

1 Like
8

Excellent advice. I have done all of that; however, there is a clear error in the results. I’m not sure what is causing it, I think it is from the Dirichlet boundary. I am going to keep digging and get back with a more complete question. That is very interesting about defining the BC before the periodic boundaries.

9

I have isolated the problem… I am having an issue where the periodic boundary conditions and the Neumann boundary conditions intersect. I assume I am going to have to go in and remove the nodes contained in the Neumann boundary from the periodic boundary? Is it possible to be in both?

10

Yes. In my opinion, you either need to remove the DOFs from the slave/master node set or define ds so that the nodes do not intersect with the periodic BCs.

11

To my surprise, that did work either. I have figured it out though! I figured out that the periodic boundaries were also getting treated as Neumann conditions. Turns out I needed to do two things.

  1. At the start when you are defining the Dirichlet (to feed into the periodic boundary condition creation), you need to make a ‘dummy’ Dirichlet BC on the boundary you intend to put the Neumann conditions, then make sure this dummy condition doesn’t go into the solve.

  2. Define the ds for the Neumann condition (as you just mentioned)

Thanks for your advice on this.

Regards
Jess