diff --git a/examples/2d/base/cookie_single.py b/examples/2d/base/cookie_single.py
new file mode 100644
index 0000000..f153e01
--- /dev/null
+++ b/examples/2d/base/cookie_single.py
@@ -0,0 +1,29 @@
+import numpy as np
+import fenics as fen
+from rrompy.hfengines.linear_problem.bidimensional import \
+ CookieEngineSingle as CES
+
+verb = 100
+
+kappa = 15. ** .5
+theta = - np.pi / 6.
+n = 30
+R = 1.
+L = np.pi
+nX = 3
+nY = 2
+
+mu0 = [25. ** .5, 1.]
+mutar = [25. ** .5, 1.5]
+
+solver = CES(kappa = kappa, theta = theta, n = n, R = R, L = L, nX = nX,
+ nY = nY, mu0 = mu0, verbosity = verb)
+uh = solver.solve(mutar)[0]
+solver.plotmesh(figsize = (7.5, 4.5))
+fen.plot(fen.project(solver.cookieIn,
+ fen.FunctionSpace(solver.V.mesh(), "DG", 0)))
+print(solver.norm(uh))
+solver.plot(uh)
+solver.outParaviewTimeDomain(uh, mutar[0], filename = 'out', folder = True,
+ forceNewFile = False)
+
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/cookie_single_pod.py
similarity index 58%
copy from examples/2d/pod/fracture_pod.py
copy to examples/2d/pod/cookie_single_pod.py
index b494062..752762f 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/cookie_single_pod.py
@@ -1,177 +1,138 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
- MembraneFractureEngine as MFE
+ CookieEngineSingle as CES
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 100
-show_sample = False
+size = 1
+show_sample = True
show_norm = True
-ignore_forcing = True
-ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
-homogeneize = False
-#homogeneize = True
-MN = 1
+MN = 15
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if size == 1: # below
- mu0 = [40 ** .5, .4]
- mutar = [45 ** .5, .4]
- murange = [[30 ** .5, .3], [50 ** .5, .5]]
-elif size == 2: # top
- mu0 = [40 ** .5, .6]
- mutar = [45 ** .5, .6]
- murange = [[30 ** .5, .5], [50 ** .5, .7]]
-elif size == 3: # interesting
- mu0 = [40 ** .5, .5]
- mutar = [45 ** .5, .5]
- murange = [[30 ** .5, .3], [50 ** .5, .7]]
-elif size == 4: # wide_low
- mu0 = [40 ** .5, .2]
- mutar = [45 ** .5, .2]
- murange = [[10 ** .5, .1], [70 ** .5, .3]]
-elif size == 5: # wide_hi
- mu0 = [40 ** .5, .8]
- mutar = [45 ** .5, .8]
- murange = [[10 ** .5, .7], [70 ** .5, .9]]
-elif size == 6: # top_zoom
- mu0 = [50 ** .5, .8]
- mutar = [55 ** .5, .8]
- murange = [[40 ** .5, .7], [60 ** .5, .9]]
-elif size == 100: # tiny
- mu0 = [32.5 ** .5, .5]
- mutar = [34 ** .5, .5]
- murange = [[30 ** .5, .3], [35 ** .5, .7]]
+ mu0 = [20 ** .5, 1. ** .5]
+ mutar = [20.5 ** .5, 1.05 ** .5]
+ murange = [[18.5 ** .5, .85 ** .5], [21.5 ** .5, 1.15 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
-H = 1.
-L = .75
-delta = .05
-n = 25
+kappa = 20. ** .5
+theta = - np.pi / 6.
+n = 40
+Rad = 1.
+L = np.pi
+nX = 2
+nY = 1
-solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
+solver = CES(kappa = kappa, theta = theta, n = n, R = Rad, L = L, nX = nX,
+ nY = nY, mu0 = mu0, verbosity = verb)
-rescaling = [lambda x: np.power(x, 2.), lambda x: x]
-rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
+rescaling = [lambda x: np.power(x, 2.)] * 2
+rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
-approx = method(solver, mu0 = mu0, approxParameters = params,
- verbosity = verb, homogeneized = homogeneize)
+approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
- import ufl
- import fenics as fen
- from rrompy.solver.fenics import affine_warping
L = mutar[1]
- y = fen.SpatialCoordinate(solver.V.mesh())[1]
- warp1, warpI1 = affine_warping(solver.V.mesh(),
- np.array([[1, 0], [0, 2. * L]]))
- warp2, warpI2 = affine_warping(solver.V.mesh(),
- np.array([[1, 0], [0, 2. - 2. * L]]))
- warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
- warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
-
- approx.plotApprox(mutar, [warp, warpI], name = 'u_app',
- homogeneized = False, what = "REAL")
- approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
- homogeneized = False, what = "REAL")
- approx.plotErr(mutar, [warp, warpI], name = 'err',
- homogeneized = False, what = "REAL")
-# approx.plotRes(mutar, [warp, warpI], name = 'res',
-# homogeneized = False, what = "REAL")
- appErr = approx.normErr(mutar, homogeneized = homogeneize)
- solNorm = approx.normHF(mutar, homogeneized = homogeneize)
- resNorm = approx.normRes(mutar, homogeneized = homogeneize)
- RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
+ approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
+ what = "REAL")
+ approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
+ approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
+# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
+ appErr = approx.normErr(mutar)
+ solNorm = approx.normHF(mutar)
+ resNorm = approx.normRes(mutar)
+ RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
- 200, [2., 1.], clip = clip)
+ 200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normErr = plotInfSet2(murange, murangeEff,
- approx, mu0, 50,
- [2., 1.], clip = clip)
+ approx, mu0, 40,
+ [2., 2.], clip = clip)
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index b494062..64f053f 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,177 +1,181 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 100
+size = 7
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
-MN = 1
+MN = 20
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
-sampling = "random"
+#sampling = "random"
if size == 1: # below
mu0 = [40 ** .5, .4]
mutar = [45 ** .5, .4]
murange = [[30 ** .5, .3], [50 ** .5, .5]]
elif size == 2: # top
mu0 = [40 ** .5, .6]
mutar = [45 ** .5, .6]
murange = [[30 ** .5, .5], [50 ** .5, .7]]
elif size == 3: # interesting
mu0 = [40 ** .5, .5]
mutar = [45 ** .5, .5]
murange = [[30 ** .5, .3], [50 ** .5, .7]]
elif size == 4: # wide_low
mu0 = [40 ** .5, .2]
mutar = [45 ** .5, .2]
murange = [[10 ** .5, .1], [70 ** .5, .3]]
elif size == 5: # wide_hi
mu0 = [40 ** .5, .8]
mutar = [45 ** .5, .8]
murange = [[10 ** .5, .7], [70 ** .5, .9]]
elif size == 6: # top_zoom
mu0 = [50 ** .5, .8]
mutar = [55 ** .5, .8]
murange = [[40 ** .5, .7], [60 ** .5, .9]]
+elif size == 7: # huge
+ mu0 = [50 ** .5, .5]
+ mutar = [55 ** .5, .5]
+ murange = [[10 ** .5, .2], [90 ** .5, .8]]
elif size == 100: # tiny
mu0 = [32.5 ** .5, .5]
mutar = [34 ** .5, .5]
murange = [[30 ** .5, .3], [35 ** .5, .7]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
H = 1.
L = .75
delta = .05
n = 25
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.), lambda x: x]
rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = mutar[1]
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normErr = plotInfSet2(murange, murangeEff,
approx, mu0, 50,
[2., 1.], clip = clip)
diff --git a/rrompy/hfengines/linear_problem/bidimensional/__init__.py b/rrompy/hfengines/linear_problem/bidimensional/__init__.py
index f972337..abf4add 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/__init__.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/__init__.py
@@ -1,34 +1,36 @@
# 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 .
#
+from .cookie_engine_single import CookieEngineSingle
from .membrane_fracture_engine import MembraneFractureEngine
from .laplace_disk_gaussian_2 import LaplaceDiskGaussian2
from .helmholtz_square_simplified_domain_problem_engine import \
HelmholtzSquareSimplifiedDomainProblemEngine
from .helmholtz_square_domain_problem_engine import \
HelmholtzSquareDomainProblemEngine
__all__ = [
+ 'CookieEngineSingle',
'MembraneFractureEngine',
'LaplaceDiskGaussian2',
'HelmholtzSquareSimplifiedDomainProblemEngine',
'HelmholtzSquareDomainProblemEngine'
]
diff --git a/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py b/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
new file mode 100644
index 0000000..5afb153
--- /dev/null
+++ b/rrompy/hfengines/linear_problem/bidimensional/cookie_engine_single.py
@@ -0,0 +1,167 @@
+# 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 .
+#
+
+import numpy as np
+import scipy.sparse as scsp
+import fenics as fen
+import ufl
+from rrompy.utilities.base.types import ScOp, List, paramVal
+from rrompy.solver.fenics import fenONE, fenZERO
+from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
+ HelmholtzProblemEngine)
+from rrompy.utilities.base import verbosityDepth
+from rrompy.utilities.poly_fitting.polynomial import (
+ hashDerivativeToIdx as hashD)
+
+__all__ = ['CookieEngineSingle']
+
+class CookieEngineSingle(HelmholtzProblemEngine):
+
+ def __init__(self, kappa:float, theta:float, n:int, R : int = 1.,
+ L : int = 2., nX : int = 1, nY : int = 1,
+ mu0 : paramVal = [12. ** .5, 1.],
+ 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 = 5
+ self.npar = 2
+ self.rescalingExp = [2., 2.]
+
+ mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(L * nX, L * nY),
+ n * nX, n * nY)
+ self.V = fen.FunctionSpace(mesh, "P", 1)
+
+ x, y = fen.SpatialCoordinate(mesh)[:]
+ cxs = np.linspace(0, L * nX, 2 * nX + 1)[1::2]
+ cys = np.linspace(0, L * nY, 2 * nY + 1)[1::2]
+ self.cookieIn = fenZERO
+ for cx in cxs:
+ for cy in cys:
+ self.cookieIn += ufl.conditional(
+ ufl.le((x-cx)**2. + (y-cy)**2., R**2.),
+ fenONE, fenZERO)
+ self.cookieOut = fenONE - self.cookieIn
+ c, s = np.cos(theta), np.sin(theta)
+ self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
+ fen.sin(kappa * (c * x + s * y))]
+
+ 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()
+ if derI <= 0 and self.As[0] is None:
+ if self.verbosity >= 20:
+ verbosityDepth("INIT", "Assembling operator term A0.",
+ timestamp = self.timestamp)
+ DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
+ self.DirichletBoundary)
+ aRe, aIm = self.diffusivity
+ hRe, hIm = self.RobinDatumH
+ termNames = ["diffusivity", "RobinDatumH"]
+ parsRe = self.iterReduceQuadratureDegree(zip(
+ [aRe, hRe],
+ [x + "Real" for x in termNames]))
+ parsIm = self.iterReduceQuadratureDegree(zip(
+ [aIm, hIm],
+ [x + "Imag" for x in termNames]))
+ a0Re = (aRe * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ + hRe * fen.dot(self.u, self.v) * self.ds(1))
+ a0Im = (aIm * fen.dot(fen.grad(self.u), fen.grad(self.v)) * fen.dx
+ + hIm * fen.dot(self.u, self.v) * self.ds(1))
+ A0Re = fen.assemble(a0Re, form_compiler_parameters = parsRe)
+ A0Im = fen.assemble(a0Im, form_compiler_parameters = parsIm)
+ DirichletBC0.apply(A0Re)
+ DirichletBC0.zero(A0Im)
+ A0ReMat = fen.as_backend_type(A0Re).mat()
+ A0ImMat = fen.as_backend_type(A0Im).mat()
+ A0Rer, A0Rec, A0Rev = A0ReMat.getValuesCSR()
+ A0Imr, A0Imc, A0Imv = A0ImMat.getValuesCSR()
+ self.As[0] = (scsp.csr_matrix((A0Rev, A0Rec, A0Rer),
+ shape = A0ReMat.size)
+ + 1.j * scsp.csr_matrix((A0Imv, A0Imc, A0Imr),
+ shape = A0ImMat.size))
+ if self.verbosity >= 20:
+ verbosityDepth("DEL", "Done assembling operator term.",
+ timestamp = self.timestamp)
+ if derI <= 1 and self.As[1] is None:
+ if self.verbosity >= 20:
+ verbosityDepth("INIT", "Assembling operator term A1.",
+ timestamp = self.timestamp)
+ 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"]))
+ a1Re = - n2Re * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
+ a1Im = - n2Im * self.cookieOut * fen.dot(self.u, self.v) * fen.dx
+ A1Re = fen.assemble(a1Re, form_compiler_parameters = parsRe)
+ A1Im = fen.assemble(a1Im, form_compiler_parameters = parsIm)
+ DirichletBC0.zero(A1Re)
+ DirichletBC0.zero(A1Im)
+ A1ReMat = fen.as_backend_type(A1Re).mat()
+ A1ImMat = fen.as_backend_type(A1Im).mat()
+ A1Rer, A1Rec, A1Rev = A1ReMat.getValuesCSR()
+ A1Imr, A1Imc, A1Imv = A1ImMat.getValuesCSR()
+ self.As[1] = (scsp.csr_matrix((A1Rev, A1Rec, A1Rer),
+ shape = A1ReMat.size)
+ + 1.j * scsp.csr_matrix((A1Imv, A1Imc, A1Imr),
+ shape = A1ImMat.size))
+ if self.verbosity >= 20:
+ verbosityDepth("DEL", "Done assembling operator term.",
+ timestamp = self.timestamp)
+ if derI <= 2 and self.As[2] is None:
+ self.As[2] = self.checkAInBounds(-1)
+ if derI <= 3 and self.As[3] is None:
+ self.As[3] = self.checkAInBounds(-1)
+ if derI <= 4 and self.As[4] is None:
+ if self.verbosity >= 20:
+ verbosityDepth("INIT", "Assembling operator term A4.",
+ timestamp = self.timestamp)
+ 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"]))
+ a4Re = - n2Re * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
+ a4Im = - n2Im * self.cookieIn * fen.dot(self.u, self.v) * fen.dx
+ A4Re = fen.assemble(a4Re, form_compiler_parameters = parsRe)
+ A4Im = fen.assemble(a4Im, form_compiler_parameters = parsIm)
+ DirichletBC0.zero(A4Re)
+ DirichletBC0.zero(A4Im)
+ A4ReMat = fen.as_backend_type(A4Re).mat()
+ A4ImMat = fen.as_backend_type(A4Im).mat()
+ A4Rer, A4Rec, A4Rev = A4ReMat.getValuesCSR()
+ A4Imr, A4Imc, A4Imv = A4ImMat.getValuesCSR()
+ self.As[4] = (scsp.csr_matrix((A4Rev, A4Rec, A4Rer),
+ shape = A4ReMat.size)
+ + 1.j * scsp.csr_matrix((A4Imv, A4Imc, A4Imr),
+ shape = A4ImMat.size))
+ if self.verbosity >= 20:
+ verbosityDepth("DEL", "Done assembling operator term.",
+ timestamp = self.timestamp)
+ return self._assembleA(mu, der, derI)
+