diff --git a/samples/__init__.py b/samples/__init__.py
index b752a17..e165fbf 100644
--- a/samples/__init__.py
+++ b/samples/__init__.py
@@ -24,24 +24,32 @@ def genMesh(stype, theta, flowdirection, saveto):
     if structure:
         salome.salome_init()
 
-        grains, cubic, rhombohedron = structure(theta)
-        fd = namedtuple("fd", ["x", "y", "z"])
-        geometry = cubic
+        #grains, cubic, rhombohedron = structure(theta)
+        #fd = namedtuple("fd", ["x", "y", "z"])
+        grains, geometry1, geometry2 = structure(theta)
+        geometry = geometry1
         
         if flowdirection == [1, 1, 1]:
-            geometry = rhombohedron
-            direction = fd([1, 1, 1], [1, -1, 1], [1, -1, -1])
+            geometry = geometry2
+            norm = [-1, -1, 1]
+            bcount = 6
+            #direction = fd([1, 1, 1], [1, -1, 1], [1, -1, -1])
         
-        else:
-            geometry = cubic
+        #else:
+        #    geometry = cubic
 
         if flowdirection == [1, 0, 0]:
-            direction = fd([1, 1, 0], [1, -1, 0], [0, 0, 1])
+            norm = [0, 0, 1]
+            bcount = 4
+            #direction = fd([1, 1, 0], [1, -1, 0], [0, 0, 1])
 
         if flowdirection == [0, 0, 1]:
-            direction = fd([0, 0, 1], [1, -1, 0], [1, 1, 0])
+            norm = [1, 1, 0]
+            bcount = 4
+            #direction = fd([0, 0, 1], [1, -1, 0], [1, 1, 0])
         
-        boundary = geometry_utils.boundaryCreate(geometry, direction, grains)
+        #boundary = geometry_utils.boundaryCreate(geometry, direction, grains)
+        boundary = geometry_utils.createBoundary(geometry, bcount, flowdirection, norm, grains)
 
         fineness = 3
         mesh = mesh_utils.meshCreate(geometry, boundary, fineness)
diff --git a/samples/bodyCenteredCubic.py b/samples/bodyCenteredCubic.py
index c36d8fe..5cfb7aa 100644
--- a/samples/bodyCenteredCubic.py
+++ b/samples/bodyCenteredCubic.py
@@ -149,7 +149,19 @@ def bodyCenteredCubic(alpha):
     # Preparation
     grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
     grains += geompy.ExtractShapes(Middle, geompy.ShapeType["SOLID"], True)
-
     grains = geompy.MakeFuseList(grains, False, False)
 
+    R = r0 / (1 - alpha)
+    C1 = 0.8
+    C2 = 0.4
+    alpha1 = 0.01
+    alpha2 = 0.13
+    
+    Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
+    R_fillet = Cf * (r0 * math.sqrt(2) - R)
+    
+    grains = geompy.MakeFilletAll(grains, R_fillet)
+    geometry1 = geompy.MakeCutList(Pore_1a, [grains], True)
+    geometry2 = geompy.MakeCutList(Pore_3c, [grains], True)
+
     return grains, Pore_1a, Pore_3c
diff --git a/samples/faceCenteredCubic.py b/samples/faceCenteredCubic.py
index fe16916..17e40a9 100644
--- a/samples/faceCenteredCubic.py
+++ b/samples/faceCenteredCubic.py
@@ -144,7 +144,19 @@ def faceCenteredCubic(alpha):
     # Preparation
     grains = geompy.ExtractShapes(Up_Down, geompy.ShapeType["SOLID"], True)
     grains += geompy.ExtractShapes(Middle_Up, geompy.ShapeType["SOLID"], True)
-
     grains = geompy.MakeFuseList(grains, False, False)
+    
+    R = r0 / (1 - alpha)
+    C1 = 0.8
+    C2 = 0.4
+    alpha1 = 0.01
+    alpha2 = 0.18
+    
+    Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
+    R_fillet = Cf * (r0 * math.sqrt(2) - R)
+    
+    grains = geompy.MakeFilletAll(grains, R_fillet)
+    geometry1 = geompy.MakeCutList(Common, [grains], True)
+    geometry2 = geompy.MakeCutList(Cut_8, [grains], True)
 
-    return grains, Common, Cut_8
+    return grains, geometry1, geometry2
diff --git a/samples/simpleCubic.py b/samples/simpleCubic.py
index dc71c69..bb621a1 100644
--- a/samples/simpleCubic.py
+++ b/samples/simpleCubic.py
@@ -113,7 +113,19 @@ def simpleCubic(alpha):
 
     # Preparation
     grains = geompy.ExtractShapes(Multi_Translation_3, geompy.ShapeType["SOLID"], True)
-
     grains = geompy.MakeFuseList(grains, False, False)
 
-    return grains, Cut_1, Cut_4
+    R = r0 / (1 - alpha)
+    C1 = 0.8
+    C2 = 0.4
+    alpha1 = 0.01
+    alpha2 = 0.28
+    
+    Cf = C1 + (C2 - C1) / (alpha2 - alpha1) * (alpha - alpha1)
+    R_fillet = Cf * (r0 * math.sqrt(2) - R)
+    
+    grains = geompy.MakeFilletAll(grains, R_fillet)
+    geometry1 = geompy.MakeCutList(Cut_1, [grains], True)
+    geometry2 = geompy.MakeCutList(Cut_4, [grains], True)
+
+    return grains, geometry1, geometry2
diff --git a/src/geometry_utils.py b/src/geometry_utils.py
index f921a61..e520631 100644
--- a/src/geometry_utils.py
+++ b/src/geometry_utils.py
@@ -41,11 +41,15 @@ def createGroup(gobj, planelist, grains, name):
     return gr
 
 
-def createBoundary(gobj, bcount, dvec, grains):
+def createBoundary(gobj, bcount, dvec, norm, grains):
+    
+    direction = np.quaternion(0, dvec[0], dvec[1], dvec[2]).normalized()
+    ax = lambda alpha: np.quaternion(
+        np.cos(alpha * 0.5),
+        np.sin(alpha * 0.5) * norm[0], 
+        np.sin(alpha * 0.5) * norm[1], 
+        np.sin(alpha * 0.5) * norm[2]) 
     
-    direction = np.quaternion(0, dvec.x[0], dvec.x[1], dvec.x[2]).normalized()
-    ax = lambda alpha: np.quaternion(np.cos(alpha * 0.5), \
-        np.sin(alpha * 0.5) * 1, np.sin(alpha * 0.5) * 1, np.sin(alpha * 0.5) * 1) 
     ang = lambda n, count: 2 * np.pi * n / count
     limit = bcount if np.mod(bcount, 2) else int(bcount / 2)
 
@@ -54,7 +58,7 @@ def createBoundary(gobj, bcount, dvec, grains):
     #
     flowvec = geompy.MakeVector(
         geompy.MakeVertex(0, 0, 0),
-        geompy.MakeVertex(dvec.x[0], dvec.x[1], dvec.x[2]))
+        geompy.MakeVertex(dvec[0], dvec[1], dvec[2]))
     symvec = []
 
     for qvec in vecs: