import numpy as np input0 = self.GetInput()points = input0.GetPoints()numPts = points.GetNumberOfPoints() gradU = input0.GetPointData().GetArray("U_gradient")U = input0.GetPointData().GetVectors("U") from vtkmodules.vtkCommonDataModel import vtkPolyDatafrom vtkmodules.vtkCommonCore import vtkPoints, vtkFloatArray, vtkIntArray outputPts = vtkPoints() # 出力用属性typeArray = vtkIntArray()typeArray.SetName("CriticalPointType") vortArray = vtkFloatArray()vortArray.SetName("vorticity") divArray = vtkFloatArray()divArray.SetName("divergence") # しきい値設定vort_eps = 1e-2 # 渦度がこれ以上ならフォーカスdiv_eps = 1e-2 # 発散のしきい値speed_eps = 1e-3 # 停留点とみなす速度 # ラベル定義LABEL_SINK = -1LABEL_SADDLE = 0LABEL_SOURCE = 1LABEL_FOCUS = 2 for i in range(numPts): g = gradU.GetTuple9(i) u = np.array(U.GetTuple3(i)) speed = np.linalg.norm(u) dudx = g[0] dudy = g[1] dvdx = g[3] dvdy = g[4] # 2D発散と渦度 div = dudx + dvdy vort = dvdx - dudy # 初期分類 label = None if abs(vort) > vort_eps and abs(div) < div_eps: label = LABEL_FOCUS elif abs(vort) < vort_eps and div > div_eps: label = LABEL_SOURCE elif abs(vort) < vort_eps and div < -div_eps: label = LABEL_SINK elif abs(div) < div_eps and abs(vort) < vort_eps and speed < speed_eps: label = LABEL_SADDLE else: continue # 特異点でないと判断 outputPts.InsertNextPoint(points.GetPoint(i)) typeArray.InsertNextValue(label) divArray.InsertNextValue(div) vortArray.InsertNextValue(vort) # 出力構築output = self.GetOutput()output.SetPoints(outputPts)output.GetPointData().AddArray(typeArray)output.GetPointData().AddArray(divArray)output.GetPointData().AddArray(vortArray)