Page MenuHomec4science
Files F121336970

fenics_plotting.py
No OneTemporary
Actions

File Metadata

Created
Thu, Jul 10, 03:28

fenics_plotting.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
import dolfin.cpp as cpp
import ufl
import fenics as fen
from rrompy.utilities.base.types import Np1D, Np2D
from .fenics_projecting import interp_project
from rrompy.utilities.exception_manager import RROMPyException
_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:
fen.ALE.move(mesh, interp_project(warping[0], mesh))
out = fen.plot(object, *args, mesh = mesh, **kwargs)
if warping is not None:
fen.ALE.move(mesh, interp_project(warping[1], mesh))
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
try:
Ainv = np.linalg.inv(A)
except np.linalg.LinAlgError as e:
raise RROMPyException(e)
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]])

Event Timeline

c4science · Help