Connector/Generator: IBM AMR (#611)

* MeshGen & IBMprep python wrapper: bug fixes

* Bug fixes

* generalize if statement - empty t

* Apply autopep8 formatting

---------

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>

Author: antjost

Cassiopee/Connector/Connector/AMR.py

@@ -245,7 +245,14 @@ def prepareAMRData(t_case, t, IBM_parameters=None, check=False, dim=3, localDir=
         if different_front_flag == False: #True is default
             Internal._rmNodesFromName(t, "TurbulentDistance")
             Internal._renameNode(t, "TurbulentDistanceForCFDComputation","TurbulentDistance")
-    if OPT:
+
+    # Catch 22: IBM prep needs dist2wall at the nodes (dist2wall@Node). blanking, ibm point location, etc. is doing using dist2wall@Node
+    #           CODA needs dist2wall at the cell centers. CODA doesnt recalculate this for IBM runs and must therefore be provided.
+    # Chosen approach: OPT=True - use only the dist2wall@Node for the IBM prep & if needed recalculate the dist2wall@Centers at the end of the function.
+    #                  Options for OPT=False have not been tested but do not seem like the appropriate solution : calcs on dist2wall@Centers + center2Node is counterintuitive
+    # TODO: Remove OPT=False options ONLY after testing with IBM Local options. F. Basile prototype did some operations with dist2wall@Centers.
+    node = Internal.getZones(t) #select cells after blanking may yields an empty zone - yields empty list or list with elements
+    if OPT and node:
         varnames = C.getVarNames(t, loc="centers")[0]
         if "TurbulentDistance" not in varnames:
             DTW._distance2Walls(t, tb2, type='ortho', signed=0, dim=dim, loc='centers')
@@ -979,19 +986,41 @@ def prepareAMRIBM(tb, levelMax, vmins, dim, IBM_parameters, toffset=None, check=
     Usage: prepareAMRIBM(tb, levelMax, vmins, dim, IBM_parameters, toffset, check, opt, octreeMode,
                          snears, dfars, loadBalancing, OutputAMRMesh, localDir, fileName, tbox, vminsTbox, tbv2, forceAlignment)"""
 
+    ## =========================
+    ## ==== Mesh Generation ====
+    ## =========================
     t_AMR = G_AMR.generateAMRMesh(tb=tb, levelMax=levelMax, vmins=vmins, dim=dim,
                                   toffset=toffset, check=check, opt=opt, octreeMode=octreeMode, localDir=localDir,
                                   snears=0.01, dfars=10, loadBalancing=loadBalancing, tbox=tbox, vminsTbox=vminsTbox, tbv2=tbv2)
 
     Cmpi.trace('AMR Mesh Dist2Walls...start', master=True)
-    if dim == 2: tb2 = T.addkplane(tb)
+    tb2  = Internal.copyTree(tb)
+    if dim == 2:tb2 = T.addkplane(tb)
+    else:
+        baseSYM = Internal.getNodesFromName1(tb2, "SYM")
+        if baseSYM:
+            tb2 = Internal.rmNodesByNameAndType(tb2, 'SYM', 'CGNSBase_t')
+            tb2 = Internal.rmNodesByNameAndType(tb2, '*_sym*', 'Zone_t')
+            tb  = Internal.copyTree(tb2)
     DTW._distance2Walls(t_AMR, tb2, type='ortho', signed=0, dim=dim, loc='centers')
     DTW._distance2Walls(t_AMR, tb2, type='ortho', signed=0, dim=dim, loc='nodes')
     del tb2
     Cmpi.trace('AMR Mesh Dist2Walls...end', master=True)
 
     if OutputAMRMesh: Cmpi.convertPyTree2File(t_AMR, localDir+'tAMRMesh.cgns')
+    ## Ncells output
+    Ncells=C.getNCells(t_AMR)
+    Ncells=Cmpi.allreduce(Ncells, op=Cmpi.SUM)
+    if Cmpi.master: print("[MESH GEN.] Number of Cells::%ge06"%(Ncells/1e06), flush=True)
+
+    ## ==================
+    ## ==== IBM Prep ====
+    ## ==================
     t_AMR = prepareAMRData(tb, t_AMR, IBM_parameters=IBM_parameters, dim=dim, check=check, localDir=localDir, forceAlignment=forceAlignment)
+    ## Ncells output
+    Ncells=C.getNCells(t_AMR)
+    Ncells=Cmpi.allreduce(Ncells, op=Cmpi.SUM)
+    if Cmpi.master: print("[IBM PREP.] Number of Cells::%ge06"%(Ncells/1e06), flush=True)
 
     if fileName is not None:
         Cmpi.convertPyTree2File(t_AMR, localDir+fileName)

Cassiopee/Generator/Generator/AMR.py

@@ -238,7 +238,8 @@ def generateListOfOffsets__(tb, snears, offsetValues=[], dim=3, opt=False, numTb
         bodies = [tbLocalTmp]; nbodies = 1
         BM = numpy.ones((1, nbodies), dtype=numpy.int32)
         t = C.newPyTree(["BASE", Internal.getZones(b)])
-        if dim == 2 or blankCellsAlgo == 'xray': X._blankCells(t, bodies, BM, blankingType='node_in', dim=dim, XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2)
+        # using the input bodies provided (tb & tbox) - safe to use blankCellsTri
+        if dim == 2: X._blankCells(t, bodies, BM, blankingType='node_in', dim=dim, XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2) #or blankCellsAlgo == 'xray'
         else: X._blankCellsTri(t, bodies, BM, blankingType='node_in')
         C._initVars(t, '{TurbulentDistance}={TurbulentDistance}*({cellN}>0.)-{TurbulentDistance}*({cellN}<1.)')
         #Cmpi.convertPyTree2File(b, 'meshForOffsetBase%d.cgns'%nBase) # DEBUG ONLY
@@ -395,6 +396,7 @@ def tagOutsideBody__(o, tbTMP, dim=3, h_target=-1., opt=False, noffsets=None, co
     # XRAYDIM1 = max(500, min(5000, XRAYDIM1)); # XRAYDIM1 = max(5000, min(50000, XRAYDIM1)); is too expensive need to find another solution
     C._initVars(to, "cellNIn", 1.)
 
+    # tbTMP - offset body - blankCells is more tested due to uncertainties on the quality of the offset
     if dim == 2 or blankCellsAlgo == 'xray': to = X.blankCells(to, bodies1, BM, blankingType='node_in',
                                                                XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM1, dim=dim,
                                                                cellNName='cellN')
@@ -424,9 +426,12 @@ def tagInsideOffset__(o, offset1=None, offset2=None, dim=3, h_target=-1., opt=Fa
     if dim==2:
         offset1 = T.addkplane(offset1)
         offset2 = T.addkplane(offset2)
-    bodies1 = [Internal.getZones(offset1)] # tag for refinement outside of offset1
+    bodies1 = [Internal.getZones(offset1)] # tag for refinement outside of offset1 (cgns base of tb or tbox)
     bodies2 = [Internal.getZones(offset2)] # tag for refinement inside of offset2
 
+    # bodies1 - tb & tbox - blankCellsTri is safe due to assure quality of the input
+    # bodies2 - offset body - blankCells is more tested due to uncertainties on the quality of the offset
+
     BM = numpy.ones((1,1), dtype=Internal.E_NpyInt)
 
     bb1 = G.bbox(offset1)
@@ -466,9 +471,9 @@ def tagInsideOffset__(o, offset1=None, offset2=None, dim=3, h_target=-1., opt=Fa
 
     if isTbox: C._initVars(to, '{cellN}=({cellNIn}<1)')
     else:
-        if dim == 2 or blankCellsAlgo == 'xray': to = X.blankCells(to, bodies1, BM, blankingType='node_in',
-                                                                   XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2, dim=dim,
-                                                                   cellNName='cellNOut')
+        if dim == 2 : to = X.blankCells(to, bodies1, BM, blankingType='node_in',
+                                        XRaydim1=XRAYDIM1, XRaydim2=XRAYDIM2, dim=dim,
+                                        cellNName='cellNOut') #or blankCellsAlgo == 'xray'
         else: to = X.blankCellsTri(to, bodies1, BM, blankingType='node_in', cellNName='cellNOut')
 
         # artificial shift the location of the boundary by depthLocal cells inwards (inside the body)