From 46715706153cc81289fa9cfbed11b046394e8955 Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 09:41:04 +0200 Subject: [PATCH 01/40] First working version + some debug --- .../mesh/processing/rotateAndTranslate.py | 111 ++++++++++++++++++ geos-mesh/tests/test_rotateAndTranslate.py | 98 ++++++++++++++++ 2 files changed, 209 insertions(+) create mode 100644 geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py create mode 100644 geos-mesh/tests/test_rotateAndTranslate.py diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py new file mode 100644 index 000000000..9a236beb1 --- /dev/null +++ b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py @@ -0,0 +1,111 @@ +from vtkmodules.numpy_interface import dataset_adapter as dsa +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from vtkmodules.vtkCommonCore import vtkPoints +import numpy as np +import logging + +def transform_mesh( input : vtkUnstructuredGrid, logger : logging.Logger | None = None ) -> vtkUnstructuredGrid: + """Apply a transformation to a mesh + + Args: + input (vtk.vtkUnstructuredGrid): input mesh + logger (logging.logger | None, optional): logger instance. Defaults to None. + + Returns: + vtk.vtkUnstructuredGrid: transformed mesh + """ + # Initialize the logger if it is empty + if not logger: + logging.basicConfig( level=logging.WARNING ) + logger = logging.getLogger( __name__ ) + + translation, rotation, pts = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy(input.GetPoints().GetData()), logger ) + + logger.info(f"Translation {translation}") + logger.info(f"Rotation {rotation}") + + output = vtkUnstructuredGrid() + (vpts:= vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) + output.SetPoints(vpts) + output.SetCells(input.GetCellTypesArray(), input.GetCellLocationsArray(), input.GetCells()) + + return output + +def local_frame(pts): + # find 3 orthogonal vectors + # we assume points are on a box + # first vector is along x axis + ori, _, _ = origin_bounding_box(pts) + assert ((ori == np.asarray([0, 0, 0])).all()) + u = pts[1] + v = u + for pt in pts[2:]:# As we assume points are on a box there should be one orthogonal to u + if (np.abs(np.dot(u, pt)) < 0.0001): + v = pt + break + + return (u, v, np.cross(u, v)) + + # utils +def origin_bounding_box( pts : np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Find the bounding box of a set of points + Args: + pts (np.ndarray): points to find the bounding box of + Returns: + tuple[np.ndarray, np.ndarray, np.ndarray]: origin, min, max of the bounding box + """ + + pt0 = pts[np.argmin(pts, axis=0), :] + pt1 = pts[np.argmax(pts, axis=0), :] + + return (pts[0], pt0, pt1) + +def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Recenter and rotate points to align with principal axes + + Args: + pts (np.ndarray): points to recenter and rotate + logger (logging.logger, optional): logger instance. + """ + + # find bounding box + org, vmin, vmax = origin_bounding_box(pts) + logger.info(f"Bounding box is {org}, {vmin}, {vmax}") + + # Transformation + translation = org + pts -= translation + + u, v, w = local_frame(pts) + logger.info(f"Local frame u {u}, v {v}, w {w}") + + # find rotation R = U sig V + rotation = np.matrix([u / np.linalg.norm(u), v / np.linalg.norm(v), w / np.linalg.norm(w)]).transpose() + logger.info(f"R {rotation}") + logger.info(f"theta {np.acos(.5 * (np.trace(rotation) - 1)) * 180 / np.pi}") + + pts = (rotation.transpose() * pts.transpose()).transpose() + pts[np.abs(pts)<1e-15] = 0. # clipping points too close to zero + logger.info(f"Un-rotated pts : {pts}") + + return translation, rotation, pts + + +# #TO DO : remove +# if __name__ == "__main__": + +# reader = vtk.vtkXMLUnstructuredGridReader() +# reader.SetFileName("/home/jfranc/codes/python/geosPythonPackages/geos-mesh/src/geos/mesh/processing/input.vtu") +# reader.Update() + +# # Transformation +# output = transform_mesh(reader.GetOutput()) + +# output.SetCells(VTK_HEXAHEDRON, reader.GetOutput().GetCells()) +# writer = vtk.vtkXMLUnstructuredGridWriter() +# writer.SetFileName("./output.vtu") +# writer.SetInputData(output) +# writer.Update() +# writer.Write() diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py new file mode 100644 index 000000000..c913015ce --- /dev/null +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -0,0 +1,98 @@ +import pytest +from dataclasses import dataclass +from typing import Generator +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, VTK_HEXAHEDRON +from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridWriter # for writing the output tmp +from vtkmodules.numpy_interface import dataset_adapter as dsa +import numpy as np +from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON + +from geos.mesh.processing.rotateAndTranslate import transform_mesh + +@dataclass( frozen=True ) +class Expected: + mesh: vtkUnstructuredGrid + +def __gen_box(Lx:int, Ly:int, Lz:int, nx:int, ny:int, nz:int)-> tuple[np.ndarray, np.ndarray]: + np.random.seed(1) + off = np.random.randn(1, 3) + pts = [] + x,y,z = np.meshgrid(np.linspace(0,Lx,nx),np.linspace(0,Ly,ny),np.linspace(0,Lz,ny)) + for i in range(x.shape[0]): + for j in range(x.shape[1]): + for k in range(x.shape[2]): + pts.append(np.asarray([x[i,j,k], y[i,j,k], z[i,j,k]])) + + return (np.asarray(pts), off) + +def __rotate_box(angles : np.ndarray, pts:np.ndarray) -> np.ndarray: + a = angles[0] + RX = np.matrix([[1., 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]]) + a = angles[1] + RY = np.matrix([[np.cos(a), 0, np.sin(a)], [0, 1, 0], [-np.sin(a), 0, np.cos(a)]]) + a = angles[2] + RZ = np.matrix([[np.cos(a), -np.sin(a), 0], [np.sin(a), np.cos(a), 0], [0, 0, 1]]) + + return np.asarray((RZ * RY * RX * pts.transpose()).transpose()) + +def __build_test_mesh() -> Generator[ Expected, None, None ]: + # generate random points in a box Lx, Ly, Lz + np.random.seed(1) + Lx, Ly, Lz = (2, 5, 8) # box size + nx, ny, nz = (10, 10, 10) # number of points in each direction + pts, off = __gen_box(Lx, Ly, Lz, nx, ny, nz) + + print(f"Offseting of {off}") + print(f"Original pts : {pts}") + angles = -2*np.pi + np.random.randn(1, 3)*np.pi # random angles in rad + print(f"angles {angles[0]}") + pts = __rotate_box(angles[0], pts) + print(f"Rotated pts : {pts}") + pts[:, 0] += off[0][0] + pts[:, 1] += off[0][1] + pts[:, 2] += off[0][2] + print(f"Translated pts : {pts}") + + # Creating multiple meshes, each time with a different angles + mesh = vtkUnstructuredGrid() + (vtps:=vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) + mesh.SetPoints( vtps ) + + ids = vtkIdList() + for i in range(nx-1): + for j in range(ny-1): + for k in range(nz-1): + hex = vtkHexahedron() + ids = hex.GetPointIds() + ids.SetId(0,i+j*nx+k*nx*ny) + ids.SetId(1,i+j*nx+k*nx*ny + 1) + ids.SetId(2,i+j*nx+k*nx*ny + nx*ny+1) + ids.SetId(3,i+j*nx+k*nx*ny + nx*ny) + ids.SetId(4,i+j*nx+k*nx*ny + nx) + ids.SetId(5,i+j*nx+k*nx*ny + nx+1) + ids.SetId(6,i+j*nx+k*nx*ny + nx*ny+nx+1) + ids.SetId(7,i+j*nx+k*nx*ny + nx*ny+nx) + mesh.InsertNextCell(VTK_HEXAHEDRON,ids) + + yield Expected( mesh=mesh ) + +@pytest.mark.parametrize( "expected", __build_test_mesh() ) +def test_reorient_polyhedron( expected: Expected ) -> None: + output_mesh = transform_mesh( expected.mesh ) + assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() + assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() + print(f"Bounds {output_mesh.GetBounds()}") + assert output_mesh.GetBounds()[0] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[2] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[4] == pytest.approx(0., abs=1e-10) + #TODO more assert but need more assumptions then + # temp + w = vtkXMLUnstructuredGridWriter() + w.SetFileName("./test_rotateAndTranslate.vtk") + w.SetInputData(output_mesh) + w.Write() + w.SetFileName("./test_rotateAndTranslate_input.vtk") + w.SetInputData(expected.mesh) + w.Write() + +if __name__ == "__main__": + pytest.main([__file__]) \ No newline at end of file From f70a9be747ea64da8234c89722a1f3d1e647eb86 Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 09:56:15 +0200 Subject: [PATCH 02/40] get rid of np.matrix --- geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py | 4 ++-- geos-mesh/tests/test_rotateAndTranslate.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py index 9a236beb1..4b464c21e 100644 --- a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py +++ b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py @@ -82,11 +82,11 @@ def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.n logger.info(f"Local frame u {u}, v {v}, w {w}") # find rotation R = U sig V - rotation = np.matrix([u / np.linalg.norm(u), v / np.linalg.norm(v), w / np.linalg.norm(w)]).transpose() + rotation = np.asarray([u / np.linalg.norm(u), v / np.linalg.norm(v), w / np.linalg.norm(w)],dtype=np.float64).transpose() logger.info(f"R {rotation}") logger.info(f"theta {np.acos(.5 * (np.trace(rotation) - 1)) * 180 / np.pi}") - pts = (rotation.transpose() * pts.transpose()).transpose() + pts = (rotation.transpose() @ pts.transpose()).transpose() pts[np.abs(pts)<1e-15] = 0. # clipping points too close to zero logger.info(f"Un-rotated pts : {pts}") diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index c913015ce..ed986f0f2 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -28,13 +28,13 @@ def __gen_box(Lx:int, Ly:int, Lz:int, nx:int, ny:int, nz:int)-> tuple[np.ndarray def __rotate_box(angles : np.ndarray, pts:np.ndarray) -> np.ndarray: a = angles[0] - RX = np.matrix([[1., 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]]) + RX = np.asarray([[1., 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]]) a = angles[1] - RY = np.matrix([[np.cos(a), 0, np.sin(a)], [0, 1, 0], [-np.sin(a), 0, np.cos(a)]]) + RY = np.asarray([[np.cos(a), 0, np.sin(a)], [0, 1, 0], [-np.sin(a), 0, np.cos(a)]]) a = angles[2] - RZ = np.matrix([[np.cos(a), -np.sin(a), 0], [np.sin(a), np.cos(a), 0], [0, 0, 1]]) + RZ = np.asarray([[np.cos(a), -np.sin(a), 0], [np.sin(a), np.cos(a), 0], [0, 0, 1]]) - return np.asarray((RZ * RY * RX * pts.transpose()).transpose()) + return np.asarray((RZ @ RY @ RX @ pts.transpose()).transpose()) def __build_test_mesh() -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz From faeb6a429c455d514796e419d09f64a882e065ef Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 10:54:22 +0200 Subject: [PATCH 03/40] Using PV Transform filter --- .../mesh/processing/rotateAndTranslate.py | 40 ++++++++++++++----- geos-mesh/tests/test_rotateAndTranslate.py | 38 +++++++++--------- 2 files changed, 49 insertions(+), 29 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py index 4b464c21e..ff343e9d1 100644 --- a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py +++ b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py @@ -1,5 +1,7 @@ from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from vtkmodules.vtkCommonTransforms import vtkTransform +from vtkmodules.vtkFiltersGeneral import vtkTransformFilter from vtkmodules.vtkCommonCore import vtkPoints import numpy as np import logging @@ -19,17 +21,26 @@ def transform_mesh( input : vtkUnstructuredGrid, logger : logging.Logger | None logging.basicConfig( level=logging.WARNING ) logger = logging.getLogger( __name__ ) - translation, rotation, pts = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy(input.GetPoints().GetData()), logger ) + translation, theta, axis = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy(input.GetPoints().GetData()), logger ) - logger.info(f"Translation {translation}") - logger.info(f"Rotation {rotation}") + logger.info(f"Theta {theta}") + logger.info(f"Axis {axis}") - output = vtkUnstructuredGrid() - (vpts:= vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) - output.SetPoints(vpts) - output.SetCells(input.GetCellTypesArray(), input.GetCellLocationsArray(), input.GetCells()) + # output = vtkUnstructuredGrid() + # (vpts:= vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) + # output.SetPoints(vpts) + # output.SetCells(input.GetCellTypesArray(), input.GetCellLocationsArray(), input.GetCells()) - return output + transform = vtkTransform() + transform.PostMultiply() # we want to apply rotation then translation + transform.RotateWXYZ(-theta, axis[0], axis[1], axis[2]) + transform.Translate(-translation[0], -translation[1], -translation[2]) + transformFilter = vtkTransformFilter() + transformFilter.SetTransform(transform) + transformFilter.SetInputData(input) + transformFilter.Update() + + return transformFilter.GetOutput() def local_frame(pts): # find 3 orthogonal vectors @@ -61,7 +72,7 @@ def origin_bounding_box( pts : np.ndarray) -> tuple[np.ndarray, np.ndarray, np.n return (pts[0], pt0, pt1) -def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.ndarray, np.ndarray, np.ndarray]: +def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.ndarray, float, np.ndarray]: """ Recenter and rotate points to align with principal axes @@ -84,13 +95,20 @@ def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.n # find rotation R = U sig V rotation = np.asarray([u / np.linalg.norm(u), v / np.linalg.norm(v), w / np.linalg.norm(w)],dtype=np.float64).transpose() logger.info(f"R {rotation}") - logger.info(f"theta {np.acos(.5 * (np.trace(rotation) - 1)) * 180 / np.pi}") + + theta = np.acos(.5 * (np.trace(rotation) - 1)) * 180 / np.pi + logger.info(f"theta {theta}") + + axis = np.asarray([rotation[2, 1] - rotation[1, 2], rotation[0, 2] - rotation[2, 0], rotation[1, 0] - rotation[0, 1]],dtype=np.float64) + axis /= np.linalg.norm(axis) + logger.info(f"axis {axis}") pts = (rotation.transpose() @ pts.transpose()).transpose() pts[np.abs(pts)<1e-15] = 0. # clipping points too close to zero logger.info(f"Un-rotated pts : {pts}") - return translation, rotation, pts + # return translation, rotation, pts + return translation, theta, axis # #TO DO : remove diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index ed986f0f2..656336baf 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -4,6 +4,7 @@ from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, VTK_HEXAHEDRON from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridWriter # for writing the output tmp +import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON @@ -15,7 +16,8 @@ class Expected: mesh: vtkUnstructuredGrid def __gen_box(Lx:int, Ly:int, Lz:int, nx:int, ny:int, nz:int)-> tuple[np.ndarray, np.ndarray]: - np.random.seed(1) + # np.random.seed(1) # for reproducibility + np.random.default_rng() off = np.random.randn(1, 3) pts = [] x,y,z = np.meshgrid(np.linspace(0,Lx,nx),np.linspace(0,Ly,ny),np.linspace(0,Lz,ny)) @@ -38,21 +40,22 @@ def __rotate_box(angles : np.ndarray, pts:np.ndarray) -> np.ndarray: def __build_test_mesh() -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz - np.random.seed(1) + # np.random.seed(1) # for reproducibility + np.random.default_rng() Lx, Ly, Lz = (2, 5, 8) # box size nx, ny, nz = (10, 10, 10) # number of points in each direction pts, off = __gen_box(Lx, Ly, Lz, nx, ny, nz) - print(f"Offseting of {off}") - print(f"Original pts : {pts}") + logging.warning(f"Offseting of {off}") + logging.warning(f"Original pts : {pts}") angles = -2*np.pi + np.random.randn(1, 3)*np.pi # random angles in rad - print(f"angles {angles[0]}") + logging.warning(f"angles {angles[0]}") pts = __rotate_box(angles[0], pts) - print(f"Rotated pts : {pts}") + logging.info(f"Rotated pts : {pts}") pts[:, 0] += off[0][0] pts[:, 1] += off[0][1] pts[:, 2] += off[0][2] - print(f"Translated pts : {pts}") + logging.info(f"Translated pts : {pts}") # Creating multiple meshes, each time with a different angles mesh = vtkUnstructuredGrid() @@ -78,21 +81,20 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: yield Expected( mesh=mesh ) @pytest.mark.parametrize( "expected", __build_test_mesh() ) -def test_reorient_polyhedron( expected: Expected ) -> None: +def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: output_mesh = transform_mesh( expected.mesh ) assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() - print(f"Bounds {output_mesh.GetBounds()}") assert output_mesh.GetBounds()[0] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[2] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[4] == pytest.approx(0., abs=1e-10) #TODO more assert but need more assumptions then # temp - w = vtkXMLUnstructuredGridWriter() - w.SetFileName("./test_rotateAndTranslate.vtk") - w.SetInputData(output_mesh) - w.Write() - w.SetFileName("./test_rotateAndTranslate_input.vtk") - w.SetInputData(expected.mesh) - w.Write() +# w = vtkXMLUnstructuredGridWriter() +# w.SetFileName("./test_rotateAndTranslate.vtu") +# w.SetInputData(output_mesh) +# w.Write() +# w.SetFileName("./test_rotateAndTranslate_input.vtu") +# w.SetInputData(expected.mesh) +# w.Write() -if __name__ == "__main__": - pytest.main([__file__]) \ No newline at end of file +# if __name__ == "__main__": +# pytest.main([__file__, '-s', '-v','--log-cli-level=WARNING']) \ No newline at end of file From dcea1429549a02a4bd0505072796d7163d138cdf Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 11:27:54 +0200 Subject: [PATCH 04/40] style --- .../mesh/processing/rotateAndTranslate.py | 90 ++++++++------- geos-mesh/tests/test_rotateAndTranslate.py | 106 ++++++++++-------- 2 files changed, 104 insertions(+), 92 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py index ff343e9d1..842b2804c 100644 --- a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py +++ b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py @@ -2,11 +2,11 @@ from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid from vtkmodules.vtkCommonTransforms import vtkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter -from vtkmodules.vtkCommonCore import vtkPoints import numpy as np import logging -def transform_mesh( input : vtkUnstructuredGrid, logger : logging.Logger | None = None ) -> vtkUnstructuredGrid: + +def transform_mesh( input: vtkUnstructuredGrid, logger: logging.Logger | None = None ) -> vtkUnstructuredGrid: """Apply a transformation to a mesh Args: @@ -21,10 +21,11 @@ def transform_mesh( input : vtkUnstructuredGrid, logger : logging.Logger | None logging.basicConfig( level=logging.WARNING ) logger = logging.getLogger( __name__ ) - translation, theta, axis = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy(input.GetPoints().GetData()), logger ) + translation, theta, axis = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( input.GetPoints().GetData() ), + logger ) - logger.info(f"Theta {theta}") - logger.info(f"Axis {axis}") + logger.info( f"Theta {theta}" ) + logger.info( f"Axis {axis}" ) # output = vtkUnstructuredGrid() # (vpts:= vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) @@ -32,80 +33,83 @@ def transform_mesh( input : vtkUnstructuredGrid, logger : logging.Logger | None # output.SetCells(input.GetCellTypesArray(), input.GetCellLocationsArray(), input.GetCells()) transform = vtkTransform() - transform.PostMultiply() # we want to apply rotation then translation - transform.RotateWXYZ(-theta, axis[0], axis[1], axis[2]) - transform.Translate(-translation[0], -translation[1], -translation[2]) + transform.PostMultiply() # we want to apply rotation then translation + transform.RotateWXYZ( -theta, axis[ 0 ], axis[ 1 ], axis[ 2 ] ) + transform.Translate( -translation[ 0 ], -translation[ 1 ], -translation[ 2 ] ) transformFilter = vtkTransformFilter() - transformFilter.SetTransform(transform) - transformFilter.SetInputData(input) + transformFilter.SetTransform( transform ) + transformFilter.SetInputData( input ) transformFilter.Update() return transformFilter.GetOutput() -def local_frame(pts): + +def local_frame( pts ): # find 3 orthogonal vectors # we assume points are on a box # first vector is along x axis - ori, _, _ = origin_bounding_box(pts) - assert ((ori == np.asarray([0, 0, 0])).all()) - u = pts[1] + ori, _, _ = origin_bounding_box( pts ) + assert ( ( ori == np.asarray( [ 0, 0, 0 ] ) ).all() ) + u = pts[ 1 ] v = u - for pt in pts[2:]:# As we assume points are on a box there should be one orthogonal to u - if (np.abs(np.dot(u, pt)) < 0.0001): + for pt in pts[ 2: ]: # As we assume points are on a box there should be one orthogonal to u + if ( np.abs( np.dot( u, pt ) ) < 0.0001 ): v = pt break - return (u, v, np.cross(u, v)) - + return ( u, v, np.cross( u, v ) ) + + # utils -def origin_bounding_box( pts : np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - Find the bounding box of a set of points +def origin_bounding_box( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: + """Find the bounding box of a set of points Args: pts (np.ndarray): points to find the bounding box of Returns: tuple[np.ndarray, np.ndarray, np.ndarray]: origin, min, max of the bounding box """ + pt0 = pts[ np.argmin( pts, axis=0 ), : ] + pt1 = pts[ np.argmax( pts, axis=0 ), : ] - pt0 = pts[np.argmin(pts, axis=0), :] - pt1 = pts[np.argmax(pts, axis=0), :] + return ( pts[ 0 ], pt0, pt1 ) - return (pts[0], pt0, pt1) -def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.ndarray, float, np.ndarray]: - """ - Recenter and rotate points to align with principal axes +def recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np.ndarray, float, np.ndarray ]: + """Recenter and rotate points to align with principal axes Args: pts (np.ndarray): points to recenter and rotate logger (logging.logger, optional): logger instance. """ - # find bounding box - org, vmin, vmax = origin_bounding_box(pts) - logger.info(f"Bounding box is {org}, {vmin}, {vmax}") + org, vmin, vmax = origin_bounding_box( pts ) + logger.info( f"Bounding box is {org}, {vmin}, {vmax}" ) # Transformation translation = org pts -= translation - u, v, w = local_frame(pts) - logger.info(f"Local frame u {u}, v {v}, w {w}") + u, v, w = local_frame( pts ) + logger.info( f"Local frame u {u}, v {v}, w {w}" ) # find rotation R = U sig V - rotation = np.asarray([u / np.linalg.norm(u), v / np.linalg.norm(v), w / np.linalg.norm(w)],dtype=np.float64).transpose() - logger.info(f"R {rotation}") + rotation = np.asarray( [ u / np.linalg.norm( u ), v / np.linalg.norm( v ), w / np.linalg.norm( w ) ], + dtype=np.float64 ).transpose() + logger.info( f"R {rotation}" ) - theta = np.acos(.5 * (np.trace(rotation) - 1)) * 180 / np.pi - logger.info(f"theta {theta}") + theta = np.acos( .5 * ( np.trace( rotation ) - 1 ) ) * 180 / np.pi + logger.info( f"theta {theta}" ) - axis = np.asarray([rotation[2, 1] - rotation[1, 2], rotation[0, 2] - rotation[2, 0], rotation[1, 0] - rotation[0, 1]],dtype=np.float64) - axis /= np.linalg.norm(axis) - logger.info(f"axis {axis}") + axis = np.asarray( [ + rotation[ 2, 1 ] - rotation[ 1, 2 ], rotation[ 0, 2 ] - rotation[ 2, 0 ], rotation[ 1, 0 ] - rotation[ 0, 1 ] + ], + dtype=np.float64 ) + axis /= np.linalg.norm( axis ) + logger.info( f"axis {axis}" ) - pts = (rotation.transpose() @ pts.transpose()).transpose() - pts[np.abs(pts)<1e-15] = 0. # clipping points too close to zero - logger.info(f"Un-rotated pts : {pts}") + pts = ( rotation.transpose() @ pts.transpose() ).transpose() + pts[ np.abs( pts ) < 1e-15 ] = 0. # clipping points too close to zero + logger.info( f"Un-rotated pts : {pts}" ) # return translation, rotation, pts return translation, theta, axis @@ -113,7 +117,7 @@ def recenter_and_rotate(pts : np.ndarray, logger : logging.Logger) -> tuple[np.n # #TO DO : remove # if __name__ == "__main__": - + # reader = vtk.vtkXMLUnstructuredGridReader() # reader.SetFileName("/home/jfranc/codes/python/geosPythonPackages/geos-mesh/src/geos/mesh/processing/input.vtu") # reader.Update() diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index 656336baf..5db8b2484 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -3,91 +3,99 @@ from typing import Generator from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, VTK_HEXAHEDRON -from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridWriter # for writing the output tmp -import logging # for debug +import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON from geos.mesh.processing.rotateAndTranslate import transform_mesh + @dataclass( frozen=True ) class Expected: mesh: vtkUnstructuredGrid -def __gen_box(Lx:int, Ly:int, Lz:int, nx:int, ny:int, nz:int)-> tuple[np.ndarray, np.ndarray]: + +def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int ) -> tuple[ np.ndarray, np.ndarray ]: # np.random.seed(1) # for reproducibility np.random.default_rng() - off = np.random.randn(1, 3) + off = np.random.randn( 1, 3 ) pts = [] - x,y,z = np.meshgrid(np.linspace(0,Lx,nx),np.linspace(0,Ly,ny),np.linspace(0,Lz,ny)) - for i in range(x.shape[0]): - for j in range(x.shape[1]): - for k in range(x.shape[2]): - pts.append(np.asarray([x[i,j,k], y[i,j,k], z[i,j,k]])) + x, y, z = np.meshgrid( np.linspace( 0, Lx, nx ), np.linspace( 0, Ly, ny ), np.linspace( 0, Lz, ny ) ) + for i in range( x.shape[ 0 ] ): + for j in range( x.shape[ 1 ] ): + for k in range( x.shape[ 2 ] ): + pts.append( np.asarray( [ x[ i, j, k ], y[ i, j, k ], z[ i, j, k ] ] ) ) + + return ( np.asarray( pts ), off ) - return (np.asarray(pts), off) -def __rotate_box(angles : np.ndarray, pts:np.ndarray) -> np.ndarray: - a = angles[0] - RX = np.asarray([[1., 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]]) - a = angles[1] - RY = np.asarray([[np.cos(a), 0, np.sin(a)], [0, 1, 0], [-np.sin(a), 0, np.cos(a)]]) - a = angles[2] - RZ = np.asarray([[np.cos(a), -np.sin(a), 0], [np.sin(a), np.cos(a), 0], [0, 0, 1]]) +def __rotate_box( angles: np.ndarray, pts: np.ndarray ) -> np.ndarray: + a = angles[ 0 ] + RX = np.asarray( [ [ 1., 0, 0 ], [ 0, np.cos( a ), -np.sin( a ) ], [ 0, np.sin( a ), np.cos( a ) ] ] ) + a = angles[ 1 ] + RY = np.asarray( [ [ np.cos( a ), 0, np.sin( a ) ], [ 0, 1, 0 ], [ -np.sin( a ), 0, np.cos( a ) ] ] ) + a = angles[ 2 ] + RZ = np.asarray( [ [ np.cos( a ), -np.sin( a ), 0 ], [ np.sin( a ), np.cos( a ), 0 ], [ 0, 0, 1 ] ] ) + + return np.asarray( ( RZ @ RY @ RX @ pts.transpose() ).transpose() ) - return np.asarray((RZ @ RY @ RX @ pts.transpose()).transpose()) def __build_test_mesh() -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz # np.random.seed(1) # for reproducibility np.random.default_rng() - Lx, Ly, Lz = (2, 5, 8) # box size - nx, ny, nz = (10, 10, 10) # number of points in each direction - pts, off = __gen_box(Lx, Ly, Lz, nx, ny, nz) - - logging.warning(f"Offseting of {off}") - logging.warning(f"Original pts : {pts}") - angles = -2*np.pi + np.random.randn(1, 3)*np.pi # random angles in rad - logging.warning(f"angles {angles[0]}") - pts = __rotate_box(angles[0], pts) - logging.info(f"Rotated pts : {pts}") - pts[:, 0] += off[0][0] - pts[:, 1] += off[0][1] - pts[:, 2] += off[0][2] - logging.info(f"Translated pts : {pts}") + Lx, Ly, Lz = ( 2, 5, 8 ) # box size + nx, ny, nz = ( 10, 10, 10 ) # number of points in each direction + pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz ) + + logging.warning( f"Offseting of {off}" ) + logging.warning( f"Original pts : {pts}" ) + angles = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad + logging.warning( f"angles {angles[0]}" ) + pts = __rotate_box( angles[ 0 ], pts ) + logging.info( f"Rotated pts : {pts}" ) + pts[ :, 0 ] += off[ 0 ][ 0 ] + pts[ :, 1 ] += off[ 0 ][ 1 ] + pts[ :, 2 ] += off[ 0 ][ 2 ] + logging.info( f"Translated pts : {pts}" ) # Creating multiple meshes, each time with a different angles mesh = vtkUnstructuredGrid() - (vtps:=vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) + ( vtps := vtkPoints() ).SetData( dsa.numpy_support.numpy_to_vtk( pts ) ) mesh.SetPoints( vtps ) ids = vtkIdList() - for i in range(nx-1): - for j in range(ny-1): - for k in range(nz-1): + for i in range( nx - 1 ): + for j in range( ny - 1 ): + for k in range( nz - 1 ): hex = vtkHexahedron() ids = hex.GetPointIds() - ids.SetId(0,i+j*nx+k*nx*ny) - ids.SetId(1,i+j*nx+k*nx*ny + 1) - ids.SetId(2,i+j*nx+k*nx*ny + nx*ny+1) - ids.SetId(3,i+j*nx+k*nx*ny + nx*ny) - ids.SetId(4,i+j*nx+k*nx*ny + nx) - ids.SetId(5,i+j*nx+k*nx*ny + nx+1) - ids.SetId(6,i+j*nx+k*nx*ny + nx*ny+nx+1) - ids.SetId(7,i+j*nx+k*nx*ny + nx*ny+nx) - mesh.InsertNextCell(VTK_HEXAHEDRON,ids) - - yield Expected( mesh=mesh ) + ids.SetId( 0, i + j * nx + k * nx * ny ) + ids.SetId( 1, i + j * nx + k * nx * ny + 1 ) + ids.SetId( 2, i + j * nx + k * nx * ny + nx * ny + 1 ) + ids.SetId( 3, i + j * nx + k * nx * ny + nx * ny ) + ids.SetId( 4, i + j * nx + k * nx * ny + nx ) + ids.SetId( 5, i + j * nx + k * nx * ny + nx + 1 ) + ids.SetId( 6, i + j * nx + k * nx * ny + nx * ny + nx + 1 ) + ids.SetId( 7, i + j * nx + k * nx * ny + nx * ny + nx ) + mesh.InsertNextCell( VTK_HEXAHEDRON, ids ) + + yield Expected( mesh=mesh ) + @pytest.mark.parametrize( "expected", __build_test_mesh() ) def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: output_mesh = transform_mesh( expected.mesh ) assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() - assert output_mesh.GetBounds()[0] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[2] == pytest.approx(0., abs=1e-10) and output_mesh.GetBounds()[4] == pytest.approx(0., abs=1e-10) + assert output_mesh.GetBounds()[ 0 ] == pytest.approx( + 0., abs=1e-10 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( + 0., abs=1e-10 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-10 ) #TODO more assert but need more assumptions then # temp + + # w = vtkXMLUnstructuredGridWriter() # w.SetFileName("./test_rotateAndTranslate.vtu") # w.SetInputData(output_mesh) @@ -97,4 +105,4 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: # w.Write() # if __name__ == "__main__": -# pytest.main([__file__, '-s', '-v','--log-cli-level=WARNING']) \ No newline at end of file +# pytest.main([__file__, '-s', '-v','--log-cli-level=WARNING']) From 0635ef1b005a206f106de64d2a61c507f80b4883 Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 11:56:36 +0200 Subject: [PATCH 05/40] more testing --- geos-mesh/tests/test_rotateAndTranslate.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index 5db8b2484..28b8f850d 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -93,7 +93,18 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: 0., abs=1e-10 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( 0., abs=1e-10 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-10 ) #TODO more assert but need more assumptions then - # temp + # test diagonal + assert np.linalg.norm( np.array( [ output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], + output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], + output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] ] ) ) == pytest.approx( + np.linalg.norm( np.array( [ 5, 8, 2 ] ) ), abs=1e-10 ) + # test aligned with axis + v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v1 = np.array( output_mesh.GetPoint( 10 ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v2 = np.array( output_mesh.GetPoint( 10 * 10 ) ) - np.array( output_mesh.GetPoint( 0 ) ) + assert np.abs( np.dot( v0, v1 ) ) < 1e-10 + assert np.abs( np.dot( v0, v2 ) ) < 1e-10 + assert np.abs( np.dot( v1, v2 ) ) < 1e-10 # w = vtkXMLUnstructuredGridWriter() From 652dd7eddbe2b4701ddf0bb91cd43836221a7ae4 Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 12:06:06 +0200 Subject: [PATCH 06/40] lint and style --- .../src/geos/mesh/processing/rotateAndTranslate.py | 8 +++++++- geos-mesh/tests/test_rotateAndTranslate.py | 10 ++++++---- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py index 842b2804c..25f7728e7 100644 --- a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py +++ b/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py @@ -44,7 +44,13 @@ def transform_mesh( input: vtkUnstructuredGrid, logger: logging.Logger | None = return transformFilter.GetOutput() -def local_frame( pts ): +def local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: + """Find a local frame for a set of points + Args: + pts (np.ndarray): points to find the local frame of + Returns: + tuple[np.ndarray, np.ndarray, np.ndarray]: three orthogonal vectors defining the local frame + """ # find 3 orthogonal vectors # we assume points are on a box # first vector is along x axis diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index 28b8f850d..22627a564 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -94,10 +94,12 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: 0., abs=1e-10 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-10 ) #TODO more assert but need more assumptions then # test diagonal - assert np.linalg.norm( np.array( [ output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], - output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], - output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] ] ) ) == pytest.approx( - np.linalg.norm( np.array( [ 5, 8, 2 ] ) ), abs=1e-10 ) + assert np.linalg.norm( + np.array( [ + output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], + output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], + output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] + ] ) ) == pytest.approx( np.linalg.norm( np.array( [ 5, 8, 2 ] ) ), abs=1e-10 ) # test aligned with axis v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) v1 = np.array( output_mesh.GetPoint( 10 ) ) - np.array( output_mesh.GetPoint( 0 ) ) From a801816497a58571485dab27365b0840a6bad3e9 Mon Sep 17 00:00:00 2001 From: jacques FRANC Date: Fri, 5 Sep 2025 13:17:04 +0200 Subject: [PATCH 07/40] yapf-ing --- geos-mesh/tests/test_rotateAndTranslate.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_rotateAndTranslate.py index 22627a564..9c2c912cf 100644 --- a/geos-mesh/tests/test_rotateAndTranslate.py +++ b/geos-mesh/tests/test_rotateAndTranslate.py @@ -116,6 +116,3 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: # w.SetFileName("./test_rotateAndTranslate_input.vtu") # w.SetInputData(expected.mesh) # w.Write() - -# if __name__ == "__main__": -# pytest.main([__file__, '-s', '-v','--log-cli-level=WARNING']) From 5b051e8eb76182a49409ba75e063b2cb8b4e47b6 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 8 Sep 2025 10:38:02 +0200 Subject: [PATCH 08/40] name conv --- ...tateAndTranslate.py => clipToMainFrame.py} | 35 +++++-------------- ...ndTranslate.py => test_clipToMainFrame.py} | 0 2 files changed, 8 insertions(+), 27 deletions(-) rename geos-mesh/src/geos/mesh/processing/{rotateAndTranslate.py => clipToMainFrame.py} (76%) rename geos-mesh/tests/{test_rotateAndTranslate.py => test_clipToMainFrame.py} (100%) diff --git a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py similarity index 76% rename from geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py rename to geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 25f7728e7..df77bda66 100644 --- a/geos-mesh/src/geos/mesh/processing/rotateAndTranslate.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -21,8 +21,8 @@ def transform_mesh( input: vtkUnstructuredGrid, logger: logging.Logger | None = logging.basicConfig( level=logging.WARNING ) logger = logging.getLogger( __name__ ) - translation, theta, axis = recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( input.GetPoints().GetData() ), - logger ) + translation, theta, axis = __recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( input.GetPoints().GetData() ), + logger ) logger.info( f"Theta {theta}" ) logger.info( f"Axis {axis}" ) @@ -44,7 +44,7 @@ def transform_mesh( input: vtkUnstructuredGrid, logger: logging.Logger | None = return transformFilter.GetOutput() -def local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: +def __local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: """Find a local frame for a set of points Args: pts (np.ndarray): points to find the local frame of @@ -54,8 +54,7 @@ def local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray # find 3 orthogonal vectors # we assume points are on a box # first vector is along x axis - ori, _, _ = origin_bounding_box( pts ) - assert ( ( ori == np.asarray( [ 0, 0, 0 ] ) ).all() ) + ori, _, _ = __origin_bounding_box( pts ) u = pts[ 1 ] v = u for pt in pts[ 2: ]: # As we assume points are on a box there should be one orthogonal to u @@ -67,7 +66,7 @@ def local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray # utils -def origin_bounding_box( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: +def __origin_bounding_box( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: """Find the bounding box of a set of points Args: pts (np.ndarray): points to find the bounding box of @@ -80,7 +79,7 @@ def origin_bounding_box( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np. return ( pts[ 0 ], pt0, pt1 ) -def recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np.ndarray, float, np.ndarray ]: +def __recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np.ndarray, float, np.ndarray ]: """Recenter and rotate points to align with principal axes Args: @@ -88,14 +87,14 @@ def recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np. logger (logging.logger, optional): logger instance. """ # find bounding box - org, vmin, vmax = origin_bounding_box( pts ) + org, vmin, vmax = __origin_bounding_box( pts ) logger.info( f"Bounding box is {org}, {vmin}, {vmax}" ) # Transformation translation = org pts -= translation - u, v, w = local_frame( pts ) + u, v, w = __local_frame( pts ) logger.info( f"Local frame u {u}, v {v}, w {w}" ) # find rotation R = U sig V @@ -119,21 +118,3 @@ def recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np. # return translation, rotation, pts return translation, theta, axis - - -# #TO DO : remove -# if __name__ == "__main__": - -# reader = vtk.vtkXMLUnstructuredGridReader() -# reader.SetFileName("/home/jfranc/codes/python/geosPythonPackages/geos-mesh/src/geos/mesh/processing/input.vtu") -# reader.Update() - -# # Transformation -# output = transform_mesh(reader.GetOutput()) - -# output.SetCells(VTK_HEXAHEDRON, reader.GetOutput().GetCells()) -# writer = vtk.vtkXMLUnstructuredGridWriter() -# writer.SetFileName("./output.vtu") -# writer.SetInputData(output) -# writer.Update() -# writer.Write() diff --git a/geos-mesh/tests/test_rotateAndTranslate.py b/geos-mesh/tests/test_clipToMainFrame.py similarity index 100% rename from geos-mesh/tests/test_rotateAndTranslate.py rename to geos-mesh/tests/test_clipToMainFrame.py From fa66f4be3e79d44c75c0762cb1ea9c67e38939e4 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 8 Sep 2025 11:20:05 +0200 Subject: [PATCH 09/40] change for full vtkTransform inheritance --- .../geos/mesh/processing/clipToMainFrame.py | 269 ++++++++++-------- geos-mesh/tests/test_clipToMainFrame.py | 18 +- 2 files changed, 168 insertions(+), 119 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index df77bda66..97365e96a 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -1,4 +1,5 @@ from vtkmodules.numpy_interface import dataset_adapter as dsa +from vtkmodules.vtkCommonCore import vtkFloatArray, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid from vtkmodules.vtkCommonTransforms import vtkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter @@ -6,115 +7,159 @@ import logging -def transform_mesh( input: vtkUnstructuredGrid, logger: logging.Logger | None = None ) -> vtkUnstructuredGrid: - """Apply a transformation to a mesh - - Args: - input (vtk.vtkUnstructuredGrid): input mesh - logger (logging.logger | None, optional): logger instance. Defaults to None. - - Returns: - vtk.vtkUnstructuredGrid: transformed mesh - """ - # Initialize the logger if it is empty - if not logger: - logging.basicConfig( level=logging.WARNING ) - logger = logging.getLogger( __name__ ) - - translation, theta, axis = __recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( input.GetPoints().GetData() ), - logger ) - - logger.info( f"Theta {theta}" ) - logger.info( f"Axis {axis}" ) - - # output = vtkUnstructuredGrid() - # (vpts:= vtkPoints()).SetData( dsa.numpy_support.numpy_to_vtk(pts) ) - # output.SetPoints(vpts) - # output.SetCells(input.GetCellTypesArray(), input.GetCellLocationsArray(), input.GetCells()) - - transform = vtkTransform() - transform.PostMultiply() # we want to apply rotation then translation - transform.RotateWXYZ( -theta, axis[ 0 ], axis[ 1 ], axis[ 2 ] ) - transform.Translate( -translation[ 0 ], -translation[ 1 ], -translation[ 2 ] ) - transformFilter = vtkTransformFilter() - transformFilter.SetTransform( transform ) - transformFilter.SetInputData( input ) - transformFilter.Update() - - return transformFilter.GetOutput() - - -def __local_frame( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: - """Find a local frame for a set of points - Args: - pts (np.ndarray): points to find the local frame of - Returns: - tuple[np.ndarray, np.ndarray, np.ndarray]: three orthogonal vectors defining the local frame - """ - # find 3 orthogonal vectors - # we assume points are on a box - # first vector is along x axis - ori, _, _ = __origin_bounding_box( pts ) - u = pts[ 1 ] - v = u - for pt in pts[ 2: ]: # As we assume points are on a box there should be one orthogonal to u - if ( np.abs( np.dot( u, pt ) ) < 0.0001 ): - v = pt - break - - return ( u, v, np.cross( u, v ) ) - - - # utils -def __origin_bounding_box( pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: - """Find the bounding box of a set of points - Args: - pts (np.ndarray): points to find the bounding box of - Returns: - tuple[np.ndarray, np.ndarray, np.ndarray]: origin, min, max of the bounding box - """ - pt0 = pts[ np.argmin( pts, axis=0 ), : ] - pt1 = pts[ np.argmax( pts, axis=0 ), : ] - - return ( pts[ 0 ], pt0, pt1 ) - - -def __recenter_and_rotate( pts: np.ndarray, logger: logging.Logger ) -> tuple[ np.ndarray, float, np.ndarray ]: - """Recenter and rotate points to align with principal axes - - Args: - pts (np.ndarray): points to recenter and rotate - logger (logging.logger, optional): logger instance. - """ - # find bounding box - org, vmin, vmax = __origin_bounding_box( pts ) - logger.info( f"Bounding box is {org}, {vmin}, {vmax}" ) - - # Transformation - translation = org - pts -= translation - - u, v, w = __local_frame( pts ) - logger.info( f"Local frame u {u}, v {v}, w {w}" ) - - # find rotation R = U sig V - rotation = np.asarray( [ u / np.linalg.norm( u ), v / np.linalg.norm( v ), w / np.linalg.norm( w ) ], - dtype=np.float64 ).transpose() - logger.info( f"R {rotation}" ) - - theta = np.acos( .5 * ( np.trace( rotation ) - 1 ) ) * 180 / np.pi - logger.info( f"theta {theta}" ) - - axis = np.asarray( [ - rotation[ 2, 1 ] - rotation[ 1, 2 ], rotation[ 0, 2 ] - rotation[ 2, 0 ], rotation[ 1, 0 ] - rotation[ 0, 1 ] - ], - dtype=np.float64 ) - axis /= np.linalg.norm( axis ) - logger.info( f"axis {axis}" ) - - pts = ( rotation.transpose() @ pts.transpose() ).transpose() - pts[ np.abs( pts ) < 1e-15 ] = 0. # clipping points too close to zero - logger.info( f"Un-rotated pts : {pts}" ) - - # return translation, rotation, pts - return translation, theta, axis +class ClipToMainFrame( vtkTransform ): + + userTranslation: vtkFloatArray | None = None # user provided translation + pts: vtkPoints + + def __init__( self, pts: vtkPoints ) -> None: + """Clip mesh to main frame + Args: userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one + + """ + super().__init__() + self.pts = pts + + def Update( self ) -> None: + """Update the transformation""" + self.__transform( self.pts, self.userTranslation ) + super().Update() + + def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: + """Set the user translation + Args: + userTranslation (vtkFloatArray): user translation to set + """ + self.userTranslation = userTranslation + + def GetUserTranslation( self ) -> vtkFloatArray | None: + """Get the user translation + Returns: + vtkFloatArray | None: user translation if set, None otherwise + """ + return self.userTranslation + + def __transform( self, pts: vtkPoints, userTranslation: vtkFloatArray | None = None ) -> None: + """Apply the transformation to the points + Args: + pts (vtkPoints): points to transform + userTranslation (vtkFloatArray | None, optional): if provided use this translation instead of the computed one. Defaults to None. + """ + translation, theta, axis = self.__recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( pts.GetData() ) ) + + if userTranslation is not None: + translation = dsa.numpy_support.vtk_to_numpy( userTranslation ) + logging.info( "Using user translation" ) + + self.PostMultiply() # we want to apply rotation then translation + self.RotateWXYZ( -theta, axis[ 0 ], axis[ 1 ], axis[ 2 ] ) + self.Translate( -translation[ 0 ], -translation[ 1 ], -translation[ 2 ] ) + + logging.info( f"Using translation {translation}" ) + logging.info( f"Theta {theta}" ) + logging.info( f"Axis {axis}" ) + + def __local_frame( self, pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: + """Find a local frame for a set of points + Args: + pts (np.ndarray): points to find the local frame of + Returns: + tuple[np.ndarray, np.ndarray, np.ndarray]: three orthogonal vectors defining the local frame + """ + # find 3 orthogonal vectors + # we assume points are on a box + # first vector is along x axis + ori, _, _ = self.__origin_bounding_box( pts ) + u = pts[ 1 ] + v = u + for pt in pts[ 2: ]: # As we assume points are on a box there should be one orthogonal to u + if ( np.abs( np.dot( u, pt ) ) < 1e-10 ): + v = pt + break + + return ( u, v, np.cross( u, v ) ) + + # utils + def __origin_bounding_box( self, pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: + """Find the bounding box of a set of points + Args: + pts (np.ndarray): points to find the bounding box of + Returns: + tuple[np.ndarray, np.ndarray, np.ndarray]: origin, min, max of the bounding box + """ + pt0 = pts[ np.argmin( pts, axis=0 ), : ] + pt1 = pts[ np.argmax( pts, axis=0 ), : ] + + print( f"pt0 {pt0}, pt1 {pt1}" ) + + return ( pts[ 0 ], pt0, pt1 ) + + def __recenter_and_rotate( self, pts: np.ndarray ) -> tuple[ np.ndarray, float, np.ndarray ]: + """Recenter and rotate points to align with principal axes + + Args: + pts (np.ndarray): points to recenter and rotate + """ + # find bounding box + org, vmin, vmax = self.__origin_bounding_box( pts ) + logging.info( f"Bounding box is {org}, {vmin}, {vmax}" ) + + # Transformation + translation = org + pts -= translation + + u, v, w = self.__local_frame( pts ) + logging.info( f"Local frame u {u}, v {v}, w {w}" ) + + # find rotation R = U sig V + rotation = np.asarray( [ u / np.linalg.norm( u ), v / np.linalg.norm( v ), w / np.linalg.norm( w ) ], + dtype=np.float64 ).transpose() + logging.info( f"R {rotation}" ) + + theta = np.acos( .5 * ( np.trace( rotation ) - 1 ) ) * 180 / np.pi + logging.info( f"theta {theta}" ) + + axis = np.asarray( [ + rotation[ 2, 1 ] - rotation[ 1, 2 ], rotation[ 0, 2 ] - rotation[ 2, 0 ], + rotation[ 1, 0 ] - rotation[ 0, 1 ] + ], + dtype=np.float64 ) + axis /= np.linalg.norm( axis ) + logging.info( f"axis {axis}" ) + + pts = ( rotation.transpose() @ pts.transpose() ).transpose() + pts[ np.abs( pts ) < 1e-15 ] = 0. # clipping points too close to zero + logging.info( f"Un-rotated pts : {pts}" ) + + # return translation, rotation, pts + return translation, theta, axis + + +class ClipToMainFrameFilter( vtkTransformFilter ): + + userTranslation: vtkFloatArray | None = None # user provided translation + + def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: + """Set the user translation + Args: + userTranslation (vtkFloatArray): user translation to set + """ + self.userTranslation = userTranslation + + def GetUserTranslation( self ) -> vtkFloatArray | None: + """Get the user translation + Returns: + vtkFloatArray | None: user translation if set, None otherwise + """ + return self.userTranslation + + def SetInputData( self, input: vtkUnstructuredGrid ) -> None: + """Set the input data and apply the transformation to it + Args: + input (vtkUnstructuredGrid): input mesh to transform + """ + super().SetInputData( input ) + clip = ClipToMainFrame( input.GetPoints() ) + clip.SetUserTranslation( self.userTranslation ) + clip.Update() + self.SetTransform( clip ) diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 9c2c912cf..6d1bd7e36 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -8,7 +8,10 @@ import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON -from geos.mesh.processing.rotateAndTranslate import transform_mesh +from geos.mesh.processing.clipToMainFrame import ClipToMainFrameFilter + +Lx, Ly, Lz = 5, 2, 8 +nx, ny, nz = 10, 10, 10 @dataclass( frozen=True ) @@ -45,8 +48,7 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz # np.random.seed(1) # for reproducibility np.random.default_rng() - Lx, Ly, Lz = ( 2, 5, 8 ) # box size - nx, ny, nz = ( 10, 10, 10 ) # number of points in each direction + pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz ) logging.warning( f"Offseting of {off}" ) @@ -86,7 +88,9 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: @pytest.mark.parametrize( "expected", __build_test_mesh() ) def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: - output_mesh = transform_mesh( expected.mesh ) + ( filter := ClipToMainFrameFilter() ).SetInputData( expected.mesh ) + filter.Update() + output_mesh = filter.GetOutput() assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() assert output_mesh.GetBounds()[ 0 ] == pytest.approx( @@ -99,11 +103,11 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] - ] ) ) == pytest.approx( np.linalg.norm( np.array( [ 5, 8, 2 ] ) ), abs=1e-10 ) + ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-10 ) # test aligned with axis v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v1 = np.array( output_mesh.GetPoint( 10 ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v2 = np.array( output_mesh.GetPoint( 10 * 10 ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v1 = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v2 = np.array( output_mesh.GetPoint( nx * ny ) ) - np.array( output_mesh.GetPoint( 0 ) ) assert np.abs( np.dot( v0, v1 ) ) < 1e-10 assert np.abs( np.dot( v0, v2 ) ) < 1e-10 assert np.abs( np.dot( v1, v2 ) ) < 1e-10 From 8db411475cd57bd8a398947cbfbab9751952f923 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 8 Sep 2025 14:33:09 +0200 Subject: [PATCH 10/40] loglevel changes and clean up --- geos-mesh/src/geos/mesh/processing/clipToMainFrame.py | 6 +----- geos-mesh/tests/test_clipToMainFrame.py | 10 +++++----- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 97365e96a..43f0c7b07 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -90,8 +90,6 @@ def __origin_bounding_box( self, pts: np.ndarray ) -> tuple[ np.ndarray, np.ndar pt0 = pts[ np.argmin( pts, axis=0 ), : ] pt1 = pts[ np.argmax( pts, axis=0 ), : ] - print( f"pt0 {pt0}, pt1 {pt1}" ) - return ( pts[ 0 ], pt0, pt1 ) def __recenter_and_rotate( self, pts: np.ndarray ) -> tuple[ np.ndarray, float, np.ndarray ]: @@ -117,7 +115,7 @@ def __recenter_and_rotate( self, pts: np.ndarray ) -> tuple[ np.ndarray, float, logging.info( f"R {rotation}" ) theta = np.acos( .5 * ( np.trace( rotation ) - 1 ) ) * 180 / np.pi - logging.info( f"theta {theta}" ) + logging.info( f"Theta {theta}" ) axis = np.asarray( [ rotation[ 2, 1 ] - rotation[ 1, 2 ], rotation[ 0, 2 ] - rotation[ 2, 0 ], @@ -125,11 +123,9 @@ def __recenter_and_rotate( self, pts: np.ndarray ) -> tuple[ np.ndarray, float, ], dtype=np.float64 ) axis /= np.linalg.norm( axis ) - logging.info( f"axis {axis}" ) pts = ( rotation.transpose() @ pts.transpose() ).transpose() pts[ np.abs( pts ) < 1e-15 ] = 0. # clipping points too close to zero - logging.info( f"Un-rotated pts : {pts}" ) # return translation, rotation, pts return translation, theta, axis diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 6d1bd7e36..e34df2b08 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -51,16 +51,16 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz ) - logging.warning( f"Offseting of {off}" ) - logging.warning( f"Original pts : {pts}" ) + logging.info( f"Offseting of {off}" ) + logging.debug( f"Original pts : {pts}" ) angles = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad - logging.warning( f"angles {angles[0]}" ) + logging.info( f"angles {angles[0]}" ) pts = __rotate_box( angles[ 0 ], pts ) - logging.info( f"Rotated pts : {pts}" ) + logging.debug( f"Rotated pts : {pts}" ) pts[ :, 0 ] += off[ 0 ][ 0 ] pts[ :, 1 ] += off[ 0 ][ 1 ] pts[ :, 2 ] += off[ 0 ][ 2 ] - logging.info( f"Translated pts : {pts}" ) + logging.debug( f"Translated pts : {pts}" ) # Creating multiple meshes, each time with a different angles mesh = vtkUnstructuredGrid() From 29d549f3b6f111a46f7684e82f96482ec14682f8 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 8 Sep 2025 14:47:39 +0200 Subject: [PATCH 11/40] file's header --- geos-mesh/src/geos/mesh/processing/clipToMainFrame.py | 3 +++ geos-mesh/tests/test_clipToMainFrame.py | 4 ++++ 2 files changed, 7 insertions(+) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 43f0c7b07..022fc9cb2 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -1,3 +1,6 @@ +# SPDX-FileContributor: Jacques Franc +# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-License-Identifier: Apache 2.0 from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkCommonCore import vtkFloatArray, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index e34df2b08..75a29b054 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -1,3 +1,7 @@ +# SPDX-FileContributor: Jacques Franc +# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-License-Identifier: Apache 2.0 +# ruff: noqa: E402 # disable Module level import not at top of file import pytest from dataclasses import dataclass from typing import Generator From 0a2f1cb05d2ca886f699015e43c82de3177a00c2 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 9 Sep 2025 11:44:33 +0200 Subject: [PATCH 12/40] using vtkLandmarkTransform and vtkOBBTree --- .../geos/mesh/processing/clipToMainFrame.py | 274 ++++++++++++------ geos-mesh/tests/test_clipToMainFrame.py | 57 +++- 2 files changed, 225 insertions(+), 106 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 022fc9cb2..121a2eb2c 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -1,33 +1,44 @@ -# SPDX-FileContributor: Jacques Franc -# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-FileContributor: Jacques Franc +# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-License-Identifier: Apache 2.0 from vtkmodules.numpy_interface import dataset_adapter as dsa -from vtkmodules.vtkCommonCore import vtkFloatArray, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid -from vtkmodules.vtkCommonTransforms import vtkTransform +from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray +from vtkmodules.vtkFiltersGeneral import vtkOBBTree +from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkDataObjectTreeIterator, vtkPolyData +from vtkmodules.vtkCommonTransforms import vtkLandmarkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter + +from geos.mesh.utils.genericHelpers import getMultiBlockBounds + import numpy as np import logging -class ClipToMainFrame( vtkTransform ): +class ClipToMainFrame( vtkLandmarkTransform ): userTranslation: vtkFloatArray | None = None # user provided translation - pts: vtkPoints + sourcePts: vtkPoints + targetPts: vtkPoints + mesh: vtkUnstructuredGrid - def __init__( self, pts: vtkPoints ) -> None: + def __init__( self, mesh: vtkUnstructuredGrid ) -> None: """Clip mesh to main frame Args: userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one """ super().__init__() - self.pts = pts + self.mesh = mesh def Update( self ) -> None: """Update the transformation""" - self.__transform( self.pts, self.userTranslation ) + self.sourcePts, self.targetPts = self.__getFramePoints( self.__getOBBTree( self.mesh ) ) + self.SetSourceLandmarks( self.sourcePts ) + self.SetTargetLandmarks( self.targetPts ) + self.SetModeToRigidBody() super().Update() + def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: """Set the user translation Args: @@ -42,99 +53,127 @@ def GetUserTranslation( self ) -> vtkFloatArray | None: """ return self.userTranslation - def __transform( self, pts: vtkPoints, userTranslation: vtkFloatArray | None = None ) -> None: - """Apply the transformation to the points - Args: - pts (vtkPoints): points to transform - userTranslation (vtkFloatArray | None, optional): if provided use this translation instead of the computed one. Defaults to None. - """ - translation, theta, axis = self.__recenter_and_rotate( dsa.numpy_support.vtk_to_numpy( pts.GetData() ) ) + def __str__(self): - if userTranslation is not None: - translation = dsa.numpy_support.vtk_to_numpy( userTranslation ) - logging.info( "Using user translation" ) - - self.PostMultiply() # we want to apply rotation then translation - self.RotateWXYZ( -theta, axis[ 0 ], axis[ 1 ], axis[ 2 ] ) - self.Translate( -translation[ 0 ], -translation[ 1 ], -translation[ 2 ] ) + return super().__str__() + f"\nUser translation: {self.userTranslation}" \ + + f"\nSource points: {self.sourcePts}" \ + + f"\nTarget points: {self.targetPts}" \ + + f"\nAngle-Axis: {self.__getAngleAxis()}" \ + + f"\nTranslation: {self._getTranslation()}" - logging.info( f"Using translation {translation}" ) - logging.info( f"Theta {theta}" ) - logging.info( f"Axis {axis}" ) - def __local_frame( self, pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: - """Find a local frame for a set of points - Args: - pts (np.ndarray): points to find the local frame of + def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: + """Get the angle and axis of the rotation Returns: - tuple[np.ndarray, np.ndarray, np.ndarray]: three orthogonal vectors defining the local frame + tuple[float, np.ndarray]: angle in degrees and axis of rotation """ - # find 3 orthogonal vectors - # we assume points are on a box - # first vector is along x axis - ori, _, _ = self.__origin_bounding_box( pts ) - u = pts[ 1 ] - v = u - for pt in pts[ 2: ]: # As we assume points are on a box there should be one orthogonal to u - if ( np.abs( np.dot( u, pt ) ) < 1e-10 ): - v = pt - break - - return ( u, v, np.cross( u, v ) ) - - # utils - def __origin_bounding_box( self, pts: np.ndarray ) -> tuple[ np.ndarray, np.ndarray, np.ndarray ]: - """Find the bounding box of a set of points - Args: - pts (np.ndarray): points to find the bounding box of + matrix = self.GetMatrix() + angle = np.arccos( ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) + if angle == 0: + return 0.0, np.array( [ 1.0, 0.0, 0.0 ] ) + rx = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) + ry = matrix.GetElement( 0, 2 ) - matrix.GetElement( 2, 0 ) + rz = matrix.GetElement( 1, 0 ) - matrix.GetElement( 0, 1 ) + r = np.array( [ rx, ry, rz ] ) + r /= np.linalg.norm( r ) + return np.degrees( angle ), r + + def _getTranslation( self ) -> np.ndarray: + """Get the translation vector Returns: - tuple[np.ndarray, np.ndarray, np.ndarray]: origin, min, max of the bounding box + np.ndarray: translation vector """ - pt0 = pts[ np.argmin( pts, axis=0 ), : ] - pt1 = pts[ np.argmax( pts, axis=0 ), : ] - - return ( pts[ 0 ], pt0, pt1 ) + matrix = self.GetMatrix() + return np.array( [ matrix.GetElement( 0, 3 ), matrix.GetElement( 1, 3 ), matrix.GetElement( 2, 3 ) ] ) - def __recenter_and_rotate( self, pts: np.ndarray ) -> tuple[ np.ndarray, float, np.ndarray ]: - """Recenter and rotate points to align with principal axes - + def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: + """Get the OBB tree of the mesh Args: - pts (np.ndarray): points to recenter and rotate + mesh (vtkUnstructuredGrid): mesh to get the OBB tree from + Returns: + vtkOBBTree: OBB tree of the mesh """ - # find bounding box - org, vmin, vmax = self.__origin_bounding_box( pts ) - logging.info( f"Bounding box is {org}, {vmin}, {vmax}" ) - - # Transformation - translation = org - pts -= translation + OBBTree = vtkOBBTree() + surfFilter = vtkDataSetSurfaceFilter() + surfFilter.SetInputData( mesh ) + surfFilter.Update() + OBBTree.SetDataSet( surfFilter.GetOutput() ) + OBBTree.BuildLocator() + pdata = vtkPolyData() + OBBTree.GenerateRepresentation( 0, pdata) + # at level 0 this should return 8 corners of the bounding box + if pdata.GetNumberOfPoints() < 3: + logging.warning( "Could not get OBB points, using bounding box points instead" ) + return self.__allpoints( mesh.GetBounds() ) - u, v, w = self.__local_frame( pts ) - logging.info( f"Local frame u {u}, v {v}, w {w}" ) + return pdata.GetPoints() - # find rotation R = U sig V - rotation = np.asarray( [ u / np.linalg.norm( u ), v / np.linalg.norm( v ), w / np.linalg.norm( w ) ], - dtype=np.float64 ).transpose() - logging.info( f"R {rotation}" ) - theta = np.acos( .5 * ( np.trace( rotation ) - 1 ) ) * 180 / np.pi - logging.info( f"Theta {theta}" ) + def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] ) -> vtkPoints: + """Get the 8 corners of the bounding box + Args: + bounds (tuple[float, float, float, float, float, float]): bounding box + Returns: + vtkPoints: 8 corners of the bounding box + """ + pts = vtkPoints() + pts.SetNumberOfPoints( 8 ) + pts.SetPoint( 0, [ bounds[ 0 ], bounds[ 2 ], bounds[ 4 ] ] ) + pts.SetPoint( 1, [ bounds[ 1 ], bounds[ 2 ], bounds[ 4 ] ] ) + pts.SetPoint( 2, [ bounds[ 1 ], bounds[ 3 ], bounds[ 4 ] ] ) + pts.SetPoint( 3, [ bounds[ 0 ], bounds[ 3 ], bounds[ 4 ] ] ) + pts.SetPoint( 4, [ bounds[ 0 ], bounds[ 2 ], bounds[ 5 ] ] ) + pts.SetPoint( 5, [ bounds[ 1 ], bounds[ 2 ], bounds[ 5 ] ] ) + pts.SetPoint( 6, [ bounds[ 1 ], bounds[ 3 ], bounds[ 5 ] ] ) + pts.SetPoint( 7, [ bounds[ 0 ], bounds[ 3 ], bounds[ 5 ] ] ) + return pts - axis = np.asarray( [ - rotation[ 2, 1 ] - rotation[ 1, 2 ], rotation[ 0, 2 ] - rotation[ 2, 0 ], - rotation[ 1, 0 ] - rotation[ 0, 1 ] - ], - dtype=np.float64 ) - axis /= np.linalg.norm( axis ) + def __getFramePoints( self, vpts : vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: + """Get the source and target points for the transformation + Args: + pts (vtkPoints): points to transform + Returns: + tuple[vtkPoints, vtkPoints]: source and target points for the transformation + """ + pts = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) + further_ix = np.argmax(np.linalg.norm(pts,axis=1)) # by default take the min point furthest from origin + org = pts[further_ix,:] + if self.userTranslation is not None: + org = dsa.numpy_support.vtk_to_numpy( self.userTranslation ) + logging.info( "Using user translation" ) + logging.info( f"Moving point {org} to origin for transformation" ) + # find 3 orthogonal vectors + # we assume points are on a box + dist_indexes = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) + # find u,v,w + v1 = pts[dist_indexes[1],:] - org + v1 /= np.linalg.norm( v1 ) + v2 = pts[dist_indexes[2],:] - org + v2 /= np.linalg.norm( v2 ) + # ensure orthogonality + v2 -= np.dot( v2, v1 ) * v1 + v2 /= np.linalg.norm( v2 ) + v3 = np.cross( v1, v2 ) + v3 /= np.linalg.norm( v3 ) - pts = ( rotation.transpose() @ pts.transpose() ).transpose() - pts[ np.abs( pts ) < 1e-15 ] = 0. # clipping points too close to zero + sourcePts = vtkPoints() + sourcePts.SetNumberOfPoints( 4 ) + sourcePts.SetPoint( 0, org ) + sourcePts.SetPoint( 1, v1 + org ) + sourcePts.SetPoint( 2, v2 + org ) + sourcePts.SetPoint( 3, v3 + org ) - # return translation, rotation, pts - return translation, theta, axis + targetPts = vtkPoints() + targetPts.SetNumberOfPoints( 4 ) + targetPts.SetPoint( 0, [ 0., 0., 0. ] ) + targetPts.SetPoint( 1, [ 1., 0., 0. ] ) + targetPts.SetPoint( 2, [ 0., 1. , 0. ] ) + targetPts.SetPoint( 3, [ 0., 0., 1. ]) + return ( sourcePts, targetPts ) class ClipToMainFrameFilter( vtkTransformFilter ): + """Filter to clip a mesh to the main frame using ClipToMainFrame class""" userTranslation: vtkFloatArray | None = None # user provided translation @@ -152,13 +191,66 @@ def GetUserTranslation( self ) -> vtkFloatArray | None: """ return self.userTranslation - def SetInputData( self, input: vtkUnstructuredGrid ) -> None: - """Set the input data and apply the transformation to it - Args: - input (vtkUnstructuredGrid): input mesh to transform - """ - super().SetInputData( input ) - clip = ClipToMainFrame( input.GetPoints() ) + def Update( self ) -> None: + """Update the transformation""" + # dispatch to ClipToMainFrame depending on input type + if isinstance( self.GetInput(), vtkMultiBlockDataSet ): + #locate reference point + idBlock = self.__locate_reference_point( self.GetInput(), self.userTranslation ) + clip = ClipToMainFrame( self.GetInput().GetBlock( idBlock - 1 ) ) + else: + clip = ClipToMainFrame( self.GetInput() ) + clip.SetUserTranslation( self.userTranslation ) clip.Update() self.SetTransform( clip ) + super().Update() + + + + def __locate_reference_point( self, input: vtkMultiBlockDataSet, userTranslation: vtkFloatArray | None ) -> int: + """Locate the block to use as reference for the transformation + Args: + input (vtkMultiBlockDataSet): input multiblock mesh + userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one + Returns: + int: index of the block to use as reference + """ + + def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, float ] ) -> bool: + """Check if a point is inside a bounding box + Args: + pt (np.ndarray): point to check + bounds (tuple[float, float, float, float, float, float]): bounding box + Returns: + bool: True if the point is inside the bounding box, False otherwise + """ + logging.info( f"Checking if point {pt} is inside bounds {bounds}" ) + return ( pt[ 0 ] >= bounds[ 0 ] and pt[ 0 ] <= bounds[ 1 ] and pt[ 1 ] >= bounds[ 2 ] and pt[ 1 ] <= bounds[ 3 ] + and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) + + #TODO (jacques) make a decorator for this + iter : vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() + iter.SetDataSet( input ) + iter.VisitOnlyLeavesOn() + iter.GoToFirstItem() + xmin, _ , ymin, _ , zmin, _ = getMultiBlockBounds( input ) + #TODO (jacques) : rewrite with a filter struct + while iter.GetCurrentDataObject() is not None: + block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) + if block.GetNumberOfPoints() > 0: + bounds = block.GetBounds() + if userTranslation is not None: + #check if user translation is inside the block + if __inside( dsa.numpy_support.vtk_to_numpy(self.userTranslation), bounds ): + logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) + return iter.GetCurrentFlatIndex() + else: + #use the lowest bounds corner as reference point + if __inside( np.asarray([xmin,ymin,zmin]) , bounds ): + logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) + return iter.GetCurrentFlatIndex() + iter.GoToNextItem() + + + diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 75a29b054..e505b9511 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -1,12 +1,13 @@ -# SPDX-FileContributor: Jacques Franc -# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-FileContributor: Jacques Franc +# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-License-Identifier: Apache 2.0 # ruff: noqa: E402 # disable Module level import not at top of file import pytest from dataclasses import dataclass from typing import Generator from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, VTK_HEXAHEDRON +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron +from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np @@ -17,7 +18,6 @@ Lx, Ly, Lz = 5, 2, 8 nx, ny, nz = 10, 10, 10 - @dataclass( frozen=True ) class Expected: mesh: vtkUnstructuredGrid @@ -98,16 +98,15 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() assert output_mesh.GetBounds()[ 0 ] == pytest.approx( - 0., abs=1e-10 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( - 0., abs=1e-10 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-10 ) - #TODO more assert but need more assumptions then + 0., abs=1e-6 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( + 0., abs=1e-6 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-6 ) # test diagonal assert np.linalg.norm( np.array( [ output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] - ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-10 ) + ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-5 ) # test aligned with axis v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) v1 = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) @@ -116,11 +115,39 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: assert np.abs( np.dot( v0, v2 ) ) < 1e-10 assert np.abs( np.dot( v1, v2 ) ) < 1e-10 + #check if input has been modified + # w = vtkXMLUnstructuredGridWriter() + # w.SetFileName("./test_rotateAndTranslate_input.vtu") + # w.SetInputData(expected.mesh) + # w.Write() + + # w.SetFileName("./test_rotateAndTranslate.vtu") + # w.SetInputData(output_mesh) + # w.Write() + +def test_rotateAndTranslate_gen(): + + reader = vtkXMLMultiBlockDataReader() + reader.SetFileName( "/data/pau901/SIM_CS/users/MargauxRaguenel/data/Gengibre/SKUA/gengibre_2faults_res100.vtm" ) + reader.Update() + input_mesh = reader.GetOutput() + ( filter := ClipToMainFrameFilter() ).SetInputData( input_mesh ) + filter.Update() + print(filter.GetTransform()) + output_mesh = filter.GetOutputDataObject( 0 ) + assert output_mesh.GetNumberOfPoints() == input_mesh.GetNumberOfPoints() + assert output_mesh.GetNumberOfCells() == input_mesh.GetNumberOfCells() + + # #check if input has been modified + # w = vtkXMLMultiBlockDataWriter() + # w.SetFileName("./test_rotateAndTranslate_input.vtm") + # w.SetInputData(input_mesh) + # w.Write() + + + # w.SetFileName("./test_rotateAndTranslate.vtm") + # w.SetInputData(output_mesh) + # w.Write() -# w = vtkXMLUnstructuredGridWriter() -# w.SetFileName("./test_rotateAndTranslate.vtu") -# w.SetInputData(output_mesh) -# w.Write() -# w.SetFileName("./test_rotateAndTranslate_input.vtu") -# w.SetInputData(expected.mesh) -# w.Write() +# if __name__ == "__main__": +# pytest.main( [ __file__, "-s", "-vvv", "--log-cli-level=INFO" ] ) From 9be24e96caa8caf674854be214165820a3c5bec8 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 9 Sep 2025 17:39:44 +0200 Subject: [PATCH 13/40] some doc and accessible tests --- geos-mesh/src/geos/mesh/processing/clipToMainFrame.py | 5 +++++ geos-mesh/tests/test_clipToMainFrame.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 121a2eb2c..50a6b04a7 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -14,6 +14,11 @@ import numpy as np import logging +__doc__ = """Module to clip a mesh to the main frame using rigid body transformation. +Methods include: + - ClipToMainFrame class to compute the transformation + - ClipToMainFrameFilter class to apply the transformation to a mesh +""" class ClipToMainFrame( vtkLandmarkTransform ): diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index e505b9511..b8eec4228 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -128,7 +128,7 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: def test_rotateAndTranslate_gen(): reader = vtkXMLMultiBlockDataReader() - reader.SetFileName( "/data/pau901/SIM_CS/users/MargauxRaguenel/data/Gengibre/SKUA/gengibre_2faults_res100.vtm" ) + reader.SetFileName( "./geos-mesh/tests/data/displacedFault.vtm" ) reader.Update() input_mesh = reader.GetOutput() ( filter := ClipToMainFrameFilter() ).SetInputData( input_mesh ) From 29ba247fafb309d55bd645ce466eb5adad0dfdc2 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 10 Sep 2025 09:15:25 +0200 Subject: [PATCH 14/40] ruff and yapf --- .../geos/mesh/processing/clipToMainFrame.py | 136 ++++++++++-------- geos-mesh/tests/test_clipToMainFrame.py | 25 ++-- 2 files changed, 86 insertions(+), 75 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py index 50a6b04a7..50ec3264b 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py @@ -20,6 +20,7 @@ - ClipToMainFrameFilter class to apply the transformation to a mesh """ + class ClipToMainFrame( vtkLandmarkTransform ): userTranslation: vtkFloatArray | None = None # user provided translation @@ -28,52 +29,55 @@ class ClipToMainFrame( vtkLandmarkTransform ): mesh: vtkUnstructuredGrid def __init__( self, mesh: vtkUnstructuredGrid ) -> None: - """Clip mesh to main frame - Args: userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one - + """Clip mesh to main frame. + + Args: + mesh (vtkUnstructuredGrid): mesh to clip. """ super().__init__() self.mesh = mesh def Update( self ) -> None: - """Update the transformation""" + """Update the transformation.""" self.sourcePts, self.targetPts = self.__getFramePoints( self.__getOBBTree( self.mesh ) ) self.SetSourceLandmarks( self.sourcePts ) self.SetTargetLandmarks( self.targetPts ) self.SetModeToRigidBody() super().Update() - def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: - """Set the user translation + """Set the user translation. + Args: - userTranslation (vtkFloatArray): user translation to set + userTranslation (vtkFloatArray): user translation to set. """ self.userTranslation = userTranslation def GetUserTranslation( self ) -> vtkFloatArray | None: - """Get the user translation + """Get the user translation. + Returns: - vtkFloatArray | None: user translation if set, None otherwise + vtkFloatArray | None: user translation if set, None otherwise. """ return self.userTranslation - def __str__(self): - + def __str__( self ) -> str: + """String representation of the transformation.""" return super().__str__() + f"\nUser translation: {self.userTranslation}" \ + f"\nSource points: {self.sourcePts}" \ + f"\nTarget points: {self.targetPts}" \ + f"\nAngle-Axis: {self.__getAngleAxis()}" \ + f"\nTranslation: {self._getTranslation()}" - def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: - """Get the angle and axis of the rotation + """Get the angle and axis of the rotation. + Returns: - tuple[float, np.ndarray]: angle in degrees and axis of rotation + tuple[float, np.ndarray]: angle in degrees and axis of rotation. """ matrix = self.GetMatrix() - angle = np.arccos( ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) + angle = np.arccos( + ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) if angle == 0: return 0.0, np.array( [ 1.0, 0.0, 0.0 ] ) rx = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) @@ -84,19 +88,22 @@ def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: return np.degrees( angle ), r def _getTranslation( self ) -> np.ndarray: - """Get the translation vector + """Get the translation vector. + Returns: - np.ndarray: translation vector + np.ndarray: translation vector. """ matrix = self.GetMatrix() return np.array( [ matrix.GetElement( 0, 3 ), matrix.GetElement( 1, 3 ), matrix.GetElement( 2, 3 ) ] ) def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: - """Get the OBB tree of the mesh + """Get the OBB tree of the mesh. + Args: - mesh (vtkUnstructuredGrid): mesh to get the OBB tree from + mesh (vtkUnstructuredGrid): mesh to get the OBB tree from. + Returns: - vtkOBBTree: OBB tree of the mesh + vtkOBBTree: OBB tree of the mesh. """ OBBTree = vtkOBBTree() surfFilter = vtkDataSetSurfaceFilter() @@ -105,7 +112,7 @@ def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: OBBTree.SetDataSet( surfFilter.GetOutput() ) OBBTree.BuildLocator() pdata = vtkPolyData() - OBBTree.GenerateRepresentation( 0, pdata) + OBBTree.GenerateRepresentation( 0, pdata ) # at level 0 this should return 8 corners of the bounding box if pdata.GetNumberOfPoints() < 3: logging.warning( "Could not get OBB points, using bounding box points instead" ) @@ -113,13 +120,14 @@ def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: return pdata.GetPoints() - def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] ) -> vtkPoints: - """Get the 8 corners of the bounding box + """Get the 8 corners of the bounding box. + Args: - bounds (tuple[float, float, float, float, float, float]): bounding box + bounds (tuple[float, float, float, float, float, float]): bounding box. + Returns: - vtkPoints: 8 corners of the bounding box + vtkPoints: 8 corners of the bounding box. """ pts = vtkPoints() pts.SetNumberOfPoints( 8 ) @@ -133,16 +141,18 @@ def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] pts.SetPoint( 7, [ bounds[ 0 ], bounds[ 3 ], bounds[ 5 ] ] ) return pts - def __getFramePoints( self, vpts : vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: - """Get the source and target points for the transformation + def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: + """Get the source and target points for the transformation. + Args: - pts (vtkPoints): points to transform + vpts (vtkPoints): points to transform. + Returns: - tuple[vtkPoints, vtkPoints]: source and target points for the transformation + tuple[vtkPoints, vtkPoints]: source and target points for the transformation. """ pts = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) - further_ix = np.argmax(np.linalg.norm(pts,axis=1)) # by default take the min point furthest from origin - org = pts[further_ix,:] + further_ix = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin + org = pts[ further_ix, : ] if self.userTranslation is not None: org = dsa.numpy_support.vtk_to_numpy( self.userTranslation ) logging.info( "Using user translation" ) @@ -151,9 +161,9 @@ def __getFramePoints( self, vpts : vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: # we assume points are on a box dist_indexes = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) # find u,v,w - v1 = pts[dist_indexes[1],:] - org + v1 = pts[ dist_indexes[ 1 ], : ] - org v1 /= np.linalg.norm( v1 ) - v2 = pts[dist_indexes[2],:] - org + v2 = pts[ dist_indexes[ 2 ], : ] - org v2 /= np.linalg.norm( v2 ) # ensure orthogonality v2 -= np.dot( v2, v1 ) * v1 @@ -164,40 +174,43 @@ def __getFramePoints( self, vpts : vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: sourcePts = vtkPoints() sourcePts.SetNumberOfPoints( 4 ) sourcePts.SetPoint( 0, org ) - sourcePts.SetPoint( 1, v1 + org ) - sourcePts.SetPoint( 2, v2 + org ) - sourcePts.SetPoint( 3, v3 + org ) + sourcePts.SetPoint( 1, v1 + org ) + sourcePts.SetPoint( 2, v2 + org ) + sourcePts.SetPoint( 3, v3 + org ) targetPts = vtkPoints() targetPts.SetNumberOfPoints( 4 ) targetPts.SetPoint( 0, [ 0., 0., 0. ] ) targetPts.SetPoint( 1, [ 1., 0., 0. ] ) - targetPts.SetPoint( 2, [ 0., 1. , 0. ] ) - targetPts.SetPoint( 3, [ 0., 0., 1. ]) + targetPts.SetPoint( 2, [ 0., 1., 0. ] ) + targetPts.SetPoint( 3, [ 0., 0., 1. ] ) return ( sourcePts, targetPts ) + class ClipToMainFrameFilter( vtkTransformFilter ): - """Filter to clip a mesh to the main frame using ClipToMainFrame class""" + """Filter to clip a mesh to the main frame using ClipToMainFrame class.""" userTranslation: vtkFloatArray | None = None # user provided translation def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: - """Set the user translation + """Set the user translation. + Args: - userTranslation (vtkFloatArray): user translation to set + userTranslation (vtkFloatArray): user translation to set. """ self.userTranslation = userTranslation def GetUserTranslation( self ) -> vtkFloatArray | None: - """Get the user translation + """Get the user translation. + Returns: - vtkFloatArray | None: user translation if set, None otherwise + vtkFloatArray | None: user translation if set, None otherwise. """ return self.userTranslation def Update( self ) -> None: - """Update the transformation""" + """Update the transformation.""" # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point @@ -211,35 +224,37 @@ def Update( self ) -> None: self.SetTransform( clip ) super().Update() - - def __locate_reference_point( self, input: vtkMultiBlockDataSet, userTranslation: vtkFloatArray | None ) -> int: - """Locate the block to use as reference for the transformation + """Locate the block to use as reference for the transformation. + Args: - input (vtkMultiBlockDataSet): input multiblock mesh - userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one + input (vtkMultiBlockDataSet): input multiblock mesh. + userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one. + Returns: - int: index of the block to use as reference + int: index of the block to use as reference. """ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, float ] ) -> bool: - """Check if a point is inside a bounding box + """Check if a point is inside a bounding box. + Args: pt (np.ndarray): point to check - bounds (tuple[float, float, float, float, float, float]): bounding box + bounds (tuple[float, float, float, float, float, float]): bounding box. + Returns: - bool: True if the point is inside the bounding box, False otherwise + bool: True if the point is inside the bounding box, False otherwise. """ logging.info( f"Checking if point {pt} is inside bounds {bounds}" ) - return ( pt[ 0 ] >= bounds[ 0 ] and pt[ 0 ] <= bounds[ 1 ] and pt[ 1 ] >= bounds[ 2 ] and pt[ 1 ] <= bounds[ 3 ] - and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) + return ( pt[ 0 ] >= bounds[ 0 ] and pt[ 0 ] <= bounds[ 1 ] and pt[ 1 ] >= bounds[ 2 ] + and pt[ 1 ] <= bounds[ 3 ] and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) #TODO (jacques) make a decorator for this - iter : vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() + iter: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() iter.SetDataSet( input ) iter.VisitOnlyLeavesOn() iter.GoToFirstItem() - xmin, _ , ymin, _ , zmin, _ = getMultiBlockBounds( input ) + xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( input ) #TODO (jacques) : rewrite with a filter struct while iter.GetCurrentDataObject() is not None: block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) @@ -247,15 +262,12 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, bounds = block.GetBounds() if userTranslation is not None: #check if user translation is inside the block - if __inside( dsa.numpy_support.vtk_to_numpy(self.userTranslation), bounds ): + if __inside( dsa.numpy_support.vtk_to_numpy( self.userTranslation ), bounds ): logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) return iter.GetCurrentFlatIndex() else: #use the lowest bounds corner as reference point - if __inside( np.asarray([xmin,ymin,zmin]) , bounds ): + if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) return iter.GetCurrentFlatIndex() iter.GoToNextItem() - - - diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index b8eec4228..a12d54d71 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -7,7 +7,7 @@ from typing import Generator from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron -from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader +from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader, vtkXMLMultiBlockDataWriter import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np @@ -18,6 +18,7 @@ Lx, Ly, Lz = 5, 2, 8 nx, ny, nz = 10, 10, 10 + @dataclass( frozen=True ) class Expected: mesh: vtkUnstructuredGrid @@ -91,7 +92,8 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: @pytest.mark.parametrize( "expected", __build_test_mesh() ) -def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: +def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: + """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" ( filter := ClipToMainFrameFilter() ).SetInputData( expected.mesh ) filter.Update() output_mesh = filter.GetOutput() @@ -125,29 +127,26 @@ def test_rotateAndTranslate_polyhedron( expected: Expected ) -> None: # w.SetInputData(output_mesh) # w.Write() -def test_rotateAndTranslate_gen(): +def test_clipToMainFrame_gen() -> None: + """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" reader = vtkXMLMultiBlockDataReader() - reader.SetFileName( "./geos-mesh/tests/data/displacedFault.vtm" ) + reader.SetFileName( "./geos-mesh/tests/data/displacedFault.vtm" ) reader.Update() input_mesh = reader.GetOutput() ( filter := ClipToMainFrameFilter() ).SetInputData( input_mesh ) filter.Update() - print(filter.GetTransform()) + print( filter.GetTransform() ) output_mesh = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == input_mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == input_mesh.GetNumberOfCells() # #check if input has been modified # w = vtkXMLMultiBlockDataWriter() - # w.SetFileName("./test_rotateAndTranslate_input.vtm") - # w.SetInputData(input_mesh) + # w.SetFileName( "./test_rotateAndTranslate_input.vtm" ) + # w.SetInputData( input_mesh ) # w.Write() - - # w.SetFileName("./test_rotateAndTranslate.vtm") - # w.SetInputData(output_mesh) + # w.SetFileName( "./test_rotateAndTranslate.vtm" ) + # w.SetInputData( output_mesh ) # w.Write() - -# if __name__ == "__main__": -# pytest.main( [ __file__, "-s", "-vvv", "--log-cli-level=INFO" ] ) From e46b174dd94b2bd7a20fa4dbe96241113f6f9ba5 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 10 Sep 2025 10:04:50 +0200 Subject: [PATCH 15/40] docs --- docs/geos_mesh_docs/processing.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/geos_mesh_docs/processing.rst b/docs/geos_mesh_docs/processing.rst index fec7f4b70..03850b769 100644 --- a/docs/geos_mesh_docs/processing.rst +++ b/docs/geos_mesh_docs/processing.rst @@ -26,6 +26,14 @@ geos.mesh.processing.SplitMesh filter -------------------------------------- .. automodule:: geos.mesh.processing.SplitMesh + :members: + :undoc-members: + :show-inheritance: + +geos.mesh.processing.ClipToMainFrame filter +-------------------------------------------- + +.. automodule:: geos.mesh.processing.ClipToMainFrame :members: :undoc-members: :show-inheritance: \ No newline at end of file From 94934511955fa1d1f4d3258a5a8914b1b55bff20 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 10 Sep 2025 10:43:31 +0200 Subject: [PATCH 16/40] fix doc and rename module --- ...{clipToMainFrame.py => ClipToMainFrame.py} | 32 +++++++++++++++++-- geos-mesh/tests/test_clipToMainFrame.py | 2 +- 2 files changed, 30 insertions(+), 4 deletions(-) rename geos-mesh/src/geos/mesh/processing/{clipToMainFrame.py => ClipToMainFrame.py} (94%) diff --git a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py similarity index 94% rename from geos-mesh/src/geos/mesh/processing/clipToMainFrame.py rename to geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 50ec3264b..363df1d48 100644 --- a/geos-mesh/src/geos/mesh/processing/clipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -1,6 +1,7 @@ -# SPDX-FileContributor: Jacques Franc -# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-License-Identifier: Apache 2.0 +# SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-FileContributor: Jacques Franc + from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray from vtkmodules.vtkFiltersGeneral import vtkOBBTree @@ -14,10 +15,35 @@ import numpy as np import logging -__doc__ = """Module to clip a mesh to the main frame using rigid body transformation. +__doc__ = """ +Module to clip a mesh to the main frame using rigid body transformation. + Methods include: - ClipToMainFrame class to compute the transformation - ClipToMainFrameFilter class to apply the transformation to a mesh + +To use it: + +.. code-block:: python + + from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameFilter + + # Filter inputs. + multiBlockDataSet: vtkMultiBlockDataSet + + # Optional inputs. + userTranslation: vtkFloatArray + + # Instantiate the filter. + filter: ClipToMainFrameFilter = ClipToMainFrameFilter() + filter.SetInputData( multiBlockDataSet ) + + # Do calculations. + if userTranslation is not None: + filter.SetUserTranslation( userTranslation ) + filter.Update() + output: vtkMultiBlockDataSet = filter.GetOutput() + """ diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index a12d54d71..fe07f9526 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -13,7 +13,7 @@ import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON -from geos.mesh.processing.clipToMainFrame import ClipToMainFrameFilter +from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameFilter Lx, Ly, Lz = 5, 2, 8 nx, ny, nz = 10, 10, 10 From 43d7e2548c7f198aefc6d578ab685fcd69015d24 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 14:04:09 +0200 Subject: [PATCH 17/40] Fix test and all quadrant --- .../geos/mesh/processing/ClipToMainFrame.py | 26 ++++++++-- geos-mesh/tests/test_clipToMainFrame.py | 50 +++++++------------ 2 files changed, 40 insertions(+), 36 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 363df1d48..da0a87aa3 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -167,6 +167,8 @@ def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] pts.SetPoint( 7, [ bounds[ 0 ], bounds[ 3 ], bounds[ 5 ] ] ) return pts + # def _translateToNegQuadrant(self, ) + def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: """Get the source and target points for the transformation. @@ -177,6 +179,12 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: tuple[vtkPoints, vtkPoints]: source and target points for the transformation. """ pts = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) + #translate pts so they always lie on the -z,-y,-x quadrant + off = np.asarray( [ + -2 * np.amax( np.abs( pts[ :, 0 ] ) ), -2 * np.amax( np.abs( pts[ :, 1 ] ) ), + -2 * np.amax( np.abs( pts[ :, 2 ] ) ) + ] ) + pts += off further_ix = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin org = pts[ further_ix, : ] if self.userTranslation is not None: @@ -188,21 +196,29 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: dist_indexes = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) # find u,v,w v1 = pts[ dist_indexes[ 1 ], : ] - org - v1 /= np.linalg.norm( v1 ) v2 = pts[ dist_indexes[ 2 ], : ] - org + v1 /= np.linalg.norm( v1 ) v2 /= np.linalg.norm( v2 ) + if np.abs( v1[ 0 ] ) > np.abs( v2[ 0 ] ): + v1, v2 = v2, v1 + # ensure orthogonality v2 -= np.dot( v2, v1 ) * v1 v2 /= np.linalg.norm( v2 ) v3 = np.cross( v1, v2 ) v3 /= np.linalg.norm( v3 ) + #reorder axis if v3 points downward + if v3[ 2 ] < 0: + v3 = -v3 + v1, v2 = v2, v1 + sourcePts = vtkPoints() sourcePts.SetNumberOfPoints( 4 ) - sourcePts.SetPoint( 0, org ) - sourcePts.SetPoint( 1, v1 + org ) - sourcePts.SetPoint( 2, v2 + org ) - sourcePts.SetPoint( 3, v3 + org ) + sourcePts.SetPoint( 0, org - off ) + sourcePts.SetPoint( 1, v1 + org - off ) + sourcePts.SetPoint( 2, v2 + org - off ) + sourcePts.SetPoint( 3, v3 + org - off ) targetPts = vtkPoints() targetPts.SetNumberOfPoints( 4 ) diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index fe07f9526..9a0f4a972 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -3,11 +3,12 @@ # SPDX-License-Identifier: Apache 2.0 # ruff: noqa: E402 # disable Module level import not at top of file import pytest +import itertools from dataclasses import dataclass -from typing import Generator +from typing import Generator, Tuple from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron -from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader, vtkXMLMultiBlockDataWriter +from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np @@ -24,10 +25,12 @@ class Expected: mesh: vtkUnstructuredGrid -def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int ) -> tuple[ np.ndarray, np.ndarray ]: +def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int, multx: int, multy: int, + multz: int ) -> tuple[ np.ndarray, np.ndarray ]: # np.random.seed(1) # for reproducibility - np.random.default_rng() - off = np.random.randn( 1, 3 ) + # np.random.default_rng() + # off = np.random.randn( 1, 3 ) + off = np.max( [ Lx, Ly, Lz ] ) * np.asarray( [ [ multx, multy, multz ] ] ) pts = [] x, y, z = np.meshgrid( np.linspace( 0, Lx, nx ), np.linspace( 0, Ly, ny ), np.linspace( 0, Lz, ny ) ) for i in range( x.shape[ 0 ] ): @@ -49,12 +52,14 @@ def __rotate_box( angles: np.ndarray, pts: np.ndarray ) -> np.ndarray: return np.asarray( ( RZ @ RY @ RX @ pts.transpose() ).transpose() ) -def __build_test_mesh() -> Generator[ Expected, None, None ]: +def __build_test_mesh( mxx: Tuple[ int ] ) -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz # np.random.seed(1) # for reproducibility np.random.default_rng() - pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz ) + #test all quadrant + multx, multy, multz = mxx + pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz, multx, multy, multz ) logging.info( f"Offseting of {off}" ) logging.debug( f"Original pts : {pts}" ) @@ -91,7 +96,8 @@ def __build_test_mesh() -> Generator[ Expected, None, None ]: yield Expected( mesh=mesh ) -@pytest.mark.parametrize( "expected", __build_test_mesh() ) +@pytest.mark.parametrize( + "expected", [ item for t in list( itertools.product( [ -1, 1 ], repeat=3 ) ) for item in __build_test_mesh( t ) ] ) def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" ( filter := ClipToMainFrameFilter() ).SetInputData( expected.mesh ) @@ -99,16 +105,18 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: output_mesh = filter.GetOutput() assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() + assert output_mesh.GetBounds()[ 0 ] == pytest.approx( - 0., abs=1e-6 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( - 0., abs=1e-6 ) and output_mesh.GetBounds()[ 4 ] == pytest.approx( 0., abs=1e-6 ) + 0., abs=1e-4 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( 0., abs=1e-4 ) and np.max( + [ np.abs( output_mesh.GetBounds()[ 4 ] ), + np.abs( output_mesh.GetBounds()[ 5 ] ) ] ) == pytest.approx( Lz, abs=1e-4 ) # test diagonal assert np.linalg.norm( np.array( [ output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] - ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-5 ) + ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-4 ) # test aligned with axis v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) v1 = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) @@ -117,16 +125,6 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: assert np.abs( np.dot( v0, v2 ) ) < 1e-10 assert np.abs( np.dot( v1, v2 ) ) < 1e-10 - #check if input has been modified - # w = vtkXMLUnstructuredGridWriter() - # w.SetFileName("./test_rotateAndTranslate_input.vtu") - # w.SetInputData(expected.mesh) - # w.Write() - - # w.SetFileName("./test_rotateAndTranslate.vtu") - # w.SetInputData(output_mesh) - # w.Write() - def test_clipToMainFrame_gen() -> None: """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" @@ -140,13 +138,3 @@ def test_clipToMainFrame_gen() -> None: output_mesh = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == input_mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == input_mesh.GetNumberOfCells() - - # #check if input has been modified - # w = vtkXMLMultiBlockDataWriter() - # w.SetFileName( "./test_rotateAndTranslate_input.vtm" ) - # w.SetInputData( input_mesh ) - # w.Write() - - # w.SetFileName( "./test_rotateAndTranslate.vtm" ) - # w.SetInputData( output_mesh ) - # w.Write() From 6a5d546a51b02ed43694cbd19f9b928ff625c4d4 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 14:44:13 +0200 Subject: [PATCH 18/40] Adressing Romain's comments --- .../geos/mesh/processing/ClipToMainFrame.py | 100 +++++------------- geos-mesh/tests/test_clipToMainFrame.py | 6 +- 2 files changed, 28 insertions(+), 78 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index da0a87aa3..524ac5019 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -3,7 +3,7 @@ # SPDX-FileContributor: Jacques Franc from vtkmodules.numpy_interface import dataset_adapter as dsa -from vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray +from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkFiltersGeneral import vtkOBBTree from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkDataObjectTreeIterator, vtkPolyData @@ -19,37 +19,31 @@ Module to clip a mesh to the main frame using rigid body transformation. Methods include: - - ClipToMainFrame class to compute the transformation - - ClipToMainFrameFilter class to apply the transformation to a mesh + - ClipToMainFrameElement class to compute the transformation + - ClipToMainFrame class to apply the transformation to a mesh To use it: .. code-block:: python - from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameFilter + from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame # Filter inputs. multiBlockDataSet: vtkMultiBlockDataSet - # Optional inputs. - userTranslation: vtkFloatArray - # Instantiate the filter. - filter: ClipToMainFrameFilter = ClipToMainFrameFilter() + filter: ClipToMainFrame = ClipToMainFrame() filter.SetInputData( multiBlockDataSet ) # Do calculations. - if userTranslation is not None: - filter.SetUserTranslation( userTranslation ) filter.Update() output: vtkMultiBlockDataSet = filter.GetOutput() """ -class ClipToMainFrame( vtkLandmarkTransform ): +class ClipToMainFrameElement( vtkLandmarkTransform ): - userTranslation: vtkFloatArray | None = None # user provided translation sourcePts: vtkPoints targetPts: vtkPoints mesh: vtkUnstructuredGrid @@ -71,26 +65,9 @@ def Update( self ) -> None: self.SetModeToRigidBody() super().Update() - def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: - """Set the user translation. - - Args: - userTranslation (vtkFloatArray): user translation to set. - """ - self.userTranslation = userTranslation - - def GetUserTranslation( self ) -> vtkFloatArray | None: - """Get the user translation. - - Returns: - vtkFloatArray | None: user translation if set, None otherwise. - """ - return self.userTranslation - def __str__( self ) -> str: """String representation of the transformation.""" - return super().__str__() + f"\nUser translation: {self.userTranslation}" \ - + f"\nSource points: {self.sourcePts}" \ + return super().__str__() + f"\nSource points: {self.sourcePts}" \ + f"\nTarget points: {self.targetPts}" \ + f"\nAngle-Axis: {self.__getAngleAxis()}" \ + f"\nTranslation: {self._getTranslation()}" @@ -187,9 +164,7 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: pts += off further_ix = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin org = pts[ further_ix, : ] - if self.userTranslation is not None: - org = dsa.numpy_support.vtk_to_numpy( self.userTranslation ) - logging.info( "Using user translation" ) + logging.info( f"Moving point {org} to origin for transformation" ) # find 3 orthogonal vectors # we assume points are on a box @@ -230,48 +205,28 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: return ( sourcePts, targetPts ) -class ClipToMainFrameFilter( vtkTransformFilter ): +class ClipToMainFrame( vtkTransformFilter ): """Filter to clip a mesh to the main frame using ClipToMainFrame class.""" - userTranslation: vtkFloatArray | None = None # user provided translation - - def SetUserTranslation( self, userTranslation: vtkFloatArray ) -> None: - """Set the user translation. - - Args: - userTranslation (vtkFloatArray): user translation to set. - """ - self.userTranslation = userTranslation - - def GetUserTranslation( self ) -> vtkFloatArray | None: - """Get the user translation. - - Returns: - vtkFloatArray | None: user translation if set, None otherwise. - """ - return self.userTranslation - def Update( self ) -> None: """Update the transformation.""" # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - idBlock = self.__locate_reference_point( self.GetInput(), self.userTranslation ) - clip = ClipToMainFrame( self.GetInput().GetBlock( idBlock - 1 ) ) + idBlock = self.__locate_reference_point( self.GetInput()) + clip = ClipToMainFrameElement( self.GetInput().GetDataSet( idBlock ) ) else: - clip = ClipToMainFrame( self.GetInput() ) + clip = ClipToMainFrameElement( self.GetInput() ) - clip.SetUserTranslation( self.userTranslation ) clip.Update() self.SetTransform( clip ) super().Update() - def __locate_reference_point( self, input: vtkMultiBlockDataSet, userTranslation: vtkFloatArray | None ) -> int: + def __locate_reference_point( self, input: vtkMultiBlockDataSet ) -> int: """Locate the block to use as reference for the transformation. Args: input (vtkMultiBlockDataSet): input multiblock mesh. - userTranslation (vtkFloatArray | None): if provided use this translation instead of the computed one. Returns: int: index of the block to use as reference. @@ -292,24 +247,19 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, and pt[ 1 ] <= bounds[ 3 ] and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) #TODO (jacques) make a decorator for this - iter: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() - iter.SetDataSet( input ) - iter.VisitOnlyLeavesOn() - iter.GoToFirstItem() + DOIterator: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() + DOIterator.SetDataSet( input ) + DOIterator.VisitOnlyLeavesOn() + DOIterator.GoToFirstItem() xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( input ) #TODO (jacques) : rewrite with a filter struct - while iter.GetCurrentDataObject() is not None: - block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( iter.GetCurrentDataObject() ) + while DOIterator.GetCurrentDataObject() is not None: + block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( DOIterator.GetCurrentDataObject() ) if block.GetNumberOfPoints() > 0: bounds = block.GetBounds() - if userTranslation is not None: - #check if user translation is inside the block - if __inside( dsa.numpy_support.vtk_to_numpy( self.userTranslation ), bounds ): - logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) - return iter.GetCurrentFlatIndex() - else: - #use the lowest bounds corner as reference point - if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): - logging.info( f"Using block {iter.GetCurrentFlatIndex()} as reference for transformation" ) - return iter.GetCurrentFlatIndex() - iter.GoToNextItem() + + #use the furthest bounds corner as reference point in the all negs quadrant + if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): + logging.info( f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) + return DOIterator.GetCurrentFlatIndex() + DOIterator.GoToNextItem() diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 9a0f4a972..403122b69 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -14,7 +14,7 @@ import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON -from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameFilter +from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame Lx, Ly, Lz = 5, 2, 8 nx, ny, nz = 10, 10, 10 @@ -100,7 +100,7 @@ def __build_test_mesh( mxx: Tuple[ int ] ) -> Generator[ Expected, None, None ]: "expected", [ item for t in list( itertools.product( [ -1, 1 ], repeat=3 ) ) for item in __build_test_mesh( t ) ] ) def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" - ( filter := ClipToMainFrameFilter() ).SetInputData( expected.mesh ) + ( filter := ClipToMainFrame() ).SetInputData( expected.mesh ) filter.Update() output_mesh = filter.GetOutput() assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() @@ -132,7 +132,7 @@ def test_clipToMainFrame_gen() -> None: reader.SetFileName( "./geos-mesh/tests/data/displacedFault.vtm" ) reader.Update() input_mesh = reader.GetOutput() - ( filter := ClipToMainFrameFilter() ).SetInputData( input_mesh ) + ( filter := ClipToMainFrame() ).SetInputData( input_mesh ) filter.Update() print( filter.GetTransform() ) output_mesh = filter.GetOutputDataObject( 0 ) From ce7b9b507d82472515ccdabfde5daf95193046b9 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 14:59:32 +0200 Subject: [PATCH 19/40] Some MaJ + conftest use --- .../geos/mesh/processing/ClipToMainFrame.py | 22 +++++++++---------- geos-mesh/tests/test_clipToMainFrame.py | 17 ++++++-------- 2 files changed, 18 insertions(+), 21 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 524ac5019..80c9a912d 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -52,7 +52,7 @@ def __init__( self, mesh: vtkUnstructuredGrid ) -> None: """Clip mesh to main frame. Args: - mesh (vtkUnstructuredGrid): mesh to clip. + mesh (vtkUnstructuredGrid): Mesh to clip. """ super().__init__() self.mesh = mesh @@ -76,7 +76,7 @@ def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: """Get the angle and axis of the rotation. Returns: - tuple[float, np.ndarray]: angle in degrees and axis of rotation. + tuple[float, np.ndarray]: Angle in degrees and axis of rotation. """ matrix = self.GetMatrix() angle = np.arccos( @@ -94,7 +94,7 @@ def _getTranslation( self ) -> np.ndarray: """Get the translation vector. Returns: - np.ndarray: translation vector. + np.ndarray: The translation vector. """ matrix = self.GetMatrix() return np.array( [ matrix.GetElement( 0, 3 ), matrix.GetElement( 1, 3 ), matrix.GetElement( 2, 3 ) ] ) @@ -103,7 +103,7 @@ def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: """Get the OBB tree of the mesh. Args: - mesh (vtkUnstructuredGrid): mesh to get the OBB tree from. + mesh (vtkUnstructuredGrid): Mesh to get the OBB tree from. Returns: vtkOBBTree: OBB tree of the mesh. @@ -127,7 +127,7 @@ def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] """Get the 8 corners of the bounding box. Args: - bounds (tuple[float, float, float, float, float, float]): bounding box. + bounds (tuple[float, float, float, float, float, float]): Bounding box. Returns: vtkPoints: 8 corners of the bounding box. @@ -150,10 +150,10 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: """Get the source and target points for the transformation. Args: - vpts (vtkPoints): points to transform. + vpts (vtkPoints): Points to transform. Returns: - tuple[vtkPoints, vtkPoints]: source and target points for the transformation. + tuple[vtkPoints, vtkPoints]: Source and target points for the transformation. """ pts = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) #translate pts so they always lie on the -z,-y,-x quadrant @@ -226,18 +226,18 @@ def __locate_reference_point( self, input: vtkMultiBlockDataSet ) -> int: """Locate the block to use as reference for the transformation. Args: - input (vtkMultiBlockDataSet): input multiblock mesh. + input (vtkMultiBlockDataSet): Input multiblock mesh. Returns: - int: index of the block to use as reference. + int: Index of the block to use as reference. """ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, float ] ) -> bool: """Check if a point is inside a bounding box. Args: - pt (np.ndarray): point to check - bounds (tuple[float, float, float, float, float, float]): bounding box. + pt (np.ndarray): Point to check + bounds (tuple[float, float, float, float, float, float]): Bounding box. Returns: bool: True if the point is inside the bounding box, False otherwise. diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 403122b69..a553d2143 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -7,8 +7,7 @@ from dataclasses import dataclass from typing import Generator, Tuple from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron -from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, vtkMultiBlockDataSet import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np @@ -126,15 +125,13 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: assert np.abs( np.dot( v1, v2 ) ) < 1e-10 -def test_clipToMainFrame_gen() -> None: +def test_clipToMainFrame_generic(dataSetTest : vtkMultiBlockDataSet) -> None: """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" - reader = vtkXMLMultiBlockDataReader() - reader.SetFileName( "./geos-mesh/tests/data/displacedFault.vtm" ) - reader.Update() - input_mesh = reader.GetOutput() - ( filter := ClipToMainFrame() ).SetInputData( input_mesh ) + + multiBlockDataSet : vtkMultiBlockDataSet = dataSetTest( "multiblock" ) + ( filter := ClipToMainFrame() ).SetInputData( multiBlockDataSet ) filter.Update() print( filter.GetTransform() ) output_mesh = filter.GetOutputDataObject( 0 ) - assert output_mesh.GetNumberOfPoints() == input_mesh.GetNumberOfPoints() - assert output_mesh.GetNumberOfCells() == input_mesh.GetNumberOfCells() + assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() + assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() From 8439b2481145bccc095cc71b6e29b2570838e257 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 15:25:45 +0200 Subject: [PATCH 20/40] minors --- .../src/geos/mesh/processing/ClipToMainFrame.py | 17 +++++++++-------- geos-mesh/tests/test_clipToMainFrame.py | 3 --- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 80c9a912d..7c004f8f4 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -4,6 +4,7 @@ from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonMath import vtkMatrix4x4 from vtkmodules.vtkFiltersGeneral import vtkOBBTree from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkMultiBlockDataSet, vtkDataObjectTreeIterator, vtkPolyData @@ -78,14 +79,14 @@ def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: Returns: tuple[float, np.ndarray]: Angle in degrees and axis of rotation. """ - matrix = self.GetMatrix() - angle = np.arccos( + matrix : vtkMatrix4x4 = self.GetMatrix() + angle : np.ndarray = np.arccos( ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) if angle == 0: return 0.0, np.array( [ 1.0, 0.0, 0.0 ] ) - rx = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) - ry = matrix.GetElement( 0, 2 ) - matrix.GetElement( 2, 0 ) - rz = matrix.GetElement( 1, 0 ) - matrix.GetElement( 0, 1 ) + rx : float = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) + ry : float = matrix.GetElement( 0, 2 ) - matrix.GetElement( 2, 0 ) + rz : float = matrix.GetElement( 1, 0 ) - matrix.GetElement( 0, 1 ) r = np.array( [ rx, ry, rz ] ) r /= np.linalg.norm( r ) return np.degrees( angle ), r @@ -99,14 +100,14 @@ def _getTranslation( self ) -> np.ndarray: matrix = self.GetMatrix() return np.array( [ matrix.GetElement( 0, 3 ), matrix.GetElement( 1, 3 ), matrix.GetElement( 2, 3 ) ] ) - def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: + def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkPoints: """Get the OBB tree of the mesh. Args: mesh (vtkUnstructuredGrid): Mesh to get the OBB tree from. Returns: - vtkOBBTree: OBB tree of the mesh. + vtkPoints: Points from the 0-level OBB tree of the mesh. Fallback on Axis Aligned Bounding Box """ OBBTree = vtkOBBTree() surfFilter = vtkDataSetSurfaceFilter() @@ -116,7 +117,7 @@ def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkOBBTree: OBBTree.BuildLocator() pdata = vtkPolyData() OBBTree.GenerateRepresentation( 0, pdata ) - # at level 0 this should return 8 corners of the bounding box + # at level 0 this should return 8 corners of the bounding box or fallback on AABB if pdata.GetNumberOfPoints() < 3: logging.warning( "Could not get OBB points, using bounding box points instead" ) return self.__allpoints( mesh.GetBounds() ) diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index a553d2143..ea266df7d 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -26,9 +26,6 @@ class Expected: def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int, multx: int, multy: int, multz: int ) -> tuple[ np.ndarray, np.ndarray ]: - # np.random.seed(1) # for reproducibility - # np.random.default_rng() - # off = np.random.randn( 1, 3 ) off = np.max( [ Lx, Ly, Lz ] ) * np.asarray( [ [ multx, multy, multz ] ] ) pts = [] x, y, z = np.meshgrid( np.linspace( 0, Lx, nx ), np.linspace( 0, Ly, ny ), np.linspace( 0, Lz, ny ) ) From 20eec19bc670718bcf3be1a2561a4d42bb5e3ef6 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 16:29:19 +0200 Subject: [PATCH 21/40] Logger --- .../geos/mesh/processing/ClipToMainFrame.py | 76 +++++++++++++++---- geos-mesh/tests/test_clipToMainFrame.py | 17 ++--- 2 files changed, 68 insertions(+), 25 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 7c004f8f4..81a9bcdfb 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -11,10 +11,12 @@ from vtkmodules.vtkCommonTransforms import vtkLandmarkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter +from geos.utils.Logger import logging, Logger, getLogger + from geos.mesh.utils.genericHelpers import getMultiBlockBounds import numpy as np -import logging +from typing import overload, Any __doc__ = """ Module to clip a mesh to the main frame using rigid body transformation. @@ -31,12 +33,19 @@ # Filter inputs. multiBlockDataSet: vtkMultiBlockDataSet + # Optional Inputs + speHandler : bool # Instantiate the filter. filter: ClipToMainFrame = ClipToMainFrame() filter.SetInputData( multiBlockDataSet ) + # Set the handler of yours (only if speHandler is True). + yourHandler: logging.Handler + filter.setLoggerHandler( yourHandler ) + # Do calculations. + filter.ComputeTransform() filter.Update() output: vtkMultiBlockDataSet = filter.GetOutput() @@ -79,14 +88,14 @@ def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: Returns: tuple[float, np.ndarray]: Angle in degrees and axis of rotation. """ - matrix : vtkMatrix4x4 = self.GetMatrix() - angle : np.ndarray = np.arccos( + matrix: vtkMatrix4x4 = self.GetMatrix() + angle: np.ndarray = np.arccos( ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) if angle == 0: return 0.0, np.array( [ 1.0, 0.0, 0.0 ] ) - rx : float = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) - ry : float = matrix.GetElement( 0, 2 ) - matrix.GetElement( 2, 0 ) - rz : float = matrix.GetElement( 1, 0 ) - matrix.GetElement( 0, 1 ) + rx: float = matrix.GetElement( 2, 1 ) - matrix.GetElement( 1, 2 ) + ry: float = matrix.GetElement( 0, 2 ) - matrix.GetElement( 2, 0 ) + rz: float = matrix.GetElement( 1, 0 ) - matrix.GetElement( 0, 1 ) r = np.array( [ rx, ry, rz ] ) r /= np.linalg.norm( r ) return np.degrees( angle ), r @@ -119,7 +128,6 @@ def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkPoints: OBBTree.GenerateRepresentation( 0, pdata ) # at level 0 this should return 8 corners of the bounding box or fallback on AABB if pdata.GetNumberOfPoints() < 3: - logging.warning( "Could not get OBB points, using bounding box points instead" ) return self.__allpoints( mesh.GetBounds() ) return pdata.GetPoints() @@ -166,7 +174,6 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: further_ix = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin org = pts[ further_ix, : ] - logging.info( f"Moving point {org} to origin for transformation" ) # find 3 orthogonal vectors # we assume points are on a box dist_indexes = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) @@ -206,22 +213,60 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: return ( sourcePts, targetPts ) +loggerTitle: str = "Clip mesh to main frame." + + class ClipToMainFrame( vtkTransformFilter ): """Filter to clip a mesh to the main frame using ClipToMainFrame class.""" - def Update( self ) -> None: + def __init__( self, speHandler: bool = False, **properties : str ) -> None: + """Initialize the ClipToMainFrame Filter with optional speHandler args and forwarding properties to main class. + + Args: + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. + properties (**kwargs): kwargs forwarded to vtkTransformFilter. + """ + super().__init__( **properties ) + # Logger. + self.logger: Logger + if not speHandler: + self.logger = getLogger( loggerTitle, True ) + else: + self.logger = logging.getLogger( loggerTitle ) + self.logger.setLevel( logging.INFO ) + + def ComputeTransform( self ) -> None: """Update the transformation.""" # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - idBlock = self.__locate_reference_point( self.GetInput()) + try: + idBlock = self.__locate_reference_point( self.GetInput() ) + except IndexError: + Logger.error("Reference point is not in the domain") + clip = ClipToMainFrameElement( self.GetInput().GetDataSet( idBlock ) ) else: clip = ClipToMainFrameElement( self.GetInput() ) clip.Update() self.SetTransform( clip ) - super().Update() + + def SetLoggerHandler( self, handler: logging.Handler ) -> None: + """Set a specific handler for the filter logger. + + In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels. + + Args: + handler (logging.Handler): The handler to add. + """ + if not self.logger.hasHandlers(): + self.logger.addHandler( handler ) + else: + self.logger.warning( + "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization." + ) def __locate_reference_point( self, input: vtkMultiBlockDataSet ) -> int: """Locate the block to use as reference for the transformation. @@ -243,17 +288,15 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, Returns: bool: True if the point is inside the bounding box, False otherwise. """ - logging.info( f"Checking if point {pt} is inside bounds {bounds}" ) + self.logger.info( f"Checking if point {pt} is inside bounds {bounds}" ) return ( pt[ 0 ] >= bounds[ 0 ] and pt[ 0 ] <= bounds[ 1 ] and pt[ 1 ] >= bounds[ 2 ] and pt[ 1 ] <= bounds[ 3 ] and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) - #TODO (jacques) make a decorator for this DOIterator: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() DOIterator.SetDataSet( input ) DOIterator.VisitOnlyLeavesOn() DOIterator.GoToFirstItem() xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( input ) - #TODO (jacques) : rewrite with a filter struct while DOIterator.GetCurrentDataObject() is not None: block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( DOIterator.GetCurrentDataObject() ) if block.GetNumberOfPoints() > 0: @@ -261,6 +304,9 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, #use the furthest bounds corner as reference point in the all negs quadrant if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): - logging.info( f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) + self.logger.info( + f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) return DOIterator.GetCurrentFlatIndex() DOIterator.GoToNextItem() + + raise IndexError diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index ea266df7d..b6b6a12ab 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -2,13 +2,13 @@ # SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-License-Identifier: Apache 2.0 # ruff: noqa: E402 # disable Module level import not at top of file +# mypy: disable-error-code="operator" import pytest import itertools from dataclasses import dataclass from typing import Generator, Tuple from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, vtkMultiBlockDataSet -import logging # for debug from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON @@ -48,7 +48,7 @@ def __rotate_box( angles: np.ndarray, pts: np.ndarray ) -> np.ndarray: return np.asarray( ( RZ @ RY @ RX @ pts.transpose() ).transpose() ) -def __build_test_mesh( mxx: Tuple[ int ] ) -> Generator[ Expected, None, None ]: +def __build_test_mesh( mxx: Tuple[ int, ...] ) -> Generator[ Expected, None, None ]: # generate random points in a box Lx, Ly, Lz # np.random.seed(1) # for reproducibility np.random.default_rng() @@ -57,16 +57,11 @@ def __build_test_mesh( mxx: Tuple[ int ] ) -> Generator[ Expected, None, None ]: multx, multy, multz = mxx pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz, multx, multy, multz ) - logging.info( f"Offseting of {off}" ) - logging.debug( f"Original pts : {pts}" ) angles = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad - logging.info( f"angles {angles[0]}" ) pts = __rotate_box( angles[ 0 ], pts ) - logging.debug( f"Rotated pts : {pts}" ) pts[ :, 0 ] += off[ 0 ][ 0 ] pts[ :, 1 ] += off[ 0 ][ 1 ] pts[ :, 2 ] += off[ 0 ][ 2 ] - logging.debug( f"Translated pts : {pts}" ) # Creating multiple meshes, each time with a different angles mesh = vtkUnstructuredGrid() @@ -92,11 +87,13 @@ def __build_test_mesh( mxx: Tuple[ int ] ) -> Generator[ Expected, None, None ]: yield Expected( mesh=mesh ) +# arg-type: ignore @pytest.mark.parametrize( "expected", [ item for t in list( itertools.product( [ -1, 1 ], repeat=3 ) ) for item in __build_test_mesh( t ) ] ) def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" ( filter := ClipToMainFrame() ).SetInputData( expected.mesh ) + filter.ComputeTransform() filter.Update() output_mesh = filter.GetOutput() assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() @@ -122,11 +119,11 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: assert np.abs( np.dot( v1, v2 ) ) < 1e-10 -def test_clipToMainFrame_generic(dataSetTest : vtkMultiBlockDataSet) -> None: +def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" - - multiBlockDataSet : vtkMultiBlockDataSet = dataSetTest( "multiblock" ) + multiBlockDataSet: vtkMultiBlockDataSet = dataSetTest( "multiblock" ) ( filter := ClipToMainFrame() ).SetInputData( multiBlockDataSet ) + filter.ComputeTransform() filter.Update() print( filter.GetTransform() ) output_mesh = filter.GetOutputDataObject( 0 ) From 633e8102f5d9eb0f92d514e06d526ce215bce223 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 16:37:58 +0200 Subject: [PATCH 22/40] placeholders --- .../src/geos/pv/plugins/PVClipToMainFrame.py | 50 +++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py new file mode 100644 index 000000000..ffc472689 --- /dev/null +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -0,0 +1,50 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Jacques Franc +# ruff: noqa: E402 # disable Module level import not at top of file +import sys +from pathlib import Path + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/util/vtkAlgorithm.py +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py + +from vtkmodules.vtkCommonDataModel import ( + vtkMultiBlockDataSet, ) + +from vtkmodules.vtkCommonCore import ( + vtkInformation, + vtkInformationVector, +) + +# update sys.path to load all GEOS Python Package dependencies +geos_pv_path: Path = Path( __file__ ).parent.parent.parent.parent.parent +sys.path.insert( 0, str( geos_pv_path / "src" ) ) +from geos.pv.utils.config import update_paths + +update_paths() + +from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame + +__doc__ = """ +Clip the input mesh to the main frame applying the correct LandmarkTransform + +To use it: + +* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVClipToMainFrame. +* Apply. + +""" + +@smproxy.filter( name="PVFillPartialArrays", label="Fill Partial Arrays" ) +@smhint.xml( '' ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( + dataTypes=[ "vtkMultiBlockDataSet" ], + composite_data_supported=True, +) +class PVClipToMainFrame( VTKPythonAlgorithmBase ): + pass \ No newline at end of file From a75723f06e83d266948eb67cb19cf6e67f3810b6 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 16:54:52 +0200 Subject: [PATCH 23/40] yapf --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 81a9bcdfb..5c6f7b46c 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -16,7 +16,6 @@ from geos.mesh.utils.genericHelpers import getMultiBlockBounds import numpy as np -from typing import overload, Any __doc__ = """ Module to clip a mesh to the main frame using rigid body transformation. @@ -219,7 +218,7 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: class ClipToMainFrame( vtkTransformFilter ): """Filter to clip a mesh to the main frame using ClipToMainFrame class.""" - def __init__( self, speHandler: bool = False, **properties : str ) -> None: + def __init__( self, speHandler: bool = False, **properties: str ) -> None: """Initialize the ClipToMainFrame Filter with optional speHandler args and forwarding properties to main class. Args: @@ -244,8 +243,8 @@ def ComputeTransform( self ) -> None: try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: - Logger.error("Reference point is not in the domain") - + self.logger.error( "Reference point is not in the domain" ) + clip = ClipToMainFrameElement( self.GetInput().GetDataSet( idBlock ) ) else: clip = ClipToMainFrameElement( self.GetInput() ) @@ -308,5 +307,5 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) return DOIterator.GetCurrentFlatIndex() DOIterator.GoToNextItem() - + raise IndexError From 7f7fea638bbb762ae75c96b63b8b9527d0f3d4b7 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 12 Sep 2025 17:07:09 +0200 Subject: [PATCH 24/40] stars not good for docs --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 5c6f7b46c..b724c4b85 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -224,7 +224,7 @@ def __init__( self, speHandler: bool = False, **properties: str ) -> None: Args: speHandler (bool, optional): True to use a specific handler, False to use the internal handler. Defaults to False. - properties (**kwargs): kwargs forwarded to vtkTransformFilter. + properties (kwargs): kwargs forwarded to vtkTransformFilter. """ super().__init__( **properties ) # Logger. From 3b27c66836e001f809fa9dad997a0ac989e0e232 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 15 Sep 2025 13:56:11 +0200 Subject: [PATCH 25/40] first working version for input: VTU --- .../src/geos/pv/plugins/PVClipToMainFrame.py | 47 +++++++++++++++++-- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index ffc472689..99ab3f729 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -13,7 +13,9 @@ ) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from vtkmodules.vtkCommonDataModel import ( - vtkMultiBlockDataSet, ) + vtkMultiBlockDataSet, + vtkUnstructuredGrid, + ) from vtkmodules.vtkCommonCore import ( vtkInformation, @@ -27,7 +29,7 @@ update_paths() -from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame +from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameElement, ClipToMainFrame __doc__ = """ Clip the input mesh to the main frame applying the correct LandmarkTransform @@ -39,12 +41,47 @@ """ -@smproxy.filter( name="PVFillPartialArrays", label="Fill Partial Arrays" ) +@smproxy.filter( name="PVClipToMainFrame", label="Clip to the main frame" ) @smhint.xml( '' ) @smproperty.input( name="Input", port_index=0 ) @smdomain.datatype( - dataTypes=[ "vtkMultiBlockDataSet" ], + dataTypes=[ "vtkMultiBlockDataSet", "vtkUnstructuredGrid" ], composite_data_supported=True, ) class PVClipToMainFrame( VTKPythonAlgorithmBase ): - pass \ No newline at end of file + + def __init__(self): + VTKPythonAlgorithmBase.__init__(self, + nInputPorts=1, + nOutputPorts=1, + outputType='vtkMultiBlockDataSet') + + self.__realFilter = ClipToMainFrame() + if not self.__realFilter.logger.hasHandlers(): + self.__realFilter.setLoggerHandler( VTKHandler() ) + + def RequestData(self, request, inInfo, outInfo) -> int: + inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData(inInfo,0,0) + outputMesh : vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData(outInfo,0) + + # struct + def logInfos(obj, logger): + logger.info(f"outputMesh has {obj.GetNumberOfPoints()} points") + logger.info(f"outputMesh has {obj.GetNumberOfCells()} cells") + logger.info(f"output type {obj.GetClassName()}") + + self.__realFilter.SetInputData(inputMesh) + self.__realFilter.logger.info(f"inInfo has {inInfo}") + logInfos(inputMesh, self.__realFilter.logger) + + self.__realFilter.ComputeTransform() + self.__realFilter.Update() + + outputMesh.SetBlock(0, self.__realFilter.GetOutput() ) + self.__realFilter.logger.info(f"outInfo has {outInfo}") + logInfos(outputMesh, self.__realFilter.logger) + + + return 1 + + From e30a4a4446088414efcf45552934ee51287d472d Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 15 Sep 2025 17:20:10 +0200 Subject: [PATCH 26/40] first working version with MBDS --- .../src/geos/pv/plugins/PVClipToMainFrame.py | 34 +++++++++++-------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index 99ab3f729..c1b37faef 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -29,7 +29,7 @@ update_paths() -from geos.mesh.processing.ClipToMainFrame import ClipToMainFrameElement, ClipToMainFrame +from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame __doc__ = """ Clip the input mesh to the main frame applying the correct LandmarkTransform @@ -45,7 +45,7 @@ @smhint.xml( '' ) @smproperty.input( name="Input", port_index=0 ) @smdomain.datatype( - dataTypes=[ "vtkMultiBlockDataSet", "vtkUnstructuredGrid" ], + dataTypes=[ "vtkMultiBlockDataSet"], composite_data_supported=True, ) class PVClipToMainFrame( VTKPythonAlgorithmBase ): @@ -54,33 +54,37 @@ def __init__(self): VTKPythonAlgorithmBase.__init__(self, nInputPorts=1, nOutputPorts=1, - outputType='vtkMultiBlockDataSet') + inputType="vtkMultiBlockDataSet", + outputType="vtkMultiBlockDataSet" + ) self.__realFilter = ClipToMainFrame() if not self.__realFilter.logger.hasHandlers(): self.__realFilter.setLoggerHandler( VTKHandler() ) - def RequestData(self, request, inInfo, outInfo) -> int: + def RequestData(self, request : vtkInformation, inInfo: list[vtkInformationVector], + outInfo:vtkInformationVector) -> int: inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData(inInfo,0,0) outputMesh : vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData(outInfo,0) - + + # outputMesh.ShallowCopy(inputMesh) # struct - def logInfos(obj, logger): - logger.info(f"outputMesh has {obj.GetNumberOfPoints()} points") - logger.info(f"outputMesh has {obj.GetNumberOfCells()} cells") - logger.info(f"output type {obj.GetClassName()}") + self.__realFilter.logger.info(f"size Info {outInfo.GetNumberOfInformationObjects()}") + def logInfos(obj, logger, header): + logger.info(f"--{header}--") + logger.info(f"mesh has {obj.GetNumberOfPoints()} points") + logger.info(f"mesh has {obj.GetNumberOfCells()} cells") + logger.info(f"type {obj.GetClassName()}") self.__realFilter.SetInputData(inputMesh) - self.__realFilter.logger.info(f"inInfo has {inInfo}") - logInfos(inputMesh, self.__realFilter.logger) + logInfos(inputMesh, self.__realFilter.logger,"input") self.__realFilter.ComputeTransform() self.__realFilter.Update() + + outputMesh.ShallowCopy( self.__realFilter.GetOutputDataObject(0) ) - outputMesh.SetBlock(0, self.__realFilter.GetOutput() ) - self.__realFilter.logger.info(f"outInfo has {outInfo}") - logInfos(outputMesh, self.__realFilter.logger) - + logInfos(outputMesh, self.__realFilter.logger,"output") return 1 From b8ffc41543090da80fcd7f7f7d7565a98430c2de Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 16 Sep 2025 11:32:17 +0200 Subject: [PATCH 27/40] last corr before PR --- .../src/geos/mesh/processing/ClipToMainFrame.py | 4 +++- geos-mesh/tests/test_clipToMainFrame.py | 1 + geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py | 13 ------------- 3 files changed, 4 insertions(+), 14 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 81a9bcdfb..b58a6ce02 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -241,13 +241,15 @@ def ComputeTransform( self ) -> None: # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point + self.logger.info(f"Processing MultiblockDataSet") try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: - Logger.error("Reference point is not in the domain") + self.logger.error("Reference point is not in the domain") clip = ClipToMainFrameElement( self.GetInput().GetDataSet( idBlock ) ) else: + self.logger.info(f"Processing untructuredGrid") clip = ClipToMainFrameElement( self.GetInput() ) clip.Update() diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index b6b6a12ab..2819d6db4 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -129,3 +129,4 @@ def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: output_mesh = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() + assert output_mesh.IsA( 'vtkMultiBlockDataSet' ) diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index c1b37faef..c7bdf0a7b 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -67,25 +67,12 @@ def RequestData(self, request : vtkInformation, inInfo: list[vtkInformationVecto inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData(inInfo,0,0) outputMesh : vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData(outInfo,0) - # outputMesh.ShallowCopy(inputMesh) # struct - self.__realFilter.logger.info(f"size Info {outInfo.GetNumberOfInformationObjects()}") - def logInfos(obj, logger, header): - logger.info(f"--{header}--") - logger.info(f"mesh has {obj.GetNumberOfPoints()} points") - logger.info(f"mesh has {obj.GetNumberOfCells()} cells") - logger.info(f"type {obj.GetClassName()}") - self.__realFilter.SetInputData(inputMesh) - logInfos(inputMesh, self.__realFilter.logger,"input") - self.__realFilter.ComputeTransform() self.__realFilter.Update() - outputMesh.ShallowCopy( self.__realFilter.GetOutputDataObject(0) ) - logInfos(outputMesh, self.__realFilter.logger,"output") - return 1 From af9f31f036176c222b0db39b415eca13750af949 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 17 Sep 2025 08:56:41 +0200 Subject: [PATCH 28/40] apply Romain's comments --- .../geos/mesh/processing/ClipToMainFrame.py | 74 +++++++++---------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index b724c4b85..806494c73 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -16,6 +16,9 @@ from geos.mesh.utils.genericHelpers import getMultiBlockBounds import numpy as np +import numpy.typing as npt + +from typing import Tuple __doc__ = """ Module to clip a mesh to the main frame using rigid body transformation. @@ -79,16 +82,15 @@ def __str__( self ) -> str: return super().__str__() + f"\nSource points: {self.sourcePts}" \ + f"\nTarget points: {self.targetPts}" \ + f"\nAngle-Axis: {self.__getAngleAxis()}" \ - + f"\nTranslation: {self._getTranslation()}" - - def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: + + f"\nTranslation: {self.__getTranslation()}" + + def __getAngleAxis( self ) -> tuple[ float, npt.NDArray[ np.double ] ]: """Get the angle and axis of the rotation. - Returns: - tuple[float, np.ndarray]: Angle in degrees and axis of rotation. + tuple[float, npt.NDArray[np.double]]: Angle in degrees and axis of rotation. """ matrix: vtkMatrix4x4 = self.GetMatrix() - angle: np.ndarray = np.arccos( + angle: np.double = np.arccos( ( matrix.GetElement( 0, 0 ) + matrix.GetElement( 1, 1 ) + matrix.GetElement( 2, 2 ) - 1 ) / 2 ) if angle == 0: return 0.0, np.array( [ 1.0, 0.0, 0.0 ] ) @@ -99,13 +101,13 @@ def __getAngleAxis( self ) -> tuple[ float, np.ndarray ]: r /= np.linalg.norm( r ) return np.degrees( angle ), r - def _getTranslation( self ) -> np.ndarray: + def __getTranslation( self ) -> npt.NDArray[ np.double ]: """Get the translation vector. Returns: - np.ndarray: The translation vector. + npt.NDArray[ np.double ]: The translation vector. """ - matrix = self.GetMatrix() + matrix: vtkMatrix4x4 = self.GetMatrix() return np.array( [ matrix.GetElement( 0, 3 ), matrix.GetElement( 1, 3 ), matrix.GetElement( 2, 3 ) ] ) def __getOBBTree( self, mesh: vtkUnstructuredGrid ) -> vtkPoints: @@ -152,7 +154,6 @@ def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] pts.SetPoint( 7, [ bounds[ 0 ], bounds[ 3 ], bounds[ 5 ] ] ) return pts - # def _translateToNegQuadrant(self, ) def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: """Get the source and target points for the transformation. @@ -163,22 +164,22 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: Returns: tuple[vtkPoints, vtkPoints]: Source and target points for the transformation. """ - pts = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) + pts : npt.NDArray[ np.double ] = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) #translate pts so they always lie on the -z,-y,-x quadrant - off = np.asarray( [ + off : npt.NDArray[ np.double ] = np.asarray( [ -2 * np.amax( np.abs( pts[ :, 0 ] ) ), -2 * np.amax( np.abs( pts[ :, 1 ] ) ), -2 * np.amax( np.abs( pts[ :, 2 ] ) ) ] ) pts += off - further_ix = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin - org = pts[ further_ix, : ] + further_ix : np.int_ = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin + org : npt.NDArray = pts[ further_ix, : ] # find 3 orthogonal vectors # we assume points are on a box - dist_indexes = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) + dist_indexes : npt.NDArray[np.double] = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) # find u,v,w - v1 = pts[ dist_indexes[ 1 ], : ] - org - v2 = pts[ dist_indexes[ 2 ], : ] - org + v1 : npt.NDArray[np.double] = pts[ dist_indexes[ 1 ], : ] - org + v2 : npt.NDArray[np.double] = pts[ dist_indexes[ 2 ], : ] - org v1 /= np.linalg.norm( v1 ) v2 /= np.linalg.norm( v2 ) if np.abs( v1[ 0 ] ) > np.abs( v2[ 0 ] ): @@ -187,7 +188,7 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: # ensure orthogonality v2 -= np.dot( v2, v1 ) * v1 v2 /= np.linalg.norm( v2 ) - v3 = np.cross( v1, v2 ) + v3 : npt.NDArray[np.double] = np.cross( v1, v2 ) v3 /= np.linalg.norm( v3 ) #reorder axis if v3 points downward @@ -222,9 +223,9 @@ def __init__( self, speHandler: bool = False, **properties: str ) -> None: """Initialize the ClipToMainFrame Filter with optional speHandler args and forwarding properties to main class. Args: - speHandler (bool, optional): True to use a specific handler, False to use the internal handler. - Defaults to False. - properties (kwargs): kwargs forwarded to vtkTransformFilter. + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. + properties (kwargs): kwargs forwarded to vtkTransformFilter. """ super().__init__( **properties ) # Logger. @@ -267,7 +268,7 @@ def SetLoggerHandler( self, handler: logging.Handler ) -> None: "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization." ) - def __locate_reference_point( self, input: vtkMultiBlockDataSet ) -> int: + def __locate_reference_point( self, multiBlockDataSet: vtkMultiBlockDataSet ) -> int: """Locate the block to use as reference for the transformation. Args: @@ -277,11 +278,10 @@ def __locate_reference_point( self, input: vtkMultiBlockDataSet ) -> int: int: Index of the block to use as reference. """ - def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, float ] ) -> bool: - """Check if a point is inside a bounding box. - + def __inside( pt: npt.NDArray[ np.double], bounds: tuple[ float, float, float, float, float, float ] ) -> bool: + """ Args: - pt (np.ndarray): Point to check + pt (npt.NDArray[np.double]): Point to check. bounds (tuple[float, float, float, float, float, float]): Bounding box. Returns: @@ -292,20 +292,20 @@ def __inside( pt: np.ndarray, bounds: tuple[ float, float, float, float, float, and pt[ 1 ] <= bounds[ 3 ] and pt[ 2 ] >= bounds[ 4 ] and pt[ 2 ] <= bounds[ 5 ] ) DOIterator: vtkDataObjectTreeIterator = vtkDataObjectTreeIterator() - DOIterator.SetDataSet( input ) + DOIterator.SetDataSet( multiBlockDataSet ) DOIterator.VisitOnlyLeavesOn() DOIterator.GoToFirstItem() - xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( input ) + xmin : float + ymin : float + zmin : float + xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( multiBlockDataSet ) while DOIterator.GetCurrentDataObject() is not None: - block: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( DOIterator.GetCurrentDataObject() ) - if block.GetNumberOfPoints() > 0: - bounds = block.GetBounds() - - #use the furthest bounds corner as reference point in the all negs quadrant - if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): - self.logger.info( - f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) - return DOIterator.GetCurrentFlatIndex() + dataSet: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( DOIterator.GetCurrentDataObject() ) + bounds: Tuple[ float, float, float, float, float, float ] = dataSet.GetBounds() + #use the furthest bounds corner as reference point in the all negs quadrant + if __inside( np.asarray( [ xmin, ymin, zmin ] ), bounds ): + self.logger.info( f"Using block {DOIterator.GetCurrentFlatIndex()} as reference for transformation" ) + return DOIterator.GetCurrentFlatIndex() DOIterator.GoToNextItem() raise IndexError From 4cdf4d96691c46a1ace5af22e7ee77d584c53d49 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 17 Sep 2025 10:28:12 +0200 Subject: [PATCH 29/40] fin comm test + mypy/ruff/yapf --- .../geos/mesh/processing/ClipToMainFrame.py | 47 +++++++++---------- geos-mesh/tests/test_clipToMainFrame.py | 39 +++++++++------ 2 files changed, 48 insertions(+), 38 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 806494c73..ecc099701 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -83,11 +83,11 @@ def __str__( self ) -> str: + f"\nTarget points: {self.targetPts}" \ + f"\nAngle-Axis: {self.__getAngleAxis()}" \ + f"\nTranslation: {self.__getTranslation()}" - + def __getAngleAxis( self ) -> tuple[ float, npt.NDArray[ np.double ] ]: """Get the angle and axis of the rotation. - tuple[float, npt.NDArray[np.double]]: Angle in degrees and axis of rotation. + tuple[float, npt.NDArray[np.double]]: Angle in degrees and axis of rotation. """ matrix: vtkMatrix4x4 = self.GetMatrix() angle: np.double = np.arccos( @@ -154,7 +154,6 @@ def __allpoints( self, bounds: tuple[ float, float, float, float, float, float ] pts.SetPoint( 7, [ bounds[ 0 ], bounds[ 3 ], bounds[ 5 ] ] ) return pts - def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: """Get the source and target points for the transformation. @@ -164,22 +163,23 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: Returns: tuple[vtkPoints, vtkPoints]: Source and target points for the transformation. """ - pts : npt.NDArray[ np.double ] = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) + pts: npt.NDArray[ np.double ] = dsa.numpy_support.vtk_to_numpy( vpts.GetData() ) #translate pts so they always lie on the -z,-y,-x quadrant - off : npt.NDArray[ np.double ] = np.asarray( [ + off: npt.NDArray[ np.double ] = np.asarray( [ -2 * np.amax( np.abs( pts[ :, 0 ] ) ), -2 * np.amax( np.abs( pts[ :, 1 ] ) ), -2 * np.amax( np.abs( pts[ :, 2 ] ) ) ] ) pts += off - further_ix : np.int_ = np.argmax( np.linalg.norm( pts, axis=1 ) ) # by default take the min point furthest from origin - org : npt.NDArray = pts[ further_ix, : ] + further_ix: np.int_ = np.argmax( np.linalg.norm( + pts, axis=1 ) ) # by default take the min point furthest from origin + org: npt.NDArray = pts[ further_ix, : ] # find 3 orthogonal vectors # we assume points are on a box - dist_indexes : npt.NDArray[np.double] = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) + dist_indexes: npt.NDArray[ np.int_ ] = np.argsort( np.linalg.norm( pts - org, axis=1 ) ) # find u,v,w - v1 : npt.NDArray[np.double] = pts[ dist_indexes[ 1 ], : ] - org - v2 : npt.NDArray[np.double] = pts[ dist_indexes[ 2 ], : ] - org + v1: npt.NDArray[ np.double ] = pts[ dist_indexes[ 1 ], : ] - org + v2: npt.NDArray[ np.double ] = pts[ dist_indexes[ 2 ], : ] - org v1 /= np.linalg.norm( v1 ) v2 /= np.linalg.norm( v2 ) if np.abs( v1[ 0 ] ) > np.abs( v2[ 0 ] ): @@ -188,7 +188,7 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: # ensure orthogonality v2 -= np.dot( v2, v1 ) * v1 v2 /= np.linalg.norm( v2 ) - v3 : npt.NDArray[np.double] = np.cross( v1, v2 ) + v3: npt.NDArray[ np.double ] = np.cross( v1, v2 ) v3 /= np.linalg.norm( v3 ) #reorder axis if v3 points downward @@ -198,10 +198,10 @@ def __getFramePoints( self, vpts: vtkPoints ) -> tuple[ vtkPoints, vtkPoints ]: sourcePts = vtkPoints() sourcePts.SetNumberOfPoints( 4 ) - sourcePts.SetPoint( 0, org - off ) - sourcePts.SetPoint( 1, v1 + org - off ) - sourcePts.SetPoint( 2, v2 + org - off ) - sourcePts.SetPoint( 3, v3 + org - off ) + sourcePts.SetPoint( 0, list( org - off ) ) + sourcePts.SetPoint( 1, list( v1 + org - off ) ) + sourcePts.SetPoint( 2, list( v2 + org - off ) ) + sourcePts.SetPoint( 3, list( v3 + org - off ) ) targetPts = vtkPoints() targetPts.SetNumberOfPoints( 4 ) @@ -254,9 +254,7 @@ def ComputeTransform( self ) -> None: self.SetTransform( clip ) def SetLoggerHandler( self, handler: logging.Handler ) -> None: - """Set a specific handler for the filter logger. - - In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels. + """Set a specific handler for the filter logger. In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels. Args: handler (logging.Handler): The handler to add. @@ -272,14 +270,15 @@ def __locate_reference_point( self, multiBlockDataSet: vtkMultiBlockDataSet ) -> """Locate the block to use as reference for the transformation. Args: - input (vtkMultiBlockDataSet): Input multiblock mesh. + multiBlockDataSet (vtkMultiBlockDataSet): Input multiblock mesh. Returns: int: Index of the block to use as reference. """ - def __inside( pt: npt.NDArray[ np.double], bounds: tuple[ float, float, float, float, float, float ] ) -> bool: - """ + def __inside( pt: npt.NDArray[ np.double ], bounds: tuple[ float, float, float, float, float, float ] ) -> bool: + """Check if a point is inside a box. + Args: pt (npt.NDArray[np.double]): Point to check. bounds (tuple[float, float, float, float, float, float]): Bounding box. @@ -295,9 +294,9 @@ def __inside( pt: npt.NDArray[ np.double], bounds: tuple[ float, float, float, f DOIterator.SetDataSet( multiBlockDataSet ) DOIterator.VisitOnlyLeavesOn() DOIterator.GoToFirstItem() - xmin : float - ymin : float - zmin : float + xmin: float + ymin: float + zmin: float xmin, _, ymin, _, zmin, _ = getMultiBlockBounds( multiBlockDataSet ) while DOIterator.GetCurrentDataObject() is not None: dataSet: vtkUnstructuredGrid = vtkUnstructuredGrid.SafeDownCast( DOIterator.GetCurrentDataObject() ) diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index b6b6a12ab..fecba2704 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -11,8 +11,11 @@ from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, vtkMultiBlockDataSet from vtkmodules.numpy_interface import dataset_adapter as dsa import numpy as np +import numpy.typing as npt from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON +from typing import List + from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame Lx, Ly, Lz = 5, 2, 8 @@ -25,9 +28,9 @@ class Expected: def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int, multx: int, multy: int, - multz: int ) -> tuple[ np.ndarray, np.ndarray ]: - off = np.max( [ Lx, Ly, Lz ] ) * np.asarray( [ [ multx, multy, multz ] ] ) - pts = [] + multz: int ) -> Tuple[ npt.NDArray[ np.double ], npt.NDArray[ np.double ] ]: + off: npt.NDArray[ np.double ] = np.max( [ Lx, Ly, Lz ] ) * np.asarray( [ [ multx, multy, multz ] ] ) + pts: List[ npt.NDArray[ np.double ] ] = [] x, y, z = np.meshgrid( np.linspace( 0, Lx, nx ), np.linspace( 0, Ly, ny ), np.linspace( 0, Lz, ny ) ) for i in range( x.shape[ 0 ] ): for j in range( x.shape[ 1 ] ): @@ -37,13 +40,16 @@ def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int, multx: int, return ( np.asarray( pts ), off ) -def __rotate_box( angles: np.ndarray, pts: np.ndarray ) -> np.ndarray: - a = angles[ 0 ] - RX = np.asarray( [ [ 1., 0, 0 ], [ 0, np.cos( a ), -np.sin( a ) ], [ 0, np.sin( a ), np.cos( a ) ] ] ) +def __rotate_box( angles: npt.NDArray[ np.double ], pts: npt.NDArray[ np.double ] ) -> npt.NDArray[ np.double ]: + a: np.double = angles[ 0 ] + RX: npt.NDArray[ np.double ] = np.asarray( [ [ 1., 0, 0 ], [ 0, np.cos( a ), -np.sin( a ) ], + [ 0, np.sin( a ), np.cos( a ) ] ] ) a = angles[ 1 ] - RY = np.asarray( [ [ np.cos( a ), 0, np.sin( a ) ], [ 0, 1, 0 ], [ -np.sin( a ), 0, np.cos( a ) ] ] ) + RY: npt.NDArray[ np.double ] = np.asarray( [ [ np.cos( a ), 0, np.sin( a ) ], [ 0, 1, 0 ], + [ -np.sin( a ), 0, np.cos( a ) ] ] ) a = angles[ 2 ] - RZ = np.asarray( [ [ np.cos( a ), -np.sin( a ), 0 ], [ np.sin( a ), np.cos( a ), 0 ], [ 0, 0, 1 ] ] ) + RZ: npt.NDArray[ np.double ] = np.asarray( [ [ np.cos( a ), -np.sin( a ), 0 ], [ np.sin( a ), + np.cos( a ), 0 ], [ 0, 0, 1 ] ] ) return np.asarray( ( RZ @ RY @ RX @ pts.transpose() ).transpose() ) @@ -54,10 +60,15 @@ def __build_test_mesh( mxx: Tuple[ int, ...] ) -> Generator[ Expected, None, Non np.random.default_rng() #test all quadrant + multx: int + multy: int + multz: int multx, multy, multz = mxx + pts: npt.NDArray[ np.double ] + off: npt.NDArray[ np.double ] pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz, multx, multy, multz ) - angles = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad + angles: npt.NDArray[ np.double ] = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad pts = __rotate_box( angles[ 0 ], pts ) pts[ :, 0 ] += off[ 0 ][ 0 ] pts[ :, 1 ] += off[ 0 ][ 1 ] @@ -95,7 +106,7 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: ( filter := ClipToMainFrame() ).SetInputData( expected.mesh ) filter.ComputeTransform() filter.Update() - output_mesh = filter.GetOutput() + output_mesh: vtkUnstructuredGrid = filter.GetOutput() assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() @@ -111,9 +122,9 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-4 ) # test aligned with axis - v0 = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v1 = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v2 = np.array( output_mesh.GetPoint( nx * ny ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v0: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v1: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) + v2: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( nx * ny ) ) - np.array( output_mesh.GetPoint( 0 ) ) assert np.abs( np.dot( v0, v1 ) ) < 1e-10 assert np.abs( np.dot( v0, v2 ) ) < 1e-10 assert np.abs( np.dot( v1, v2 ) ) < 1e-10 @@ -126,6 +137,6 @@ def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: filter.ComputeTransform() filter.Update() print( filter.GetTransform() ) - output_mesh = filter.GetOutputDataObject( 0 ) + output_mesh: vtkMultiBlockDataSet = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() From 3d4d09c5ff611825552644b9ee10fab19a4c7d41 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 17 Sep 2025 10:53:45 +0200 Subject: [PATCH 30/40] wip --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index ecc099701..7bbe904a4 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -271,7 +271,7 @@ def __locate_reference_point( self, multiBlockDataSet: vtkMultiBlockDataSet ) -> Args: multiBlockDataSet (vtkMultiBlockDataSet): Input multiblock mesh. - + Returns: int: Index of the block to use as reference. """ From f6c38343e9cb9fbcde2f64d353dbce0d7bd3e3fd Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 17 Sep 2025 13:43:44 +0200 Subject: [PATCH 31/40] Romain's trick to get both Unstructured and MBDS --- .../geos/mesh/processing/ClipToMainFrame.py | 6 +- .../src/geos/pv/plugins/PVClipToMainFrame.py | 79 +++++++++++++------ 2 files changed, 59 insertions(+), 26 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 45703eb9c..4fe1ff76e 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -241,7 +241,7 @@ def ComputeTransform( self ) -> None: # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - self.logger.info(f"Processing MultiblockDataSet") + self.logger.info( "Processing MultiblockDataSet" ) try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: @@ -249,7 +249,7 @@ def ComputeTransform( self ) -> None: clip = ClipToMainFrameElement( self.GetInput().GetDataSet( idBlock ) ) else: - self.logger.info(f"Processing untructuredGrid") + self.logger.info( "Processing untructuredGrid" ) clip = ClipToMainFrameElement( self.GetInput() ) clip.Update() @@ -273,7 +273,7 @@ def __locate_reference_point( self, multiBlockDataSet: vtkMultiBlockDataSet ) -> Args: multiBlockDataSet (vtkMultiBlockDataSet): Input multiblock mesh. - + Returns: int: Index of the block to use as reference. """ diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index c7bdf0a7b..269039a82 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -7,15 +7,15 @@ from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, -) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/util/vtkAlgorithm.py +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/util/vtkAlgorithm.py from paraview.detail.loghandler import ( # type: ignore[import-not-found] VTKHandler, -) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from vtkmodules.vtkCommonDataModel import ( vtkMultiBlockDataSet, - vtkUnstructuredGrid, - ) + vtkUnstructuredGrid, +) from vtkmodules.vtkCommonCore import ( vtkInformation, @@ -41,38 +41,71 @@ """ + @smproxy.filter( name="PVClipToMainFrame", label="Clip to the main frame" ) @smhint.xml( '' ) @smproperty.input( name="Input", port_index=0 ) @smdomain.datatype( - dataTypes=[ "vtkMultiBlockDataSet"], + dataTypes=[ "vtkMultiBlockDataSet", "vtkUnstructuredGrid" ], composite_data_supported=True, ) class PVClipToMainFrame( VTKPythonAlgorithmBase ): - def __init__(self): - VTKPythonAlgorithmBase.__init__(self, - nInputPorts=1, - nOutputPorts=1, - inputType="vtkMultiBlockDataSet", - outputType="vtkMultiBlockDataSet" - ) - + def __init__( self ) -> None: + """Init motherclass, filter and logger.""" + VTKPythonAlgorithmBase.__init__( self, + nInputPorts=1, + nOutputPorts=1, + inputType="vtkDataObject", + outputType="vtkDataObject" ) + self.__realFilter = ClipToMainFrame() if not self.__realFilter.logger.hasHandlers(): self.__realFilter.setLoggerHandler( VTKHandler() ) - - def RequestData(self, request : vtkInformation, inInfo: list[vtkInformationVector], - outInfo:vtkInformationVector) -> int: - inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData(inInfo,0,0) - outputMesh : vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData(outInfo,0) - + + #ensure I/O consistency + def RequestDataObject( self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inData = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None + if outData is None or ( not outData.IsA( inData.GetClassName() ) ): + outData = inData.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] + + def RequestData( self, request: vtkInformation, inInfo: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. Apply ClipToMainFrame filter. + + Args: + request (vtkInformation): Request + inInfo (list[vtkInformationVector]): Input objects + outInfo (vtkInformationVector): Output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData( inInfo, 0, 0 ) + outputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData( outInfo, 0 ) + + print( inputMesh.GetClassName() ) + # struct - self.__realFilter.SetInputData(inputMesh) + self.__realFilter.SetInputData( inputMesh ) self.__realFilter.ComputeTransform() self.__realFilter.Update() - outputMesh.ShallowCopy( self.__realFilter.GetOutputDataObject(0) ) + outputMesh.ShallowCopy( self.__realFilter.GetOutputDataObject( 0 ) ) + print( outputMesh.GetClassName() ) return 1 - - From 5f65ffaf3617c6613f07cbb3882edb6e64da4b76 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 18 Sep 2025 14:51:28 +0200 Subject: [PATCH 32/40] few fixes --- docs/geos_pv_docs/processing.rst | 6 ++++++ geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 11 +++++------ geos-mesh/tests/test_clipToMainFrame.py | 4 +--- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/docs/geos_pv_docs/processing.rst b/docs/geos_pv_docs/processing.rst index f5fcaf243..c2797cd06 100644 --- a/docs/geos_pv_docs/processing.rst +++ b/docs/geos_pv_docs/processing.rst @@ -18,3 +18,9 @@ PVSplitMesh ---------------------------------- .. automodule:: geos.pv.plugins.PVSplitMesh + + +PVClipToMainFrame +---------------------------------- + +.. automodule:: geos.pv.plugins.PVSClipToMainFrame diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 4fe1ff76e..4c394142d 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -1,6 +1,10 @@ # SPDX-License-Identifier: Apache 2.0 # SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-FileContributor: Jacques Franc +from typing import Tuple + +import numpy as np +import numpy.typing as npt from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkCommonCore import vtkPoints @@ -11,15 +15,10 @@ from vtkmodules.vtkCommonTransforms import vtkLandmarkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter -from geos.utils.Logger import logging, Logger, getLogger +from geos.utils.Logger import logging, Logger, getLogger from geos.mesh.utils.genericHelpers import getMultiBlockBounds -import numpy as np -import numpy.typing as npt - -from typing import Tuple - __doc__ = """ Module to clip a mesh to the main frame using rigid body transformation. diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 31414feb6..40cd3f9f2 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -6,7 +6,7 @@ import pytest import itertools from dataclasses import dataclass -from typing import Generator, Tuple +from typing import Generator, Tuple, List from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, vtkMultiBlockDataSet from vtkmodules.numpy_interface import dataset_adapter as dsa @@ -14,8 +14,6 @@ import numpy.typing as npt from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON -from typing import List - from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame Lx, Ly, Lz = 5, 2, 8 From a470718d09a37cddd92b7cce6439d43c207d67e4 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 18 Sep 2025 15:15:00 +0200 Subject: [PATCH 33/40] typo --- docs/geos_pv_docs/processing.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/geos_pv_docs/processing.rst b/docs/geos_pv_docs/processing.rst index c2797cd06..7655ede61 100644 --- a/docs/geos_pv_docs/processing.rst +++ b/docs/geos_pv_docs/processing.rst @@ -23,4 +23,4 @@ PVSplitMesh PVClipToMainFrame ---------------------------------- -.. automodule:: geos.pv.plugins.PVSClipToMainFrame +.. automodule:: geos.pv.plugins.PVClipToMainFrame From 809991736279d68baa1b727175a085fbab9d94ec Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 18 Sep 2025 15:15:00 +0200 Subject: [PATCH 34/40] infortunate duplication removed --- docs/geos_pv_docs/processing.rst | 2 +- geos-mesh/tests/test_clipToMainFrame.py | 142 +----------------------- 2 files changed, 2 insertions(+), 142 deletions(-) diff --git a/docs/geos_pv_docs/processing.rst b/docs/geos_pv_docs/processing.rst index c2797cd06..7655ede61 100644 --- a/docs/geos_pv_docs/processing.rst +++ b/docs/geos_pv_docs/processing.rst @@ -23,4 +23,4 @@ PVSplitMesh PVClipToMainFrame ---------------------------------- -.. automodule:: geos.pv.plugins.PVSClipToMainFrame +.. automodule:: geos.pv.plugins.PVClipToMainFrame diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index 8db2a71c3..bd019ba40 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -138,144 +138,4 @@ def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: output_mesh: vtkMultiBlockDataSet = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() -# SPDX-FileContributor: Jacques Franc -# SPEDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies -# SPDX-License-Identifier: Apache 2.0 -# ruff: noqa: E402 # disable Module level import not at top of file -# mypy: disable-error-code="operator" -import pytest -import itertools -from dataclasses import dataclass -from typing import Generator, Tuple, List -from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkHexahedron, vtkMultiBlockDataSet -from vtkmodules.numpy_interface import dataset_adapter as dsa -import numpy as np -import numpy.typing as npt -from vtkmodules.util.vtkConstants import VTK_HEXAHEDRON - -from geos.mesh.processing.ClipToMainFrame import ClipToMainFrame - -Lx, Ly, Lz = 5, 2, 8 -nx, ny, nz = 10, 10, 10 - - -@dataclass( frozen=True ) -class Expected: - mesh: vtkUnstructuredGrid - - -def __gen_box( Lx: int, Ly: int, Lz: int, nx: int, ny: int, nz: int, multx: int, multy: int, - multz: int ) -> Tuple[ npt.NDArray[ np.double ], npt.NDArray[ np.double ] ]: - off: npt.NDArray[ np.double ] = np.max( [ Lx, Ly, Lz ] ) * np.asarray( [ [ multx, multy, multz ] ] ) - pts: List[ npt.NDArray[ np.double ] ] = [] - x, y, z = np.meshgrid( np.linspace( 0, Lx, nx ), np.linspace( 0, Ly, ny ), np.linspace( 0, Lz, ny ) ) - for i in range( x.shape[ 0 ] ): - for j in range( x.shape[ 1 ] ): - for k in range( x.shape[ 2 ] ): - pts.append( np.asarray( [ x[ i, j, k ], y[ i, j, k ], z[ i, j, k ] ] ) ) - - return ( np.asarray( pts ), off ) - - -def __rotate_box( angles: npt.NDArray[ np.double ], pts: npt.NDArray[ np.double ] ) -> npt.NDArray[ np.double ]: - a: np.double = angles[ 0 ] - RX: npt.NDArray[ np.double ] = np.asarray( [ [ 1., 0, 0 ], [ 0, np.cos( a ), -np.sin( a ) ], - [ 0, np.sin( a ), np.cos( a ) ] ] ) - a = angles[ 1 ] - RY: npt.NDArray[ np.double ] = np.asarray( [ [ np.cos( a ), 0, np.sin( a ) ], [ 0, 1, 0 ], - [ -np.sin( a ), 0, np.cos( a ) ] ] ) - a = angles[ 2 ] - RZ: npt.NDArray[ np.double ] = np.asarray( [ [ np.cos( a ), -np.sin( a ), 0 ], [ np.sin( a ), - np.cos( a ), 0 ], [ 0, 0, 1 ] ] ) - - return np.asarray( ( RZ @ RY @ RX @ pts.transpose() ).transpose() ) - - -def __build_test_mesh( mxx: Tuple[ int, ...] ) -> Generator[ Expected, None, None ]: - # generate random points in a box Lx, Ly, Lz - # np.random.seed(1) # for reproducibility - np.random.default_rng() - - #test all quadrant - multx: int - multy: int - multz: int - multx, multy, multz = mxx - pts: npt.NDArray[ np.double ] - off: npt.NDArray[ np.double ] - pts, off = __gen_box( Lx, Ly, Lz, nx, ny, nz, multx, multy, multz ) - - angles: npt.NDArray[ np.double ] = -2 * np.pi + np.random.randn( 1, 3 ) * np.pi # random angles in rad - pts = __rotate_box( angles[ 0 ], pts ) - pts[ :, 0 ] += off[ 0 ][ 0 ] - pts[ :, 1 ] += off[ 0 ][ 1 ] - pts[ :, 2 ] += off[ 0 ][ 2 ] - - # Creating multiple meshes, each time with a different angles - mesh = vtkUnstructuredGrid() - ( vtps := vtkPoints() ).SetData( dsa.numpy_support.numpy_to_vtk( pts ) ) - mesh.SetPoints( vtps ) - - ids = vtkIdList() - for i in range( nx - 1 ): - for j in range( ny - 1 ): - for k in range( nz - 1 ): - hex = vtkHexahedron() - ids = hex.GetPointIds() - ids.SetId( 0, i + j * nx + k * nx * ny ) - ids.SetId( 1, i + j * nx + k * nx * ny + 1 ) - ids.SetId( 2, i + j * nx + k * nx * ny + nx * ny + 1 ) - ids.SetId( 3, i + j * nx + k * nx * ny + nx * ny ) - ids.SetId( 4, i + j * nx + k * nx * ny + nx ) - ids.SetId( 5, i + j * nx + k * nx * ny + nx + 1 ) - ids.SetId( 6, i + j * nx + k * nx * ny + nx * ny + nx + 1 ) - ids.SetId( 7, i + j * nx + k * nx * ny + nx * ny + nx ) - mesh.InsertNextCell( VTK_HEXAHEDRON, ids ) - - yield Expected( mesh=mesh ) - - -# arg-type: ignore -@pytest.mark.parametrize( - "expected", [ item for t in list( itertools.product( [ -1, 1 ], repeat=3 ) ) for item in __build_test_mesh( t ) ] ) -def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: - """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" - ( filter := ClipToMainFrame() ).SetInputData( expected.mesh ) - filter.ComputeTransform() - filter.Update() - output_mesh: vtkUnstructuredGrid = filter.GetOutput() - assert output_mesh.GetNumberOfPoints() == expected.mesh.GetNumberOfPoints() - assert output_mesh.GetNumberOfCells() == expected.mesh.GetNumberOfCells() - - assert output_mesh.GetBounds()[ 0 ] == pytest.approx( - 0., abs=1e-4 ) and output_mesh.GetBounds()[ 2 ] == pytest.approx( 0., abs=1e-4 ) and np.max( - [ np.abs( output_mesh.GetBounds()[ 4 ] ), - np.abs( output_mesh.GetBounds()[ 5 ] ) ] ) == pytest.approx( Lz, abs=1e-4 ) - # test diagonal - assert np.linalg.norm( - np.array( [ - output_mesh.GetBounds()[ 1 ] - output_mesh.GetBounds()[ 0 ], - output_mesh.GetBounds()[ 3 ] - output_mesh.GetBounds()[ 2 ], - output_mesh.GetBounds()[ 5 ] - output_mesh.GetBounds()[ 4 ] - ] ) ) == pytest.approx( np.linalg.norm( np.array( [ Lx, Ly, Lz ] ) ), abs=1e-4 ) - # test aligned with axis - v0: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( 1 ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v1: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( nx ) ) - np.array( output_mesh.GetPoint( 0 ) ) - v2: npt.NDArray[ np.double ] = np.array( output_mesh.GetPoint( nx * ny ) ) - np.array( output_mesh.GetPoint( 0 ) ) - assert np.abs( np.dot( v0, v1 ) ) < 1e-10 - assert np.abs( np.dot( v0, v2 ) ) < 1e-10 - assert np.abs( np.dot( v1, v2 ) ) < 1e-10 - - -def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: - """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" - multiBlockDataSet: vtkMultiBlockDataSet = dataSetTest( "multiblock" ) - ( filter := ClipToMainFrame() ).SetInputData( multiBlockDataSet ) - filter.ComputeTransform() - filter.Update() - print( filter.GetTransform() ) - output_mesh: vtkMultiBlockDataSet = filter.GetOutputDataObject( 0 ) - assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() - assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() - assert output_mesh.IsA( 'vtkMultiBlockDataSet' ) + assert output_mesh.IsA('vtkMultiBlockDataSet') From 02598c7572cd484fdc0efb91b7d2066b7ffa0b26 Mon Sep 17 00:00:00 2001 From: Jacques Franc <49998870+jafranc@users.noreply.github.com> Date: Thu, 18 Sep 2025 16:25:13 +0200 Subject: [PATCH 35/40] yapf --- geos-mesh/tests/test_clipToMainFrame.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index bd019ba40..40cd3f9f2 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -138,4 +138,4 @@ def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: output_mesh: vtkMultiBlockDataSet = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() - assert output_mesh.IsA('vtkMultiBlockDataSet') + assert output_mesh.IsA( 'vtkMultiBlockDataSet' ) From a684fb15d0d035171c5ca7f853602ffb6196b6c1 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 19 Sep 2025 11:33:26 +0200 Subject: [PATCH 36/40] Ensure compat w PV12.X and python <3.10 --- geos-mesh/src/geos/mesh/utils/genericHelpers.py | 2 +- geos-mesh/tests/test_clipToMainFrame.py | 13 ++++++++----- geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py | 5 +++-- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/geos-mesh/src/geos/mesh/utils/genericHelpers.py b/geos-mesh/src/geos/mesh/utils/genericHelpers.py index de0624fd9..481b391ff 100644 --- a/geos-mesh/src/geos/mesh/utils/genericHelpers.py +++ b/geos-mesh/src/geos/mesh/utils/genericHelpers.py @@ -36,7 +36,7 @@ def to_vtk_id_list( data: List[ int ] ) -> vtkIdList: return result -def vtk_iter( vtkContainer: vtkIdList | vtkCellTypes ) -> Iterator[ Any ]: +def vtk_iter( vtkContainer: Union[ vtkIdList, vtkCellTypes ] ) -> Iterator[ Any ]: """Utility function transforming a vtk "container" into an iterable. Args: diff --git a/geos-mesh/tests/test_clipToMainFrame.py b/geos-mesh/tests/test_clipToMainFrame.py index bd019ba40..12f2c1940 100644 --- a/geos-mesh/tests/test_clipToMainFrame.py +++ b/geos-mesh/tests/test_clipToMainFrame.py @@ -74,8 +74,9 @@ def __build_test_mesh( mxx: Tuple[ int, ...] ) -> Generator[ Expected, None, Non # Creating multiple meshes, each time with a different angles mesh = vtkUnstructuredGrid() - ( vtps := vtkPoints() ).SetData( dsa.numpy_support.numpy_to_vtk( pts ) ) - mesh.SetPoints( vtps ) + vpts = vtkPoints() + vpts.SetData( dsa.numpy_support.numpy_to_vtk( pts ) ) + mesh.SetPoints( vpts ) ids = vtkIdList() for i in range( nx - 1 ): @@ -101,7 +102,8 @@ def __build_test_mesh( mxx: Tuple[ int, ...] ) -> Generator[ Expected, None, Non "expected", [ item for t in list( itertools.product( [ -1, 1 ], repeat=3 ) ) for item in __build_test_mesh( t ) ] ) def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: """Test the ClipToMainFrameFilter on a rotated and translated box hexa mesh.""" - ( filter := ClipToMainFrame() ).SetInputData( expected.mesh ) + filter = ClipToMainFrame() + filter.SetInputData( expected.mesh ) filter.ComputeTransform() filter.Update() output_mesh: vtkUnstructuredGrid = filter.GetOutput() @@ -131,11 +133,12 @@ def test_clipToMainFrame_polyhedron( expected: Expected ) -> None: def test_clipToMainFrame_generic( dataSetTest: vtkMultiBlockDataSet ) -> None: """Test the ClipToMainFrameFilter on a MultiBlockDataSet.""" multiBlockDataSet: vtkMultiBlockDataSet = dataSetTest( "multiblock" ) - ( filter := ClipToMainFrame() ).SetInputData( multiBlockDataSet ) + filter = ClipToMainFrame() + filter.SetInputData( multiBlockDataSet ) filter.ComputeTransform() filter.Update() print( filter.GetTransform() ) output_mesh: vtkMultiBlockDataSet = filter.GetOutputDataObject( 0 ) assert output_mesh.GetNumberOfPoints() == multiBlockDataSet.GetNumberOfPoints() assert output_mesh.GetNumberOfCells() == multiBlockDataSet.GetNumberOfCells() - assert output_mesh.IsA('vtkMultiBlockDataSet') + assert output_mesh.IsA( 'vtkMultiBlockDataSet' ) diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index 269039a82..f58657387 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -4,6 +4,7 @@ # ruff: noqa: E402 # disable Module level import not at top of file import sys from pathlib import Path +from typing import Union from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, @@ -96,8 +97,8 @@ def RequestData( self, request: vtkInformation, inInfo: list[ vtkInformationVect Returns: int: 1 if calculation successfully ended, 0 otherwise. """ - inputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetInputData( inInfo, 0, 0 ) - outputMesh: vtkMultiBlockDataSet | vtkUnstructuredGrid = self.GetOutputData( outInfo, 0 ) + inputMesh: Union[ vtkMultiBlockDataSet, vtkUnstructuredGrid ] = self.GetInputData( inInfo, 0, 0 ) + outputMesh: Union[ vtkMultiBlockDataSet, vtkUnstructuredGrid ] = self.GetOutputData( outInfo, 0 ) print( inputMesh.GetClassName() ) From 2d2f7620149dbb6509578532734f98bf967e56eb Mon Sep 17 00:00:00 2001 From: jacques franc Date: Fri, 19 Sep 2025 14:42:47 +0200 Subject: [PATCH 37/40] Romain's comments --- .../src/geos/mesh/processing/ClipToMainFrame.py | 3 ++- .../src/geos/pv/plugins/PVClipToMainFrame.py | 17 +++++++---------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index d821a4b78..98881e2cd 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -2,6 +2,7 @@ # SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-FileContributor: Jacques Franc from typing import Tuple +import logging import numpy as np import numpy.typing as npt @@ -15,7 +16,7 @@ from vtkmodules.vtkCommonTransforms import vtkLandmarkTransform from vtkmodules.vtkFiltersGeneral import vtkTransformFilter -from geos.utils.Logger import ( logging, Logger, getLogger ) +from geos.utils.Logger import ( Logger, getLogger ) from geos.mesh.utils.genericHelpers import getMultiBlockBounds __doc__ = """ diff --git a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py index f58657387..8c0cfdb2d 100644 --- a/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py +++ b/geos-pv/src/geos/pv/plugins/PVClipToMainFrame.py @@ -60,9 +60,9 @@ def __init__( self ) -> None: inputType="vtkDataObject", outputType="vtkDataObject" ) - self.__realFilter = ClipToMainFrame() - if not self.__realFilter.logger.hasHandlers(): - self.__realFilter.setLoggerHandler( VTKHandler() ) + self._realFilter = ClipToMainFrame( speHandler=True ) + if not self._realFilter.logger.hasHandlers(): + self._realFilter.SetLoggerHandler( VTKHandler() ) #ensure I/O consistency def RequestDataObject( self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], @@ -100,13 +100,10 @@ def RequestData( self, request: vtkInformation, inInfo: list[ vtkInformationVect inputMesh: Union[ vtkMultiBlockDataSet, vtkUnstructuredGrid ] = self.GetInputData( inInfo, 0, 0 ) outputMesh: Union[ vtkMultiBlockDataSet, vtkUnstructuredGrid ] = self.GetOutputData( outInfo, 0 ) - print( inputMesh.GetClassName() ) - # struct - self.__realFilter.SetInputData( inputMesh ) - self.__realFilter.ComputeTransform() - self.__realFilter.Update() - outputMesh.ShallowCopy( self.__realFilter.GetOutputDataObject( 0 ) ) - print( outputMesh.GetClassName() ) + self._realFilter.SetInputData( inputMesh ) + self._realFilter.ComputeTransform() + self._realFilter.Update() + outputMesh.ShallowCopy( self._realFilter.GetOutputDataObject( 0 ) ) return 1 From 7f5076a98a3353cf3268214a594e49ab9c525697 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 24 Sep 2025 17:21:06 +0200 Subject: [PATCH 38/40] Message from Romain --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 98881e2cd..6390c80a4 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -240,7 +240,7 @@ def ComputeTransform( self ) -> None: # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - self.logger.info( "Processing MultiblockDataSet" ) + self.logger.info( "Processing MultiblockDataSet.\n Please make sure your MultiblockDataSet is owning a vtkUnstructured grid as a leaf." ) try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: From ff8260d2d28984206e4ed040241623e6159eab95 Mon Sep 17 00:00:00 2001 From: Jacques Franc <49998870+jafranc@users.noreply.github.com> Date: Wed, 24 Sep 2025 18:24:44 +0200 Subject: [PATCH 39/40] Update ClipToMainFrame.py --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 6390c80a4..a5cc1b90b 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -240,7 +240,9 @@ def ComputeTransform( self ) -> None: # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - self.logger.info( "Processing MultiblockDataSet.\n Please make sure your MultiblockDataSet is owning a vtkUnstructured grid as a leaf." ) + self.logger.info( + "Processing MultiblockDataSet.\n Please make sure your MultiblockDataSet is owning a vtkUnstructured grid as a leaf." + ) try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: From 7fa7bdeb3f8bba66c458cf80115f2717724d5e47 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 25 Sep 2025 08:32:53 +0200 Subject: [PATCH 40/40] yapfing --- geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py index 6390c80a4..bc24b0349 100644 --- a/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py +++ b/geos-mesh/src/geos/mesh/processing/ClipToMainFrame.py @@ -240,7 +240,9 @@ def ComputeTransform( self ) -> None: # dispatch to ClipToMainFrame depending on input type if isinstance( self.GetInput(), vtkMultiBlockDataSet ): #locate reference point - self.logger.info( "Processing MultiblockDataSet.\n Please make sure your MultiblockDataSet is owning a vtkUnstructured grid as a leaf." ) + self.logger.info( + "Processing MultiblockDataSet.\n Please make sure your MultiblockDataSet is owning a vtkUnstructured grid as a leaf." + ) try: idBlock = self.__locate_reference_point( self.GetInput() ) except IndexError: