Page MenuHomec4science
Files F63113540

membrane_fracture_engine3.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 17, 20:58

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 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.numerical import (nextDerivativeIndices,
hashDerivativeToIdx as hashD)
from rrompy.solver.fenics import fenics2Sparse
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._affinePoly = False
self.nAs = 6
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 buildA(self):
"""Build terms of operator of linear system."""
if self.As[1] is None or self.As[2] is None or self.As[5] is None:
thy2 = [(1.,), ('x', '()', 1, '*', 4., '-', 2. * self._H),
('x', '()', 1, '**', 2., '*', 6., '-',
('x', '()', 1, '*', 6. * self._H), '+', self._H ** 2.),
('x', '()', 1, '**', 2., '*', 4., '-',
('x', '()', 1, '*', 3. * self._H), '+', self._H ** 2.,
'*', ('x', '()', 1), '*', 2.),
('x', '()', 1, '*', -1., '+', self._H, '*', ('x', '()', 1),
'**', 2.)]
if self.As[3] is None or self.As[4] is None or self.As[5] is None:
thz2 = [(1.,), ('x', '()', 2, '*', 4., '-', 2. * self._W),
('x', '()', 2, '**', 2., '*', 6., '-',
('x', '()', 2, '*', 6. * self._W), '+', self._W ** 2.),
('x', '()', 2, '**', 2., '*', 4., '-',
('x', '()', 2, '*', 3. * self._W), '+', self._W ** 2.,
'*', ('x', '()', 2), '*', 2.),
('x', '()', 2, '*', -1., '+', self._W, '*', ('x', '()', 2),
'**', 2.)]
if 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)
self.thAs[0] = self.getMonomialSingleWeight([0])
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[1] is None:
vbMng(self, "INIT", "Assembling operator term A1.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = ((self._aboveIndicator + self._belowSidesIndicator)
* 16. * self._L ** 2. * self.u.dx(0) * self.v.dx(0) * fen.dx)
self.As[1] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thz0 = [(1.,), ('x', '()', 2, '*', 2.), ('x', '()', 2, '**', 2.)]
self.thAs[1] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
for idx in idxs:
dy, dz = 4 - idx[1], 2 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[1] += [thy2[dy] + ('*') + (thz0[dz],)]
else:
self.thAs[1] += [(0.,)]
self.thAs[1] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[2] is None:
vbMng(self, "INIT", "Assembling operator term A2.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = (4. * self._delta ** 2. * self._belowCenterIndicator
* self.u.dx(0) * self.v.dx(0) * fen.dx)
self.As[2] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thz1 = [(1.,), ('x', '()', 2, '*', 2., '-', 2. * self._W),
('x', '()', 2, '*', -1., '+', self._W, '**', 2.)]
self.thAs[2] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
for idx in idxs:
dy, dz = 4 - idx[1], 2 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[2] += [thy2[dy] + ('*') + (thz1[dz],)]
else:
self.thAs[2] += [(0.,)]
self.thAs[2] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[3] is None:
vbMng(self, "INIT", "Assembling operator term A3.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = (self._H ** 2. * self._aboveIndicator
* self.u.dx(1) * self.v.dx(1) * fen.dx)
self.As[3] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thy1 = [(1.,), ('x', '()', 1, '*', 2., '-', 2. * self._H),
(self._H, '-', ('x', '()', 1), '**', 2.)]
self.thAs[3] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
for idx in idxs:
dy, dz = 2 - idx[1], 4 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[3] += [thy1[dy] + ('*') + (thz2[dz],)]
else:
self.thAs[3] += [(0.,)]
self.thAs[3] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[4] is None:
vbMng(self, "INIT", "Assembling operator term A4.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = ((self._belowSidesIndicator + self._belowCenterIndicator)
* self._H ** 2. * self.u.dx(1) * self.v.dx(1) * fen.dx)
self.As[4] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thy0 = [(1.,), ('x', '()', 1, '*', 2.), ('x', '()', 1, '**', 2.)]
self.thAs[4] = []
idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
for idx in idxs:
dy, dz = 2 - idx[1], 4 - idx[2]
if idx[0] == 0 and dy >= 0 and dz >= 0:
self.thAs[4] += [thy0[dy] + ('*') + (thz2[dz],)]
else:
self.thAs[4] += [(0.,)]
self.thAs[4] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)
if self.As[5] is None:
vbMng(self, "INIT", "Assembling operator term A5.", 20)
DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
self.DirichletBoundary)
ajRe = - 4. * self.u * self.v * fen.dx
self.As[5] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
thx = [(1.,), ('x', '()', 0)]
self.thAs[5] = []
idxs = nextDerivativeIndices([], 3, hashD([1, 4, 4]) + 1)
for idx in idxs:
dx, dy, dz = 1 - idx[0], 4 - idx[1], 4 - idx[2]
if dx >= 0 and dy >= 0 and dz >= 0:
self.thAs[5] += [thx[dx] + ('*') + (thy2[dy],)
+ ('*') + (thz2[dz],)]
else:
self.thAs[5] += [(0.,)]
self.thAs[5] += [None]
vbMng(self, "DEL", "Done assembling operator term.", 20)

Event Timeline

c4science · Help