Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/hfengines/linear_problem/tridimensional/membrane_fracture_engine3.py
1818294aa86adevelop

Last Change
Raw File

membrane_fracture_engine3.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 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__ = ['MembraneFractureEngine3']
def gen_mesh_fracture(W, delta, H, n):
n = 2 * (n // 2) + 1
xL = np.linspace(-.5 * W, - .5 * delta, n // 2 + 1)
xR = - xL[::-1]
nEff = int(np.ceil(delta / (xL[1] - xL[0])))
xC = np.linspace(-.5 * delta, .5 * delta, nEff + 1)[1:-1]
yA = np.linspace(-.5 * H, .5 * H, n)
yL = np.linspace(-.5 * H, 0., n // 2 + 1)
editor = fen.MeshEditor()
mesh = fen.Mesh()
editor.open(mesh, "triangle", 2, 2)
editor.init_vertices((2 * n + nEff - 1) * (n // 2 + 1))
editor.init_cells(2 * (2 * n - 2 + nEff) * (n // 2))
for j, y in enumerate(yA):
for i, x in enumerate(xL):
editor.add_vertex(i + j * (n // 2 + 1), np.array([x, y]))
for j in range(n - 1):
for i in range(n // 2):
editor.add_cell(2 * (i + j * (n // 2)),
np.array([i + j * (n // 2 + 1), i + 1 + j * (n // 2 + 1),
i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
editor.add_cell(2 * (i + j * (n // 2)) + 1,
np.array([i + 1 + j * (n // 2 + 1), i + 1 + (j + 1) * (n // 2 + 1),
i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
vBase1, cBase1 = n * (n // 2 + 1), 2 * (n - 1) * (n // 2)
for j, y in enumerate(yA):
for i, x in enumerate(xR):
editor.add_vertex(vBase1 + i + j * (n // 2 + 1), np.array([x, y]))
for j in range(n - 1):
for i in range(n // 2):
editor.add_cell(cBase1 + 2 * (i + j * (n // 2)),
np.array([vBase1 + i + j * (n // 2 + 1), vBase1 + i + 1 + j * (n // 2 + 1),
vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
editor.add_cell(cBase1 + 2 * (i + j * (n // 2)) + 1,
np.array([vBase1 + i + 1 + j * (n // 2 + 1),
vBase1 + i + 1 + (j + 1) * (n // 2 + 1),
vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
vBase2, cBase2 = 2 * n * (n // 2 + 1), 4 * (n - 1) * (n // 2)
for j, y in enumerate(yL):
for i, x in enumerate(xC):
editor.add_vertex(vBase2 + i + j * (nEff - 1), np.array([x, y]))
for j in range(n // 2):
for i in range(nEff - 2):
editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)),
np.array([vBase2 + i + j * (nEff - 1), vBase2 + i + 1 + j * (nEff - 1),
vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)) + 1,
np.array([vBase2 + i + 1 + j * (nEff - 1),
vBase2 + i + 1 + (j + 1) * (nEff - 1),
vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
if nEff == 1:
for j in range(n // 2):
editor.add_cell(cBase2 + 2 * j,
np.array([(j + 1) * (n // 2 + 1) - 1, vBase1 + j * (n // 2 + 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
editor.add_cell(cBase2 + 2 * j + 1,
np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
else:
cBase3 = 2 * (2 * n + nEff - 4) * (n // 2)
for j in range(n // 2):
editor.add_cell(cBase3 + 2 * j,
np.array([(j + 1) * (n // 2 + 1) - 1, vBase2 + j * (nEff - 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
editor.add_cell(cBase3 + 2 * j + 1,
np.array([vBase2 + j * (nEff - 1), vBase2 + (j + 1) * (nEff - 1),
(j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
cBase4 = 2 * (2 * n + nEff - 3) * (n // 2)
for j in range(n // 2):
editor.add_cell(cBase4 + 2 * j,
np.array([vBase2 + (j + 1) * (nEff - 1) - 1, vBase1 + j * (n // 2 + 1),
vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
editor.add_cell(cBase4 + 2 * j + 1,
np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
editor.close()
return mesh
class MembraneFractureEngine3(HelmholtzProblemEngine):
def __init__(self, mu0 : paramVal = [20. ** .5, .6, .08], 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 = 206
self.npar = 3
self._H = H
self._delta = delta
self._L = L
self._W = 2 * L + delta
self.rescalingExp = [2., 1., 1.]
mesh = gen_mesh_fracture(self._W, delta, H, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
and x[0] >= - .5 * delta and x[0] <= .5 * delta)
self.DirichletBoundary = "REST"
x, y = fen.SpatialCoordinate(mesh)[:]
self._aboveIndicator = ufl.conditional(ufl.gt(y, 0.), fenONE, fenZERO)
self._belowCenterIndicator = ufl.conditional(
ufl.And(ufl.ge(x, - .5 * delta),
ufl.le(x, .5 * delta)),
fenONE, fenZERO)
self._belowSidesIndicator = (fenONE - self._aboveIndicator
- self._belowCenterIndicator)
self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
- .5 * ((x + .4 * L + .5 * delta) / (.1 * L)) ** 2.
) * self._belowSidesIndicator, 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()
spKx, spKy = [None] * 2, [None] * 2
spM = None
def getstiffnessblockj(direc:int, j:int):
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
if direc == 0:
if j == 0:
fact = 16. * self._L ** 2.
else:
fact = 4. * self._delta ** 2.
else:
fact = self._H ** 2.
if direc == 0:
if j == 0:
fact *= (self._aboveIndicator + self._belowSidesIndicator)
else:
fact *= self._belowCenterIndicator
else:
if j == 0:
fact *= self._aboveIndicator
else:
fact *= (self._belowSidesIndicator
+ self._belowCenterIndicator)
ajRe = fact * self.u.dx(direc) * self.v.dx(direc) * fen.dx
return fenics2Sparse(ajRe, {}, DirichletBC0, 0)
def getmass():
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = - 4. * self.u * self.v * fen.dx
return fenics2Sparse(ajRe, {}, DirichletBC0, 0)
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 = fenZERO * fen.dot(self.u, self.v) * fen.dx
self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 7 and self.As[7] is None:
vbMng(self, "INIT", "Assembling operator term A7.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[7] = self._W ** 2. * self._H ** 2. * spKx[1]
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)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[9] = self._W ** 2. * self._H ** 2. * spKy[0]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 16 and self.As[16] is None:
vbMng(self, "INIT", "Assembling operator term A16.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[16] = - 2. * self._W ** 2. * self._H * spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 17 and self.As[17] is None:
vbMng(self, "INIT", "Assembling operator term A17.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[17] = - 2. * self._W * self._H ** 2. * spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 18 and self.As[18] is None:
vbMng(self, "INIT", "Assembling operator term A18.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[18] = - 2. * self._W ** 2. * self._H * spKy[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)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[19] = - 2. * self._W * self._H ** 2. * spKy[0]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 30 and self.As[30] is None:
vbMng(self, "INIT", "Assembling operator term A30.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[30] = self._W ** 2. * spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 31 and self.As[31] is None:
vbMng(self, "INIT", "Assembling operator term A31.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[31] = 4. * self._W * self._H * spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 32 and self.As[32] is None:
vbMng(self, "INIT", "Assembling operator term A32.", 20)
if spKx[0] is None: spKx[0] = getstiffnessblockj(0, 0)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
if spKy[1] is None: spKy[1] = getstiffnessblockj(1, 1)
self.As[32] = (self._H ** 2. * (spKx[0] + spKx[1])
+ self._W ** 2. * (spKy[0] + spKy[1]))
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 33 and self.As[33] is None:
vbMng(self, "INIT", "Assembling operator term A33.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[33] = 4. * self._W * self._H * spKy[0]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 34 and self.As[34] is None:
vbMng(self, "INIT", "Assembling operator term A34.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[34] = self._H ** 2. * spKy[0]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 47 and self.As[47] is None:
vbMng(self, "INIT", "Assembling operator term A47.", 20)
if spM is None: spM = getmass()
self.As[47] = self._W ** 2. * self._H ** 2. * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 51 and self.As[51] is None:
vbMng(self, "INIT", "Assembling operator term A51.", 20)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[51] = - 2. * self._W * spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 52 and self.As[52] is None:
vbMng(self, "INIT", "Assembling operator term A52.", 20)
if spKx[0] is None: spKx[0] = getstiffnessblockj(0, 0)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[52] = - 2. * self._H * (spKx[0] + spKx[1])
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 53 and self.As[53] is None:
vbMng(self, "INIT", "Assembling operator term A53.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
if spKy[1] is None: spKy[1] = getstiffnessblockj(1, 1)
self.As[53] = - 2. * self._W * (spKy[0] + spKy[1])
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 54 and self.As[54] is None:
vbMng(self, "INIT", "Assembling operator term A54.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
self.As[54] = - 2. * self._H * spKy[0]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 73 and self.As[73] is None:
vbMng(self, "INIT", "Assembling operator term A73.", 20)
if spM is None: spM = getmass()
self.As[73] = - 2. * self._W ** 2. * self._H * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 74 and self.As[74] is None:
vbMng(self, "INIT", "Assembling operator term A74.", 20)
if spM is None: spM = getmass()
self.As[74] = - 2. * self._W * self._H ** 2. * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 79 and self.As[79] is None:
vbMng(self, "INIT", "Assembling operator term A79.", 20)
if spKx[0] is None: spKx[0] = getstiffnessblockj(0, 0)
if spKx[1] is None: spKx[1] = getstiffnessblockj(0, 1)
self.As[79] = spKx[0] + spKx[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 81 and self.As[81] is None:
vbMng(self, "INIT", "Assembling operator term A81.", 20)
if spKy[0] is None: spKy[0] = getstiffnessblockj(1, 0)
if spKy[1] is None: spKy[1] = getstiffnessblockj(1, 1)
self.As[81] = spKy[0] + spKy[1]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 107 and self.As[107] is None:
vbMng(self, "INIT", "Assembling operator term A107.", 20)
if spM is None: spM = getmass()
self.As[107] = self._W ** 2. * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 108 and self.As[108] is None:
vbMng(self, "INIT", "Assembling operator term A108.", 20)
if spM is None: spM = getmass()
self.As[108] = 4. * self._W * self._H * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 109 and self.As[109] is None:
vbMng(self, "INIT", "Assembling operator term A109.", 20)
if spM is None: spM = getmass()
self.As[109] = self._H ** 2. * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 151 and self.As[151] is None:
vbMng(self, "INIT", "Assembling operator term A151.", 20)
if spM is None: spM = getmass()
self.As[151] = - 2. * self._W * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 152 and self.As[152] is None:
vbMng(self, "INIT", "Assembling operator term A152.", 20)
if spM is None: spM = getmass()
self.As[152] = - 2. * self._H * spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
if derI <= 205 and self.As[205] is None:
vbMng(self, "INIT", "Assembling operator term A205.", 20)
if spM is None: spM = getmass()
self.As[205] = spM
vbMng(self, "DEL", "Done assembling operator term.", 20)
for j in range(derI, self.nAs):
if self.As[j] is None:
self.As[j] = self.checkAInBounds(-1)
return self._assembleA(mu, der, derI)
c4science · Help