Merge branch 'dev' into 'main'

Merged dev branch

See merge request jcfergus/puma-dev!14

python/pumapy/__init__.py

@@ -7,10 +7,19 @@
 THE UNITED STATES GOVERNMENT DISCLAIMS ALL WARRANTIES AND LIABILITIES REGARDING THIRD PARTY COMPUTER SOFTWARE, DATA, OR DOCUMENTATION, IF SAID THIRD PARTY COMPUTER SOFTWARE, DATA, OR DOCUMENTATION IS PRESENT IN THE NASA SOFTWARE AND/OR TECHNICAL DATA, AND DISTRIBUTES IT "AS IS."
 RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED STATES GOVERNMENT AND ITS CONTRACTORS AND SUBCONTRACTORS, AND SHALL INDEMNIFY AND HOLD HARMLESS THE UNITED STATES GOVERNMENT AND ITS CONTRACTORS AND SUBCONTRACTORS FOR ANY LIABILITIES, DEMANDS, DAMAGES, EXPENSES OR LOSSES THAT MAY ARISE FROM RECIPIENT'S USE OF THE SOFTWARE AND/OR TECHNICAL DATA, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON, OR RESULTING FROM, THE USE THEREOF.
 IF RECIPIENT FURTHER RELEASES OR DISTRIBUTES THE NASA SOFTWARE AND/OR TECHNICAL DATA, RECIPIENT AGREES TO OBTAIN THIS IDENTICAL WAIVER OF CLAIMS, INDEMNIFICATION AND HOLD HARMLESS, AGREEMENT WITH ANY ENTITIES THAT ARE PROVIDED WITH THE SOFTWARE AND/OR TECHNICAL DATA.
-"""
 
-""" pumapy
 Root directory for the pumapy package.
+Please refer to this publication for a detailed software architecture explanation:
+
+@article{ferguson2021update,
+  title={Update 3.0 to “puma: The porous microstructure analysis software”,(pii: s2352711018300281)},
+  author={Ferguson, Joseph C and Semeraro, Federico and Thornton, John M and Panerai, Francesco and Borner, Arnaud and Mansour, Nagi N},
+  journal={SoftwareX},
+  volume={15},
+  pages={100775},
+  year={2021},
+  publisher={Elsevier}
+}
 """
 
 
@@ -22,15 +31,16 @@
 # - git push nasa main
 # - git tag -a v$(python setup.py --version) -m 'INPUT DESCRIPTION'
 # - gh release create v$(python setup.py --version) --target main
-__version__ = "3.1.5"
+__version__ = "3.1.6"
 
 
 # utilities
 from pumapy.utilities.workspace import Workspace
 from pumapy.utilities.timer import Timer
-from pumapy.utilities.property_maps import IsotropicConductivityMap, AnisotropicConductivityMap, ElasticityMap
-from pumapy.utilities.boundary_conditions import ConductivityBC, ElasticityBC
+from pumapy.physicsmodels.property_maps import IsotropicConductivityMap, AnisotropicConductivityMap, ElasticityMap
+from pumapy.physicsmodels.boundary_conditions import ConductivityBC, ElasticityBC
 from pumapy.utilities.example_files import path_to_example_file, list_example_files
+from pumapy.utilities.generic_checks import estimate_max_memory, set_random_seed
 
 # input/output
 from pumapy.io.input import import_3Dtiff, import_bin, import_weave_vtu, import_vti
@@ -47,7 +57,7 @@
 from pumapy.materialproperties.orientation import compute_orientation_st, compute_angular_differences
 from pumapy.materialproperties.conductivity import compute_thermal_conductivity, compute_electrical_conductivity
 from pumapy.materialproperties.tortuosity import compute_continuum_tortuosity
-from pumapy.materialproperties.elasticity import compute_elasticity, compute_stress_analysis
+from pumapy.materialproperties.elasticity import compute_elasticity, compute_stress_analysis, get_E_nu_from_elasticity
 from pumapy.materialproperties.radiation import compute_radiation, compute_extinction_coefficients
 from pumapy.materialproperties.permeability import compute_permeability
 
@@ -72,7 +82,7 @@
 from pumapy.visualization.slicer import plot_slices, compare_slices
 
 # segmentation
-from pumapy.segmentation.porespace import identify_porespace, fill_closed_pores
+from pumapy.segmentation.ccl import identify_porespace, fill_closed_pores, remove_rbms
 
 # global settings
 settings = {"log_location": 'logs'}

python/pumapy/io/input.py

@@ -223,7 +223,7 @@ def import_weave_vtu(filename, from_texgen_gui=False):
     yarn_index = vtk_to_numpy(vtkobject.GetCellData().GetArray(0)) + 1
     ws = Workspace.from_array(yarn_index.reshape(int(dims[0]), int(dims[1]), int(dims[2]), order="F"))
 
-    if vtkobject.GetCellData().GetNumberOfArrays() == 2:
+    if vtkobject.GetCellData().GetNumberOfArrays() > 2:
         if from_texgen_gui:
             # ORIGINAL TEXGEN (GUI in Windows)
             # Number Of Arrays: 6

python/pumapy/materialproperties/conductivity.py

@@ -1,6 +1,6 @@
 from pumapy.physicsmodels.isotropic_conductivity import IsotropicConductivity
 from pumapy.physicsmodels.mpfa_conductivity import AnisotropicConductivity
-from pumapy.utilities.property_maps import IsotropicConductivityMap, AnisotropicConductivityMap
+from pumapy.physicsmodels.property_maps import IsotropicConductivityMap, AnisotropicConductivityMap
 
 
 def compute_thermal_conductivity(workspace, cond_map, direction, side_bc='s', prescribed_bc=None, tolerance=1e-4,
@@ -16,12 +16,12 @@ def compute_thermal_conductivity(workspace, cond_map, direction, side_bc='s', pr
         :param side_bc: side boundary conditions can be symmetric ('s'), periodic ('p') or dirichlet ('d')
         :type side_bc: string
         :param prescribed_bc: 3D array holding dirichlet BC.
-        :type prescribed_bc: pumapy.ConductivityBC
+        :type prescribed_bc: pumapy.ConductivityBC or None
         :param tolerance: tolerance for iterative solver
         :type tolerance: float
         :param maxiter: maximum Iterations for solver
         :type maxiter: int
-        :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres', 'direct'
+        :param solver_type: solver type, options: 'bicgstab' (default), 'cg', 'gmres', 'direct'
         :type solver_type: string
         :param display_iter: display iterations and residual
         :type display_iter: bool
@@ -32,13 +32,23 @@ def compute_thermal_conductivity(workspace, cond_map, direction, side_bc='s', pr
 
         :Example:
         >>> import pumapy as puma
-        >>> ws_fiberform = puma.import_3Dtiff(puma.path_to_example_file("100_fiberform.tif"), 1.3e-6)
+        >>> ws_fiberform = puma.import_3Dtiff(puma.path_to_example_file("200_fiberform.tif"), 1.3e-6)
+        >>> ws_fiberform.rescale(0.5, segmented=False)
+        >>>
+        >>> # Conductivity with Isotropic local phases
         >>> cond_map = puma.IsotropicConductivityMap()
         >>> cond_map.add_material((0, 89), 0.0257)
         >>> cond_map.add_material((90, 255), 12)
-        >>> k_eff_x, T_x, q_x = puma.compute_thermal_conductivity(ws_fiberform, cond_map, 'x', 's', tolerance=1e-2, solver_type='cg')
-        >>> print("Effective thermal conductivity tensor:")
-        >>> print(k_eff_x)
+        >>> k_eff_x, T_x, q_x = puma.compute_thermal_conductivity(ws_fiberform, cond_map, 'x', 's')
+        >>>
+        >>> # Conductivity with Anisotropic local phases
+        >>> puma.compute_orientation_st(ws_fiberform, cutoff=(90, 255))
+        >>> cond_map = puma.AnisotropicConductivityMap()
+        >>> # conductivity of the void phase to be 0.0257 (air at STP)
+        >>> cond_map.add_isotropic_material((0, 89), 0.0257)
+        >>> # axial fiber conductivity of 12, radial fiber conductivity of 0.7
+        >>> cond_map.add_material_to_orient((90, 255), 12., 0.7)
+        >>> k_eff_z, T_z, q_z = puma.compute_thermal_conductivity(ws_fiberform, cond_map, 'z', 's')
     """
     if isinstance(cond_map, IsotropicConductivityMap):
         solver = IsotropicConductivity(workspace, cond_map, direction, side_bc, prescribed_bc, tolerance, maxiter,
@@ -70,12 +80,12 @@ def compute_electrical_conductivity(workspace, cond_map, direction, side_bc='p',
     :param side_bc: side boundary conditions can be symmetric ('s'), periodic ('p') or dirichlet ('d')
     :type side_bc: string
     :param prescribed_bc: 3D array holding dirichlet BC
-    :type prescribed_bc: pumapy.ConductivityBC
+    :type prescribed_bc: pumapy.ConductivityBC or None
     :param tolerance: tolerance for iterative solver
     :type tolerance: float
     :param maxiter: maximum Iterations for solver
     :type maxiter: int
-    :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres', 'direct'
+    :param solver_type: solver type, options: 'bicgstab' (default), 'cg', 'gmres', 'direct'
     :type solver_type: string
     :param display_iter: display iterations and residual
     :type display_iter: bool
@@ -92,8 +102,6 @@ def compute_electrical_conductivity(workspace, cond_map, direction, side_bc='p',
     >>> cond_map.add_material((0, 89), 0.0257)
     >>> cond_map.add_material((90, 255), 12)
     >>> k_eff_x, P_x, q_x = puma.compute_electrical_conductivity(ws_fiberform, cond_map, 'x', 's', tolerance=1e-2, solver_type='cg')
-    >>> print("Effective electrical conductivity tensor:")
-    >>> print(k_eff_x)
     """
     return compute_thermal_conductivity(workspace, cond_map, direction, side_bc, prescribed_bc, tolerance, maxiter,
                                         solver_type, display_iter, print_matrices)

python/pumapy/materialproperties/elasticity.py

@@ -1,10 +1,13 @@
 from pumapy.physicsmodels.mpsa_elasticity import Elasticity
-from pumapy.utilities.property_maps import ElasticityMap
+from pumapy.physicsmodels.property_maps import ElasticityMap
+from scipy.optimize import fsolve
+import numpy as np
 
 
 def compute_elasticity(workspace, elast_map, direction, side_bc='p', prescribed_bc=None, tolerance=1e-4,
-                       maxiter=10000, solver_type='bicgstab', display_iter=True, print_matrices=(0, 0, 0, 0, 0)):
-    """ Compute the thermal conductivity (N.B. 0 material ID in workspace refers to air unless otherwise specified)
+                       maxiter=100000, solver_type='gmres', display_iter=True, print_matrices=(0, 0, 0, 0, 0)):
+    """ Compute the effective elasticity coefficient
+        (N.B. 0 material ID in workspace refers to air unless otherwise specified)
 
         :param workspace: domain
         :type workspace: pumapy.Workspace
@@ -15,12 +18,12 @@ def compute_elasticity(workspace, elast_map, direction, side_bc='p', prescribed_
         :param side_bc: side boundary conditions can be symmetric ('s'), periodic ('p'), dirichlet ('d') or free ('f')
         :type side_bc: string
         :param prescribed_bc: 3D array holding dirichlet BC
-        :type prescribed_bc: pumapy.ElasticityBC
+        :type prescribed_bc: pumapy.ElasticityBC or None
         :param tolerance: tolerance for iterative solver
         :type: tolerance: float
         :param maxiter: maximum Iterations for solver
         :type maxiter: int
-        :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres', 'direct'
+        :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres' (default), 'direct'
         :type solver_type: string
         :param display_iter: display iterations and residual
         :type display_iter: bool
@@ -54,8 +57,9 @@ def compute_elasticity(workspace, elast_map, direction, side_bc='p', prescribed_
 
 
 def compute_stress_analysis(workspace, elast_map, prescribed_bc, side_bc='p', tolerance=1e-4,
-                            maxiter=10000, solver_type='bicgstab', display_iter=True, print_matrices=(0, 0, 0, 0, 0)):
-    """ Compute the thermal conductivity (N.B. 0 material ID in workspace refers to air unless otherwise specified)
+                            maxiter=100000, solver_type='gmres', display_iter=True, print_matrices=(0, 0, 0, 0, 0)):
+    """ Compute stress analysis
+        (N.B. 0 material ID in workspace refers to air unless otherwise specified)
 
         :param workspace: domain
         :type workspace: pumapy.Workspace
@@ -69,7 +73,7 @@ def compute_stress_analysis(workspace, elast_map, prescribed_bc, side_bc='p', to
         :type tolerance: float
         :param maxiter: maximum Iterations for solver
         :type maxiter: int
-        :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres', 'direct'
+        :param solver_type: solver type, options: 'bicgstab', 'cg', 'gmres' (default), 'direct'
         :type solver_type: string
         :param display_iter: display iterations and residual
         :type display_iter: bool
@@ -85,10 +89,10 @@ def compute_stress_analysis(workspace, elast_map, prescribed_bc, side_bc='p', to
         >>> elast_map = puma.ElasticityMap()
         >>> elast_map.add_isotropic_material((1, 1), 200, 0.3)
         >>> elast_map.add_isotropic_material((2, 2), 400, 0.1)
-        >>> bc = puma.ElasticityBC.from_workspace(ws)
-        >>> bc[0] = 0  # hold x -ve face
-        >>> bc[-1, :, :, 0] = 1   # displace x +ve face by 1 in x direction
-        >>> bc[-1, :, :, 1:] = 0  # hold x +ve face in y and z directions
+        >>> bc = puma.ElasticityBC(ws)
+        >>> bc.dirichlet[0] = 0  # hold x -ve face
+        >>> bc.dirichlet[-1, :, :, 0] = 1   # displace x +ve face by 1 in x direction
+        >>> bc.dirichlet[-1, :, :, 1:] = 0  # hold x +ve face in y and z directions
         >>> u, s, t = puma.compute_stress_analysis(ws, elast_map, bc, side_bc='f', solver_type="direct")
         >>> puma.render_volume(u[:, :, :, 1], cmap='jet')  # displacement magnitude in y direction
     """
@@ -104,3 +108,39 @@ def compute_stress_analysis(workspace, elast_map, prescribed_bc, side_bc='p', to
     solver.compute()
     solver.log_output()
     return solver.u, solver.s, solver.t
+
+
+def get_E_nu_from_elasticity(C11, C12, C13, C22, C23, C33):
+    """ Compute Young's moduli E1, E2, E3 and Poisson's ratios nu12, nu23, nu31 from symmetric elastic stiffness tensor
+        :param C11: elasticity tensor component
+        :type C11: float
+        :param C12: elasticity tensor component
+        :type C12: float
+        :param C13: elasticity tensor component
+        :type C13: float
+        :param C22: elasticity tensor component
+        :type C22: float
+        :param C23: elasticity tensor component
+        :type C23: float
+        :param C33: elasticity tensor component
+        :type C33: float
+        :return: Young's moduli E1, E2, E3 and Poisson's ratios nu12, nu23, nu31
+        :rtype: (float, float, float, float, float, float)
+    """
+    def eqs(unknowns, *cs):
+        E1, E2, E3, nu12, nu23, nu31 = unknowns
+        c11, c12, c31, c22, c23, c33 = cs
+        nu21 = nu12 * E2 / E1
+        nu32 = nu23 * E3 / E2
+        nu13 = nu31 * E1 / E3
+        c0 = (1 - nu12 * nu21 - nu23 * nu32 - nu13 * nu31 - 2 * nu21 * nu32 * nu13) / (E1 * E2 * E3)
+        return (c11 - (1 - nu23 * nu32) / (c0 * E2 * E3),
+                c22 - (1 - nu13 * nu31) / (c0 * E1 * E3),
+                c33 - (1 - nu21 * nu12) / (c0 * E2 * E1),
+                c12 - (nu12 + nu32 * nu13) / (c0 * E1 * E3),
+                c23 - (nu23 + nu21 * nu13) / (c0 * E1 * E2),
+                c31 - (nu31 + nu21 * nu32) / (c0 * E2 * E3))
+
+    s = fsolve(eqs, np.full(6, 0.1), args=(C11, C12, C13, C22, C23, C33))
+    [print(i, j) for i, j in zip(["E1", "E2", "E3", "nu12", "nu23", "nu31"], s)]
+    return s

python/pumapy/materialproperties/orientation.py

@@ -1,3 +1,22 @@
+"""
+Further explained in publication:
+Semeraro, F., Ferguson, J.C., Panerai, F., King, R.J. and Mansour, N.N., 2020.
+Anisotropic analysis of fibrous and woven materials part 1: Estimation of local orientation.
+Computational Materials Science, 178, p.109631.
+(https://www.sciencedirect.com/science/article/abs/pii/S0927025620301221)
+
+Please cite using this BibTex:
+@article{semeraro2020anisotropic,
+  title={Anisotropic analysis of fibrous and woven materials part 1: Estimation of local orientation},
+  author={Semeraro, Federico and Ferguson, Joseph C and Panerai, Francesco and King, Robert J and Mansour, Nagi N},
+  journal={Computational Materials Science},
+  volume={178},
+  pages={109631},
+  year={2020},
+  publisher={Elsevier}
+}
+"""
+
 import sys
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter

python/pumapy/materialproperties/permeability.py

@@ -14,9 +14,10 @@ def compute_permeability(workspace, solid_cutoff, tol=1e-8, maxiter=10000, solve
         :type maxiter: int
         :param solver_type: solver type, options: 'minres' (default), 'direct', 'cg', 'bicgstab'
         :type solver_type: string
+        :param display_iter: display iteration in iterative solver
         :type display_iter: bool
         :return: effective permeability (3x3 matrix), velocity and pressure fields for x, y, z directions (u_x, p_x, etc)
-        :rtype: (numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray)
+        :rtype: numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray
 
         :Example:
         >>> import pumapy as puma

python/pumapy/materialproperties/radiation.py

@@ -1,4 +1,4 @@
-from pumapy.utilities.raycasting import RayCasting
+from pumapy.physicsmodels.raycasting import RayCasting
 from scipy.optimize import curve_fit
 from matplotlib import pyplot as plt
 import numpy as np

python/pumapy/materialproperties/tortuosity.py

@@ -1,5 +1,5 @@
 from pumapy.materialproperties.volumefraction import compute_volume_fraction
-from pumapy.utilities.property_maps import IsotropicConductivityMap
+from pumapy.physicsmodels.property_maps import IsotropicConductivityMap
 from pumapy.physicsmodels.isotropic_conductivity import IsotropicConductivity
 
 
@@ -16,12 +16,12 @@ def compute_continuum_tortuosity(workspace, cutoff, direction, side_bc='p', pres
         :param side_bc: side boundary conditions (string) can be symmetric ('s'), periodic ('p') or dirichlet ('d')
         :type side_bc: string
         :param prescribed_bc: 3D array holding dirichlet BC
-        :type prescribed_bc: pumapy.ConductivityBC
+        :type prescribed_bc: pumapy.ConductivityBC or None
         :param tolerance: tolerance for iterative solver
         :type tolerance: float
         :param maxiter: maximum Iterations for solver
         :type maxiter: int
-        :param solver_type: solver type, options: 'cg', 'bicgstab', 'direct'
+        :param solver_type: solver type, options: 'cg' (default), 'bicgstab', 'direct'
         :type solver_type: string
         :param display_iter: display iterations and residual
         :type display_iter: bool