Page MenuHomec4science
Files F61737829

membrane_fracture_engine.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 8, 16:04

membrane_fracture_engine.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import fenics as fen
import mshr, ufl
from rrompy.utilities.base.types import ScOp, List, paramVal
from rrompy.solver.fenics import fenZERO, fenONE
from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
HelmholtzProblemEngine)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenics2Sparse
__all__ = ['MembraneFractureEngine']
class MembraneFractureEngine(HelmholtzProblemEngine):
def __init__(self, mu0 : paramVal = [20. ** .5, .6], H : float = 1.,
L : float = .75, delta : float = .05, n : int = 50,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.nAs = 20
self.npar = 2
self.H = H
self.rescalingExp = [2., 1.]
domain = (mshr.Rectangle(fen.Point(0., - H / 2.),
fen.Point(2. * L + delta, H / 2.))
- mshr.Rectangle(fen.Point(L, 0.),
fen.Point(L + delta, H / 2.)))
mesh = mshr.generate_mesh(domain, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
and x[0] >= L
and x[0] <= L + delta)
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
self._belowIndicator = ufl.conditional(ufl.le(y, 0.), fenONE, fenZERO)
self._aboveIndicator = fenONE - self._belowIndicator
self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
- .5 * ((x - .6 * L) / (.1 * L)) ** 2.
) * self._belowIndicator, fenZERO]
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
self.autoSetDS()
for j in [1, 3, 4, 6, 7, 10, 11, 12, 15, 16, 17, 18]:
if derI <= j and self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
if derI <= 0 and self.As[0] is None:
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a0Re = (self.H ** 4 / 4. * self._aboveIndicator
* fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx)
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 2 and self.As[2] is None:
vbMng(self, "INIT", "Assembling operator term A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a2Re = (- self.H ** 3 / 2. * self._aboveIndicator
* fen.dot(self.u.dx(1), self.v.dx(1)) * fen.dx)
self.As[2] = fenics2Sparse(a2Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 5 and self.As[5] is None:
vbMng(self, "INIT", "Assembling operator term A6.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a5Re = self.H ** 2 * (fen.dot(self.u.dx(0), self.v.dx(0))
+ .25 * fen.dot(self.u.dx(1), self.v.dx(1))) * fen.dx
self.As[5] = fenics2Sparse(a5Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 8 and self.As[8] is None:
vbMng(self, "INIT", "Assembling operator term A8.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a8Re = - self.H ** 2. * n2Re * fen.dot(self.u, self.v) * fen.dx
a8Im = - self.H ** 2. * n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[8] = (fenics2Sparse(a8Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a8Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 9 and self.As[9] is None:
vbMng(self, "INIT", "Assembling operator term A9.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a9Re = - 2. * self.H * fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[9] = fenics2Sparse(a9Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 13 and self.As[13] is None:
vbMng(self, "INIT", "Assembling operator term A13.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a13Re = 2. * self.H * n2Re * fen.dot(self.u, self.v) * fen.dx
a13Im = 2. * self.H * n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[13] = (fenics2Sparse(a13Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a13Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 14 and self.As[14] is None:
vbMng(self, "INIT", "Assembling operator term A14.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
a14Re = fen.dot(self.u.dx(0), self.v.dx(0)) * fen.dx
self.As[14] = fenics2Sparse(a14Re, {}, DirichletBC0, 0)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 19 and self.As[19] is None:
vbMng(self, "INIT", "Assembling operator term A19.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
nRe, nIm = self.refractionIndex
n2Re, n2Im = nRe * nRe - nIm * nIm, 2 * nRe * nIm
parsRe = self.iterReduceQuadratureDegree(zip([n2Re],
["refractionIndexSquaredReal"]))
parsIm = self.iterReduceQuadratureDegree(zip([n2Im],
["refractionIndexSquaredImag"]))
a19Re = - n2Re * fen.dot(self.u, self.v) * fen.dx
a19Im = - n2Im * fen.dot(self.u, self.v) * fen.dx
self.As[19] = (fenics2Sparse(a19Re, parsRe, DirichletBC0, 0)
+ 1.j * fenics2Sparse(a19Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
return self._assembleA(mu, der, derI)

Event Timeline

c4science · Help