fix nondimensionals

ntrfc/math/vectorcalc.py

@@ -209,6 +209,16 @@ def ellipsoidVol(sig):
 
 
 def vecProjection(direction, vector):
+    """
+    Calculate the projection of a vector onto a direction vector.
+
+    Parameters:
+    direction (list or numpy array): The direction vector onto which the projection will be calculated.
+    vector (list or numpy array): The vector to be projected onto the direction vector.
+
+    Returns:
+    projection (numpy array): A vector representing the projection of the input vector onto the direction vector.
+    """
     unitDir = vecDir(direction)
     return np.dot(vector, unitDir) * unitDir
 

ntrfc/meshquality/nondimensionals.py

@@ -24,8 +24,8 @@ from ntrfc.math.vectorcalc import vecAbs, vecAbs_list, vecProjection, vecAngle,
 
 def cellDirections(cellUMean, wallNorm):
     x = unitvec(cellUMean)  # mainDirection
-    z = unitvec(wallNorm)  # tangential
-    y = unitvec(np.cross(x, z))  # spanwise
+    y = unitvec(wallNorm)  # tangential
+    z = unitvec(np.cross(x, y))  # spanwise
     return np.array([x, y, z])
 
 
@@ -67,20 +67,28 @@ def cellSpans(solutionMesh, calcFrom):
     return cell_spans
 
 
-def getWalluTaus(solutionMesh, surfaceMesh, mu_0, rhofieldname, velfieldname):
-    surfaceNormals = surfaceMesh.extract_surface().compute_normals()
-    cell_centers = solutionMesh["cellCenters"]
 
-    relatedpointonwallids = [surfaceNormals.find_closest_point(i) for i in cell_centers]
-    relatedpointonwall = [surfaceNormals.points[i] for i in relatedpointonwallids]
-    relatedpointonwall_samples = pv.PolyData(relatedpointonwall).sample(surfaceMesh)
-    graduwall = relatedpointonwall_samples["gradient"].reshape(relatedpointonwall_samples.number_of_points, 3, 3)
-    flowdirections = unitvec_list(solutionMesh[velfieldname])
-    gradUyWs = vecAbs_list(
-        [np.dot(flowdirections[i], graduwall[i]) for i in range(relatedpointonwall_samples.number_of_points)])
-    rhoWs = relatedpointonwall_samples[rhofieldname]
+def getWalluTaus(solutionMesh, mu_0, rhofieldname, velfieldname):
+    """
+    Calculate the wall shear stress at various points on a surface in a CFD simulation.
+
+    Parameters:
+    solutionMesh (pv.PolyData): A mesh that represents the solution of the CFD simulation.
+    surfaceMesh (pv.PolyData): A mesh that represents the surface on which the wall shear stress is to be calculated.
+    mu_0 (float): The dynamic viscosity of the fluid.
+    rhofieldname (str): The name of the field that contains the density of the fluid.
+    velfieldname (str): The name of the field that contains the velocity of the fluid.
+
+    Returns:
+    u_taus (np.ndarray): An array containing the wall shear stress at each point on the surface.
+    """
 
-    tauWs = gradUyWs * mu_0 * rhoWs
+
+    graduwall = solutionMesh[f"grad_{velfieldname}"].reshape(solutionMesh.number_of_cells,3,3)
+    wallnormals = solutionMesh["wallNormal"]
+    dudy = vecAbs_list([np.dot(graduwall[i], wallnormals[i]) for i in range(solutionMesh.number_of_cells)])
+    rhoWs = solutionMesh[rhofieldname]
+    tauWs = dudy * mu_0
     u_taus = (tauWs / rhoWs) ** 0.5
 
     return u_taus
@@ -115,6 +123,13 @@ def constructWallMesh(surfaces):
 #                                                                               #
 #################################################################################
 
+
+def compute_scalar_gradient(mesh,arrayname):
+    mesh = mesh.compute_derivative(scalars=arrayname)
+    mesh[f"grad_{arrayname}"] = mesh["gradient"].reshape((mesh.number_of_cells,3,3))
+    return mesh
+
+
 def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield, mu_0):
     """
     :param volmesh: pyvista-vtk object
@@ -127,31 +142,27 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
 
     print("constructing surfacemesh from wall meshes ...")
     surfaceMesh = constructWallMesh(surfaces)
-    surfaceMeshcopy = surfaceMesh.copy()
-    surfaceMeshcopy.clear_data()
-
-    print("preparing processData from meshes")
-
-    volmesh = volmesh.compute_derivative(scalars=use_velfield)
-    volmesh_walladjacent = volmesh.extract_cells(volmesh.find_containing_cell(volmesh.extract_surface().points))
-    surfaceMeshcopy = surfaceMeshcopy.sample(volmesh_walladjacent)
-    surfaceMesh["gradient"] = surfaceMeshcopy["gradient"]
-    volmesh_walladjacent["cellCenters"] = volmesh_walladjacent.cell_centers().points
-
     print("calculating wall-normal vectors...")
     surfacenormals_surface = surfaceMesh.extract_surface().compute_normals()
 
+    print("finding walladjacent cells")
+    volmesh = compute_scalar_gradient(volmesh,use_velfield)
+    walladjacentids = volmesh.find_containing_cell(surfacenormals_surface.points)
+    volmesh_walladjacent = volmesh.extract_cells(walladjacentids)
+
+    volmesh_walladjacent["cellCenters"] = volmesh_walladjacent.cell_centers().points
     volmesh_walladjacent["wallNormal"] = [
         surfacenormals_surface.point_data["Normals"][surfacenormals_surface.find_closest_point(i)]
         for i in volmesh_walladjacent.points]
-    print("calculating cell spans from WallNormals and CellEdges...")
+
+    print("calculating cell spans from WallNormals...")
     spanS = cellSpans(volmesh_walladjacent, use_velfield)
     volmesh_walladjacent["xSpan"] = np.array([i[0] for i in spanS])  # calculate cell span in flow direction
     volmesh_walladjacent["ySpan"] = np.array([i[1] for i in spanS])  # calculate cell span in wall normal direction
     volmesh_walladjacent["zSpan"] = np.array([i[2] for i in spanS])  # calculate cell span in span direction
 
     print("calculating wall-shear and friction-velocity")
-    volmesh_walladjacent["uTaus"] = getWalluTaus(volmesh_walladjacent, surfaceMesh, mu_0, use_rhofield, use_velfield)
+    volmesh_walladjacent["uTaus"] = getWalluTaus(volmesh_walladjacent, mu_0, use_rhofield, use_velfield)
 
     print("calculating grid spacing")
     gridSpacings = gridSpacing(mu_0, volmesh_walladjacent)
@@ -161,3 +172,4 @@ def calc_dimensionless_gridspacing(volmesh, surfaces, use_velfield, use_rhofield
     volmesh_walladjacent["DeltaZPlus"] = gridSpacings[2]
 
     return volmesh_walladjacent
+

tests/meshquality/test_ntrfc_nondimensionals.py

@@ -20,9 +20,9 @@ def test_calc_dimensionless_gridspacing():
         velocity = 2
         gradient = velocity / height
         # analytical solution
-        span_x = length / (nodes - 1)
-        span_y = width / (nodes - 1)
-        span_z = height / (nodes - 1)
+        span_x = width / (nodes - 1)
+        span_y = height / (nodes - 1)
+        span_z = length / (nodes - 1)
         wallshearstress = mu_0 * gradient
         wallshearvelocity = np.sqrt(wallshearstress / rho)
         deltaxplus = wallshearvelocity * span_x / mu_0
@@ -36,7 +36,13 @@ def test_calc_dimensionless_gridspacing():
         grid = pv.StructuredGrid(x, y, z)
         # scale the mesh
         grid.points /= nodes - 1
-        grid.points *= np.array([length, height, width])
+        grid.points *= np.array([width, height, length])
+
+        # define velocityfield
+        bounds = grid.bounds
+        min_z = bounds[4]
+        grid["U"] = [gradient * (grid.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(grid.number_of_cells)]
+        grid["rho"] = np.ones(grid.number_of_cells)
         # extract surface
         surface = grid.extract_surface()
         facecellids = surface.surface_indices()
@@ -50,19 +56,16 @@ def test_calc_dimensionless_gridspacing():
                 upperwallids.append(faceid)
         lowerwall = surface.extract_cells(lowerwallids)
         upperwall = surface.extract_cells(upperwallids)
-        # define velocityfield
-        grid["U"] = np.array([gradient, 0, 0]) * grid.cell_centers().points
-        lowerwall["U"] = np.array([gradient, 0, 0]) * lowerwall.cell_centers().points
-        upperwall["U"] = np.array([gradient, 0, 0]) * upperwall.cell_centers().points
-        grid["rho"] = np.ones(grid.number_of_cells)
+        lowerwall["U"] =  [gradient * (lowerwall.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(lowerwall.number_of_cells)]
+        upperwall["U"] = [gradient * (upperwall.cell_centers().points[::,1][i]-min_z) * np.array([0,0,1]) for i in range(upperwall.number_of_cells)]
         lowerwall["rho"] = np.ones(lowerwall.number_of_cells) * rho
         upperwall["rho"] = np.ones(upperwall.number_of_cells) * rho
         dimless_gridspacings = calc_dimensionless_gridspacing(grid, [lowerwall, upperwall], "U", "rho", mu_0)
-        assert all(np.isclose(deltaxplus, dimless_gridspacings[
+        assert all(np.isclose(deltazplus, dimless_gridspacings[
             "DeltaXPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in x direction not accurate"
         assert all(np.isclose(deltayplus, dimless_gridspacings[
             "DeltaYPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in y direction not accurate"
-        assert all(np.isclose(deltazplus, dimless_gridspacings[
+        assert all(np.isclose(deltaxplus, dimless_gridspacings[
             "DeltaZPlus"], rtol=0.05)), "calc_dimensionelss_gridspcing in z direction not accurate"
 
     runtest(height=1, length=1, width=1)