From b6ed2c64168b68e9237c2fb87faee4ecc85ac35d Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 27 Nov 2025 14:35:32 +0100 Subject: [PATCH 01/11] start filter diff --- .../generic_processing_tools/DiffFields.py | 115 ++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py new file mode 100644 index 000000000..771b586ab --- /dev/null +++ b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py @@ -0,0 +1,115 @@ +# SPDX-License-Identifier: Apache 2.0 +# SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies +# SPDX-FileContributor: Jacques Franc +import logging + +import numpy as np +import numpy.typing as npt +from functools import partial + +from vtkmodules.numpy_interface import dataset_adapter as dsa +from vtkmodules.vtkFiltersCore import vtkAppendFilter +from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataObject +from geos.utils.Logger import ( Logger, getLogger ) + +__doc__ = """ +Module to diff(operate) on fields between two meshes composing on AppendFilter + +""" + +loggerTitle : str = "DiffFieldsFilter" + +class DiffFieldsFilter: + + def __init__(self): + """ Difference of Fields""" + self.appendFilter = vtkAppendFilter() + self.fieldnames = set() + self.logger = getLogger( loggerTitle ) + + def setMeshes(self, mesh1: vtkDataObject, mesh2: vtkDataObject): + """ Setting meshes to diff.""" + if isinstance(mesh1, vtkMultiBlockDataSet): + #discard partial for now + fieldnames_1 = self._displayFields_mb(mesh1) + else: + fieldnames_1 = self._displayFields(mesh1) + + if isinstance(mesh2, vtkMultiBlockDataSet): + #discard partial for now + fieldnames_2 = self._displayFields_mb(mesh2) + else: + fieldnames_2 = self._displayFields(mesh2) + + self.fieldnames = fieldnames_1.intersection(fieldnames_2) + self.logger.info(f"Fields have been dropped as unique to one of the meshes: {fieldnames_1.difference(fieldnames_2)}") + + #rename fields + cfields = mesh1.GetCellData() + ptable1 = self.getPtable(mesh1) + for i in range(cfields.GetNumberOfArrays()): + cfields.GetArray(i).SetName(f"{cfields.GetArrayName(i)}_0") + self.appendFilter.AddInputData(mesh1) + # + cfields = mesh2.GetCellData() + for i in range(cfields.GetNumberOfArrays()): + cfields.GetArray(i).SetName(f"{cfields.GetArrayName(i)}_1") + #eventually append + self.appendFilter.AddInputData(self._repartition(mesh2, ptable1)) + + + def _display_fields(self,mesh): + cfields = mesh.GetCellData() + namelist = set() + for i in range(0,cfields.GetNumberOfArrays()): + namelist.add(cfields.GetArrayName(i)) + + return namelist + + def display_fields_mb(self, ugrid): + it = ugrid.NewIterator() + namelist = set() + namelist.update(self.display_fields(ugrid.GetDataSet(it))) + return namelist + + def _computeL2(self, f, callback = partial(np.linalg.norm(ord=np.inf))): + s = f.shape + #loop + + self.logger.info("L2 norm for fields {name} is {value}") + pass + + def _computeLmax(self): + #loop + self.logger.info("L2 norm for fields {name} is {value}") + pass + + def _computeL1(self): + + pass + + def _repartition(mesh, ptable): + + pass + + def Update(self): + self.appendFilter.Update() + # size the output numpy array + try: + nv = fields.GetCellData().GetArray("ghostRank").GetNumberOfValues() + except: + nv = 0 + it = self.fields.NewIterator() + while not it.IsDoneWithTraversal(): + nv += self.fields.GetDataSet(it).GetCellData().GetArray("ghostRank").GetNumberOfValues() + it.GoToNextItem() + + for fielname in self.fieldnames: + self.appendFilter.Get + + pass + + def getOutput(self) -> vtkDataObject: + return self.appendFilter.GetOutput() + + \ No newline at end of file From cd1f6f5f76498ae9eb58e7d894f38b38547418c0 Mon Sep 17 00:00:00 2001 From: Jacques Franc <49998870+jafranc@users.noreply.github.com> Date: Mon, 15 Dec 2025 10:24:10 +0000 Subject: [PATCH 02/11] in process transcription to filter --- .../generic_processing_tools/DiffFields.py | 421 ++++++++++++++++-- 1 file changed, 379 insertions(+), 42 deletions(-) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py index 771b586ab..b63273ac9 100644 --- a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py +++ b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py @@ -7,13 +7,14 @@ import numpy.typing as npt from functools import partial +from vtk.util import numpy_support from vtkmodules.numpy_interface import dataset_adapter as dsa from vtkmodules.vtkFiltersCore import vtkAppendFilter from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataObject from geos.utils.Logger import ( Logger, getLogger ) __doc__ = """ -Module to diff(operate) on fields between two meshes composing on AppendFilter +Module to diff(operate) on fields between two meshes composing based on globalIds indirection """ @@ -44,72 +45,408 @@ def setMeshes(self, mesh1: vtkDataObject, mesh2: vtkDataObject): self.fieldnames = fieldnames_1.intersection(fieldnames_2) self.logger.info(f"Fields have been dropped as unique to one of the meshes: {fieldnames_1.difference(fieldnames_2)}") - #rename fields - cfields = mesh1.GetCellData() - ptable1 = self.getPtable(mesh1) - for i in range(cfields.GetNumberOfArrays()): - cfields.GetArray(i).SetName(f"{cfields.GetArrayName(i)}_0") - self.appendFilter.AddInputData(mesh1) - # - cfields = mesh2.GetCellData() - for i in range(cfields.GetNumberOfArrays()): - cfields.GetArray(i).SetName(f"{cfields.GetArrayName(i)}_1") - #eventually append - self.appendFilter.AddInputData(self._repartition(mesh2, ptable1)) - - - def _display_fields(self,mesh): - cfields = mesh.GetCellData() - namelist = set() - for i in range(0,cfields.GetNumberOfArrays()): - namelist.add(cfields.GetArrayName(i)) + +#displays + def _display_cfields(self,fields,namelist): + print("Cell Fields available are :\n") + cfields = fields.GetCellData() + for i in range((off:=len(namelist)),off+cfields.GetNumberOfArrays()): + if cfields.GetArrayName(i - len(namelist)): + print(str(i)+": "+cfields.GetArrayName(i-off)) + namelist.append(cfields.GetArrayName(i)) + else: + print(f"fields {i} is undefined") + + return namelist + + def _display_pfields(self,fields,namelist): + print("Point Fields available are :\n") + pfields = fields.GetPointData() + for i in range(len(namelist),len(namelist) + pfields.GetNumberOfArrays()): + print(str(i)+": *"+pfields.GetArrayName(i)) + namelist.append('*'+pfields.GetArrayName(i)) + + return namelist + + def display_fields(self, fields): + namelist = [] + self._display_pfields(fields,namelist) + self._display_cfields(fields,namelist) return namelist def display_fields_mb(self, ugrid): it = ugrid.NewIterator() - namelist = set() - namelist.update(self.display_fields(ugrid.GetDataSet(it))) + namelist = [] + namelist.extend(self.display_fields(ugrid.GetDataSet(it))) return namelist - def _computeL2(self, f, callback = partial(np.linalg.norm(ord=np.inf))): - s = f.shape - #loop +# - self.logger.info("L2 norm for fields {name} is {value}") - pass - - def _computeLmax(self): + def _computeL(self, f, callback = partial(np.linalg.norm(ord=np.inf))): + """ compute by default inf norm """ + s = f.shape #loop - self.logger.info("L2 norm for fields {name} is {value}") - pass + sp = fp.shape + s = f.shape + for i in range(0,s[1]): + n = callback( f[:,i,i1]-f[:,i,i2]) + print(self.flist[i]+": "+str(n)+" ") + for i in range(0,sp[1]): + n = callback( fp[:,i,i1]-fp[:,i,i2]) + print(self.flist[i]+": "+str(n)+" ") + self.logger.info("Lmax norm for fields {name} is {value}") def _computeL1(self): - - pass + # writer.SetDataModeToAscii() + #mesh = vtk.vtkUnstructuredGrid() + postfix = "" - def _repartition(mesh, ptable): + print(self.flist) + for i,fname in enumerate(self.flist): + try: + arr =numpy_support.numpy_to_vtk(np.abs( f[:,i,i1]-f[:,i,i2] )) + arr.SetName("d"+fname) + self.fields.GetCellData().AddArray(arr) + except: + it = self.fields.NewIterator() + start = 0 + while not it.IsDoneWithTraversal(): + #scalar fill only + if fname[0] == '*': + l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) + arr = numpy_support.numpy_to_vtk(np.abs( fp[l2g,i,i1]-fp[l2g,i,i2] )) + arr.SetName("d"+fname) + self.fields.GetDataSet(it).GetPointData().AddArray(arr) + else: + l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) + arr = numpy_support.numpy_to_vtk(np.abs( f[l2g-nbp,i,i1]-f[l2g-nbp,i,i2] )) + arr.SetName("d"+fname) + self.fields.GetDataSet(it).GetCellData().AddArray(arr) + it.GoToNextItem() - pass - def Update(self): - self.appendFilter.Update() + def _getdata(self): # size the output numpy array + nv =0 + nbp = 0 + npp = 0 + try: nv = fields.GetCellData().GetArray("ghostRank").GetNumberOfValues() + npp = fields.GetPointData().GetArray("ghostRank").GetNumberOfValues() + nbp = np.max( numpy_support.vtk_to_numpy( fields.GetPointData().GetArray("localToGlobalMap") ) ) except: - nv = 0 it = self.fields.NewIterator() while not it.IsDoneWithTraversal(): nv += self.fields.GetDataSet(it).GetCellData().GetArray("ghostRank").GetNumberOfValues() + npp += self.fields.GetDataSet(it).GetPointData().GetArray("ghostRank").GetNumberOfValues() + nbp = np.max( numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) ) it.GoToNextItem() - for fielname in self.fieldnames: - self.appendFilter.Get + cnf = 0 + ncc = 0 + npcc = 0 + for ifields in self.olist: + try: + if ifields[0] == '*': + npcc = fields.GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() + else: + ncc = fields.GetCellData().GetArray(ifields).GetNumberOfComponents() + except: + it = fields.NewIterator() + if ifields[0] == '*': + npcc = fields.GetDataSet(it).GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() + else: + ncc = fields.GetDataSet(it).GetCellData().GetArray(ifields).GetNumberOfComponents() - pass + + ncnf = ncnf + ncc + npcnf = npcnf + npcc + + f = np.zeros(shape=(nv,ncnf,len(filelist)),dtype='float') + fp = np.zeros(shape=(npp,npcnf,len(filelist)),dtype='float') + print(f"nv {nv} ncnf {ncnf} nb {nbp}") + + i = 0 # for file loop + for fileid in filelist: + self.reader.SetFileName(fileid) + print(fileid) + self.reader.Update() + # fields = self.reader.GetOutput() + j = 0 # for field loop + nc = 0 + for nfields in self.olist: + try: + if nfields[0] == '*': + field = self.fields.GetPointData().GetArray(nfields[1:]) + else: + field = self.fields.GetCellData().GetArray(nfields) + + nc = field.GetNumberOfComponents() + if nfields[0] == '*': + f[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(nv,nc) + else: + fp[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(npb,nc) + except: + it = self.fields.NewIterator() + start = 0 + while not it.IsDoneWithTraversal(): + # use localToGlobalMap + if nfields[0] == '*': + field = self.fields.GetDataSet(it).GetPointData().GetArray(nfields[1:]) + nt = field.GetNumberOfValues() + nc = field.GetNumberOfComponents() + l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) + fp[l2g,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) + else: + field = self.fields.GetDataSet(it).GetCellData().GetArray(nfields) + nt = field.GetNumberOfValues() + nc = field.GetNumberOfComponents() + l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) + f[l2g-nbp,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) + it.GoToNextItem() + + if nc>1 & i ==0 : + self.flist[j:j+1] = [ self.flist[j]+"_"+str(k) for k in range(0,nc) ] + j = j + nc + i = i+1 + + return fp,f,nbp def getOutput(self) -> vtkDataObject: - return self.appendFilter.GetOutput() + return self.GetOutput() + +# ##### +# class diff_visu: + +# def __init__(self, fname): +# self.t = fname[-1] +# self.filelist = fname +# print(self.filelist) + +# self.extension = fname[0].split('.')[-1] +# #TODO check extension for all files +# if self.extension == "vtk": +# self.reader = vtk.vtkUnstructuredGridReader() +# self.writer = vtk.vtkUnstructuredGridWriter() +# self.reader.SetFileName(fname[0]) +# self.reader.Update() +# self.fields = self.reader.GetOutput() +# namelist = self.display_fields(self.fields)# as indexes change between time 0 and others +# elif self.extension == "vtr" : +# self.reader = vtk.vtkRectilinearGridReader() +# self.writer = vtk.vtkRectilinearGridWriter() +# self.reader.SetFileName(fname[0]) +# self.reader.Update() +# self.fields = self.reader.GetOutput() +# namelist = self.display_fields(self.fields)# as indexes change between time 0 and others +# elif self.extension == "vtu": +# self.reader = vtk.vtkXMLUnstructuredGridReader() +# self.writer = vtk.vtkXMLUnstructuredGridWriter() +# self.reader.SetFileName(fname[0]) +# self.reader.Update() +# self.fields = self.reader.GetOutput() +# namelist = self.display_fields(self.fields)# as indexes change between time 0 and others +# elif self.extension == "vtm": +# self.reader = vtk.vtkXMLMultiBlockDataReader() +# self.writer = vtk.vtkXMLMultiBlockDataWriter() +# self.reader.SetFileName(fname[0]) +# self.reader.Update() +# self.fields = self.reader.GetOutput() +# namelist = self.display_fields_mb(self.fields)# as indexes change between time 0 and others + +# else: +# raise NotImplementedError + +# #self.data_d = self.len_*[vtk.vtkFloatArray()] +# prs = input("number to diff ?\n") +# # debug +# # prs = '22 23' +# olist = [ namelist[int(item)] for item in prs.split() ] +# print(olist) +# self.flist = olist.copy() + +# fp, f ,nbp = self.extract_data(self.filelist,olist) +# self.write_report(fp,f,0,1)#diff first and second +# self.write_vizdiff(fp,f,0,1,nbp) + + +# #displays +# def _display_cfields(self,fields,namelist): +# print("Cell Fields available are :\n") +# cfields = fields.GetCellData() +# for i in range((off:=len(namelist)),off+cfields.GetNumberOfArrays()): +# if cfields.GetArrayName(i - len(namelist)): +# print(str(i)+": "+cfields.GetArrayName(i-off)) +# namelist.append(cfields.GetArrayName(i)) +# else: +# print(f"fields {i} is undefined") + +# return namelist + +# def _display_pfields(self,fields,namelist): +# print("Point Fields available are :\n") +# pfields = fields.GetPointData() +# for i in range(len(namelist),len(namelist) + pfields.GetNumberOfArrays()): +# print(str(i)+": *"+pfields.GetArrayName(i)) +# namelist.append('*'+pfields.GetArrayName(i)) + +# return namelist + +# def display_fields(self, fields): +# namelist = [] +# self._display_pfields(fields,namelist) +# self._display_cfields(fields,namelist) + +# return namelist + +# def display_fields_mb(self, ugrid): +# it = ugrid.NewIterator() +# namelist = [] +# namelist.extend(self.display_fields(ugrid.GetDataSet(it))) +# return namelist + + +# #extract + +# def extract_data(self,filelist,olist): +# self.reader.SetFileName(filelist[0]) +# self.reader.Update() +# fields = self.reader.GetOutput() +# nv =0 +# nbp = 0 +# npp = 0 + +# #number of cells +# try: +# nv = fields.GetCellData().GetArray("ghostRank").GetNumberOfValues() +# npp = fields.GetPointData().GetArray("ghostRank").GetNumberOfValues() +# nbp = np.max( numpy_support.vtk_to_numpy( fields.GetPointData().GetArray("localToGlobalMap") ) ) +# except: +# it = self.fields.NewIterator() +# while not it.IsDoneWithTraversal(): +# nv += self.fields.GetDataSet(it).GetCellData().GetArray("ghostRank").GetNumberOfValues() +# npp += self.fields.GetDataSet(it).GetPointData().GetArray("ghostRank").GetNumberOfValues() +# nbp = np.max( numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) ) +# it.GoToNextItem() + + + +# ncnf = 0 +# npcnf = 0 +# ncc = 0 +# npcc = 0 +# for ifields in olist: +# try: +# if ifields[0] == '*': +# npcc = fields.GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() +# else: +# ncc = fields.GetCellData().GetArray(ifields).GetNumberOfComponents() +# except: +# it = fields.NewIterator() +# if ifields[0] == '*': +# npcc = fields.GetDataSet(it).GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() +# else: +# ncc = fields.GetDataSet(it).GetCellData().GetArray(ifields).GetNumberOfComponents() + + +# ncnf = ncnf + ncc +# npcnf = npcnf + npcc + +# f = np.zeros(shape=(nv,ncnf,len(filelist)),dtype='float') +# fp = np.zeros(shape=(npp,npcnf,len(filelist)),dtype='float') +# print(f"nv {nv} ncnf {ncnf} nb {nbp}") + +# i = 0 # for file loop +# for fileid in filelist: +# self.reader.SetFileName(fileid) +# print(fileid) +# self.reader.Update() +# # fields = self.reader.GetOutput() +# j = 0 # for field loop +# nc = 0 +# for nfields in olist: +# try: +# if nfields[0] == '*': +# field = self.fields.GetPointData().GetArray(nfields[1:]) +# else: +# field = self.fields.GetCellData().GetArray(nfields) + +# nc = field.GetNumberOfComponents() +# if nfields[0] == '*': +# f[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(nv,nc) +# else: +# fp[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(npb,nc) +# except: +# it = self.fields.NewIterator() +# start = 0 +# while not it.IsDoneWithTraversal(): +# # use localToGlobalMap +# if nfields[0] == '*': +# field = self.fields.GetDataSet(it).GetPointData().GetArray(nfields[1:]) +# nt = field.GetNumberOfValues() +# nc = field.GetNumberOfComponents() +# l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) +# fp[l2g,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) +# else: +# field = self.fields.GetDataSet(it).GetCellData().GetArray(nfields) +# nt = field.GetNumberOfValues() +# nc = field.GetNumberOfComponents() +# l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) +# f[l2g-nbp,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) +# it.GoToNextItem() + +# if nc>1 & i ==0 : +# self.flist[j:j+1] = [ self.flist[j]+"_"+str(k) for k in range(0,nc) ] +# j = j + nc +# i = i+1 + +# return fp,f,nbp + +# def write_report(self,fp,f,i1,i2): +# sp = fp.shape +# s = f.shape +# print(s) +# for i in range(0,s[1]): +# n = np.linalg.norm( f[:,i,i1]-f[:,i,i2], np.inf ) +# print(self.flist[i]+": "+str(n)+" ") +# for i in range(0,sp[1]): +# n = np.linalg.norm( fp[:,i,i1]-fp[:,i,i2], np.inf ) +# print(self.flist[i]+": "+str(n)+" ") + + +# def write_vizdiff(self,fp,f,i1,i2, nbp): +# # writer.SetDataModeToAscii() +# #mesh = vtk.vtkUnstructuredGrid() +# postfix = "" + +# print(self.flist) +# for i,fname in enumerate(self.flist): +# try: +# arr =numpy_support.numpy_to_vtk(np.abs( f[:,i,i1]-f[:,i,i2] )) +# arr.SetName("d"+fname) +# self.fields.GetCellData().AddArray(arr) +# except: +# it = self.fields.NewIterator() +# start = 0 +# while not it.IsDoneWithTraversal(): +# #scalar fill only +# if fname[0] == '*': +# l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) +# arr = numpy_support.numpy_to_vtk(np.abs( fp[l2g,i,i1]-fp[l2g,i,i2] )) +# arr.SetName("d"+fname) +# self.fields.GetDataSet(it).GetPointData().AddArray(arr) +# else: +# l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) +# arr = numpy_support.numpy_to_vtk(np.abs( f[l2g-nbp,i,i1]-f[l2g-nbp,i,i2] )) +# arr.SetName("d"+fname) +# self.fields.GetDataSet(it).GetCellData().AddArray(arr) +# it.GoToNextItem() + + +# self.writer.SetFileName( "diff_field_"+postfix+"."+self.extension ) +# self.writer.SetInputData(self.fields) +# self.writer.Write() \ No newline at end of file From 0ea6dcf885f8c4e76d6e66c968bc5039bf551a39 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Wed, 17 Dec 2025 18:54:25 +0100 Subject: [PATCH 03/11] Refactor the script to compute L1Diff --- .../generic_processing_tools/DiffFields.py | 489 +++++++++++------- 1 file changed, 295 insertions(+), 194 deletions(-) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py index b63273ac9..342facb04 100644 --- a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py +++ b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py @@ -1,19 +1,22 @@ # SPDX-License-Identifier: Apache 2.0 # SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies # SPDX-FileContributor: Jacques Franc +# ruff: noqa: E402 # disable Module level import not at top of file import logging import numpy as np import numpy.typing as npt -from functools import partial +from typing_extensions import Self, Any +# from functools import partial -from vtk.util import numpy_support -from vtkmodules.numpy_interface import dataset_adapter as dsa -from vtkmodules.vtkFiltersCore import vtkAppendFilter -from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataObject +from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataSet + +from geos.mesh.utils.arrayModifiers import createAttribute +from geos.mesh.utils.arrayHelpers import ( getAttributeSet, getNumberOfComponents, getArrayInObject ) +from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten from geos.utils.Logger import ( Logger, getLogger ) -__doc__ = """ +__doc__ = """ Module to diff(operate) on fields between two meshes composing based on globalIds indirection """ @@ -22,200 +25,299 @@ class DiffFieldsFilter: - def __init__(self): - """ Difference of Fields""" - self.appendFilter = vtkAppendFilter() - self.fieldnames = set() - self.logger = getLogger( loggerTitle ) - - def setMeshes(self, mesh1: vtkDataObject, mesh2: vtkDataObject): - """ Setting meshes to diff.""" - if isinstance(mesh1, vtkMultiBlockDataSet): - #discard partial for now - fieldnames_1 = self._displayFields_mb(mesh1) - else: - fieldnames_1 = self._displayFields(mesh1) - - if isinstance(mesh2, vtkMultiBlockDataSet): - #discard partial for now - fieldnames_2 = self._displayFields_mb(mesh2) + def __init__( + self: Self, + computeL2Diff: bool = False, + speHandler: bool = False, + ) -> None: + """ Difference of Fields. + + Args: + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. + """ + self.listMeshes: list[ vtkMultiBlockDataSet | vtkDataSet ] = [] + self.nbCells: int = 0 + self.nbPoints: int = 0 + self.idPointMax: int = 0 + + self.dicSharedAttributes: dict[ bool, set[ str ] ] = {} + self.dicAttributesToCompare: dict[ bool, set[ str ] ] = {} + self.dicAttributesDiffNames: dict[ bool, list[ str ] ] = {} + self.cellsAttributesArray: npt.NDArray[ np.float32 ] = np.array( [] ) + self.pointsAttributesArray: npt.NDArray[ np.float32 ] = np.array( [] ) + + self.computeL2Diff: bool = computeL2Diff + + self.outputMesh: vtkMultiBlockDataSet | vtkDataSet = vtkDataSet() + + # Logger. + self.logger: Logger + if not speHandler: + self.logger = getLogger( loggerTitle, True ) else: - fieldnames_2 = self._displayFields(mesh2) - - self.fieldnames = fieldnames_1.intersection(fieldnames_2) - self.logger.info(f"Fields have been dropped as unique to one of the meshes: {fieldnames_1.difference(fieldnames_2)}") - - -#displays - def _display_cfields(self,fields,namelist): - print("Cell Fields available are :\n") - cfields = fields.GetCellData() - for i in range((off:=len(namelist)),off+cfields.GetNumberOfArrays()): - if cfields.GetArrayName(i - len(namelist)): - print(str(i)+": "+cfields.GetArrayName(i-off)) - namelist.append(cfields.GetArrayName(i)) - else: - print(f"fields {i} is undefined") - - return namelist - - def _display_pfields(self,fields,namelist): - print("Point Fields available are :\n") - pfields = fields.GetPointData() - for i in range(len(namelist),len(namelist) + pfields.GetNumberOfArrays()): - print(str(i)+": *"+pfields.GetArrayName(i)) - namelist.append('*'+pfields.GetArrayName(i)) - - return namelist - - def display_fields(self, fields): - namelist = [] - self._display_pfields(fields,namelist) - self._display_cfields(fields,namelist) - - return namelist - - def display_fields_mb(self, ugrid): - it = ugrid.NewIterator() - namelist = [] - namelist.extend(self.display_fields(ugrid.GetDataSet(it))) - return namelist - -# - - def _computeL(self, f, callback = partial(np.linalg.norm(ord=np.inf))): - """ compute by default inf norm """ - s = f.shape - #loop - sp = fp.shape - s = f.shape - for i in range(0,s[1]): - n = callback( f[:,i,i1]-f[:,i,i2]) - print(self.flist[i]+": "+str(n)+" ") - for i in range(0,sp[1]): - n = callback( fp[:,i,i1]-fp[:,i,i2]) - print(self.flist[i]+": "+str(n)+" ") - self.logger.info("Lmax norm for fields {name} is {value}") - - def _computeL1(self): - # writer.SetDataModeToAscii() - #mesh = vtk.vtkUnstructuredGrid() - postfix = "" - - print(self.flist) - for i,fname in enumerate(self.flist): - try: - arr =numpy_support.numpy_to_vtk(np.abs( f[:,i,i1]-f[:,i,i2] )) - arr.SetName("d"+fname) - self.fields.GetCellData().AddArray(arr) - except: - it = self.fields.NewIterator() - start = 0 - while not it.IsDoneWithTraversal(): - #scalar fill only - if fname[0] == '*': - l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) - arr = numpy_support.numpy_to_vtk(np.abs( fp[l2g,i,i1]-fp[l2g,i,i2] )) - arr.SetName("d"+fname) - self.fields.GetDataSet(it).GetPointData().AddArray(arr) - else: - l2g = numpy_support.numpy_to_vtk( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) - arr = numpy_support.numpy_to_vtk(np.abs( f[l2g-nbp,i,i1]-f[l2g-nbp,i,i2] )) - arr.SetName("d"+fname) - self.fields.GetDataSet(it).GetCellData().AddArray(arr) - it.GoToNextItem() + self.logger = logging.getLogger( loggerTitle ) + self.logger.setLevel( logging.INFO ) + def setLoggerHandler( self: Self, handler: logging.Handler ) -> None: + """Set a specific handler for the filter logger. - def _getdata(self): - # size the output numpy array - nv =0 - nbp = 0 - npp = 0 - - try: - nv = fields.GetCellData().GetArray("ghostRank").GetNumberOfValues() - npp = fields.GetPointData().GetArray("ghostRank").GetNumberOfValues() - nbp = np.max( numpy_support.vtk_to_numpy( fields.GetPointData().GetArray("localToGlobalMap") ) ) - except: - it = self.fields.NewIterator() - while not it.IsDoneWithTraversal(): - nv += self.fields.GetDataSet(it).GetCellData().GetArray("ghostRank").GetNumberOfValues() - npp += self.fields.GetDataSet(it).GetPointData().GetArray("ghostRank").GetNumberOfValues() - nbp = np.max( numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) ) - it.GoToNextItem() - - cnf = 0 - ncc = 0 - npcc = 0 - for ifields in self.olist: - try: - if ifields[0] == '*': - npcc = fields.GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() - else: - ncc = fields.GetCellData().GetArray(ifields).GetNumberOfComponents() - except: - it = fields.NewIterator() - if ifields[0] == '*': - npcc = fields.GetDataSet(it).GetPointData().GetArray(ifields[1:]).GetNumberOfComponents() + 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: + # This warning does not count for the number of warning created during the application of the filter. + self.logger.warning( "The logger already has an handler, to use yours set the argument 'speHandler' to True" + " during the filter initialization." ) + + def setMeshes( + self: Self, + listMeshes: list[ vtkMultiBlockDataSet | vtkDataSet ], + ) -> None: + """Setter of the two meshes with the attributes to compare. + + Setting the two meshes will automatically compute the dictionary with the shared attribute per localization. + + Args: + listMeshes (list[vtkMultiBlockDataSet | vtkDataSet]): The list of the meshes to compare. + + Raises: + TypeError: The meshes do not have the same type. + ValueError: The meshes do not have the same cells or points number or datasets indexes or the number of meshes is to small. + """ + if len( listMeshes ) < 2: + raise ValueError ( "The list of meshes must contain at least two meshes." ) + + meshTypeRef: str = listMeshes[ 0 ].GetClassName() + for mesh in listMeshes[ 1: ]: + if mesh.GetClassName() != meshTypeRef: + raise TypeError( "All the meshes must have the same type." ) + + nbDataSet: int + idPointMax: int + nbCellsRef: int + nbPointsRef: int + if isinstance( listMeshes[ 0 ], vtkDataSet ): + nbDataSet = 1 + idPointMax = np.max( getArrayInObject( listMeshes[ 0 ], "localToGlobalMap", True ) ) + nbCellsRef = listMeshes[ 0 ].GetNumberOfCells() + nbPointsRef = listMeshes[ 0 ].GetNumberOfPoints() + for mesh in listMeshes[ 1: ]: + if mesh.GetNumberOfCells() != nbCellsRef: + raise ValueError( "All the meshes must have the same number of cells." ) + if mesh.GetNumberOfPoints() != nbPointsRef: + raise ValueError( "All the meshes must have the same number of points." ) + elif isinstance( listMeshes[ 0 ], vtkMultiBlockDataSet ): + blockElementIndexesFlattenRef: list[ int ] = getBlockElementIndexesFlatten( listMeshes[ 0 ] ) + nbDataSet = len( blockElementIndexesFlattenRef ) + for mesh in listMeshes[ 1: ]: + if getBlockElementIndexesFlatten( mesh ) != blockElementIndexesFlattenRef: + raise ValueError( "All the meshes do not have the same blocks indexes.") + + for blockId in blockElementIndexesFlattenRef: + datasetTypeRef: str = listMeshes[ 0 ].GetDataSet( blockId ).GetClassName() + for mesh in listMeshes[ 1: ]: + if mesh.GetDataSet( blockId ).GetClassName() != datasetTypeRef: + raise TypeError( "The datasets with the same flatten index of all the meshes must have the same type.") + + datasetRef: vtkDataSet = vtkDataSet.SafeDownCast( listMeshes[ 0 ].GetDataSet( blockId ) ) + idPointMax = np.max( getArrayInObject( datasetRef, "localToGlobalMap", True ) ) + nbCellsRef = datasetRef.GetNumberOfCells() + nbPointsRef = datasetRef.GetNumberOfPoints() + for mesh in listMeshes[ 1: ]: + datasetTest: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( blockId ) ) + if datasetTest.GetNumberOfCells() != nbCellsRef: + raise ValueError( "The datasets with the same flatten index of all the meshes must have the same number of cells." ) + if datasetTest.GetNumberOfPoints() != nbPointsRef: + raise ValueError( "The datasets with the same flatten index of all the meshes must have the same number of points." ) + else: + raise TypeError( "The meshes must be inherited from vtkMultiBlockDataSet or vtkDataSet.") + + self.listMeshes = listMeshes + self.idPointMax = idPointMax + self.nbCells = nbCellsRef * nbDataSet + self.nbPoints = nbPointsRef * nbDataSet + self.computeSharedAttributes() + self.outputMesh = listMeshes[ 0 ].NewInstance() + self.outputMesh.ShallowCopy( listMeshes[ 0 ] ) + + return + + def computeSharedAttributes( self: Self ) -> None: + """Compute the dictionary with the shared attributes per localization between the two meshes. + + Keys of the dictionary are the attribute localization and the value are the shared attribute per localization. + """ + for piece in [ True, False ]: + sharedAttributes: set[ str ] = getAttributeSet( self.listMeshes[ 0 ], piece ) + for mesh in self.listMeshes[ 1: ]: + sharedAttributes.update( getAttributeSet( mesh, piece ) ) + if sharedAttributes != set(): + self.dicSharedAttributes[ piece ] = sharedAttributes + + return + + def getSharedAttribute( self: Self ) -> dict[ bool, set[ str ] ]: + """Getter of the dictionary with the shared attributes per localization.""" + return self.dicSharedAttributes + + def logSharedAttributeInfo( self: Self ) -> None: + """Log the shared attributes per localization.""" + if self.dicSharedAttributes == {}: + self.logger.warning( "The two meshes do not share any attribute." ) + else: + for piece, sharedAttributes in self.dicSharedAttributes.items(): + self.logger.info( f"Shared attributes on { piece } are { sharedAttributes }." ) + + return + + def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ bool, set[ str ] ] ) -> None: + """Setter of the dictionary with the shared attribute per localization to compare. + + Args: + dicAttributesToCompare (dict[bool, set[ str ]]): The dictionary with the attributes to compare per localization. + + Raises: + ValueError: At least one attribute to compare is not a shared attribute. + """ + for piece, sharedAttributesToCompare in dicAttributesToCompare.items(): + if not sharedAttributesToCompare.issubset( self.dicSharedAttributes[ piece ] ): + wrongAttributes: set[ str ] = sharedAttributesToCompare.difference( self.dicSharedAttributes[ piece ] ) + raise ValueError( f"The attributes to compare { wrongAttributes } are not shared attributes.") + + for piece in dicAttributesToCompare: + self.dicAttributesDiffNames[ piece ] = [] + + nbCellsComponents = 0 + nbPointsComponents = 0 + for piece, sharedAttributesToCompare in dicAttributesToCompare.items(): + for attributeName in sharedAttributesToCompare: + self.dicAttributesDiffNames[ piece ].append( f"diff_{ attributeName }" ) + if piece: + nbPointsComponents += getNumberOfComponents( self.listMeshes[ 0 ], attributeName, piece ) else: - ncc = fields.GetDataSet(it).GetCellData().GetArray(ifields).GetNumberOfComponents() - - - ncnf = ncnf + ncc - npcnf = npcnf + npcc - - f = np.zeros(shape=(nv,ncnf,len(filelist)),dtype='float') - fp = np.zeros(shape=(npp,npcnf,len(filelist)),dtype='float') - print(f"nv {nv} ncnf {ncnf} nb {nbp}") - - i = 0 # for file loop - for fileid in filelist: - self.reader.SetFileName(fileid) - print(fileid) - self.reader.Update() - # fields = self.reader.GetOutput() - j = 0 # for field loop - nc = 0 - for nfields in self.olist: - try: - if nfields[0] == '*': - field = self.fields.GetPointData().GetArray(nfields[1:]) + nbCellsComponents += getNumberOfComponents( self.listMeshes[ 0 ], attributeName, piece ) + + self.dicAttributesToCompare = dicAttributesToCompare + self.cellsAttributesArray = np.zeros( shape=( self.nbCells, nbCellsComponents, len( self.listMeshes ) ), dtype=np.float32 ) + self.pointsAttributesArray = np.zeros( shape=( self.nbPoints, nbPointsComponents, len( self.listMeshes ) ), dtype=np.float32 ) + + return + + def getDicAttributesToCompare( self: Self ) -> dict[ bool, set[ str ] ]: + """Getter of the dictionary of the attribute to compare per localization.""" + return self.dicAttributesToCompare + + def applyFilter( self: Self ) -> None: + """Apply the diffFieldsFilter.""" + self.logger.info( f"Apply filter { self.logger.name }." ) + + if self.listMeshes == []: + raise ValueError( "Set a list of meshes to compare." ) + + if self.dicAttributesToCompare == {}: + raise ValueError( "Set the attribute to compare per localization." ) + + self.computePointsAndCellsAttributesArrays() + self.computeL1() + # if self.computeL2Diff: + # self.computeL2() + + self.logger.info( f"The filter { self.logger.name } succeed." ) + + return + + def computePointsAndCellsAttributesArrays( self: Self ) -> None: + """Compute one array per localization with all the values of all the attributes to compare.""" + for piece, sharedAttributesToCompare in self.dicAttributesToCompare.items(): + idComponents: int = 0 + for attributeName in sharedAttributesToCompare: + arrayAttributeData: npt.NDArray[ Any ] + nbAttributeComponents: int + for idMesh, mesh in enumerate( self.listMeshes ): + if isinstance( mesh, vtkDataSet ): + arrayAttributeData = getArrayInObject( mesh, attributeName, piece ) + nbAttributeComponents = getNumberOfComponents( mesh, attributeName, piece ) + if piece: + self.pointsAttributesArray[ :, idComponents : idComponents + nbAttributeComponents, idMesh ] = arrayAttributeData.reshape( self.nbPoints, nbAttributeComponents ) + else: + self.cellsAttributesArray[ :, idComponents : idComponents + nbAttributeComponents, idMesh ] = arrayAttributeData.reshape( self.nbCells, nbAttributeComponents ) else: - field = self.fields.GetCellData().GetArray(nfields) - - nc = field.GetNumberOfComponents() - if nfields[0] == '*': - f[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(nv,nc) + it = mesh.NewIterator() + while not it.IsDoneWithTraversal(): + # use localToGlobalMap + dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( it ) ) + arrayAttributeData = getArrayInObject( dataset, attributeName, piece ) + nbAttributeComponents = getNumberOfComponents( dataset, attributeName, piece ) + lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) + if piece: + nbPoints: int = dataset.GetNumberOfPoints() + self.pointsAttributesArray[ lToG, idComponents : idComponents + nbAttributeComponents, idMesh ] += arrayAttributeData.reshape( int( nbPoints / nbAttributeComponents ), nbAttributeComponents ) + else: + nbCells: int = dataset.GetNumberOfCells() + self.cellsAttributesArray[ lToG - self.idPointMax, idComponents : idComponents + nbAttributeComponents, idMesh ] += arrayAttributeData.reshape( int( nbCells / nbAttributeComponents ), nbAttributeComponents ) + it.GoToNextItem() + + if nbAttributeComponents > 1: + self.dicAttributesDiffNames[ piece ][ idComponents : idComponents + 1 ] = str( [ self.dicAttributesDiffNames[ piece ][ idComponents ] + "_component" + str( idAttributeComponent ) for idAttributeComponent in range( nbAttributeComponents ) ] ) + + idComponents += nbAttributeComponents + + return + + def computeL1( self: Self ) -> None: + """Compute the L1 diff for all the wanted attribute and create attribute with it on the output mesh.""" + for attributeId, dicItems in enumerate( self.dicAttributesDiffNames.items() ): + piece: bool + attributeDiffName: list[ str ] + piece, attributeDiffName = dicItems + attributeArray: npt.NDArray[ Any ] + if isinstance( self.listMeshes[ 0 ], vtkDataSet ): + if piece: + attributeArray = np.abs( self.pointsAttributesArray[ :, attributeId, 0 ] - self.pointsAttributesArray[ :, attributeId, 1 ] ) + else: + attributeArray = np.abs( self.cellsAttributesArray[ :, attributeId, 0 ] - self.cellsAttributesArray[ :, attributeId, 1 ] ) + createAttribute( self.outputMesh, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) + else: + it = self.outputMesh.NewIterator() + while not it.IsDoneWithTraversal(): + dataset: vtkDataSet = vtkDataSet.SafeDownCast( self.outputMesh.GetDataSet( it ) ) + lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) + if piece: + attributeArray = np.abs( self.cellsAttributesArray[ lToG, attributeId, 0 ] - self.cellsAttributesArray[ lToG, attributeId, 1 ] ) else: - fp[:,j:j+nc,i] = numpy_support.vtk_to_numpy(field).reshape(npb,nc) - except: - it = self.fields.NewIterator() - start = 0 - while not it.IsDoneWithTraversal(): - # use localToGlobalMap - if nfields[0] == '*': - field = self.fields.GetDataSet(it).GetPointData().GetArray(nfields[1:]) - nt = field.GetNumberOfValues() - nc = field.GetNumberOfComponents() - l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetPointData().GetArray("localToGlobalMap") ) - fp[l2g,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) - else: - field = self.fields.GetDataSet(it).GetCellData().GetArray(nfields) - nt = field.GetNumberOfValues() - nc = field.GetNumberOfComponents() - l2g = numpy_support.vtk_to_numpy( self.fields.GetDataSet(it).GetCellData().GetArray("localToGlobalMap") ) - f[l2g-nbp,j:j+nc,i] += numpy_support.vtk_to_numpy(field).reshape(int(nt/nc),nc) - it.GoToNextItem() - - if nc>1 & i ==0 : - self.flist[j:j+1] = [ self.flist[j]+"_"+str(k) for k in range(0,nc) ] - j = j + nc - i = i+1 - - return fp,f,nbp + attributeArray = np.abs( self.cellsAttributesArray[ lToG - self.idPointMax, attributeId, 0 ] - self.cellsAttributesArray[ lToG - self.idPointMax, attributeId, 1 ] ) + createAttribute( dataset, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) + it.GoToNextItem() - def getOutput(self) -> vtkDataObject: - return self.GetOutput() + return + + # def computeL2(self, f, callback = partial(np.linalg.norm(ord=np.inf))): + # """ compute by default inf norm """ + # s = f.shape + # #loop + # sp = fp.shape + # s = f.shape + # for i in range(0,s[1]): + # n = callback( f[:,i,i1]-f[:,i,i2]) + # print(self.flist[i]+": "+str(n)+" ") + # for i in range(0,sp[1]): + # n = callback( fp[:,i,i1]-fp[:,i,i2]) + # print(self.flist[i]+": "+str(n)+" ") + # self.logger.info("Lmax norm for fields {name} is {value}") + + def getOutput( self: Self ) -> vtkMultiBlockDataSet | vtkDataSet: + """Return the mesh with the computed diff as attributes for the wanted attributes. + + Returns: + (vtkMultiBlockDataSet | vtkDataSet): The mesh with the computed attributes diff. + """ + return self.outputMesh # ##### # class diff_visu: @@ -449,4 +551,3 @@ def getOutput(self) -> vtkDataObject: # self.writer.SetInputData(self.fields) # self.writer.Write() - \ No newline at end of file From 1b3dc366db0973bc6a7d8b0b58a7b88633cc7632 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Mon, 22 Dec 2025 14:18:28 +0100 Subject: [PATCH 04/11] test diffFields --- geos-processing/tests/test_DiffMesh.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 geos-processing/tests/test_DiffMesh.py diff --git a/geos-processing/tests/test_DiffMesh.py b/geos-processing/tests/test_DiffMesh.py new file mode 100644 index 000000000..640c79a07 --- /dev/null +++ b/geos-processing/tests/test_DiffMesh.py @@ -0,0 +1,12 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Romain Baville +# 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 + +from typing import Any +from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet + +from geos.processing.generic_processing_tools.DiffFields import DiffFields From 75198ebb67ddc2a955fd2bccc45928390a227b53 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 18:50:16 +0100 Subject: [PATCH 05/11] Clean variables names and add the doc --- .../generic_processing_tools/DiffFields.py | 265 +++++++++--------- 1 file changed, 138 insertions(+), 127 deletions(-) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py index 342facb04..ddd5f8af3 100644 --- a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py +++ b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py @@ -17,39 +17,80 @@ from geos.utils.Logger import ( Logger, getLogger ) __doc__ = """ -Module to diff(operate) on fields between two meshes composing based on globalIds indirection +Attributes Diff is a vtk that compute L1 and L2 differences between attributes shared by two identical meshes. +Input meshes cans be vtkDataSet or vtkMultiBlockDataSet. + +To use the filter: + +.. code-block:: python + + import logging + from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataSet + + from geos.processing.generic_processing_tools.AttributesDiff import AttributesDiff + + # Filter inputs: + computeL2Diff: bool # defaults to False + speHandler: bool # defaults to False + + # Instantiate the filter: + attributesDiffFilter: AttributesDiff = AttributesDiff( computeL2Diff, speHandler ) + + # Set the handler of yours (only if speHandler is True): + yourHandler: logging.Handler + attributesDiffFilter.setLoggerHandler( yourHandler ) + + # Set the meshes to compare: + listMeshes: list[ vtkMultiBlockDataSet | vtkDataSet ] + attributesDiffFilter.setMeshes( listMeshes ) + + # Log the shared attributes info (optional): + attributesDiffFilter.logSharedAttributeInfo() + + # Get the shared attributes (optional): + dicSharedAttributes: dict[ bool, set[ str ] ] + dicSharedAttributes = attributesDiffFilter.getDicSharedAttribute() + + # Set the attributes to compare: + dicAttributesToCompare: dict[ bool, set[ str ] ] + attributesDiffFilter.setDicAttributesToCompare( dicAttributesToCompare ) + + # Do calculations: + attributesDiffFilter.applyFilter() + + # Get the mesh with the computed attributes differences as new attributes: + outputMesh: vtkMultiBlockDataSet | vtkDataSet + outputMesh = attributesDiffFilter.getOutput() """ -loggerTitle : str = "DiffFieldsFilter" +loggerTitle : str = "Attributes Diff" -class DiffFieldsFilter: +class AttributesDiff: def __init__( self: Self, computeL2Diff: bool = False, speHandler: bool = False, ) -> None: - """ Difference of Fields. + """Compute differences (L1 and L2) between two identical meshes attribute. Args: + computeL2Diff (bool, optional): True to compute the L2 difference speHandler (bool, optional): True to use a specific handler, False to use the internal handler. Defaults to False. """ self.listMeshes: list[ vtkMultiBlockDataSet | vtkDataSet ] = [] - self.nbCells: int = 0 - self.nbPoints: int = 0 - self.idPointMax: int = 0 + self.dicNbElements: dict[ bool, int ] = {} self.dicSharedAttributes: dict[ bool, set[ str ] ] = {} self.dicAttributesToCompare: dict[ bool, set[ str ] ] = {} self.dicAttributesDiffNames: dict[ bool, list[ str ] ] = {} - self.cellsAttributesArray: npt.NDArray[ np.float32 ] = np.array( [] ) - self.pointsAttributesArray: npt.NDArray[ np.float32 ] = np.array( [] ) + self.dicAttributesArray: dict[ bool, npt.NDArray[ np.float32 ] ] = {} self.computeL2Diff: bool = computeL2Diff - self.outputMesh: vtkMultiBlockDataSet | vtkDataSet = vtkDataSet() + self.outputMesh: vtkMultiBlockDataSet | vtkDataSet = vtkMultiBlockDataSet() # Logger. self.logger: Logger @@ -90,79 +131,64 @@ def setMeshes( TypeError: The meshes do not have the same type. ValueError: The meshes do not have the same cells or points number or datasets indexes or the number of meshes is to small. """ - if len( listMeshes ) < 2: - raise ValueError ( "The list of meshes must contain at least two meshes." ) - - meshTypeRef: str = listMeshes[ 0 ].GetClassName() - for mesh in listMeshes[ 1: ]: - if mesh.GetClassName() != meshTypeRef: - raise TypeError( "All the meshes must have the same type." ) - - nbDataSet: int - idPointMax: int - nbCellsRef: int - nbPointsRef: int - if isinstance( listMeshes[ 0 ], vtkDataSet ): - nbDataSet = 1 - idPointMax = np.max( getArrayInObject( listMeshes[ 0 ], "localToGlobalMap", True ) ) - nbCellsRef = listMeshes[ 0 ].GetNumberOfCells() - nbPointsRef = listMeshes[ 0 ].GetNumberOfPoints() - for mesh in listMeshes[ 1: ]: - if mesh.GetNumberOfCells() != nbCellsRef: - raise ValueError( "All the meshes must have the same number of cells." ) - if mesh.GetNumberOfPoints() != nbPointsRef: - raise ValueError( "All the meshes must have the same number of points." ) - elif isinstance( listMeshes[ 0 ], vtkMultiBlockDataSet ): - blockElementIndexesFlattenRef: list[ int ] = getBlockElementIndexesFlatten( listMeshes[ 0 ] ) - nbDataSet = len( blockElementIndexesFlattenRef ) - for mesh in listMeshes[ 1: ]: - if getBlockElementIndexesFlatten( mesh ) != blockElementIndexesFlattenRef: - raise ValueError( "All the meshes do not have the same blocks indexes.") - - for blockId in blockElementIndexesFlattenRef: - datasetTypeRef: str = listMeshes[ 0 ].GetDataSet( blockId ).GetClassName() - for mesh in listMeshes[ 1: ]: - if mesh.GetDataSet( blockId ).GetClassName() != datasetTypeRef: - raise TypeError( "The datasets with the same flatten index of all the meshes must have the same type.") - - datasetRef: vtkDataSet = vtkDataSet.SafeDownCast( listMeshes[ 0 ].GetDataSet( blockId ) ) - idPointMax = np.max( getArrayInObject( datasetRef, "localToGlobalMap", True ) ) - nbCellsRef = datasetRef.GetNumberOfCells() - nbPointsRef = datasetRef.GetNumberOfPoints() - for mesh in listMeshes[ 1: ]: - datasetTest: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( blockId ) ) - if datasetTest.GetNumberOfCells() != nbCellsRef: - raise ValueError( "The datasets with the same flatten index of all the meshes must have the same number of cells." ) - if datasetTest.GetNumberOfPoints() != nbPointsRef: - raise ValueError( "The datasets with the same flatten index of all the meshes must have the same number of points." ) + if len( listMeshes ) != 2: + raise ValueError ( "The list of meshes must contain two meshes." ) + + if listMeshes[ 0 ].GetClassName() != listMeshes[ 1 ].GetClassName(): + raise TypeError( "The meshes must have the same type." ) + + dicMeshesMaxElementId: dict[ bool, list[ int ] ] = { False: [ 0, 0 ], True: [ 0, 0 ], } + if isinstance( self.listMeshes[ 0 ], vtkDataSet ): + for meshId, mesh in enumerate( listMeshes ): + for piece in dicMeshesMaxElementId: + dicMeshesMaxElementId[ piece ][ meshId ] = np.max( getArrayInObject( mesh, "localToGlobalMap", piece ) ) + elif isinstance( self.listMeshes[ 0 ], vtkMultiBlockDataSet ): + setDatasetType: set[ str ] = set() + for meshId, mesh in enumerate( listMeshes ): + listMeshBlockId: list[ int ] = getBlockElementIndexesFlatten( mesh ) + for meshBlockId in listMeshBlockId: + setDatasetType.add( mesh.GetDataSet( meshBlockId ).GetClassName() ) + dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( meshBlockId ) ) + for piece in dicMeshesMaxElementId: + dicMeshesMaxElementId[ piece ][ meshId ] = max( dicMeshesMaxElementId[ piece ][ meshId ], np.max( getArrayInObject( dataset, "localToGlobalMap", piece ) ) ) + if piece: + a= dataset.GetNumberOfPoints() + else: + a= dataset.GetNumberOfCells() + if a - len( getArrayInObject( dataset, "localToGlobalMap", piece ) ) != 0: + print("r") + if len( setDatasetType ) != 1: + raise TypeError( "All datasets of the meshes must have the same type." ) else: raise TypeError( "The meshes must be inherited from vtkMultiBlockDataSet or vtkDataSet.") + for piece, listMeshMaxElementId in dicMeshesMaxElementId.items(): + if listMeshMaxElementId[ 0 ] != listMeshMaxElementId[ 1 ]: + raise ValueError( f"The total number of { piece } in the meshes must be the same." ) + self.listMeshes = listMeshes - self.idPointMax = idPointMax - self.nbCells = nbCellsRef * nbDataSet - self.nbPoints = nbPointsRef * nbDataSet - self.computeSharedAttributes() + self.dicNbElements[ True ] = dicMeshesMaxElementId[ True ][ 0 ] + 1 + self.dicNbElements[ False ] = dicMeshesMaxElementId[ False ][ 0 ] + 1 self.outputMesh = listMeshes[ 0 ].NewInstance() self.outputMesh.ShallowCopy( listMeshes[ 0 ] ) + self._computeDicSharedAttributes() return - def computeSharedAttributes( self: Self ) -> None: + def _computeDicSharedAttributes( self: Self ) -> None: """Compute the dictionary with the shared attributes per localization between the two meshes. Keys of the dictionary are the attribute localization and the value are the shared attribute per localization. """ for piece in [ True, False ]: - sharedAttributes: set[ str ] = getAttributeSet( self.listMeshes[ 0 ], piece ) - for mesh in self.listMeshes[ 1: ]: - sharedAttributes.update( getAttributeSet( mesh, piece ) ) - if sharedAttributes != set(): - self.dicSharedAttributes[ piece ] = sharedAttributes + setSharedAttributes: set[ str ] = getAttributeSet( self.listMeshes[ 0 ], piece ) + setSharedAttributes.update( getAttributeSet( self.listMeshes[ 1 ], piece ) ) + if setSharedAttributes != set(): + self.dicSharedAttributes[ piece ] = setSharedAttributes return - def getSharedAttribute( self: Self ) -> dict[ bool, set[ str ] ]: + def getDicSharedAttribute( self: Self ) -> dict[ bool, set[ str ] ]: """Getter of the dictionary with the shared attributes per localization.""" return self.dicSharedAttributes @@ -185,27 +211,30 @@ def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ bool, s Raises: ValueError: At least one attribute to compare is not a shared attribute. """ - for piece, sharedAttributesToCompare in dicAttributesToCompare.items(): - if not sharedAttributesToCompare.issubset( self.dicSharedAttributes[ piece ] ): - wrongAttributes: set[ str ] = sharedAttributesToCompare.difference( self.dicSharedAttributes[ piece ] ) + for piece, setSharedAttributesToCompare in dicAttributesToCompare.items(): + if not setSharedAttributesToCompare.issubset( self.dicSharedAttributes[ piece ] ): + wrongAttributes: set[ str ] = setSharedAttributesToCompare.difference( self.dicSharedAttributes[ piece ] ) raise ValueError( f"The attributes to compare { wrongAttributes } are not shared attributes.") - for piece in dicAttributesToCompare: - self.dicAttributesDiffNames[ piece ] = [] - - nbCellsComponents = 0 - nbPointsComponents = 0 - for piece, sharedAttributesToCompare in dicAttributesToCompare.items(): - for attributeName in sharedAttributesToCompare: - self.dicAttributesDiffNames[ piece ].append( f"diff_{ attributeName }" ) - if piece: - nbPointsComponents += getNumberOfComponents( self.listMeshes[ 0 ], attributeName, piece ) + dicNbComponents: dict[ bool, int ] = {} + dicAttributesDiffNames: dict[ bool, list[ str ] ] = {} + dicAttributesArray: dict[ bool, npt.NDArray[ np.float32 ] ] = {} + for piece, setSharedAttributesToCompare in dicAttributesToCompare.items(): + dicNbComponents[ piece ] = 0 + dicAttributesDiffNames[ piece ] = [] + for attributeName in setSharedAttributesToCompare: + nbAttributeComponents = getNumberOfComponents( self.outputMesh, attributeName, piece ) + dicNbComponents[ piece ] += nbAttributeComponents + diffAttributeName: str = f"diff_{ attributeName }" + if nbAttributeComponents > 1: + dicAttributesDiffNames[ piece ].extend( [ diffAttributeName + "_component" + str( idAttributeComponent ) for idAttributeComponent in range( nbAttributeComponents ) ] ) else: - nbCellsComponents += getNumberOfComponents( self.listMeshes[ 0 ], attributeName, piece ) + dicAttributesDiffNames[ piece ].append( diffAttributeName ) + dicAttributesArray[ piece ] = np.zeros( shape=( self.dicNbElements[ piece ], dicNbComponents[ piece ], 2 ), dtype=np.float32 ) + self.dicAttributesArray = dicAttributesArray self.dicAttributesToCompare = dicAttributesToCompare - self.cellsAttributesArray = np.zeros( shape=( self.nbCells, nbCellsComponents, len( self.listMeshes ) ), dtype=np.float32 ) - self.pointsAttributesArray = np.zeros( shape=( self.nbPoints, nbPointsComponents, len( self.listMeshes ) ), dtype=np.float32 ) + self.dicAttributesDiffNames = dicAttributesDiffNames return @@ -213,6 +242,10 @@ def getDicAttributesToCompare( self: Self ) -> dict[ bool, set[ str ] ]: """Getter of the dictionary of the attribute to compare per localization.""" return self.dicAttributesToCompare + def getDicAttributesDiffNames( self: Self ) -> dict[ bool, list[ str ] ]: + """Getter of the dictionary with the name of the attribute created with the calculated attributes diff.""" + return self.dicAttributesDiffNames + def applyFilter( self: Self ) -> None: """Apply the diffFieldsFilter.""" self.logger.info( f"Apply filter { self.logger.name }." ) @@ -223,8 +256,8 @@ def applyFilter( self: Self ) -> None: if self.dicAttributesToCompare == {}: raise ValueError( "Set the attribute to compare per localization." ) - self.computePointsAndCellsAttributesArrays() - self.computeL1() + self._computeDicAttributesArray() + self._computeL1() # if self.computeL2Diff: # self.computeL2() @@ -232,68 +265,46 @@ def applyFilter( self: Self ) -> None: return - def computePointsAndCellsAttributesArrays( self: Self ) -> None: - """Compute one array per localization with all the values of all the attributes to compare.""" + def _computeDicAttributesArray( self: Self ) -> None: + """Compute the dictionary with one array per localization with all the values of all the attributes to compare.""" for piece, sharedAttributesToCompare in self.dicAttributesToCompare.items(): idComponents: int = 0 for attributeName in sharedAttributesToCompare: arrayAttributeData: npt.NDArray[ Any ] nbAttributeComponents: int - for idMesh, mesh in enumerate( self.listMeshes ): + for meshId, mesh in enumerate( self.listMeshes ): if isinstance( mesh, vtkDataSet ): arrayAttributeData = getArrayInObject( mesh, attributeName, piece ) nbAttributeComponents = getNumberOfComponents( mesh, attributeName, piece ) - if piece: - self.pointsAttributesArray[ :, idComponents : idComponents + nbAttributeComponents, idMesh ] = arrayAttributeData.reshape( self.nbPoints, nbAttributeComponents ) - else: - self.cellsAttributesArray[ :, idComponents : idComponents + nbAttributeComponents, idMesh ] = arrayAttributeData.reshape( self.nbCells, nbAttributeComponents ) + self.dicAttributesArray[ piece ][ :, idComponents : idComponents + nbAttributeComponents, meshId ] = arrayAttributeData.reshape( self.dicNbElements[ piece ], nbAttributeComponents ) else: - it = mesh.NewIterator() - while not it.IsDoneWithTraversal(): - # use localToGlobalMap - dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( it ) ) + listMeshBlockId: list[ int ] = getBlockElementIndexesFlatten( mesh ) + for meshBlockId in listMeshBlockId: + dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( meshBlockId ) ) arrayAttributeData = getArrayInObject( dataset, attributeName, piece ) nbAttributeComponents = getNumberOfComponents( dataset, attributeName, piece ) lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) - if piece: - nbPoints: int = dataset.GetNumberOfPoints() - self.pointsAttributesArray[ lToG, idComponents : idComponents + nbAttributeComponents, idMesh ] += arrayAttributeData.reshape( int( nbPoints / nbAttributeComponents ), nbAttributeComponents ) - else: - nbCells: int = dataset.GetNumberOfCells() - self.cellsAttributesArray[ lToG - self.idPointMax, idComponents : idComponents + nbAttributeComponents, idMesh ] += arrayAttributeData.reshape( int( nbCells / nbAttributeComponents ), nbAttributeComponents ) - it.GoToNextItem() - - if nbAttributeComponents > 1: - self.dicAttributesDiffNames[ piece ][ idComponents : idComponents + 1 ] = str( [ self.dicAttributesDiffNames[ piece ][ idComponents ] + "_component" + str( idAttributeComponent ) for idAttributeComponent in range( nbAttributeComponents ) ] ) + self.dicAttributesArray[ piece ][ lToG, idComponents : idComponents + nbAttributeComponents, meshId ] = arrayAttributeData.reshape( len( lToG ), nbAttributeComponents ) idComponents += nbAttributeComponents return - def computeL1( self: Self ) -> None: + def _computeL1( self: Self ) -> None: """Compute the L1 diff for all the wanted attribute and create attribute with it on the output mesh.""" - for attributeId, dicItems in enumerate( self.dicAttributesDiffNames.items() ): - piece: bool - attributeDiffName: list[ str ] - piece, attributeDiffName = dicItems - attributeArray: npt.NDArray[ Any ] - if isinstance( self.listMeshes[ 0 ], vtkDataSet ): - if piece: - attributeArray = np.abs( self.pointsAttributesArray[ :, attributeId, 0 ] - self.pointsAttributesArray[ :, attributeId, 1 ] ) + for piece in self.dicAttributesDiffNames: + for attributeId, attributeDiffName in enumerate( self.dicAttributesDiffNames[ piece ] ): + attributeArray: npt.NDArray[ Any ] + if isinstance( self.listMeshes[ 0 ], vtkDataSet ): + attributeArray = np.abs( self.dicAttributesArray[ piece ][ :, attributeId, 0 ] - self.dicAttributesArray[ piece ][ :, attributeId, 1 ] ) + createAttribute( self.outputMesh, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) else: - attributeArray = np.abs( self.cellsAttributesArray[ :, attributeId, 0 ] - self.cellsAttributesArray[ :, attributeId, 1 ] ) - createAttribute( self.outputMesh, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) - else: - it = self.outputMesh.NewIterator() - while not it.IsDoneWithTraversal(): - dataset: vtkDataSet = vtkDataSet.SafeDownCast( self.outputMesh.GetDataSet( it ) ) - lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) - if piece: - attributeArray = np.abs( self.cellsAttributesArray[ lToG, attributeId, 0 ] - self.cellsAttributesArray[ lToG, attributeId, 1 ] ) - else: - attributeArray = np.abs( self.cellsAttributesArray[ lToG - self.idPointMax, attributeId, 0 ] - self.cellsAttributesArray[ lToG - self.idPointMax, attributeId, 1 ] ) - createAttribute( dataset, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) - it.GoToNextItem() + listBlockId: list[ int ] = getBlockElementIndexesFlatten( self.outputMesh ) + for BlockId in listBlockId: + dataset: vtkDataSet = vtkDataSet.SafeDownCast( self.outputMesh.GetDataSet( BlockId ) ) + lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) + attributeArray = np.abs( self.dicAttributesArray[ piece ][ lToG, attributeId, 0 ] - self.dicAttributesArray[ piece ][ lToG, attributeId, 1 ] ) + createAttribute( dataset, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) return From 47bcd9fdd84d6d89d0bb412301d8e5223a8baf8a Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 18:55:22 +0100 Subject: [PATCH 06/11] fix setMeshes --- .../geos/processing/generic_processing_tools/DiffFields.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py index ddd5f8af3..2abe072e1 100644 --- a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py +++ b/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py @@ -138,11 +138,11 @@ def setMeshes( raise TypeError( "The meshes must have the same type." ) dicMeshesMaxElementId: dict[ bool, list[ int ] ] = { False: [ 0, 0 ], True: [ 0, 0 ], } - if isinstance( self.listMeshes[ 0 ], vtkDataSet ): + if isinstance( listMeshes[ 0 ], vtkDataSet ): for meshId, mesh in enumerate( listMeshes ): for piece in dicMeshesMaxElementId: dicMeshesMaxElementId[ piece ][ meshId ] = np.max( getArrayInObject( mesh, "localToGlobalMap", piece ) ) - elif isinstance( self.listMeshes[ 0 ], vtkMultiBlockDataSet ): + elif isinstance( listMeshes[ 0 ], vtkMultiBlockDataSet ): setDatasetType: set[ str ] = set() for meshId, mesh in enumerate( listMeshes ): listMeshBlockId: list[ int ] = getBlockElementIndexesFlatten( mesh ) From b132c9fd538fb8886e2f6f99cd578cb09f14bf5e Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 18:57:06 +0100 Subject: [PATCH 07/11] Rename the AttributeDiff file --- .../generic_processing_tools/{DiffFields.py => AttributesDiff.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename geos-processing/src/geos/processing/generic_processing_tools/{DiffFields.py => AttributesDiff.py} (100%) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py b/geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py similarity index 100% rename from geos-processing/src/geos/processing/generic_processing_tools/DiffFields.py rename to geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py From 2f8e82e1491f0082e4c83351af524bbb3eec133f Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 18:58:39 +0100 Subject: [PATCH 08/11] Update the doc --- docs/geos_processing_docs/generic_processing_tools.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/geos_processing_docs/generic_processing_tools.rst b/docs/geos_processing_docs/generic_processing_tools.rst index 66bd22d02..261e6dcb9 100644 --- a/docs/geos_processing_docs/generic_processing_tools.rst +++ b/docs/geos_processing_docs/generic_processing_tools.rst @@ -12,6 +12,13 @@ geos.processing.generic_processing_tools.AttributeMapping filter :undoc-members: :show-inheritance: +geos.processing.generic_processing_tools.AttributesDiff filter +---------------------------------------------------------------- + +.. automodule:: geos.processing.generic_processing_tools.AttributesDiff + :members: + :undoc-members: + :show-inheritance: geos.processing.generic_processing_tools.CreateConstantAttributePerRegion filter -------------------------------------------------------------------------------- From debf4ce854d430100501d0786a059768e0c5992a Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 19:03:17 +0100 Subject: [PATCH 09/11] Update the test of AttributesDiff --- geos-processing/tests/test_DiffMesh.py | 36 ++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/geos-processing/tests/test_DiffMesh.py b/geos-processing/tests/test_DiffMesh.py index 640c79a07..f51dee3f2 100644 --- a/geos-processing/tests/test_DiffMesh.py +++ b/geos-processing/tests/test_DiffMesh.py @@ -5,8 +5,40 @@ # ruff: noqa: E402 # disable Module level import not at top of file # mypy: disable-error-code="operator" import pytest +import numpy as np from typing import Any -from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet +from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataSet -from geos.processing.generic_processing_tools.DiffFields import DiffFields +from geos.processing.generic_processing_tools.AttributesDiff import AttributesDiff +from geos.mesh.utils.arrayHelpers import getArrayInObject +from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten + + +@pytest.mark.parametrize( "mesh1Name, mesh2Name", [ +] ) +def test_AttributesDiff( + dataSetTest: Any, + mesh1Name: str, + mesh2Name: str, +) -> None: + """Test the filter AttributesDiff.""" + mesh1: vtkDataSet | vtkMultiBlockDataSet = dataSetTest( mesh1Name ) + mesh2: vtkDataSet | vtkMultiBlockDataSet = dataSetTest( mesh2Name ) + + AttributesDiffFilter: AttributesDiff = AttributesDiff() + AttributesDiffFilter.setMeshes( [ mesh1, mesh2 ] ) + AttributesDiffFilter.logSharedAttributeInfo() + dicAttributesToCompare: dict[ bool, set[ str ] ] = { False: set( [ "elementCenter", "localToGlobalMap" ] ), True: set( [ "localToGlobalMap" ] ) } + AttributesDiffFilter.setDicAttributesToCompare( dicAttributesToCompare ) + AttributesDiffFilter.applyFilter() + mesh: vtkDataSet | vtkMultiBlockDataSet = mesh1.NewInstance() + mesh.ShallowCopy( AttributesDiffFilter.getOutput() ) + dicAttributesDiffNames: dict[ bool, set[ str ] ] = AttributesDiffFilter.getDicAttributesDiffNames() + listFlattenIndexes = getBlockElementIndexesFlatten( mesh ) + for it in listFlattenIndexes: + dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( it ) ) + for piece, listDiffAttributesName in dicAttributesDiffNames.items(): + for diffAttributeName in listDiffAttributesName: + test = getArrayInObject( dataset, diffAttributeName, piece ) + assert ( test == np.zeros( test.shape ) ).all() From 7ce39620cecd694a7cdd5c10e92cea6af101a154 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Tue, 23 Dec 2025 19:04:34 +0100 Subject: [PATCH 10/11] Rename the test file of AttributesDiff filter --- .../tests/{test_DiffMesh.py => test_AttributesDiff.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename geos-processing/tests/{test_DiffMesh.py => test_AttributesDiff.py} (100%) diff --git a/geos-processing/tests/test_DiffMesh.py b/geos-processing/tests/test_AttributesDiff.py similarity index 100% rename from geos-processing/tests/test_DiffMesh.py rename to geos-processing/tests/test_AttributesDiff.py From 6317757d2dc96b936ca0e0e2ab6a96db35e84e45 Mon Sep 17 00:00:00 2001 From: Romain Baville Date: Wed, 24 Dec 2025 14:05:14 +0100 Subject: [PATCH 11/11] Change bool with enum for the attribute localization --- .../AttributesDiff.py | 57 +++++++++---------- geos-processing/tests/test_AttributesDiff.py | 7 ++- 2 files changed, 30 insertions(+), 34 deletions(-) diff --git a/geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py b/geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py index 2abe072e1..1af11aad4 100644 --- a/geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py +++ b/geos-processing/src/geos/processing/generic_processing_tools/AttributesDiff.py @@ -1,6 +1,6 @@ # SPDX-License-Identifier: Apache 2.0 # SPDX-FileCopyrightText: Copyright 2023-2025 TotalEnergies -# SPDX-FileContributor: Jacques Franc +# SPDX-FileContributor: Jacques Franc, Romain Baville # ruff: noqa: E402 # disable Module level import not at top of file import logging @@ -15,6 +15,7 @@ from geos.mesh.utils.arrayHelpers import ( getAttributeSet, getNumberOfComponents, getArrayInObject ) from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten from geos.utils.Logger import ( Logger, getLogger ) +from geos.utils.pieceEnum import Piece __doc__ = """ Attributes Diff is a vtk that compute L1 and L2 differences between attributes shared by two identical meshes. @@ -27,8 +28,8 @@ import logging from vtkmodules.vtkCommonDataModel import vtkMultiBlockDataSet, vtkDataSet - from geos.processing.generic_processing_tools.AttributesDiff import AttributesDiff + from geos.utils.pieceEnum import Piece # Filter inputs: computeL2Diff: bool # defaults to False @@ -49,11 +50,11 @@ attributesDiffFilter.logSharedAttributeInfo() # Get the shared attributes (optional): - dicSharedAttributes: dict[ bool, set[ str ] ] + dicSharedAttributes: dict[ Piece, set[ str ] ] dicSharedAttributes = attributesDiffFilter.getDicSharedAttribute() # Set the attributes to compare: - dicAttributesToCompare: dict[ bool, set[ str ] ] + dicAttributesToCompare: dict[ Piece, set[ str ] ] attributesDiffFilter.setDicAttributesToCompare( dicAttributesToCompare ) # Do calculations: @@ -81,12 +82,12 @@ def __init__( Defaults to False. """ self.listMeshes: list[ vtkMultiBlockDataSet | vtkDataSet ] = [] - self.dicNbElements: dict[ bool, int ] = {} + self.dicNbElements: dict[ Piece, int ] = {} - self.dicSharedAttributes: dict[ bool, set[ str ] ] = {} - self.dicAttributesToCompare: dict[ bool, set[ str ] ] = {} - self.dicAttributesDiffNames: dict[ bool, list[ str ] ] = {} - self.dicAttributesArray: dict[ bool, npt.NDArray[ np.float32 ] ] = {} + self.dicSharedAttributes: dict[ Piece, set[ str ] ] = {} + self.dicAttributesToCompare: dict[ Piece, set[ str ] ] = {} + self.dicAttributesDiffNames: dict[ Piece, list[ str ] ] = {} + self.dicAttributesArray: dict[ Piece, npt.NDArray[ np.float32 ] ] = {} self.computeL2Diff: bool = computeL2Diff @@ -137,7 +138,7 @@ def setMeshes( if listMeshes[ 0 ].GetClassName() != listMeshes[ 1 ].GetClassName(): raise TypeError( "The meshes must have the same type." ) - dicMeshesMaxElementId: dict[ bool, list[ int ] ] = { False: [ 0, 0 ], True: [ 0, 0 ], } + dicMeshesMaxElementId: dict[ Piece, list[ int ] ] = { Piece.CELLS: [ 0, 0 ], Piece.POINTS: [ 0, 0 ], } if isinstance( listMeshes[ 0 ], vtkDataSet ): for meshId, mesh in enumerate( listMeshes ): for piece in dicMeshesMaxElementId: @@ -151,12 +152,6 @@ def setMeshes( dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( meshBlockId ) ) for piece in dicMeshesMaxElementId: dicMeshesMaxElementId[ piece ][ meshId ] = max( dicMeshesMaxElementId[ piece ][ meshId ], np.max( getArrayInObject( dataset, "localToGlobalMap", piece ) ) ) - if piece: - a= dataset.GetNumberOfPoints() - else: - a= dataset.GetNumberOfCells() - if a - len( getArrayInObject( dataset, "localToGlobalMap", piece ) ) != 0: - print("r") if len( setDatasetType ) != 1: raise TypeError( "All datasets of the meshes must have the same type." ) else: @@ -164,11 +159,11 @@ def setMeshes( for piece, listMeshMaxElementId in dicMeshesMaxElementId.items(): if listMeshMaxElementId[ 0 ] != listMeshMaxElementId[ 1 ]: - raise ValueError( f"The total number of { piece } in the meshes must be the same." ) + raise ValueError( f"The total number of { piece.value } in the meshes must be the same." ) self.listMeshes = listMeshes - self.dicNbElements[ True ] = dicMeshesMaxElementId[ True ][ 0 ] + 1 - self.dicNbElements[ False ] = dicMeshesMaxElementId[ False ][ 0 ] + 1 + self.dicNbElements[ Piece.CELLS ] = dicMeshesMaxElementId[ Piece.CELLS ][ 0 ] + 1 + self.dicNbElements[ Piece.POINTS ] = dicMeshesMaxElementId[ Piece.POINTS ][ 0 ] + 1 self.outputMesh = listMeshes[ 0 ].NewInstance() self.outputMesh.ShallowCopy( listMeshes[ 0 ] ) self._computeDicSharedAttributes() @@ -180,7 +175,7 @@ def _computeDicSharedAttributes( self: Self ) -> None: Keys of the dictionary are the attribute localization and the value are the shared attribute per localization. """ - for piece in [ True, False ]: + for piece in [ Piece.POINTS, Piece.CELLS ]: setSharedAttributes: set[ str ] = getAttributeSet( self.listMeshes[ 0 ], piece ) setSharedAttributes.update( getAttributeSet( self.listMeshes[ 1 ], piece ) ) if setSharedAttributes != set(): @@ -188,7 +183,7 @@ def _computeDicSharedAttributes( self: Self ) -> None: return - def getDicSharedAttribute( self: Self ) -> dict[ bool, set[ str ] ]: + def getDicSharedAttribute( self: Self ) -> dict[ Piece, set[ str ] ]: """Getter of the dictionary with the shared attributes per localization.""" return self.dicSharedAttributes @@ -198,15 +193,15 @@ def logSharedAttributeInfo( self: Self ) -> None: self.logger.warning( "The two meshes do not share any attribute." ) else: for piece, sharedAttributes in self.dicSharedAttributes.items(): - self.logger.info( f"Shared attributes on { piece } are { sharedAttributes }." ) + self.logger.info( f"Shared attributes on { piece.value } are { sharedAttributes }." ) return - def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ bool, set[ str ] ] ) -> None: + def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ Piece, set[ str ] ] ) -> None: """Setter of the dictionary with the shared attribute per localization to compare. Args: - dicAttributesToCompare (dict[bool, set[ str ]]): The dictionary with the attributes to compare per localization. + dicAttributesToCompare (dict[Piece, set[str]]): The dictionary with the attributes to compare per localization. Raises: ValueError: At least one attribute to compare is not a shared attribute. @@ -216,9 +211,9 @@ def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ bool, s wrongAttributes: set[ str ] = setSharedAttributesToCompare.difference( self.dicSharedAttributes[ piece ] ) raise ValueError( f"The attributes to compare { wrongAttributes } are not shared attributes.") - dicNbComponents: dict[ bool, int ] = {} - dicAttributesDiffNames: dict[ bool, list[ str ] ] = {} - dicAttributesArray: dict[ bool, npt.NDArray[ np.float32 ] ] = {} + dicNbComponents: dict[ Piece, int ] = {} + dicAttributesDiffNames: dict[ Piece, list[ str ] ] = {} + dicAttributesArray: dict[ Piece, npt.NDArray[ np.float32 ] ] = {} for piece, setSharedAttributesToCompare in dicAttributesToCompare.items(): dicNbComponents[ piece ] = 0 dicAttributesDiffNames[ piece ] = [] @@ -238,11 +233,11 @@ def setDicAttributesToCompare( self: Self, dicAttributesToCompare: dict[ bool, s return - def getDicAttributesToCompare( self: Self ) -> dict[ bool, set[ str ] ]: + def getDicAttributesToCompare( self: Self ) -> dict[ Piece, set[ str ] ]: """Getter of the dictionary of the attribute to compare per localization.""" return self.dicAttributesToCompare - def getDicAttributesDiffNames( self: Self ) -> dict[ bool, list[ str ] ]: + def getDicAttributesDiffNames( self: Self ) -> dict[ Piece, list[ str ] ]: """Getter of the dictionary with the name of the attribute created with the calculated attributes diff.""" return self.dicAttributesDiffNames @@ -297,14 +292,14 @@ def _computeL1( self: Self ) -> None: attributeArray: npt.NDArray[ Any ] if isinstance( self.listMeshes[ 0 ], vtkDataSet ): attributeArray = np.abs( self.dicAttributesArray[ piece ][ :, attributeId, 0 ] - self.dicAttributesArray[ piece ][ :, attributeId, 1 ] ) - createAttribute( self.outputMesh, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) + createAttribute( self.outputMesh, attributeArray, attributeDiffName, piece=piece, logger=self.logger ) else: listBlockId: list[ int ] = getBlockElementIndexesFlatten( self.outputMesh ) for BlockId in listBlockId: dataset: vtkDataSet = vtkDataSet.SafeDownCast( self.outputMesh.GetDataSet( BlockId ) ) lToG: npt.NDArray[ Any ] = getArrayInObject( dataset, "localToGlobalMap", piece ) attributeArray = np.abs( self.dicAttributesArray[ piece ][ lToG, attributeId, 0 ] - self.dicAttributesArray[ piece ][ lToG, attributeId, 1 ] ) - createAttribute( dataset, attributeArray, attributeDiffName, onPoints=piece, logger=self.logger ) + createAttribute( dataset, attributeArray, attributeDiffName, piece=piece, logger=self.logger ) return diff --git a/geos-processing/tests/test_AttributesDiff.py b/geos-processing/tests/test_AttributesDiff.py index f51dee3f2..f1c873e64 100644 --- a/geos-processing/tests/test_AttributesDiff.py +++ b/geos-processing/tests/test_AttributesDiff.py @@ -13,8 +13,9 @@ from geos.processing.generic_processing_tools.AttributesDiff import AttributesDiff from geos.mesh.utils.arrayHelpers import getArrayInObject from geos.mesh.utils.multiblockHelpers import getBlockElementIndexesFlatten +from geos.utils.pieceEnum import Piece - +# TODO: Create meshes for test @pytest.mark.parametrize( "mesh1Name, mesh2Name", [ ] ) def test_AttributesDiff( @@ -29,12 +30,12 @@ def test_AttributesDiff( AttributesDiffFilter: AttributesDiff = AttributesDiff() AttributesDiffFilter.setMeshes( [ mesh1, mesh2 ] ) AttributesDiffFilter.logSharedAttributeInfo() - dicAttributesToCompare: dict[ bool, set[ str ] ] = { False: set( [ "elementCenter", "localToGlobalMap" ] ), True: set( [ "localToGlobalMap" ] ) } + dicAttributesToCompare: dict[ Piece, set[ str ] ] = { Piece.CELLS: set( [ "elementCenter", "localToGlobalMap" ] ), Piece.POINTS: set( [ "localToGlobalMap" ] ) } AttributesDiffFilter.setDicAttributesToCompare( dicAttributesToCompare ) AttributesDiffFilter.applyFilter() mesh: vtkDataSet | vtkMultiBlockDataSet = mesh1.NewInstance() mesh.ShallowCopy( AttributesDiffFilter.getOutput() ) - dicAttributesDiffNames: dict[ bool, set[ str ] ] = AttributesDiffFilter.getDicAttributesDiffNames() + dicAttributesDiffNames: dict[ Piece, set[ str ] ] = AttributesDiffFilter.getDicAttributesDiffNames() listFlattenIndexes = getBlockElementIndexesFlatten( mesh ) for it in listFlattenIndexes: dataset: vtkDataSet = vtkDataSet.SafeDownCast( mesh.GetDataSet( it ) )