Ms homo

setup.py

@@ -41,7 +41,7 @@
     keywords=[
         # eg: 'keyword1', 'keyword2', 'keyword3',
     ],
-    
+
     packages=['bgem', 'bgem.polygons', 'bgem.bspline', 'bgem.gmsh', 'bgem.external', 'bgem.geometry', 'bgem.stochastic',
               'bgem.core', 'bgem.upscale'], #setuptools.find_packages(where='src'),
     package_dir={'': 'src'},

src/bgem/core/flow_call.py

@@ -57,12 +57,12 @@ def __init__(self, process: subprocess.CompletedProcess, stdout: File, stderr: F
         self.process = process
         self.stdout = stdout
         self.stderr = stderr
-        with workdir(output_dir):
-            self.log = File("flow123.0.log")
-            # TODO: flow ver 4.0 unify output file names
-            self.hydro = EquationOutput("flow", "water")
-            self.solute = EquationOutput("solute", "mass")
-            self.mechanic = EquationOutput("mechanics", "mechanics")
+        #with workdir(output_dir):
+        self.log = File("flow123.0.log")
+        # TODO: flow ver 4.0 unify output file names
+        self.hydro = EquationOutput("flow", "water")
+        self.solute = EquationOutput("solute", "mass")
+        self.mechanic = EquationOutput("mechanics", "mechanics")
 
     @property
     def success(self):
@@ -88,10 +88,20 @@ def check_conv_reasons(self):
                     continue
         return True
 
+
 #@memoize
 def _prepare_inputs(file_in, params):
+    import pathlib
     in_dir, template = os.path.split(file_in)
-    root = template.removesuffix(".yaml").removesuffix("_tmpl")
+
+    root = pathlib.Path(template).with_suffix(".yaml")
+    try:
+        root = pathlib.Path(root).with_suffix("_tmpl")
+    except:
+        pass
+
+    root = str(root)
+    #root = template.removesuffix(".yaml").removesuffix("_tmpl")
     template_path = Path(file_in).rename(Path(in_dir) / (root + "_tmpl.yaml"))
     #suffix = "_tmpl.yaml"
     #assert template[-len(suffix):] == suffix
@@ -100,6 +110,7 @@ def _prepare_inputs(file_in, params):
     main_input, used_params =  substitute_placeholders(str(template_path), str(main_input), params)
     return main_input
 
+
 #@memoize
 def _flow_subprocess(arguments, main_input):
     filebase, ext = os.path.splitext(os.path.basename(main_input.path))

src/bgem/gmsh/gmsh.py

@@ -530,8 +530,6 @@ def import_shapes(self, fileName, highestDimOnly=True):
         self._need_synchronize = True
         return ObjectSet(self, shapes, [Region.default_region[dim] for dim, _ in shapes])
 
-
-
     def synchronize(self):
         """
         Not clear which actions requires synchronization. Seems that it should be called after calculation of
@@ -631,8 +629,8 @@ def _assign_physical_groups(self, obj):
 
         # set physical groups
         for reg, tags in reg_to_tags.values():
-            reg._gmsh_id = gmsh.model.addPhysicalGroup(reg.dim, tags, tag=-1)
-            gmsh.model.setPhysicalName(reg.dim, reg._gmsh_id, reg.name)
+            gmsh_id = gmsh.model.addPhysicalGroup(reg.dim, tags, tag=-1)
+            gmsh.model.setPhysicalName(reg.dim, gmsh_id, reg.name)
 
 
     def _set_mesh_step(self, obj: 'ObjectSet'):

src/bgem/stochastic/fr_set.py

@@ -992,8 +992,110 @@ def transform_mat(self):
         :return: Shape (N, 3, 3).
         """
         trans_mat = (self.rotation_mat[:, :, :]).copy()
-        trans_mat[:, :, 0] *=  (self.radius[:, 0])[:, None]
-        trans_mat[:, :, 1] *=  (self.radius[:, 1])[:, None]
+        trans_mat[:, :, 0] *= (self.radius[:, 0])[:, None]
+        trans_mat[:, :, 1] *= (self.radius[:, 1])[:, None]
+        return trans_mat
+
+    @staticmethod
+    def get_rotation_mat(normal, shape_axis):
+        """
+        Rotate and scale matrices for the fractures. The full transform involves 'self.center' as well:
+        ambient_space_points = self.center + self.transform_mat @ local_fr_points[:, None, :]
+
+        The Z- local axis is transformed to normal N (assumed unit).
+        The shape_axis S = [sx, sy, 0] must be rotated  to SS by the rotation that rotates Z -> N.
+        Then SS is transformation of the local X axis.
+        The transformation of the Y axis is then computed by the corss product.
+
+        Let's compute S':
+        1. Z -> N rotation unit axis K = [-Ny, Nx, 0] / Nxy
+        2. K . S = Sx Ny - Sy Nx
+        3. follow Rodriguez formula proof, split S into part parallel (p) with K and ortogonal (o) to K
+           Sp = (K.S) K
+           So = S - Sp
+        4. In the plane perpendicular to K,
+           we have vertical component giving: cos(th) = Nz
+           and horizontal component giving sin(th) = Nxy = sqrt(Nx^2 + Ny^2)
+        5. We rotate So by angle th:
+           SSo[z] = -|So| Nxy *sng(Nx)
+           SSo[x,y] = (So [x,y]) Nz
+        6. SSp = Sp
+        7. Sum:
+           SS = SSo + SSp :
+           SSx = Spx + (Nz)Sox
+           SSy = Spy + (Nz)Soy
+           SSz = -|So| Nxy *sng(Nx)
+        Finally, the third vector of the rotated bases is   cross(N, S')
+        TODO: find a better vector representation of the rotations, allowing faster construction of the transfrom matrix
+        :return: Shape (N, 3, 3).
+        """
+        N = normal
+        assert np.allclose(np.linalg.norm(N, axis=1), 1)
+        Nxy = np.stack([-N[:, 1], N[:, 0]], axis=1)
+        norm_Nxy = np.linalg.norm(Nxy, axis=1)
+        K = Nxy / norm_Nxy[:, None]
+        arg_small = np.argwhere(norm_Nxy < 1e-13)[:, 0]
+        K[arg_small, :] = np.array([1, 0], dtype=float)
+        K_dot_S = shape_axis[:, None, :] @ K[:, :, None]
+        S_p = K_dot_S[:, 0, :] * K
+        S_o = shape_axis - S_p
+        cos_th = N[:, 2:3]
+        SS_xy = S_p + cos_th * S_o
+        # ?? th is in (0, pi) -> sin(th) always positive
+        #pos_nx = np.logical_xor(N[:, 0] > 0, N[:, 2] < 0)
+        pos_nx = (N[:, None, 0:2] @ shape_axis[:, :, None])[:, 0, 0] > 0
+        sin_th = norm_Nxy
+        sin_th[pos_nx] = - sin_th[pos_nx]
+        SS_z = np.linalg.norm(S_o, axis=1) * sin_th
+
+        # Construct the rotated X axis SS vector, shape (N, 3)
+        SS = np.concatenate([
+            SS_xy,
+            SS_z[:, None]
+        ], axis=1)
+        scaled_trans_x = SS     #* self.radius[:, 0:1]
+        scaled_trans_y = np.cross(N, SS, axis=1)    #* self.radius[:, 1:2]
+        trans_z = N
+        rot_mat = np.stack([scaled_trans_x, scaled_trans_y, trans_z], axis=2)
+        return rot_mat
+
+    @staticmethod
+    def transform_mat_static(normal, shape_axis, radius):
+        """
+        Rotate and scale matrices for the fractures. The full transform involves 'self.center' as well:
+        ambient_space_points = self.center + self.transform_mat @ local_fr_points[:, None, :]
+
+        The Z- local axis is transformed to normal N (assumed unit).
+        The shape_axis S = [sx, sy, 0] must be rotated  to SS by the rotation that rotates Z -> N.
+        Then SS is transformation of the local X axis.
+        The transformation of the Y axis is then computed by the corss product.
+
+        Let's compute S':
+        1. Z -> N rotation unit axis K = [-Ny, Nx, 0] / Nxy
+        2. K . S = Sx Ny - Sy Nx
+        3. follow Rodriguez formula proof, split S into part parallel (p) with K and ortogonal (o) to K
+           Sp = (K.S) K
+           So = S - Sp
+        4. In the plane perpendicular to K,
+           we have vertical component giving: cos(th) = Nz
+           and horizontal component giving sin(th) = Nxy = sqrt(Nx^2 + Ny^2)
+        5. We rotate So by angle th:
+           SSo[z] = -|So| Nxy *sng(Nx)
+           SSo[x,y] = (So [x,y]) Nz
+        6. SSp = Sp
+        7. Sum:
+           SS = SSo + SSp :
+           SSx = Spx + (Nz)Sox
+           SSy = Spy + (Nz)Soy
+           SSz = -|So| Nxy *sng(Nx)
+        Finally, the third vector of the rotated bases is   cross(N, S')
+        TODO: find a better vector representation of the rotations, allowing faster construction of the transfrom matrix
+        :return: Shape (N, 3, 3).
+        """
+        rotation_mat = FractureSet.get_rotation_mat(normal, shape_axis)
+        trans_mat = (rotation_mat[:, :, :]).copy()
+        trans_mat[:, :, 0] *= (radius[:, 0])[:, None]
+        trans_mat[:, :, 1] *= (radius[:, 1])[:, None]
         return trans_mat
 
     @fn.cached_property

src/bgem/upscale/voxelize.py

@@ -414,7 +414,6 @@ class FracturedMedia:
     fr_conductivity: np.ndarray     # shape (n_fractures,)
     conductivity: float
 
-
     @staticmethod
     def fracture_cond_params(dfn :FractureSet, unit_cross_section, bulk_conductivity):
         # unit_cross_section = 1e-4
@@ -427,8 +426,8 @@ def fracture_cond_params(dfn :FractureSet, unit_cross_section, bulk_conductivity
         # fr cond r=10 ~ 0.8
         fr_r = np.array([fr.r for fr in dfn])
         fr_cross_section = unit_cross_section * fr_r
-        fr_cond = permeability_to_conductivity * permeability_factor * fr_r ** 2
-        fr_cond = np.full_like(fr_r, 10)
+        fr_cond = permeability_to_conductivity * permeability_factor * fr_cross_section ** 2
+        #fr_cond = np.full_like(fr_r, 10)
         return FracturedMedia(dfn, fr_cross_section, fr_cond, bulk_conductivity)
 
     @classmethod