Page MenuHomec4science
Files F61613045

anisotropic_square_engine.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 7, 19:56

anisotropic_square_engine.py
View Options

import numpy as np
import fenics as fen
import mshr
import ufl
from rrompy.hfengines.fenics_engines import HelmholtzProblemEngine
from rrompy.solver.fenics import fenONE, fenZERO, fenics2Sparse
from rrompy.parameter import parameterMap as pMap
class AnisotropicSquareEngine(HelmholtzProblemEngine):
def __init__(self, k2:float, L2:float, n:int):
super().__init__(mu0 = [k2, L2])
self._affinePoly = True
self.nAs = 3
self.npar = 2
self.parameterMap = pMap(1., 2)
self.V = fen.FunctionSpace(fen.UnitSquareMesh(n, n), "P", 1)
eps = 1e-6
self.DirichletBoundary = lambda x, on_boundary: (on_boundary
and x[1] < eps)
self.NeumannBoundary = "REST"
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
self.NeumannDatum = ufl.conditional(ufl.ge(y, 1. - eps),
fen.cos(np.pi * x), fenZERO)
self.forcingTerm = ufl.conditional(ufl.ge(y, .5), fenONE, fenZERO) * (
5 * ufl.conditional(ufl.lt(x, .1), fenONE, fenZERO)
- 5 * ufl.conditional(ufl.And(ufl.gt(x, .2), ufl.lt(x, .3)),
fenONE, fenZERO)
+ 10 * ufl.conditional(ufl.And(ufl.gt(x, .45), ufl.lt(x, .55)),
fenONE, fenZERO)
- 5 * ufl.conditional(ufl.And(ufl.gt(x, .7), ufl.lt(x, .8)),
fenONE, fenZERO)
+ 5 * ufl.conditional(ufl.gt(x, .9), fenONE, fenZERO))
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0 = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[0] = fenics2Sparse(a0, {}, DirichletBC0, 1)
if self.As[1] is None:
self.autoSetDS()
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a1 = - fen.dot(self.u, self.v) * fen.dx
self.As[1] = fenics2Sparse(a1, {}, DirichletBC0, 0)
if self.As[2] is None:
self.autoSetDS()
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a2 = fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx
self.As[2] = fenics2Sparse(a2, {}, DirichletBC0, 0)
def AnisotropicSquareEnginePoles(L2:float, k2min:float, k2max:float):
poles = []
for alpha in np.arange(np.ceil((k2max) ** .5 / np.pi)):
p = (np.pi * alpha) ** 2.
pkmin = np.ceil(max(0., (k2min - p) * 4 / L2) ** .5 / np.pi)
pkmin += 1 - (pkmin % 2)
pkmax = np.floor(max(0., (k2max - p) * 4 / L2) ** .5 / np.pi)
for beta in np.arange(pkmin, pkmax + 1, 2):
poles += [p + L2 * (np.pi * beta / 2.) ** 2.]
return np.unique(poles)

Event Timeline

c4science · Help