progressing meshing-implementation

examples/gwk_compressor/Snakefile

 from snakemake.utils import validate
-from ntrfc.utils.filehandling.read_datafiles import yaml_dict_read
-from ntrfc.utils.filehandling.read_mesh import load_mesh
+from ntrfc.utils.filehandling.datafiles import yaml_dict_read
+from ntrfc.utils.filehandling.mesh import load_mesh
 import numpy as np
 
 #todo move this snakepipe to according dir in sourcecode
@@ -20,39 +20,65 @@ rule profile_from_pointcloud:
 
         generate_profile_pointcloud_geometry(config["blade"],basedir)
 """
+
+"""
+rule dumdidum:
+
+    #snakemake Wasd1193.csv
+    #1193 --> step für das Erzeugen von parameter (float, etc.)
+
+    output:
+        "Wasd{parameter_x}.csv"
+    run:
+        print(wildcards.parameter_x)
 """
+
 rule create_2d_domain:
     input:
-        "01_profile/sortedPoints.txt", "01_profile/psPoly.vtk", "01_profile/ssPoly.vtk", "01_profile/midsPoly.vtk", "01_profile/geometry_paras.yml"
+#        sortedPoints_fpath="01_profile/sortedPoints.txt",
+        psPoly_fpath="01_profile/psPoly.vtk",
+        ssPoly_fpath="01_profile/ssPoly.vtk",
+        midsPoly_fpath="01_profile/midsPoly.vtk",
+        geometry_paras_fpath="01_profile/geometry_paras.yml"
     output:
-        "02_boundaries/inlet_2d.vtk","02_boundaries/outlet_2d.vtk","02_boundaries/y_upper_2d.vtk","02_boundaries/y_lower_2d.vtk",
-        "02_boundaries/domain.pdf"
+        "02_domainboundaries/inlet_2d.vtk","02_domainboundaries/outlet_2d.vtk",
+        "02_domainboundaries/y_upper_2d.vtk","02_domainboundaries/y_lower_2d.vtk",
+        "02_domainboundaries/domain.pdf"
+    params:
+        geometry=config["geometry"]
     run:
         from ntrfc.preprocessing.mesh_creation.domain_creation import create_2d_domain
 
-        sortedPoints_fpath = input[0]
-        psPoly_fpath = input[1]
-        ssPoly_fpath = input[2]
-        midsPoly_fpath = input[3]
-        geometry_paras_fpath = input[4]
+        #todo: dateieinlesen nicht im rule, sondern in der funktion
+        psPoly = load_mesh(input.psPoly_fpath)
+        ssPoly = load_mesh(input.ssPoly_fpath)
+        midsPoly = load_mesh(input.midsPoly_fpath)
+        geometry_paras = yaml_dict_read(input.geometry_paras_fpath)
 
-        psPoly = load_mesh(psPoly_fpath)
-        ssPoly = load_mesh(ssPoly_fpath)
-        midsPoly = load_mesh(midsPoly_fpath)
-        geometry_paras = yaml_dict_read(geometry_paras_fpath)
+        create_2d_domain(params["geometry"],basedir,midsPoly,ssPoly,psPoly,geometry_paras)
 
-        create_2d_domain(config["geometry"],basedir,midsPoly,ssPoly,psPoly,geometry_paras)
 """
+rule create_cascademesh:
+    input:
+    output:
+    run:
+"""
+
 rule create_blockboundaries:
     input:
         "01_profile/sortedPoints.txt", "01_profile/psPoly.vtk", "01_profile/ssPoly.vtk",
-        "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"
+        "01_profile/geometry_paras.yml" ,"01_profile/midsPoly.vtk" ,"02_domainboundaries/inlet_2d.vtk",
+        "02_domainboundaries/outlet_2d.vtk", "02_domainboundaries/y_upper_2d.vtk", "02_domainboundaries/y_lower_2d.vtk"
 
     output:
-        "02_domain/ogridline_2d.vtk"
+        "03_meshgeometry/blocks.geom",
+        "03_meshgeometry/domain.geom",
+        "03_meshgeometry/blocklines.pdf"
+
     run:
         from ntrfc.preprocessing.mesh_creation.domain_creation import create_2d_blocklines
+
+        #todo: params einpflegen, files in funktion laden
         sortedPoints_fpath = input[0]
         psPoly_fpath = input[1]
         ssPoly_fpath = input[2]
@@ -75,6 +101,22 @@ rule create_blockboundaries:
 
         create_2d_blocklines(config["blade"],basedir,sortedPoints,psPoly,ssPoly,geometry_paras,midsPoly,inlet,outlet,y_upper,y_lower)
 
+"""
+rule run_igg_cascademesh:
+    input:
+        "02_domainboundaries/inlet_2d.vtk", "02_domainboundaries/outlet_2d.vtk",
+        "02_domainboundaries/y_upper_2d.vtk", "02_domainboundaries/y_lower_2d.vtk",
+        "03_blockboundaries/ogridline_2d.vtk","03_blockboundaries/ogrid_oulet.vtk",
+        "03_blockboundaries/ogrid_inlet.vtk","03_blockboundaries/te_ogrid.vtk","03_blockboundaries/le_ogrid.vtk",
+        "03_blockboundaries/yupper_ogrid_x1.vtk","03_blockboundaries/ylower_ogrid_x1.vtk","03_blockboundaries/yupper_ogrid_x0.vtk",
+        "03_blockboundaries/ylower_ogrid_x0.vtk","03_blockboundaries/ss_x1_ogrid_line.vtk","03_blockboundaries/ss_x0_ogrid_line.vtk",
+        "03_blockboundaries/ps_x1_ogrid_line.vtk", "03_blockboundaries/ps_x0_ogrid_line.vtk"
+    output:
+        "04_meshgeometry/blocks.geom",
+        "04_meshgeometry/domain.geom"
+    run:
+        from ntrf.preprocessing.mesh_creation.domain_creation import create_igg_geometry
+"""
 """
 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]
-    blade_shift : 0.015                                 #Shift Domain to have probes inside domain [m]
+    blade_shift : 0.005                                 #Shift Domain to have probes inside domain [m]

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
-from ntrfc.utils.math.vectorcalc import closest_node_index, vecAbs
+from ntrfc.utils.math.vectorcalc import closest_node_index, vecAbs, vecDir
+from ntrfc.utils.filehandling.tecplot import writeTecplot1DFile
 import numpy as np
 import pyvista as pv
 import os
@@ -44,7 +45,7 @@ def create_2d_domain(geosettings, basedir, midsPoly, ssPoly, psPoly, geometry_pa
 
     outletPoly = pv.Line(*outlet_pts)
 
-    domain_dir = os.path.join(basedir, "02_boundaries")
+    domain_dir = os.path.join(basedir, "02_domainboundaries")
     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)
@@ -83,7 +84,8 @@ def create_2d_blocklines(geo_settings, basedir, sortedPoints, psPoly, ssPoly, ge
     ogridhelpersurface.points += ogrid_size * ogridhelpersurface.point_normals
     ogridline = ogridhelpersurface.slice(normal="z")
     ogridline.points[:, 2] = 0
-    ogridline.save(os.path.join(basedir, "ogridline_2d.vtk"))
+    domain_dir = os.path.join(basedir,"03_meshgeometry")
+
 
     blademids = midsPoly.copy()
     # blademids.points += pitch/2
@@ -120,25 +122,79 @@ def create_2d_blocklines(geo_settings, basedir, sortedPoints, psPoly, ssPoly, ge
                                    geometry_paras["beta_meta_02"],
                                    inlet.bounds[0],outlet.bounds[0],0)
 
-    mspPoly = pv.PolyData(np.stack([msp_xx,msp_yy,np.zeros(len(msp_yy))]).T)
-    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(blademids)
-
-    p.add_mesh(mspPoly)
-    p.add_mesh(ogridline)
-    p.add_mesh(ps_x0_ogrid_line)
-    p.add_mesh(ps_x1_ogrid_line)
-    p.add_mesh(ss_x0_ogrid_line)
-    p.add_mesh(ss_x1_ogrid_line)
-    p.add_mesh(ylower_ogrid_x0)
-    p.add_mesh(yupper_ogrid_x0)
-    p.add_mesh(ylower_ogrid_x1)
-    p.add_mesh(yupper_ogrid_x1)
-    p.show(cpos=[0, 0, 1])
+    mspPoly = lines_from_points(np.stack([msp_xx,msp_yy,np.zeros(len(msp_yy))]).T)
+
+    le_ogrid = pv.Line(midsPoly.points[0],midsPoly.points[0]-vecDir(midsPoly.points[1]-midsPoly.points[0])*ogrid_size)
+    te_ogrid = pv.Line(midsPoly.points[-1],midsPoly.points[-1]-vecDir(midsPoly.points[-2]-midsPoly.points[-1])*ogrid_size)
+    ogrid_inlet_dist = vecAbs(mspPoly.points[0]-le_ogrid.points[-1])
+    ogrid_outlet_dist = vecAbs(mspPoly.points[-1]-te_ogrid.points[-1])
+    ogrid_inlet = pv.Line(le_ogrid.points[-1],le_ogrid.points[-1]-vecDir(midsPoly.points[1]-midsPoly.points[0])*ogrid_inlet_dist)
+    ogrid_oulet = pv.Line(te_ogrid.points[-1],te_ogrid.points[-1]-vecDir(midsPoly.points[-2]-midsPoly.points[-1])*ogrid_outlet_dist)
+    """
+    ogridline.save(os.path.join(domain_dir, "ogridline_2d.vtk"))
+
+
+    ogrid_oulet.save(os.path.join(domain_dir, "ogrid_oulet.vtk"))
+    ogrid_inlet.save(os.path.join(domain_dir, "ogrid_inlet.vtk"))
+    te_ogrid.save(os.path.join(domain_dir, "te_ogrid.vtk"))
+    le_ogrid.save(os.path.join(domain_dir, "le_ogrid.vtk"))
+    yupper_ogrid_x1.save(os.path.join(domain_dir, "yupper_ogrid_x1.vtk"))
+    ylower_ogrid_x1.save(os.path.join(domain_dir, "ylower_ogrid_x1.vtk"))
+    yupper_ogrid_x0.save(os.path.join(domain_dir, "yupper_ogrid_x0.vtk"))
+    ylower_ogrid_x0.save(os.path.join(domain_dir, "ylower_ogrid_x0.vtk"))
+    ss_x1_ogrid_line.save(os.path.join(domain_dir, "ss_x1_ogrid_line.vtk"))
+    ss_x0_ogrid_line.save(os.path.join(domain_dir, "ss_x0_ogrid_line.vtk"))
+    ps_x1_ogrid_line.save(os.path.join(domain_dir, "ps_x1_ogrid_line.vtk"))
+    ps_x0_ogrid_line.save(os.path.join(domain_dir, "ps_x0_ogrid_line.vtk"))
+    """
+    plt.figure()
+    plt.plot(inlet.points[::, 0], inlet.points[::, 1], color="#6c3376")
+    plt.plot(outlet.points[::, 0], outlet.points[::, 1], color="#FF2211")
+    plt.plot(y_upper.points[::, 0], y_upper.points[::, 1], color="#FF22CC")
+    plt.plot(y_lower.points[::, 0], y_lower.points[::, 1], color="#FF22CC")
+    plt.plot(bladeline.points[::, 0], bladeline.points[::, 1], color="#6c3376")
+    plt.plot(le_ogrid.points[::, 0], le_ogrid.points[::, 1], color="#FF2211")
+    plt.plot(te_ogrid.points[::, 0], te_ogrid.points[::, 1], color="#FF22CC")
+    plt.plot(ogrid_inlet.points[::, 0], ogrid_inlet.points[::, 1], color="#FF22CC")
+    plt.plot(ogrid_oulet.points[::, 0], ogrid_oulet.points[::, 1], color="#6c3376")
+    plt.plot(ogridline.points[::, 0], ogridline.points[::, 1], color="#FF2211")
+    plt.plot(ps_x0_ogrid_line.points[::, 0], ps_x0_ogrid_line.points[::, 1], color="#FF22CC")
+    plt.plot(ps_x1_ogrid_line.points[::, 0], ps_x1_ogrid_line.points[::, 1], color="#FF22CC")
+    plt.plot(ss_x0_ogrid_line.points[::, 0], ss_x0_ogrid_line.points[::, 1], color="#FF22CC")
+    plt.plot(ss_x1_ogrid_line.points[::, 0], ss_x1_ogrid_line.points[::, 1], color="#FF22CC")
+    plt.plot(ylower_ogrid_x0.points[::, 0], ylower_ogrid_x0.points[::, 1], color="#FF22CC")
+    plt.plot(yupper_ogrid_x0.points[::, 0], yupper_ogrid_x0.points[::, 1], color="#FF22CC")
+    plt.plot(ylower_ogrid_x1.points[::, 0], ylower_ogrid_x1.points[::, 1], color="#FF22CC")
+    plt.plot(yupper_ogrid_x1.points[::, 0], yupper_ogrid_x1.points[::, 1], color="#FF22CC")
+    plt.xlabel('x')
+    plt.ylabel('y')
+    plt.savefig(os.path.join(domain_dir, 'blocklines.pdf'))
+
+
+    writeTecplot1DFile(os.path.join(domain_dir,"domain.geom"),['x', 'z'],
+                       ["inlet", "outlet", "y_upper", "y_lower", "psPoly", "ssPoly"],
+                       [[inlet.points[::,0],inlet.points[::,1]],[outlet.points[::,0],outlet.points[::,1]],
+                        [y_upper.points[::,0],y_upper.points[::,1]],[y_lower.points[::,0],y_lower.points[::,1]],
+                        [psPoly.points[::,0],psPoly.points[::,1]],[ssPoly.points[::,0],ssPoly.points[::,1]]],
+                       "domainboundaries")
+
+    writeTecplot1DFile(os.path.join(domain_dir,"blocks.geom"),['x', 'z'],
+                       ["le_ogrid","te_ogrid","ogrid_inlet","ogrid_oulet","ogridline",
+                       "ps_x0_ogrid_line","ps_x1_ogrid_line","ss_x0_ogrid_line","ss_x1_ogrid_line",
+                       "ylower_ogrid_x0","yupper_ogrid_x0","ylower_ogrid_x1","yupper_ogrid_x1"],
+                       [[le_ogrid.points[::,0],le_ogrid.points[::,1]],
+                        [te_ogrid.points[::,0],te_ogrid.points[::,1]],
+                        [ogrid_inlet.points[::,0],ogrid_inlet.points[::,1]],
+                        [ogrid_oulet.points[::,0],ogrid_oulet.points[::,1]],
+                        [ogridline.points[::,0],ogridline.points[::,1]],
+                        [ps_x0_ogrid_line.points[::,0],ps_x0_ogrid_line.points[::,1]],
+                        [ps_x1_ogrid_line.points[::,0],ps_x1_ogrid_line.points[::,1]],
+                        [ss_x0_ogrid_line.points[::,0],ss_x0_ogrid_line.points[::,1]],
+                        [ss_x1_ogrid_line.points[::,0],ss_x1_ogrid_line.points[::,1]],
+                        [ylower_ogrid_x0.points[::,0],ylower_ogrid_x0.points[::,1]],
+                        [yupper_ogrid_x0.points[::,0],yupper_ogrid_x0.points[::,1]],
+                        [ylower_ogrid_x1.points[::,0],ylower_ogrid_x1.points[::,1]]],
+                       "blockboundaries")
+
+
     return 0

ntrfc/utils/filehandling/tecplot.py

+import os
+
+
+def writeTecplot1DFile(output_path, var_string, zone_string, values, title):
+    # var_string: namen der variablen als liste ['U','p']
+    # zone_string: namen der zonen als liste ['saugseute','druckseite']
+    # values: erster index der liste steht fuer zone, dann folgen die listen der eigentlichen variablen
+    # Beispiel: [[[10,11,10],[10000,11000,12000]],[[10,11,10],[10000,11000,12000]]]
+    data = open(os.path.join(output_path), 'w')
+    data.write('TITLE     ="' + title + '"\n')
+    var = 'VARIABLES = '
+    for i in range(len(var_string)):
+        if i < len(var_string) - 1:
+            var = var + '"' + var_string[i] + '", '
+        else:
+            var = var + '"' + var_string[i] + '"\n'
+
+    data.write(var)
+
+    for i in range(len(values)):
+        data.write('ZONE T="' + zone_string[i] + '",I=' + str(len(values[i][0])) + '\n')
+
+        # data.write('ZONE  T= "'+zone_string[i]+'"\n')
+        for j in range(len(values[i][0])):
+            line_string = ''
+            for k in range(len(values[i])):
+                if k < len(values[i]) - 1:
+                    line_string = line_string + str(values[i][k][j]) + '\t'
+                else:
+                    line_string = line_string + str(values[i][k][j]) + '\n'
+            data.write(line_string)
+
+    data.close()
+
+
+def writeTecplot2D3DFile(filename, X, Y, Z, vars):
+    def pad(s, width):
+        s2 = s
+        while len(s2) < width:
+            s2 = ' ' + s2
+        if s2[0] != ' ':
+            s2 = ' ' + s2
+        if len(s2) > width:
+            s2 = s2[:width]
+        return s2
+
+    def varline(vars, id, fw):
+        s = ""
+        for v in vars:
+            s = s + pad(str(v[1][id]), fw)
+        s = s + '\n'
+        return s
+
+    fw = 10  # field width
+
+    f = open(filename, "wt")
+
+    f.write('Variables="X","Y"')
+    if len(Z) > 0:
+        f.write(',"Z"')
+    for v in vars:
+        f.write(',"%s"' % v[0])
+    f.write('\n\n')
+
+    f.write('Zone I=' + pad(str(len(X)), 6) + ',J=' + pad(str(len(Y)), 6))
+    if len(Z) > 0:
+        f.write(',K=' + pad(str(len(Z)), 6))
+    f.write(', F=POINT\n')
+
+    if len(Z) > 0:
+        id = 0
+        for k in range(len(Z)):
+            for j in range(len(Y)):
+                for i in range(len(X)):
+                    f.write(pad(str(X[i]), fw) + pad(str(Y[j]), fw) + pad(str(Z[k]), fw))
+                    f.write(varline(vars, id, fw))
+                    id = id + 1
+    else:
+        id = 0
+        for j in range(len(Y)):
+            for i in range(len(X)):
+                f.write(pad(str(X[i]), fw) + pad(str(Y[j]), fw))
+                f.write(varline(vars, id, fw))
+                id = id + 1
+
+    f.close()
+
+
+def openTecplotFile(path):
+
+    values = []
+    var = []
+    zones = []
+
+    zone_bool = -1
+    with open(path, 'r') as f:
+        for line in f:
+            if line.startswith('VARIABLES'):
+                line_string = line.replace('\n', '').split('=')[-1].split(',')
+                for string in line_string:
+                    var.append(string.replace('"', ''))
+
+            if line.startswith('ZONE'):
+
+                zone_bool = zone_bool + 1
+
+                zones.append(line.replace('\n', '').split('=')[-1].replace('"', ''))
+                list_01 = []
+                for i in range(len(var)):
+                    list_01.append([])
+
+                values.append(list_01)
+
+            if line.startswith('ZONE') or line.startswith('VARIABLES') or line.startswith('TITLE'):
+                pass
+
+            else:
+                line_values = line.split()
+                for i in range(len(line_values)):
+                    values[zone_bool][i].append(float(line_values[i]))
+
+    return values
+
+
+def getTecplotFileVarNames(path):
+    var = []
+
+    with open(path, 'r') as f:
+        for line in f:
+            if line.startswith('VARIABLES'):
+                line = line.replace(' ', '')
+                line_string = line.replace('\n', '').split('=')[-1].split(',')
+
+                for string in line_string:
+                    var.append(string.replace('"', ''))
+
+    return var
+
+
+def getTecplotFileZoneNames(path):
+    var = []
+
+    with open(path, 'r') as f:
+        for line in f:
+            if line.startswith('ZONE') or line.startswith('Zone') or line.startswith('zone'):
+                line = line.replace(' ', '')
+                line_string = line.replace('\n', '').split('"')[1]
+
+                var.append(line_string)
+
+    return var

tests/test_ntrfc.py

@@ -21,7 +21,7 @@ def test_yamlDictRead(tmpdir):
     """
     tests if yaml is returning a known dictionary
     """
-    from ntrfc.utils.filehandling.read_datafiles import yaml_dict_read
+    from ntrfc.utils.filehandling.datafiles import yaml_dict_read
 
     test_file = tmpdir / "test.yaml"
     with open(test_file, "w") as handle:
@@ -33,7 +33,7 @@ def test_yamlDictWrite(tmpdir):
     """
     tests if yaml is writing and returning a known dictionary
     """
-    from ntrfc.utils.filehandling.read_datafiles import yaml_dict_read, write_yaml_dict
+    from ntrfc.utils.filehandling.datafiles import yaml_dict_read, write_yaml_dict
 
     test_file = tmpdir / "test.yaml"
     test_dict = {"test_key": True}
@@ -62,7 +62,7 @@ def test_loadmesh_vtk(tmpdir):
     """
     tests if a vtk mesh can be read and Density is translated to rho
     """
-    from ntrfc.utils.filehandling.read_mesh import load_mesh
+    from ntrfc.utils.filehandling.mesh import load_mesh
 
     test_file = tmpdir / "tmp.vtk"
     mesh = pv.Box()
@@ -84,28 +84,28 @@ def test_surface_distance():
 
 
 def test_cgnsReader():
-    from ntrfc.utils.filehandling.read_mesh import cgnsReader
+    from ntrfc.utils.filehandling.mesh import cgnsReader
     # todo fill
     a = cgnsReader
     return 0
 
 
 def test_vtkUnstructuredGridReader():
-    from ntrfc.utils.filehandling.read_mesh import vtkUnstructuredGridReader
+    from ntrfc.utils.filehandling.mesh import vtkUnstructuredGridReader
     # todo fill
     a = vtkUnstructuredGridReader
     return 0
 
 
 def test_vtkFLUENTReader():
-    from ntrfc.utils.filehandling.read_mesh import vtkFLUENTReader
+    from ntrfc.utils.filehandling.mesh import vtkFLUENTReader
     # todo fill
     a = vtkFLUENTReader
     return 0
 
 
 def test_pickle_operations(tmpdir):
-    from ntrfc.utils.filehandling.read_datafiles import write_pickle, read_pickle, write_pickle_protocolzero
+    from ntrfc.utils.filehandling.datafiles import write_pickle, read_pickle, write_pickle_protocolzero
 
     fname = tmpdir / "test.pkl"
     dict = {"test": 1}