Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/Python_FEniCS/FEniCS_snippets.py
6e6c171cf128master

Last Change
Raw File

FEniCS_snippets.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import sympy as sp
import fenics as fen
PI = np.pi
def SquareHomogeneousBubble(kappa:float, theta:float, n:int):
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI, PI), n, n)
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 16/PI**4 * x * y * (x - PI) * (y - PI)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
sp.printing.ccode(sp.simplify(sp.im(fex)))]
return "ALL", mesh, forcingTerm
def SquareTransmissionDirichlet(nT:float, nB:float, kappa:float,
theta:float, n:int, deg : int = 4):
mesh = fen.RectangleMesh(fen.Point(-PI/2, -PI/2), fen.Point(PI/2, PI/2),
n, n)
nIdx = fen.Expression('x[1] >= 0 ? nT : nB', nT = nT, nB = nB, degree = deg)
dx, dy = np.cos(theta), np.sin(theta)
Kx = kappa * nB * dx
Ky = kappa * (nT**2. - (nB * dx)**2. + 0.j)**.5
T = 2 * kappa * nB * dy / (Ky + kappa * nB * dy)
x, y = sp.symbols('x[0] x[1]', real=True)
uT = T * sp.exp(1.j * (Kx*x + Ky*y))
uB = ( sp.exp(1.j * kappa * nB * (dx*x + dy*y))
+ (T - 1) * sp.exp(1.j * kappa * nB * (dx*x - dy*y)))
uRe = fen.Expression('x[1]>=0 ? {} : {}'.format(
sp.printing.ccode(sp.re(uT)), sp.printing.ccode(sp.re(uB))),
degree = deg)
uIm = fen.Expression('x[1]>=0 ? {} : {}'.format(
sp.printing.ccode(sp.im(uT)), sp.printing.ccode(sp.im(uB))),
degree = deg)
return "ALL", mesh, nIdx, (uRe, uIm)
def SquareScatteringTB(kappa:float, theta:float, n:int):
def Dboundary(x, on_boundary):
return on_boundary and (fen.near(x[0], 0) or fen.near(x[0], PI))
mesh = fen.RectangleMesh(fen.Point(0,0), fen.Point(PI, PI), n, n)
x, y = sp.symbols('x[0] x[1]', real=True)
wex = 4/PI**2 * x * (PI - x)
phiex = kappa * (x * np.cos(theta) + y * np.sin(theta))
uex = wex * sp.exp(-1.j * phiex)
fex = - uex.diff(x, 2) - uex.diff(y, 2) - kappa**2 * uex
forcingTerm = [sp.printing.ccode(sp.simplify(sp.re(fex))),
sp.printing.ccode(sp.simplify(sp.im(fex)))]
return Dboundary, "rest", mesh, forcingTerm
c4science · Help