3D Periodic Boundary Conditions dofmap error- MixedElement

Errors
1

Hi FEniCS family,
I am getting error in creating function space using periodic boundary conditions. Please find the below MWE.

from dolfin import *
import numpy as np
a1 = 1
## Create mesh and define function space
mesh = BoxMesh(Point(-0.5, -0.5, -0.5), Point(0.5,0.5,0.5), a1, a1, a1)


class PeriodicBoundary(SubDomain):
    
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two slave edges
        return bool ((near(x[0], -a1/2) or near(x[1], -a1/2) or near(x[2], -a1/2)) and 
            (not ((near(x[0], a1/2) and near(x[2], a1/2)) or 
                  (near(x[0], a1/2) and near(x[1], a1/2)) or
                  (near(x[1], a1/2) and near(x[2], a1/2)))) and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
    	#### define mapping for a single point in the box, such that 3 mappings are required
        if near(x[0], a1/2) and near(x[1], a1/2) and near(x[2],a1/2):
            y[0] = x[0] - a1
            y[1] = x[1] - a1
            y[2] = x[2] - a1
        ##### define mapping for edges in the box, such that mapping in 2 Cartesian coordinates are required
        elif near(x[0], a1/2) and near(x[2], a1/2):
            y[0] = x[0] - a1
            y[1] = x[1] 
            y[2] = x[2] - a1      
        elif near(x[1], a1/2) and near(x[2], a1/2):
            y[0] = x[0] 
            y[1] = x[1] - a1
            y[2] = x[2] - a1
        elif near(x[0], a1/2) and near(x[1], a1/2):
            y[0] = x[0] - a1
            y[1] = x[1] - a1
            y[2] = x[2]         
        #### right maps to left: left/right is defined as the x-direction
        elif near(x[0], a1/2):
            y[0] = x[0] - a1
            y[1] = x[1]
            y[2] = x[2]
        ### back maps to front: front/back is defined as the y-direction    
        elif near(x[1], a1/2):
            y[0] = x[0]
            y[1] = x[1] - a1
            y[2] = x[2] 
        #### top maps to bottom: top/bottom is defined as the z-direction        
        elif near(x[2], a1/2):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - a1

Ve = VectorElement("CG", mesh.ufl_cell(), 1,dim=3)
Re = VectorElement("R", mesh.ufl_cell(), 0,dim=3)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary())

I am getting the following error

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_166/3967937147.py in <module>
      1 Ve = VectorElement("CG", mesh.ufl_cell(), 1,dim=3)
      2 Re = VectorElement("R", mesh.ufl_cell(), 0,dim=3)
----> 3 W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary())

/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/functionspace.py in __init__(self, *args, **kwargs)
     29                 pass
     30             elif len(args) == 2:
---> 31                 self._init_from_ufl(*args, **kwargs)
     32             else:
     33                 self._init_convenience(*args, **kwargs)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/functionspace.py in _init_from_ufl(self, mesh, element, constrained_domain)
     50             dolfin_dofmap = cpp.fem.DofMap(ufc_dofmap, mesh)
     51         else:
---> 52             dolfin_dofmap = cpp.fem.DofMap(ufc_dofmap, mesh,
     53                                            constrained_domain)
     54 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to periodic boundary mapping.
*** Reason:  Need to set coordinate 0 in sub_domain.map.
*** Where:   This error was encountered inside PeriodicBoundaryComputation.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

The error,I guess, is due to the MixedElement in the FunctionSpace. When I changed my FunctionSpace to
V = FunctionSpace(mesh, 'Lagrange', 3, constrained_domain=PeriodicBoundary())

I got no error. In the example of periodic homogenizationof FEniCS, there was no error in using MixedElement. I have only changed the 2D terms to 3D using the following resources.

Can you help me to remove this error?
Any help is appreciated.

2

Hi, thanks for posting this problem. I have a similar problem on legacy fences (2019.1.0), where the mixed function space trigger an error. It seems like this problem is unique to the 3D periodic boundary. Have you found a fix?

3

Try this. It worked for me.


class PeriodicBoundary(SubDomain):
        def __init__(self, vertices, tolerance=DOLFIN_EPS):
            SubDomain.__init__(self, tolerance)
            self.tol = tolerance
            self.vv = vertices
            self.a1 = self.vv[1,:]-self.vv[4,:] # first vector generating periodicity
            self.a2 = self.vv[2,:]-self.vv[5,:] # second vector generating periodicity 
            self.a3 = self.vv[3,:]-self.vv[6,:] # third vector generating periodicity       
        def inside(self, x, on_boundary):     

            """return True if on left, bottom, or back faces & 
            not on one of the top, front or right faces or related edges/ vertices)"""

            # faces
            left = near(x[0],self.vv[4,0]) 
            bottom = near(x[1],self.vv[5,1]) 
            back = near(x[2],self.vv[6,2])
            right = near(x[0],self.vv[1,0])
            top = near(x[1],self.vv[2,1])
            front = near(x[2],self.vv[3,2])

            # line-segments (bottom 4; top 4; vertical 4) edges
            bottom_front = bottom and front 
            bottom_right = bottom and right 

            top_left = top and left 
            top_back = top and back 

            left_front = left and front 
            right_back = right and back 

            return bool((left or back or bottom) and 
                        (not( (top_left) or (left_front) or (top_back) or (right_back) or (bottom_right) or (bottom_front))) and on_boundary)

        def map(self, x, y):
            """ Mapping the right boundary to left, front to back and top to bottom"""

            # faces
            right = near(x[0],self.vv[1,0])
            top = near(x[1],self.vv[2,1])
            front = near(x[2],self.vv[3,2])

            # line-segments 
            top_right = top and right 
            top_front = top and front 
            right_front = right and front 
            point_6 = right and front and top 

            if point_6:
                y[0] = x[0] - (self.a1[0] + self.a2[0] + self.a3[0])
                y[1] = x[1] - (self.a1[1] + self.a2[1] + self.a3[1])
                y[2] = x[2] - (self.a1[2] + self.a2[2] + self.a3[2])
            elif top_right:
                y[0] = x[0] - (self.a1[0] + self.a2[0])
                y[1] = x[1] - (self.a1[1] + self.a2[1])
                y[2] = x[2] - (self.a1[2] + self.a2[2])
            elif top_front:
                y[0] = x[0] - (self.a2[0] + self.a3[0])
                y[1] = x[1] - (self.a2[1] + self.a3[1])
                y[2] = x[2] - (self.a2[2] + self.a3[2])
            elif right_front: 
                y[0] = x[0] - (self.a1[0] + self.a3[0])
                y[1] = x[1] - (self.a1[1] + self.a3[1])
                y[2] = x[2] - (self.a1[2] + self.a3[2])
            elif right:
                y[0] = x[0] - (self.a1[0])
                y[1] = x[1] - (self.a1[1])
                y[2] = x[2] - (self.a1[2])
            elif front:
                y[0] = x[0] - (self.a3[0])
                y[1] = x[1] - (self.a3[1])
                y[2] = x[2] - (self.a3[2])
            elif top:
                y[0] = x[0] - (self.a2[0])
                y[1] = x[1] - (self.a2[1])
                y[2] = x[2] - (self.a2[2])
            else: 
                y[0] = -1. 
                y[1] = -1. 
                y[2] = -1. 


W = FunctionSpace(mesh, MixedElement([Ve,Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))