diff --git a/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py b/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py index 008d773ff..76b6f079c 100644 --- a/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py +++ b/Cassiopee/Apps/test/FastIBMWireMeshModel_t1.py @@ -45,7 +45,6 @@ t,tc = X_IBM.prepareIBMData(tb , None , None , tbox=tboffset, snears=snears , dfars=dfars , vmin=vmin, check=False , frontType=1) - test.testT(t , 1) test.testT(tc, 2) diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 2730715f1..72205fe06 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -224,27 +224,23 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N #=================== # STEP 0 : GET FILAMENT BODIES #=================== - [filamentBases, isFilamentOnly, isOrthoProjectFirst, tb, tbFilament]=D_IBM.determineClosedSolidFilament__(tb) - isWireModel=False - for z in Internal.getZones(t_case): - soldef = Internal.getNodeFromName(z,'.Solver#define') - if soldef is not None: - ibctype = Internal.getNodeFromName(soldef,'ibctype') - if ibctype is not None: - if Internal.getValue(ibctype) == "wiremodel": - isWireModel=True - break - if isFilamentOnly: tbLocal=tbFilament - else: tbLocal=Internal.merge([tb,tbFilament]) + tb, tbFilament = D_IBM.determineClosedSolidFilament__(tb) + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + #=================== # STEP 1 : GENERATE MESH #=================== + if isFilamentOnly: tbLocal=tbFilament + elif tbFilament: tbLocal=Internal.merge([tb,tbFilament]) + else: tbLocal=tb + + if t_in is None: if verbose: pt0 = python_time.time(); printTimeAndMemory__('generate Cartesian mesh', time=-1) test.printMem("Info: prepareIBMData: generate Cartesian mesh [start]") t = G_IBM.generateIBMMesh(tbLocal, vmin=vmin, snears=snears, dimPb=dimPb, dfars=dfars, tbox=tbox, - snearsf=snearsf, check=check, to=to, ext=depth+1, - expand=expand, dfarDir=dfarDir, octreeMode=octreeMode) + snearsf=snearsf, check=check, to=to, ext=depth+1, + expand=expand, dfarDir=dfarDir, octreeMode=octreeMode) Internal._rmNodesFromName(tb,"SYM") if balancing and Cmpi.size > 1: _redispatch__(t=t) @@ -263,7 +259,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute wall distance', time=-1) _dist2wallIBM(t, tb, dimPb=dimPb, frontType=frontType, Reynolds=Reynolds, yplus=yplus, Lref=Lref, correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42, - filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) + tbFilament=tbFilament) if verbose: printTimeAndMemory__('compute wall distance', time=python_time.time()-pt0) #=================== @@ -274,8 +270,10 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, - filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament, - isWireModel=isWireModel) + tbFilament=tbFilament) + ##Keep for now to check + #filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament, + #isWireModel=isWireModel) Cmpi.barrier() _redispatch__(t=t) if verbose: printTimeAndMemory__('blank by IBC bodies', time=python_time.time()-pt0) @@ -295,9 +293,9 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N # STEP 4 : BUILD FRONT #=================== if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1) - t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, + t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, tb=tb, dimPb=dimPb, frontType=frontType, cartesian=cartesian, twoFronts=twoFronts, check=check, - isWireModel=isWireModel) + tbFilament=tbFilament) if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0) #=================== @@ -307,8 +305,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N _setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, cartesian=cartesian, twoFronts=twoFronts, check=check, - isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst, isFilamentOnly=isFilamentOnly, - frontWMM=frontWMM) + tbFilament=tbFilament, frontWMM=frontWMM) if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0) #=================== @@ -317,7 +314,7 @@ def prepareIBMData(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tinit=N if verbose: pt0 = python_time.time(); printTimeAndMemory__('initialize and clean', time=-1) t, tc, tc2 = initializeIBM(t, tc, tb, tinit=tinit, tbCurvi=tbCurvi, dimPb=dimPb, twoFronts=twoFronts, - isWireModel=isWireModel, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly) + tbFilament=tbFilament) if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts) @@ -393,13 +390,12 @@ def _redispatch__(t=None, tc=None, tc2=None, twoFronts=False): # OUT: (optional) centers:TurbulentDistance_body%i fields #========================================================================= def dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1., - correctionMultiCorpsF42=False, heightMaxF42=-1., filamentBases=[], isFilamentOnly=False, - tbFilament=None): + correctionMultiCorpsF42=False, heightMaxF42=-1., tbFilament=None): """Compute the wall distance for IBM pre-processing.""" tp = Internal.copyRef(t) _dist2wallIBM(t, tb, dimPb=dimPb, frontType=frontType, Reynolds=Reynolds, yplus=yplus, Lref=Lref, correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42, - filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) + tbFilament=tbFilament) return tp def _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb): @@ -435,9 +431,11 @@ def _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb): return None def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1., - correctionMultiCorpsF42=False, heightMaxF42=-1., - filamentBases=[], isFilamentOnly=False, tbFilament=None): + correctionMultiCorpsF42=False, heightMaxF42=-1., tbFilament=None): """Compute the wall distance for IBM pre-processing.""" + + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + if dimPb == 2: # Set CoordinateZ to a fixed value z0 = Internal.getNodeFromType2(t, "Zone_t") @@ -452,7 +450,7 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. tbsave = tb2 - if filamentBases and not isFilamentOnly: + if tbFilament and not isFilamentOnly: C._initVars(t,'{centers:TurbulentDistanceSolid}={centers:TurbulentDistance}') if dimPb ==2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz) _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb) @@ -548,7 +546,9 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. #========================================================================= def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., - filamentBases=[], isFilamentOnly=False, tbFilament=None): + tbFilament=None): + + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) snear_min = 10e10 for z in Internal.getZones(tb): @@ -568,7 +568,7 @@ def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e C._initVars(t,'{centers:cellN}={centers:cellNIBC}') - if filamentBases: + if tbFilament: C._initVars(t,'{centers:cellNFil}={centers:cellN}') C._initVars(t,'{centers:cellNFilWMM}={centers:cellN}') @@ -687,7 +687,7 @@ def findIsoFront(cellNFront, Dist_1, Dist_2): if wallAdaptF42 is None: X._maximizeBlankedCells(t, depth=2, addGC=False) else: print("Info: blankingIBM: blankingF42 cannot operate with a local modeling height") - if filamentBases: + if tbFilament: tbFilamentWMM = [] for z in Internal.getZones(tbFilament): soldef = Internal.getNodeFromName(z,'.Solver#define') @@ -754,7 +754,7 @@ def findIsoFront(cellNFront, Dist_1, Dist_2): cellNFil[i,j]=2 cellN[i,j]=2 C._rmVars(t,['centers:TurbulentDistanceFilament','centers:TurbulentDistanceFilamentWMM']) - if filamentBases: C._rmVars(t,['centers:TurbulentDistanceSolid']) + if tbFilament: C._rmVars(t,['centers:TurbulentDistanceSolid']) if isFilamentOnly: C._initVars(t,'{centers:cellNFilWMM}={centers:cellNFil}') if not isFilamentOnly: _removeBlankedGrids(t, loc='centers') @@ -762,20 +762,23 @@ def findIsoFront(cellNFront, Dist_1, Dist_2): def blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., twoFronts=False, correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., - filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False): + tbFilament=None): """Blank the computational tree by IBC bodies for IBM pre-processing.""" tp = Internal.copyRef(t) _blankingIBM(t, tb, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, - filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament, isWireModel=isWireModel) + tbFilament=tbFilament) return tp def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., twoFronts=False, correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., - filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False): + tbFilament=None): """Blank the computational tree by IBC bodies for IBM pre-processing.""" + + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + if dimPb == 2: T._addkplane(tb) T._contract(tb, (0,0,0), (1,0,0), (0,1,0), 0.01) @@ -791,7 +794,7 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, Reynolds=Reynolds, yplus=yplus, Lref=Lref, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, - filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) + tbFilament=tbFilament) C._initVars(t, '{centers:cellNIBC}={centers:cellN}') @@ -945,9 +948,12 @@ def _pushBackImageFront2__(t, tc, tbbc, cartesian=False): return None -def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False, - isWireModel=False): +def buildFrontIBM(t, tc, tb=None, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False, + tbFilament=None): """Build the IBM front for IBM pre-processing.""" + + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + tbbc = Cmpi.createBBoxTree(tc) interpDataType = 0 if cartesian else 1 @@ -1023,20 +1029,22 @@ def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, #========================================================================= def setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False, - isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, frontWMM=None): + tbFilament=None, frontWMM=None): """Compute the transfer coefficients and data for IBM pre-processing.""" tp = Internal.copyRef(t) _setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, cartesian=cartesian, twoFronts=twoFronts, check=check, - isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst, - frontWMM=frontWMM) + tbFilament=tbFilament, frontWMM=frontWMM) return tp def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False, - isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, isFilamentOnly=False, frontWMM=None): + tbFilament=None, frontWMM=None): """Compute the transfer coefficients and data for IBM pre-processing.""" + + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + tbbc = Cmpi.createBBoxTree(tc) interpDataType = 0 if cartesian else 1 @@ -1057,7 +1065,7 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy else: tb_local= Internal.merge([tb,tbFilament]) res = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_local, tfront=front, frontType=frontType, cellNName='cellNIBC', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref, - isOrthoFirst=isOrthoProjectFirst) + isOrthoFirst=isFilamentOnly) if twoFronts: res2 = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb, tfront=front2, frontType=frontType, cellNName='cellNIBC', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref) @@ -1070,11 +1078,11 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy Internal._rmNode(tb_filament_localWMM,zlocal) res2 = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_filament_localWMM, tfront=frontWMM, frontType=frontType, cellNName='cellNFilWMM', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref, - isWireModel=isWireModel, isOrthoFirst=isOrthoProjectFirst) + isWireModel=isWireModel, isOrthoFirst=isFilamentOnly) restmp = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb_filament_localWMM, tfront=frontWMM, frontType=frontType, cellNName='cellNFilWMM', depth=depth, IBCType=IBCType, Reynolds=Reynolds, - yplus=yplus, Lref=Lref, isOrthoFirst=isOrthoProjectFirst) + yplus=yplus, Lref=Lref, isOrthoFirst=isFilamentOnly) for j in range(3): ##delete in res @@ -1325,8 +1333,10 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy # OUT: updated t, tc # OUT: (optional) new connectivity tree tc2 #========================================================================= -def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None, filamentBases=[], isFilamentOnly=False): +def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None): + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + filteredBC=['outpress', 'inj', 'slip', 'overlap'] recompute=False tb2 = Internal.copyRef(tb) @@ -1351,7 +1361,7 @@ def _recomputeDistForViscousWall__(t, tb, tbCurvi=None, dimPb=3, tbFilament=None tbsave = tb2 - if filamentBases and not isFilamentOnly: + if tbFilament and not isFilamentOnly: if dimPb == 2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz) tbFilamentnoWMM = [] tbFilamentWMM = [] @@ -1440,10 +1450,11 @@ def _tInitialize__(t, tinit=None, model='NSTurbulent', isWireModel=False): C._initVars(z,'{centers:'+v_local+'_WM}=0.') return None -def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False, isWireModel=False, - tbFilament=None, filamentBases=[], isFilamentOnly=False): +def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False, tbFilament=None): """Initialize the computational and connectivity trees for IBM pre-processing.""" + isFilamentOnly, isWireModel = D_IBM.localWMMFlags__(tb, tbFilament) + model = Internal.getNodeFromName(tb, 'GoverningEquations') if model is None: raise ValueError('initializeIBM: GoverningEquations is missing in input geometry tree.') model = Internal.getValue(model) @@ -1453,7 +1464,7 @@ def initializeIBM(t, tc, tb, tinit=None, tbCurvi=None, dimPb=3, twoFronts=False, ibctypes = list(set(Internal.getValue(ibc) for ibc in ibctypes)) if model != 'Euler': - _recomputeDistForViscousWall__(t, tb, tbCurvi=tbCurvi, dimPb=dimPb, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly ) + _recomputeDistForViscousWall__(t, tb, tbCurvi=tbCurvi, dimPb=dimPb, tbFilament=tbFilament) tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None diff --git a/Cassiopee/Geom/Geom/IBM.py b/Cassiopee/Geom/Geom/IBM.py index c9b948bc0..95be68b36 100644 --- a/Cassiopee/Geom/Geom/IBM.py +++ b/Cassiopee/Geom/Geom/IBM.py @@ -542,21 +542,46 @@ def determineClosedSolidFilament__(tb): filamentBases = [] isFilamentOnly= False - len_tb = len(Internal.getBases(tb)) for b in Internal.getBases(tb): if "IBCFil" in b[0]:filamentBases.append(b[0]) - if len(filamentBases) == len_tb:isFilamentOnly=True + if len(filamentBases) == len(Internal.getBases(tb)):isFilamentOnly=True isOrthoProjectFirst = isFilamentOnly - ## if tb has both a closed solid and filaments + ## assume isFilamentOnly=True tbFilament = Internal.copyTree(tb) + + ## if tb has 1) no filaments or 2) filaments & closed bodies if not isFilamentOnly: - tbFilament = [] - for b in filamentBases: - node_local = Internal.getNodeFromNameAndType(tb, b, 'CGNSBase_t') - tbFilament.append(node_local) - Internal._rmNode(tb,node_local) - isOrthoProjectFirst = True - tbFilament = C.newPyTree(tbFilament); - return [filamentBases, isFilamentOnly, isOrthoProjectFirst, tb, tbFilament] + if len(filamentBases)==0: + tbFilament=None + else: + ##tb : only closed bodies + ##tbFilament: filament bodies + tbFilament = [] + for b in filamentBases: + node_local = Internal.getNodeFromNameAndType(tb, b, 'CGNSBase_t') + tbFilament.append(node_local) + Internal._rmNode(tb,node_local) + isOrthoProjectFirst = True + tbFilament = C.newPyTree(tbFilament); + + return tb, tbFilament + + +def localWMMFlags__(tb,tbFilament): + isFilamentOnly=False + isWireModel =False + + if tbFilament: + if len(Internal.getBases(tbFilament))==len(Internal.getBases(tb)): + isFilamentOnly=True + for z in Internal.getZones(tbFilament): + soldef = Internal.getNodeFromName(z,'.Solver#define') + if soldef is not None: + ibctype = Internal.getNodeFromName(soldef,'ibctype') + if ibctype is not None: + if Internal.getValue(ibctype) == "wiremodel": + isWireModel=True + break + return isFilamentOnly,isWireModel