diff --git a/examples/2d/base/fracture.py b/examples/2d/base/fracture.py
index dac585f..a03d6b9 100644
--- a/examples/2d/base/fracture.py
+++ b/examples/2d/base/fracture.py
@@ -1,18 +1,36 @@
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
verb = 100
mu0 = [45. ** .5, .6]
H = 1.
L = .75
delta = .05
-n = 500
+n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
uh = solver.solve(mu0)[0]
solver.plotmesh(figsize = (7.5, 4.5))
print(solver.norm(uh))
solver.plot(uh, what = 'REAL', figsize = (8, 5))
-solver.plot(solver.residual(uh, mu0)[0], 'res',
+solver.plot(solver.residual(uh, mu0)[0], name = 'res',
what = 'REAL', figsize = (8, 5))
+##
+import numpy as np
+import ufl
+import fenics as fen
+from rrompy.solver.fenics import affine_warping
+L = mu0[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)
+
+
+solver.plotmesh([warp, warpI], figsize = (7.5, 4.5))
+solver.plot(uh, [warp, warpI], what = 'REAL', figsize = (8, 5))
+solver.plot(solver.residual(uh, mu0)[0], [warp, warpI], name = 'res',
+ what = 'REAL', figsize = (8, 5))
diff --git a/rrompy/hfengines/base/problem_engine_base.py b/rrompy/hfengines/base/problem_engine_base.py
index fa33798..5aa9a4f 100644
--- a/rrompy/hfengines/base/problem_engine_base.py
+++ b/rrompy/hfengines/base/problem_engine_base.py
@@ -1,330 +1,333 @@
# 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 os import path, mkdir
import fenics as fen
import numpy as np
from matplotlib import pyplot as plt
from copy import deepcopy as copy
from rrompy.utilities.base.types import (Np1D, strLst, FenFunc, Tuple, List,
paramVal)
from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
-from rrompy.solver.fenics import L2NormMatrix
+from rrompy.solver.fenics import L2NormMatrix, fenplot
from .boundary_conditions import BoundaryConditions
from .matrix_engine_base import MatrixEngineBase
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['ProblemEngineBase']
class ProblemEngineBase(MatrixEngineBase):
"""
Generic solver for parametric problems.
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
"""
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(verbosity = verbosity, timestamp = timestamp)
self.BCManager = BoundaryConditions("Dirichlet")
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.mu0BC = np.nan
self.degree_threshold = degree_threshold
self.npar = 0
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
self.resetAs()
self.resetbs()
if not type(V).__name__ == 'FunctionSpace':
raise RROMPyException("V type not recognized.")
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
def spacedim(self):
return self.V.dim()
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = L2NormMatrix(self.V)
def liftDirichletData(self, mu : paramVal = []) -> Np1D:
"""Lift Dirichlet datum."""
mu = self.checkParameter(mu)
if mu != self.mu0BC:
self.mu0BC = copy(mu)
try:
liftRe = fen.interpolate(self.DirichletDatum[0], self.V)
except:
liftRe = fen.project(self.DirichletDatum[0], self.V)
try:
liftIm = fen.interpolate(self.DirichletDatum[1], self.V)
except:
liftIm = fen.project(self.DirichletDatum[1], self.V)
self.liftedDirichletDatum = (np.array(liftRe.vector())
+ 1.j * np.array(liftIm.vector()))
return self.liftedDirichletDatum
def reduceQuadratureDegree(self, fun:FenFunc, name:str):
"""Check whether to reduce compiler parameters to degree threshold."""
if not np.isinf(self.degree_threshold):
from ufl.algorithms.estimate_degrees import (
estimate_total_polynomial_degree as ETPD)
try:
deg = ETPD(fun)
except:
return False
if deg > self.degree_threshold:
if self.verbosity >= 15:
verbosityDepth("MAIN", ("Reducing quadrature degree from "
"{} to {} for {}.").format(
deg,
self.degree_threshold,
name),
timestamp = self.timestamp)
return True
return False
def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]):
"""
Iterate reduceQuadratureDegree over list and define reduce compiler
parameters.
"""
if funsNames is not None:
for fun, name in funsNames:
if self.reduceQuadratureDegree(fun, name):
return {"quadrature_degree" : self.degree_threshold}
return {}
- def plot(self, u:Np1D, name : str = "u", save : str = None,
- what : strLst = 'all', saveFormat : str = "eps",
- saveDPI : int = 100, show : bool = True, **figspecs):
+ def plot(self, u:Np1D, warping : List[callable] = None, name : str = "u",
+ save : str = None, what : strLst = 'all',
+ saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
+ **figspecs):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
+ warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = (13. * len(what) / 4, 3)
subplotcode = 100 + len(what) * 10
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
uAb = fen.Function(self.V)
uAb.vector().set_local(np.abs(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uAb, title = "|{0}|".format(name))
+ p = fenplot(uAb, warping = warping, title = "|{0}|".format(name))
plt.colorbar(p)
if 'PHASE' in what:
uPh = fen.Function(self.V)
uPh.vector().set_local(np.angle(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uPh, title = "phase({0})".format(name))
+ p = fenplot(uPh, warping = warping,
+ title = "phase({0})".format(name))
plt.colorbar(p)
if 'REAL' in what:
uRe = fen.Function(self.V)
uRe.vector().set_local(np.real(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uRe, title = "Re({0})".format(name))
+ p = fenplot(uRe, warping = warping, title = "Re({0})".format(name))
plt.colorbar(p)
if 'IMAG' in what:
uIm = fen.Function(self.V)
uIm.vector().set_local(np.imag(u))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uIm, title = "Im({0})".format(name))
+ p = fenplot(uIm, warping = warping, title = "Im({0})".format(name))
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
- def plotmesh(self, name : str = "Mesh", save : str = None,
- saveFormat : str = "eps", saveDPI : int = 100,
- show : bool = True, **figspecs):
+ def plotmesh(self, warping : List[callable] = None, name : str = "Mesh",
+ save : str = None, saveFormat : str = "eps",
+ saveDPI : int = 100, show : bool = True, **figspecs):
"""
Do a nice plot of the mesh.
Args:
u: numpy complex array with function dofs.
+ warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
plt.figure(**figspecs)
- fen.plot(self.V.mesh())
+ fenplot(self.V.mesh(), warping = warping)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
def outParaview(self, u:Np1D, name : str = "u", filename : str = "out",
time : float = 0., what : strLst = 'all',
forceNewFile : bool = True, folder : bool = False,
filePW = None):
"""
Output complex-valued function with given dofs to ParaView file.
Args:
u: numpy complex array with function dofs.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
time(optional): Timestamp.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
filePW(optional): Fenics File entity (for time series).
"""
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if filePW is None:
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
if what == ['MESH']:
filePW << (self.V.mesh(), time)
if 'ABS' in what:
uAb = fen.Function(self.V, name = "{}_ABS".format(name))
uAb.vector().set_local(np.abs(u))
filePW << (uAb, time)
if 'PHASE' in what:
uPh = fen.Function(self.V, name = "{}_PHASE".format(name))
uPh.vector().set_local(np.angle(u))
filePW << (uPh, time)
if 'REAL' in what:
uRe = fen.Function(self.V, name = "{}_REAL".format(name))
uRe.vector().set_local(np.real(u))
filePW << (uRe, time)
if 'IMAG' in what:
uIm = fen.Function(self.V, name = "{}_IMAG".format(name))
uIm.vector().set_local(np.imag(u))
filePW << (uIm, time)
return filePW
def outParaviewTimeDomain(self, u:Np1D, omega:float,
timeFinal : float = None,
periodResolution : int = 20, name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folder : bool = False):
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
Args:
u: numpy complex array with function dofs.
omega: frequency.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
folder(optional): Whether to create an additional folder layer.
"""
if folder:
if not path.exists(filename + "/"):
mkdir(filename)
idxpath = filename.rfind("/")
filename += "/" + filename[idxpath + 1 :]
if forceNewFile:
filePW = fen.File(getNewFilename(filename, "pvd"))
else:
filePW = fen.File("{}.pvd".format(filename))
omega = np.abs(omega)
t = 0.
dt = 2. * np.pi / omega / periodResolution
if timeFinal is None: timeFinal = 2. * np.pi / omega - dt
for j in range(int(np.ceil(timeFinal / dt)) + 1):
ut = fen.Function(self.V, name = name)
ut.vector().set_local(np.real(u) * np.cos(omega * t)
+ np.imag(u) * np.sin(omega * t))
filePW << (ut, t)
t += dt
return filePW
-
diff --git a/rrompy/hfengines/base/vector_problem_engine_base.py b/rrompy/hfengines/base/vector_problem_engine_base.py
index 8fba867..0e41e6a 100644
--- a/rrompy/hfengines/base/vector_problem_engine_base.py
+++ b/rrompy/hfengines/base/vector_problem_engine_base.py
@@ -1,197 +1,206 @@
# 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 fenics as fen
import numpy as np
from matplotlib import pyplot as plt
-from rrompy.utilities.base.types import Np1D, strLst
+from rrompy.utilities.base.types import Np1D, List, strLst
from rrompy.utilities.base import purgeList, getNewFilename
+from rrompy.solver.fenics import fenplot
from .problem_engine_base import ProblemEngineBase
__all__ = ['VectorProblemEngineBase']
class VectorProblemEngineBase(ProblemEngineBase):
"""
Generic solver for parametric vector problems.
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real FE space.
u: Generic trial functions for variational form evaluation.
v: Generic test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
liftedDirichletDatum: Dofs of Dirichlet datum lifting.
mu0BC: Mu value of last Dirichlet datum lifting.
degree_threshold: Threshold for ufl expression interpolation degree.
"""
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self.V = fen.VectorFunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.npar = 0
- def plot(self, u:Np1D, name : str = "u", save : str = None,
- what : strLst = 'all', saveFormat : str = "eps",
- saveDPI : int = 100, show : bool = True, **figspecs):
+ def plot(self, u:Np1D, warping : List[callable] = None, name : str = "u",
+ save : str = None, what : strLst = 'all',
+ saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
+ **figspecs):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
+ warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = (13. * max(len(what), 1) / 4, 3)
if len(what) > 0:
for j in range(self.V.num_sub_spaces()):
subplotcode = 100 + len(what) * 10
II = self.V.sub(j).dofmap().dofs()
Vj = self.V.sub(j).collapse()
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
uAb = fen.Function(Vj)
uAb.vector().set_local(np.abs(u[II]))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uAb, title = "|{}_comp{}|".format(name, j))
+ p = fenplot(uAb, warping = warping,
+ title = "|{}_comp{}|".format(name, j))
plt.colorbar(p)
if 'PHASE' in what:
uPh = fen.Function(Vj)
uPh.vector().set_local(np.angle(u[II]))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uPh, title = "phase({}_comp{})".format(name,
- j))
+ p = fenplot(uPh, warping = warping,
+ title = "phase({}_comp{})".format(name, j))
plt.colorbar(p)
if 'REAL' in what:
uRe = fen.Function(Vj)
uRe.vector().set_local(np.real(u[II]))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uRe, title = "Re({}_comp{})".format(name, j))
+ p = fenplot(uRe, warping = warping,
+ title = "Re({}_comp{})".format(name, j))
plt.colorbar(p)
if 'IMAG' in what:
uIm = fen.Function(Vj)
uIm.vector().set_local(np.imag(u[II]))
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
- p = fen.plot(uIm, title = "Im({}_comp{})".format(name, j))
+ p = fenplot(uIm, warping = warping,
+ title = "Im({}_comp{})".format(name, j))
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_comp{}_fig_".format(save, j),
saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
try:
if len(what) > 1:
figspecs['figsize'] = (2. / len(what) * figspecs['figsize'][0],
figspecs['figsize'][1])
elif len(what) == 0:
figspecs['figsize'] = (2. * figspecs['figsize'][0],
figspecs['figsize'][1])
if len(what) == 0 or 'ABS' in what or 'REAL' in what:
uVRe = fen.Function(self.V)
uVRe.vector().set_local(np.real(u))
plt.figure(**figspecs)
plt.jet()
- p = fen.plot(uVRe, title = "{}_Re".format(name),
- mode = "displacement")
+ p = fenplot(uVRe, warping = warping,
+ title = "{}_Re".format(name),
+ mode = "displacement")
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_disp_Re_fig_".format(save),
saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
if 'ABS' in what or 'IMAG' in what:
uVIm = fen.Function(self.V)
uVIm.vector().set_local(np.imag(u))
plt.figure(**figspecs)
plt.jet()
- p = fen.plot(uVIm, title = "{}_Im".format(name),
- mode = "displacement")
+ p = fenplot(uVIm, warping = warping,
+ title = "{}_Im".format(name),
+ mode = "displacement")
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_disp_Im_fig_".format(save, j),
saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
except:
pass
- def plotmesh(self, name : str = "Mesh", save : str = None,
- saveFormat : str = "eps", saveDPI : int = 100,
- show : bool = True, **figspecs):
+ def plotmesh(self, warping : List[callable] = None, name : str = "Mesh",
+ save : str = None, saveFormat : str = "eps",
+ saveDPI : int = 100, show : bool = True, **figspecs):
"""
Do a nice plot of the mesh.
Args:
u: numpy complex array with function dofs.
+ warping(optional): Domain warping functions.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
plt.figure(**figspecs)
- fen.plot(self.V.mesh())
+ fenplot(self.V.mesh(), warping = warping)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
if show:
plt.show()
plt.close()
diff --git a/rrompy/sampling/base/pod_engine.py b/rrompy/sampling/base/pod_engine.py
index 1df7cf3..82a8699 100644
--- a/rrompy/sampling/base/pod_engine.py
+++ b/rrompy/sampling/base/pod_engine.py
@@ -1,149 +1,124 @@
# 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
from copy import deepcopy as copy
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, HFEng, sampList
from rrompy.sampling import sampleList
-from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['PODEngine']
class PODEngine:
"""
POD engine for general matrix orthogonalization.
"""
def __init__(self, HFEngine:HFEng):
self.HFEngine = HFEngine
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
- def GS(self, a:Np1D, Q:sampList, n : int = None,
- aA:Np1D = None, QA:sampList = None) -> Tuple[Np1D, Np1D, Np1D]:
+ def GS(self, a:Np1D, Q:sampList, n : int = -1) -> Tuple[Np1D, Np1D, bool]:
"""
Compute 1 Gram-Schmidt step with given projector.
Args:
a: vector to be projected;
Q: orthogonal projection matrix;
n: number of columns of Q to be considered;
- aA: augmented components of vector to be projected;
- QA: augmented components of projection matrix.
Returns:
Resulting normalized vector, coefficients of a wrt the updated
- basis.
+ basis, whether computation is ill-conditioned.
"""
- if n is None:
+ if n == -1:
n = Q.shape[1]
- if aA is None != QA is None:
- raise RROMPyException(("Either both or none of augmented "
- "components must be provided."))
- r = np.zeros((n + 1,), dtype = a.dtype)
+ r = np.zeros((n + 1,), dtype = Q.dtype)
if n > 0:
Q = Q[: n]
for j in range(2): # twice is enough!
nu = self.HFEngine.innerProduct(a, Q)
a = a - Q.dot(nu)
- if aA is not None:
- aA = aA - QA.dot(nu)
r[:-1] = r[:-1] + nu.flatten()
r[-1] = self.HFEngine.norm(a)
- if np.isclose(np.abs(r[-1]), 0.):
+ ill_cond = False
+ if np.isclose(np.abs(r[-1]) / np.linalg.norm(r), 0.):
+ ill_cond = True
r[-1] = 1.
a = a / r[-1]
- if aA is not None:
- aA = aA / r[-1]
- return a, r, aA
+ return a, r, ill_cond
- def QRGramSchmidt(self, A:sampList,
+ def generalizedQR(self, A:sampList, Q0 : sampList = None,
only_R : bool = False) -> Tuple[sampList, Np2D]:
"""
- Compute QR decomposition of a matrix through Gram-Schmidt method.
-
- Args:
- A: matrix to be decomposed;
- only_R(optional): whether to skip reconstruction of Q; defaults to
- False.
-
- Returns:
- Resulting orthogonal and upper-triangular factors.
- """
- N = A.shape[1]
- Q = copy(A)
- R = np.zeros((N, N), dtype = A.dtype)
- for k in range(N):
- Q[k], R[: k + 1, k], _ = self.GS(A[k], Q, k)
- if only_R:
- return R
- return Q, R
-
- def QRHouseholder(self, A:sampList, Q0 : sampList = None,
- only_R : bool = False) -> Tuple[sampList, Np2D]:
- """
- Compute QR decomposition of a matrix through Householder method.
+ Compute generalized QR decomposition of a matrix through Householder
+ method.
Args:
A: matrix to be decomposed;
Q0(optional): initial orthogonal guess for Q; defaults to random;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting (orthogonal and )upper-triangular factor(s).
"""
- N = A.shape[1]
+ Nh, N = A.shape
B = copy(A)
V = copy(A)
R = np.zeros((N, N), dtype = A.dtype)
if Q0 is None:
Q = sampleList(np.zeros(A.shape, dtype = A.dtype)
+ np.random.randn(*(A.shape)))
else:
Q = copy(Q0)
for k in range(N):
- if Q0 is None:
- Q[k], _, _ = self.GS(Q[k], Q, k)
+ if k <= Nh:
+ if Q0 is None:
+ illC = True
+ while illC:
+ Q[k], _, illC = self.GS(np.random.randn(Nh), Q, k)
+ else:
+ Q[k] = np.zeros(Nh, dtype = Q.dtype)
a = B[k]
R[k, k] = self.HFEngine.norm(a)
alpha = self.HFEngine.innerProduct(a, Q[k])
if np.isclose(np.abs(alpha), 0.): s = 1.
else: s = - alpha / np.abs(alpha)
Q[k] = s * Q[k]
V[k], _, _ = self.GS(R[k, k] * Q[k] - a, Q, k)
J = np.arange(k + 1, N)
vtB = self.HFEngine.innerProduct(B[J], V[k])
B.data[:, J] -= 2 * np.outer(V[k], vtB)
R[k, J] = self.HFEngine.innerProduct(B[J], Q[k])
B.data[:, J] -= np.outer(Q[k], R[k, J])
if only_R:
return R
for k in range(N - 1, -1, -1):
J = list(range(k, N))
vtQ = self.HFEngine.innerProduct(Q[J], V[k])
Q.data[:, J] -= 2 * np.outer(V[k], vtQ)
return Q, R
diff --git a/rrompy/sampling/linear_problem/__init__.py b/rrompy/sampling/linear_problem/__init__.py
index 397644d..8326efd 100644
--- a/rrompy/sampling/linear_problem/__init__.py
+++ b/rrompy/sampling/linear_problem/__init__.py
@@ -1,29 +1,27 @@
# 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 .sampling_engine_linear import SamplingEngineLinear
from .sampling_engine_linear_pod import SamplingEngineLinearPOD
-from .sampling_engine_linear_pod_naive import SamplingEngineLinearPODNaive
__all__ = [
'SamplingEngineLinear',
- 'SamplingEngineLinearPOD',
- 'SamplingEngineLinearPODNaive'
+ 'SamplingEngineLinearPOD'
]
diff --git a/rrompy/sampling/linear_problem/sampling_engine_linear.py b/rrompy/sampling/linear_problem/sampling_engine_linear.py
index d7a6412..cdd7246 100644
--- a/rrompy/sampling/linear_problem/sampling_engine_linear.py
+++ b/rrompy/sampling/linear_problem/sampling_engine_linear.py
@@ -1,123 +1,118 @@
# 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 copy import deepcopy as copy
import numpy as np
from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityDepth
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.poly_fitting.polynomial import nextDerivativeIndices
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList
__all__ = ['SamplingEngineLinear']
class SamplingEngineLinear(SamplingEngineBase):
"""HERE"""
def preprocesssamples(self, idxs:Np1D) -> sampList:
if self.samples is None or len(self.samples) == 0: return
return self.samples(idxs)
def postprocessu(self, u:sampList, overwrite : bool = False) -> Np1D:
return copy(u)
def postprocessuBulk(self, u:sampList) -> sampList:
return copy(u)
def lastSampleManagement(self):
pass
def _getSampleConcurrence(self, mu:paramVal, previous:Np1D,
homogeneized : bool = False) -> sampList:
if len(previous) >= len(self._derIdxs):
self._derIdxs += nextDerivativeIndices(self._derIdxs,
self.HFEngine.npar,
len(previous) + 1 - len(self._derIdxs))
derIdx = self._derIdxs[len(previous)]
mu = checkParameter(mu, self.HFEngine.npar)
samplesOld = self.preprocesssamples(previous)
RHS = self.HFEngine.b(mu, derIdx, homogeneized = homogeneized)
for j, derP in enumerate(self._derIdxs[: len(previous)]):
diffP = [x - y for (x, y) in zip(derIdx, derP)]
if np.all([x >= 0 for x in diffP]):
RHS -= self.HFEngine.A(mu, diffP).dot(samplesOld[j])
return self.solveLS(mu, RHS = RHS, homogeneized = homogeneized)
def nextSample(self, mu : paramVal = [], overwrite : bool = False,
homogeneized : bool = False,
lastSample : bool = True) -> Np1D:
mu = checkParameter(mu, self.HFEngine.npar)
ns = self.nsamples
muidxs = self.mus.findall(mu[0])
if len(muidxs) > 0:
u = self._getSampleConcurrence(mu, np.sort(muidxs), homogeneized)
else:
u = self.solveLS(mu, homogeneized = homogeneized)
u = self.postprocessu(u, overwrite = overwrite)
if overwrite:
self.samples[ns] = u
self.mus[ns] = mu[0]
else:
if ns == 0:
self.samples = sampleList(u)
else:
self.samples.append(u)
self.mus.append(mu)
self.nsamples += 1
if lastSample: self.lastSampleManagement()
return u
def iterSample(self, mus:paramList,
homogeneized : bool = False) -> sampList:
mus, _ = checkParameterList(mus, self.HFEngine.npar)
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting sampling iterations.",
timestamp = self.timestamp)
n = len(mus)
if n <= 0:
raise RROMPyException(("Number of samples must be positive."))
self.resetHistory()
if self.allowRepeatedSamples:
- if self.verbosity >= 7:
- verbosityDepth("MAIN", "Computing sample {}/{}.".format(1, n),
- timestamp = self.timestamp)
- u = self.nextSample(mus[0], homogeneized = homogeneized,
- lastSample = (n == 1))
- if n > 1:
- self.preallocateSamples(u, mus[0], n)
- for j in range(1, n):
- if self.verbosity >= 7:
- verbosityDepth("MAIN", ("Computing sample "
- "{}/{}.").format(j + 1, n),
- timestamp = self.timestamp)
- self.nextSample(mus[j], overwrite = True,
- homogeneized = homogeneized,
- lastSample = (n == j + 1))
+ for j in range(n):
+ if self.verbosity >= 7:
+ verbosityDepth("MAIN", ("Computing sample "
+ "{} / {}.").format(j + 1, n),
+ timestamp = self.timestamp)
+ self.nextSample(mus[j], overwrite = (j > 0),
+ homogeneized = homogeneized,
+ lastSample = (n == j + 1))
+ if j == 0:
+ self.preallocateSamples(self.samples[0], mus[0], n)
else:
self.samples = self.postprocessuBulk(self.solveLS(mus,
homogeneized = homogeneized))
self.mus = copy(mus)
- self.nsamples = len(mus)
+ self.nsamples = n
if self.verbosity >= 5:
verbosityDepth("DEL", "Finished sampling iterations.",
timestamp = self.timestamp)
return self.samples
diff --git a/rrompy/sampling/linear_problem/sampling_engine_linear_pod.py b/rrompy/sampling/linear_problem/sampling_engine_linear_pod.py
index 30b71aa..bcfea23 100644
--- a/rrompy/sampling/linear_problem/sampling_engine_linear_pod.py
+++ b/rrompy/sampling/linear_problem/sampling_engine_linear_pod.py
@@ -1,73 +1,84 @@
# 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 copy import deepcopy as copy
from rrompy.sampling.base.pod_engine import PODEngine
from .sampling_engine_linear import SamplingEngineLinear
from rrompy.utilities.base.types import Np1D, paramVal, sampList
+from rrompy.utilities.base import verbosityDepth
from rrompy.sampling import sampleList
__all__ = ['SamplingEngineLinearPOD']
class SamplingEngineLinearPOD(SamplingEngineLinear):
"""HERE"""
def resetHistory(self):
super().resetHistory()
self.samples_full = None
self.RPOD = None
def popSample(self):
if hasattr(self, "nsamples") and self.nsamples > 1:
self.RPOD = self.RPOD[: -1, : -1]
self.samples_full.pop()
super().popSample()
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
self.PODEngine = PODEngine(self._HFEngine)
def preprocesssamples(self, idxs:Np1D) -> sampList:
if self.samples_full is None or len(self.samples_full) == 0: return
return self.samples_full(idxs)
def postprocessu(self, u:sampList, overwrite : bool = False) -> Np1D:
ns = self.nsamples
if overwrite:
self.samples_full[ns] = copy(u)
else:
if ns == 0:
self.samples_full = sampleList(u)
else:
self.samples_full.append(u)
return u
+ def postprocessuBulk(self, u:sampList) -> sampList:
+ self.samples_full = copy(u)
+ if self.verbosity >= 10:
+ verbosityDepth("INIT", "Starting orthogonalization.",
+ timestamp = self.timestamp)
+ u, self.RPOD = self.PODEngine.generalizedQR(self.samples_full)
+ if self.verbosity >= 10:
+ verbosityDepth("DEL", "Done orthogonalizing.",
+ timestamp = self.timestamp)
+ return u
+
def lastSampleManagement(self):
- self.samples, self.RPOD = self.PODEngine.QRHouseholder(
- self.samples_full)
+ self.samples = self.postprocessuBulk(self.samples_full)
def preallocateSamples(self, u:Np1D, mu:paramVal, n:int):
super().preallocateSamples(u, mu, n)
self.samples_full.reset((u.shape[0], n), u.dtype)
self.samples_full[0] = u
diff --git a/rrompy/sampling/linear_problem/sampling_engine_linear_pod_naive.py b/rrompy/sampling/linear_problem/sampling_engine_linear_pod_naive.py
deleted file mode 100644
index 83dddae..0000000
--- a/rrompy/sampling/linear_problem/sampling_engine_linear_pod_naive.py
+++ /dev/null
@@ -1,93 +0,0 @@
-# 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
-from rrompy.sampling.base.pod_engine import PODEngine
-from .sampling_engine_linear import SamplingEngineLinear
-from rrompy.utilities.base.types import Np1D, paramVal, sampList
-from rrompy.utilities.base import verbosityDepth
-
-__all__ = ['SamplingEngineLinearPODNaive']
-
-class SamplingEngineLinearPODNaive(SamplingEngineLinear):
- """HERE"""
-
- def resetHistory(self):
- super().resetHistory()
- self.RPOD = None
-
- def popSample(self):
- if hasattr(self, "nsamples") and self.nsamples > 1:
- self.RPOD = self.RPOD[: -1, : -1]
- super().popSample()
-
- @property
- def HFEngine(self):
- """Value of HFEngine. Its assignment resets history."""
- return self._HFEngine
- @HFEngine.setter
- def HFEngine(self, HFEngine):
- self._HFEngine = HFEngine
- self.resetHistory()
- self.PODEngine = PODEngine(self._HFEngine)
-
- def preprocesssamples(self, idxs:Np1D) -> sampList:
- idxMax = np.max(idxs) + 1
- sampleBase = super().preprocesssamples(np.arange(idxMax))
- RPODBase = self.RPOD[: idxMax, idxs]
- return sampleBase.dot(RPODBase)
-
- def postprocessu(self, u:sampList, overwrite : bool = False) -> Np1D:
- if self.verbosity >= 10:
- verbosityDepth("INIT", "Starting orthogonalization.",
- timestamp = self.timestamp)
- ns = self.nsamples
- if ns == 0:
- u, r, _ = self.PODEngine.GS(u, np.empty((0, 0)))
- r = r[0]
- else:
- u, r, _ = self.PODEngine.GS(u, self.samples(np.arange(ns)), ns)
- if overwrite:
- self.RPOD[: ns + 1, ns] = r
- else:
- if ns == 0:
- self.RPOD = r.reshape((1, 1))
- else:
- self.RPOD=np.block([[ self.RPOD, r[:-1, None]],
- [np.zeros((1, ns)), r[-1]]])
- if self.verbosity >= 10:
- verbosityDepth("DEL", "Done orthogonalizing.",
- timestamp = self.timestamp)
- return u
-
- def postprocessuBulk(self, u:sampList) -> sampList:
- if self.verbosity >= 10:
- verbosityDepth("INIT", "Starting orthogonalization.",
- timestamp = self.timestamp)
- u, self.RPOD = self.PODEngine.QRHouseholder(u)
- if self.verbosity >= 10:
- verbosityDepth("DEL", "Done orthogonalizing.",
- timestamp = self.timestamp)
- return u
-
- def preallocateSamples(self, u:Np1D, mu:paramVal, n:int):
- super().preallocateSamples(u, mu, n)
- r = self.RPOD
- self.RPOD = np.zeros((n, n), dtype = u.dtype)
- self.RPOD[0, 0] = r[0, 0]
-
diff --git a/rrompy/solver/fenics/__init__.py b/rrompy/solver/fenics/__init__.py
index 1e1a506..93cca6a 100644
--- a/rrompy/solver/fenics/__init__.py
+++ b/rrompy/solver/fenics/__init__.py
@@ -1,38 +1,41 @@
# 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 .fenics_constants import (fenZERO, fenZEROS, fenONE, fenONES,
bdrTrue, bdrFalse)
from .fenics_norms import (L2NormMatrix, H1NormMatrix, Hminus1NormMatrix,
elasticNormMatrix, elasticDualNormMatrix)
+from .fenics_plotting import fenplot, affine_warping
__all__ = [
'fenZERO',
'fenZEROS',
'fenONE',
'fenONES',
'bdrTrue',
'bdrFalse',
'L2NormMatrix',
'H1NormMatrix',
'Hminus1NormMatrix',
'elasticNormMatrix',
- 'elasticDualNormMatrix'
+ 'elasticDualNormMatrix',
+ 'fenplot',
+ 'affine_warping'
]
diff --git a/rrompy/solver/fenics/fenics_plotting.py b/rrompy/solver/fenics/fenics_plotting.py
new file mode 100644
index 0000000..eedbe3e
--- /dev/null
+++ b/rrompy/solver/fenics/fenics_plotting.py
@@ -0,0 +1,91 @@
+# 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 dolfin.cpp as cpp
+import ufl
+import fenics as fen
+from rrompy.utilities.base.types import Np1D, Np2D
+
+_all_plottable_types = (cpp.function.Function, cpp.function.Expression,
+ cpp.mesh.Mesh, cpp.fem.DirichletBC,
+ cpp.mesh.MeshFunctionBool, cpp.mesh.MeshFunctionInt,
+ cpp.mesh.MeshFunctionDouble,
+ cpp.mesh.MeshFunctionSizet)
+
+__all__ = ['fenplot', 'affine_warping']
+
+def fenplot(object, *args, warping = None, **kwargs):
+ "See dolfin.common.plot for more details."
+ mesh = kwargs.pop('mesh', None)
+ if isinstance(object, cpp.mesh.Mesh):
+ if mesh is not None and mesh.id() != object.id():
+ raise RuntimeError("Got different mesh in plot object and keyword "
+ "argument")
+ mesh = object
+ if mesh is None:
+ if isinstance(object, cpp.function.Function):
+ mesh = object.function_space().mesh()
+ elif hasattr(object, "mesh"):
+ mesh = object.mesh()
+ if not isinstance(object, _all_plottable_types):
+ from dolfin.fem.projection import project
+ try:
+ cpp.log.info("Object cannot be plotted directly, projecting to "
+ "piecewise linears.")
+ object = project(object, mesh=mesh)
+ mesh = object.function_space().mesh()
+ object = object._cpp_object
+ except Exception as e:
+ msg = "Don't know how to plot given object:\n %s\n" \
+ "and projection failed:\n %s" % (str(object), str(e))
+ raise RuntimeError(msg)
+
+ if warping is not None:
+ try:
+ disp = fen.interpolate(warping[0],
+ fen.VectorFunctionSpace(mesh, "Lagrange", 1))
+ except:
+ disp = fen.project(warping[0],
+ fen.VectorFunctionSpace(mesh, "Lagrange", 1))
+ fen.ALE.move(mesh, disp)
+ out = fen.plot(object, *args, mesh = mesh, **kwargs)
+ if warping is not None:
+ try:
+ disp = fen.interpolate(warping[1],
+ fen.VectorFunctionSpace(mesh, "Lagrange", 1))
+ except:
+ disp = fen.project(warping[1],
+ fen.VectorFunctionSpace(mesh, "Lagrange", 1))
+ fen.ALE.move(mesh, disp)
+ return out
+
+def affine_warping(mesh, A:Np2D, b : Np1D = None):
+ coords = fen.SpatialCoordinate(mesh)[:]
+ ndim = mesh.topology().dim()
+ if b is None: b = [0.] * ndim
+ assert A.shape[0] == ndim and A.shape[1] == ndim and len(b) == ndim
+ Ainv = np.linalg.inv(A)
+ warp = [- 1. * c for c in coords]
+ warpInv = [- 1. * c for c in coords]
+ for i in range(ndim):
+ warp[i] = warp[i] + b[i]
+ for j in range(ndim):
+ warp[i] = warp[i] + A[i, j] * coords[j]
+ warpInv[i] = warpInv[i] + Ainv[i, j] * (coords[j] - b[j])
+ return tuple([ufl.as_vector(tuple(w)) for w in [warp, warpInv]])