Add the shift of points as a vector field.

programming/mesh/cgns/add_cell_quality.py

@@ -41,17 +41,8 @@ def measure(mesh_file: str) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
 
     # read the cgns file
     cgns, zone, zone_size = wrapper.getUniqueZone(args.mesh)
-    # get the existing FlowSolution_t<CellCenter> or create a new one if not found
-    cell_data = None
-    solutions = wrapper.getChildrenByType(zone, 'FlowSolution_t')
-    for solution in solutions:
-        location = wrapper.getUniqueChildByType(solution, 'GridLocation_t')
-        location_str = wrapper.arr2str(wrapper.getNodeData(location))
-        if location_str == 'CellCenter':
-            cell_data = solution
-    # add the cell qualities as fields of cell data
-    if cell_data is None:
-        cell_data = cgl.newFlowSolution(zone, 'CellQualities', 'CellCenter')
+    cell_data = wrapper.getSolutionByLocation(zone,
+        'CellCenter', 'CellQualities')
     cgl.newDataArray(cell_data, 'MinDetJac', min_det_jac)
     cgl.newDataArray(cell_data, 'MaxDetJac', max_det_jac)
     cgl.newDataArray(cell_data, 'Volume', volume)

programming/mesh/cgns/shift_points.py → programming/mesh/cgns/add_point_shift.py

@@ -13,7 +13,7 @@
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(
         prog = f'python3 {sys.argv[0]}',
-        description='Shift the points in a given mesh to their given coordinates.')
+        description='Add the shift of each point as a vector field.')
     parser.add_argument('--folder', type=str, help='the working folder containing input and output files')
     parser.add_argument('--input', type=str, help='the CGNS file of the mesh to be shifted')
     parser.add_argument('--cad_points', type=str, help='the NPY file of points on CAD surfaces')
@@ -29,17 +29,10 @@
     mesh_zone_size = wrapper.getNodeData(mesh_zone)
     n_node = mesh_zone_size[0][0]
     n_cell = mesh_zone_size[0][1]
-    if args.verbose:
-        print(f'n_node = {n_node}, n_cell = {n_cell}')
-    mesh_points = wrapper.getChildrenByType(
-        wrapper.getUniqueChildByType(mesh_zone, 'GridCoordinates_t'), 'DataArray_t')
-    mesh_points_x, mesh_points_y, mesh_points_z = mesh_points[X][1], mesh_points[Y][1], mesh_points[Z][1]
-    assert n_node == len(mesh_points_x) == len(mesh_points_y) == len(mesh_points_z)
+    print(f'n_node = {n_node}, n_cell = {n_cell}')
 
-    # backup the original coords
-    old_points_x = np.array(mesh_points_x[:])
-    old_points_y = np.array(mesh_points_y[:])
-    old_points_z = np.array(mesh_points_z[:])
+    # read the original coords
+    _, old_x, old_y, old_z = wrapper.readPoints(mesh_zone, mesh_zone_size)
 
     # load the shifted coords
     cad_points = np.load(args.cad_points)
@@ -48,25 +41,28 @@
     stl_points = np.load(args.stl_points)
     assert (n_node, 3) == stl_points.shape
 
-    # update the coords
+    # get the new coords
+    new_x = np.ndarray(old_x.shape)
+    new_y = np.ndarray(old_y.shape)
+    new_z = np.ndarray(old_z.shape)
     for i_node in range(n_node):
         if args.verbose:
             print(f'shifting node {i_node} / {n_node}')
-        x = mesh_points_x[i_node]
+        x = old_x[i_node]
         if x < args.x_cut:
-            mesh_points_x[i_node] = stl_points[i_node, X]
-            mesh_points_y[i_node] = stl_points[i_node, Y]
-            mesh_points_z[i_node] = stl_points[i_node, Z]
+            new_x[i_node] = stl_points[i_node, X]
+            new_y[i_node] = stl_points[i_node, Y]
+            new_z[i_node] = stl_points[i_node, Z]
         else:
-            mesh_points_x[i_node] = cad_points[i_node, X]
-            mesh_points_y[i_node] = cad_points[i_node, Y]
-            mesh_points_z[i_node] = cad_points[i_node, Z]
+            new_x[i_node] = cad_points[i_node, X]
+            new_y[i_node] = cad_points[i_node, Y]
+            new_z[i_node] = cad_points[i_node, Z]
 
     # write the shifts as a field of point data
-    point_data = cgl.newFlowSolution(mesh_zone, 'FlowSolutionCellHelper', 'Vertex')
-    cgl.newDataArray(point_data, 'ShiftX', mesh_points_x - old_points_x)
-    cgl.newDataArray(point_data, 'ShiftY', mesh_points_y - old_points_y)
-    cgl.newDataArray(point_data, 'ShiftZ', mesh_points_z - old_points_z)
+    point_data = wrapper.getSolutionByLocation(mesh_zone, 'Vertex', 'PointShift')
+    cgl.newDataArray(point_data, 'ShiftX', new_x - old_x)
+    cgl.newDataArray(point_data, 'ShiftY', new_y - old_y)
+    cgl.newDataArray(point_data, 'ShiftZ', new_z - old_z)
 
     # write the shifted mesh
     output = f'{args.folder}/shifted.cgns'

programming/mesh/cgns/wrapper.py

@@ -113,5 +113,22 @@ def removeSolutionsByLocation(zone, location):
             cgu.removeChildByName(zone, getNodeName(solution))
 
 
+def getSolutionByLocation(zone, location_given, solution_name):
+    """Get an existing FlowSolution_t with a given GridLocation_t or create a new one if not found.
+    """
+    the_solution = None
+    solutions = getChildrenByType(zone, 'FlowSolution_t')
+    for solution in solutions:
+        location = getUniqueChildByType(solution, 'GridLocation_t')
+        location_found = arr2str(getNodeData(location))
+        if location_found == location_given:
+            the_solution = solution
+            break
+    # add the cell qualities as fields of cell data
+    if the_solution is None:
+        the_solution = cgl.newFlowSolution(zone, solution_name, location_given)
+    return the_solution
+
+
 if __name__ == '__main__':
     printInfo(sys.argv[1])