diff --git a/rrompy/hfengines/base/fenics_engine_base.py b/rrompy/hfengines/base/fenics_engine_base.py
index 4000a61..c11b2c9 100644
--- a/rrompy/hfengines/base/fenics_engine_base.py
+++ b/rrompy/hfengines/base/fenics_engine_base.py
@@ -1,513 +1,515 @@
# 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/>.
#
from os import path, mkdir
import fenics as fen
import numpy as np
from matplotlib import pyplot as plt
from .numpy_engine_base import NumpyEngineBase, checknports
from rrompy.utilities.base.types import (Np1D, strLst, FenFunc, Tuple, List,
FigHandle)
from rrompy.utilities.base.data_structures import purgeList, getNewFilename
from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.solver.fenics import L2NormMatrix, fenplot, interp_project
+from rrompy.solver.fenics import (L2NormMatrix, fenplot, interp_project,
+ serializeFunctionSpace)
from .boundary_conditions import BoundaryConditions
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.parallel import (masterCore, bcast, indicesScatter,
- listGather)
+ listGather, updateSerialIndex)
__all__ = ['FenicsEngineBase', 'FenicsEngineBaseTensorized']
def plottingBaseFen(u, fig, V, what, nRows, subplotidx, warping, name,
colorbar, fenplotArgs):
if 'ABS' in what:
uAb = fen.Function(V)
uAb.vector().set_local(np.abs(u))
subplotidx = subplotidx + 1
ax = fig.add_subplot(nRows, len(what), subplotidx)
p = fenplot(uAb, warping = warping, title = "|{}|".format(name),
**fenplotArgs)
if colorbar: fig.colorbar(p, ax = ax)
if 'PHASE' in what:
uPh = fen.Function(V)
uPh.vector().set_local(np.angle(u))
subplotidx = subplotidx + 1
ax = fig.add_subplot(nRows, len(what), subplotidx)
p = fenplot(uPh, warping = warping, title = "phase({})".format(name),
**fenplotArgs)
if colorbar: fig.colorbar(p, ax = ax)
if 'REAL' in what:
uRe = fen.Function(V)
uRe.vector().set_local(np.real(u))
subplotidx = subplotidx + 1
ax = fig.add_subplot(nRows, len(what), subplotidx)
p = fenplot(uRe, warping = warping, title = "Re({})".format(name),
**fenplotArgs)
if colorbar: fig.colorbar(p, ax = ax)
if 'IMAG' in what:
uIm = fen.Function(V)
uIm.vector().set_local(np.imag(u))
subplotidx = subplotidx + 1
ax = fig.add_subplot(nRows, len(what), subplotidx)
p = fenplot(uIm, warping = warping, title = "Im({})".format(name),
**fenplotArgs)
if colorbar: fig.colorbar(p, ax = ax)
class FenicsEngineBase(NumpyEngineBase):
"""Generic solver for parametric fenics problems."""
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.degree_threshold = degree_threshold
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
if not type(V).__name__ == 'FunctionSpace':
raise RROMPyException("V type not recognized.")
self.dsToBeSet = True
+ if V.mesh().mpi_comm().size > 1: V = serializeFunctionSpace(V)
self._V = V
self.u = fen.TrialFunction(V)
self.v = fen.TestFunction(V)
@property
def spacedim(self):
if hasattr(self, "_V"): return self.V.dim()
return super().spacedim
def autoSetDS(self):
"""Set FEniCS boundary measure based on boundary function handles."""
if self.dsToBeSet:
vbMng(self, "INIT", "Initializing boundary measures.", 20)
mesh = self.V.mesh()
NB = self.NeumannBoundary
RB = self.RobinBoundary
boundary_markers = fen.MeshFunction("size_t", mesh,
mesh.topology().dim() - 1)
NB.mark(boundary_markers, 0)
RB.mark(boundary_markers, 1)
self.ds = fen.Measure("ds", domain = mesh,
subdomain_data = boundary_markers)
self.dsToBeSet = False
vbMng(self, "DEL", "Done assembling boundary measures.", 20)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self.energyNormMatrix = L2NormMatrix(self.V)
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
self.energyNormDualMatrix = self.energyNormMatrix
def liftDirichletData(self) -> Np1D:
"""Lift Dirichlet datum."""
if not hasattr(self, "_liftedDirichletDatum"):
liftRe = interp_project(self.DirichletDatum[0], self.V)
liftIm = interp_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:
vbMng(self, "MAIN",
("Reducing quadrature degree from {} to {} for "
"{}.").format(deg, self.degree_threshold, name), 15)
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, warping : List[callable] = None,
is_state : bool = False, name : str = "u", save : str = None,
what : strLst = 'all', forceNewFile : bool = True,
saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
colorMap : str = "jet", fenplotArgs : dict = {},
**figspecs) -> Tuple[FigHandle, str]:
"""
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.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
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'.
forceNewFile(optional): Whether to create new output file.
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.
colorMap(optional): Pyplot colormap. Defaults to 'jet'.
fenplotArgs(optional): Optional arguments for fenplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename and figure handle.
"""
if not is_state and not self.isCEye:
return super().plot(u, warping, False, name, save, what,
forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
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
out = None
if masterCore():
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = plt.figaspect(1. / len(what))
fig = plt.figure(**figspecs)
plt.set_cmap(colorMap)
plottingBaseFen(u, fig, self.V, what, 1, 0, warping, name,
self.V.mesh().geometric_dimension() > 1,
fenplotArgs)
plt.tight_layout()
if save is not None:
save = save.strip()
if forceNewFile:
fileOut = getNewFilename("{}_fig_".format(save),
saveFormat)
else:
fileOut = "{}_fig.{}".format(save, saveFormat)
fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show: plt.show()
out = fig if fileOut is None else (fig, fileOut)
return bcast(out)
def plotmesh(self, warping : List[callable] = None, name : str = "Mesh",
save : str = None, forceNewFile : bool = True,
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, fenplotArgs : dict = {},
**figspecs) -> Tuple[FigHandle, str]:
"""
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.
forceNewFile(optional): Whether to create new output file.
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.
fenplotArgs(optional): Optional arguments for fenplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename and figure handle.
"""
out = None
if masterCore():
fig = plt.figure(**figspecs)
fenplot(self.V.mesh(), warping = warping, **fenplotArgs)
plt.tight_layout()
if save is not None:
save = save.strip()
if forceNewFile:
fileOut = getNewFilename("{}_msh_".format(save),
saveFormat)
else:
fileOut = "{}_msh.{}".format(save, saveFormat)
fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show: plt.show()
out = fig if fileOut is None else (fig, fileOut)
return bcast(out)
def outParaview(self, u:Np1D, warping : List[callable] = None,
is_state : bool = False, name : str = "u",
filename : str = "out", time : float = 0.,
what : strLst = 'all', forceNewFile : bool = True,
folder : bool = False, filePW = None) -> str:
"""
Output complex-valued function with given dofs to ParaView file.
Args:
u: numpy complex array with function dofs.
warping(optional): Domain warping functions.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
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).
Returns:
Output filename.
"""
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
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
filePW = None
if masterCore():
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 warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[0], self.V.mesh()))
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)
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[1], self.V.mesh()))
return bcast(filePW)
def outParaviewTimeDomain(self, u:Np1D, omega:float,
warping : List[callable] = None,
is_state : bool = False,
timeFinal : float = None,
periodResolution : int = 20, name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folder : bool = False) -> str:
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
Args:
u: numpy complex array with function dofs.
omega: frequency.
warping(optional): Domain warping functions.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
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.
Returns:
Output filename.
"""
if not is_state and not self.isCEye:
raise RROMPyException(("Cannot output to Paraview non-state "
"object."))
filePW = None
if masterCore():
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
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[0], self.V.mesh()))
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
if warping is not None:
fen.ALE.move(self.V.mesh(),
interp_project(warping[1], self.V.mesh()))
return bcast(filePW)
class FenicsEngineBaseTensorized(FenicsEngineBase):
"""The number of tensorized dimensions should be assigned to nports."""
def plot(self, u:Np1D, warping : List[callable] = None,
is_state : bool = False, name : str = "u", save : str = None,
what : strLst = 'all', forceNewFile : bool = True,
saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
colorMap : str = "jet", fenplotArgs : dict = {},
**figspecs) -> Tuple[FigHandle, str]:
"""
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.
is_state(optional): whether given u is value before multiplication
by c. Defaults to False.
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'.
forceNewFile(optional): Whether to create new output file.
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.
colorMap(optional): Pyplot colormap. Defaults to 'jet'.
fenplotArgs(optional): Optional arguments for fenplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename and figure handle.
"""
nP = checknports(self)
if not is_state and not self.isCEye:
return super().plot(u.reshape(-1, nP), warping, False, name, save,
what, forceNewFile, saveFormat, saveDPI, show,
colorMap, fenplotArgs, **figspecs)
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
- filePW = None
+ out = None
if masterCore():
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = plt.figaspect(1. / len(what))
figspecs['figsize'][1] *= nP
fig = plt.figure(**figspecs)
plt.set_cmap(colorMap)
for i in range(nP):
plottingBaseFen(u[i :: nP], fig, self.V, what, nP,
i * len(what), warping,
"{}_port{}".format(name, i + 1),
self.V.mesh().geometric_dimension() > 1,
fenplotArgs)
plt.tight_layout()
if save is not None:
save = save.strip()
if forceNewFile:
fileOut = getNewFilename("{}_fig_".format(save),
saveFormat)
else:
fileOut = "{}_fig.{}".format(save, saveFormat)
fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show: plt.show()
out = fig if fileOut is None else (fig, fileOut)
return bcast(out)
def outParaview(self, u:Np1D, *args, **kwargs) -> List[str]:
nP = checknports(self)
idx = indicesScatter(nP)[0]
filesOut = []
if len(idx) > 0:
for j in idx:
filesOut += super().outParaview(u[j :: nP], *args, **kwargs)
filesOut = listGather(filesOut)
if filesOut[0] is None: return None
return filesOut
def outParaviewTimeDomain(self, u:Np1D, *args, **kwargs) -> List[str]:
nP = checknports(self)
idx = indicesScatter(nP)[0]
filesOut = []
if len(idx) > 0:
for j in idx:
filesOut += super().outParaviewTimeDomain(u[j :: nP], *args,
**kwargs)
filesOut = listGather(filesOut)
if filesOut[0] is None: return None
return filesOut
diff --git a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
index 6c19311..2a589b8 100644
--- a/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
+++ b/rrompy/hfengines/vector_linear_problem/linear_elasticity_problem_engine.py
@@ -1,293 +1,289 @@
# 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
from rrompy.hfengines.base.linear_affine_engine import LinearAffineEngine
from rrompy.hfengines.base.vector_fenics_engine_base import \
VectorFenicsEngineBase
from rrompy.utilities.base.types import paramVal
from rrompy.solver.fenics import (fenZERO, fenZEROS, fenONE,
elasticNormMatrix, elasticDualNormMatrix)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.parameter import checkParameter
from rrompy.solver.fenics import fenics2Sparse, fenics2Vector
__all__ = ['LinearElasticityProblemEngine']
class LinearElasticityProblemEngine(LinearAffineEngine,
VectorFenicsEngineBase):
"""
Solver for generic linear elasticity problems.
- div(lambda_ * div(u) * I + 2 * mu_ * epsilon(u)) = f in \Omega
u = u0 on \Gamma_D
\partial_nu = g1 on \Gamma_N
\partial_nu + h u = g2 on \Gamma_R
Attributes:
verbosity: Verbosity level.
BCManager: Boundary condition manager.
V: Real vector FE space.
u: Generic vector trial functions for variational form evaluation.
v: Generic vector test functions for variational form evaluation.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product over
V.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product between Riesz representers V-V.
degree_threshold: Threshold for ufl expression interpolation degree.
lambda_: Value of lambda_.
mu_: Value of mu_.
forcingTerm: Value of f.
DirichletDatum: Value of u0.
NeumannDatum: Value of g1.
RobinDatumG: Value of g2.
RobinDatumH: Value of h.
DirichletBoundary: Function handle to \Gamma_D.
NeumannBoundary: Function handle to \Gamma_N.
RobinBoundary: Function handle to \Gamma_R.
ds: Boundary measure 2-tuple (resp. for Neumann and Robin boundaries).
dsToBeSet: Whether ds needs to be set.
"""
_energyDualNormCompress = None
def __init__(self, mu0 : paramVal = [], degree_threshold : int = np.inf,
verbosity : int = 10, timestamp : bool = True):
super().__init__(degree_threshold = degree_threshold,
verbosity = verbosity, timestamp = timestamp)
self._affinePoly = True
self.lambda_ = fenONE
self.mu_ = fenONE
self.mu0 = checkParameter(mu0)
self.npar = self.mu0.shape[1]
- self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
- self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
- self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
- self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumH = fenZERO
@property
def V(self):
"""Value of V."""
return self._V
@V.setter
def V(self, V):
VectorFenicsEngineBase.V.fset(self, V)
self.forcingTerm = fenZEROS(self.V.mesh().topology().dim())
self.DirichletDatum = fenZEROS(self.V.mesh().topology().dim())
self.NeumannDatum = fenZEROS(self.V.mesh().topology().dim())
self.RobinDatumG = fenZEROS(self.V.mesh().topology().dim())
self.dsToBeSet = True
@property
def lambda_(self):
"""Value of lambda_."""
return self._lambda_
@lambda_.setter
def lambda_(self, lambda_):
self.resetAs()
if not isinstance(lambda_, (list, tuple,)):
lambda_ = [lambda_, fenZERO]
self._lambda_ = lambda_
@property
def mu_(self):
"""Value of mu_."""
return self._mu_
@mu_.setter
def mu_(self, mu_):
self.resetAs()
if not isinstance(mu_, (list, tuple,)):
mu_ = [mu_, fenZERO]
self._mu_ = mu_
@property
def forcingTerm(self):
"""Value of f."""
return self._forcingTerm
@forcingTerm.setter
def forcingTerm(self, forcingTerm):
self.resetbs()
if not isinstance(forcingTerm, (list, tuple,)):
forcingTerm = [forcingTerm,
fenZEROS(self.V.mesh().topology().dim())]
self._forcingTerm = forcingTerm
@property
def DirichletDatum(self):
"""Value of u0."""
return self._DirichletDatum
@DirichletDatum.setter
def DirichletDatum(self, DirichletDatum):
self.resetbs()
if not isinstance(DirichletDatum, (list, tuple,)):
DirichletDatum = [DirichletDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._DirichletDatum = DirichletDatum
@property
def NeumannDatum(self):
"""Value of g1."""
return self._NeumannDatum
@NeumannDatum.setter
def NeumannDatum(self, NeumannDatum):
self.resetbs()
if not isinstance(NeumannDatum, (list, tuple,)):
NeumannDatum = [NeumannDatum,
fenZEROS(self.V.mesh().topology().dim())]
self._NeumannDatum = NeumannDatum
@property
def RobinDatumG(self):
"""Value of g2."""
return self._RobinDatumG
@RobinDatumG.setter
def RobinDatumG(self, RobinDatumG):
self.resetbs()
if not isinstance(RobinDatumG, (list, tuple,)):
RobinDatumG = [RobinDatumG,
fenZEROS(self.V.mesh().topology().dim())]
self._RobinDatumG = RobinDatumG
@property
def RobinDatumH(self):
"""Value of h."""
return self._RobinDatumH
@RobinDatumH.setter
def RobinDatumH(self, RobinDatumH):
self.resetAs()
if not isinstance(RobinDatumH, (list, tuple,)):
RobinDatumH = [RobinDatumH, fenZERO]
self._RobinDatumH = RobinDatumH
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self.BCManager.DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self.resetAs()
self.resetbs()
self.BCManager.DirichletBoundary = DirichletBoundary
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self.BCManager.NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.NeumannBoundary = NeumannBoundary
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self.BCManager.RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self.resetAs()
self.resetbs()
self.dsToBeSet = True
self.BCManager.RobinBoundary = RobinBoundary
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
vbMng(self, "INIT", "Assembling energy matrix.", 20)
self.energyNormMatrix = elasticNormMatrix(self.V, self.lambda_[0],
self.mu_[0])
vbMng(self, "DEL", "Done assembling energy matrix.", 20)
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
vbMng(self, "INIT", "Assembling energy dual matrix.", 20)
self.energyNormDualMatrix = elasticDualNormMatrix(
self.V, self.lambda_[0], self.mu_[0],
compressRank = self._energyDualNormCompress)
vbMng(self, "DEL", "Done assembling energy dual matrix.", 20)
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling operator term A0.", 20)
DirichletBC0 = fen.DirichletBC(self.V,
fenZEROS(self.V.mesh().topology().dim()),
self.DirichletBoundary)
lambda_Re, lambda_Im = self.lambda_
mu_Re, mu_Im = self.mu_
hRe, hIm = self.RobinDatumH
termNames = ["lambda_", "mu_", "RobinDatumH"]
parsRe = self.iterReduceQuadratureDegree(zip(
[lambda_Re, mu_Re, hRe],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip(
[lambda_Im, mu_Re, hIm],
[x + "Imag" for x in termNames]))
epsilon = lambda u: 0.5 * (fen.grad(u) + fen.nabla_grad(u))
sigma = lambda u, l_, m_: (
l_ * fen.div(u) * fen.Identity(u.geometric_dimension())
+ 2. * m_ * epsilon(u))
a0Re = (fen.inner(sigma(self.u, lambda_Re, mu_Re),
epsilon(self.v)) * fen.dx
+ hRe * fen.inner(self.u, self.v) * self.ds(1))
a0Im = (fen.inner(sigma(self.u, lambda_Im, mu_Im),
epsilon(self.v)) * fen.dx
+ hIm * fen.inner(self.u, self.v) * self.ds(1))
self.As[0] = (fenics2Sparse(a0Re, parsRe, DirichletBC0, 1)
+ 1.j * fenics2Sparse(a0Im, parsIm, DirichletBC0, 0))
vbMng(self, "DEL", "Done assembling operator term.", 20)
def buildb(self):
"""Build terms of operator of linear system."""
if self.thbs[0] is None:
self.thbs = self.getMonomialWeights(self.nbs)
if self.bs[0] is None:
self.autoSetDS()
vbMng(self, "INIT", "Assembling forcing term b0.", 20)
u0Re, u0Im = self.DirichletDatum
fRe, fIm = self.forcingTerm
g1Re, g1Im = self.NeumannDatum
g2Re, g2Im = self.RobinDatumG
termNames = ["forcingTerm", "NeumannDatum", "RobinDatumG"]
parsRe = self.iterReduceQuadratureDegree(zip([fRe, g1Re, g2Re],
[x + "Real" for x in termNames]))
parsIm = self.iterReduceQuadratureDegree(zip([fIm, g1Im, g2Im],
[x + "Imag" for x in termNames]))
L0Re = (fen.inner(fRe, self.v) * fen.dx
+ fen.inner(g1Re, self.v) * self.ds(0)
+ fen.inner(g2Re, self.v) * self.ds(1))
L0Im = (fen.inner(fIm, self.v) * fen.dx
+ fen.inner(g1Im, self.v) * self.ds(0)
+ fen.inner(g2Im, self.v) * self.ds(1))
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
DBCI = fen.DirichletBC(self.V, u0Im, self.DirichletBoundary)
self.bs[0] = (fenics2Vector(L0Re, parsRe, DBCR, 1)
+ 1.j * fenics2Vector(L0Im, parsIm, DBCI, 1))
vbMng(self, "DEL", "Done assembling forcing term.", 20)
diff --git a/rrompy/solver/fenics/__init__.py b/rrompy/solver/fenics/__init__.py
index 6936f3c..5043f0b 100644
--- a/rrompy/solver/fenics/__init__.py
+++ b/rrompy/solver/fenics/__init__.py
@@ -1,48 +1,51 @@
# 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/>.
#
from .fenics_constants import (fenZERO, fenZEROS, fenONE, fenONES,
bdrTrue, bdrFalse)
from .fenics_la import fenics2Sparse, fenics2Vector
+from .fenics_serial import serializeMesh, serializeFunctionSpace
from .fenics_norms import (L2NormMatrix, L2InverseNormMatrix, H1NormMatrix,
Hminus1NormMatrix, elasticNormMatrix,
elasticDualNormMatrix)
from .fenics_plotting import fenplot, affine_warping
from .fenics_projecting import interp_project
__all__ = [
'fenZERO',
'fenZEROS',
'fenONE',
'fenONES',
'bdrTrue',
'bdrFalse',
'fenics2Sparse',
'fenics2Vector',
+ 'serializeMesh',
+ 'serializeFunctionSpace',
'L2NormMatrix',
'L2InverseNormMatrix',
'H1NormMatrix',
'Hminus1NormMatrix',
'elasticNormMatrix',
'elasticDualNormMatrix',
'fenplot',
'affine_warping',
'interp_project'
]
diff --git a/rrompy/solver/fenics/fenics_serial.py b/rrompy/solver/fenics/fenics_serial.py
new file mode 100644
index 0000000..89b1160
--- /dev/null
+++ b/rrompy/solver/fenics/fenics_serial.py
@@ -0,0 +1,49 @@
+# 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/>.
+#
+
+from os import remove
+import fenics as fen
+from rrompy.utilities.base.data_structures import getNewFilename
+from rrompy.utilities.parallel import SELF
+
+__all__ = ['serializeMesh', 'serializeFunctionSpace']
+
+def serializeMesh(mesh):
+ mesh_comm = mesh.mpi_comm()
+ if mesh_comm.size == 1: return mesh
+ if mesh_comm.Get_rank() == 0:
+ filename = getNewFilename("mesh", "h5", False)
+ else:
+ filename = None
+ filename = mesh_comm.bcast(filename, root = 0)
+ fout = fen.HDF5File(mesh_comm, filename, "w")
+ fout.write(mesh, "mesh")
+ fout.close()
+ del fout
+ fin = fen.HDF5File(SELF, filename, "r")
+ meshS = fen.Mesh(SELF)
+ fin.read(meshS, "mesh", False)
+ fin.close()
+ del fin
+ mesh_comm.Barrier()
+ if mesh_comm.Get_rank() == 0: remove(filename)
+ return meshS
+
+def serializeFunctionSpace(space):
+ return fen.FunctionSpace(serializeMesh(space.mesh()), space.ufl_element(),
+ constrained_domain = space.dofmap().constrained_domain)
diff --git a/rrompy/utilities/base/verbosity_depth.py b/rrompy/utilities/base/verbosity_depth.py
index 9dae790..a3c674e 100644
--- a/rrompy/utilities/base/verbosity_depth.py
+++ b/rrompy/utilities/base/verbosity_depth.py
@@ -1,96 +1,87 @@
# 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/>.
#
from copy import deepcopy as copy
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ["verbosityDepth", "verbosityManager"]
from datetime import datetime
def getTimestamp() -> str:
time = datetime.now().strftime("%H:%M:%S.%f")
return "\x1b[42m{}\x1b[0m".format(time)
-def updateVerbosityCheckpoint(vdtype:str) -> str:
+def updateVerbosityCheckpoint(vctype:int) -> str:
global RROMPy_verbosity_checkpoint, RROMPy_verbosity_buffer
- assert isinstance(vdtype, str)
- vdtype = vdtype.upper()
- if vdtype not in ["UP", "DOWN"]:
- raise RROMPyException("Verbosity checkpoint type not recognized.")
if "RROMPy_verbosity_checkpoint" not in globals():
RROMPy_verbosity_checkpoint = 0
- if vdtype == "UP":
- RROMPy_verbosity_checkpoint += 1
- if "RROMPy_verbosity_buffer" not in globals():
- RROMPy_verbosity_buffer = ""
- else:
- if "RROMPy_verbosity_checkpoint" in globals():
- RROMPy_verbosity_checkpoint -= 1
- if RROMPy_verbosity_checkpoint <= 0:
- assert "RROMPy_verbosity_buffer" in globals()
- buffer = copy(RROMPy_verbosity_buffer)
- del RROMPy_verbosity_buffer
- return buffer
+ RROMPy_verbosity_checkpoint += vctype
+ if "RROMPy_verbosity_buffer" not in globals():
+ RROMPy_verbosity_buffer = ""
+ if RROMPy_verbosity_checkpoint <= 0:
+ buffer = copy(RROMPy_verbosity_buffer)
+ del RROMPy_verbosity_buffer
+ return buffer
return None
def verbosityDepth(vdtype:str, message:str, end : str = "\n",
timestamp : bool = True):
global RROMPy_verbosity_depth, RROMPy_verbosity_checkpoint, \
RROMPy_verbosity_buffer
assert isinstance(vdtype, str)
vdtype = vdtype.upper()
if vdtype not in ["INIT", "MAIN", "DEL"]:
raise RROMPyException("Verbosity depth type not recognized.")
if "RROMPy_verbosity_checkpoint" not in globals():
RROMPy_verbosity_checkpoint = 0
out = "{} ".format(getTimestamp()) if timestamp else ""
if vdtype == "INIT":
if "RROMPy_verbosity_depth" not in globals():
RROMPy_verbosity_depth = 0
RROMPy_verbosity_depth += 1
out += "│" * (RROMPy_verbosity_depth - 1)
out += "┌"
else:
assert "RROMPy_verbosity_depth" in globals()
if vdtype == "MAIN":
out += "│" * (RROMPy_verbosity_depth - 1)
out += "├"
elif vdtype == "DEL":
RROMPy_verbosity_depth -= 1
out += "│" * RROMPy_verbosity_depth
out += "└"
if RROMPy_verbosity_depth <= 0: del RROMPy_verbosity_depth
from rrompy.utilities.parallel import poolRank, poolSize, masterCore
if message != "" and masterCore():
if RROMPy_verbosity_checkpoint and poolSize() > 1:
poolrk = "{{\x1b[34m{}\x1b[0m}}".format(poolRank())
else:
poolrk = ""
msg = "{}{}{}{}".format(out, poolrk, message, end)
if RROMPy_verbosity_checkpoint:
assert "RROMPy_verbosity_buffer" in globals()
RROMPy_verbosity_buffer += msg
else:
print(msg, end = "", flush = True)
return
def verbosityManager(object, vdtype:str, message:str, vlvl : int = 0,
end : str = "\n"):
if object.verbosity >= vlvl:
return verbosityDepth(vdtype, message, end, object.timestamp)
diff --git a/rrompy/utilities/parallel/__init__.py b/rrompy/utilities/parallel/__init__.py
index 3724fe9..d0d2cf2 100644
--- a/rrompy/utilities/parallel/__init__.py
+++ b/rrompy/utilities/parallel/__init__.py
@@ -1,43 +1,46 @@
# 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/>.
#
-from .base import (_MPI_IS_LOADED, COMM, poolRank, poolSize, masterCore,
- forceSerialExecution, allowParallelExecution, forcedSerial)
+from .base import (_MPI_IS_LOADED, COMM, SELF, poolRank, poolSize, masterCore,
+ updateSerialIndex, forceSerialExecution,
+ allowParallelExecution, forcedSerial)
from .partition import partitionIndices
from .broadcast import bcast
from .scatter import indicesScatter, listScatter
from .gather import listGather, matrixGather
__all__ = [
'_MPI_IS_LOADED',
'COMM',
+ 'SELF',
'poolRank',
'poolSize',
'masterCore',
+ 'updateSerialIndex',
'forceSerialExecution',
'allowParallelExecution',
'forcedSerial',
'partitionIndices',
'bcast',
'indicesScatter',
'listScatter',
'listGather',
'matrixGather'
]
diff --git a/rrompy/utilities/parallel/base.py b/rrompy/utilities/parallel/base.py
index df1f319..d8b2e4c 100644
--- a/rrompy/utilities/parallel/base.py
+++ b/rrompy/utilities/parallel/base.py
@@ -1,63 +1,67 @@
# 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/>.
#
from rrompy.utilities.base.verbosity_depth import updateVerbosityCheckpoint
try:
from mpi4py import MPI
_MPI_IS_LOADED = True
- COMM = MPI.COMM_WORLD
+ COMM, SELF = MPI.COMM_WORLD, MPI.COMM_SELF
except ImportError:
_MPI_IS_LOADED = False
- COMM = None
+ COMM, SELF = None, None
-__all__ = ['_MPI_IS_LOADED', 'COMM', 'poolRank', 'poolSize', 'masterCore',
- 'forceSerialExecution', 'allowParallelExecution', 'forcedSerial']
+__all__ = ['_MPI_IS_LOADED', 'COMM', 'SELF', 'poolRank', 'poolSize',
+ 'masterCore', 'updateSerialIndex', 'forceSerialExecution',
+ 'allowParallelExecution', 'forcedSerial']
def poolRank() -> int:
if _MPI_IS_LOADED: return COMM.Get_rank()
return 0
def poolSize() -> int:
if _MPI_IS_LOADED: return COMM.Get_size()
return 1
def masterCore() -> bool:
return forcedSerial() or poolRank() == 0
-def forceSerialExecution(verbosity : bool = False):
+def updateSerialIndex(n : int = 0, delta : bool = True):
global RROMPy_force_serial
if "RROMPy_force_serial" not in globals(): RROMPy_force_serial = 0
- RROMPy_force_serial += 1
- if verbosity: updateVerbosityCheckpoint("UP")
+ if delta: RROMPy_force_serial += n
+ else: RROMPy_force_serial = n
+ if RROMPy_force_serial <= 0: del RROMPy_force_serial
+
+def forceSerialExecution(verbosity : bool = False):
+ updateSerialIndex(1)
+ if verbosity: updateVerbosityCheckpoint(1)
def allowParallelExecution(verbosity : bool = False):
- global RROMPy_force_serial
- if "RROMPy_force_serial" in globals(): RROMPy_force_serial -= 1
- if RROMPy_force_serial <= 0: del RROMPy_force_serial
+ updateSerialIndex(-1)
if verbosity:
- buffer = updateVerbosityCheckpoint("DOWN")
+ buffer = updateVerbosityCheckpoint(-1)
if buffer is not None:
msg = COMM.gather(buffer, root = 0)
if masterCore():
out = ""
for ms in msg: out += ms
print(out, end = "", flush = True)
def forcedSerial():
global RROMPy_force_serial
return "RROMPy_force_serial" in globals() and RROMPy_force_serial > 0