fixes and progress

examples/gwk_compressor/Snakefile

 from snakemake.utils import validate
-import pyvista as pv
 from ntrfc.utils.filehandling.read_datafiles import yaml_dict_read
 from ntrfc.utils.filehandling.read_mesh import load_mesh
+import numpy as np
 
 #todo move this snakepipe to according dir in sourcecode
 
@@ -9,23 +9,25 @@ configfile : "casesettings.yaml"
 validate(config,"config.schema.yaml")
 basedir = workflow.current_basedir
 
+
 rule profile_from_pointcloud:
     input:
         "00_resources/profile_pointcloud.txt"
     output:
-        "01_geometry/sortedPoints.txt", "01_geometry/psPoly.vtk", "01_geometry/ssPoly.vtk", "01_geometry/midsPoly.vtk", "01_geometry/geometry_paras.yml"
+        "01_profile/sortedPoints.txt", "01_profile/psPoly.vtk", "01_profile/ssPoly.vtk", "01_profile/midsPoly.vtk", "01_profile/geometry_paras.yml",
+        "01_profile/profile.pdf"
     run:
         from ntrfc.preprocessing.geometry_creation.profile_pointcloud import generate_profile_pointcloud_geometry
 
         generate_profile_pointcloud_geometry(config["blade"],basedir)
 
+"""
 
 rule create_2d_domain:
     input:
-        "01_geometry/sortedPoints.txt", "01_geometry/psPoly.vtk", "01_geometry/ssPoly.vtk", "01_geometry/midsPoly.vtk", "01_geometry/geometry_paras.yml"
+        "01_profile/sortedPoints.txt", "01_profile/psPoly.vtk", "01_profile/ssPoly.vtk", "01_profile/midsPoly.vtk", "01_profile/geometry_paras.yml"
     output:
-        "02_domain/inlet_2d.vtk","02_domain/outlet_2d.vtk","02_domain/y_upper_2d.vtk","02_domain/y_lower_2d.vtk",
-        "02_domain/ogridline_2d.vtk"
+        "02_boundaries/inlet_2d.vtk","02_boundaries/outlet_2d.vtk","02_boundaries/y_upper_2d.vtk","02_boundaries/y_lower_2d.vtk"
     run:
         from ntrfc.preprocessing.mesh_creation.domain_creation import create_2d_domain
 
@@ -40,9 +42,37 @@ rule create_2d_domain:
         midsPoly = load_mesh(midsPoly_fpath)
         geometry_paras = yaml_dict_read(geometry_paras_fpath)
 
-        create_2d_domain(config["blade"],basedir,midsPoly,ssPoly,psPoly)
+        create_2d_domain(config["geometry"],basedir,midsPoly,ssPoly,psPoly,geometry_paras)
+"""
+"""
+rule create_blockboundaries:
+    input:
+        "01_profile/sortedPoints.txt", "01_profile/geometry_paras.yml" ,"01_profile/midsPoly.vtk" ,"02_boundaries/inlet_2d.vtk",
+        "02_boundaries/outlet_2d.vtk", "02_boundaries/y_upper_2d.vtk", "02_boundaries/y_lower_2d.vtk"
+
+    output:
+        "02_domain/inlet_2d.vtk","02_domain/outlet_2d.vtk","02_domain/y_upper_2d.vtk","02_domain/y_lower_2d.vtk",
+        "02_domain/ogridline_2d.vtk"
+    run:
+        from ntrfc.preprocessing.mesh_creation.domain_creation import create_2d_blocklines
+        sortedPoints_fpath = input[0]
+        geometry_paras_fpath = input[1]
+        midsPoly_fpath = input[2]
+        inlet_fpath = input[3]
+        outlet_fpath = input[4]
+        y_upper_fpath = input[5]
+        y_lower_paras_fpath = input[6]
 
+        sortedPoints = np.loadtxt(os.path.join(basedir,sortedPoints_fpath))
+        geometry_paras = yaml_dict_read(os.path.join(basedir,geometry_paras_fpath))
+        midsPoly = load_mesh(os.path.join(basedir,midsPoly_fpath))
+        inlet = load_mesh(os.path.join(basedir,inlet_fpath))
+        outlet = load_mesh(os.path.join(basedir,outlet_fpath))
+        y_upper = load_mesh(os.path.join(basedir,y_upper_fpath))
+        y_lower = load_mesh(os.path.join(basedir,y_lower_paras_fpath))
 
+        create_2d_blocklines(config["blade"],basedir,sortedPoints,geometry_paras,midsPoly,inlet,outlet,y_upper,y_lower)
+"""
 """
 rule profile_from_naca:
 #    input:

examples/gwk_compressor/casesettings.yaml

@@ -8,4 +8,4 @@ geometry :
     pitch : 0.0765                                #Pitch [m]
     x_inlet : -0.13                                 #Shift Inlet along x [m]
     x_outlet : 0.26                                  #Shift Outlet along x [m]
-    shift_domain : 0.015                                 #Shift Domain to have probes inside domain [m]
+    blade_shift : 0.015                                 #Shift Domain to have probes inside domain [m]

ntrfc/database/naca_airfoil_creator.py

@@ -256,6 +256,9 @@ def naca(number, n, finite_TE = False, half_cosine_spacing = True):
     else:
         raise Exception
 
+    #delete last point as it is defined twice
+    X = X[:-1]
+    Y = Y[:-1]
     """
     fix for trailing_edge finite_TE
     """

ntrfc/preprocessing/geometry_creation/profile_pointcloud.py

@@ -4,6 +4,7 @@ from ntrfc.utils.geometry.pointcloud_methods import extract_geo_paras
 
 import numpy as np
 import os
+import matplotlib.pyplot as plt
 
 def generate_profile_pointcloud_geometry(settings,basedir):
     ptcloud_profile = settings["ptcloud_profile"]
@@ -36,9 +37,19 @@ def generate_profile_pointcloud_geometry(settings,basedir):
                     "camber_angle":float(camber_angle)}
 
 
-    geo_dir = os.path.join(basedir,"01_geometry")
+    geo_dir = os.path.join(basedir,"01_profile")
     np.savetxt(os.path.join(geo_dir,"sortedPoints.txt"), sortedPoints)
     psPoly.save(os.path.join(geo_dir,"psPoly.vtk"),False)
     ssPoly.save(os.path.join(geo_dir,"ssPoly.vtk"),False)
     midsPoly.save(os.path.join(geo_dir,"midsPoly.vtk"),False)
     write_yaml_dict(os.path.join(geo_dir,"geometry_paras.yml"),geometry_paras)
+
+    #todo: here, we should also plot the geometry-parameters (chord, angles, ...)
+    plt.figure()
+    plt.plot(psPoly.points[::,0],psPoly.points[::,1], color="#6c3376")
+    plt.plot(ssPoly.points[::,0],ssPoly.points[::,1], color="#FF2211")
+    plt.plot(midsPoly.points[::,0],midsPoly.points[::,1], color="#FF22CC")
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.savefig(os.path.join(geo_dir,'profile.pdf'))
+

ntrfc/preprocessing/mesh_creation/domain_creation.py

 from ntrfc.utils.geometry.pointcloud_methods import calcMidPassageStreamLine
+from ntrfc.utils.pyvista_utils.line import lines_from_points
+import numpy as np
+import pyvista as pv
+import os
 
-def create_2d_domain(settings,basedir,midsPoly,ssPoly,psPoly):
-
+def create_2d_domain(geosettings, basedir, midsPoly, ssPoly, psPoly, geometry_paras):
+    beta_meta_01 = geometry_paras["beta_meta_01"]
+    beta_meta_02 = geometry_paras["beta_meta_02"]
+    x_inlet = geosettings["x_inlet"]
+    x_outlet = geosettings["x_outlet"]
+    pitch = geosettings["pitch"]
+    blade_shift = geosettings["blade_shift"]
     x_mids = midsPoly.points[::, 0]
     y_mids = midsPoly.points[::, 1]
     x_ss = ssPoly.points[::, 0]
@@ -31,3 +40,30 @@ def create_2d_domain(settings,basedir,midsPoly,ssPoly,psPoly):
                            per_y_upper.points[per_y_upper.points[::, 0].argmax()]])
 
     outletPoly = pv.Line(*outlet_pts)
+
+    domain_dir = os.path.join(basedir,"02_boundaries")
+    inletPoly.save(os.path.join(domain_dir,"inlet_2d.vtk"),False)
+    outletPoly.save(os.path.join(domain_dir,"outlet_2d.vtk"),False)
+    per_y_upper.save(os.path.join(domain_dir,"y_upper_2d.vtk"),False)
+    per_y_lower.save(os.path.join(domain_dir,"y_lower_2d.vtk"),False)
+
+def create_2d_blocklines(geo_settings,basedir,sortedPoints,geometry_paras,midsPoly,inlet,outlet,y_upper,y_lower):
+    bladeline = lines_from_points(sortedPoints)
+    bladehelpersurface = bladeline.extrude((0, 0, 0.01))
+    ogridhelpersurface = bladehelpersurface.compute_normals()
+    #todo: replace .002 with a propper generic scaling (using the length of the profile and the distance to a periodic as reference)
+    ogridhelpersurface.points += .002 * ogridhelpersurface.point_normals
+    ogridline=ogridhelpersurface.slice(normal="z")
+    ogridline.points[:, 2] = 0
+
+    p=pv.Plotter()
+    p.add_mesh(inlet)
+    p.add_mesh(inlet)
+    p.add_mesh(y_upper)
+    p.add_mesh(y_lower)
+    p.add_mesh(bladeline)
+
+    p.add_mesh(ogridline)
+
+    p.show()
+    return 0

ntrfc/utils/geometry/pointcloud_methods.py

@@ -178,7 +178,7 @@ def extract_vk_hk(sortedPoly, verbose=False):
         midsPoly = midline_from_sides(ind_1, ind_2, sortedPoly.points, psPoly, ssPoly)
         arclengths = midsPoly.compute_arc_length()["arc_length"]
         midslength = sum(arclengths)
-        return midslength
+        return midsPoly.length
 
     xs, ys = sortedPoly.points[::, 0], sortedPoly.points[::, 1]
     ind_1, ind_2 = calc_largedistant_idx(xs, ys)
@@ -285,9 +285,9 @@ def extract_geo_paras(points, alpha, verbose=False):
     vk_tangent = np.stack((xmids[0] - xmids[1], ymids[0] - ymids[1], 0)).T
     hk_tangent = np.stack((xmids[-2] - xmids[-1], ymids[-2] - ymids[-1], 0)).T
     camber = np.stack((xmids[0] - xmids[-1], ymids[0] - ymids[-1], 0)).T[::-1]
-    beta_leading = vecAngle(vk_tangent, np.array([0, 1, 0])) / np.pi * 180
-    beta_trailing = vecAngle(hk_tangent, np.array([0, 1, 0])) / np.pi * 180
-    camber_angle = vecAngle(camber, np.array([0, 1, 0])) / np.pi * 180
+    beta_leading = vecAngle(vk_tangent, np.array([1, 0, 0])) / np.pi * 180
+    beta_trailing = vecAngle(hk_tangent, np.array([1, 0, 0])) / np.pi * 180
+    camber_angle = vecAngle(camber, np.array([1, 0, 0])) / np.pi * 180
 
     if verbose:
         p = pv.Plotter()

tests/test_ntrfc_geometry.py

@@ -85,7 +85,9 @@ def test_extract_vk_hk(verbose=False):
     # d1,d2,d3,d4 = np.random.randint(0,9),np.random.randint(0,9),np.random.randint(0,9),np.random.randint(0,9)
     # digitstring = str(d1)+str(d2)+str(d3)+str(d4)
     # manifold problems with other profiles with veronoi-mid and other unoptimized code. therefor tests only 0009
-    X, Y = naca("6409", res, finite_TE=False, half_cosine_spacing=True)
+    #todo: currently we cant test half_cosine_spacing profiles, as the resolution is too good for extract_vk_hk, this must be changed!
+    #todo: we should test finite_TE profiles with round TE. the current implementation is not validated for this
+    X, Y = naca("6506", res, finite_TE=False, half_cosine_spacing=False)
     ind_hk_test = 0
     ind_vk_test = res