How to implement Bloch-Floquet Boundary condition?

1

Hello!

I’m working on solving the eigenvalue problem of the Laplace operator with Bloch-Floquet boundary conditions.
I already solved it with periodic boundary conditions:

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

Sub domain for Periodic boundary condition

class PeriodicBoundary(SubDomain):

# Left and bottom boundary is "target domain" G
def inside(self, x, on_boundary):
    return bool((near(x[0], (-Laenge/2), tol) or near(x[1], (-Weite/2), tol)) and 
        (not ((near(x[0], (-Laenge/2)) and near(x[1], (Weite/2))) or 
              (near(x[0], (Laenge/2)) and near(x[1], (-Weite/2))))) and on_boundary)

def map(self, x, y):
    if near(x[0], (Laenge/2)):
        y[0] = x[0] - Laenge
        y[1] = x[1] 
    else:   # near(x[1], Weite/2)
        y[0] = x[0]
        y[1] = x[1] - Weite

Create periodic boundary condition

pbc = PeriodicBoundary()

Create mesh and define function space

Laenge = 10.0
Weite = 10.0
tol = 1E-14
alpha = 1.0
domain = Rectangle(Point(-(Laenge/2), -(Weite/2)), Point(Laenge/2, Weite/2)) - Circle(Point(0,0), alpha, 100)
mesh = generate_mesh(domain, 20)
V = FunctionSpace(mesh, “Lagrange”, 1, constrained_domain=pbc)

Unfortunately I have no idea how to implement the Bloch-Floquet boundary condition.

Can anyone help me with this?

Thanks!

2

I am interested in this too. Is there a way to invoke a periodic BC with a phase factor \exp(\Gamma l) ?

For example, in solving the Maxwell system on a 1-D periodic structure, we have two planes \partial\Omega_1, \partial\Omega_2 connected to each other via the relationship:
\mathbf{n}\times\mathbf{E}\vert_{\partial\Omega_1} \exp(-\Gamma l) = \mathbf{n}\times\mathbf{E}\vert_{\partial\Omega_2}
where l is the unit cell length and \Gamma is an eigenvalue that depends on the frequency of operation. One can find a dispersion relationship \Gamma (k).

I see two problems:

  1. The weighted periodic boundary condition.
  2. Properly setting up the eigenvalue problem so it is tractable using Fenics.

The first is most important, since a search algorithm could be used to find \Gamma (not the best option). Generally, for a 2D or 3D structure, the mode set would be reflected in the computation of \Gamma.

3

My experience with Bloch boundary conditions comes from phononic materials, so I apologize that all of my references are in the context of elastodynamics. With that being said, I think there are three approaches you could use:

  1. Substitute the Bloch solution \boldsymbol{E}(\boldsymbol{x}) = \bar{\boldsymbol{E}} (\boldsymbol{x}) e^{i \boldsymbol{k} \cdot \boldsymbol{x}} into Maxwell’s equations and expand the spatial derivatives to rewrite the equations in terms of \bar{\boldsymbol{E}}(\boldsymbol{x}). The boundary conditions can then be expressed in terms of standard periodic boundary conditions \bar{\boldsymbol{E}}(\boldsymbol{x}) = \bar{\boldsymbol{E}}(\boldsymbol{x}+\boldsymbol{R}\boldsymbol{n}) (as opposed to a Floquet-Bloch boundary condition), since \bar{\boldsymbol{E}}(\boldsymbol{x}) is strictly periodic (not Floquet-periodic). The spatial derivatives of the Bloch propagator e^{i \boldsymbol{k} \cdot \boldsymbol{x}} give rise to additional terms in the weak form, see e.g. Eq. 2.13 of this paper for the weak form in the context of mechanics.
  2. If you use DOLFINx, I think you could probably use @dokken’s dolfinx_mpc to apply the Floquet-Bloch boundary condition via a multi-point constraint, i.e. by setting the scale equal to e^{i \boldsymbol{k} \cdot \boldsymbol{x}} in create_periodic_constraint_geometrical. However, I’ve never attempted this myself.
  3. Use the formulation of #1 to set up the inverse eigenvalue problem k(\Gamma), a la Collet et al. Here, again, the appropriate boundary conditions are the standard periodic boundary conditions, and the Floquet-Bloch propagator is not needed.
2 Likes
4

Thanks for that. You have given me some food for thought. The first option is what is often described in classic EM texts and may be the way to go. Cheers!

2 Likes
5

Did you manage to implement that? I am having the same issue now… So providing very minimal working code would be super helpful.

1 Like
6

Hello Ekrem,
Yes, I implemented it and it seems to work. Note that the code is for the older version of Fenics (2019). It should be easily adaptable to the new DolfinX version. The brief article describing numerical experiments as well as the code is here.

Code is also found at Github.
Have fun!

1 Like