Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/9_active_remeshing/active_remeshing_engine.py
0cbbfffdd8f5master

Last Change
Raw File

active_remeshing_engine.py
View Options

import numpy as np
import ufl
import fenics as fen
import mshr
from rrompy.utilities.base.decorators import (pivot_affine_construct,
mupivot_independent)
from rrompy.hfengines.fenics_engines import HelmholtzProblemEngine
from rrompy.utilities.numerical.hash_derivative import (
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenZERO, fenONE, fenics2Vector
from rrompy.parameter import parameterMap as pMap
from rrompy.utilities.exception_manager import RROMPyException
class ActiveRemeshingEngine(HelmholtzProblemEngine):
def __init__(self, z0:float, t0:float, n:int):
super().__init__(mu0 = [z0, t0])
self._affinePoly = False
self._nMesh = n
self.meshGen(t0)
self.parameterMap = pMap(1., 2)
self.DirichletBoundary = lambda x, on_boundary: (on_boundary
and np.abs(x[0]) >= .75 - 1e-5
or np.abs(x[1]) >= .5 - 1e-5)
self.NeumannBoundary = "REST"
self.cutOffPolesIMin, self.cutOffPolesIMax = -1e-2, 1e-2
def meshGen(self, t:float):
t = np.real(t)
if (not hasattr(self, "_tMesh") or not np.isclose(self._tMesh, t)):
e, self._tMesh = .01, t
tipx, tipy = .5 * np.sin(t), .5 - .5 * np.cos(t)
tiplx, tiply = tipx + .5 * e * np.cos(t), tipy + .5 * e * np.sin(t)
tiprx, tipry = tipx - .5 * e * np.cos(t), tipy - .5 * e * np.sin(t)
basx, basy = - e * np.sin(t), .5 + e * np.cos(t)
baslx, basly = basx + .5 * e * np.cos(t), basy + .5 * e * np.sin(t)
basrx, basry = basx - .5 * e * np.cos(t), basy - .5 * e * np.sin(t)
mesh = mshr.generate_mesh(
mshr.Rectangle(fen.Point(-.75, -.5), fen.Point(.75, .5))
- mshr.Polygon([fen.Point(tiprx, tipry),
fen.Point(tiplx, tiply),
fen.Point(baslx, basly),
fen.Point(basrx, basry)])
- mshr.Circle(fen.Point(tipx, tipy), .5 * e), self._nMesh)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.As, self._C = [None] * 2, None
self.autoSetDS()
if hasattr(self, "_energyNormMatrix"):
del self._energyNormMatrix
if hasattr(self, "_energyNormDualMatrix"):
del self._energyNormDualMatrix
def getForcingTerm(self, mu = []):
mu = self.checkParameter(mu)
self.meshGen(mu(0, 1))
x, y = fen.SpatialCoordinate(self.V.mesh())[:]
rightZone = .1875**-2 * ufl.conditional(ufl.And(
ufl.And(ufl.ge(x, -.5625), ufl.le(x, -.375)),
ufl.And(ufl.ge(y, .125), ufl.le(y, .3125))),
fenONE, fenZERO)
return rightZone, fenZERO
@pivot_affine_construct
def A(self, mu = [], der = 0):
derI = hashD(der) if hasattr(der, "__len__") else der
if derI > 0: raise Exception("Derivatives not implemented.")
mu = self.checkParameter(mu)
self.meshGen(mu(0, 1))
return HelmholtzProblemEngine.A(self, mu, der)
@pivot_affine_construct
def b(self, mu = [], der = 0):
derI = hashD(der) if hasattr(der, "__len__") else der
if derI > 0: raise Exception("Derivatives not implemented.")
if self.thbs[0] is None: self.thbs = self.getMonomialWeights(self.nbs)
fen0 = self.getForcingTerm(mu)[0] * self.v * fen.dx
DBC = fen.DirichletBC(self.V, fenZERO, self.DirichletBoundary)
self.bs = [fenics2Vector(fen0, {}, DBC, 1)]
return HelmholtzProblemEngine.b(self, mu, der)
@mupivot_independent
def C(self, mu):
mu = self.checkParameterList(mu)
if not np.all(np.isclose(mu(1), mu(0, 1))):
raise RROMPyException(("Simultaneous evaluation of C on multiple "
"meshes not supported."))
self.meshGen(mu(0, 1))
if self._C is None:
self._C = fenics2Vector(self.v * fen.dx, {}).reshape(1, -1) / 1.5
return self._C
c4science · Help