Hi,
I am currently facing some problems with complex meshes and periodic boundary conditions in fenics 2019.1: as long as I am doing 2D simulations, my implementation of periodic boundary conditions works fine - when I switch to 3D, I always find points on the boundary of my domains for which the periodicity mapping did not work.
One strange thing I noticed is that - during the checkup whether a point is inside the periodic subdomain or how the mapping works - the points that are passed to these methods are totally different from the ones that I can find by, e.g. checking the function spaces dof coordinates. The ones that are passed to the inside() and map() methods are really on the domain boundaries - if I check the dof coordinates, the number of coordinates on the boundary decreases. Does anybody know where that comes from? I am using simple Langrangian elements!
For the minimal code example below, I generated a gmsh mesh using this tool: https://github.com/NEFM-TUDresden/gmshModel. The meshes from the package work fine, if I am not using periodic BCs and I double-checked them to be periodic. The example I am using in the code is the hexagonal unit cell example from the repository. From what I have noticed, periodicity errors are only present for complex meshes - as here with inclusions on the boundary - and never in simple meshes.
To check for periodicity I am simply loading the mesh, generating a periodic domain for which the mapping seems to work fine and then I am creating a function on a FiniteElement function space (Lagrangian, first order). I have filled the vector of the function f with random numbers, before I applied the BCs explicitely. Afterwards, I am evaluating the values of the function f at the coordinates that I have logged during the mapping within the periodic subdomain. Comparing “master” and “slave” results with this method, leads to considerable errors that indicate that the periodicity mapping somehow is not working correctly, which is what I am experiencing in many of my simulations.
##############################################################################
# TEST FILE FOR THE IMPLEMENTATION OF PERIODIC BOUNDARY CONDITIONS IN FENICS #
##############################################################################
# load required libraries
#########################
from dolfin import * # dolfin for fenics fe capabilities
import numpy as np # numpy for fast array computations
import copy as cp
import matplotlib.pyplot as plt # matplotlib for figures
# Optimization options for the form compiler
############################################
parameters["form_compiler"]["cpp_optimize"] = True # optimize automatically generated cpp code
parameters["allow_extrapolation"] = True # must-do: otherwise this example would not work
###############################
# Function to read XDMF files #
###############################
def readXDMF(fileName,comm=MPI.comm_world,boundaries=False,subdomains=False):
# create empty mesh object
mesh = Mesh()
# read mesh
with XDMFFile(comm, fileName+"_domain.xdmf") as infile: # temporarily open meshFile (assume that mesh is stores as "file"_domain.xdmf)
infile.read(mesh) # read mesh
# get mesh dimension
dim=mesh.topology().dim()
# read boundary data if necessary
if boundaries: # boundary data is stored in mesh file
boundariesMVC = MeshValueCollection("size_t",mesh,dim=dim) # create mesh value collection
with XDMFFile(comm,fileName+"_boundaries.xdmf") as infile: # temporarily open boundary meshFile (assume that mesh is stored as "file"_boundaries.xdmf)
infile.read(boundariesMVC,'boundaries') # read boundary mesh
boundaryMarkers = cpp.mesh.MeshFunctionSizet(mesh,boundariesMVC) # create MeshFunction for boundary data
# read subdomain data if necessary
if subdomains:
subdomainsMVC = MeshValueCollection("size_t",mesh,dim=dim) # create mesh value collection
with XDMFFile(comm,fileName+"_domain.xdmf") as infile: # temporarily open boundary meshFile (assume that mesh is stored as "file"_domain.xdmf)
infile.read(subdomainsMVC,'subdomains') # read boundary mesh
subdomainMarkers = cpp.mesh.MeshFunctionSizet(mesh,subdomainsMVC) # create MeshFunction for subdomain data
return mesh, boundaryMarkers, subdomainMarkers
########################
# PeriodicDomain class #
########################
class PeriodicDomain(SubDomain):
def __init__(self,mesh,boundingBox,periodicAxes=None,tolerance=DOLFIN_EPS):
# calculate relevant information from bounding box array
axesMin=np.array(boundingBox[0])
axesMax=np.array(boundingBox[1])
relevantAxes=np.flatnonzero(axesMax-axesMin)
# initialize parent class
super().__init__()
# set relevant information for further calculations
self.relevantAxes=relevantAxes # relevant calculation axes
self.periodicAxes=periodicAxes # periodic axes
self.axesMin=axesMin # minimum axes coordinates (lower, left, back point of RVE)
self.axesMax=axesMax # maximum axes coordinates (upper, right, front point of RVE)
self.tolerance=tolerance # allowed tolerance
self.mappedMaster=[] # list of mapped master coordinates (for debugging)
self.mappedSlave=[] # list of mapped slave coordinates (for debugging)
##########################################################################
# Overwrite inside method of SubDomain class to define master subdomains #
##########################################################################
def inside(self,x,on_boundary):
# check if point x belongs to AT LEAST ONE master but NO slave entities
xInMasterEntity=False # initialze xInMasterEntity flag to be False by default
if on_boundary: # check if x is a boundary coordinate
for axis in self.periodicAxes: # loop over all periodic axes
if abs(x[axis]-self.axesMax[axis]) < self.tolerance: # -> x is part of the current axis slave entity
xInMasterEntity=False # ->-> set master entity membership information
break # ->-> stop further checks: x cannot be part of a master entity
elif abs(x[axis]-self.axesMin[axis]) < self.tolerance: # -> x is part of the current axis master entity
xInMasterEntity=True # ->-> set master entity membership information
return xInMasterEntity
#################################################################################
# Overwrite map method of SubDomain class to define the mapping master -> slave #
#################################################################################
def map(self,x,y):
xInSlaveEntity=False # initialize xInSlaveEntity flag to be False by default
for axis in self.relevantAxes: # loop over all axes with nonzero extent
if axis in self.periodicAxes: # current axis is supposed to be periodic
if abs(x[axis]-self.axesMax[axis]) < self.tolerance: # -> point is on slave entity for the current axis
y[axis]=self.axesMin[axis] # ->-> set coordinate to master nodes coordinate
xInSlaveEntity=True # ->-> update slave entity membership flag
else: # -> point is not on slave entity for the current axis
y[axis]=x[axis] # ->-> leave coordinate unchanged
else: # current axis is not supposed to be periodic
y[axis]=x[axis] # -> leave coordinate unchanged
# debugging information
if xInSlaveEntity:
self.mappedMaster.append(cp.deepcopy(y)) # add y to list mapped master coordinates
self.mappedSlave.append(cp.deepcopy(x)) # add x to list of mapped slave coordinates
# relevant simulation parameters:
#################################
# domain information
meshFile = "HexagonalCell" # mesh file name
domainBB=[[0, 0, 0],[6, 6*3**(0.5), 6*3**(0.5)]] # domain bounding box (lower front left and upper back right points)
# domain information (mesh, function spaces, ...)
#################################################
# import mesh
mesh, boundaries, subdomains = readXDMF(meshFile, MPI.comm_world, boundaries=True, subdomains=True) #read mesh
# generate subdomains for periodicity constraints and fixed center point
periodicDomain=PeriodicDomain(mesh,domainBB,periodicAxes=[0,1,2],tolerance=1e-4)# periodicity constraints
# define element and function spaces
V = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,V,constrained_domain=periodicDomain)
# boundary conditions
BC=DirichletBC(W,Constant(1.5),periodicDomain)
f = Function(W)
f.vector()[:]=np.random.rand(f.vector().size())
BC.apply(f.vector())
# check error in periodicity
periodicityError = np.zeros(len(periodicDomain.mappedMaster))
for bndPt in range(0,len(periodicDomain.mappedMaster)):
periodicityError[bndPt]=f(periodicDomain.mappedSlave[bndPt])-f(periodicDomain.mappedMaster[bndPt])
plt.plot(periodicityError)
plt.show()
Maybe someone of you has made similar observations or has a hint on what I might be doing wrong…
To make it a little easier:
The first thing that would help me while trying to find what is going wrong here is an answer to the question why the function spaces dof coordinates (W.tabulate_dof_coordinates) for Langrangian finite elements are not the same as the mesh coordinates - especially on the boundary where I find dofs that are supposed to be on the boundary (and they are, when I check mesh.coordinates()) but apparently are not.
