diff --git a/.gitignore b/.gitignore index 8062809..87b636f 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,6 @@ example/bar2vtk/result/* # testing .coverage + +# Documentation +docs/_build/* diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/api/requirements.txt b/docs/api/requirements.txt new file mode 100644 index 0000000..06f6e99 --- /dev/null +++ b/docs/api/requirements.txt @@ -0,0 +1,2 @@ +sphinx-rtd-theme +autodocsumm diff --git a/docs/api/vtkpytools.barfiletools.rst b/docs/api/vtkpytools.barfiletools.rst new file mode 100644 index 0000000..4e431f6 --- /dev/null +++ b/docs/api/vtkpytools.barfiletools.rst @@ -0,0 +1,20 @@ +Bar File tools (:mod:`barfiletools`) +================================================ + +Creating VTM files (:mod:`barfiletools.bar2vtk`) +------------------------------------------------------------- + +.. automodule:: vtkpytools.barfiletools.bar2vtk + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: + +Post-Processing VTM files (:mod:`barfiletools.data`) +----------------------------------------------------------- + +.. automodule:: vtkpytools.barfiletools.data + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: diff --git a/docs/api/vtkpytools.bl.rst b/docs/api/vtkpytools.bl.rst new file mode 100644 index 0000000..8270536 --- /dev/null +++ b/docs/api/vtkpytools.bl.rst @@ -0,0 +1,9 @@ +Boundary Layer tools (:mod:`bl`) +============================================================= + +.. automodule:: vtkpytools.bl + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: + diff --git a/docs/api/vtkpytools.common.rst b/docs/api/vtkpytools.common.rst new file mode 100644 index 0000000..cdf5a64 --- /dev/null +++ b/docs/api/vtkpytools.common.rst @@ -0,0 +1,7 @@ +Common Utilities module (:mod:`common`) +========================================= + +.. automodule:: vtkpytools.common + :members: + :autosummary: + :autosummary-nosignatures: diff --git a/docs/api/vtkpytools.gridtools2d.rst b/docs/api/vtkpytools.gridtools2d.rst new file mode 100644 index 0000000..fcba06b --- /dev/null +++ b/docs/api/vtkpytools.gridtools2d.rst @@ -0,0 +1,9 @@ +Tools for 2D Grids (:mod:`gridtools2d`) +=========================================== + +.. automodule:: vtkpytools.gridtools2d.core + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: + :show-inheritance: diff --git a/docs/api/vtkpytools.gridtools3d.rst b/docs/api/vtkpytools.gridtools3d.rst new file mode 100644 index 0000000..ab538d3 --- /dev/null +++ b/docs/api/vtkpytools.gridtools3d.rst @@ -0,0 +1,9 @@ +Tools for 3D Grids (:mod:`gridtools3d`) +========================================= + +.. automodule:: vtkpytools.gridtools3d.core + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: + :show-inheritance: diff --git a/docs/api/vtkpytools.numtools.rst b/docs/api/vtkpytools.numtools.rst new file mode 100644 index 0000000..76c7d9f --- /dev/null +++ b/docs/api/vtkpytools.numtools.rst @@ -0,0 +1,8 @@ +Numerical Tools (:mod:`numtools`) +================================== + +.. automodule:: vtkpytools.numtools + :members: + :undoc-members: + :autosummary: + :autosummary-nosignatures: diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..5b0cd33 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,66 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +import sphinx_rtd_theme +sys.path.insert(0, os.path.abspath('../')) + + +# -- Project information ----------------------------------------------------- + +project = 'vtkpytools' +copyright = '2021, James Wright' +author = 'James Wright' + +# The full version, including alpha/beta/rc tags +release = '1.1.0' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'autodocsumm', + 'sphinx.ext.autodoc', + 'sphinx.ext.autosummary', + 'sphinx.ext.napoleon', + 'sphinx.ext.mathjax', + 'sphinx_rtd_theme', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +autodoc_member_order = 'bysource' + +default_role = 'autolink' diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..1c32c1c --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,29 @@ +Welcome to vtkpytools's documentation! +====================================== + +Subpackages +----------- + +.. toctree:: + :maxdepth: 4 + + api/vtkpytools.barfiletools + api/vtkpytools.gridtools2d + api/vtkpytools.gridtools3d + +Modules +-------------------- + +.. toctree:: + :maxdepth: 4 + + api/vtkpytools.bl + api/vtkpytools.common + api/vtkpytools.numtools + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..2119f51 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..06f6e99 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,2 @@ +sphinx-rtd-theme +autodocsumm diff --git a/tests/test_numtools.py b/tests/test_numtools.py index 518c8c6..0044a8c 100644 --- a/tests/test_numtools.py +++ b/tests/test_numtools.py @@ -31,19 +31,71 @@ def test_calcStrainRate(twoFullTensors): assert np.allclose(strainrate, solution) assert strainrate.shape[0] == twoFullTensors.shape[0] +class Test_symmetric2FullTensor(): -def test_symmetric2FullTensor(symAndFullTensor): - (full, sym) = symAndFullTensor + @staticmethod + def test_func2d(symAndFullTensor): + full, sym = symAndFullTensor - full_test = vpt.symmetric2FullTensor(sym) - assert np.allclose(full_test, full) + full_test = vpt.symmetric2FullTensor(sym) + assert np.allclose(full_test, full) + @staticmethod + def test_func_Nd(symAndFullTensor): + """Test (n,m,6) -> (n,m,3,3)""" + full, sym = symAndFullTensor -def test_full2SymmetricTensor(symAndFullTensor): - (full, sym) = symAndFullTensor + sym = sym.reshape((1,-1,6)) + full = full.reshape((1,-1,3,3)) + full_test = vpt.symmetric2FullTensor(sym) + assert full_test.shape == full.shape + assert np.allclose(full_test, full) - sym_test = vpt.full2SymmetricTensor(full) - assert np.allclose(sym_test, sym) + sym = sym.reshape((-1,1,6)) + full = full.reshape((-1,1,3,3)) + full_test = vpt.symmetric2FullTensor(sym) + assert full_test.shape == full.shape + assert np.allclose(full_test, full) + + @staticmethod + def test_inputCheck(): + array = np.empty((6,4)) + with pytest.raises(Exception): + vpt.symmetric2FullTensor(array) + + +class Test_full2SymmetricTensor(): + + @staticmethod + def test_func3d(symAndFullTensor): + full, sym = symAndFullTensor + + sym_test = vpt.full2SymmetricTensor(full) + assert sym_test.shape == sym.shape + assert np.allclose(sym_test, sym) + + @staticmethod + def test_func_Nd(symAndFullTensor): + """Test (n,m,3,3) -> (n,m,6)""" + full, sym = symAndFullTensor + + sym = sym.reshape((1,-1,6)) + full = full.reshape((1,-1,3,3)) + sym_test = vpt.full2SymmetricTensor(full) + assert sym_test.shape == sym.shape + assert np.allclose(sym_test, sym) + + sym = sym.reshape((-1,1,6)) + full = full.reshape((-1,1,3,3)) + sym_test = vpt.full2SymmetricTensor(full) + assert sym_test.shape == sym.shape + assert np.allclose(sym_test, sym) + + @staticmethod + def test_inputCheck(): + array = np.empty((3,3,4)) + with pytest.raises(Exception): + vpt.full2SymmetricTensor(array) @pytest.mark.parametrize('offset,expected', [(1, [0, 1] ), @@ -53,7 +105,7 @@ def test_pwlinRoots(offset, expected): y = np.array([-1, 1, -1]) roots = vpt.pwlinRoots(x, y + offset) - assert(np.allclose(roots, expected)) + assert np.allclose(roots, expected) @pytest.mark.parametrize('kwargs,expected', [({'dx': 2.4}, [0, 1, 2, 4, 6.4, 8]), @@ -63,8 +115,8 @@ def test_seriesDiffLimiter(kwargs, expected): series = np.array([0, 1, 2, 4, 8]) result = vpt.seriesDiffLimiter(series, **kwargs) - assert(result.size == len(expected)) - assert(np.allclose(result, expected)) + assert result.size == len(expected) + assert np.allclose(result, expected) @pytest.mark.parametrize('kwargs', [{}, {'dx': 2.0, 'magnitude': 4.0}]) diff --git a/vtkpytools/__init__.py b/vtkpytools/__init__.py index 3fefedb..ce70e81 100644 --- a/vtkpytools/__init__.py +++ b/vtkpytools/__init__.py @@ -1,3 +1,4 @@ +from ._version import __version__ from .barfiletools import * from .gridtools3d import * from .gridtools2d import * @@ -8,3 +9,4 @@ pwlinRoots, seriesDiffLimiter) from .bl import (integratedVortBLThickness, sampleAlongVectors, delta_vortInt, delta_velInt, delta_percent) +from .misc import * diff --git a/vtkpytools/barfiletools/bar2vtk.py b/vtkpytools/barfiletools/bar2vtk.py index 47186e0..becfce5 100644 --- a/vtkpytools/barfiletools/bar2vtk.py +++ b/vtkpytools/barfiletools/bar2vtk.py @@ -1,13 +1,16 @@ +"""Functions related to the creation of the VTK files from bar* files """ + import vtk, os, argparse, datetime, platform, warnings import pytomlpp import pyvista as pv import numpy as np from pathlib import Path, PurePath -from .data import binaryVelbar, binaryStsbar, calcReynoldsStresses, compute_vorticity -from ..common import globFile +from .data import binaryVelbar, binaryStsbar, calcReynoldsStresses +from ..common import globFile, readBinaryArray from .._version import __version__ + def bar2vtk_main(args=None): """Function that runs the "binary" bar2vtk""" @@ -35,10 +38,13 @@ def bar2vtk_main(args=None): val[i] = Path(item) else: bar2vtkargs[key] = Path(val) + else: + raise RuntimeError tomlMetadata = bar2vtk_function(**bar2vtkargs, returnTomlMetadata=True) tomlReceipt(bar2vtkargs, tomlMetadata) + def bar2vtk_parse(args=None): GeneralDescription="""Tool for putting velbar and stsbar data onto a vtm grid.""" @@ -105,6 +111,9 @@ class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescri cliparser.add_argument('-a', '--ascii', help='Read *bar files as ASCII', action='store_true') cliparser.add_argument('--velbar', help='Path to velbar file(s)', type=Path, nargs='+', default=[]) cliparser.add_argument('--stsbar', help='Path to stsbar file(s)', type=Path, nargs='+', default=[]) + cliparser.add_argument('--consrvstress', + help='Calculate Reynolds stress assuming conservative stsbar data', + type=bool, nargs=1, default=False) # Toml Parser Setup tomlparser = subparser.add_parser('toml', description=TomlDescription, formatter_class=CustomFormatter, @@ -117,10 +126,11 @@ class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescri return args + def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ ts0: int=-1, new_file_prefix: str='', outpath: Path=None, \ velonly=False, debug=False, asciidata=False, \ - velbar=[], stsbar=[], returnTomlMetadata=False): + velbar=[], stsbar=[], returnTomlMetadata=False, consrvstress=False): """Convert velbar and stsbar files into 2D vtk files See bar2vtk_commandline help documentation for more information. @@ -129,7 +139,7 @@ def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ ---------- blankvtmfile : Path Path to blank VTM file. Data is loaded onto this file and then saved to - different file. VTM file should contain 'grid' and 'wall' blocks. + different file. VTM file should contain ``grid`` and ``wall`` blocks. barfiledir : Path Directory to where the velbar and stsbar files are located timestep : str @@ -143,22 +153,29 @@ def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ (Default: -1). new_file_prefix : str If given, the newly created file will take the form of - '{new_file_prefix}_{timestep}.vtm'. If not given, the file prefix will - be the same as blankvtmfile. (Default: ''). + ``{new_file_prefix}_{timestep}.vtm``. If not given, the file prefix will + be the same as blankvtmfile. (Default: ``''``). outpath : Path Path where the new file should be written. If not given, the new file will be placed in the current working directory. velonly : bool - Do not include the stsbar file. (Default: False) + Do not include the stsbar file. (Default: ``False``) debug : bool Adds raw stsbar array as data field. This is purely for debugging - purposes. (Default: False) + purposes. (Default: ``False``) asciidata : bool - Whether file paths are in ASCII format. (Default: False) + Whether file paths are in ASCII format. (Default: ``False``) velbar : List of Path - Path(s) to velbar files. If doing time windows, must have two Paths. (Default: []) + Path(s) to velbar files. If doing time windows, must have two Paths. (Default: ``[]``) stsbar : List of Path Path(s) to stsbar files. If doing time windows, must have two Paths. (Default: []) + returnTomlMetadata : bool + Whether to return the metadata required for writing a toml reciept file. (Default: False) + consrvstress : bool + Whether the data in the stsbar file is conservative or not. This + affects both how the stsbar file is read and how the Reynolds stresses + are calculated from them. See calcReynoldsStresses() for more + information. (Default: False) """ ## ---- Process/check script arguments @@ -189,7 +206,12 @@ def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ vtmPath = (outpath if outpath else os.getcwd()) / vtmName velbarReader = np.loadtxt if asciidata else binaryVelbar - stsbarReader = np.loadtxt if asciidata else binaryStsbar + if asciidata: + stsbarReader = np.loadtxt + elif consrvstress: + stsbarReader = lambda path: readBinaryArray(path, ncols=9) + else: + stsbarReader = binaryStsbar ## ---- Loading data arrays if '-' in timestep and ts0 == -1: @@ -211,7 +233,7 @@ def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ grid['Velocity'] = velbarArray[:,1:4] if not velonly: - ReyStrTensor = calcReynoldsStresses(stsbarArray, velbarArray) + ReyStrTensor = calcReynoldsStresses(stsbarArray, velbarArray, consrvstress) grid['ReynoldsStress'] = ReyStrTensor if debug and not velonly: @@ -236,6 +258,7 @@ def bar2vtk_function(blankvtmfile: Path, barfiledir: Path, timestep: str, \ return tomlMetadata + def getBarData(_bar: list, timestep_str: str, barfiledir: Path, _barReader, ts0: int, globname: str): """Get array of data from bar2vtk arguments""" @@ -266,6 +289,7 @@ def getBarData(_bar: list, timestep_str: str, barfiledir: Path, _barReader, return _barArray, _barPaths + def blankToml(tomlfilepath: Path, returndict=False): """Write blank toml file to tomlfilepath""" tomldict = { 'arguments': { @@ -280,12 +304,14 @@ def blankToml(tomlfilepath: Path, returndict=False): 'asciidata': False, 'velbar': [], 'stsbar': [], + 'consrvstress': False, }} with tomlfilepath.open('w') as file: pytomlpp.dump(tomldict, file) if returndict: return tomldict + def tomlReceipt(args: dict, tomlMetadata: dict): """Creates a receipt of the created file @@ -298,7 +324,7 @@ def tomlReceipt(args: dict, tomlMetadata: dict): """ convertArray = [(PurePath, lambda x: x.as_posix()), - (type(None), lambda x: '') + (type(None), lambda _: '') ] _convertArray2TomlTypes(args, convertArray) @@ -331,6 +357,7 @@ def tomlReceipt(args: dict, tomlMetadata: dict): with tomlPath.open(mode='w') as file: pytomlpp.dump(tomldict, file) + def _convertArray2TomlTypes(array, convertArray: list): """Recursively convert objects in array to toml compatible objects based on convertArray diff --git a/vtkpytools/barfiletools/data.py b/vtkpytools/barfiletools/data.py index 9a4fee5..ce4c11b 100644 --- a/vtkpytools/barfiletools/data.py +++ b/vtkpytools/barfiletools/data.py @@ -4,11 +4,14 @@ from ..common import readBinaryArray, Profile, vCutter from ..numtools import makeRotationTensor import warnings +from os import PathLike +from typing import Union, Optional +from typing_extensions import Literal -def binaryVelbar(velbar_path) -> np.ndarray: +def binaryVelbar(velbar_path: Union[str, PathLike]) -> np.ndarray: """Get velbar array from binary file. - Wrapping around vpt.readBinaryArray. Assumes that the number of columns in + Wrapping around `vpt.readBinaryArray`. Assumes that the number of columns in the array is 5. Parameters @@ -18,21 +21,30 @@ def binaryVelbar(velbar_path) -> np.ndarray: """ return readBinaryArray(velbar_path, 5) -def binaryStsbar(stsbar_path) -> np.ndarray: +def binaryStsbar(stsbar_path: Union[str, PathLike]) -> np.ndarray: """Get stsbar array from binary file. - Wrapping around vpt.readBinaryArray. Assumes that the number of columns in + Wrapping around `vpt.readBinaryArray`. Assumes that the number of columns in the array is 6. Parameters ---------- - stsbar_path : Path + stsbar_path : PathLike Path to stsbar file. """ return readBinaryArray(stsbar_path, 6) -def calcReynoldsStresses(stsbar_array, velbar_array, conservative_stresses=False) -> np.ndarray: - """Calculate Reynolds Stresses from velbar and stsbar data. +def calcReynoldsStresses(stsbar_array: np.ndarray, velbar_array: np.ndarray, + conservative_stresses: bool=False) -> np.ndarray: + r"""Calculate Reynolds Stresses from velbar and stsbar data. + + Calculates via: + + .. math:: \langle u_i' u_j' \rangle = \langle u_i u_j \rangle - \langle u_i + \rangle \langle u_j \rangle + + where :math:`u_i` are the instantaneous velocity components and + :math:`u_i'` is the fluctuating velocity component. Parameters ---------- @@ -41,16 +53,14 @@ def calcReynoldsStresses(stsbar_array, velbar_array, conservative_stresses=False velbar_array : ndarray Array of velbar values conservative_stresses : bool - Whether the stsbar file used the - 'Conservative Stresses' option (default:False) + Whether the stsbar file used the 'Conservative Stresses' option + (Default: ``False``) Returns ------- numpy.ndarray """ if conservative_stresses: - warnings.warn("Calculation of Reynolds Stresses when using the " - "'Conservative Stress' option for stsbar has not been validated.") ReyStrTensor = np.empty((stsbar_array.shape[0], 6)) ReyStrTensor[:,0] = stsbar_array[:,3] - stsbar_array[:,0]**2 ReyStrTensor[:,1] = stsbar_array[:,4] - stsbar_array[:,1]**2 @@ -71,7 +81,7 @@ def calcReynoldsStresses(stsbar_array, velbar_array, conservative_stresses=False return ReyStrTensor -def calcWallShearGradient(wall) -> np.ndarray: +def calcWallShearGradient(wall: pv.DataSet) -> np.ndarray: """Calcuate the shear gradient at the wall Wall shear gradient is defined as the gradient of the velocity tangent to @@ -79,11 +89,11 @@ def calcWallShearGradient(wall) -> np.ndarray: multiply by mu. If the wall normal unit vector is taken to be n and the - gradient of velocity is e_ij, then the wall shear gradient is: + gradient of velocity is :math:`e_{ij}`, then the wall shear gradient is: - (delta_ik - n_k n_i) n_j e_kj + .. math:: (\delta_{ik} - n_k n_i) n_j e_{kj} - where n_j e_kj is the "traction" vector at the wall. See post at + where :math:`n_j e_{kj}` is the "traction" vector at the wall. See post at https://www.jameswright.xyz/post/wall_shear_gradient_from_velocity_gradient/ for derivation of method. @@ -111,17 +121,18 @@ def calcWallShearGradient(wall) -> np.ndarray: return wall_shear_gradient -def calcCf(wall, Uref, nu, rho, plane_normal='XY') -> np.ndarray: - """Calcuate the Coefficient of Friction of the wall +def calcCf(wall: pv.DataSet, Uref: float, nu: float, rho: float, + plane_normal: Literal['XY', 'XZ', 'YZ'] ='XY') -> np.ndarray: + r"""Calcuate the Coefficient of Friction of the wall Uses vpt.calcWallShearGradient to get du/dn, then uses input values to - calculate Cf using: + calculate :math:`C_f` using: - C_f = T_w / (0.5 * rho * Uref^2) + .. math:: C_f = \frac{\tau_w}{0.5 * \rho * U_\mathrm{ref}^2} Parameters ---------- - wall : pv.PolyData + wall : pv.DataSet Wall cells and points Uref : float Reference velocity @@ -131,7 +142,7 @@ def calcCf(wall, Uref, nu, rho, plane_normal='XY') -> np.ndarray: Density plane_normal : {'XY', 'XZ', 'YZ'}, optional Plane that the wall lies on. The shear stress vector will be projected - onto it. (default: 'XY') + onto it. (default: ``'XY'``) Returns ------- @@ -156,6 +167,9 @@ def calcCf(wall, Uref, nu, rho, plane_normal='XY') -> np.ndarray: streamwise_vectors = np.array([wall['Normals'][:,2], np.zeros_like(-wall['Normals'][:,0]), -wall['Normals'][:,0]]).T + else: + raise RuntimeError("'plane_normal' must be either 'xy', 'xz', or 'yz'. " + f"Instead, given {plane_normal}") # Project tangential gradient vector onto the chosen plane using n x (T_w x n) Tw = mu * np.cross(plane_normal[None,:], @@ -186,33 +200,34 @@ def compute_vorticity(dataset, scalars, vorticity_name='vorticity'): alg.Update() return pv.filters._get_output(alg) -def sampleDataBlockProfile(dataBlock, line_walldists, pointid=None, - cutterobj=None, normal=None) -> Profile: +def sampleDataBlockProfile(dataBlock: pv.MultiBlock, line_walldists: np.ndarray, + pointid: Optional[int] = None, + cutterobj: Optional[vtk.vtkPlane] = None, + normal: Optional[np.ndarray]=None) -> Profile: """Sample data block over a wall-normal profile - Given a dataBlock containing a 'grid' and 'wall' block, this will return a - PolyData object that samples 'grid' at the wall distances specified in - line_walldists. This assumes that the 'wall' block has a field named - 'Normals' containing the wall-normal vectors. + Given a dataBlock containing a ``grid`` and ``wall`` block, this will + return a PolyData object that samples ``grid`` at the wall distances + specified in line_walldists. This assumes that the ``wall`` block has a + field named ``Normals`` containing the wall-normal vectors. The location of the profile is defined by either the index of a point in - the 'wall' block or by specifying a vtk implicit function (such as - vtk.vtkPlane) that intersects the 'wall' object. The latter uses the - vtk.vtkCutter filter to determine the intersection. + the ``wall`` block or by specifying a vtk implicit function (such as + :class:`vtk.vtkPlane`) that intersects the ``wall`` object. The latter uses + the :class:`vtk.vtkCutter` filter to determine the intersection. Parameters ---------- dataBlock : pv.MultiBlock - MultiBlock containing the 'grid' and 'wall' objects + MultiBlock containing the ``grid`` and ``wall`` objects line_walldists : numpy.ndarray The locations normal to the wall that should be sampled and returned. Locations are expected to be in order. pointid : int, optional - Index of the point in 'wall' where the profile should be taken. - (default: None) + Index of the point in ``wall`` where the profile should be taken. cutterobj : vtk.vtkPlane, optional VTK object that defines the profile location via intersection with the - 'wall' + ``wall`` normal : numpy.ndarray, optional If given, use this vector as the wall normal. @@ -225,8 +240,6 @@ def sampleDataBlockProfile(dataBlock, line_walldists, pointid=None, if 'Normals' not in wall.array_names: raise RuntimeError('The wall object must have a "Normals" field present.') - if not (isinstance(pointid, int) ^ bool(cutterobj) ): #xnor - raise RuntimeError('Must provide either pointid or cutterobj.') if isinstance(pointid, int): wallnormal = wall['Normals'][pointid,:] if normal is None else normal @@ -238,8 +251,9 @@ def sampleDataBlockProfile(dataBlock, line_walldists, pointid=None, sample_line = pv.lines_from_points(sample_points) sample_line = sample_line.sample(dataBlock['grid']) sample_line['WallDistance'] = line_walldists + sample_line = Profile(sample_line) - else: + elif cutterobj: cutterout = vCutter(wall, cutterobj) if cutterout.points.shape[0] != 1: raise RuntimeError('vCutter resulted in {:d} points instead of 1.'.format( @@ -257,9 +271,13 @@ def sampleDataBlockProfile(dataBlock, line_walldists, pointid=None, sample_line = Profile(sample_line) sample_line.setWallDataFromPolyDataPoint(cutterout) + else: + raise RuntimeError('Must provide either pointid or cutterobj.') + return sample_line -def wallAlignRotationTensor(wallnormal, cart_normal, plane='xy') -> np.ndarray: +def wallAlignRotationTensor(wallnormal: np.ndarray, cart_normal: Union[list, np.ndarray], + plane: Literal['xy', 'xz', 'yz'] = 'xy') -> np.ndarray: """Create rotation tensor for wall alligning quantities For 2D xy plane and streamwise x, cart_normal should be the y unit vector. @@ -271,13 +289,16 @@ def wallAlignRotationTensor(wallnormal, cart_normal, plane='xy') -> np.ndarray: cart_normal : numpy.ndarray The Cartesian equivalent to the wall normal. If the wall were "flat", this would be the wall normal vector. - plane : str - 2D plane to rotate on, from 'xy', 'xz', 'yz'. (Default is 'xy') + plane : str, default: 'xy' + 2D plane to rotate on, from 'xy', 'xz', 'yz'. """ if plane.lower() == 'xy': rotation_axis = np.array([0, 0, 1]) elif plane.lower() == 'xz': rotation_axis = np.array([0, 1, 0]) elif plane.lower() == 'yz': rotation_axis = np.array([1, 0, 0]) + else: + raise RuntimeError("'plane' must be either 'xy', 'xz', or 'yz'. " + f"Instead, given {plane}") theta = np.arccos(np.dot(wallnormal, cart_normal)) # Determine whether clockwise or counter-clockwise rotation diff --git a/vtkpytools/bl.py b/vtkpytools/bl.py index d8ed600..6506677 100644 --- a/vtkpytools/bl.py +++ b/vtkpytools/bl.py @@ -3,14 +3,15 @@ from .barfiletools import * from scipy.integrate import cumtrapz from numpy import ndarray +import warnings def sampleAlongVectors(dataBlock, sample_dists, vectors, locations) -> pv.PolyData: - """Sample dataBlock at sample_dists away from locations in vectors direction + """Sample dataBlock at ``sample_dists`` away from ``locations`` in ``vectors`` direction - Each nth location in 'locations' corresponds with the nth vector in - 'vectors'. At each location, a sample is take at each sample_dist in - 'sample_dists' in the direction of the location's corresponding vector. + Each nth location in ``locations`` corresponds with the nth vector in + ``vectors``. At each location, a sample is take at each sample_dist in + ``sample_dists`` in the direction of the location's corresponding vector. In other words, a "profile" is taken at each location. The vector defines the direction the profile is taken (with respect to the location itself) @@ -19,25 +20,26 @@ def sampleAlongVectors(dataBlock, sample_dists, vectors, locations) -> pv.PolyDa Parameters ---------- dataBlock : pyvista.MultiBlock - dataBlock from where data will be interpolated. Must contain 'grid' and 'wall' VTK objects. - sample_dists : (S) ndarray - The distances away from location along the vectors - vectors : (L, 3) ndarray + dataBlock from where data will be interpolated. Must contain ``grid`` + and ``wall`` VTK objects. + sample_dists : ndarray + The distances away from location along the vectors. Shape: (`nsamples`) + vectors : ndarray Unit vectors that define the direction the samples should be taken away - from location. - locations : (L, 3) ndarray - Coordinates from which samples should start from + from location. Shape: (`nlocations`, 3) + locations : ndarray + Coordinates from which samples should start from. Shape: (`nlocations`, 3) Returns ------- samples: pv.PolyData - Contains L*S data points in S-major order (ie. [:S] contains the all + Contains `nlocations`*`nsamples` data points in S-major order (ie. [:S] contains the all the samples associated with location[0]) """ nlocations = locations.shape[0] - nprofilesamples = sample_dists.size - ntotsamples = nlocations*nprofilesamples + samplesperprofile = sample_dists.size + ntotsamples = nlocations*samplesperprofile profiles = np.einsum('i,jk->jik', sample_dists, vectors) profiles += locations[:, None, :] @@ -51,65 +53,67 @@ def sampleAlongVectors(dataBlock, sample_dists, vectors, locations) -> pv.PolyDa def delta_vortInt(vorticity, wall_distance, nwallpnts: int, displace=False, momentum=False, returnUvort=False, Uedge=None) -> dict: - """Calculate vorticity-integrated BL thicknesses (momentum and displacement) + r"""Calculate vorticity-integrated BL thicknesses (momentum and displacement) Based on eqns 3.1-3.4 in "Numerical study of turbulent separation bubbles with varying pressure gradient and Reynolds number" Coleman et. al. 2018 First, we define U as the cumulative integral of vorticity from the wall: - U_i(y) = \int_0^y vorticity_i(j) dj + .. math:: U(y) = \int_0^y \omega(n) dn - where i indexes the wall point and y is distance from the wall. This is - Uvort (which is returned by setting 'returnUvort' to True) and can be used - to find delta_percent. + where :math:`y` is distance from the wall. This is `Uvort` (which is + returned by setting ``returnUvort = True``) and can be used to find + `delta_percent`: + + .. code-block:: python result = vpt.delta_vortInt(...., returnUvort=True) delta_percent_vorticity = delta_percent(U=result['Uvort'], ....) - Displacement thickness, d_i, is defined as: + Displacement thickness, :math:`\delta^*`, is defined as: - d_i = -1/Uedge_i \int y * vorticity_i(y) dy + .. math:: \delta^* = -1/U_e \int y * \omega(y) dy - over the full profile height, with i indexing the wall point. + over the full profile height for each wall point. :math:`U_e` is the + velocity at the edge of the boundary layer (see `Uedge` argument). - Momentum thickness, m_i, is defined as: + Momentum thickness, :math:`\delta_\theta`, is defined as: - m_i = -2/Uedge_i^2 \int [y * U_i(y) * vorticity_i(y)]dy - d_i + .. math:: \delta_\theta = -2/U_e^2 \int [y * U(y) * \omega(y)]dy - \delta^* - over the full profile height, with i indexing the wall point. Note that U_i - is the vorticity-defined velocity defined above. + over the full profile height for each wall point. Parameters ---------- - vorticity : [nwallpnts*nprofilesamples] ndarray - Spanwise component of vorticity of the flow. - wall_distance : [nwallpnts*nprofilesamples] ndarray - Distance to wall for all the sample points + vorticity : ndarray + Spanwise component of vorticity of the flow. Shape: (`nwallpnts` * `samplesperprofile`) + wall_distance : ndarray + Distance to wall for all the sample points Shape: (`nwallpnts` * `samplesperprofile`) nwallpnts : int - Number of wall locations used in the sampling of U. The size of U and - wall_distance must be evenly divisible by nwallpnts. + Number of wall locations used in the sampling of vorticity. The size of + `vorticity` and `wall_distance` must be evenly divisible by `nwallpnts`. displace : bool, optional Whether to calculate displacement thickness momentum : bool, optional Whether to calcualte momentum thickness - Uedge : [nwallpnts] ndarray, optional + Uedge : ndarray, optional Sets the values for the edge (vorticity-integrated) velocity used in calculating the boundary layer height. If not given, the last sample point for each wall point - profile will be used (Default: None). + profile will be used. Shape: (`nwallpnts`) Returns ------- Dictionary with the following items optionally defined: - - delta_displace: ndarray, optional - Displacement thickness. Not passed if displace=False - - delta_momentum: ndarray, optional - Momentum thickness. Not passed if momentum=False - - Uvort: ndarray, optional - Vorticity-integrated velocity. Not passed if returnUvort=False + delta_displace : ndarray, optional + Displacement thickness. Not passed if ``displace = False``. Shape: + (`nwallpnts`) + delta_momentum : ndarray, optional + Momentum thickness. Not passed if ``momentum = False``. Shape: + (`nwallpnts`) + Uvort : ndarray, optional + Vorticity-integrated velocity. Not passed if ``returnUvort = False``. + Shape: (`nwallpnts`) """ @@ -141,48 +145,48 @@ def delta_vortInt(vorticity, wall_distance, nwallpnts: int, displace=False, def delta_velInt(U, wall_distance, nwallpnts: int, displace: bool = False, momentum: bool = False, Uedge=None) -> dict: - """Calculate velocity-integrated BL thicknesses (momentum and displacement) + r"""Calculate velocity-integrated BL thicknesses - Displacement thickness, d_i, defined as integrating: + Displacement thickness, :math:`\delta^*`, defined as integrating: - (1 - U_i(wall_distance)/Uedge_i) + .. math:: \int_0 (1 - U(y)/U_e) dy - over the full profile height, with i indexing the wall point. + over the full profile height for each wall point. - Momentum thickness, m_i, defined as integrating: + Momentum thickness, :math:`\delta_\theta`, defined as integrating: - (1 - U_i(wall_distance)/Uedge_i) * (U_i(wall_distance)/Uedge_i) + .. math:: \int_0 (1 - U(y)/U_e) * (U(y)/U_e) dy - over the full profile height, with i indexing the wall point. + over the full profile height for each wall point. Parameters ---------- - U : [nwallpnts*nprofilesamples] ndarray + U : ndarray Quantity to base boundary layer height on (generally streamwise - velocity) - wall_distance : [nwallpnts*nprofilesamples] ndarray - Distance to wall for all the sample points + velocity). Shape: (`nwallpnts` * `samplesperprofile`) + wall_distance : ndarray + Distance to wall for all the sample points Shape: + (`nwallpnts` * `samplesperprofile`) nwallpnts : int - Number of wall locations used in the sampling of U. The size of U and - wall_distance must be evenly divisible by nwallpnts. + Number of wall locations used in the sampling of `U`. The size of + `U` and `wall_distance` must be evenly divisible by `nwallpnts`. displace : bool, optional Whether to calculate displacement thickness momentum : bool, optional - Whether to calcualte momentum thickness - Uedge : [nwallpnts] ndarray, optional + Whether to calculate momentum thickness + Uedge : ndarray, optional Sets the values for the edge velocity used in calculating the boundary layer height. If not given, the last sample point for each wall point - profile will be used (Default: None). + profile will be used. Shape: (`nwallpnts`) Returns ------- Dictionary with the following items optionally defined: + delta_displace : ndarray, optional + Displacement thickness. Not passed if ``displace = False``. Shape: (`nwallpnts`) - delta_displace: ndarray, optional - Displacement thickness. Not passed if displace=False - - delta_momentum: ndarray, optional - Momentum thickness. Not passed if momentum=False + delta_momentum : ndarray, optional + Momentum thickness. Not passed if ``momentum = False``. Shape: (`nwallpnts`) """ if U.size % nwallpnts != 0: @@ -208,35 +212,39 @@ def delta_velInt(U, wall_distance, nwallpnts: int, def delta_percent(U, wall_distance, nwallpnts: int, percent: float, Uedge=None) -> ndarray: - """Calculate the boundary layer height based on percentage of U + r"""Calculate the boundary layer height based on percentage of `U` - Define the percent boundary layer thickness as h_i such that + Define the percent boundary layer thickness as :math:`\delta` such that - U_i(h_i) = percent*Uedge_i + .. math:: U(\delta) = percent*U_e - for U_i(wall_distance) and i indexing every wall point. + for :math:`U(y) \quad \forall y \in` `wall_distance` and :math:`U_e` the + edge velocity (see ``Uedge``). Linear interpolation between values of + :math:`y \in` `wall_distance` is done to determine :math:`\delta`. This is + repeated for each wall point. Parameters ---------- - U : [nwallpnts*nprofilesamples] ndarray + U : ndarray Quantity to base boundary layer height on (generally streamwise - velocity) - wall_distance : [nwallpnts*nprofilesamples] ndarray - Distance to wall for all the sample points + velocity). Shape: (`nwallpnts`*`samplesperprofile`) + wall_distance : ndarray + Distance to wall for all the sample points Shape: (`nwallpnts` * `samplesperprofile`) nwallpnts : int - Number of wall locations used in the sampling of U. The size of U and - wall_distance must be evenly divisible by nwallpnts. + Number of wall locations used in the sampling of `U`. The size of + `U` and `wall_distance` must be evenly divisible by `nwallpnts`. percent : float The percentage of the "edge" value that defines the boundary layer height. - Uedge : [nwallpnts] ndarray, optional + Uedge : ndarray, optional Sets the values for the edge velocity used in calculating the boundary layer height. If not given, the last sample point for each wall point - profile will be used (Default: None). + profile will be used. Shape: (`nwallpnts`) Returns ------- - delta_percent: [nwallpnts] ndarray + delta_percent : ndarray + Shape: (`nwallpnts`) """ if U.size % nwallpnts != 0: @@ -257,8 +265,17 @@ def delta_percent(U, wall_distance, nwallpnts: int, percent: float, Uedge=None) def integratedVortBLThickness(vorticity, wall_distance, delta_percent=0.995, delta_displace=False, delta_momentum=False) -> dict: - """(DEPRECATED) Computes vorticity-integrated BL thickness for a profile""" - # Use eqns 3.1-3.4 in Numerical study of turbulent separation bubbles with varying pressure gradient and Reynolds number + """(DEPRECATED) Computes vorticity-integrated BL thickness for a profile + + .. deprecated:: + Use `delta_vortInt` with `sampleAlongVectors` instead. Only kept to not + break previous scripts. + """ + + warnings.warn("This function is deprecated. Use 'delta_vortInt' with " + "'sampleAlongVectors' instead", FutureWarning) + # Use eqns 3.1-3.4 in Numerical study of turbulent separation bubbles + # with varying pressure gradient and Reynolds number U = -cumtrapz(vorticity, wall_distance, initial=0) Ue = U[-1] diff --git a/vtkpytools/common.py b/vtkpytools/common.py index 1caee2f..a591850 100644 --- a/vtkpytools/common.py +++ b/vtkpytools/common.py @@ -4,6 +4,8 @@ from scipy.io import FortranFile from pathlib import Path import re +from typing import Union, Optional +from os import PathLike def unstructuredToPoly(unstructured_grid): """Convert vtk.UnstructruedGrid to vtk.PolyData""" @@ -19,7 +21,7 @@ def orderPolyDataLine(polydata): strip.Update() return pv.wrap(strip.GetOutput()) -def vCutter(input_data, cut_function): +def vCutter(input_data: pv.PolyData, cut_function): """Returns the intersection of input_data and cut_function Wrapper around vtkCutter filter. Output contains interpolated data from @@ -48,13 +50,24 @@ class Profile(pv.PolyData): Use case is for storage of wall local data (boundary layer metrics, Cf etc.) with profiles that correspond to that wall local data. + Attributes + ---------- + walldata : dict + Dictionary of data from the wall corresponding to the profile location. + Generated using `setWallDataFromPolyDataPoint` """ walldata = {} - def setWallDataFromPolyDataPoint(self, PolyPoint): + def setWallDataFromPolyDataPoint(self, PolyPoint: pv.PolyData): """Set walldata attribute from PolyData Point Primary use case is the using the output of vtkpytools.vCutter() + + Parameters + ---------- + PolyPoint : pv.PolyPoint + PolyData object containing one point. The point arrays attached to + the points are converted to a `dict` and set to `walldata`. """ if PolyPoint.n_points != 1: raise RuntimeError('Profile should only have 1 wallpoint, {:d} given'.format( @@ -62,12 +75,11 @@ def setWallDataFromPolyDataPoint(self, PolyPoint): self.walldata = dict(PolyPoint.point_arrays) self.walldata['Point'] = PolyPoint.points -def readBinaryArray(path, ncols) -> np.ndarray: +def readBinaryArray(path: Union[str, PathLike], ncols: int) -> np.ndarray: """Get array from Fortran binary file. Parameters ---------- - path : Path Path to Fortran binary array. ncols : uint @@ -80,7 +92,7 @@ def readBinaryArray(path, ncols) -> np.ndarray: return array -def globFile(globstring, path: Path, regex=False) -> Path: +def globFile(globstring, path: Path, regex=False) -> Optional[Path]: """ Glob for one file in directory, then return. If it finds more than one file matching the globstring, it will error out. diff --git a/vtkpytools/gridtools2d/core.py b/vtkpytools/gridtools2d/core.py index 95b163f..58f4050 100644 --- a/vtkpytools/gridtools2d/core.py +++ b/vtkpytools/gridtools2d/core.py @@ -79,20 +79,22 @@ def form2DGrid(coords_array, connectivity_array=None) -> pv.UnstructuredGrid: return grid def computeEdgeNormals(edges, domain_point) -> pv.PolyData: - """Compute the normals of the edge assuming coplanar to XY plane + r"""Compute the normals of the edge assuming coplanar to XY plane Loops through every line (or cell) in the edges to calculate it's normal vector. Then it will ensure that the vectors face towards the inside of the domain using the following: - 1) Create vector from point on line segment to domain_point called "insideVector" + #. Create vector from point on line segment to ``domain_point`` + :math:`\mathbf{a}` - 2) Determine whether the current normal vector points in the direction as - the insideVector using a dot product. + #. Determine whether the current normal vector points in the direction as + the :math:`\mathbf{a}` using a dot product. - 3) Reverse the wall normal vector if it points outside the domain + #. Reverse the wall normal vector if it points outside the domain - n = n * dot(n, insideVector / |insideVector|) + .. math:: \mathbf{n} = \mathbf{n} * \frac{\mathbf{n} \cdot + \mathbf{a}}{\vert \mathbf{a} \vert} Parameters ---------- @@ -104,7 +106,7 @@ def computeEdgeNormals(edges, domain_point) -> pv.PolyData: """ normals = np.zeros((edges.n_cells, 3)) - i = 0 + # Indices of 2 points forming line cell indices = edges.lines.reshape((edges.n_cells, 3))[:,1:3] pnts1 = edges.points[indices[:,0], :] @@ -112,7 +114,7 @@ def computeEdgeNormals(edges, domain_point) -> pv.PolyData: normals[:, 0:2] = np.array([-(pnts1[:,1] - pnts2[:,1]), (pnts1[:,0] - pnts2[:,0])]).T - inside_vector = domain_point - pnts1 + inside_vector = domain_point - pnts1 # \mathbf{a} vector # Dot product of all the normal vectors with the inside-pointing vector inOrOut = np.einsum('ij,ij->i',normals,inside_vector) diff --git a/vtkpytools/gridtools3d/core.py b/vtkpytools/gridtools3d/core.py index a4c8435..9de6002 100644 --- a/vtkpytools/gridtools3d/core.py +++ b/vtkpytools/gridtools3d/core.py @@ -1,6 +1,5 @@ import numpy as np import vtk -from pathlib import Path import pyvista as pv @@ -13,12 +12,12 @@ def form3DGrid(coords_array, connectivity_array) -> pv.UnstructuredGrid: Parameters ---------- coords_array : numpy.ndarray - Coordinates of the points. Shape in (nPoints,2). + Coordinates of the points. Shape: (nPoints, 2) connectivity_array : numpy.ndarray - Connectivity array of the cells. Shape is (nCells,PointsPerCell). The - index of points should start from 0. The type of cell will be inferred - based on the number of points per cell. Mixed meshes inferred by - repeated node IDs in a single element. + Connectivity array of the cells. The index of points should start from + 0. The type of cell will be inferred based on the number of points per + cell. Mixed meshes inferred by repeated node IDs in a single element. + Shape: (nCells, PointsPerCell) Returns ------- diff --git a/vtkpytools/misc/__init__.py b/vtkpytools/misc/__init__.py new file mode 100644 index 0000000..683617c --- /dev/null +++ b/vtkpytools/misc/__init__.py @@ -0,0 +1 @@ +from . import spectra diff --git a/vtkpytools/misc/spectra.py b/vtkpytools/misc/spectra.py new file mode 100644 index 0000000..b57cfea --- /dev/null +++ b/vtkpytools/misc/spectra.py @@ -0,0 +1,46 @@ +import numpy as np +import re +from typing import Tuple, Optional, Sequence, Dict +from pathlib import Path + + +def energySpectra1D(u: np.ndarray, L: float, periodic: bool = True) -> Tuple[np.ndarray, np.ndarray]: + r"""Calculate energy spectra for a 1D signal + + Given :math:`u_n,\ n=0\...N-1` samples of a signal, return the energy + associated a frequency (bin) :math:`k`. This is done using a + real-restricted FFT. + + .. math:: \int_0^\pi + + Parameters + ---------- + u : np.ndarray + Array of signal(s). If a 2D array is given, then the signals are + assumed to be continuous in the last index. ie. ``u[0, :]`` is a vector + of a continuous signal. + L : float + This is the length scale used to scale the frequency values + periodic : bool + Whether the continuous signal samples contain the periodic endpoint. + + Returns + ------- + Tuple[np.ndarray, np.ndarray] + + """ + if periodic: u = u[..., 0:-1] # Remove periodic data point + + n = u.shape[-1] + + utrans = np.fft.rfft(u)/n + energy = 0.5*(utrans*np.conj(utrans)) + + ks = np.fft.rfftfreq(n, L)*2*np.pi + + # double the contributions of non-zero wavemodes + # rfft gives only half of the spectrum (for real input, it is symmetric + # about k_0) + energy[1:] *= 2 + + return ks, energy.astype(float) diff --git a/vtkpytools/numtools.py b/vtkpytools/numtools.py index dcf1fd0..bc40a16 100644 --- a/vtkpytools/numtools.py +++ b/vtkpytools/numtools.py @@ -129,26 +129,36 @@ def seriesDiffLimiter(series: np.ndarray, dx=None, magnitude=None) -> np.ndarray def symmetric2FullTensor(tensor_array) -> np.ndarray: - """ Turn (n, 6) shape array of tensor entries into (n, 3, 3) + """Turn (..., 6) shape array of tensor entries into (..., 3, 3) Assumed that symmtetric entires are in XX YY ZZ XY XZ YZ order.""" + if tensor_array.shape[-1] != 6: + raise ValueError('The array is not in symmetric form! ' + f'The last axis must be size 6, not {tensor_array.shape[-1]}.') + shaped_tensors = np.array([ - tensor_array[:,0], tensor_array[:,3], tensor_array[:,4], - tensor_array[:,3], tensor_array[:,1], tensor_array[:,5], - tensor_array[:,4], tensor_array[:,5], tensor_array[:,2] + tensor_array[...,0], tensor_array[...,3], tensor_array[...,4], + tensor_array[...,3], tensor_array[...,1], tensor_array[...,5], + tensor_array[...,4], tensor_array[...,5], tensor_array[...,2] ]).T - return shaped_tensors.reshape(shaped_tensors.shape[0], 3, 3) + return shaped_tensors.reshape(*tensor_array.shape[:-1], 3, 3) def full2SymmetricTensor(tensor_array) -> np.ndarray: - """ Turn (n, 3, 3) shape array of tensor entries into (n, 6) + """Turn (..., 3, 3) shape array of tensor entries into (..., 6) - Symmetric entires are in XX YY ZZ XY XZ YZ order.""" - return np.array([ - tensor_array[:,0,0], tensor_array[:,1,1], tensor_array[:,2,2], - tensor_array[:,0,1], tensor_array[:,0,2], tensor_array[:,1,2] - ]).T + Symmetric entires are in XX YY ZZ XY XZ YZ order. + """ + if tensor_array.shape[-2:] != (3,3): + raise ValueError('The array is not in full form! ' + f'The last two axis must be shape (3,3), not {tensor_array.shape[-2:]}.') + + shaped_tensors = np.empty((*tensor_array.shape[:-2], 6)) + indices = [(0,0), (1,1), (2,2), (0,1), (0,2), (1,2)] + for i, ind in enumerate(indices): + shaped_tensors[...,i] = tensor_array[...,ind[0], ind[1]] + return shaped_tensors def makeRotationTensor(rotation_axis, theta) -> np.ndarray: