Quaternions: fix

run.py

@@ -33,7 +33,7 @@ def createTasks():
 
     for structure in structures:
         if structure == "simpleCubic":
-            theta = [c * 0.01 for c in range(1, 14)]
+            theta = [c * 0.01 for c in range(1, 29 + 1)]
 
         elif structure == "faceCenteredCubic":
             theta = []

samples/__init__.py

@@ -37,31 +37,35 @@ def genMesh(stype, theta, flowdirection, saveto):
 
             # initial angle
             angle = math.pi / 6
-            vec = Quaternion(0, flowdirection[0], flowdirection[1], flowdirection[2])
+            #vec = Quaternion(0, norm[0], norm[1], norm[2])
-            ax = Quaternion(
+            #ax = Quaternion(
-                math.cos(angle * 0.5), 
+            #    math.cos(angle * 0.5), 
-                math.sin(angle * 0.5) * norm[0], 
+            #    math.sin(angle * 0.5) * flowdirection[0], 
-                math.sin(angle * 0.5) * norm[1], 
+            #    math.sin(angle * 0.5) * flowdirection[1], 
-                math.sin(angle * 0.5) * norm[2])
+            #    math.sin(angle * 0.5) * flowdirection[2])
-            qvec = (ax * vec * ax.inverse).vector
+            #qvec = (ax * vec * ax.inverse).vector
-            norm = [qvec[0], qvec[1], qvec[2]]
+            #normvec = [qvec[0], qvec[1], qvec[2]]
+            normvec = Quaternion(axis = flowdirection, angle = angle).rotate(norm)
+            direction = [1, 1, 1]
             #direction = fd([1, 1, 1], [1, -1, 1], [1, -1, -1])
         
         #else:
         #    geometry = cubic
 
         if flowdirection == [1, 0, 0]:
-            norm = [0, 0, 1]
+            normvec = [0, 0, 1]
             bcount = 4
+            direction = [1, 1, 0]
             #direction = fd([1, 1, 0], [1, -1, 0], [0, 0, 1])
 
         if flowdirection == [0, 0, 1]:
-            norm = [1, 1, 0]
+            normvec = [1, 1, 0]
             bcount = 4
+            direction = [0, 0, 1]
             #direction = fd([0, 0, 1], [1, -1, 0], [1, 1, 0])
         
         #boundary = geometry_utils.boundaryCreate(geometry, direction, grains)
-        boundary = geometry_utils.createBoundary(geometry, bcount, flowdirection, norm, grains)
+        boundary = geometry_utils.createBoundary(geometry, bcount, direction, normvec, grains)
 
         fineness = 3
         mesh = mesh_utils.meshCreate(geometry, boundary, fineness)

samples/bodyCenteredCubic.py

@@ -127,22 +127,21 @@ def bodyCenteredCubic(alpha):
     Pore_3b = geompy.MakeCutList(Cut_3b, [Middle])
     Pore_3c = geompy.MakeCutList(Cut_3c, [Middle])
 
-    geompy.addToStudy( Sphere_1, 'Sphere_'+str(i+1) )
+    geompy.addToStudy( Sphere_1, 'Sphere_' )
-    geompy.addToStudy( Down, 'Down_'+str(i+1) )
+    geompy.addToStudy( Down, 'Down_' )
-    geompy.addToStudy( Up, 'Up_'+str(i+1) )
+    geompy.addToStudy( Up, 'Up_' )
-    geompy.addToStudy( Up_Down, 'Up_Down_'+str(i+1) )
+    geompy.addToStudy( Up_Down, 'Up_Down_' )
-    geompy.addToStudy( Cut_1a, 'Cut_1a_'+str(i+1) )
+    geompy.addToStudy( Cut_1a, 'Cut_1a_' )
-    geompy.addToStudy( Cut_1b, 'Cut_1b_'+str(i+1) )
+    geompy.addToStudy( Cut_1b, 'Cut_1b_' )
-    geompy.addToStudy( Middle, 'Middle_'+str(i+1) )
+    geompy.addToStudy( Middle, 'Middle_' )
-    geompy.addToStudy( Pore_1a, 'Pore_1a_'+str(i+1) )
+    geompy.addToStudy( Pore_1a, 'Pore_1a_' )
-    geompy.addToStudy( Pore_1b, 'Pore_1b_'+str(i+1) )
+    geompy.addToStudy( Pore_1b, 'Pore_1b_' )
-    geompy.addToStudy( Pore_2a, 'Pore_2a_'+str(i+1) )
+    geompy.addToStudy( Pore_2a, 'Pore_2a_' )
-    geompy.addToStudy( Pore_2b, 'Pore_2b_'+str(i+1) )
+    geompy.addToStudy( Pore_2b, 'Pore_2b_' )
-    geompy.addToStudy( Pore_3a, 'Pore_3a_'+str(i+1) )
+    geompy.addToStudy( Pore_3a, 'Pore_3a_' )
-    geompy.addToStudy( Pore_3b, 'Pore_3b_'+str(i+1) )
+    geompy.addToStudy( Pore_3b, 'Pore_3b_' )
-    geompy.addToStudy( Pore_3c, 'Pore_3c_'+str(i+1) )
+    geompy.addToStudy( Pore_3c, 'Pore_3c_' )
 
-    i = i + 1
     if salome.sg.hasDesktop():
         salome.sg.updateObjBrowser()
 

samples/faceCenteredCubic.py

@@ -125,19 +125,18 @@ def faceCenteredCubic(alpha):
     Cut_7 = geompy.MakeCutList(Extrusion_4, [Up_Down], True)
     Cut_8 = geompy.MakeCutList(Cut_7, [Middle_Up], True)
 
-    #geompy.addToStudy( Sphere_1, 'Sphere_'+str(i+1) )
+    #geompy.addToStudy( Sphere_1, 'Sphere_' )
-    geompy.addToStudy( Down, 'Down_'+str(i+1) )
+    geompy.addToStudy( Down, 'Down_' )
-    geompy.addToStudy( Up_Down, 'Up_Down_'+str(i+1) )
+    geompy.addToStudy( Up_Down, 'Up_Down_' )
-    geompy.addToStudy( Cut_1, 'Cut_1_'+str(i+1) )
+    geompy.addToStudy( Cut_1, 'Cut_1_' )
-    geompy.addToStudy( Middle_Up, 'Middle_Up_'+str(i+1) )
+    geompy.addToStudy( Middle_Up, 'Middle_Up_' )
-    geompy.addToStudy( Cut_2, 'Pore_'+str(i+1) )
+    geompy.addToStudy( Cut_2, 'Pore_' )
-    geompy.addToStudy( Common, 'Pore_2_'+str(i+1) )
+    geompy.addToStudy( Common, 'Pore_2_' )
-    geompy.addToStudy( Pore_3, 'Pore_3_'+str(i+1) )
+    geompy.addToStudy( Pore_3, 'Pore_3_' )
-    geompy.addToStudy( Cut_4, 'Pore_3a_'+str(i+1) )
+    geompy.addToStudy( Cut_4, 'Pore_3a_' )
-    geompy.addToStudy( Cut_6, 'Pore_3b_'+str(i+1) )
+    geompy.addToStudy( Cut_6, 'Pore_3b_' )
-    geompy.addToStudy( Cut_8, 'Pore_3c_'+str(i+1) )
+    geompy.addToStudy( Cut_8, 'Pore_3c_' )
 
-    i = i + 1
     if salome.sg.hasDesktop():
         salome.sg.updateObjBrowser()
 

samples/simpleCubic.py

@@ -99,15 +99,14 @@ def simpleCubic(alpha):
     Cut_4 = geompy.MakeCutList(Extrusion_3, [Multi_Translation_3])
 
     geompy.addToStudy( Sphere_1, 'Sphere_' )
-    geompy.addToStudy( Multi_Translation_2, 'Multi-Translation_2_'+str(i+1)  )
+    geompy.addToStudy( Multi_Translation_2, 'Multi-Translation_2_'  )
-    geompy.addToStudy( Multi_Translation_3, 'Multi-Translation_3_'+str(i+1)  )
+    geompy.addToStudy( Multi_Translation_3, 'Multi-Translation_3_'  )
-    geompy.addToStudy( Cut_1, 'Pore1_'+str(i+1) )
+    geompy.addToStudy( Cut_1, 'Pore1_' )
-    geompy.addToStudy( Cut_2, 'Pore2_'+str(i+1) )
+    geompy.addToStudy( Cut_2, 'Pore2_' )
-    geompy.addToStudy( Cut_3, 'Pore3_'+str(i+1) )
+    geompy.addToStudy( Cut_3, 'Pore3_' )
-    geompy.addToStudy( Cut_V, 'Cut_V_'+str(i+1) )
+    geompy.addToStudy( Cut_V, 'Cut_V_' )
-    geompy.addToStudy( Cut_4, 'Pore4_'+str(i+1) )
+    geompy.addToStudy( Cut_4, 'Pore4_' )
 
-    i = i + 1
     if salome.sg.hasDesktop():
         salome.sg.updateObjBrowser()
 

src/foam_utils.py

@@ -6,10 +6,6 @@ from datetime import timedelta
 
 def application(name, case, log=False, args=[], parallel=False):
     
-
-    #if log:
-    #    logfile = open("{}/{}.log".format(case, name), "a")
-    
     mpirun = []
     if parallel:
         mpirun = ["mpirun", "-np", "4", "--oversubscribe"]
@@ -18,7 +14,6 @@ def application(name, case, log=False, args=[], parallel=False):
     logging.info("Running '{}'".format(" ".join(cmd)))
    
     with subprocess.Popen(cmd, 
-        #shell = True,
         stdout = subprocess.PIPE, 
         stderr = subprocess.PIPE) as p, \
         open("{}/{}.log".format(case, name), "wb") as logfile:

src/geometry_utils.py

@@ -1,6 +1,3 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
 import GEOM
 from salome.geom import geomBuilder
 geompy = geomBuilder.New()
@@ -14,19 +11,19 @@ import numpy as np
 def getGeom():
     return geompy
 
-def rotate(gobj, ang):
+#def rotate(gobj, ang):
-    x = geompy.MakeVectorDXDYDZ(1, 0, 0)
+#    x = geompy.MakeVectorDXDYDZ(1, 0, 0)
-    y = geompy.MakeVectorDXDYDZ(0, 1, 0)
+#    y = geompy.MakeVectorDXDYDZ(0, 1, 0)
-    z = geompy.MakeVectorDXDYDZ(0, 0, 1)
+#    z = geompy.MakeVectorDXDYDZ(0, 0, 1)
 
     # yaw
-    rotated = geompy.MakeRotation(gobj, z, ang[2])
+#    rotated = geompy.MakeRotation(gobj, z, ang[2])
     # pitch
-    rotated = geompy.MakeRotation(rotated, y, ang[1])
+#    rotated = geompy.MakeRotation(rotated, y, ang[1])
     # roll
-    rotated = geompy.MakeRotation(rotated, x, ang[0])
+#    rotated = geompy.MakeRotation(rotated, x, ang[0])
 
-    return rotated
+#    return rotated
 
 def createGroup(gobj, planelist, grains, name):
     gr = geompy.CreateGroup(gobj, geompy.ShapeType["FACE"], name)
@@ -43,17 +40,26 @@ def createGroup(gobj, planelist, grains, name):
 
 def createBoundary(gobj, bcount, dvec, norm, grains):
     
-    direction = Quaternion(0, dvec[0], dvec[1], dvec[2]).normalized
+    #normvec = Quaternion(0, norm[0], norm[1], norm[2]).normalised
-    ax = lambda alpha: Quaternion(
+    #ax = lambda alpha: Quaternion(
-        np.cos(alpha * 0.5),
+    #    np.cos(alpha * 0.5),
-        np.sin(alpha * 0.5) * norm[0], 
+    #    np.sin(alpha * 0.5) * dvec[0], 
-        np.sin(alpha * 0.5) * norm[1], 
+    #    0, 0).normalised
-        np.sin(alpha * 0.5) * norm[2]) 
+        #np.sin(alpha * 0.5) * dvec[1], 
+        #np.sin(alpha * 0.5) * dvec[2]).normalised 
     
     ang = lambda n, count: 2 * np.pi * n / count
     limit = bcount if np.mod(bcount, 2) else int(bcount / 2)
 
-    vecs = [ ax(ang(n, bcount)) * direction * ax(ang(n, bcount)).inverse for n in range(limit) ]
+    #vecs = [ ax(ang(n, bcount)) * normvec * ax(ang(n, bcount)).inverse for n in range(limit) ]
+    vecs = [ Quaternion(axis = dvec, angle = ang(n, bcount)).rotate(norm) for n in range(limit) ]
+    
+    logging.info("""createBoundary: 
+    flow direction:\t{}
+    side boundaries:\t{}
+    normal direction:\t{}
+    angles:\t{}
+    side directions:\t{}""".format(dvec, bcount, norm, [ ang(n, bcount) / (2 * np.pi) * 360 for n in range(limit) ], [ v for v in vecs ]))
 
     #
     flowvec = geompy.MakeVector(
@@ -61,8 +67,8 @@ def createBoundary(gobj, bcount, dvec, norm, grains):
         geompy.MakeVertex(dvec[0], dvec[1], dvec[2]))
     symvec = []
 
-    for qvec in vecs:
+    for vec in vecs:
-        vec = qvec.vector
+        #vec = qvec.vector
         symvec.append(geompy.MakeVector(
             geompy.MakeVertex(0, 0, 0),
             geompy.MakeVertex(vec[0], vec[1], vec[2])))
@@ -73,6 +79,7 @@ def createBoundary(gobj, bcount, dvec, norm, grains):
     planes = geompy.MakeCompound(planes)
     planes = geompy.MakeCutList(planes, [grains], False)
     planes = geompy.ExtractShapes(planes, geompy.ShapeType["FACE"], False)
+    #print("planes: ", len(planes))
 
     inletplanes = []
     outletplanes = []
@@ -82,6 +89,7 @@ def createBoundary(gobj, bcount, dvec, norm, grains):
         nvec = geompy.GetNormal(plane)
         
         fwang = round(geompy.GetAngle(nvec, flowvec), 0)
+        #print("fwang = ", fwang)
 
         if fwang == 0:
             inletplanes.append(plane)
@@ -90,13 +98,26 @@ def createBoundary(gobj, bcount, dvec, norm, grains):
             outletplanes.append(plane)
 
         for n in range(len(symvec)):
-            sang = round(geompy.GetAngle(nvec, svec[n]), 0)
+            sang = round(geompy.GetAngle(nvec, symvec[n]), 0)
+            #print("sang = ", sang, "\n")
 
             if sang == 0:
-                symetryplanes[n][0] = plane
+                if symetryplanes[n][0] == None:
+                    symetryplanes[n][0] = []
+
+                symetryplanes[n][0].append(plane)
 
             elif sang == 180:
-                symetryplanes[n][1] = plane
+                if symetryplanes[n][1] == None:
+                    symetryplanes[n][1] = []
+
+                symetryplanes[n][1].append(plane)
+
+    logging.info("""createBoundary: 
+    planes:\t{}
+    inlet planes:\t{}
+    outlet planes:\t{}
+    side planes:\t{}""".format(len(planes), inletplanes, outletplanes, symetryplanes))
 
     #
     boundary = {}

src/mesh_utils.py

@@ -1,6 +1,3 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
 import SMESH
 from salome.smesh import smeshBuilder
 smesh = smeshBuilder.New()
@@ -10,6 +7,7 @@ import logging
 def getSmesh():
     return smesh
 
+
 def meshCreate(gobj, boundary, fineness, viscousLayers=None):
     """
     Creates a mesh from a geometry.
@@ -89,6 +87,7 @@ def meshCreate(gobj, boundary, fineness, viscousLayers=None):
 
     return mesh
 
+
 def meshCompute(mobj):
     """Compute the mesh."""
     status = mobj.Compute()
@@ -103,6 +102,7 @@ def meshCompute(mobj):
     logging.info("""meshCompute:
     status:\t{}""".format(msg))
 
+
 def meshExport(mobj, path):
     """
     Export the mesh in a file in UNV format.
@@ -110,10 +110,9 @@ def meshExport(mobj, path):
     Parameters:
         path (string): full path to the expected directory.
     """
-    exportpath = path #os.path.join(path, "{}.unv".format(mobj.name))
 
     try:
-        mobj.ExportUNV(exportpath)
+        mobj.ExportUNV(path)
 
         logging.info("""meshExport:
         format:\t{}""".format("unv"))