diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml deleted file mode 100644 index d3686289..00000000 --- a/.github/workflows/ci.yml +++ /dev/null @@ -1,50 +0,0 @@ -name: "CI" - -on: push - -jobs: - build: - runs-on: ${{ matrix.os }} - strategy: - matrix: - os: [ubuntu-latest] - python-version: ["3.10"] - - steps: - - uses: actions/checkout@v4 - - - name: Set up Python - uses: actions/setup-python@v5 - with: - python-version: ${{ matrix.python-version }} - - - name: Setup Conda - uses: s-weigand/setup-conda@v1 - with: - python-version: ${{ matrix.python-version }} - update-conda: true - conda-channels: conda-forge, anaconda - - - name: Install pythonocc-core (Conda) - run: conda install --yes -c conda-forge -c anaconda pythonocc-core - shell: bash - - - name: Install pip dependencies - run: | - python -m pip install --upgrade pip - python -m pip install numpy scipy matplotlib pytest pytest-cov - python -m pip install smithers[vtk] - python -m pip install .[test] - - - name: Run tests with pytest - env: - CODACY_API_TOKEN: ${{ secrets.CODACY_API_TOKEN }} - run: | - if [ -z "$CODACY_API_TOKEN" ]; then - python -m pytest - else - python -m pytest --cov-report term --cov-report xml:cobertura.xml --cov=pygem - curl -s https://coverage.codacy.com/get.sh -o CodacyCoverageReporter.sh - chmod +x CodacyCoverageReporter.sh - ./CodacyCoverageReporter.sh report -r cobertura.xml -t $CODACY_API_TOKEN - fi \ No newline at end of file diff --git a/.github/workflows/cd.yml b/.github/workflows/deployer.yml similarity index 58% rename from .github/workflows/cd.yml rename to .github/workflows/deployer.yml index 81bd0825..aefb7483 100644 --- a/.github/workflows/cd.yml +++ b/.github/workflows/deployer.yml @@ -1,4 +1,4 @@ -name: "Continuous Deployment" +name: "Deployer" on: push: @@ -7,8 +7,6 @@ on: jobs: - ############################################################################# - # Create and online deployment of the documentation ######################### docs: ####################################################################### runs-on: ubuntu-latest steps: @@ -21,13 +19,14 @@ jobs: python-version: "3.10" - name: Install Dependencies (conda and pip) - shell: bash + shell: bash -el {0} run: | - conda install --yes -c conda-forge pythonocc-core + conda install --yes -c conda-forge pythonocc-core pip python -m pip install --upgrade pip python -m pip install .[docs] - name: Build Documentation + shell: bash -el {0} run: | make html working-directory: docs/ @@ -38,9 +37,7 @@ jobs: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: ./docs/build/html allow_empty_commit: true - - ############################################################################# - ## Create a public "Release" on the Github page ############################# + release_github: ############################################################# runs-on: ubuntu-latest permissions: @@ -50,3 +47,27 @@ jobs: - uses: ncipollo/release-action@v1 with: token: ${{ secrets.GITHUB_TOKEN }} + + pypi: ####################################################################### + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Setup Python + uses: actions/setup-python@v5 + with: + python-version: "3.10" + + - name: Install build + run: >- + python -m pip install build + + - name: Build a binary wheel and a source tarball + run: >- + python -m build --sdist --wheel --outdir dist/ . + + - name: Publish distribution to PyPI + if: startsWith(github.ref, 'refs/tags') + uses: pypa/gh-action-pypi-publish@release/v1 + with: + password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/master_cleaner.yml b/.github/workflows/master_cleaner.yml new file mode 100644 index 00000000..66106786 --- /dev/null +++ b/.github/workflows/master_cleaner.yml @@ -0,0 +1,29 @@ +name: Master Cleaner + +on: + push: + branches: + - master + +jobs: + formatter: + name: runner / black + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - uses: psf/black@stable + with: + src: "./pygem" + + - name: Create Pull Request + uses: peter-evans/create-pull-request@v3 + with: + token: ${{ secrets.GITHUB_TOKEN }} + title: "Format Python code with psf/black push" + commit-message: ":art: Format Python code with psf/black" + body: | + There appear to be some python formatting errors in ${{ github.sha }}. This pull request + uses the [psf/black](https://github.com/psf/black) formatter to fix these issues. + base: ${{ github.head_ref }} + branch: actions/black diff --git a/.github/workflows/monthly-tagger.yml b/.github/workflows/monthly-tagger.yml new file mode 100644 index 00000000..299794b4 --- /dev/null +++ b/.github/workflows/monthly-tagger.yml @@ -0,0 +1,49 @@ +name: "Monthly Tagger" + +on: + schedule: + - cron: '20 2 1 * *' + +jobs: + + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [windows-latest, macos-latest, ubuntu-latest] + python-version: ["3.9", "3.10", "3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: ${{ matrix.python-version }} + channels: conda-forge + - name: Install Python dependencies + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core + python -m pip install --upgrade pip + python -m pip install .[test] + - name: Test with pytest + shell: bash -el {0} + run: | + python -m pytest + + monthly_tag: + runs-on: ubuntu-latest + needs: test + steps: + - uses: actions/checkout@v4 + with: + token: ${{ secrets.PAT_PYGEM_PUSH }} + + - name: Create and push the tag + run: | + VERS="post$(date +%y%m)" + git config --global user.name 'Monthly Tag bot' + git config --global user.email 'mtbot@noreply.github.com' + git commit --allow-empty -m "monthly version $VERS" + git tag -a "v$VERS" -m "Monthly version $VERS" + git push origin "v$VERS" diff --git a/.github/workflows/tester.yml b/.github/workflows/tester.yml new file mode 100644 index 00000000..fee0afa4 --- /dev/null +++ b/.github/workflows/tester.yml @@ -0,0 +1,111 @@ +name: "Testing Pull Request" + +on: + pull_request: + branches: + - "master" + - "dev" + +jobs: + unittests: ################################################################# + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [windows-latest, macos-latest, ubuntu-latest] + python-version: ["3.9", "3.10", "3.11", "3.12"] + steps: + - uses: actions/checkout@v4 + + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: ${{ matrix.python-version }} + channels: conda-forge + + - name: Installing packages + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core pip + python -m pip install --upgrade pip + python -m pip install .[test] + + - name: Test with pytest + shell: bash -el {0} + run: | + python -m pytest + + linter: #################################################################### + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Run Black formatter (check mode) + uses: psf/black@stable + with: + src: "./pygem" + + testdocs: ################################################################## + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: "3.10" + channels: conda-forge + + - name: Install Python dependencies + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core pip + python -m pip install --upgrade pip + python -m pip install .[docs] + + - name: Build Documentation + shell: bash -el {0} + run: | + make html SPHINXOPTS+='-W' + working-directory: docs/ + + coverage: ################################################################## + runs-on: ubuntu-latest + permissions: + contents: write + pull-requests: write + + steps: + - uses: actions/checkout@v4 + + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: "3.10" + channels: conda-forge + + - name: Install Python dependencies + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core pip + python -m pip install --upgrade pip + python -m pip install pytest pytest-cov + python -m pip install .[test] + + - name: Generate coverage report + shell: bash -el {0} + run: | + python -m pytest --cov-report term --cov-report xml:cobertura.xml --cov=pygem + + - name: Produce the coverage report + if: github.event.pull_request.head.repo.full_name == github.repository + uses: insightsengineering/coverage-action@v2 + with: + path: ./cobertura.xml + threshold: 80.0 + fail: true + publish: ${{ github.event.pull_request.head.repo.full_name == github.repository }} + coverage-summary-title: "Code Coverage Summary" diff --git a/.github/workflows/testing_pr.yml b/.github/workflows/testing_pr.yml deleted file mode 100644 index 80f5196b..00000000 --- a/.github/workflows/testing_pr.yml +++ /dev/null @@ -1,38 +0,0 @@ -name: "Testing Pull Request" - -on: - pull_request: - branches: - - "master" - -jobs: - build: - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - os: [windows-latest, macos-latest, ubuntu-latest] - python-version: ["3.9", "3.10", "3.11", "3.12"] - steps: - - uses: actions/checkout@v2 - - - - name: Installing conda - uses: conda-incubator/setup-miniconda@v3 - - with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - channels: conda-forge - - - name: Installing packages - shell: bash -el {0} - run: | - conda install --yes -c conda-forge pythonocc-core - python -m pip install --upgrade pip - python -m pip install .[test] - - - name: Test with pytest - shell: bash -el {0} - run: | - python -m pytest diff --git a/.github/workflows/tutorial_exporter.yml b/.github/workflows/tutorial_exporter.yml new file mode 100644 index 00000000..547e0f6f --- /dev/null +++ b/.github/workflows/tutorial_exporter.yml @@ -0,0 +1,143 @@ +name: "Export Tutorials" + +on: + workflow_dispatch: + push: + branches: + - "dev" + - "master" + paths: + - 'tutorials/**/*.ipynb' + +jobs: + # run on push + export_tutorials_on_push: + if: ${{ github.event_name == 'push' }} + permissions: write-all + runs-on: ubuntu-latest + env: + TUTORIAL_TIMEOUT: 1200s + steps: + - uses: actions/checkout@v4 + + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: "3.10" + channels: conda-forge + + - name: Install dependencies + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core pip + python -m pip install --upgrade pip -e .[tut] black[jupyter] ipykernel + + - name: Setup FFmpeg + uses: FedericoCarboni/setup-ffmpeg@v2 + + - id: files + uses: jitterbit/get-changed-files@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + format: space-delimited + + - name: Configure git + run: | + git config user.name "github-actions[bot]" + git config user.email 41898282+github-actions[bot]@users.noreply.github.com + + - name: Run formatter + shell: bash -el {0} + run: black tutorials/ + + - name: Export tutorials to .py and .html + shell: bash -el {0} + run: | + set -x + for file in ${{ steps.files.outputs.all }}; do + if [[ $file == *.ipynb ]]; then + filename=$(basename $file) + pyfilename=$(echo ${filename%?????})py + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert $file --to python --output $pyfilename --output-dir=$(dirname $file) + htmlfilename=$(echo ${filename%?????} | sed -e 's/-//g')html + htmldir="docs/source"/$(echo ${file%??????????????} | sed -e 's/-//g') + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert --execute $file --to html --output $htmlfilename --output-dir=$htmldir + fi + done + set +x + + - uses: benjlevesque/short-sha@v2.1 + id: short-sha + + - name: Create Pull Request + uses: peter-evans/create-pull-request@v5.0.2 + with: + labels: maintenance + title: Export tutorial changed in ${{ steps.short-sha.outputs.sha }} + branch: export-tutorial-${{ steps.short-sha.outputs.sha }} + base: ${{ github.head_ref }} + commit-message: export tutorials changed in ${{ steps.short-sha.outputs.sha }} + delete-branch: true + + # run on workflow_dispatch + export_tutorials_workflow_dispatch: + if: ${{ github.event_name == 'workflow_dispatch' }} + permissions: write-all + runs-on: ubuntu-latest + env: + TUTORIAL_TIMEOUT: 1200s + steps: + - uses: actions/checkout@v4 + + - name: Installing conda + uses: conda-incubator/setup-miniconda@v3 + with: + auto-update-conda: true + python-version: "3.10" + channels: conda-forge + + - name: Install dependencies + shell: bash -el {0} + run: | + conda install --yes -c conda-forge pythonocc-core pip + python -m pip install --upgrade pip -e .[tut] black[jupyter] ipykernel + + - name: Setup FFmpeg + uses: FedericoCarboni/setup-ffmpeg@v2 + + - name: Configure git + run: | + git config user.name "github-actions[bot]" + git config user.email 41898282+github-actions[bot]@users.noreply.github.com + + - name: Run formatter + shell: bash -el {0} + run: black tutorials/ + + - name: Export all tutorials to .py and .html + shell: bash -el {0} + run: | + set -x + for file in $(find tutorials -type f -name "*.ipynb"); do + filename=$(basename $file) + pyfilename="${filename%.ipynb}.py" + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert $file --to python --output $pyfilename --output-dir=$(dirname $file) + htmlfilename="${filename%.ipynb}.html" + htmldir="docs/source"/$(dirname $file) + timeout --signal=SIGKILL $TUTORIAL_TIMEOUT python -Xfrozen_modules=off -m jupyter nbconvert --execute $file --to html --output $htmlfilename --output-dir=$htmldir + done + set +x + + - uses: benjlevesque/short-sha@v2.1 + id: short-sha + + - name: Create Pull Request + uses: peter-evans/create-pull-request@v5.0.2 + with: + labels: maintenance + title: Export tutorial changed in ${{ steps.short-sha.outputs.sha }} + branch: export-tutorial-${{ steps.short-sha.outputs.sha }} + base: ${{ github.head_ref }} + commit-message: export tutorials changed in ${{ steps.short-sha.outputs.sha }} + delete-branch: true diff --git a/pygem/__init__.py b/pygem/__init__.py index f8bf8fad..c6d94388 100644 --- a/pygem/__init__.py +++ b/pygem/__init__.py @@ -1,6 +1,7 @@ """ PyGeM init """ + __all__ = [ "deformation", "ffd", @@ -8,9 +9,7 @@ "idw", "rbf_factory", "custom_deformation", - "cffd" - "bffd" - "vffd" + "cffd" "bffd" "vffd", ] from .deformation import Deformation @@ -22,4 +21,4 @@ from .custom_deformation import CustomDeformation from .meta import * from .bffd import BFFD -from .vffd import VFFD \ No newline at end of file +from .vffd import VFFD diff --git a/pygem/bffd.py b/pygem/bffd.py index 4a86cefb..bb72bb9f 100644 --- a/pygem/bffd.py +++ b/pygem/bffd.py @@ -1,12 +1,14 @@ from pygem.cffd import CFFD import numpy as np + + class BFFD(CFFD): - ''' + """ Class that handles the Barycenter Free Form Deformation on the mesh points. - + :param list n_control_points: number of control points in the x, y, and z direction. Default is [2, 2, 2]. - + :cvar numpy.ndarray box_length: dimension of the FFD bounding box, in the x, y and z direction (local coordinate system). :cvar numpy.ndarray box_origin: the x, y and z coordinates of the origin of @@ -21,8 +23,8 @@ class BFFD(CFFD): z, normalized with the box length z. :cvar callable fun: it defines the F of the constraint F(x)=c. Default is the constant 1 function. :cvar numpy.ndarray fixval: it defines the c of the constraint F(x)=c. Default is 1. - :cvar numpy.ndarray mask: a boolean tensor that tells to the class - which control points can be moved, and in what direction, to enforce the constraint. + :cvar numpy.ndarray mask: a boolean tensor that tells to the class + which control points can be moved, and in what direction, to enforce the constraint. The tensor has shape (n_x,n_y,n_z,3), where the last dimension indicates movement on x,y,z respectively. Default is all true. @@ -36,7 +38,7 @@ class BFFD(CFFD): >>> bffd.adjust_control_points(original_mesh_points[:-4]) >>> assert np.isclose(np.linalg.norm(bffd.fun(bffd.ffd(original_mesh_points[:-4])) - b), np.array([0.])) >>> new_mesh_points = bffd.ffd(original_mesh_points) - ''' + """ def __init__(self, fixval=None, n_control_points=None, ffd_mask=None): super().__init__(fixval, None, n_control_points, ffd_mask, None) @@ -46,7 +48,6 @@ def linfun(x): self.fun = linfun self.fixval = fixval - self.fun_mask = np.array([[True, False, False], [False, True, False], - [False, False, True]]) - - + self.fun_mask = np.array( + [[True, False, False], [False, True, False], [False, False, True]] + ) diff --git a/pygem/cad/__init__.py b/pygem/cad/__init__.py index f7d46a10..4b15c271 100644 --- a/pygem/cad/__init__.py +++ b/pygem/cad/__init__.py @@ -5,10 +5,10 @@ try: import OCC except ModuleNotFoundError as e: - print('\nOCC not found, but required to deal with CAD files') - print('Install it using:') - print('\tconda install -c conda-forge pythonocc-core=7.4.0') - print('or visit https://github.com/tpaviot/pythonocc-core for more info\n') + print("\nOCC not found, but required to deal with CAD files") + print("Install it using:") + print("\tconda install -c conda-forge pythonocc-core=7.4.0") + print("or visit https://github.com/tpaviot/pythonocc-core for more info\n") raise e diff --git a/pygem/cad/cad_deformation.py b/pygem/cad/cad_deformation.py index b527e95e..83496596 100644 --- a/pygem/cad/cad_deformation.py +++ b/pygem/cad/cad_deformation.py @@ -1,38 +1,53 @@ """ Module for generic deformation for CAD file. """ + import os import numpy as np from itertools import product -from OCC.Core.TopoDS import (TopoDS_Shape, TopoDS_Compound, - TopoDS_Face, TopoDS_Wire) +from OCC.Core.TopoDS import ( + TopoDS_Shape, + TopoDS_Compound, + TopoDS_Face, + TopoDS_Wire, +) from OCC.Core.BRep import BRep_Builder from OCC.Core.TopExp import TopExp_Explorer -from OCC.Core.TopAbs import (TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE) +from OCC.Core.TopAbs import TopAbs_EDGE, TopAbs_FACE, TopAbs_WIRE from OCC.Core.TopTools import TopTools_ListOfShape -from OCC.Core.BRepBuilderAPI import (BRepBuilderAPI_MakeFace, - BRepBuilderAPI_MakeWire, - BRepBuilderAPI_MakeEdge, - BRepBuilderAPI_NurbsConvert) +from OCC.Core.BRepBuilderAPI import ( + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_MakeWire, + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_NurbsConvert, +) from OCC.Core.BRep import BRep_Tool from OCC.Core.Geom import Geom_BSplineCurve, Geom_BSplineSurface -from OCC.Core.GeomConvert import (geomconvert_SurfaceToBSplineSurface, - geomconvert, - GeomConvert_CompCurveToBSplineCurve) +from OCC.Core.GeomConvert import ( + geomconvert_SurfaceToBSplineSurface, + geomconvert, + GeomConvert_CompCurveToBSplineCurve, +) from OCC.Core.gp import gp_Pnt from OCC.Core.BRepTools import breptools from OCC.Core.IFSelect import IFSelect_RetDone from OCC.Core.Interface import Interface_Static_SetCVal -from OCC.Core.STEPControl import (STEPControl_Writer, STEPControl_Reader, - STEPControl_AsIs) -from OCC.Core.IGESControl import (IGESControl_Writer, IGESControl_Reader, - IGESControl_Controller) - - -class CADDeformation(): +from OCC.Core.STEPControl import ( + STEPControl_Writer, + STEPControl_Reader, + STEPControl_AsIs, +) +from OCC.Core.IGESControl import ( + IGESControl_Writer, + IGESControl_Reader, + IGESControl_Controller, +) + + +class CADDeformation: """ Base class for applting deformation to CAD geometries. - + :param int u_knots_to_add: the number of knots to add to the NURBS surfaces along `u` direction before the deformation. This parameter is useful whenever the gradient of the imposed deformation present spatial scales @@ -63,9 +78,9 @@ class CADDeformation(): default value only if the input file scale is significantly different form mm, making some of the aforementioned operations fail. Default is 0.0001. - + :cvar int u_knots_to_add: the number of knots to add to the NURBS surfaces - along `u` direction before the deformation. + along `u` direction before the deformation. :cvar int v_knots_to_add: the number of knots to add to the NURBS surfaces along `v` direction before the deformation. :cvar int t_knots_to_add: the number of knots to add to the NURBS curves @@ -74,11 +89,14 @@ class CADDeformation(): operations of the procedure (joining wires in a single curve before deformation and placing new poles on curves and surfaces). """ - def __init__(self, - u_knots_to_add=0, - v_knots_to_add=0, - t_knots_to_add=0, - tolerance=1e-4): + + def __init__( + self, + u_knots_to_add=0, + v_knots_to_add=0, + t_knots_to_add=0, + tolerance=1e-4, + ): self.u_knots_to_add = u_knots_to_add self.v_knots_to_add = v_knots_to_add self.t_knots_to_add = t_knots_to_add @@ -87,27 +105,27 @@ def __init__(self, @staticmethod def read_shape(filename): """ - Static method to load the `TopoDS_Shape` from a file. - Supported extensions are: ".iges", ".step". - - :param str filename: the name of the file containing the geometry. - - Example: - - >>> from pygem.cad import CADDeformation as CAD - >>> shape = CAD.read_shape('example.iges') - """ + Static method to load the `TopoDS_Shape` from a file. + Supported extensions are: ".iges", ".step". + + :param str filename: the name of the file containing the geometry. + + Example: + + >>> from pygem.cad import CADDeformation as CAD + >>> shape = CAD.read_shape('example.iges') + """ file_extension = os.path.splitext(filename)[1] av_readers = { - '.step': STEPControl_Reader, - '.iges': IGESControl_Reader, + ".step": STEPControl_Reader, + ".iges": IGESControl_Reader, } reader_class = av_readers.get(file_extension) if reader_class is None: - raise ValueError('Unable to open the input file') + raise ValueError("Unable to open the input file") reader = reader_class() return_reader = reader.ReadFile(filename) @@ -126,52 +144,53 @@ def read_shape(filename): @staticmethod def write_shape(filename, shape): """ - Static method to save a `TopoDS_Shape` object to a file. - Supported extensions are: ".iges", ".step". - - :param str filename: the name of the file where the shape will be saved. - - Example: - - >>> from pygem.cad import CADDeformation as CAD - >>> CAD.read_shape('example.step', my_shape) - """ + Static method to save a `TopoDS_Shape` object to a file. + Supported extensions are: ".iges", ".step". + + :param str filename: the name of the file where the shape will be saved. + + Example: + + >>> from pygem.cad import CADDeformation as CAD + >>> CAD.read_shape('example.step', my_shape) + """ + def write_iges(filename, shape): - """ IGES writer """ + """IGES writer""" IGESControl_Controller.Init() writer = IGESControl_Writer() writer.AddShape(shape) writer.Write(filename) def write_step(filename, shape): - """ STEP writer """ + """STEP writer""" step_writer = STEPControl_Writer() # Changes write schema to STEP standard AP203 # It is considered the most secure standard for STEP. # *According to PythonOCC documentation (http://www.pythonocc.org/) Interface_Static_SetCVal("write.step.schema", "AP203") - Interface_Static_SetCVal('write.surfacecurve.mode','0') + Interface_Static_SetCVal("write.surfacecurve.mode", "0") step_writer.Transfer(shape, STEPControl_AsIs) step_writer.Write(filename) file_extension = os.path.splitext(filename)[1] av_writers = { - '.step': write_step, - '.iges': write_iges, + ".step": write_step, + ".iges": write_iges, } writer = av_writers.get(file_extension) if writer is None: - raise ValueError('Unable to open the output file') + raise ValueError("Unable to open the output file") writer(filename, shape) def _bspline_surface_from_face(self, face): """ Private method that takes a TopoDS_Face and transforms it into a Bspline_Surface. - - :param TopoDS_Face face: the TopoDS_Face to be converted + + :param TopoDS_Face face: the TopoDS_Face to be converted :rtype: Geom_BSplineSurface """ if not isinstance(face, TopoDS_Face): @@ -233,7 +252,7 @@ def _enrich_curve_knots(self, bsp_curve): """ Private method that adds `self.t_knots_to_add` poles to the passed curve. - + :param Geom_BSplineCurve bsp_curve: the NURBS curve to enrich """ if not isinstance(bsp_curve, Geom_BSplineCurve): @@ -245,9 +264,12 @@ def _enrich_curve_knots(self, bsp_curve): # end parameter of composite curve last_param = bsp_curve.LastParameter() for i in range(self.t_knots_to_add): - bsp_curve.InsertKnot(first_param+ \ - i*(last_param-first_param)/self.t_knots_to_add, 1, \ - self.tolerance) + bsp_curve.InsertKnot( + first_param + + i * (last_param - first_param) / self.t_knots_to_add, + 1, + self.tolerance, + ) def _enrich_surface_knots(self, bsp_surface): """ @@ -265,18 +287,24 @@ def _enrich_surface_knots(self, bsp_surface): bounds = bsp_surface.Bounds() for i in range(self.u_knots_to_add): - bsp_surface.InsertUKnot(bounds[0]+ \ - i*(bounds[1]-bounds[0])/self.u_knots_to_add, 1, self.tolerance) + bsp_surface.InsertUKnot( + bounds[0] + i * (bounds[1] - bounds[0]) / self.u_knots_to_add, + 1, + self.tolerance, + ) for i in range(self.v_knots_to_add): - bsp_surface.InsertVKnot(bounds[2]+ \ - i*(bounds[3]-bounds[2])/self.v_knots_to_add, 1, self.tolerance) + bsp_surface.InsertVKnot( + bounds[2] + i * (bounds[3] - bounds[2]) / self.v_knots_to_add, + 1, + self.tolerance, + ) def _pole_get_components(self, pole): - """ Extract component from gp_Pnt """ + """Extract component from gp_Pnt""" return pole.X(), pole.Y(), pole.Z() def _pole_set_components(self, components): - """ Return a gp_Pnt with the passed components """ + """Return a gp_Pnt with the passed components""" assert len(components) == 3 return gp_Pnt(*components) @@ -349,7 +377,7 @@ def __call__(self, obj, dst=None): geometry. If `obj` is a filename, the method deforms the geometry contained in the file and writes the deformed shape to `dst` (which has to be set). - + :param obj: the input geometry. :type obj: str or TopoDS_Shape :param str dst: if `obj` is a string containing the input filename, @@ -359,7 +387,8 @@ def __call__(self, obj, dst=None): if isinstance(obj, str): # if a input filename is passed if dst is None: raise ValueError( - 'Source file is provided, but no destination specified') + "Source file is provided, but no destination specified" + ) shape = self.read_shape(obj) elif isinstance(obj, TopoDS_Shape): shape = obj @@ -367,7 +396,7 @@ def __call__(self, obj, dst=None): else: raise TypeError - #create compound to store modified faces + # create compound to store modified faces compound_builder = BRep_Builder() compound = TopoDS_Compound() compound_builder.MakeCompound(compound) @@ -397,8 +426,8 @@ def __call__(self, obj, dst=None): # curves (actually, the WIRES) that define the bounds of the # surface and TRIM the surface with them, to obtain the new FACE - #we now start really looping on the wires - #we will create a single curve joining all the edges in the wire + # we now start really looping on the wires + # we will create a single curve joining all the edges in the wire # the curve must be a bspline curve so we need to make conversions # through all the way @@ -431,8 +460,9 @@ def __call__(self, obj, dst=None): # edge (to be converted to wire) obtained from the modified # Bspline curve - modified_composite_edge = \ - BRepBuilderAPI_MakeEdge(composite_curve).Edge() + modified_composite_edge = BRepBuilderAPI_MakeEdge( + composite_curve + ).Edge() # modified edge is added to shapes_list shapes_list.Append(modified_composite_edge) @@ -457,8 +487,9 @@ def __call__(self, obj, dst=None): # so once we finished looping on all the wires to modify them, # we first use the only outer one to trim the surface # face builder object - face_maker = BRepBuilderAPI_MakeFace(bspline_surface, - outer_wires[0]) + face_maker = BRepBuilderAPI_MakeFace( + bspline_surface, outer_wires[0] + ) # and then add all other inner wires for the holes for inner_wire in inner_wires: diff --git a/pygem/cad/custom_deformation.py b/pygem/cad/custom_deformation.py index 09e46fe4..c6d8910f 100644 --- a/pygem/cad/custom_deformation.py +++ b/pygem/cad/custom_deformation.py @@ -6,6 +6,7 @@ from pygem import CustomDeformation as OriginalCustomDeformation from .cad_deformation import CADDeformation + class CustomDeformation(CADDeformation, OriginalCustomDeformation): """ Class to perform a custom deformation to the CAD geometries. @@ -45,17 +46,17 @@ class CustomDeformation(CADDeformation, OriginalCustomDeformation): default value only if the input file scale is significantly different form mm, making some of the aforementioned operations fail. Default is 0.0001. - + :cvar int u_knots_to_add: the number of knots to add to the NURBS surfaces - along `u` direction before the deformation. + along `u` direction before the deformation. :cvar int v_knots_to_add: the number of knots to add to the NURBS surfaces along `v` direction before the deformation. :cvar int t_knots_to_add: the number of knots to add to the NURBS curves before the deformation. :cvar float tolerance: the tolerance involved in several internal operations of the procedure (joining wires in a single curve before - deformation and placing new poles on curves and surfaces). - + deformation and placing new poles on curves and surfaces). + :Example: >>> from pygem.cad import CustomDeformation @@ -64,15 +65,20 @@ class CustomDeformation(CADDeformation, OriginalCustomDeformation): >>> deform = CustomDeformation(move) >>> deform('original_shape.iges', dst='new_shape.iges') """ - def __init__(self, - func, - u_knots_to_add=0, - v_knots_to_add=0, - t_knots_to_add=0, - tolerance=1e-4): + + def __init__( + self, + func, + u_knots_to_add=0, + v_knots_to_add=0, + t_knots_to_add=0, + tolerance=1e-4, + ): OriginalCustomDeformation.__init__(self, func) - CADDeformation.__init__(self, - u_knots_to_add=u_knots_to_add, - v_knots_to_add=v_knots_to_add, - t_knots_to_add=t_knots_to_add, - tolerance=tolerance) + CADDeformation.__init__( + self, + u_knots_to_add=u_knots_to_add, + v_knots_to_add=v_knots_to_add, + t_knots_to_add=t_knots_to_add, + tolerance=tolerance, + ) diff --git a/pygem/cad/ffd.py b/pygem/cad/ffd.py index f3ffcdb3..bea7ea80 100644 --- a/pygem/cad/ffd.py +++ b/pygem/cad/ffd.py @@ -45,6 +45,7 @@ from pygem import FFD as OriginalFFD from .cad_deformation import CADDeformation + class FFD(CADDeformation, OriginalFFD): """ Class that handles the Free Form Deformation on CAD geometries. @@ -139,16 +140,20 @@ class FFD(CADDeformation, OriginalFFD): >>> # is equivalent to >>> new_shape = ffd(CADDeformation.read_shape(input_cad_file_name)) """ - def __init__(self, - n_control_points=None, - u_knots_to_add=0, - v_knots_to_add=0, - t_knots_to_add=0, - tolerance=1e-4): - OriginalFFD.__init__(self, - n_control_points=n_control_points) - CADDeformation.__init__(self, - u_knots_to_add=u_knots_to_add, - v_knots_to_add=v_knots_to_add, - t_knots_to_add=t_knots_to_add, - tolerance=tolerance) + + def __init__( + self, + n_control_points=None, + u_knots_to_add=0, + v_knots_to_add=0, + t_knots_to_add=0, + tolerance=1e-4, + ): + OriginalFFD.__init__(self, n_control_points=n_control_points) + CADDeformation.__init__( + self, + u_knots_to_add=u_knots_to_add, + v_knots_to_add=v_knots_to_add, + t_knots_to_add=t_knots_to_add, + tolerance=tolerance, + ) diff --git a/pygem/cad/idw.py b/pygem/cad/idw.py index 76e616de..903afbeb 100644 --- a/pygem/cad/idw.py +++ b/pygem/cad/idw.py @@ -41,6 +41,7 @@ from pygem import IDW as OriginalIDW from .cad_deformation import CADDeformation + class IDW(CADDeformation, OriginalIDW): """ Class that perform the Inverse Distance Weighting (IDW) on CAD geometries. @@ -84,7 +85,7 @@ class IDW(CADDeformation, OriginalIDW): default value only if the input file scale is significantly different form mm, making some of the aforementioned operations fail. Default is 0.0001. - + :cvar int power: the power parameter. The default value is 2. :cvar numpy.ndarray original_control_points: it is an (*n_control_points*, *3*) array with the coordinates of the original @@ -135,20 +136,27 @@ class IDW(CADDeformation, OriginalIDW): >>> modified_cad_file_name = "output.iges" >>> idw(input_cad_file_name, modified_cad_file_name) """ - def __init__(self, - original_control_points=None, - deformed_control_points=None, - power=2, - u_knots_to_add=0, - v_knots_to_add=0, - t_knots_to_add=0, - tolerance=1e-4): - OriginalIDW.__init__(self, - original_control_points=original_control_points, - deformed_control_points=deformed_control_points, - power=power) - CADDeformation.__init__(self, - u_knots_to_add=u_knots_to_add, - v_knots_to_add=v_knots_to_add, - t_knots_to_add=t_knots_to_add, - tolerance=tolerance) + + def __init__( + self, + original_control_points=None, + deformed_control_points=None, + power=2, + u_knots_to_add=0, + v_knots_to_add=0, + t_knots_to_add=0, + tolerance=1e-4, + ): + OriginalIDW.__init__( + self, + original_control_points=original_control_points, + deformed_control_points=deformed_control_points, + power=power, + ) + CADDeformation.__init__( + self, + u_knots_to_add=u_knots_to_add, + v_knots_to_add=v_knots_to_add, + t_knots_to_add=t_knots_to_add, + tolerance=tolerance, + ) diff --git a/pygem/cad/igeshandler.py b/pygem/cad/igeshandler.py index 6ad6ac48..67200140 100644 --- a/pygem/cad/igeshandler.py +++ b/pygem/cad/igeshandler.py @@ -1,8 +1,12 @@ """ Derived module from filehandler.py to handle iges and igs files. """ -from OCC.Core.IGESControl import (IGESControl_Reader, IGESControl_Writer, - IGESControl_Controller_Init) + +from OCC.Core.IGESControl import ( + IGESControl_Reader, + IGESControl_Writer, + IGESControl_Controller_Init, +) from OCC.Core.IFSelect import IFSelect_RetDone from pygem.cad import NurbsHandler @@ -31,7 +35,7 @@ class IgesHandler(NurbsHandler): def __init__(self): super(IgesHandler, self).__init__() - self.extensions = ['.iges', '.igs'] + self.extensions = [".iges", ".igs"] def load_shape_from_file(self, filename): """ diff --git a/pygem/cad/nurbshandler.py b/pygem/cad/nurbshandler.py index 22c5f126..dfb90280 100644 --- a/pygem/cad/nurbshandler.py +++ b/pygem/cad/nurbshandler.py @@ -4,29 +4,48 @@ File handling operations (reading/writing) must be implemented in derived classes. """ + import os import numpy as np from OCC.Core.BRep import BRep_Tool, BRep_Builder, BRep_Tool_Curve from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh from OCC.Core.BRepAlgo import brepalgo_IsValid from OCC.Core.BRepBuilderAPI import ( - BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeFace, - BRepBuilderAPI_NurbsConvert, BRepBuilderAPI_MakeWire, BRepBuilderAPI_Sewing) + BRepBuilderAPI_MakeEdge, + BRepBuilderAPI_MakeFace, + BRepBuilderAPI_NurbsConvert, + BRepBuilderAPI_MakeWire, + BRepBuilderAPI_Sewing, +) from OCC.Core.BRepOffsetAPI import BRepOffsetAPI_FindContigousEdges from OCC.Display.SimpleGui import init_display -from OCC.Core.GeomConvert import (geomconvert_SurfaceToBSplineSurface, - geomconvert_CurveToBSplineCurve) +from OCC.Core.GeomConvert import ( + geomconvert_SurfaceToBSplineSurface, + geomconvert_CurveToBSplineCurve, +) from OCC.Core.gp import gp_Pnt, gp_XYZ from OCC.Core.Precision import precision_Confusion from OCC.Core.ShapeAnalysis import ShapeAnalysis_WireOrder from OCC.Core.ShapeFix import ShapeFix_ShapeTolerance, ShapeFix_Shell from OCC.Core.StlAPI import StlAPI_Writer from OCC.Core.TColgp import TColgp_Array1OfPnt, TColgp_Array2OfPnt -from OCC.Core.TopAbs import (TopAbs_FACE, TopAbs_EDGE, TopAbs_WIRE, TopAbs_FORWARD, - TopAbs_SHELL) +from OCC.Core.TopAbs import ( + TopAbs_FACE, + TopAbs_EDGE, + TopAbs_WIRE, + TopAbs_FORWARD, + TopAbs_SHELL, +) from OCC.Core.TopExp import TopExp_Explorer, topexp -from OCC.Core.TopoDS import (topods_Face, TopoDS_Compound, topods_Shell, topods_Edge, - topods_Wire, topods, TopoDS_Shape) +from OCC.Core.TopoDS import ( + topods_Face, + TopoDS_Compound, + topods_Shell, + topods_Edge, + topods_Wire, + topods, + TopoDS_Shape, +) from matplotlib import pyplot from mpl_toolkits import mplot3d from stl import mesh @@ -69,7 +88,8 @@ def _check_infile_instantiation(self): """ if not self.shape or not self.infile: raise RuntimeError( - 'You can not write a file without having parsed one.') + "You can not write a file without having parsed one." + ) def load_shape_from_file(self, filename): """ @@ -77,9 +97,10 @@ def load_shape_from_file(self, filename): Not implemented, it has to be implemented in subclasses. """ - raise NotImplementedError('Subclass must implement abstract method' - '{}.load_shape_from_file'.format( - self.__class__.__name__)) + raise NotImplementedError( + "Subclass must implement abstract method" + "{}.load_shape_from_file".format(self.__class__.__name__) + ) def parse(self, filename): """ @@ -119,26 +140,30 @@ def parse(self, filename): n_poles_u = occ_face.NbUPoles() n_poles_v = occ_face.NbVPoles() control_polygon_coordinates = np.zeros( - shape=(n_poles_u * n_poles_v, 3)) + shape=(n_poles_u * n_poles_v, 3) + ) # cycle over the poles to get their coordinates i = 0 for pole_u_direction in range(n_poles_u): for pole_v_direction in range(n_poles_v): control_point_coordinates = occ_face.Pole( - pole_u_direction + 1, pole_v_direction + 1) + pole_u_direction + 1, pole_v_direction + 1 + ) control_polygon_coordinates[i, :] = [ control_point_coordinates.X(), control_point_coordinates.Y(), - control_point_coordinates.Z() + control_point_coordinates.Z(), ] i += 1 # pushing the control points coordinates to the mesh_points array # (used for FFD) mesh_points = np.append( - mesh_points, control_polygon_coordinates, axis=0) + mesh_points, control_polygon_coordinates, axis=0 + ) control_point_position.append( - control_point_position[-1] + n_poles_u * n_poles_v) + control_point_position[-1] + n_poles_u * n_poles_v + ) n_faces += 1 faces_explorer.Next() @@ -195,19 +220,20 @@ def write(self, mesh_points, filename, tolerance=None): for pole_u_direction in range(n_poles_u): for pole_v_direction in range(n_poles_v): control_point_coordinates = mesh_points[ - i + control_point_position[n_faces], :] + i + control_point_position[n_faces], : + ] point_xyz = gp_XYZ(*control_point_coordinates) gp_point = gp_Pnt(point_xyz) - occ_face.SetPole(pole_u_direction + 1, pole_v_direction + 1, - gp_point) + occ_face.SetPole( + pole_u_direction + 1, pole_v_direction + 1, gp_point + ) i += 1 # construct the deformed wire for the trimmed surfaces wire_maker = BRepBuilderAPI_MakeWire() tol = ShapeFix_ShapeTolerance() - brep = BRepBuilderAPI_MakeFace(occ_face, - self.tolerance).Face() + brep = BRepBuilderAPI_MakeFace(occ_face, self.tolerance).Face() brep_face = BRep_Tool.Surface(brep) # cycle on the edges @@ -219,7 +245,8 @@ def write(self, mesh_points, filename, tolerance=None): # evaluating the new edge: same (u,v) coordinates, but # different (x,y,x) ones edge_phis_coordinates_aux = BRepBuilderAPI_MakeEdge( - edge_uv_coordinates[0], brep_face) + edge_uv_coordinates[0], brep_face + ) edge_phis_coordinates = edge_phis_coordinates_aux.Edge() tol.SetTolerance(edge_phis_coordinates, self.tolerance) wire_maker.Add(edge_phis_coordinates) @@ -229,8 +256,7 @@ def write(self, mesh_points, filename, tolerance=None): wire = wire_maker.Wire() # trimming the surfaces - brep_surf = BRepBuilderAPI_MakeFace(occ_face, - wire).Shape() + brep_surf = BRepBuilderAPI_MakeFace(occ_face, wire).Shape() compound_builder.Add(compound, brep_surf) n_faces += 1 faces_explorer.Next() @@ -301,8 +327,9 @@ def parse_face(topo_face): bspline_converter = BRepBuilderAPI_NurbsConvert(edge) bspline_converter.Perform(edge) bspline_tshape_edge = bspline_converter.Shape() - h_geom_edge = BRep_Tool_Curve( - topods_Edge(bspline_tshape_edge))[0] + h_geom_edge = BRep_Tool_Curve(topods_Edge(bspline_tshape_edge))[ + 0 + ] h_bspline_edge = geomconvert_CurveToBSplineCurve(h_geom_edge) bspline_geom_edge = h_bspline_edge @@ -316,11 +343,11 @@ def parse_face(topo_face): for i in range(1, nb_poles + 1): ctrlpt = edge_ctrlpts.Value(i) ctrlpt_position = np.array( - [[ctrlpt.Coord(1), - ctrlpt.Coord(2), - ctrlpt.Coord(3)]]) + [[ctrlpt.Coord(1), ctrlpt.Coord(2), ctrlpt.Coord(3)]] + ) points_single_edge = np.append( - points_single_edge, ctrlpt_position, axis=0) + points_single_edge, ctrlpt_position, axis=0 + ) mesh_points_edge.append(points_single_edge) @@ -347,11 +374,11 @@ def parse_face(topo_face): for indice_v_direction in range(1, nb_v + 1): ctrlpt = ctrlpts.Value(indice_u_direction, indice_v_direction) ctrlpt_position = np.array( - [[ctrlpt.Coord(1), - ctrlpt.Coord(2), - ctrlpt.Coord(3)]]) + [[ctrlpt.Coord(1), ctrlpt.Coord(2), ctrlpt.Coord(3)]] + ) mesh_points_face = np.append( - mesh_points_face, ctrlpt_position, axis=0) + mesh_points_face, ctrlpt_position, axis=0 + ) return mesh_points_face, mesh_points_edge @@ -442,8 +469,10 @@ def write_edge(points_edge, topo_edge): nb_cpt = bspline_edge_curve.NbPoles() # check consistency if points_edge.shape[0] != nb_cpt: - raise ValueError("Input control points do not have not have the " - "same number as the geometric edge!") + raise ValueError( + "Input control points do not have not have the " + "same number as the geometric edge!" + ) else: for i in range(1, nb_cpt + 1): @@ -482,8 +511,10 @@ def write_face(self, points_face, list_points_edge, topo_face, toledge): nb_v = bsurface.NbVPoles() # check consistency if points_face.shape[0] != nb_u * nb_v: - raise ValueError("Input control points do not have not have the " - "same number as the geometric face!") + raise ValueError( + "Input control points do not have not have the " + "same number as the geometric face!" + ) # cycle on the face points indice_cpt = 0 @@ -500,7 +531,8 @@ def write_face(self, points_face, list_points_edge, topo_face, toledge): # cycle on the wires face_wires_explorer = TopExp_Explorer( - topo_nurbsface.Oriented(TopAbs_FORWARD), TopAbs_WIRE) + topo_nurbsface.Oriented(TopAbs_FORWARD), TopAbs_WIRE + ) ind_edge_total = 0 while face_wires_explorer.More(): @@ -510,7 +542,8 @@ def write_face(self, points_face, list_points_edge, topo_face, toledge): # cycle on the edges ind_edge = 0 wire_explorer_edge = TopExp_Explorer( - twire.Oriented(TopAbs_FORWARD), TopAbs_EDGE) + twire.Oriented(TopAbs_FORWARD), TopAbs_EDGE + ) # check edges order on the wire mode3d = True tolerance_edges = toledge @@ -522,7 +555,8 @@ def write_face(self, points_face, list_points_edge, topo_face, toledge): while wire_explorer_edge.More(): tedge = topods_Edge(wire_explorer_edge.Current()) new_bspline_tedge = self.write_edge( - list_points_edge[ind_edge_total], tedge) + list_points_edge[ind_edge_total], tedge + ) deformed_edges.append(new_bspline_tedge) analyzer = topexp() @@ -555,7 +589,8 @@ def write_face(self, points_face, list_points_edge, topo_face, toledge): new_bspline_twire.Add(new_edge_toadd) else: deformed_edge_revers = deformed_edges[ - np.abs(deformed_edge_i) - 1] + np.abs(deformed_edge_i) - 1 + ] stol.SetTolerance(deformed_edge_revers, toledge) new_bspline_twire.Add(deformed_edge_revers) if new_bspline_twire.Error() != 0: @@ -634,7 +669,8 @@ def write_shape(self, l_shells, filename, tol): if self.check_topo == 0: # cycle on shells (multiple objects) shape_shells_explorer = TopExp_Explorer( - self.shape.Oriented(TopAbs_FORWARD), TopAbs_SHELL) + self.shape.Oriented(TopAbs_FORWARD), TopAbs_SHELL + ) ishell = 0 while shape_shells_explorer.More(): @@ -646,13 +682,17 @@ def write_shape(self, l_shells, filename, tol): # cycle on faces faces_explorer = TopExp_Explorer( - per_shell.Oriented(TopAbs_FORWARD), TopAbs_FACE) + per_shell.Oriented(TopAbs_FORWARD), TopAbs_FACE + ) iface = 0 while faces_explorer.More(): topoface = topods.Face(faces_explorer.Current()) - newface = self.write_face(l_shells[ishell][iface][0], - l_shells[ishell][iface][1], - topoface, tol) + newface = self.write_face( + l_shells[ishell][iface][0], + l_shells[ishell][iface][1], + topoface, + tol, + ) # add face to compound compound_builder.Add(comp, newface) @@ -665,8 +705,8 @@ def write_shape(self, l_shells, filename, tol): global_compound_builder.Add(global_comp, new_shell) # TODO - #print("Shell {0} of type {1} Processed ".format(ishell, itype)) - #print "==============================================" + # print("Shell {0} of type {1} Processed ".format(ishell, itype)) + # print "==============================================" ishell += 1 shape_shells_explorer.Next() @@ -680,12 +720,14 @@ def write_shape(self, l_shells, filename, tol): # cycle on faces faces_explorer = TopExp_Explorer( - self.shape.Oriented(TopAbs_FORWARD), TopAbs_FACE) + self.shape.Oriented(TopAbs_FORWARD), TopAbs_FACE + ) iface = 0 while faces_explorer.More(): topoface = topods.Face(faces_explorer.Current()) - newface = self.write_face(l_shells[0][iface][0], - l_shells[0][iface][1], topoface, tol) + newface = self.write_face( + l_shells[0][iface][0], l_shells[0][iface][1], topoface, tol + ) # add face to compound compound_builder.Add(comp, newface) @@ -709,9 +751,11 @@ def write_shape_to_file(self, shape, filename): Not implemented, it has to be implemented in subclasses. """ - raise NotImplementedError(\ - "Subclass must implement abstract method " +\ - self.__class__.__name__ + ".write_shape_to_file") + raise NotImplementedError( + "Subclass must implement abstract method " + + self.__class__.__name__ + + ".write_shape_to_file" + ) def plot(self, plot_file=None, save_fig=False): """ @@ -740,44 +784,54 @@ def plot(self, plot_file=None, save_fig=False): stl_mesh = BRepMesh_IncrementalMesh(shape, 0.01) stl_mesh.Perform() - f = stl_writer.Write(shape, 'aux_figure.stl') + f = stl_writer.Write(shape, "aux_figure.stl") # Create a new plot figure = pyplot.figure() axes = mplot3d.Axes3D(figure) # Load the STL files and add the vectors to the plot - stl_mesh = mesh.Mesh.from_file('aux_figure.stl') - os.remove('aux_figure.stl') + stl_mesh = mesh.Mesh.from_file("aux_figure.stl") + os.remove("aux_figure.stl") axes.add_collection3d( - mplot3d.art3d.Poly3DCollection(stl_mesh.vectors / 1000)) + mplot3d.art3d.Poly3DCollection(stl_mesh.vectors / 1000) + ) # Get the limits of the axis and center the geometry - max_dim = np.array([\ - np.max(stl_mesh.vectors[:, :, 0]) / 1000,\ - np.max(stl_mesh.vectors[:, :, 1]) / 1000,\ - np.max(stl_mesh.vectors[:, :, 2]) / 1000]) - min_dim = np.array([\ - np.min(stl_mesh.vectors[:, :, 0]) / 1000,\ - np.min(stl_mesh.vectors[:, :, 1]) / 1000,\ - np.min(stl_mesh.vectors[:, :, 2]) / 1000]) + max_dim = np.array( + [ + np.max(stl_mesh.vectors[:, :, 0]) / 1000, + np.max(stl_mesh.vectors[:, :, 1]) / 1000, + np.max(stl_mesh.vectors[:, :, 2]) / 1000, + ] + ) + min_dim = np.array( + [ + np.min(stl_mesh.vectors[:, :, 0]) / 1000, + np.min(stl_mesh.vectors[:, :, 1]) / 1000, + np.min(stl_mesh.vectors[:, :, 2]) / 1000, + ] + ) max_lenght = np.max(max_dim - min_dim) - axes.set_xlim(\ - -.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2,\ - .6 * max_lenght + (max_dim[0] + min_dim[0]) / 2) - axes.set_ylim(\ - -.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2,\ - .6 * max_lenght + (max_dim[1] + min_dim[1]) / 2) - axes.set_zlim(\ - -.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2,\ - .6 * max_lenght + (max_dim[2] + min_dim[2]) / 2) + axes.set_xlim( + -0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + 0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + ) + axes.set_ylim( + -0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + 0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + ) + axes.set_zlim( + -0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + 0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + ) # Show the plot to the screen if not save_fig: pyplot.show() else: - figure.savefig(plot_file.split('.')[0] + '.png') + figure.savefig(plot_file.split(".")[0] + ".png") return figure diff --git a/pygem/cad/rbf.py b/pygem/cad/rbf.py index 3f5c310f..9709332a 100644 --- a/pygem/cad/rbf.py +++ b/pygem/cad/rbf.py @@ -24,7 +24,7 @@ is defines as follows .. math:: - \\mathcal{M}(\\boldsymbol{x}) = p(\\boldsymbol{x}) + + \\mathcal{M}(\\boldsymbol{x}) = p(\\boldsymbol{x}) + \\sum_{i=1}^{\\mathcal{N}_C} \\gamma_i \\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}} \\|) @@ -60,11 +60,12 @@ from pygem import RBF as OriginalRBF from .cad_deformation import CADDeformation + class RBF(CADDeformation, OriginalRBF): """ Class that handles the Radial Basis Functions interpolation on CAD geometries. - + :param numpy.ndarray original_control_points: it is an (*n_control_points*, *3*) array with the coordinates of the original interpolation control points before the deformation. The default is the @@ -83,7 +84,7 @@ class RBF(CADDeformation, OriginalRBF): basis functions. For details see the class :class:`~pygem.radialbasis.RBF`. The default value is 0.5. :param dict extra_parameter: the additional parameters that may be passed to - the kernel function. Default is None. + the kernel function. Default is None. :param int u_knots_to_add: the number of knots to add to the NURBS surfaces along `u` direction before the deformation. This parameter is useful whenever the gradient of the imposed deformation present spatial scales @@ -114,7 +115,7 @@ class RBF(CADDeformation, OriginalRBF): default value only if the input file scale is significantly different form mm, making some of the aforementioned operations fail. Default is 0.0001. - + :cvar numpy.ndarray weights: the matrix formed by the weights corresponding to the a-priori selected N control points, associated to the basis functions and c and Q terms that describe the polynomial of order one @@ -131,9 +132,9 @@ class RBF(CADDeformation, OriginalRBF): :cvar float radius: the scaling parameter that affects the shape of the basis functions. :cvar dict extra_parameter: the additional parameters that may be passed to - the kernel function. + the kernel function. :cvar int u_knots_to_add: the number of knots to add to the NURBS surfaces - along `u` direction before the deformation. + along `u` direction before the deformation. :cvar int v_knots_to_add: the number of knots to add to the NURBS surfaces along `v` direction before the deformation. :cvar int t_knots_to_add: the number of knots to add to the NURBS curves @@ -152,24 +153,31 @@ class RBF(CADDeformation, OriginalRBF): >>> modified_cad_file_name = "output.iges" >>> rbf(input_cad_file_name, modified_cad_file_name) """ - def __init__(self, - original_control_points=None, - deformed_control_points=None, - func='gaussian_spline', - radius=0.5, - extra_parameter=None, - u_knots_to_add=0, - v_knots_to_add=0, - t_knots_to_add=0, - tolerance=1e-4): - OriginalRBF.__init__(self, - original_control_points=original_control_points, - deformed_control_points=deformed_control_points, - func=func, - radius=radius, - extra_parameter=extra_parameter) - CADDeformation.__init__(self, - u_knots_to_add=u_knots_to_add, - v_knots_to_add=v_knots_to_add, - t_knots_to_add=t_knots_to_add, - tolerance=tolerance) + + def __init__( + self, + original_control_points=None, + deformed_control_points=None, + func="gaussian_spline", + radius=0.5, + extra_parameter=None, + u_knots_to_add=0, + v_knots_to_add=0, + t_knots_to_add=0, + tolerance=1e-4, + ): + OriginalRBF.__init__( + self, + original_control_points=original_control_points, + deformed_control_points=deformed_control_points, + func=func, + radius=radius, + extra_parameter=extra_parameter, + ) + CADDeformation.__init__( + self, + u_knots_to_add=u_knots_to_add, + v_knots_to_add=v_knots_to_add, + t_knots_to_add=t_knots_to_add, + tolerance=tolerance, + ) diff --git a/pygem/cad/stephandler.py b/pygem/cad/stephandler.py index 7a43d7af..239d7a58 100644 --- a/pygem/cad/stephandler.py +++ b/pygem/cad/stephandler.py @@ -1,6 +1,7 @@ """ Derived module from nurbshandler.py to handle step and stp files. """ + from OCC.Core.IFSelect import IFSelect_RetDone from OCC.Core.Interface import Interface_Static_SetCVal from OCC.Core.STEPControl import STEPControl_Writer, STEPControl_Reader @@ -34,7 +35,7 @@ class StepHandler(NurbsHandler): def __init__(self): super(StepHandler, self).__init__() self._control_point_position = None - self.extensions = ['.step', '.stp'] + self.extensions = [".step", ".stp"] def load_shape_from_file(self, filename): """ diff --git a/pygem/cffd.py b/pygem/cffd.py index 80337fec..1d864b6f 100644 --- a/pygem/cffd.py +++ b/pygem/cffd.py @@ -2,8 +2,8 @@ Utilities for performing Constrained Free Form Deformation (CFFD). :Theoretical Insight: - - It performs Free Form Deformation while trying to enforce a costraint of the form F(x)=c. + + It performs Free Form Deformation while trying to enforce a costraint of the form F(x)=c. The constraint is enforced exactly (up to numerical errors) if and only if the function is linear. For details on Free Form Deformation see the mother class. @@ -20,7 +20,7 @@ class CFFD(FFD): :param list n_control_points: number of control points in the x, y, and z direction. Default is [2, 2, 2]. - :param string mode: it can be ``affine`` or ``triaffine``. The first option is for the F that are affine in all the coordinates of the points. + :param string mode: it can be ``affine`` or ``triaffine``. The first option is for the F that are affine in all the coordinates of the points. The second one is for functions that are F in the coordinates of the points. The first option implies the second, but is optimal for that class of functions. :cvar numpy.ndarray box_length: dimension of the FFD bounding box, in the x, y and z direction (local coordinate system). @@ -36,11 +36,11 @@ class CFFD(FFD): z, normalized with the box length z. :cvar callable fun: it defines the F of the constraint F(x)=c. Default is the constant 1 function. :cvar numpy.ndarray fixval: it defines the c of the constraint F(x)=c. Default is 1. - :cvar numpy.ndarray ffd_mask: a boolean tensor that tells to the class - which control points can be moved, and in what direction, to enforce the constraint. + :cvar numpy.ndarray ffd_mask: a boolean tensor that tells to the class + which control points can be moved, and in what direction, to enforce the constraint. The tensor has shape (n_x,n_y,n_z,3), where the last dimension indicates movement on x,y,z respectively. Default is all true. - :cvar numpy.ndarray fun_mask: a boolean tensor that tells to the class + :cvar numpy.ndarray fun_mask: a boolean tensor that tells to the class on which axis which constraint depends on. The tensor has shape (n_cons,3), where the last dimension indicates dependency on on x,y,z respectively. Default is all true. It used only in the triaffine mode. @@ -58,18 +58,16 @@ class CFFD(FFD): >>> new_mesh_points = cffd.ffd(original_mesh_points) """ - def __init__(self, - fixval, - fun, - n_control_points=None, - ffd_mask=None, - fun_mask=None): + + def __init__( + self, fixval, fun, n_control_points=None, ffd_mask=None, fun_mask=None + ): super().__init__(n_control_points) if ffd_mask is None: - self.ffd_mask = np.full((*self.n_control_points, 3), - True, - dtype=bool) + self.ffd_mask = np.full( + (*self.n_control_points, 3), True, dtype=bool + ) else: self.ffd_mask = ffd_mask @@ -82,60 +80,66 @@ def __init__(self, self.fun_mask = fun_mask def adjust_control_points(self, src_pts): - ''' + """ Adjust the FFD control points such that fun(ffd(src_pts))=fixval - - :param np.ndarray src_pts: the points whose deformation we want to be + + :param np.ndarray src_pts: the points whose deformation we want to be constrained. :rtype: None. - ''' + """ hyper_param = self.fun_mask.copy().astype(float) hyper_param = hyper_param / np.sum(hyper_param, axis=1) mask_bak = self.ffd_mask.copy() fixval_bak = self.fixval.copy() diffvolume = self.fixval - self.fun(self.ffd(src_pts)) for i in range(3): - self.ffd_mask = np.full((*self.n_control_points, 3), - False, - dtype=bool) + self.ffd_mask = np.full( + (*self.n_control_points, 3), False, dtype=bool + ) self.ffd_mask[:, :, :, i] = mask_bak[:, :, :, i].copy() - self.fixval = self.fun( - self.ffd(src_pts)) + hyper_param[:, i] * (diffvolume) + self.fixval = self.fun(self.ffd(src_pts)) + hyper_param[:, i] * ( + diffvolume + ) saved_parameters = self._save_parameters() - indices = np.arange(np.prod(self.n_control_points) * - 3)[self.ffd_mask.reshape(-1)] - A, b = self._compute_linear_map(src_pts, saved_parameters.copy(), - indices) + indices = np.arange(np.prod(self.n_control_points) * 3)[ + self.ffd_mask.reshape(-1) + ] + A, b = self._compute_linear_map( + src_pts, saved_parameters.copy(), indices + ) A = A[self.fun_mask[:, i].reshape(-1), :] b = b[self.fun_mask[:, i].reshape(-1)] d = A @ saved_parameters[indices] + b fixval = self.fixval[self.fun_mask[:, i].reshape(-1)] - deltax = np.linalg.multi_dot([ - A.T, - np.linalg.inv(np.linalg.multi_dot([A, A.T])), (fixval - d) - ]) + deltax = np.linalg.multi_dot( + [ + A.T, + np.linalg.inv(np.linalg.multi_dot([A, A.T])), + (fixval - d), + ] + ) saved_parameters[indices] = saved_parameters[indices] + deltax self._load_parameters(saved_parameters) self.ffd_mask = mask_bak.copy() self.fixval = fixval_bak.copy() def ffd(self, src_pts): - ''' + """ Performs Classic Free Form Deformation. - + :param np.ndarray src_pts: the points to deform. :return: the deformed points. :rtype: numpy.ndarray - ''' + """ return super().__call__(src_pts) def _save_parameters(self): - ''' + """ Saves the FFD control points in an array of shape [n_x,ny,nz,3]. :return: the FFD control points in an array of shape [n_x,ny,nz,3]. :rtype: numpy.ndarray - ''' + """ tmp = np.zeros([*self.n_control_points, 3]) tmp[:, :, :, 0] = self.array_mu_x tmp[:, :, :, 1] = self.array_mu_y @@ -143,45 +147,49 @@ def _save_parameters(self): return tmp.reshape(-1) def _load_parameters(self, tmp): - ''' + """ Loads the FFD control points from an array of shape [n_x,ny,nz,3]. :param np.ndarray tmp: the array of FFD control points. :rtype: None - ''' + """ tmp = tmp.reshape(*self.n_control_points, 3) self.array_mu_x = tmp[:, :, :, 0] self.array_mu_y = tmp[:, :, :, 1] self.array_mu_z = tmp[:, :, :, 2] def _compute_linear_map(self, src_pts, saved_parameters, indices): - ''' + """ Computes the coefficient and the intercept of the linear map from the control points to the output. - + :param np.ndarray src_pts: the points to deform. :param np.ndarray saved_parameters: the array of FFD control points. :return: a tuple containing the coefficient and the intercept. :rtype: tuple(np.ndarray,np.ndarray) - ''' + """ n_indices = len(indices) inputs = np.zeros([n_indices + 1, n_indices + 1]) outputs = np.zeros([n_indices + 1, self.fixval.shape[0]]) np.random.seed(0) - for i in range(n_indices + - 1): ##now we generate the interpolation points + for i in range( + n_indices + 1 + ): ##now we generate the interpolation points tmp = np.random.rand(1, n_indices) tmp = tmp.reshape(1, -1) - inputs[i] = np.hstack([tmp, np.ones( - (tmp.shape[0], 1))]) #dependent variable + inputs[i] = np.hstack( + [tmp, np.ones((tmp.shape[0], 1))] + ) # dependent variable saved_parameters[indices] = tmp self._load_parameters( saved_parameters - ) #loading the depent variable as a control point + ) # loading the depent variable as a control point def_pts = super().__call__( - src_pts) #computing the deformation with the dependent variable - outputs[i] = self.fun(def_pts) #computing the independent variable - sol = np.linalg.lstsq(inputs, outputs, - rcond=None) #computation of the linear map - A = sol[0].T[:, :-1] #coefficient - b = sol[0].T[:, -1] #intercept - return A, b \ No newline at end of file + src_pts + ) # computing the deformation with the dependent variable + outputs[i] = self.fun(def_pts) # computing the independent variable + sol = np.linalg.lstsq( + inputs, outputs, rcond=None + ) # computation of the linear map + A = sol[0].T[:, :-1] # coefficient + b = sol[0].T[:, -1] # intercept + return A, b diff --git a/pygem/custom_deformation.py b/pygem/custom_deformation.py index d938099a..7a99e5b8 100644 --- a/pygem/custom_deformation.py +++ b/pygem/custom_deformation.py @@ -1,10 +1,12 @@ """ Module for a custom deformation. """ + import numpy as np from pygem import Deformation + class CustomDeformation(Deformation): """ Class to perform a custom deformation to the mesh points. diff --git a/pygem/deformation.py b/pygem/deformation.py index 7cfb13ce..a3dca9b2 100644 --- a/pygem/deformation.py +++ b/pygem/deformation.py @@ -3,14 +3,15 @@ """ from abc import ABC, abstractmethod - + + class Deformation(ABC): """ Abstract class for generic deformation. This class should be inherited for the development of new deformation techniques. """ - + @abstractmethod def __init__(self, value): pass diff --git a/pygem/elmerhandler.py b/pygem/elmerhandler.py index 11569176..5dc6d0ff 100644 --- a/pygem/elmerhandler.py +++ b/pygem/elmerhandler.py @@ -1,6 +1,7 @@ """ Derived module from filehandler.py to handle ElmerFEM files. """ + import numpy as np import pygem.filehandler as fh @@ -17,7 +18,7 @@ class ElmerHandler(fh.FileHandler): def __init__(self): super(ElmerHandler, self).__init__() - self.extensions = ['.nodes'] + self.extensions = [".nodes"] def parse(self, filename): """ @@ -25,7 +26,7 @@ def parse(self, filename): the coordinates. :param string filename: name of the input file. - + :return: mesh_points: it is a `n_points`-by-3 matrix containing the coordinates of the points of the mesh :rtype: numpy.ndarray @@ -41,17 +42,17 @@ def parse(self, filename): i = 0 n_points = 0 - with open(self.infile, 'r') as input_file: + with open(self.infile, "r") as input_file: for line in input_file: n_points += 1 mesh_points = np.zeros(shape=(n_points, 3)) - with open(self.infile, 'r') as input_file: - + with open(self.infile, "r") as input_file: + i = 0 for line in input_file: - numbers = line.split() #[n1 p x y z] -> [x y z] - del numbers[0:2] + numbers = line.split() # [n1 p x y z] -> [x y z] + del numbers[0:2] j = 0 for number in numbers: @@ -82,16 +83,25 @@ def write(self, mesh_points, filename): n_points = mesh_points.shape[0] i = 0 - with open(self.infile, 'r') as input_file, open(self.outfile, - 'w') as output_file: + with ( + open(self.infile, "r") as input_file, + open(self.outfile, "w") as output_file, + ): for line in input_file: - numbers = line.split() #[n1 p x y z] - - output_file.write(numbers[0] + ' ' +numbers[1] + ' ' \ - + str(mesh_points[i][0]) + ' ' + str(mesh_points[i][1]) \ - + ' ' + str(mesh_points[i][2])) + numbers = line.split() # [n1 p x y z] + + output_file.write( + numbers[0] + + " " + + numbers[1] + + " " + + str(mesh_points[i][0]) + + " " + + str(mesh_points[i][1]) + + " " + + str(mesh_points[i][2]) + ) i += 1 if i != n_points: - output_file.write('\n') - \ No newline at end of file + output_file.write("\n") diff --git a/pygem/ffd.py b/pygem/ffd.py index 9b40ead1..a6bd9522 100644 --- a/pygem/ffd.py +++ b/pygem/ffd.py @@ -8,7 +8,7 @@ *Sederberg, Thomas W., and Scott R. Parry. "Free-form deformation of solid geometric models." ACM SIGGRAPH computer graphics 20.4 (1986): 151-160*. It consists in three different step: - + - Mapping the physical domain to the reference one with map :math:`\\boldsymbol{\\psi}`. In the code it is named *transformation*. @@ -35,11 +35,12 @@ where :math:`\\mathsf{b}_{lmn}` are Bernstein polynomials. We improve the traditional version by allowing a rotation of the FFD lattice in order to give more flexibility to the tool. - + You can try to add more shapes to the lattice to allow more and more involved transformations. """ + try: import configparser as configparser except ImportError: @@ -59,7 +60,7 @@ class FFD(Deformation): :param list n_control_points: number of control points in the x, y, and z direction. Default is [2, 2, 2]. - + :cvar numpy.ndarray box_length: dimension of the FFD bounding box, in the x, y and z direction (local coordinate system). :cvar numpy.ndarray box_origin: the x, y and z coordinates of the origin of @@ -86,14 +87,15 @@ class FFD(Deformation): >>> 'tests/test_datasets/meshpoints_sphere_orig.npy') >>> new_mesh_points = ffd(original_mesh_points) """ + reference_frame = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]]) def __init__(self, n_control_points=None): - self.conversion_unit = 1. # TODO: unused at the moment + self.conversion_unit = 1.0 # TODO: unused at the moment - self.box_length = np.array([1., 1., 1.]) - self.box_origin = np.array([0., 0., 0.]) - self.rot_angle = np.array([0., 0., 0.]) + self.box_length = np.array([1.0, 1.0, 1.0]) + self.box_origin = np.array([0.0, 0.0, 0.0]) + self.rot_angle = np.array([0.0, 0.0, 0.0]) self.array_mu_x = None self.array_mu_y = None @@ -148,9 +150,10 @@ def T(self): :rtype: callable """ + def T_mapping(points): - (n_rows, n_cols) = points.shape - (dim_n_mu, dim_m_mu, dim_t_mu) = self.array_mu_x.shape + n_rows, n_cols = points.shape + dim_n_mu, dim_m_mu, dim_t_mu = self.array_mu_x.shape # Initialization. In order to exploit the contiguity in memory the # following are transposed @@ -160,37 +163,41 @@ def T_mapping(points): shift_points = np.zeros((n_cols, n_rows)) # TODO check no-loop implementation - #bernstein_x = ( + # bernstein_x = ( # np.power(mesh_points[:, 0][:, None], range(dim_n_mu)) * # np.power(1 - mesh_points[:, 0][:, None], range(dim_n_mu-1, -1, -1)) * # special.binom(np.array([dim_n_mu-1]*dim_n_mu), np.arange(dim_n_mu)) - #) + # ) for i in range(0, dim_n_mu): aux1 = np.power((1 - points[:, 0]), dim_n_mu - 1 - i) aux2 = np.power(points[:, 0], i) - bernstein_x[i, :] = (special.binom(dim_n_mu - 1, i) * - np.multiply(aux1, aux2)) + bernstein_x[i, :] = special.binom( + dim_n_mu - 1, i + ) * np.multiply(aux1, aux2) for i in range(0, dim_m_mu): aux1 = np.power((1 - points[:, 1]), dim_m_mu - 1 - i) aux2 = np.power(points[:, 1], i) - bernstein_y[i, :] = special.binom(dim_m_mu - 1, - i) * np.multiply(aux1, aux2) + bernstein_y[i, :] = special.binom( + dim_m_mu - 1, i + ) * np.multiply(aux1, aux2) for i in range(0, dim_t_mu): aux1 = np.power((1 - points[:, 2]), dim_t_mu - 1 - i) aux2 = np.power(points[:, 2], i) - bernstein_z[i, :] = special.binom(dim_t_mu - 1, - i) * np.multiply(aux1, aux2) + bernstein_z[i, :] = special.binom( + dim_t_mu - 1, i + ) * np.multiply(aux1, aux2) - aux_x = 0. - aux_y = 0. - aux_z = 0. + aux_x = 0.0 + aux_y = 0.0 + aux_z = 0.0 for j in range(0, dim_m_mu): for k in range(0, dim_t_mu): - bernstein_yz = np.multiply(bernstein_y[j, :], - bernstein_z[k, :]) + bernstein_yz = np.multiply( + bernstein_y[j, :], bernstein_z[k, :] + ) for i in range(0, dim_n_mu): aux = np.multiply(bernstein_x[i, :], bernstein_yz) aux_x += aux * self.array_mu_x[i, j, k] @@ -212,9 +219,11 @@ def rotation_matrix(self): :rtype: numpy.ndarray """ - return angles2matrix(np.radians(self.rot_angle[2]), - np.radians(self.rot_angle[1]), - np.radians(self.rot_angle[0])) + return angles2matrix( + np.radians(self.rot_angle[2]), + np.radians(self.rot_angle[1]), + np.radians(self.rot_angle[0]), + ) @property def position_vertices(self): @@ -223,10 +232,12 @@ def position_vertices(self): :rtype: numpy.ndarray """ - return self.box_origin + np.vstack([ - np.zeros((1, 3)), - self.rotation_matrix.dot(np.diag(self.box_length)).T - ]) + return self.box_origin + np.vstack( + [ + np.zeros((1, 3)), + self.rotation_matrix.dot(np.diag(self.box_length)).T, + ] + ) def reset_weights(self): """ @@ -235,7 +246,7 @@ def reset_weights(self): self.array_mu_x.fill(0.0) self.array_mu_y.fill(0.0) self.array_mu_z.fill(0.0) - + def fit_to_points(self, points, padding_factor=0.05): """ Automatically set the FFD box_origin and box_length based on the input points. @@ -255,8 +266,7 @@ def fit_to_points(self, points, padding_factor=0.05): self.box_length = side_lengths + 2 * padding self.box_origin = min_coords - padding - - def read_parameters(self, filename='parameters.prm'): + def read_parameters(self, filename="parameters.prm"): """ Reads in the parameters file and fill the self structure. @@ -274,46 +284,49 @@ def read_parameters(self, filename='parameters.prm'): config = configparser.RawConfigParser() config.read(filename) - self.n_control_points[0] = config.getint('Box info', - 'n control points x') - self.n_control_points[1] = config.getint('Box info', - 'n control points y') - self.n_control_points[2] = config.getint('Box info', - 'n control points z') + self.n_control_points[0] = config.getint( + "Box info", "n control points x" + ) + self.n_control_points[1] = config.getint( + "Box info", "n control points y" + ) + self.n_control_points[2] = config.getint( + "Box info", "n control points z" + ) - self.box_length[0] = config.getfloat('Box info', 'box length x') - self.box_length[1] = config.getfloat('Box info', 'box length y') - self.box_length[2] = config.getfloat('Box info', 'box length z') + self.box_length[0] = config.getfloat("Box info", "box length x") + self.box_length[1] = config.getfloat("Box info", "box length y") + self.box_length[2] = config.getfloat("Box info", "box length z") - self.box_origin[0] = config.getfloat('Box info', 'box origin x') - self.box_origin[1] = config.getfloat('Box info', 'box origin y') - self.box_origin[2] = config.getfloat('Box info', 'box origin z') + self.box_origin[0] = config.getfloat("Box info", "box origin x") + self.box_origin[1] = config.getfloat("Box info", "box origin y") + self.box_origin[2] = config.getfloat("Box info", "box origin z") - self.rot_angle[0] = config.getfloat('Box info', 'rotation angle x') - self.rot_angle[1] = config.getfloat('Box info', 'rotation angle y') - self.rot_angle[2] = config.getfloat('Box info', 'rotation angle z') + self.rot_angle[0] = config.getfloat("Box info", "rotation angle x") + self.rot_angle[1] = config.getfloat("Box info", "rotation angle y") + self.rot_angle[2] = config.getfloat("Box info", "rotation angle z") self.array_mu_x = np.zeros(self.n_control_points) self.array_mu_y = np.zeros(self.n_control_points) self.array_mu_z = np.zeros(self.n_control_points) - mux = config.get('Parameters weights', 'parameter x') - muy = config.get('Parameters weights', 'parameter y') - muz = config.get('Parameters weights', 'parameter z') + mux = config.get("Parameters weights", "parameter x") + muy = config.get("Parameters weights", "parameter y") + muz = config.get("Parameters weights", "parameter z") - for line in mux.split('\n'): + for line in mux.split("\n"): values = np.array(line.split()) self.array_mu_x[tuple(map(int, values[0:3]))] = float(values[3]) - for line in muy.split('\n'): + for line in muy.split("\n"): values = line.split() self.array_mu_y[tuple(map(int, values[0:3]))] = float(values[3]) - for line in muz.split('\n'): + for line in muz.split("\n"): values = line.split() self.array_mu_z[tuple(map(int, values[0:3]))] = float(values[3]) - def write_parameters(self, filename='parameters.prm'): + def write_parameters(self, filename="parameters.prm"): """ This method writes a parameters file (.prm) called `filename` and fills it with all the parameters class members. @@ -324,124 +337,151 @@ def write_parameters(self, filename='parameters.prm'): raise TypeError("filename must be a string") output_string = "" - output_string += '\n[Box info]\n' - output_string += '# This section collects all the properties of the' - output_string += ' FFD bounding box.\n' - - output_string += '\n# n control points indicates the number of control' - output_string += ' points in each direction (x, y, z).\n' - output_string += '# For example, to create a 2 x 3 x 2 grid, use the' - output_string += ' following: n control points: 2, 3, 2\n' - output_string += 'n control points x: ' + str( - self.n_control_points[0]) + '\n' - output_string += 'n control points y: ' + str( - self.n_control_points[1]) + '\n' - output_string += 'n control points z: ' + str( - self.n_control_points[2]) + '\n' - - output_string += '\n# box length indicates the length of the FFD ' - output_string += 'bounding box along the three canonical directions ' - output_string += '(x, y, z).\n' - - output_string += '# It uses the local coordinate system.\n' - output_string += '# For example to create a 2 x 1.5 x 3 meters box ' - output_string += 'use the following: box length: 2.0, 1.5, 3.0\n' - - output_string += 'box length x: ' + str(self.box_length[0]) + '\n' - output_string += 'box length y: ' + str(self.box_length[1]) + '\n' - output_string += 'box length z: ' + str(self.box_length[2]) + '\n' - - output_string += '\n# box origin indicates the x, y, and z coordinates ' - output_string += 'of the origin of the FFD bounding box. That is ' - output_string += 'center of\n' - - output_string += '# rotation of the bounding box. It corresponds to ' - output_string += 'the point coordinates with position [0][0][0].\n' + output_string += "\n[Box info]\n" + output_string += "# This section collects all the properties of the" + output_string += " FFD bounding box.\n" + + output_string += "\n# n control points indicates the number of control" + output_string += " points in each direction (x, y, z).\n" + output_string += "# For example, to create a 2 x 3 x 2 grid, use the" + output_string += " following: n control points: 2, 3, 2\n" + output_string += ( + "n control points x: " + str(self.n_control_points[0]) + "\n" + ) + output_string += ( + "n control points y: " + str(self.n_control_points[1]) + "\n" + ) + output_string += ( + "n control points z: " + str(self.n_control_points[2]) + "\n" + ) + + output_string += "\n# box length indicates the length of the FFD " + output_string += "bounding box along the three canonical directions " + output_string += "(x, y, z).\n" + + output_string += "# It uses the local coordinate system.\n" + output_string += "# For example to create a 2 x 1.5 x 3 meters box " + output_string += "use the following: box length: 2.0, 1.5, 3.0\n" + + output_string += "box length x: " + str(self.box_length[0]) + "\n" + output_string += "box length y: " + str(self.box_length[1]) + "\n" + output_string += "box length z: " + str(self.box_length[2]) + "\n" + + output_string += "\n# box origin indicates the x, y, and z coordinates " + output_string += "of the origin of the FFD bounding box. That is " + output_string += "center of\n" + + output_string += "# rotation of the bounding box. It corresponds to " + output_string += "the point coordinates with position [0][0][0].\n" output_string += '# See section "Parameters weights" for more ' - output_string += 'details.\n' - output_string += '# For example, if the origin is equal to 0., 0., 0., ' - output_string += 'use the following: box origin: 0., 0., 0.\n' - - output_string += 'box origin x: ' + str(self.box_origin[0]) + '\n' - output_string += 'box origin y: ' + str(self.box_origin[1]) + '\n' - output_string += 'box origin z: ' + str(self.box_origin[2]) + '\n' - - output_string += '\n# rotation angle indicates the rotation angle ' - output_string += 'around the x, y, and z axis of the FFD bounding box ' - output_string += 'in degrees.\n' - - output_string += '# The rotation is done with respect to the box ' - output_string += 'origin.\n' - output_string += '# For example, to rotate the box by 2 deg along ' - output_string += 'the z ' - output_string += 'direction, use the following: rotation angle: ' - output_string += '0., 0., 2.\n' - - output_string += 'rotation angle x: ' + str(self.rot_angle[0]) + '\n' - output_string += 'rotation angle y: ' + str(self.rot_angle[1]) + '\n' - output_string += 'rotation angle z: ' + str(self.rot_angle[2]) + '\n' - - output_string += '\n\n[Parameters weights]\n' - output_string += '# This section describes the weights of the FFD ' - output_string += 'control points.\n' - - output_string += '# We adopt the following convention:\n' - output_string += '# For example with a 2x2x2 grid of control points we ' - output_string += 'have to fill a 2x2x2 matrix of weights.\n' - - output_string += '# If a weight is equal to zero you can discard the ' - output_string += 'line since the default is zero.\n' - - output_string += '#\n' - output_string += '# | x index | y index | z index | weight |\n' - output_string += '# --------------------------------------\n' - output_string += '# | 0 | 0 | 0 | 1.0 |\n' - output_string += '# | 0 | 1 | 1 | 0.0 | --> you ' - output_string += 'can erase this line without effects\n' - output_string += '# | 0 | 1 | 0 | -2.1 |\n' - output_string += '# | 0 | 0 | 1 | 3.4 |\n' - - output_string += '\n# parameter x collects the displacements along x, ' - output_string += 'normalized with the box length x.' - - output_string += '\nparameter x:' + output_string += "details.\n" + output_string += "# For example, if the origin is equal to 0., 0., 0., " + output_string += "use the following: box origin: 0., 0., 0.\n" + + output_string += "box origin x: " + str(self.box_origin[0]) + "\n" + output_string += "box origin y: " + str(self.box_origin[1]) + "\n" + output_string += "box origin z: " + str(self.box_origin[2]) + "\n" + + output_string += "\n# rotation angle indicates the rotation angle " + output_string += "around the x, y, and z axis of the FFD bounding box " + output_string += "in degrees.\n" + + output_string += "# The rotation is done with respect to the box " + output_string += "origin.\n" + output_string += "# For example, to rotate the box by 2 deg along " + output_string += "the z " + output_string += "direction, use the following: rotation angle: " + output_string += "0., 0., 2.\n" + + output_string += "rotation angle x: " + str(self.rot_angle[0]) + "\n" + output_string += "rotation angle y: " + str(self.rot_angle[1]) + "\n" + output_string += "rotation angle z: " + str(self.rot_angle[2]) + "\n" + + output_string += "\n\n[Parameters weights]\n" + output_string += "# This section describes the weights of the FFD " + output_string += "control points.\n" + + output_string += "# We adopt the following convention:\n" + output_string += "# For example with a 2x2x2 grid of control points we " + output_string += "have to fill a 2x2x2 matrix of weights.\n" + + output_string += "# If a weight is equal to zero you can discard the " + output_string += "line since the default is zero.\n" + + output_string += "#\n" + output_string += "# | x index | y index | z index | weight |\n" + output_string += "# --------------------------------------\n" + output_string += "# | 0 | 0 | 0 | 1.0 |\n" + output_string += "# | 0 | 1 | 1 | 0.0 | --> you " + output_string += "can erase this line without effects\n" + output_string += "# | 0 | 1 | 0 | -2.1 |\n" + output_string += "# | 0 | 0 | 1 | 3.4 |\n" + + output_string += "\n# parameter x collects the displacements along x, " + output_string += "normalized with the box length x." + + output_string += "\nparameter x:" offset = 1 for i in range(0, self.n_control_points[0]): for j in range(0, self.n_control_points[1]): for k in range(0, self.n_control_points[2]): - output_string += offset * ' ' + str(i) + ' ' + str( - j) + ' ' + str(k) + ' ' + str( - self.array_mu_x[i][j][k]) + '\n' + output_string += ( + offset * " " + + str(i) + + " " + + str(j) + + " " + + str(k) + + " " + + str(self.array_mu_x[i][j][k]) + + "\n" + ) offset = 13 - output_string += '\n# parameter y collects the displacements along y, ' - output_string += 'normalized with the box length y.' + output_string += "\n# parameter y collects the displacements along y, " + output_string += "normalized with the box length y." - output_string += '\nparameter y:' + output_string += "\nparameter y:" offset = 1 for i in range(0, self.n_control_points[0]): for j in range(0, self.n_control_points[1]): for k in range(0, self.n_control_points[2]): - output_string += offset * ' ' + str(i) + ' ' + str( - j) + ' ' + str(k) + ' ' + str( - self.array_mu_y[i][j][k]) + '\n' + output_string += ( + offset * " " + + str(i) + + " " + + str(j) + + " " + + str(k) + + " " + + str(self.array_mu_y[i][j][k]) + + "\n" + ) offset = 13 - output_string += '\n# parameter z collects the displacements along z, ' - output_string += 'normalized with the box length z.' + output_string += "\n# parameter z collects the displacements along z, " + output_string += "normalized with the box length z." - output_string += '\nparameter z:' + output_string += "\nparameter z:" offset = 1 for i in range(0, self.n_control_points[0]): for j in range(0, self.n_control_points[1]): for k in range(0, self.n_control_points[2]): - output_string += offset * ' ' + str(i) + ' ' + str( - j) + ' ' + str(k) + ' ' + str( - self.array_mu_z[i][j][k]) + '\n' + output_string += ( + offset * " " + + str(i) + + " " + + str(j) + + " " + + str(k) + + " " + + str(self.array_mu_z[i][j][k]) + + "\n" + ) offset = 13 - with open(filename, 'w') as f: + with open(filename, "w") as f: f.write(output_string) def __str__(self): @@ -450,16 +490,16 @@ def __str__(self): for debugging. """ string = "" - string += 'conversion_unit = {}\n'.format(self.conversion_unit) - string += 'n_control_points = {}\n\n'.format(self.n_control_points) - string += 'box_length = {}\n'.format(self.box_length) - string += 'box_origin = {}\n'.format(self.box_origin) - string += 'rot_angle = {}\n'.format(self.rot_angle) - string += '\narray_mu_x =\n{}\n'.format(self.array_mu_x) - string += '\narray_mu_y =\n{}\n'.format(self.array_mu_y) - string += '\narray_mu_z =\n{}\n'.format(self.array_mu_z) - string += '\nrotation_matrix = \n{}\n'.format(self.rotation_matrix) - string += '\nposition_vertices = {}\n'.format(self.position_vertices) + string += "conversion_unit = {}\n".format(self.conversion_unit) + string += "n_control_points = {}\n\n".format(self.n_control_points) + string += "box_length = {}\n".format(self.box_length) + string += "box_origin = {}\n".format(self.box_origin) + string += "rot_angle = {}\n".format(self.rot_angle) + string += "\narray_mu_x =\n{}\n".format(self.array_mu_x) + string += "\narray_mu_y =\n{}\n".format(self.array_mu_y) + string += "\narray_mu_z =\n{}\n".format(self.array_mu_z) + string += "\nrotation_matrix = \n{}\n".format(self.rotation_matrix) + string += "\nposition_vertices = {}\n".format(self.position_vertices) return string def control_points(self, deformed=True): @@ -480,21 +520,23 @@ def control_points(self, deformed=True): y_coords, x_coords, z_coords = np.meshgrid(y, x, z) box_points = np.array( - [x_coords.ravel(), - y_coords.ravel(), - z_coords.ravel()]) + [x_coords.ravel(), y_coords.ravel(), z_coords.ravel()] + ) if deformed: - box_points += np.array([ - self.array_mu_x.ravel() * self.box_length[0], - self.array_mu_y.ravel() * self.box_length[1], - self.array_mu_z.ravel() * self.box_length[2] - ]) + box_points += np.array( + [ + self.array_mu_x.ravel() * self.box_length[0], + self.array_mu_y.ravel() * self.box_length[1], + self.array_mu_z.ravel() * self.box_length[2], + ] + ) n_rows = box_points.shape[1] box_points = np.dot(self.rotation_matrix, box_points) + np.transpose( - np.tile(self.box_origin, (n_rows, 1))) + np.tile(self.box_origin, (n_rows, 1)) + ) return box_points.T @@ -518,7 +560,7 @@ def reflect(self, axis=0, in_place=True): Default is True. :return: a new object with the same parameters and the reflected lattice if `in_place` is False, otherwise NoneType. - + :Example: >>> ffd.reflect(axis=0, in_place=True) # irreversible @@ -528,17 +570,20 @@ def reflect(self, axis=0, in_place=True): # check axis value if axis not in (0, 1, 2): raise ValueError( - "The axis has to be 0, 1, or 2. Current value {}.".format(axis)) + "The axis has to be 0, 1, or 2. Current value {}.".format(axis) + ) # check that the plane of symmetry is undeformed - if (axis == 0 and np.count_nonzero(self.array_mu_x[-1, :, :]) != 0) or ( - axis == 1 and np.count_nonzero(self.array_mu_y[:, -1, :]) != 0 - ) or (axis == 2 and np.count_nonzero(self.array_mu_z[:, :, -1]) != 0): + if ( + (axis == 0 and np.count_nonzero(self.array_mu_x[-1, :, :]) != 0) + or (axis == 1 and np.count_nonzero(self.array_mu_y[:, -1, :]) != 0) + or (axis == 2 and np.count_nonzero(self.array_mu_z[:, :, -1]) != 0) + ): raise RuntimeError( - "If you want to reflect the FFD bounding box along axis " + \ - "{} you can not diplace the control ".format(axis) + \ - "points in the symmetry plane along that axis." - ) + "If you want to reflect the FFD bounding box along axis " + + "{} you can not diplace the control ".format(axis) + + "points in the symmetry plane along that axis." + ) if in_place is False: self = copy.deepcopy(self) @@ -559,54 +604,61 @@ def reflect(self, axis=0, in_place=True): # we append along the given axis all the displacements reflected # and in the reverse order - self.array_mu_x = np.append(self.array_mu_x, - reflection[0] * - np.flip(self.array_mu_x, axis)[indeces], - axis=axis) - self.array_mu_y = np.append(self.array_mu_y, - reflection[1] * - np.flip(self.array_mu_y, axis)[indeces], - axis=axis) - self.array_mu_z = np.append(self.array_mu_z, - reflection[2] * - np.flip(self.array_mu_z, axis)[indeces], - axis=axis) + self.array_mu_x = np.append( + self.array_mu_x, + reflection[0] * np.flip(self.array_mu_x, axis)[indeces], + axis=axis, + ) + self.array_mu_y = np.append( + self.array_mu_y, + reflection[1] * np.flip(self.array_mu_y, axis)[indeces], + axis=axis, + ) + self.array_mu_z = np.append( + self.array_mu_z, + reflection[2] * np.flip(self.array_mu_z, axis)[indeces], + axis=axis, + ) if in_place is False: return self - - def __call__(self, src_pts): """ This method performs the FFD to `src_pts` and return the deformed points. - + :param numpy.ndarray src_pts: the array of dimensions (*n_points*, *3*) containing the points to deform. The points have to be arranged by row. :return: the deformed points :rtype: numpy.ndarray (with shape = (*n_points*, *3*)) """ + def is_inside(pts, boundaries): - """ + """ Check is `pts` is inside the ranges provided by `boundaries`. """ - return np.all(np.logical_and(pts >= boundaries[0], - pts <= boundaries[1]), - axis=1) + return np.all( + np.logical_and(pts >= boundaries[0], pts <= boundaries[1]), + axis=1, + ) # map to the reference domain src_reference_frame_pts = self.psi(src_pts - self.box_origin) # apply deformation for all the pts in the unit cube - index_pts_inside = is_inside(src_reference_frame_pts, - np.array([[0., 0., 0.], [1., 1., 1.]])) + index_pts_inside = is_inside( + src_reference_frame_pts, + np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]), + ) shifted_reference_frame_pts = self.T( - src_reference_frame_pts[index_pts_inside]) + src_reference_frame_pts[index_pts_inside] + ) # map to the physical domain - shifted_pts = self.inverse_psi( - shifted_reference_frame_pts) + self.box_origin + shifted_pts = ( + self.inverse_psi(shifted_reference_frame_pts) + self.box_origin + ) dst_pts = src_pts.copy() dst_pts[index_pts_inside] = shifted_pts diff --git a/pygem/filehandler.py b/pygem/filehandler.py index f56811e5..b8568b8f 100644 --- a/pygem/filehandler.py +++ b/pygem/filehandler.py @@ -3,11 +3,15 @@ .. warning:: This module will be deprecated in next releases. Follow updates on - https://github.com/mathLab for news about file handling. + https://github.com/mathLab for news about file handling. """ + import os import warnings -warnings.warn("This module will be deprecated in next releases", DeprecationWarning) + +warnings.warn( + "This module will be deprecated in next releases", DeprecationWarning +) class FileHandler(object): @@ -32,8 +36,10 @@ def parse(self, *args): Not implemented, it has to be implemented in subclasses. """ raise NotImplementedError( - "Subclass must implement abstract method " \ - + self.__class__.__name__ + ".parse") + "Subclass must implement abstract method " + + self.__class__.__name__ + + ".parse" + ) def write(self, *args): """ @@ -42,8 +48,10 @@ def write(self, *args): Not implemented, it has to be implemented in subclasses. """ raise NotImplementedError( - "Subclass must implement abstract method " \ - + self.__class__.__name__ + ".write") + "Subclass must implement abstract method " + + self.__class__.__name__ + + ".write" + ) def _check_extension(self, filename): """ @@ -55,9 +63,11 @@ def _check_extension(self, filename): __, file_ext = os.path.splitext(filename) if file_ext not in self.extensions: raise ValueError( - 'The input file does not have the proper extension. \ - It is {0!s}, instead of {1!s}.'.format(file_ext, - self.extensions)) + "The input file does not have the proper extension. \ + It is {0!s}, instead of {1!s}.".format( + file_ext, self.extensions + ) + ) @staticmethod def _check_filename_type(filename): @@ -69,7 +79,8 @@ def _check_filename_type(filename): """ if not isinstance(filename, str): raise TypeError( - 'The given filename ({0!s}) must be a string'.format(filename)) + "The given filename ({0!s}) must be a string".format(filename) + ) def _check_infile_instantiation(self): """ @@ -80,4 +91,5 @@ def _check_infile_instantiation(self): """ if not self.infile: raise RuntimeError( - "You can not write a file without having parsed one.") + "You can not write a file without having parsed one." + ) diff --git a/pygem/idw.py b/pygem/idw.py index a7f979c8..3c77bfbf 100644 --- a/pygem/idw.py +++ b/pygem/idw.py @@ -36,8 +36,10 @@ :math:`\\mathrm{x}_i` and :math:`p` is a power parameter, typically equal to 2. """ + import os import numpy as np + try: import configparser as configparser except ImportError: @@ -60,7 +62,7 @@ class IDW(Deformation): (*n_control_points*, *3*) array with the coordinates of the interpolation control points after the deformation. The default is the vertices of the unit cube. - + :cvar int power: the power parameter. The default value is 2. :cvar numpy.ndarray original_control_points: it is an (*n_control_points*, *3*) array with the coordinates of the original @@ -86,26 +88,43 @@ class IDW(Deformation): >>> idw.read_parameters('tests/test_datasets/parameters_idw_cube.prm') >>> new_mesh_points = idw(mesh_points.T) """ - def __init__(self, - original_control_points=None, - deformed_control_points=None, - power=2): + + def __init__( + self, + original_control_points=None, + deformed_control_points=None, + power=2, + ): if original_control_points is None: - self.original_control_points = np.array([[0., 0., 0.], [0., 0., 1.], - [0., 1., 0.], [1., 0., 0.], - [0., 1., 1.], [1., 0., 1.], - [1., 1., 0.], [1., 1., - 1.]]) + self.original_control_points = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + ] + ) else: self.original_control_points = original_control_points if deformed_control_points is None: - self.deformed_control_points = np.array([[0., 0., 0.], [0., 0., 1.], - [0., 1., 0.], [1., 0., 0.], - [0., 1., 1.], [1., 0., 1.], - [1., 1., 0.], [1., 1., - 1.]]) + self.deformed_control_points = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + ] + ) else: self.deformed_control_points = deformed_control_points @@ -116,8 +135,9 @@ def __call__(self, src_pts): This method performs the deformation of the mesh points. After the execution it sets `self.modified_mesh_points`. """ + def distance(u, v): - """ Norm of u - v """ + """Norm of u - v""" return np.linalg.norm(u - v, ord=self.power) # Compute displacement of the control points @@ -130,12 +150,14 @@ def distance(u, v): # not zero, otherwise 1.0 where distance is zero. weights = np.zeros(dist.shape) for i, d in enumerate(dist): - weights[i] = 1. / d if d.all() else np.where(d == 0.0, 1.0, 0.0) + weights[i] = 1.0 / d if d.all() else np.where(d == 0.0, 1.0, 0.0) - offset = np.array([ - np.sum(displ * wi[:, np.newaxis] / np.sum(wi), axis=0) - for wi in weights - ]) + offset = np.array( + [ + np.sum(displ * wi[:, np.newaxis] / np.sum(wi), axis=0) + for wi in weights + ] + ) return src_pts + offset @@ -146,23 +168,25 @@ def read_parameters(self, filename): :param string filename: parameters file to be read in. """ if not isinstance(filename, str): - raise TypeError('filename must be a string') + raise TypeError("filename must be a string") if not os.path.isfile(filename): - raise IOError('filename does not exist') + raise IOError("filename does not exist") config = configparser.RawConfigParser() config.read(filename) - self.power = config.getint('Inverse Distance Weighting', 'power') + self.power = config.getint("Inverse Distance Weighting", "power") - ctrl_points = config.get('Control points', 'original control points') + ctrl_points = config.get("Control points", "original control points") self.original_control_points = np.array( - [line.split() for line in ctrl_points.split('\n')], dtype=float) + [line.split() for line in ctrl_points.split("\n")], dtype=float + ) - defo_points = config.get('Control points', 'deformed control points') + defo_points = config.get("Control points", "deformed control points") self.deformed_control_points = np.array( - [line.split() for line in defo_points.split('\n')], dtype=float) + [line.split() for line in defo_points.split("\n")], dtype=float + ) def write_parameters(self, filename): """ @@ -188,20 +212,22 @@ def write_parameters(self, filename): output_string += "original control points: " output_string += ( - ' '.join(map(str, self.original_control_points[0])) + "\n") + " ".join(map(str, self.original_control_points[0])) + "\n" + ) for points in self.original_control_points[1:]: - output_string += 25 * ' ' + ' '.join(map(str, points)) + "\n" + output_string += 25 * " " + " ".join(map(str, points)) + "\n" output_string += "\n" output_string += "# deformed control points collects the coordinates\n" output_string += "# of the interpolation control points after the\n" output_string += "# deformation.\n" output_string += "deformed control points: " output_string += ( - ' '.join(map(str, self.original_control_points[0])) + "\n") + " ".join(map(str, self.original_control_points[0])) + "\n" + ) for points in self.deformed_control_points[1:]: - output_string += 25 * ' ' + ' '.join(map(str, points)) + "\n" + output_string += 25 * " " + " ".join(map(str, points)) + "\n" - with open(filename, 'w') as f: + with open(filename, "w") as f: f.write(output_string) def __str__(self): @@ -209,10 +235,10 @@ def __str__(self): This method prints all the IDW parameters on the screen. Its purpose is for debugging. """ - string = '' - string += 'p = {}\n'.format(self.power) - string += '\noriginal_control_points =\n' - string += '{}\n'.format(self.original_control_points) - string += '\ndeformed_control_points =\n' - string += '{}\n'.format(self.deformed_control_points) + string = "" + string += "p = {}\n".format(self.power) + string += "\noriginal_control_points =\n" + string += "{}\n".format(self.original_control_points) + string += "\ndeformed_control_points =\n" + string += "{}\n".format(self.deformed_control_points) return string diff --git a/pygem/khandler.py b/pygem/khandler.py index 1061ea78..d3d25537 100644 --- a/pygem/khandler.py +++ b/pygem/khandler.py @@ -1,7 +1,8 @@ """ Derived module from filehandler.py to handle LS-DYNA keyword (.k) files. """ -import re + +import re import numpy as np import pygem.filehandler as fh @@ -18,7 +19,7 @@ class KHandler(fh.FileHandler): def __init__(self): super(KHandler, self).__init__() - self.extensions = ['.k'] + self.extensions = [".k"] def parse(self, filename): """ @@ -39,29 +40,32 @@ def parse(self, filename): mesh_points = [] node_indicator = False - with open(self.infile, 'r') as input_file: + with open(self.infile, "r") as input_file: for num, line in enumerate(input_file): - expression = re.compile(r'(.+?)(?:,|$)') + expression = re.compile(r"(.+?)(?:,|$)") expression = expression.findall(line) if line.startswith("$"): continue - if line.startswith('*NODE'): + if line.startswith("*NODE"): node_indicator = True continue - if line.startswith('*ELEMENT'): + if line.startswith("*ELEMENT"): break if not node_indicator: pass else: if len(expression) == 1: - expression = re.findall(r'\S+', expression[0]) - l = [float(expression[1]), float(expression[2]), - float(expression[3])] + expression = re.findall(r"\S+", expression[0]) + l = [ + float(expression[1]), + float(expression[2]), + float(expression[3]), + ] mesh_points.append(l) mesh_points = np.array(mesh_points) @@ -85,22 +89,22 @@ def write(self, mesh_points, filename): i = 0 node_indicator = False - with open(self.outfile, 'w') as output_file: - with open(self.infile, 'r') as input_file: + with open(self.outfile, "w") as output_file: + with open(self.infile, "r") as input_file: for _, line in enumerate(input_file): - get_num = re.findall(r'[-+]?[0-9]*\.?[0-9]+', line) + get_num = re.findall(r"[-+]?[0-9]*\.?[0-9]+", line) # Write header files - if line.startswith('$'): + if line.startswith("$"): output_file.write(line) continue # Change the node indicator if you find the elements section - if line.startswith('*ELEMENT'): + if line.startswith("*ELEMENT"): node_indicator = False # Change the nodes indicator if you find the nodes section - if line.startswith('*NODE'): + if line.startswith("*NODE"): node_indicator = True output_file.write(line) continue @@ -114,10 +118,12 @@ def write(self, mesh_points, filename): split_line = line.split(" ") # Format the data into correct format - data = [int(get_num[0]), - '{:.10f}'.format(float(mesh_points[i][0])), - '{:.10f}'.format(float(mesh_points[i][1])), - '{:.10f}'.format(float(mesh_points[i][2]))] + data = [ + int(get_num[0]), + "{:.10f}".format(float(mesh_points[i][0])), + "{:.10f}".format(float(mesh_points[i][1])), + "{:.10f}".format(float(mesh_points[i][2])), + ] comma_seperator = False pointer = 0 @@ -128,13 +134,20 @@ def write(self, mesh_points, filename): if value[len(value) - 1] == ",": comma_seperator = True - new_str = value.replace(value[:-1], str(data[pointer])) + new_str = value.replace( + value[:-1], str(data[pointer]) + ) split_line[index] = new_str else: - new_str = value.replace(value, str(data[pointer])) + new_str = value.replace( + value, str(data[pointer]) + ) split_line[index] = new_str - if float(data[pointer]) < 0 and not comma_seperator: + if ( + float(data[pointer]) < 0 + and not comma_seperator + ): del split_line[index - 1] pointer += 1 diff --git a/pygem/mdpahandler.py b/pygem/mdpahandler.py index ac341933..e1036a25 100644 --- a/pygem/mdpahandler.py +++ b/pygem/mdpahandler.py @@ -1,6 +1,7 @@ """ Derived module from khandler.py to handle Kratos mesh (.mdpa) files. """ + import numpy as np import pygem.filehandler as fh @@ -17,7 +18,7 @@ class MdpaHandler(fh.FileHandler): def __init__(self): super(MdpaHandler, self).__init__() - self.extensions = ['.mdpa'] + self.extensions = [".mdpa"] def parse(self, filename): """ @@ -36,15 +37,15 @@ def parse(self, filename): self.infile = filename index = -9 mesh_points = [] - with open(self.infile, 'r') as input_file: + with open(self.infile, "r") as input_file: for num, line in enumerate(input_file): - if line.startswith('Begin Nodes'): + if line.startswith("Begin Nodes"): index = num if num == index + 1: - if line.startswith('End Nodes'): + if line.startswith("End Nodes"): break else: - line = line.replace('D', 'E') + line = line.replace("D", "E") li = [] for t in line.split()[1:]: try: @@ -72,19 +73,23 @@ def write(self, mesh_points, filename): self.outfile = filename index = -9 i = 0 - with open(self.outfile, 'w') as output_file: - with open(self.infile, 'r') as input_file: + with open(self.outfile, "w") as output_file: + with open(self.infile, "r") as input_file: for num, line in enumerate(input_file): - if line.startswith('Begin Nodes'): + if line.startswith("Begin Nodes"): index = num if num == index + 1: - if line.startswith('End Nodes'): + if line.startswith("End Nodes"): index = -9 else: - line = (" {:6d} {:23.16E} {:23.16E} {:23.16E}\n" - .format(i+1, mesh_points[i][0], - mesh_points[i][1], - mesh_points[i][2])) + line = ( + " {:6d} {:23.16E} {:23.16E} {:23.16E}\n".format( + i + 1, + mesh_points[i][0], + mesh_points[i][1], + mesh_points[i][2], + ) + ) i += 1 index = num output_file.write(line) diff --git a/pygem/meta.py b/pygem/meta.py index e4f97ee7..87a40e85 100644 --- a/pygem/meta.py +++ b/pygem/meta.py @@ -12,7 +12,7 @@ def get_current_year(): - """ Return current year """ + """Return current year""" from datetime import datetime return datetime.now().year diff --git a/pygem/openfhandler.py b/pygem/openfhandler.py index d21485ab..6af7cbf9 100644 --- a/pygem/openfhandler.py +++ b/pygem/openfhandler.py @@ -3,12 +3,16 @@ .. warning:: This module will be deprecated in next releases. Follow updates on - https://github.com/mathLab for news about file handling. + https://github.com/mathLab for news about file handling. """ + import numpy as np import pygem.filehandler as fh import warnings -warnings.warn("This module will be deprecated in next releases", DeprecationWarning) + +warnings.warn( + "This module will be deprecated in next releases", DeprecationWarning +) class OpenFoamHandler(fh.FileHandler): @@ -23,7 +27,7 @@ class OpenFoamHandler(fh.FileHandler): def __init__(self): super(OpenFoamHandler, self).__init__() - self.extensions = [''] + self.extensions = [""] def parse(self, filename): """ @@ -31,7 +35,7 @@ def parse(self, filename): the coordinates. :param string filename: name of the input file. - + :return: mesh_points: it is a `n_points`-by-3 matrix containing the coordinates of the points of the mesh :rtype: numpy.ndarray @@ -47,14 +51,14 @@ def parse(self, filename): nrow = 0 i = 0 - with open(self.infile, 'r') as input_file: + with open(self.infile, "r") as input_file: for line in input_file: nrow += 1 if nrow == 19: n_points = int(line) mesh_points = np.zeros(shape=(n_points, 3)) if 20 < nrow < 21 + n_points: - line = line[line.index("(") + 1:line.rindex(")")] + line = line[line.index("(") + 1 : line.rindex(")")] j = 0 for number in line.split(): mesh_points[i][j] = float(number) @@ -85,14 +89,23 @@ def write(self, mesh_points, filename): n_points = mesh_points.shape[0] nrow = 0 i = 0 - with open(self.infile, 'r') as input_file, open(self.outfile, - 'w') as output_file: + with ( + open(self.infile, "r") as input_file, + open(self.outfile, "w") as output_file, + ): for line in input_file: nrow += 1 if 20 < nrow < 21 + n_points: - output_file.write('(' + str(mesh_points[i][0]) + ' ' + str( - mesh_points[i][1]) + ' ' + str(mesh_points[i][2]) + ')') - output_file.write('\n') + output_file.write( + "(" + + str(mesh_points[i][0]) + + " " + + str(mesh_points[i][1]) + + " " + + str(mesh_points[i][2]) + + ")" + ) + output_file.write("\n") i += 1 else: output_file.write(line) diff --git a/pygem/rbf.py b/pygem/rbf.py index cb1f6fff..4741f9b4 100644 --- a/pygem/rbf.py +++ b/pygem/rbf.py @@ -24,7 +24,7 @@ is defines as follows .. math:: - \\mathcal{M}(\\boldsymbol{x}) = p(\\boldsymbol{x}) + + \\mathcal{M}(\\boldsymbol{x}) = p(\\boldsymbol{x}) + \\sum_{i=1}^{\\mathcal{N}_C} \\gamma_i \\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}} \\|) @@ -55,8 +55,10 @@ Wendland :math:`C^2` basis and Polyharmonic splines all defined and implemented below. """ + import os import numpy as np + try: import configparser as configparser except ImportError: @@ -134,25 +136,27 @@ class RBF(Deformation): # Precision mapping DTYPE_MAP = { - 'fp16': np.float16, - 'float16': np.float16, - 'fp32': np.float32, - 'float32': np.float32, - 'fp64': np.float64, - 'float64': np.float64, - 'fp96': np.float96 if hasattr(np, 'float96') else np.float64, - 'float96': np.float96 if hasattr(np, 'float96') else np.float64, - 'fp128': np.float128 if hasattr(np, 'float128') else np.float64, - 'float128': np.float128 if hasattr(np, 'float128') else np.float64, + "fp16": np.float16, + "float16": np.float16, + "fp32": np.float32, + "float32": np.float32, + "fp64": np.float64, + "float64": np.float64, + "fp96": np.float96 if hasattr(np, "float96") else np.float64, + "float96": np.float96 if hasattr(np, "float96") else np.float64, + "fp128": np.float128 if hasattr(np, "float128") else np.float64, + "float128": np.float128 if hasattr(np, "float128") else np.float64, } - def __init__(self, - original_control_points=None, - deformed_control_points=None, - func='gaussian_spline', - radius=0.5, - extra_parameter=None, - dtype='fp64'): + def __init__( + self, + original_control_points=None, + deformed_control_points=None, + func="gaussian_spline", + radius=0.5, + extra_parameter=None, + dtype="fp64", + ): # Parse and set dtype with platform check if isinstance(dtype, str): @@ -164,23 +168,23 @@ def __init__(self, ) # Check for fp128 fallback - if dtype_lower in ['fp128', 'float128']: - if not hasattr(np, 'float128'): + if dtype_lower in ["fp128", "float128"]: + if not hasattr(np, "float128"): warnings.warn( "fp128/float128 is not supported on this platform. " "Automatically falling back to fp64. " "For true quad-precision, consider using Linux platform.", - RuntimeWarning + RuntimeWarning, ) # Check for fp96 fallback - if dtype_lower in ['fp96', 'float96']: - if not hasattr(np, 'float96'): + if dtype_lower in ["fp96", "float96"]: + if not hasattr(np, "float96"): warnings.warn( "fp96/float96 is not supported on this platform. " "Automatically falling back to fp64. " "For higher precision consider using 'fp128' (if available) ", - RuntimeWarning + RuntimeWarning, ) self._dtype = self.DTYPE_MAP[dtype_lower] @@ -191,27 +195,48 @@ def __init__(self, self.radius = radius if original_control_points is None: - self.original_control_points = np.array([[0., 0., 0.], [0., 0., 1.], - [0., 1., 0.], [1., 0., 0.], - [0., 1., 1.], [1., 0., 1.], - [1., 1., 0.], [1., 1., 1.]], - dtype=self._dtype) + self.original_control_points = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + ], + dtype=self._dtype, + ) else: - self.original_control_points = np.asarray(original_control_points, dtype=self._dtype) + self.original_control_points = np.asarray( + original_control_points, dtype=self._dtype + ) if deformed_control_points is None: - self.deformed_control_points = np.array([[0., 0., 0.], [0., 0., 1.], - [0., 1., 0.], [1., 0., 0.], - [0., 1., 1.], [1., 0., 1.], - [1., 1., 0.], [1., 1., 1.]], - dtype=self._dtype) + self.deformed_control_points = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + [1.0, 1.0, 1.0], + ], + dtype=self._dtype, + ) else: - self.deformed_control_points = np.asarray(deformed_control_points, dtype=self._dtype) + self.deformed_control_points = np.asarray( + deformed_control_points, dtype=self._dtype + ) self.extra = extra_parameter if extra_parameter else dict() - self.weights = self._get_weights(self.original_control_points, - self.deformed_control_points) + self.weights = self._get_weights( + self.original_control_points, self.deformed_control_points + ) @property def n_control_points(self): @@ -243,7 +268,7 @@ def basis(self, func): elif isinstance(func, str): self.__basis = RBFFactory(func) else: - raise TypeError('`func` is not valid.') + raise TypeError("`func` is not valid.") def _get_weights(self, X, Y): """ @@ -280,11 +305,17 @@ def _get_weights(self, X, Y): rhs = np.zeros((size, dim), dtype=self._dtype) rhs[:npts, :] = Y - solve_dtype = np.float64 if self._dtype not in (np.float32, np.float64) else self._dtype - weights = np.linalg.solve(H.astype(solve_dtype), rhs.astype(solve_dtype)).astype(self._dtype) + solve_dtype = ( + np.float64 + if self._dtype not in (np.float32, np.float64) + else self._dtype + ) + weights = np.linalg.solve( + H.astype(solve_dtype), rhs.astype(solve_dtype) + ).astype(self._dtype) return weights.astype(self._dtype) - def read_parameters(self, filename='parameters_rbf.prm'): + def read_parameters(self, filename="parameters_rbf.prm"): """ Reads in the parameters file and fill the self structure. @@ -292,7 +323,7 @@ def read_parameters(self, filename='parameters_rbf.prm'): parameters_rbf.prm. """ if not isinstance(filename, str): - raise TypeError('filename must be a string') + raise TypeError("filename must be a string") # Checks if the parameters file exists. If not it writes the default # class into filename. It consists in the vetices of a cube of side one @@ -304,29 +335,33 @@ def read_parameters(self, filename='parameters_rbf.prm'): config = configparser.RawConfigParser() config.read(filename) - rbf_settings = dict(config.items('Radial Basis Functions')) + rbf_settings = dict(config.items("Radial Basis Functions")) - self.basis = rbf_settings.pop('basis function') - self.radius = float(rbf_settings.pop('radius')) + self.basis = rbf_settings.pop("basis function") + self.radius = float(rbf_settings.pop("radius")) self.extra = {k: eval(v) for k, v in rbf_settings.items()} - ctrl_points = config.get('Control points', 'original control points') - lines = ctrl_points.split('\n') + ctrl_points = config.get("Control points", "original control points") + lines = ctrl_points.split("\n") self.original_control_points = np.array( - list(map(lambda x: x.split(), lines)), dtype=self._dtype) + list(map(lambda x: x.split(), lines)), dtype=self._dtype + ) - mod_points = config.get('Control points', 'deformed control points') - lines = mod_points.split('\n') + mod_points = config.get("Control points", "deformed control points") + lines = mod_points.split("\n") self.deformed_control_points = np.array( - list(map(lambda x: x.split(), lines)), dtype=self._dtype) + list(map(lambda x: x.split(), lines)), dtype=self._dtype + ) if len(lines) != self.n_control_points: - raise TypeError("The number of control points must be equal both in" - "the 'original control points' and in the 'deformed" - "control points' section of the parameters file" - "({0!s})".format(filename)) - - def write_parameters(self, filename='parameters_rbf.prm'): + raise TypeError( + "The number of control points must be equal both in" + "the 'original control points' and in the 'deformed" + "control points' section of the parameters file" + "({0!s})".format(filename) + ) + + def write_parameters(self, filename="parameters_rbf.prm"): """ This method writes a parameters file (.prm) called `filename` and fills it with all the parameters class members. Default value is @@ -338,57 +373,67 @@ def write_parameters(self, filename='parameters_rbf.prm'): raise TypeError("filename must be a string") output_string = "" - output_string += '\n[Radial Basis Functions]\n' - output_string += '# This section describes the radial basis functions' - output_string += ' shape.\n' - - output_string += '\n# basis funtion is the name of the basis functions' - output_string += ' to use in the transformation. The functions\n' - output_string += '# implemented so far are: gaussian_spline,' - output_string += ' multi_quadratic_biharmonic_spline,\n' - output_string += '# inv_multi_quadratic_biharmonic_spline,' - output_string += ' thin_plate_spline, beckert_wendland_c2_basis,' - output_string += ' polyharmonic_spline.\n' - output_string += '# For a comprehensive list with details see the' - output_string += ' class RBF.\n' - output_string += 'basis function: {}\n'.format('gaussian_spline') - - output_string += '\n# radius is the scaling parameter r that affects' - output_string += ' the shape of the basis functions. See the' - output_string += ' documentation\n' - output_string += '# of the class RBF for details.\n' - output_string += 'radius: {}\n'.format(str(self.radius)) - - output_string += '\n\n[Control points]\n' - output_string += '# This section describes the RBF control points.\n' - - output_string += '\n# original control points collects the coordinates' - output_string += ' of the interpolation control points before the' - output_string += ' deformation.\n' - - output_string += 'original control points:' + output_string += "\n[Radial Basis Functions]\n" + output_string += "# This section describes the radial basis functions" + output_string += " shape.\n" + + output_string += "\n# basis funtion is the name of the basis functions" + output_string += " to use in the transformation. The functions\n" + output_string += "# implemented so far are: gaussian_spline," + output_string += " multi_quadratic_biharmonic_spline,\n" + output_string += "# inv_multi_quadratic_biharmonic_spline," + output_string += " thin_plate_spline, beckert_wendland_c2_basis," + output_string += " polyharmonic_spline.\n" + output_string += "# For a comprehensive list with details see the" + output_string += " class RBF.\n" + output_string += "basis function: {}\n".format("gaussian_spline") + + output_string += "\n# radius is the scaling parameter r that affects" + output_string += " the shape of the basis functions. See the" + output_string += " documentation\n" + output_string += "# of the class RBF for details.\n" + output_string += "radius: {}\n".format(str(self.radius)) + + output_string += "\n\n[Control points]\n" + output_string += "# This section describes the RBF control points.\n" + + output_string += "\n# original control points collects the coordinates" + output_string += " of the interpolation control points before the" + output_string += " deformation.\n" + + output_string += "original control points:" offset = 1 for i in range(0, self.n_control_points): - output_string += offset * ' ' + str( - self.original_control_points[i][0]) + ' ' + str( - self.original_control_points[i][1]) + ' ' + str( - self.original_control_points[i][2]) + '\n' + output_string += ( + offset * " " + + str(self.original_control_points[i][0]) + + " " + + str(self.original_control_points[i][1]) + + " " + + str(self.original_control_points[i][2]) + + "\n" + ) offset = 25 - output_string += '\n# deformed control points collects the coordinates' - output_string += ' of the interpolation control points after the' - output_string += ' deformation.\n' + output_string += "\n# deformed control points collects the coordinates" + output_string += " of the interpolation control points after the" + output_string += " deformation.\n" - output_string += 'deformed control points:' + output_string += "deformed control points:" offset = 1 for i in range(0, self.n_control_points): - output_string += offset * ' ' + str( - self.deformed_control_points[i][0]) + ' ' + str( - self.deformed_control_points[i][1]) + ' ' + str( - self.deformed_control_points[i][2]) + '\n' + output_string += ( + offset * " " + + str(self.deformed_control_points[i][0]) + + " " + + str(self.deformed_control_points[i][1]) + + " " + + str(self.deformed_control_points[i][2]) + + "\n" + ) offset = 25 - with open(filename, 'w') as f: + with open(filename, "w") as f: f.write(output_string) def __str__(self): @@ -396,14 +441,14 @@ def __str__(self): This method prints all the RBF parameters on the screen. Its purpose is for debugging. """ - string = '' - string += 'basis function = {}\n'.format(self.basis) - string += 'radius = {}\n'.format(self.radius) - string += 'extra_parameter = {}\n'.format(self.extra) - string += '\noriginal control points =\n' - string += '{}\n'.format(self.original_control_points) - string += '\ndeformed control points =\n' - string += '{}\n'.format(self.deformed_control_points) + string = "" + string += "basis function = {}\n".format(self.basis) + string += "radius = {}\n".format(self.radius) + string += "extra_parameter = {}\n".format(self.extra) + string += "\noriginal control points =\n" + string += "{}\n".format(self.original_control_points) + string += "\ndeformed control points =\n" + string += "{}\n".format(self.deformed_control_points) return string def plot_points(self, filename=None): @@ -415,27 +460,34 @@ def plot_points(self, filename=None): on the specified `filename`. Default is None. """ fig = plt.figure(1) - axes = fig.add_subplot(111, projection='3d') - orig = axes.scatter(self.original_control_points[:, 0], - self.original_control_points[:, 1], - self.original_control_points[:, 2], - c='blue', - marker='o') - defor = axes.scatter(self.deformed_control_points[:, 0], - self.deformed_control_points[:, 1], - self.deformed_control_points[:, 2], - c='red', - marker='x') - - axes.set_xlabel('X axis') - axes.set_ylabel('Y axis') - axes.set_zlabel('Z axis') - - plt.legend((orig, defor), ('Original', 'Deformed'), - scatterpoints=1, - loc='lower left', - ncol=2, - fontsize=10) + axes = fig.add_subplot(111, projection="3d") + orig = axes.scatter( + self.original_control_points[:, 0], + self.original_control_points[:, 1], + self.original_control_points[:, 2], + c="blue", + marker="o", + ) + defor = axes.scatter( + self.deformed_control_points[:, 0], + self.deformed_control_points[:, 1], + self.deformed_control_points[:, 2], + c="red", + marker="x", + ) + + axes.set_xlabel("X axis") + axes.set_ylabel("Y axis") + axes.set_zlabel("Z axis") + + plt.legend( + (orig, defor), + ("Original", "Deformed"), + scatterpoints=1, + loc="lower left", + ncol=2, + fontsize=10, + ) # Show the plot to the screen if filename is None: @@ -448,8 +500,9 @@ def compute_weights(self): This method compute the weights according to the `original_control_points` and `deformed_control_points` arrays. """ - self.weights = self._get_weights(self.original_control_points, - self.deformed_control_points) + self.weights = self._get_weights( + self.original_control_points, self.deformed_control_points + ) def __call__(self, src_pts): """ @@ -459,13 +512,16 @@ def __call__(self, src_pts): src = np.asarray(src_pts, dtype=self._dtype) self.compute_weights() - H = np.zeros((src.shape[0], self.n_control_points + 3 + 1), - dtype=self._dtype) + H = np.zeros( + (src.shape[0], self.n_control_points + 3 + 1), dtype=self._dtype + ) dists = cdist(src, self.original_control_points).astype(self._dtype) basis_block = self.basis(dists, self.radius, **self.extra) - H[:, :self.n_control_points] = np.asarray(basis_block, dtype=self._dtype) + H[:, : self.n_control_points] = np.asarray( + basis_block, dtype=self._dtype + ) H[:, self.n_control_points] = self._dtype(1.0) H[:, -3:] = src diff --git a/pygem/rbf_factory.py b/pygem/rbf_factory.py index 8441c0d8..eaad6717 100644 --- a/pygem/rbf_factory.py +++ b/pygem/rbf_factory.py @@ -1,10 +1,11 @@ """ Factory class for radial basis functions """ + import numpy as np -class classproperty(): +class classproperty: def __init__(self, f): self.f = f @@ -12,18 +13,19 @@ def __get__(self, obj, owner): return self.f(owner) -class RBFFactory(): +class RBFFactory: """ Factory class that spawns the radial basis functions. :Example: - + >>> from pygem import RBFFactory >>> import numpy as np >>> x = np.linspace(0, 1) >>> for fname in RBFFactory.bases: >>> y = RBFFactory(fname)(x) """ + @staticmethod def gaussian_spline(X, r=1): """ @@ -89,11 +91,11 @@ def thin_plate_spline(X, r=1, k=2): With k=2 the function is "radius free", that means independent of radius value. :param numpy.ndarray X: the norm x in the formula above. - :param float r: the parameter r in the formula above. + :param float r: the parameter r in the formula above. :param float k: the parameter k in the formula above. - + :return: result: the result of the formula above. - :rtype: float + :rtype: float """ arg = X / r result = np.power(arg, k) @@ -106,7 +108,7 @@ def beckert_wendland_c2_basis(X, r=1): It implements the following formula: .. math:: - \\varphi(\\boldsymbol{x}) = + \\varphi(\\boldsymbol{x}) = \\left( 1 - \\frac{\\boldsymbol{x}}{r}\\right)^4 + \\left( 4 \\frac{ \\boldsymbol{x} }{r} + 1 \\right) @@ -158,9 +160,11 @@ def polyharmonic_spline(X, r=1, k=2): return np.power(r_sc, k) # k even - result = np.where(r_sc < 1, - np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc)), - np.power(r_sc, k) * np.log(r_sc)) + result = np.where( + r_sc < 1, + np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc)), + np.power(r_sc, k) * np.log(r_sc), + ) return result ############################################################################ @@ -171,14 +175,12 @@ def polyharmonic_spline(X, r=1, k=2): ## ## ############################################################################ __bases = { - 'gaussian_spline': gaussian_spline.__func__, - 'multi_quadratic_biharmonic_spline': - multi_quadratic_biharmonic_spline.__func__, - 'inv_multi_quadratic_biharmonic_spline': - inv_multi_quadratic_biharmonic_spline.__func__, - 'thin_plate_spline': thin_plate_spline.__func__, - 'beckert_wendland_c2_basis': beckert_wendland_c2_basis.__func__, - 'polyharmonic_spline': polyharmonic_spline.__func__ + "gaussian_spline": gaussian_spline.__func__, + "multi_quadratic_biharmonic_spline": multi_quadratic_biharmonic_spline.__func__, + "inv_multi_quadratic_biharmonic_spline": inv_multi_quadratic_biharmonic_spline.__func__, + "thin_plate_spline": thin_plate_spline.__func__, + "beckert_wendland_c2_basis": beckert_wendland_c2_basis.__func__, + "polyharmonic_spline": polyharmonic_spline.__func__, } def __new__(self, fname): @@ -190,7 +192,8 @@ def __new__(self, fname): raise NameError( """The name of the basis function in the parameters file is not correct or not implemented. Check the documentation for - all the available functions.""") + all the available functions.""" + ) @classproperty def bases(self): diff --git a/pygem/stlhandler.py b/pygem/stlhandler.py index e2219183..02441604 100644 --- a/pygem/stlhandler.py +++ b/pygem/stlhandler.py @@ -3,8 +3,9 @@ .. warning:: This module will be deprecated in next releases. Follow updates on - https://github.com/mathLab for news about file handling. + https://github.com/mathLab for news about file handling. """ + import numpy as np from mpl_toolkits import mplot3d import matplotlib.pyplot as plt @@ -12,7 +13,10 @@ import pygem.filehandler as fh import vtk import warnings -warnings.warn("This module will be deprecated in next releases", DeprecationWarning) + +warnings.warn( + "This module will be deprecated in next releases", DeprecationWarning +) class StlHandler(fh.FileHandler): @@ -27,7 +31,7 @@ class StlHandler(fh.FileHandler): def __init__(self): super(StlHandler, self).__init__() - self.extensions = ['.stl'] + self.extensions = [".stl"] def parse(self, filename): """ @@ -35,7 +39,7 @@ def parse(self, filename): coordinates. :param string filename: name of the input file. - + :return: mesh_points: it is a `n_points`-by-3 matrix containing the coordinates of the points of the mesh :rtype: numpy.ndarray @@ -54,8 +58,9 @@ def parse(self, filename): mesh_points = np.zeros([n_points, 3]) for i in range(n_points): - mesh_points[i][0], mesh_points[i][1], mesh_points[i][ - 2] = data.GetPoint(i) + mesh_points[i][0], mesh_points[i][1], mesh_points[i][2] = ( + data.GetPoint(i) + ) return mesh_points @@ -109,7 +114,7 @@ def plot(self, plot_file=None, save_fig=False): :param string plot_file: the stl filename you want to plot. :param bool save_fig: a flag to save the figure in png or not. If True the plot is not shown. The default value is False. - + :return: figure: matlplotlib structure for the figure of the chosen geometry :rtype: matplotlib.pyplot.figure @@ -136,35 +141,40 @@ def plot(self, plot_file=None, save_fig=False): for j in range(0, 3): cell = data.GetCell(i).GetPointId(j) vtx[i][j][0], vtx[i][j][1], vtx[i][j][2] = points.GetPoint( - int(cell)) + int(cell) + ) tri = a3.art3d.Poly3DCollection([vtx[i]]) - tri.set_color('b') - tri.set_edgecolor('k') + tri.set_color("b") + tri.set_edgecolor("k") axes.add_collection3d(tri) ## Get the limits of the axis and center the geometry max_dim = np.array( - [np.max(vtx[:, :, 0]), - np.max(vtx[:, :, 1]), - np.max(vtx[:, :, 2])]) + [np.max(vtx[:, :, 0]), np.max(vtx[:, :, 1]), np.max(vtx[:, :, 2])] + ) min_dim = np.array( - [np.min(vtx[:, :, 0]), - np.min(vtx[:, :, 1]), - np.min(vtx[:, :, 2])]) + [np.min(vtx[:, :, 0]), np.min(vtx[:, :, 1]), np.min(vtx[:, :, 2])] + ) max_lenght = np.max(max_dim - min_dim) - axes.set_xlim(-.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, - .6 * max_lenght + (max_dim[0] + min_dim[0]) / 2) - axes.set_ylim(-.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, - .6 * max_lenght + (max_dim[1] + min_dim[1]) / 2) - axes.set_zlim(-.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, - .6 * max_lenght + (max_dim[2] + min_dim[2]) / 2) + axes.set_xlim( + -0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + 0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + ) + axes.set_ylim( + -0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + 0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + ) + axes.set_zlim( + -0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + 0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + ) # Show the plot to the screen if not save_fig: plt.show() else: - figure.savefig(plot_file.split('.')[0] + '.png') + figure.savefig(plot_file.split(".")[0] + ".png") return figure diff --git a/pygem/unvhandler.py b/pygem/unvhandler.py index bd772b5d..425f5d94 100644 --- a/pygem/unvhandler.py +++ b/pygem/unvhandler.py @@ -3,12 +3,16 @@ .. warning:: This module will be deprecated in next releases. Follow updates on - https://github.com/mathLab for news about file handling. + https://github.com/mathLab for news about file handling. """ + import numpy as np import pygem.filehandler as fh import warnings -warnings.warn("This module will be deprecated in next releases", DeprecationWarning) + +warnings.warn( + "This module will be deprecated in next releases", DeprecationWarning +) class UnvHandler(fh.FileHandler): @@ -23,7 +27,7 @@ class UnvHandler(fh.FileHandler): def __init__(self): super(UnvHandler, self).__init__() - self.extensions = ['.unv'] + self.extensions = [".unv"] def parse(self, filename): """ @@ -44,15 +48,15 @@ def parse(self, filename): index = -9 mesh_points = [] - with open(self.infile, 'r') as input_file: + with open(self.infile, "r") as input_file: for num, line in enumerate(input_file): - if line.startswith(' 2411'): + if line.startswith(" 2411"): index = num if num == index + 2: - if line.startswith(' -1'): + if line.startswith(" -1"): break else: - line = line.replace('D', 'E') + line = line.replace("D", "E") l = [] for t in line.split(): try: @@ -83,20 +87,22 @@ def write(self, mesh_points, filename): index = -9 i = 0 - with open(self.outfile, 'w') as output_file: - with open(self.infile, 'r') as input_file: + with open(self.outfile, "w") as output_file: + with open(self.infile, "r") as input_file: for num, line in enumerate(input_file): - if line.startswith(' 2411'): + if line.startswith(" 2411"): index = num if num == index + 2: - if line.startswith(' -1'): + if line.startswith(" -1"): index = -9 output_file.write(line) else: for j in range(0, 3): - output_file.write(3 * ' ' + '{:.16E}'.format( - mesh_points[i][j])) - output_file.write('\n') + output_file.write( + 3 * " " + + "{:.16E}".format(mesh_points[i][j]) + ) + output_file.write("\n") i += 1 index = num else: diff --git a/pygem/utils.py b/pygem/utils.py index 6bfb066b..8ece7a51 100644 --- a/pygem/utils.py +++ b/pygem/utils.py @@ -2,6 +2,7 @@ Utilities for the affine transformations of the bounding box of the Free Form Deformation. """ + import math from functools import reduce import numpy as np @@ -44,17 +45,20 @@ def angles2matrix(rot_z=0, rot_y=0, rot_x=0): cos = math.cos(rot_z) sin = math.sin(rot_z) rot_matrix.append( - np.array([cos, -sin, 0, sin, cos, 0, 0, 0, 1]).reshape((3, 3))) + np.array([cos, -sin, 0, sin, cos, 0, 0, 0, 1]).reshape((3, 3)) + ) if rot_y: cos = math.cos(rot_y) sin = math.sin(rot_y) rot_matrix.append( - np.array([cos, 0, sin, 0, 1, 0, -sin, 0, cos]).reshape((3, 3))) + np.array([cos, 0, sin, 0, 1, 0, -sin, 0, cos]).reshape((3, 3)) + ) if rot_x: cos = math.cos(rot_x) sin = math.sin(rot_x) rot_matrix.append( - np.array([1, 0, 0, 0, cos, -sin, 0, sin, cos]).reshape((3, 3))) + np.array([1, 0, 0, 0, cos, -sin, 0, sin, cos]).reshape((3, 3)) + ) if rot_matrix: return reduce(np.dot, rot_matrix[::-1]) return np.eye(3) @@ -90,14 +94,15 @@ def fit_affine_transformation(points_start, points_end): dim = len(points_start[0]) if len(points_start) < dim: raise RuntimeError( - "Too few starting points => under-determined system.") + "Too few starting points => under-determined system." + ) def pad_column_ones(x): - """ Add right column of 1.0 to the given 2D numpy array """ + """Add right column of 1.0 to the given 2D numpy array""" return np.hstack([x, np.ones((x.shape[0], 1))]) def unpad_column(x): - """ Remove last column to the given 2D numpy array """ + """Remove last column to the given 2D numpy array""" return x[:, :-1] def transform(src): @@ -109,13 +114,11 @@ def transform(src): A, res, rank, _ = np.linalg.lstsq(X, Y, rcond=None) # TODO add check condition number - #if np.linalg.cond(A) >= 1 / sys.float_info.epsilon: + # if np.linalg.cond(A) >= 1 / sys.float_info.epsilon: # raise RuntimeError( # "Error: singular matrix. Points are probably coplanar.") return unpad_column( - np.dot( - pad_column_ones(np.atleast_2d(src)), - A) - ).reshape(shape) + np.dot(pad_column_ones(np.atleast_2d(src)), A) + ).reshape(shape) return transform diff --git a/pygem/vffd.py b/pygem/vffd.py index c464ec25..660bdbd9 100644 --- a/pygem/vffd.py +++ b/pygem/vffd.py @@ -1,12 +1,14 @@ from pygem.cffd import CFFD import numpy as np + + class VFFD(CFFD): - ''' + """ Class that handles the Volumetric Free Form Deformation on the mesh points. - + :param list n_control_points: number of control points in the x, y, and z direction. Default is [2, 2, 2]. - :param string mode: it can be ``affine`` or ``triaffine``. The first option is for the F that are affine in all the coordinates of the points. + :param string mode: it can be ``affine`` or ``triaffine``. The first option is for the F that are affine in all the coordinates of the points. The second one is for functions that are F in the coordinates of the points. The first option implies the second, but is optimal for that class of functions. :cvar numpy.ndarray box_length: dimension of the FFD bounding box, in the x, y and z direction (local coordinate system). @@ -22,11 +24,11 @@ class VFFD(CFFD): z, normalized with the box length z. :cvar callable fun: it defines the F of the constraint F(x)=c. Default is the constant 1 function. :cvar numpy.ndarray fixval: it defines the c of the constraint F(x)=c. Default is 1. - :cvar numpy.ndarray ffd_mask: a boolean tensor that tells to the class - which control points can be moved, and in what direction, to enforce the constraint. + :cvar numpy.ndarray ffd_mask: a boolean tensor that tells to the class + which control points can be moved, and in what direction, to enforce the constraint. The tensor has shape (n_x,n_y,n_z,3), where the last dimension indicates movement on x,y,z respectively. Default is all true. - :cvar numpy.ndarray fun_mask: a boolean tensor that tells to the class + :cvar numpy.ndarray fun_mask: a boolean tensor that tells to the class on which axis which constraint depends on. The tensor has shape (n_cons,3), where the last dimension indicates dependency on on x,y,z respectively. Default is all true. It used only in the triaffine mode. @@ -45,7 +47,8 @@ class VFFD(CFFD): >>> new_mesh_points = vffd(original_mesh_points) >>> assert np.isclose(np.linalg.norm(vffd.fun(new_mesh_points) - b), np.array([0.]), atol=1e-07) - ''' + """ + def __init__(self, triangles, fixval, n_control_points=None, ffd_mask=None): super().__init__(fixval, None, n_control_points, ffd_mask, None) @@ -63,5 +66,3 @@ def _volume(x, triangles): x = x.reshape(-1, 3) mesh = x[triangles] return np.array([np.sum(np.linalg.det(mesh))]) - - diff --git a/pygem/vtkhandler.py b/pygem/vtkhandler.py index 3fcb61ab..41a7ba8b 100644 --- a/pygem/vtkhandler.py +++ b/pygem/vtkhandler.py @@ -3,15 +3,19 @@ .. warning:: This module will be deprecated in next releases. Follow updates on - https://github.com/mathLab for news about file handling. + https://github.com/mathLab for news about file handling. """ + import numpy as np import matplotlib.pyplot as plt import mpl_toolkits.mplot3d as a3 import vtk import pygem.filehandler as fh import warnings -warnings.warn("This module will be deprecated in next releases", DeprecationWarning) + +warnings.warn( + "This module will be deprecated in next releases", DeprecationWarning +) class VtkHandler(fh.FileHandler): @@ -26,7 +30,7 @@ class VtkHandler(fh.FileHandler): def __init__(self): super(VtkHandler, self).__init__() - self.extensions = ['.vtk'] + self.extensions = [".vtk"] def parse(self, filename): """ @@ -59,8 +63,9 @@ def parse(self, filename): mesh_points = np.zeros([n_points, 3]) for i in range(n_points): - mesh_points[i][0], mesh_points[i][1], mesh_points[i][ - 2] = data.GetPoint(i) + mesh_points[i][0], mesh_points[i][1], mesh_points[i][2] = ( + data.GetPoint(i) + ) return mesh_points @@ -108,7 +113,7 @@ def plot(self, plot_file=None, save_fig=False): :param string plot_file: the vtk filename you want to plot. :param bool save_fig: a flag to save the figure in png or not. If True the plot is not shown. - + :return: figure: matlplotlib structure for the figure of the chosen geometry :rtype: matplotlib.pyplot.figure @@ -135,35 +140,40 @@ def plot(self, plot_file=None, save_fig=False): for j in range(0, 3): cell = data.GetCell(i).GetPointId(j) vtx[i][j][0], vtx[i][j][1], vtx[i][j][2] = points.GetPoint( - int(cell)) + int(cell) + ) tri = a3.art3d.Poly3DCollection([vtx[i]]) - tri.set_color('b') - tri.set_edgecolor('k') + tri.set_color("b") + tri.set_edgecolor("k") axes.add_collection3d(tri) ## Get the limits of the axis and center the geometry max_dim = np.array( - [np.max(vtx[:, :, 0]), - np.max(vtx[:, :, 1]), - np.max(vtx[:, :, 2])]) + [np.max(vtx[:, :, 0]), np.max(vtx[:, :, 1]), np.max(vtx[:, :, 2])] + ) min_dim = np.array( - [np.min(vtx[:, :, 0]), - np.min(vtx[:, :, 1]), - np.min(vtx[:, :, 2])]) + [np.min(vtx[:, :, 0]), np.min(vtx[:, :, 1]), np.min(vtx[:, :, 2])] + ) max_lenght = np.max(max_dim - min_dim) - axes.set_xlim(-.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, - .6 * max_lenght + (max_dim[0] + min_dim[0]) / 2) - axes.set_ylim(-.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, - .6 * max_lenght + (max_dim[1] + min_dim[1]) / 2) - axes.set_zlim(-.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, - .6 * max_lenght + (max_dim[2] + min_dim[2]) / 2) + axes.set_xlim( + -0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + 0.6 * max_lenght + (max_dim[0] + min_dim[0]) / 2, + ) + axes.set_ylim( + -0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + 0.6 * max_lenght + (max_dim[1] + min_dim[1]) / 2, + ) + axes.set_zlim( + -0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + 0.6 * max_lenght + (max_dim[2] + min_dim[2]) / 2, + ) # Show the plot to the screen if not save_fig: plt.show() else: - figure.savefig(plot_file.split('.')[0] + '.png') + figure.savefig(plot_file.split(".")[0] + ".png") return figure