FEniCS: задача с периодической границей в 3D

Question

Я хочу разрешить следующую систему:

Сильная проблема Неймана

в жидком домене квадрата, который является объединением двух доменов Omega_fluid и Omega_solid, последний представляет собой набор из 2- или 3D евклидовых чаш, пересекающихся с квадратом.

Пример составной ячейки

Дело в том, что результирующий вектор khi, а точнее его градиент, является непериодическим в 3D, если я вычисляю его с использованием FEniCS со следующим кодом:

## Usual python packages

import numpy as np
import pylab as pl
from pylab import *
import matplotlib.pyplot as plt
from math import sqrt

## FEniCS

from fenics import *
from dolfin import *
from mshr import *

## A periodic boundary defined as a SubDomain

tol=1e-10

xinf=0.0
yinf=0.0
zinf=0.0
xsup=1.0
ysup=1.0
zsup=1.0

dimension=3

class PeriodicBoundary(SubDomain):
 # Left boundary is "target domain" G
 def inside(self, x, on_boundary):
  return on_boundary and not(near(x[0],xsup,tol) or near(x[1],ysup,tol) or near(x[2],zsup,tol))## a special thank to Arnold Douglas
 # Map right boundary (H) to left boundary (G)
 def map(self, x, y):
  for i in range(dimension):
   if near(x[i],1.0,tol):
    y[i]=0.0
   else:
    y[i]=x[i]

## Creates a mesh for the unit square minus a circular inclusion

c_x=0.5
c_y=0.5
c_z=0.5
cen=[c_x,c_y,c_z]
r=0.35

res=15

box=Box(Point(-1,-1,-1),Point(2,2,2))
domain=box
# Creates a periodic inclusion, useful if cen != [0.5,0.5]
l_sph=[]
for i in range(-1,2):
 for j in range(-1,2):
  for k in range(-1,2):
   l_sph.append(Sphere(Point(cen[0]+i,cen[1]+j,cen[2]+k),r))
for sph_per in l_sph:
 domain=domain-sph_per
domain=domain*Box(Point(0,0,0),Point(1,1,1))

# A first mesh for the periodic domain
mesh_s_r=generate_mesh(domain,res)

## Refines the mesh at the edges of the square and of the inclusion
# not relevant

## Computes khi with a Finite Element Method

# A vector function space for the weak problem
V=VectorFunctionSpace(mesh_s_r, 'P', 2, constrained_domain=PeriodicBoundary())
# A SubDomain for the inclusion boundary
l_cen=[]
for i in range(-1,2):
 for j in range(-1,2):
  for k in range(-1,2):
   l_cen.append([cen[0]+i,cen[1]+j,cen[2]+k])
class periodic_inclusion(SubDomain):
 def inside(self,x,on_boundary):
  return (on_boundary and any([between((x[0]-c[0]), (-r-tol, r+tol)) for c in l_cen]) and any([between((x[1]-c[1]), (-r-tol, r+tol)) for c in l_cen]) and any([between((x[2]-c[2]), (-r-tol, r+tol)) for c in l_cen]))
## Marks for the domain boundaries
Gamma_sf = periodic_inclusion()
boundaries = MeshFunction("size_t", mesh_s_r, mesh_s_r.topology().dim()-1)
boundaries.set_all(1)
Gamma_sf.mark(boundaries, 7)
ds = Measure("ds")(subdomain_data=boundaries)
num_ff=1
num_front_sphere=7
# We solve the problem
normale = FacetNormal(mesh_s_r)
nb_noeuds=V.dim()
u = TrialFunction(V)
v = TestFunction(V)
a=tr(dot((grad(u)).T, grad(v)))*dx
L=-dot(normale,v)*ds(num_front_sphere)
##
u=Function(V)
solve(a==L,u)
# Mean value of khi
moy_u_x=assemble(u[0]*dx)/(1-4/3*pi*r**3)
moy_u_y=assemble(u[1]*dx)/(1-4/3*pi*r**3)
moy_u_z=assemble(u[2]*dx)/(1-4/3*pi*r**3)
moy=Function(V)
moy=Constant((moy_u_x,moy_u_y,moy_u_z))
khi=project(u-moy,V)

## Plots with matplotlib

for b in range(0,3):
 p_gk=plot(-grad(khi)[:,b])
 #plt.colorbar(p_gk)
 #plt.savefig("inc_s"+str(c_x)+str(c_y)+str(c_z)+"m_dkhi"+"_d"+str(b+1)+".png")
 plt.show(p_gk)
 plt.close()

Следующее векторное поле, которое является векторной компонентой Grad khi и, следовательно, не имеет расходимости, является сильно непериодическим.

Эта ситуация, по сути, характерна для третьего измерения: тот же код для 2D дает периодическое векторное поле, как на следующих рисунках:

Значения khi, черного и синего, на противоположных сторонах ячейки

Первый компонент вектора хи

Второй компонент вектора хи

Мы видим, что класс PeriodicBoundary() вызывается, когда пространство векторной функции определено, а сетка не должна быть периодической. Оправдывает ли это невыполнение периодического граничного условия? Однако это не объясняет того, что мы находим в двумерном случае: тангенциальная составляющая khi идеально совпадает на противоположных границах, нормальная составляющая совпадает, за исключением нескольких точек границы, где изменяется знак - я напечатал значения нормальный компонент, чтобы гарантировать, что это соответствует, где участок говорит "ноль".

Спасибо за помощь,

Антуан Моро (Университет Ла-Рошели)