Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/2_double_slit/double_slit_engine.py
6c4eee04e0dfmaster

Last Change
Raw File

double_slit_engine.py
View Options

import numpy as np
import ufl
import fenics as fen
import mshr
from rrompy.utilities.base.decorators import (affine_construct,
nonaffine_construct)
from rrompy.hfengines.fenics_engines import ScatteringProblemEngine
from rrompy.utilities.numerical.hash_derivative import (
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenZERO, fenics2Vector
class DoubleSlitEngine(ScatteringProblemEngine):
def __init__(self, k0:float, n:int):
super().__init__(mu0 = [k0])
self._affinePoly = False
delta, eps = .1, .01
mesh = mshr.generate_mesh(
mshr.Circle(fen.Point(0., 0.), 5.)
- mshr.Rectangle(fen.Point(-5., -delta), fen.Point(-.75, delta))
- mshr.Rectangle(fen.Point(-.5, -delta), fen.Point(.5, delta))
- mshr.Rectangle(fen.Point(.75, -delta), fen.Point(5., delta)), n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.DirichletBoundary = lambda x, on_boundary: (on_boundary
and np.abs(x[1]) <= delta
and np.abs(x[1]) > delta - eps)
self.NeumannBoundary = lambda x, on_boundary: (on_boundary
and np.abs(x[1]) <= delta - eps)
self.RobinBoundary = "REST"
def getDirichletValues(self, mu = []):
mu = self.checkParameter(mu)
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
c, s = .5, - .5 * 3. ** .5
muR, muI = np.real(mu[0])[0], np.imag(mu[0])[0]
mod = - muI * (c * x + s * y)
angle = muR * (c * x + s * y)
DR = fen.exp(mod) * fen.cos(angle)
DI = fen.exp(mod) * fen.sin(angle)
DR = ufl.conditional(ufl.ge(y, 0), DR, fenZERO)
DI = ufl.conditional(ufl.ge(y, 0), DI, fenZERO)
return DR, DI
@affine_construct
def A(self, mu = [], der = 0):
return ScatteringProblemEngine.A(self, mu, der)
@nonaffine_construct
def b(self, mu = [], der = 0):
derI = hashD(der) if hasattr(der, "__len__") else der
if derI < 0: return self.baselineb()
if derI > 0: raise Exception("Derivatives not implemented.")
fen0 = fen.inner(fenZERO, self.v) * fen.dx
DR, DI = self.getDirichletValues(mu)
DBCR = fen.DirichletBC(self.V, DR, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, DI, self.DirichletBoundary)
return (fenics2Vector(fen0, {}, DBCR, 1)
+ 1.j * fenics2Vector(fen0, {}, DBCI, 1))
c4science · Help