Mod: improved cli commands

Mod: new database entries (porosity, volumeCell, ...) + some fixes
Fix: foamClean function now cleans numbered directories

anisotropy/core/cli.py

@@ -97,7 +97,7 @@ def init(path):
     from anisotropy import env
     from anisotropy.core.main import Database
 
-    if not os.path.exist(path) or not os.path.isdir(path):
+    if not os.path.exists(path) or not os.path.isdir(path):
         click.echo(f"Cannot find directory { path }")
 
         return
@@ -154,12 +154,12 @@ def update(force, params, path):
 
 
     model = Anisotropy()
-    model.db = Database(env["db_name"], env["db_path"]) 
+    database = Database(env["db_name"], env["db_path"]) 
         
     click.echo("Configuring database ...")
-    model.db.setup()
+    database.setup()
 
-    if model.db.isempty() or update:
+    if database.isempty() or update:
         paramsAll = model.loadFromScratch(env["CONFIG"])
 
         if args.get("type"):
@@ -172,7 +172,7 @@ def update(force, params, path):
             paramsAll = [ entry for entry in paramsAll if args["theta"] == entry["structure"]["theta"] ]
 
         for entry in paramsAll:
-            model.db.update(entry)
+            database.update(entry)
 
         click.echo("{} entries was updated.".format(len(paramsAll)))
 
@@ -188,7 +188,7 @@ def update(force, params, path):
 )
 @click.option(
     "-s", "--stage", "stage", 
-    type = click.Choice(["all", "mesh", "flow"]), 
+    type = click.Choice(["all", "mesh", "flow", "postProcessing"]), 
     default = "all",
     help = "Current computation stage"
 )
@@ -219,7 +219,7 @@ def update(force, params, path):
 def compute(stage, nprocs, force, params, path):
     from anisotropy import env
     from anisotropy.core.main import Anisotropy, Database, logger
-    from anisotropy.core.utils import setupLogger, timer, parallel
+    from anisotropy.core.utils import setupLogger, parallel
 
     env.update(
         LOG = os.path.join(path, "logs"),
@@ -236,70 +236,20 @@ def compute(stage, nprocs, force, params, path):
 
     ###
     logger.info("Writing pid ...")
+    pidpath = os.path.join(path, "anisotropy.pid")
 
-    with open(os.path.join(path, "anisotropy.pid"), "w") as io:
+    with open(pidpath, "w") as io:
         io.write(str(os.getpid()))
 
     ###
-    model = Anisotropy()
+    #   Preparations
-    model.db = Database(env["db_name"], env["db_path"]) 
+    ##
+    database = Database(env["db_name"], env["db_path"]) 
         
     logger.info("Loading database ...")
-    model.db.setup()
+    database.setup()
 
-    ###
+    params = database.loadGeneral(
-    def computeCase(stage, type, direction, theta):
-        case = Anisotropy()
-        case.load(type, direction, theta)
-        case.evalParams()
-        case.update()
-
-        logger.info(f"Case: type = { type }, direction = { direction }, theta = { theta }")
-        logger.info(f"Stage: { stage }")
-
-        if stage == "all" or stage == "mesh":
-            if not case.params.get("meshresult", {}).get("meshStatus") == "Done" or force:
-                (out, err, returncode), elapsed = timer(case.computeMesh)(path)
-
-                if out: logger.info(out)
-                if err: logger.error(err)
-
-                case.load(type, direction, theta)
-
-                if case.params.get("meshresult"):
-                    case.params["meshresult"]["meshCalculationTime"] = elapsed
-                    case.update()
-
-                if returncode:
-                    logger.error("Mesh computation failed. Skipping flow computation ...")
-                    return
-
-            else:
-                logger.info("Mesh exists. Skipping ...")
-
-        if stage == "all" or stage == "flow":
-            if not case.params.get("flowresult", {}).get("flowStatus") == "Done" or force:
-                (out, err, returncode), elapsed = timer(case.computeFlow)(path)
-
-                if out: logger.info(out)
-                if err: logger.error(err)
-
-                case.load(type, direction, theta)
-
-                if case.params.get("flowresult"):
-                    case.params["flowresult"]["flowCalculationTime"] = elapsed
-                    case.update()
-            
-                if returncode:
-                    logger.error("Flow computation failed.")
-                    return
-
-            else:
-                logger.info("Flow exists. Skipping ...")
-
-
-    ###
-    params = model.db.loadGeneral(
         args.get("type"), 
         args.get("direction"), 
         args.get("theta")
@@ -309,8 +259,76 @@ def compute(stage, nprocs, force, params, path):
     for p in params:
         s = p["structure"]
 
-        queueargs.append((stage, s["type"], s["direction"], s["theta"]))
+        queueargs.append((s["type"], s["direction"], s["theta"]))
     
+    ###
+    #   Wrap function
+    ##
+    def computeCase(type, direction, theta):
+        case = Anisotropy()
+        case.db = database
+        case.load(type, direction, theta)
+        case.evalParams()
+        case.update()
+
+        logger.info(f"Case: type = { type }, direction = { direction }, theta = { theta }")
+        logger.info(f"Stage mode: { stage }")
+
+        if stage in ["mesh", "all"]:
+            case.load(type, direction, theta)
+
+            if not case.params["meshresult"]["meshStatus"] == "Done" or force:
+                logger.info("Current stage: mesh")
+                out, err, returncode = case.computeMesh(path)
+
+                if out: logger.info(out)
+                if err: logger.error(err)
+                if returncode:
+                    logger.error("Mesh computation failed. Skipping flow computation ...")
+
+                    return
+
+            else:
+                logger.info("Mesh exists. Skipping ...")
+
+        if stage in ["flow", "all"]:
+            case.load(type, direction, theta)
+
+            if not case.params["flowresult"]["flowStatus"] == "Done" or force:
+                logger.info("Current stage: flow")
+                out, err, returncode = case.computeFlow(path)
+
+                if out: logger.info(out)
+                if err: logger.error(err)
+                if returncode:
+                    logger.error("Flow computation failed.")
+
+                    return
+
+            else:
+                logger.info("Flow exists. Skipping ...")
+
+        if stage in ["postProcessing", "all"]:
+            case.load(type, direction, theta)
+
+            if case.params["meshresult"]["meshStatus"] == "Done":
+                logger.info("Current stage: mesh postProcessing")
+                case.porosity()
+
+            else:
+                logger.warning("Cannot compute mesh post processing values.")
+            
+            if case.params["flowresult"]["flowStatus"] == "Done":
+                logger.info("Current stage: flow postProcessing")
+                case.flowRate()
+
+            else:
+                logger.warning("Cannot compute flow post processing values.")
+
+
+    ###
+    #   Run
+    ##
     if nprocs == 1:
         for pos, qarg in enumerate(queueargs):
             computeCase(*qarg)
@@ -318,6 +336,12 @@ def compute(stage, nprocs, force, params, path):
     else:
         parallel(nprocs, queueargs, computeCase)
 
+    if os.path.exists(pidpath):
+        logger.info("Removing pid ...")
+        os.remove(pidpath)
+
+    logger.info("Computation done.")
+
 
 @anisotropy.command(
     help = "Kill process by pid file"
@@ -369,7 +393,7 @@ def computemesh(root, type, direction, theta, path):
     ##
     import os, sys
     
-    pyversion = "{}.{}".format(*sys.version_info[:2])
+    pyversion = "{}.{}".format(3, 9) #(*sys.version_info[:2])
     sys.path.extend([
         root,
         os.path.join(root, "env/lib/python{}/site-packages".format(pyversion)),
@@ -493,14 +517,14 @@ def show(params, path, printlist, export, fields, output):
     else:
         tables.append(df)
 
-    fig, ax = plt.subplots(nrows = 1, ncols = 1)
+    if output == "plot":
+        fig, ax = plt.subplots(nrows = 1, ncols = 1)
 
-    for table in tables:
+        for table in tables:
-        table.plot(table.keys()[0], table.keys()[1], ax = ax, style = "o")
+            table.plot(table.keys()[0], table.keys()[1], ax = ax, style = "o")
 
-    plt.legend()
+        plt.legend()
-    plt.grid()
+        plt.grid()
-    #plt.show()
     
     if export:
         supported = ["csv", "jpg"]
@@ -537,19 +561,19 @@ def show(params, path, printlist, export, fields, output):
 #   CLI entry
 ##
 if __name__ == "__main__":
-    try:
+    #try:
         anisotropy()
 
-    except KeyboardInterrupt:
+    #except KeyboardInterrupt:
-        click.echo("Interrupted!")
+    #    click.echo("Interrupted!")
 
-    finally:
+    #finally:
-        from anisotropy.salomepl.utils import SalomeManager
+    #    from anisotropy.salomepl.utils import SalomeManager
-        click.echo("Exiting ...")
+    #    click.echo("Exiting ...")
 
-        if os.path.exists("anisotropy.pid"):
+    #    if os.path.exists("anisotropy.pid"):
-            os.remove("anisotropy.pid")
+    #        os.remove("anisotropy.pid")
 
-        SalomeManager().killall()
+    #    SalomeManager().killall()
 
-        sys.exit(0)
+    #    sys.exit(0)

anisotropy/core/main.py

@@ -16,7 +16,7 @@ from anisotropy import (
     __version__, env,
     openfoam
 )
-from anisotropy.core.utils import setupLogger, timer
+from anisotropy.core.utils import setupLogger, Timer
 from anisotropy.core.database import Database
 from anisotropy import salomepl
 import anisotropy.salomepl.utils
@@ -38,7 +38,7 @@ class Anisotropy(object):
         """Constructor method"""
 
         self.env = env
-        self.db = Database(self.env["db_name"], self.env["db_path"])
+        self.db = None #Database(self.env["db_name"], self.env["db_path"])
         self.params = []
 
 
@@ -63,8 +63,8 @@ class Anisotropy(object):
     def version():
         """Returns versions of all used main programs
 
-        :return: Versions joined by next line symbol
+        :return: 
-        :rtype: str
+            Versions joined by next line symbol
         """
         versions = {
             "anisotropy": __version__,
@@ -86,8 +86,8 @@ class Anisotropy(object):
     def loadFromScratch(self, configpath: str = None) -> list:
         """Loads parameters from configuration file and expands special values
 
-        :return: List of dicts with parameters
+        :return: 
-        :rtype: list
+            List of dicts with parameters
         """
         config = configpath or self.env["CONFIG"]
 
@@ -131,8 +131,10 @@ class Anisotropy(object):
                         ),
                         "mesh": mesh,
                         "submesh": deepcopy(entry["submesh"]),
+                        "meshresult": dict(),
                         "flow": deepcopy(entry["flow"]),
-                        "flowapproximation": deepcopy(entry["flowapproximation"])
+                        "flowapproximation": deepcopy(entry["flowapproximation"]),
+                        "flowresult": dict(),
                     }
 
                     # For `type = fixedValue` only
@@ -251,13 +253,33 @@ class Anisotropy(object):
         manager = salomepl.utils.SalomeManager()
         casepath = self.getCasePath(path)
 
-        return manager.execute(
+        self.params["meshresult"]["meshStatus"] = "Computing"
+        self.update()
+        timer = Timer()
+
+        out, err, returncode = manager.execute(
             scriptpath, 
             *salomeargs, 
             timeout = self.env["salome_timeout"],
             root = self.env["ROOT"],
             logpath = casepath
         )
+        self.load(p["type"], p["direction"], p["theta"])
+
+        if not returncode:
+            self.params["meshresult"].update(
+                meshStatus = "Done",
+                meshCalculationTime = timer.elapsed()
+            )
+
+        else:
+            self.params["meshresult"].update(
+                meshStatus = "Failed"
+            )
+
+        self.update()
+
+        return out, err, returncode
 
 
     def genmesh(self, path):
@@ -280,7 +302,8 @@ class Anisotropy(object):
             bodyCentered = BodyCentered,
             faceCentered = FaceCentered
         )[p["structure"]["type"]]
-        shape, groups = structure(**p["structure"]).build()
+        shapeGeometry = structure(**p["structure"])
+        shape, groups = shapeGeometry.build()
 
         [length, surfaceArea, volume] = geompy.BasicProperties(shape, theTolerance = 1e-06)
 
@@ -312,6 +335,7 @@ class Anisotropy(object):
         if mp["viscousLayers"]:
             mesh.ViscousLayers(**mp, faces = faces)
 
+        #   Submesh
         smp = p["submesh"]
 
         for submesh in smp:
@@ -333,7 +357,7 @@ class Anisotropy(object):
         ###
         #   Results
         ##
-        p["meshresult"] = dict()
+        #p["meshresult"] = dict()
 
         if not returncode:
             mesh.removePyramids()
@@ -349,9 +373,9 @@ class Anisotropy(object):
 
             meshStats = mesh.stats()
             p["meshresult"].update(
-                status = "Done",
                 surfaceArea = surfaceArea,
                 volume = volume,
+                volumeCell = shapeGeometry.volumeCell,
                 **meshStats
             )
             self.update()
@@ -360,9 +384,9 @@ class Anisotropy(object):
             logger.error(err)
 
             p["meshresult"].update(
-                status = "Failed",
                 surfaceArea = surfaceArea,
-                volume = volume
+                volume = volume,
+                volumeCell = shapeGeometry.volumeCell
             )
             self.update()
 
@@ -370,23 +394,30 @@ class Anisotropy(object):
     def computeFlow(self, path):
         """Computes a flow on mesh via OpenFOAM
 
-        :return: Process output, error messages and returncode
+        :return: 
-        :rtype: tuple(str, str, int)
+            Process output, error messages and returncode
         """
         ###
         #   Case preparation
         ##
         foamCase = [ "0", "constant", "system" ]
+        #self.params["flowresult"] = dict()
+        self.params["flowresult"]["flowStatus"] = "Computing"
+        self.update()
+        timer = Timer()
 
         flow = self.params["flow"]
         flowapproximation = self.params["flowapproximation"]
 
         # ISSUE: ideasUnvToFoam cannot import mesh with '-case' flag so 'os.chdir' for that
-        casePath = self.getCasePath()
+        casePath = self.getCasePath(path)
 
         if not os.path.exists(casePath):
-            logger.warning(f"Cannot find case path. Skipping computation ...\n\t{ casePath }")
+            err = f"Cannot find case path { casePath }"
-            return "", "", 1
+            self.params["flowresult"]["flowStatus"] = "Failed"
+            self.update()
+
+            return "", err, 1
 
         os.chdir(casePath)
         openfoam.foamClean()
@@ -401,14 +432,22 @@ class Anisotropy(object):
         #   Mesh manipulations
         ##
         if not os.path.exists("mesh.unv"):
-            logger.error(f"missed 'mesh.unv'")
             os.chdir(path or self.env["ROOT"])
-            return "", "", 1
+
+            err = f"Missed 'mesh.unv'"
+            self.params["flowresult"]["flowStatus"] = "Failed"
+            self.update()
+
+            return "", err, 1
 
         out, err, returncode = openfoam.ideasUnvToFoam("mesh.unv")
 
         if returncode:
             os.chdir(path or self.env["ROOT"])
+
+            self.params["flowresult"]["flowStatus"] = "Failed"
+            self.update()
+            
             return out, err, returncode
         
         openfoam.createPatch(dictfile = "system/createPatchDict")
@@ -424,10 +463,9 @@ class Anisotropy(object):
             "1 (wall)"
         )
         
-        out = openfoam.checkMesh()
+        out, err, returncode = openfoam.checkMesh()
         
-        if out:
+        if out: logger.warning(out)
-            logger.info(out)
         
         openfoam.transformPoints(flow["scale"])
         
@@ -442,7 +480,7 @@ class Anisotropy(object):
             str(flow["transportProperties"]["nu"])
         )
 
-        #openfoam.decomposePar()
+        # openfoam.decomposePar()
 
         openfoam.renumberMesh()
 
@@ -499,33 +537,44 @@ class Anisotropy(object):
         
         out, err, returncode = openfoam.simpleFoam()
 
-        ###
-        #   Results
-        ##
-        self.params["flowresult"] = dict()
-
         if not returncode:
-            postProcessing = "postProcessing/flowRatePatch(name=outlet)/0/surfaceFieldValue.dat"
+            self.params["flowresult"]["flowCalculationTime"] = timer.elapsed()
-
+            self.params["flowresult"]["flowStatus"] = "Done"
-            with open(os.path.join(casePath, postProcessing), "r") as io:
-                lastLine = io.readlines()[-1]
-                flowRate = float(lastLine.replace(" ", "").replace("\n", "").split("\t")[1])
-                
-            self.params["flowresult"].update(
-                status = "Done",
-                flowRate = flowRate
-            )
 
         else:
-            self.params["flowresult"].update(
+            self.params["flowresult"]["flowStatus"] = "Failed"
-                status = "Failed"
-            )
-
 
         self.update()
-
+        os.chdir(path or self.env["ROOT"])
-        os.chdir(self.env["ROOT"])
         
         return out, str(err, "utf-8"), returncode
 
 
+    def flowRate(self):
+        casePath = self.getCasePath()
+        foamPostProcessing = "postProcessing/flowRatePatch(name=outlet)/0/surfaceFieldValue.dat"
+        path = os.path.join(casePath, foamPostProcessing)
+
+        if not os.path.exists(path):
+            logger.warning(f"Unable to compute flow rate. Missed { path }")
+
+            return
+
+        with open(path, "r") as io:
+            lastLine = io.readlines()[-1]
+            flowRate = float(lastLine.replace(" ", "").replace("\n", "").split("\t")[1])
+
+        self.params["flowresult"]["flowRate"] = flowRate
+        self.update()
+
+        return flowRate
+
+
+    def porosity(self):
+        mr = self.params["meshresult"] 
+        fr = self.params["flowresult"]
+
+        fr["porosity"] = mr["volume"] / mr["volumeCell"]
+        self.update()
+
+        return fr["porosity"]

anisotropy/core/models.py

@@ -129,6 +129,7 @@ class MeshResult(BaseModel):
     
     surfaceArea = FloatField(null = True)
     volume = FloatField(null = True)
+    volumeCell = FloatField(null = True)
     
     elements = IntegerField(null = True)
     edges = IntegerField(null = True)

anisotropy/core/utils.py

@@ -205,6 +205,13 @@ def timer(func: FunctionType) -> (tuple, float):
     return inner
 
 
+class Timer(object):
+    def __init__(self):
+        self.start = time.monotonic()
+
+    def elapsed(self):
+        return time.monotonic() - self.start
+
 
 def queue(cmd, qin, qout, *args):
     

anisotropy/openfoam/meshManipulation.py

@@ -22,7 +22,7 @@ def transformPoints(scale, case: str = None):
 
 
 def checkMesh(case: str = None) -> str:
-    application("checkMesh", "-allGeometry", "-allTopology", case = case, stderr = True)
+    _, err, returncode = application("checkMesh", "-allGeometry", "-allTopology", case = case, stderr = True)
     out = ""
 
     with open("checkMesh.log", "r") as io:
@@ -34,7 +34,7 @@ def checkMesh(case: str = None) -> str:
         if warnings:
             out = "checkMesh:\n\t{}".format("\n\t".join(warnings))
 
-    return out
+    return out, err, returncode
 
 
 def renumberMesh(case: str = None):

anisotropy/openfoam/utils.py

@@ -4,12 +4,13 @@
 
 import os
 import shutil
+from .application import application
 
 def version() -> str:
     return os.environ["WM_PROJECT_VERSION"]
 
 
-def foamClean(case: str = None):
+def foamCleanCustom(case: str = None):
     rmDirs = ["0", "constant", "system", "postProcessing", "logs"]
     rmDirs.extend([ "processor{}".format(n) for n in range(os.cpu_count()) ])
     path = case if case else ""
@@ -18,6 +19,15 @@ def foamClean(case: str = None):
         if os.path.exists(os.path.join(path, d)):
             shutil.rmtree(os.path.join(path, d))
 
+def foamClean(case: str = None):
+    rmDirs = ["0", "constant", "system"]
+    path = case if case else ""
+
+    for d in rmDirs:
+        if os.path.exists(os.path.join(path, d)):
+            shutil.rmtree(os.path.join(path, d))
+
+    application("foamCleanTutorials", useMPI = False, case = case, stderr = True)
 
 def uniform(value) -> str:
     if type(value) == list or type(value) == tuple:

anisotropy/samples/bodyCentered.py

@@ -15,6 +15,7 @@ class BodyCentered(object):
         self.radius = kwargs.get("radius", self.r0 / (1 - self.theta)) 
         self.filletsEnabled = kwargs.get("filletsEnabled", False)
         self.fillets = kwargs.get("fillets", 0)
+        self.volumeCell = None
 
 
     def build(self):
@@ -72,6 +73,7 @@ class BodyCentered(object):
                 vecflow = geompy.GetNormal(inletface)
                 poreCell = geompy.MakePrismVecH(inletface, vecflow, zh)
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
             
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
@@ -127,6 +129,8 @@ class BodyCentered(object):
             vecflow = geompy.GetNormal(inletface)
             poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
+
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
 

anisotropy/samples/faceCentered.py

@@ -15,6 +15,7 @@ class FaceCentered(object):
         self.radius = kwargs.get("radius", self.r0 / (1 - self.theta)) 
         self.filletsEnabled = kwargs.get("filletsEnabled", False)
         self.fillets = kwargs.get("fillets", 0)
+        self.volumeCell = None
 
 
     def build(self):
@@ -72,6 +73,7 @@ class FaceCentered(object):
                 vecflow = geompy.GetNormal(inletface)
                 poreCell = geompy.MakePrismVecH(inletface, vecflow, 2 * zh)
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
             
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
@@ -127,6 +129,8 @@ class FaceCentered(object):
             vecflow = geompy.GetNormal(inletface)
             poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
+
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
 

anisotropy/samples/simple.py

@@ -15,6 +15,7 @@ class Simple(object):
         self.radius = kwargs.get("radius", self.r0 / (1 - self.theta)) 
         self.filletsEnabled = kwargs.get("filletsEnabled", False)
         self.fillets = kwargs.get("fillets", 0)
+        self.volumeCell = None
 
 
     def build(self):
@@ -68,6 +69,7 @@ class Simple(object):
                 vecflow = geompy.GetNormal(inletface)
                 poreCell = geompy.MakePrismVecH(inletface, vecflow, height)
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
                 
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)
@@ -120,6 +122,8 @@ class Simple(object):
             vecflow = geompy.GetNormal(inletface)
             poreCell = geompy.MakePrismVecH(inletface, vecflow, self.L * sqrt(3))
 
+            [_, _, self.volumeCell] = geompy.BasicProperties(poreCell, theTolerance = 1e-06)
+
             inletface = geompy.MakeScaleTransform(inletface, oo, scale)
             poreCell = geompy.MakeScaleTransform(poreCell, oo, scale)