diff --git a/Cassiopee/Apps/Apps/Fast/IBM.py b/Cassiopee/Apps/Apps/Fast/IBM.py index 5354b2854..be44b6ea8 100644 --- a/Cassiopee/Apps/Apps/Fast/IBM.py +++ b/Cassiopee/Apps/Apps/Fast/IBM.py @@ -43,6 +43,13 @@ __IBCNameServer__={} + +## ============================= June 2024 ================================= +## No further developments in this file. All IBM developments need to be +## performed in the IBM.py files in each Cassiopee module (e.g.$CASSIOPEE/Cassiopee/Connector/Connector/IBM.py) +## Any developments in this file will be REJECTED to be merged (pull request) in the official/upstream ONERA Cassiopee. +## Please address any questions to: Christophe Benoit, Benjamin Constant, Antoine Jost, and Stephanie Peron. + def getIBCDName(proposedName): global __IBCNameServer__ (ibcname,__IBCNameServer__) = C.getUniqueName(proposedName, __IBCNameServer__) diff --git a/Cassiopee/Connector/Connector/IBM.py b/Cassiopee/Connector/Connector/IBM.py index 1c14e172b..380c50d28 100644 --- a/Cassiopee/Connector/Connector/IBM.py +++ b/Cassiopee/Connector/Connector/IBM.py @@ -12,6 +12,7 @@ import Initiator.PyTree as I import Converter.Distributed as Distributed import Generator.IBMmodelHeight as G_IBM_Height +import Geom.IBM as D_IBM import Transform.PyTree as T import Converter.Internal as Internal import Compressor.PyTree as Compressor @@ -184,7 +185,7 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin IBCType=1, verbose=True, check=False, balancing=False, distribute=False, twoFronts=False, cartesian=False, yplus=100., Lref=1., correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., - tbOneOver=None): + tbOneOver=None, generateCartesianMeshOnly=False): import Generator.IBM as G_IBM import time as python_time @@ -192,6 +193,9 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin if isinstance(t_case, str): tb = C.convertFile2PyTree(t_case) else: tb = Internal.copyTree(t_case) + if isinstance(tc_out, str): fileoutpre = tc_out.split('/') + else: fileoutpre = ['.','template.cgns'] + refstate = Internal.getNodeFromName(tb, 'ReferenceState') flowEqn = Internal.getNodeFromName(tb, 'FlowEquationSet') @@ -220,6 +224,21 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin if any(ibc in ['Musker', 'MuskerMob', 'Mafzal', 'Log', 'TBLE', 'TBLE_FULL'] for ibc in ibctypes): raise ValueError("prepareIBMDataPara: governing equations (Euler) not consistent with ibc types %s"%(ibctypes)) + #=================== + # 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]) #=================== # STEP 1 : GENERATE MESH #=================== @@ -227,14 +246,20 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin if t_in is None: if verbose: pt0 = python_time.time(); printTimeAndMemory__('generate Cartesian mesh', time=-1) test.printMem("Info: prepareIBMDataPara: generate Cartesian mesh [start]") - t = G_IBM.generateIBMMeshPara(tb, vmin=vmin, snears=snears, dimPb=dimPb, dfar=dfar, dfarList=dfarList, tbox=tbox, + t = G_IBM.generateIBMMeshPara(tbLocal, vmin=vmin, snears=snears, dimPb=dimPb, dfar=dfar, dfarList=dfarList, tbox=tbox, snearsf=snearsf, check=check, to=to, ext=depth+1, expand=expand, dfarDir=dfarDir, check_snear=False, mode=mode, - tbOneOver=tbOneOver, listF1save=listF1save) + tbOneOver=tbOneOver, listF1save=listF1save, fileoutpre=fileoutpre) Internal._rmNodesFromName(tb,"SYM") if balancing and Cmpi.size > 1: _redispatch__(t=t) if verbose: printTimeAndMemory__('generate Cartesian mesh', time=python_time.time()-pt0) + if generateCartesianMeshOnly: + fileoutpre[-1]='CartMesh.cgns' + fileout = '/'.join(fileoutpre) + Cmpi.convertPyTree2File(t,fileout) + if Cmpi.rank == 0: print("Wrote CartMesh.cgns...Exiting Prep...Bye") + exit() else: t = t_in @@ -248,17 +273,20 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin #=================== 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) + correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42, + filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) if verbose: printTimeAndMemory__('compute wall distance', time=python_time.time()-pt0) - + #=================== # STEP 3 : BLANKING IBM #=================== if verbose: pt0 = python_time.time(); printTimeAndMemory__('blank by IBC bodies', time=-1) _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, listF1save=listF1save) + Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, + heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save, + 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) @@ -278,8 +306,9 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin # STEP 4 : BUILD FRONT #=================== if verbose: pt0 = python_time.time(); printTimeAndMemory__('build IBM front', time=-1) - t, tc, front, front2 = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, - cartesian=cartesian, twoFronts=twoFronts, check=check) + t, tc, front, front2, frontWMM = buildFrontIBM(t, tc, dimPb=dimPb, frontType=frontType, + cartesian=cartesian, twoFronts=twoFronts, check=check, + isWireModel=isWireModel, fileoutpre=fileoutpre) if verbose: printTimeAndMemory__('build IBM front', time=python_time.time()-pt0) #=================== @@ -287,15 +316,18 @@ def prepareIBMDataPara(t_case, t_out, tc_out, t_in=None, to=None, tbox=None, tin #=================== if verbose: pt0 = python_time.time(); printTimeAndMemory__('compute interpolation data (IBM)', time=-1) _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) + Reynolds=Reynolds, yplus=yplus, Lref=Lref, + cartesian=cartesian, twoFronts=twoFronts, check=check, + isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst, isFilamentOnly=isFilamentOnly, + frontWMM=frontWMM, fileoutpre=fileoutpre) if verbose: printTimeAndMemory__('compute interpolation data (IBM)', time=python_time.time()-pt0) #=================== # STEP 6 : INIT IBM #=================== if verbose: pt0 = python_time.time(); printTimeAndMemory__('initialize and clean', time=-1) - t, tc, tc2 = initializeIBM(t, tc, tb, tinit=tinit, dimPb=dimPb, twoFronts=twoFronts) + t, tc, tc2 = initializeIBM(t, tc, tb, tinit=tinit, dimPb=dimPb, twoFronts=twoFronts, + isWireModel=isWireModel, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly) if distribute and Cmpi.size > 1: _redispatch__(t=t, tc=tc, tc2=tc2, twoFronts=twoFronts) @@ -367,15 +399,50 @@ 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.): + correctionMultiCorpsF42=False, heightMaxF42=-1., filamentBases=[], isFilamentOnly=False, + 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) + correctionMultiCorpsF42=correctionMultiCorpsF42, heightMaxF42=heightMaxF42, + filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) return tp +def _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb): + tbFilamentnoWMM = [] + tbFilamentWMM = [] + for z in Internal.getZones(tb2): + soldef = Internal.getNodeFromName(z,'.Solver#define') + ibctype = Internal.getNodeFromName(soldef,'ibctype') + if Internal.getValue(ibctype) != "wiremodel": + tbFilamentnoWMM.append(z) + else: + tbFilamentWMM.append(z) + + if tbFilamentnoWMM: + tbFilamentnoWMM = C.newPyTree(['BasenoWMM', tbFilamentnoWMM]) + tbsave = Internal.merge([tbsave,tbFilamentnoWMM]) + C._initVars(t,'{centers:TurbulentDistance}=1e06') + DTW._distance2Walls(t, tbFilamentnoWMM, type='ortho', signed=0, dim=dimPb, loc='centers') + C._initVars(t,'{centers:TurbulentDistanceFilament}={centers:TurbulentDistance}') + C._initVars(t,'{centers:TurbulentDistance}=minimum({centers:TurbulentDistanceSolid},{centers:TurbulentDistanceFilament})') + + if tbFilamentWMM: + C._initVars(t,'{centers:TurbulentDistance}=1e06') + DTW._distance2Walls(t, tbFilamentWMM, type='ortho', signed=0, dim=dimPb, loc='centers') + C._initVars(t,'{centers:TurbulentDistanceFilamentWMM}={centers:TurbulentDistance}') + C._initVars(t,'{centers:TurbulentDistance}=minimum({centers:TurbulentDistanceSolid},{centers:TurbulentDistanceFilamentWMM})') + if tbFilamentnoWMM: + C._initVars(t,'{centers:TurbulentDistance}=minimum({centers:TurbulentDistance},{centers:TurbulentDistanceFilament})') + + C._initVars(t,"{centers:TurbulentDistanceSolid}=({centers:TurbulentDistanceSolid}>1e03)*0+({centers:TurbulentDistanceSolid}<1e03)*{centers:TurbulentDistanceSolid}") + C._initVars(t,"{centers:TurbulentDistanceFilament}=({centers:TurbulentDistanceFilament}>1e03)*0+({centers:TurbulentDistanceFilament}<1e03)*{centers:TurbulentDistanceFilament}") + C._initVars(t,"{centers:TurbulentDistanceFilamentWMM}=({centers:TurbulentDistanceFilamentWMM}>1e03)*0+({centers:TurbulentDistanceFilamentWMM}<1e03)*{centers:TurbulentDistanceFilamentWMM}") + return None + def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1., - correctionMultiCorpsF42=False, heightMaxF42=-1.): + correctionMultiCorpsF42=False, heightMaxF42=-1., + filamentBases=[], isFilamentOnly=False, tbFilament=None): """Compute the wall distance for IBM pre-processing.""" if dimPb == 2: # Set CoordinateZ to a fixed value @@ -387,9 +454,20 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. else: tb2 = tb - # Compute distance to bodies DTW._distance2Walls(t, tb2, type='ortho', signed=0, dim=dimPb, loc='centers') + tbsave = tb2 + + if filamentBases and not isFilamentOnly: + C._initVars(t,'{centers:TurbulentDistanceSolid}={centers:TurbulentDistance}') + if dimPb ==2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz) + _dist2wallIBMFilamentWMM__(t, tb2, tbsave, dimPb) + elif isFilamentOnly: + C._initVars(t,"{centers:TurbulentDistanceSolid}=0") + C._initVars(t,"{centers:TurbulentDistanceFilament}={centers:TurbulentDistance}") + C._initVars(t,"{centers:TurbulentDistanceFilamentWMM}={centers:TurbulentDistance}") + + # Compute turbulentdistance wrt each body that is not a sym plan (centers:TurbulentDistance_bodyX) if correctionMultiCorpsF42 and frontType == 42: print("correctionMultiCorpsF42") @@ -467,12 +545,16 @@ def _dist2wallIBM(t, tb, dimPb=3, frontType=1, Reynolds=1.e6, yplus=100, Lref=1. # IN: wallAdaptF42 (cloud of IB target points with yplus information): # use a previous computation to adapt the positioning of IB target points around the geometry according to a target yplus (F42) # IN: heightMaxF42 (float): maximum modeling height for the location of IB target points around the geometry (F42) +# IN: listF1save: list of zones that will have an F1 front +# IN: filamentBases: list of bases that are IBC filaments +# IN: isFilamentOnly: boolean on whether there is only a IBC filament +# IN: tbFilament: PyTree of the IBC filaments # OUT: centers:cellN, centers:cellNIBC, centers:cellNChim, centers:cellNFront fields # OUT: (optional) centers:cellNIBC_2, centers:cellNFront_2 fields for second image points #========================================================================= 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., - listF1save=[]): + listF1save=[], filamentBases=[], isFilamentOnly=False, tbFilament=None): snear_min = 10e10 for z in Internal.getZones(tb): @@ -484,13 +566,17 @@ def _blankingIBM__(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e if snearl is not None: snear_min = min(snear_min,snearl) snear_min = Cmpi.allreduce(snear_min, op=Cmpi.MIN) - _blankByIBCBodies(t, tb, 'centers', dimPb) + if not isFilamentOnly: _blankByIBCBodies(t, tb, 'centers', dimPb) C._initVars(t, '{centers:cellNIBC}={centers:cellN}') - _signDistance(t) + if not isFilamentOnly: _signDistance(t) C._initVars(t,'{centers:cellN}={centers:cellNIBC}') + + if filamentBases: + C._initVars(t,'{centers:cellNFil}={centers:cellN}') + C._initVars(t,'{centers:cellNFilWMM}={centers:cellN}') # determination des pts IBC if frontType != 42: @@ -605,26 +691,101 @@ 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") - _removeBlankedGrids(t, loc='centers') + if filamentBases: + tbFilamentWMM = [] + for z in Internal.getZones(tbFilament): + soldef = Internal.getNodeFromName(z,'.Solver#define') + ibctype = Internal.getNodeFromName(soldef,'ibctype') + if Internal.getValue(ibctype) == "wiremodel": + tbFilamentWMM.append(z) + maxy=C.getMaxValue(tbFilamentWMM, ['CoordinateY']); + miny=C.getMinValue(tbFilamentWMM, ['CoordinateY']); + + if dimPb == 3: + maxz=C.getMaxValue(tbFilamentWMM, ['CoordinateZ']); + minz=C.getMinValue(tbFilamentWMM, ['CoordinateZ']); + + for z in Internal.getZones(t): + if C.getMaxValue(z, 'centers:TurbulentDistanceFilament')> 1e-05: + sol = Internal.getNodeByName(z,"FlowSolution#Centers") + cellNFil= Internal.getNodeByName(sol,'cellNFil')[1] + cellN = Internal.getNodeByName(sol,'cellN')[1] + dist = Internal.getNodeByName(sol,'TurbulentDistanceFilament')[1] + ycord = Internal.getNodeByName(z,'CoordinateY')[1] + h = abs(C.getValue(z,'CoordinateX',4)-C.getValue(z,'CoordinateX',5)) + sh = numpy.shape(dist) + if dimPb == 3: + for k in range(sh[2]): + for j in range(sh[1]): + for i in range(sh[0]): + valy=0.5*(ycord[i,j,k]+ycord[i,j+1,k]) + if dist[i,j,k] 1e-05: + sol = Internal.getNodeByName(z,"FlowSolution#Centers") + cellNFil= Internal.getNodeByName(sol,'cellNFilWMM')[1] + cellN = Internal.getNodeByName(sol,'cellN')[1] + dist = Internal.getNodeByName(sol,'TurbulentDistanceFilamentWMM')[1] + ycord = Internal.getNodeByName(z,'CoordinateY')[1] + h = abs(C.getValue(z,'CoordinateX',4)-C.getValue(z,'CoordinateX',5)) + sh = numpy.shape(dist) + if dimPb == 3: + zcord = Internal.getNodeByName(z,'CoordinateZ')[1] + for k in range(sh[2]): + for j in range(sh[1]): + for i in range(sh[0]): + valy=0.5*(ycord[i,j,k]+ycord[i,j+1,k ]) + valz=0.5*(zcord[i,j,k]+zcord[i,j ,k+1]) + if dist[i,j,k]miny and \ + valzminz : + cellNFil[i,j,k]=2 + cellN[i,j,k]=2 + else: + for j in range(sh[1]): + for i in range(sh[0]): + valy=0.5*(ycord[i,j,0]+ycord[i,j+1,0]) + if dist[i,j,0]miny and cellN[i,j]==1: + cellNFil[i,j]=2 + cellN[i,j]=2 + C._rmVars(t,['centers:TurbulentDistanceFilament','centers:TurbulentDistanceFilamentWMM']) + if filamentBases: C._rmVars(t,['centers:TurbulentDistanceSolid']) + if isFilamentOnly: C._initVars(t,'{centers:cellNFilWMM}={centers:cellNFil}') + + if not isFilamentOnly: _removeBlankedGrids(t, loc='centers') return None 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., listF1save=[]): + correctionMultiCorpsF42=False, blankingF42=False, wallAdaptF42=None, heightMaxF42=-1., listF1save=[], + filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False): """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, listF1save=listF1save) + Reynolds=Reynolds, yplus=yplus, Lref=Lref, twoFronts=twoFronts, + heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save, + filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament, isWireModel=isWireModel) 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., - listF1save=[]): + listF1save=[], filamentBases=[], isFilamentOnly=False, tbFilament=None, isWireModel=False): """Blank the computational tree by IBC bodies for IBM pre-processing.""" if dimPb == 2: T._addkplane(tb) T._contract(tb, (0,0,0), (1,0,0), (0,1,0), 0.01) + if tbFilament: + T._addkplane(tbFilament) + T._contract(tbFilament, (0,0,0), (1,0,0), (0,1,0), 0.01) X._applyBCOverlaps(t, depth=depth, loc='centers', val=2, cellNName='cellN') C._initVars(t,'{centers:cellNChim}={centers:cellN}') @@ -633,7 +794,8 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, _blankingIBM__(t, tb, dimPb=dimPb, frontType=frontType, IBCType=IBCType, depth=depth, Reynolds=Reynolds, yplus=yplus, Lref=Lref, heightMaxF42=heightMaxF42, correctionMultiCorpsF42=correctionMultiCorpsF42, - wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save) + wallAdaptF42=wallAdaptF42, blankingF42=blankingF42, listF1save=listF1save, + filamentBases=filamentBases, isFilamentOnly=isFilamentOnly, tbFilament=tbFilament) C._initVars(t, '{centers:cellNIBC}={centers:cellN}') @@ -649,7 +811,11 @@ def _blankingIBM(t, tb, dimPb=3, frontType=1, IBCType=1, depth=2, Reynolds=1.e6, Internal.__FlowSolutionCenters__) else: C._initVars(t,'{centers:cellNFront}=logical_and({centers:cellNIBC}>0.5, {centers:cellNIBC}<1.5)') - + + if isWireModel: + C._initVars(t,'{centers:cellNFrontFilWMM}={centers:cellNFilWMM}*({centers:cellNFilWMM}>0.5)+1*({centers:cellNFilWMM}<0.5)') + C._initVars(t,'{centers:cellNFrontFilWMM}=logical_and({centers:cellNFrontFilWMM}>0.5, {centers:cellNFrontFilWMM}<1.5)') + for z in Internal.getZones(t): if twoFronts: epsilon_dist = abs(C.getValue(z,'CoordinateX',1)-C.getValue(z,'CoordinateX',0)) @@ -783,7 +949,8 @@ def _pushBackImageFront2__(t, tc, tbbc, cartesian=False): return None -def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False): +def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, check=False, + isWireModel=False, fileoutpre='./'): """Build the IBM front for IBM pre-processing.""" tbbc = Cmpi.createBBoxTree(tc) @@ -824,11 +991,26 @@ def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, else: front2 = None + if isWireModel: + frontWMM = getIBMFront(tc, 'cellNFrontFilWMM', dim=dimPb, frontType=frontType) + frontWMM = gatherFront(frontWMM) + else: + frontWMM = None + if check and Cmpi.rank == 0: - C.convertPyTree2File(front, 'front.cgns') - if twoFronts: C.convertPyTree2File(front2, 'front2.cgns') + fileoutpre[-1]='front.cgns' + fileout = '/'.join(fileoutpre) + C.convertPyTree2File(front, fileout) + if twoFronts: + fileoutpre[-1]='front2.cgns' + fileout = '/'.join(fileoutpre) + C.convertPyTree2File(front2, fileout) + if isWireModel: + fileoutpre[-1]='frontWMM.cgns' + fileout = '/'.join(fileoutpre) + C.convertPyTree2File(frontWMM, fileout) - return t, tc, front, front2 + return t, tc, front, front2, frontWMM #========================================================================= # Compute the transfer coefficients and data for IBM pre-processing. @@ -850,16 +1032,20 @@ def buildFrontIBM(t, tc, dimPb=3, frontType=1, cartesian=False, twoFronts=False, # OUT: (optional) 2_IBCD* zones inside tc #========================================================================= 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): + yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False, + isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, frontWMM=None, fileoutpre='./'): """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) + Reynolds=Reynolds, yplus=yplus, Lref=Lref, + cartesian=cartesian, twoFronts=twoFronts, check=check, + isWireModel=isWireModel, tbFilament=tbFilament, isOrthoProjectFirst=isOrthoProjectFirst, + frontWMM=frontWMM, fileoutpre=fileoutpre) 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): + yplus=100, Lref=1., cartesian=False, twoFronts=False, check=False, + isWireModel=False, tbFilament=None, isOrthoProjectFirst=False, isFilamentOnly=False, frontWMM=None, fileoutpre='./'): """Compute the transfer coefficients and data for IBM pre-processing.""" tbbc = Cmpi.createBBoxTree(tc) @@ -873,14 +1059,47 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy nbZonesIBC = len(zonesRIBC) if nbZonesIBC == 0: res = [{},{},{}] - if twoFronts: res2 = [{},{},{}] + if twoFronts or isWireModel: res2 = [{},{},{}] else: - res = getAllIBMPoints(zonesRIBC, loc='centers',tb=tb, tfront=front, frontType=frontType, - cellNName='cellNIBC', depth=depth, IBCType=IBCType, Reynolds=Reynolds, yplus=yplus, Lref=Lref) + tb_local = tb + if tbFilament: + if isFilamentOnly:tb_local= tbFilament + 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) 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) - + if isWireModel: + tb_filament_localWMM=Internal.copyTree(tbFilament) + for z in Internal.getZones(tbFilament): + ibctype = Internal.getNodeFromName(Internal.getNodeFromName(z,'.Solver#define'),'ibctype') + if Internal.getValue(ibctype) != "wiremodel": + zlocal=Internal.getNodeFromName(tb_filament_localWMM,z[0]) + 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) + + 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) + + for j in range(3): + ##delete in res + item_del=[] + for ii in res[j]: + if "140" in ii: item_del.append(ii) + for ii in item_del: del res[j][ii] + ##detele in res2 + item_del=[] + for ii in res2[j]: + if "140" not in ii: item_del.append(ii) + for ii in item_del: del res2[j][ii] + ##add to res + for ii in restmp[j]: res[j][ii] = restmp[j][ii] + # cleaning C._rmVars(tc,['cellNChim','cellNIBC','TurbulentDistance','cellNFront']) # dans t, il faut cellNChim et cellNIBCDnr pour recalculer le cellN a la fin @@ -889,6 +1108,11 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy front = None if twoFronts: front2 = None + if isWireModel: + C._rmVars(tc,['cellNFrontFilWMM']) + C._rmVars(t,['cellNFrontFilWMM']) + frontWMM = None + graph = {}; datas = {} procDict = Cmpi.getProcDict(tc) @@ -900,7 +1124,7 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy dictOfInterpPtsByIBCType = res[2] interDictIBM={} - if twoFronts: + if twoFronts or isWireModel: dictOfCorrectedPtsByIBCType2 = res2[0] dictOfWallPtsByIBCType2 = res2[1] dictOfInterpPtsByIBCType2 = res2[2] @@ -917,6 +1141,13 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy allInterpPts = dictOfInterpPtsByIBCType[ibcTypeL] for nozr in range(nbZonesIBC): if allCorrectedPts[nozr] != []: + ##[AJ] Keep temporarily + ##if ibcTypeL=="140#filament": + ## nlen = numpy.shape(allInterpPts[nozr][1])[1] + ## save2file = numpy.zeros((nlen,2),dtype=float) + ## save2file[:,0]=allInterpPts[nozr][1][0][:] + ## save2file[:,1]=allInterpPts[nozr][1][1][:] + ## numpy.savetxt('allInterpPts.txt', save2file, delimiter=',') # X is an array zrname = zonesRIBC[nozr][0] interpPtsBB = Generator.BB(allInterpPts[nozr]) for z in zones: @@ -929,13 +1160,20 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy if zrname not in interDictIBM: interDictIBM[zrname]=[zname] else: if zname not in interDictIBM[zrname]: interDictIBM[zrname].append(zname) - if twoFronts: + if twoFronts or isWireModel: for ibcTypeL in dictOfCorrectedPtsByIBCType2: allCorrectedPts2 = dictOfCorrectedPtsByIBCType2[ibcTypeL] allWallPts2 = dictOfWallPtsByIBCType2[ibcTypeL] allInterpPts2 = dictOfInterpPtsByIBCType2[ibcTypeL] for nozr in range(nbZonesIBC): if allCorrectedPts2[nozr] != []: + ##[AJ] Keep temporarily + ##if ibcTypeL=="140#filament": + ## nlen = numpy.shape(allInterpPts2[nozr][1])[1] + ## save2file = numpy.zeros((nlen,2),dtype=float) + ## save2file[:,0]=allInterpPts2[nozr][1][0][:] + ## save2file[:,1]=allInterpPts2[nozr][1][1][:] + ## numpy.savetxt('allInterpPts2.txt', save2file, delimiter=',') # X is an array zrname = zonesRIBC[nozr][0] interpPtsBB2 = Generator.BB(allInterpPts2[nozr]) for z in zones: @@ -1052,7 +1290,7 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy dictOfWallPtsByIBCType = None dictOfInterpPtsByIBCType = None interDictIBM = None - if twoFronts: + if twoFronts or isWireModel: dictOfCorrectedPtsByIBCType2 = None dictOfWallPtsByIBCType2 = None dictOfInterpPtsByIBCType2 = None @@ -1078,8 +1316,14 @@ def _setInterpDataIBM(t, tc, tb, front, front2=None, dimPb=3, frontType=1, IBCTy C._rmVars(t, varsRM) if check: - extractIBMInfo(tc, IBCNames="IBCD_*", filename_out='IBMInfo.cgns') - if twoFronts: extractIBMInfo(tc, IBCNames="2_IBCD_*", filename_out='IBMInfo2.cgns') + fileoutpre[-1]='IBMInfo.cgns' + fileout = '/'.join(fileoutpre) + extractIBMInfo(tc, IBCNames="IBCD_*", filename_out=fileout) + + if twoFronts: + fileoutpre[-1]='IBMInfo2.cgns' + fileout = '/'.join(fileoutpre) + extractIBMInfo(tc, IBCNames="2_IBCD_*", filename_out=fileout) return None @@ -1094,7 +1338,7 @@ 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, dimPb=3): +def _recomputeDistForViscousWall__(t, tb, dimPb=3, tbFilament=None, filamentBases=[], isFilamentOnly=False): for z in Internal.getZones(tb): ibc = Internal.getNodeFromName(z, 'ibctype') @@ -1110,13 +1354,30 @@ def _recomputeDistForViscousWall__(t, tb, dimPb=3): else: tb2 = tb - DTW._distance2Walls(t, tb2, type='ortho', signed=0, dim=dimPb, loc='centers') + tbsave = tb2 + + if filamentBases and not isFilamentOnly: + if dimPb ==2: tb2 = C.initVars(tbFilament, 'CoordinateZ', dz) + tbFilamentnoWMM = [] + tbFilamentWMM = [] + for z in Internal.getZones(tb2): + soldef = Internal.getNodeFromName(z,'.Solver#define') + ibctype = Internal.getNodeFromName(soldef,'ibctype') + if Internal.getValue(ibctype) != "wiremodel": + tbFilamentnoWMM.append(z) + else: + tbFilamentWMM.append(z) + + if tbFilamentnoWMM: + tbFilamentnoWMM = C.newPyTree(['BasenoWMM', tbFilamentnoWMM]) + tbsave = Internal.merge([tbsave,tbFilamentnoWMM]) + + DTW._distance2Walls(t, tbsave, type='ortho', signed=0, dim=dimPb, loc='centers') return None -def _tcInitialize__(tc, tc2=None, twoFronts=False, ibctypes=[]): - import Geom.IBM as D_IBM - +def _tcInitialize__(tc, tc2=None, twoFronts=False, ibctypes=[], isWireModel=False): + if twoFronts: Internal._rmNodesByName(tc2, 'IBCD*') Internal._rmNodesByName(tc, '2_IBCD*') @@ -1158,6 +1419,12 @@ def _tcInitialize__(tc, tc2=None, twoFronts=False, ibctypes=[]): # if twoFronts: # Internal._createUniqueChild(solverIBC, 'isgradP' , 'DataArray_t', 'True') # Internal._createUniqueChild(solverIBC, 'alphaGrad' , 'DataArray_t', 0) + if isWireModel: + Internal._createUniqueChild(solverIBC, 'isWireModel' , 'DataArray_t', 'True') + Internal._createUniqueChild(solverIBC, 'DeltaVWire' , 'DataArray_t', 0) + Internal._createUniqueChild(solverIBC, 'KWire' , 'DataArray_t', 0) + Internal._createUniqueChild(solverIBC, 'DiameterWire', 'DataArray_t', 0) + Internal._createUniqueChild(solverIBC, 'CtWire' , 'DataArray_t', 0) if 'TBLE' in ibctypes or 'TBLE_FULL' in ibctypes: Internal._createUniqueChild(solverIBC, 'isTBLE' , 'DataArray_t', 'True') Internal._createUniqueChild(solverIBC, 'alphaGrad' , 'DataArray_t', 0) @@ -1165,14 +1432,21 @@ def _tcInitialize__(tc, tc2=None, twoFronts=False, ibctypes=[]): return None -def _tInitialize__(t, tinit=None, model='NSTurbulent'): +def _tInitialize__(t, tinit=None, model='NSTurbulent', isWireModel=False): if tinit is None: I._initConst(t, loc='centers') else: t = Pmpi.extractMesh(tinit, t, mode='accurate') if model != "Euler": C._initVars(t, 'centers:ViscosityEddy', 0.) + if isWireModel: + vars_wm = ['Density','VelocityX','VelocityY','VelocityZ','Temperature'] + if model == 'NSTurbulent':vars_wm.append('TurbulentSANuTilde') + for z in Internal.getZones(t): + for v_local in vars_wm: + C._initVars(z,'{centers:'+v_local+'_WM}=0.') return None -def initializeIBM(t, tc, tb, tinit=None, dimPb=3, twoFronts=False): +def initializeIBM(t, tc, tb, tinit=None, dimPb=3, twoFronts=False, isWireModel=False, + tbFilament=None, filamentBases=[], isFilamentOnly=False): """Initialize the computational and connectivity trees for IBM pre-processing.""" model = Internal.getNodeFromName(tb, 'GoverningEquations') @@ -1185,12 +1459,22 @@ def initializeIBM(t, tc, tb, tinit=None, dimPb=3, twoFronts=False): if model != 'Euler': if any(ibc in ['outpress', 'inj', 'slip'] for ibc in ibctypes): - _recomputeDistForViscousWall__(t, tb, dimPb=dimPb) + _recomputeDistForViscousWall__(t, tb, dimPb=dimPb, tbFilament=tbFilament, filamentBases=filamentBases, isFilamentOnly=isFilamentOnly) - tc2 = Internal.copyTree(tc) if twoFronts else None - - _tcInitialize__(tc, tc2=tc2, twoFronts=twoFronts, ibctypes=ibctypes) - _tInitialize__(t, tinit=tinit, model=model) + tc2 = Internal.copyTree(tc) if twoFronts or isWireModel else None + + if isWireModel: + Internal._rmNodesByName(tc2, 'IBCD*') + Internal._rmNodesByName(tc, '2_IBCD*') + D_IBM._transformTc2(tc2) + tc2 = Internal.rmNodesByName(tc2, 'ID*') + NewIBCD= 141 + _changeNameIBCD__(tc2,NewIBCD) + tc = Internal.merge([tc,tc2]) + tc2 = None + + _tcInitialize__(tc, tc2=tc2, twoFronts=twoFronts, ibctypes=ibctypes, isWireModel=isWireModel) + _tInitialize__(t, tinit=tinit, model=model, isWireModel=isWireModel) return t, tc, tc2 diff --git a/Cassiopee/Connector/test/buildFrontIBMPT_t1.py b/Cassiopee/Connector/test/buildFrontIBMPT_t1.py index 448661314..d589d2632 100755 --- a/Cassiopee/Connector/test/buildFrontIBMPT_t1.py +++ b/Cassiopee/Connector/test/buildFrontIBMPT_t1.py @@ -19,7 +19,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2) XIBM._blankingIBM(t, tb, dimPb=2) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) test.testT(t, 1) test.testT(tc, 11) @@ -36,7 +36,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) test.testT(t, 2) test.testT(tc, 21) @@ -55,7 +55,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, correctionMultiCorpsF42=True, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, correctionMultiCorpsF42=True, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) test.testT(t, 3) test.testT(tc, 31) @@ -72,7 +72,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42, blankingF42=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) test.testT(t, 4) test.testT(tc, 41) @@ -89,6 +89,6 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42, twoFronts=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2, twoFronts=True) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2, twoFronts=True) test.testT(t, 5) test.testT(tc, 51) diff --git a/Cassiopee/Connector/test/buildFrontIBMPT_t2.py b/Cassiopee/Connector/test/buildFrontIBMPT_t2.py index c6e83105b..aa8832584 100755 --- a/Cassiopee/Connector/test/buildFrontIBMPT_t2.py +++ b/Cassiopee/Connector/test/buildFrontIBMPT_t2.py @@ -20,7 +20,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3) XIBM._blankingIBM(t, tb, dimPb=3) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3) test.testT(t, 1) test.testT(tc, 11) @@ -37,7 +37,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) test.testT(t, 2) test.testT(tc, 21) @@ -56,7 +56,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, correctionMultiCorpsF42=True, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, correctionMultiCorpsF42=True, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) test.testT(t, 3) test.testT(tc, 31) @@ -73,7 +73,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42, blankingF42=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) test.testT(t, 4) test.testT(tc, 41) @@ -90,6 +90,6 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42, twoFronts=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42, twoFronts=True) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42, twoFronts=True) test.testT(t, 5) test.testT(tc, 51) diff --git a/Cassiopee/Connector/test/setInterpDataIBMPT_t1.py b/Cassiopee/Connector/test/setInterpDataIBMPT_t1.py index 1ba751cbc..631093f27 100755 --- a/Cassiopee/Connector/test/setInterpDataIBMPT_t1.py +++ b/Cassiopee/Connector/test/setInterpDataIBMPT_t1.py @@ -22,7 +22,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2) XIBM._blankingIBM(t, tb, dimPb=2) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=2, frontType=1) test.testT(t, 1) test.testT(tc, 11) @@ -42,7 +42,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=2, frontType=42) test.testT(t, 2) test.testT(tc, 21) @@ -64,7 +64,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, correctionMultiCorpsF42=True, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, correctionMultiCorpsF42=True, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42) XIBM._setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=2, frontType=42) test.testT(t, 3) test.testT(tc, 31) @@ -84,7 +84,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42, blankingF42=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=2, frontType=42) test.testT(t, 4) test.testT(tc, 41) @@ -104,7 +104,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=2, frontType=42) XIBM._blankingIBM(t, tb, dimPb=2, frontType=42, twoFronts=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42, twoFronts=True) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=2, frontType=42, twoFronts=True) XIBM._setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=2, frontType=42, twoFronts=True) test.testT(t, 5) test.testT(tc, 51) diff --git a/Cassiopee/Connector/test/setInterpDataIBMPT_t2.py b/Cassiopee/Connector/test/setInterpDataIBMPT_t2.py index c82db4f1a..7d0edb066 100755 --- a/Cassiopee/Connector/test/setInterpDataIBMPT_t2.py +++ b/Cassiopee/Connector/test/setInterpDataIBMPT_t2.py @@ -20,7 +20,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3) XIBM._blankingIBM(t, tb, dimPb=3) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=3, frontType=1) test.testT(t, 1) test.testT(tc, 11) @@ -38,7 +38,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=3, frontType=42) test.testT(t, 2) test.testT(tc, 21) @@ -58,7 +58,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, correctionMultiCorpsF42=True, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, correctionMultiCorpsF42=True, frontType=42) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=3, frontType=42) test.testT(t, 3) test.testT(tc, 31) @@ -76,7 +76,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42, blankingF42=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42) XIBM._setInterpDataIBM(t, tc, tb, front, dimPb=3, frontType=42) test.testT(t, 4) test.testT(tc, 41) @@ -94,7 +94,7 @@ XIBM._dist2wallIBM(t, tb, dimPb=3, frontType=42) XIBM._blankingIBM(t, tb, dimPb=3, frontType=42, twoFronts=True) tc = C.node2Center(t) -t, tc, front, front2 = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42, twoFronts=True) +t, tc, front, front2, frontWMM = XIBM.buildFrontIBM(t, tc, dimPb=3, frontType=42, twoFronts=True) XIBM._setInterpDataIBM(t, tc, tb, front, front2=front2, dimPb=3, frontType=42, twoFronts=True) test.testT(t, 5) test.testT(tc, 51) diff --git a/Cassiopee/Generator/Generator/IBM.py b/Cassiopee/Generator/Generator/IBM.py index c47497bbf..5bb7f1219 100644 --- a/Cassiopee/Generator/Generator/IBM.py +++ b/Cassiopee/Generator/Generator/IBM.py @@ -364,7 +364,7 @@ def _addBCsForSymmetry(t, bbox=None, dimPb=3, dir_sym=0, X_SYM=0., depth=2): def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[], tbox=None, snearsf=None, check=True, to=None, ext=2, expand=3, dfarDir=0, check_snear=False, mode=0, - tbOneOver=None, listF1save=[]): + tbOneOver=None, listF1save=[], fileoutpre='./'): import KCore.test as test # list of dfars if dfarList == []: @@ -387,7 +387,9 @@ def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[] else: snearsf.append(1.) fileout = None - if check: fileout = 'octree.cgns' + if check: + fileoutpre[-1]='octree.cgns' + fileout = '/'.join(fileoutpre) # Octree identical on all procs if to is not None: if isinstance(to, str): @@ -398,9 +400,9 @@ def generateIBMMeshPara(tb, vmin=15, snears=None, dimPb=3, dfar=10., dfarList=[] parento = None else: o = buildOctree(tb, snears=snears, snearFactor=1., dfar=dfar, dfarList=dfarList, - to=to, tbox=tbox, snearsf=snearsf, - dimPb=dimPb, vmin=vmin, fileout=None, rank=Cmpi.rank, - expand=expand, dfarDir=dfarDir, mode=mode) + to=to, tbox=tbox, snearsf=snearsf, + dimPb=dimPb, vmin=vmin, fileout=None, rank=Cmpi.rank, + expand=expand, dfarDir=dfarDir, mode=mode) if Cmpi.rank==0 and check: C.convertPyTree2File(o,fileout) # build parent octree 3 levels higher @@ -732,6 +734,7 @@ def buildOctree(tb, snears=None, snearFactor=1., dfar=10., dfarList=[], to=None, def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): + tbSolid = Internal.rmNodesByName(tb, 'IBCFil*') boxes = [] for b in Internal.getBases(tbox): boxes.append(Internal.getNodesFromType1(b, 'Zone_t')) @@ -752,7 +755,7 @@ def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): snearl = Internal.getValue(snearl) snearsf.append(snearl*(vmin-1)) - + to = C.newPyTree(['Base', o]) end = 0 G._getVolumeMap(to) @@ -762,7 +765,7 @@ def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): while end == 0: # Do not refine inside obstacles C._initVars(to, 'centers:cellN', 1.) - to = X_IBM.blankByIBCBodies(to, tb, 'centers', dim) + to = X_IBM.blankByIBCBodies(to, tbSolid, 'centers', dim) C._initVars(to, '{centers:cellNBody}={centers:cellN}') nob = 0 C._initVars(to, 'centers:indicator', 0.) @@ -771,7 +774,6 @@ def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): C._initVars(to,'centers:cellN',1.) tboxl = C.newPyTree(['BOXLOC']); tboxl[2][1][2] = box to = X_IBM.blankByIBCBodies(to, tboxl, 'centers', dim) - fact = 1.1 while C.getMinValue(to, 'centers:cellN') == 1 and fact < 10.: print("Info: addRefinementZones: tbox too small - increase tbox by fact = %2.1f"%(fact)) @@ -782,10 +784,11 @@ def addRefinementZones(o, tb, tbox, snearsf, vmin, dim): C._initVars(to,'{centers:indicator}=({centers:indicator}>0.)+({centers:indicator}<1.)*logical_and({centers:cellN}<0.001, {centers:vol}>%g)'%volmin2) nob += 1 - + + end = 1 C._initVars(to,'{centers:indicator}={centers:indicator}*({centers:cellNBody}>0.)*({centers:vol}>%g)'%volmin0) - + if C.getMaxValue(to, 'centers:indicator') == 1.: end = 0 # Maintien du niveau de raffinement le plus fin @@ -830,9 +833,10 @@ def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6 if parento is None: ## Rectilinear mesh modifications + ZONEStbOneOver = None + ZONEStbOneOverTmp = None if tbOneOver: listSavetbOneOverZones = [] - ZONEStbOneOver = None ##Add zones to list of TbOneOver list & ##Remove zone from orig list of zones @@ -853,8 +857,8 @@ def octree2StructLoc__(o, parento=None, vmin=21, ext=0, optimized=0, sizeMax=4e6 zones = T.mergeCart(zones,sizeMax=sizeMax) else: ## Rectilinear mesh modifications - ZONEStbOneOver = None - + ZONEStbOneOver = None + ZONEStbOneOverTmp = None eps = 1.e-10 bbo = G.bbox(parento[0])# 1st octant lower left side xmeano=bbo[3]; ymeano=bbo[4]; zmeano=bbo[5] diff --git a/Cassiopee/Geom/Geom/IBM.py b/Cassiopee/Geom/Geom/IBM.py index c78fb05be..deb436967 100644 --- a/Cassiopee/Geom/Geom/IBM.py +++ b/Cassiopee/Geom/Geom/IBM.py @@ -505,4 +505,36 @@ def _transformTc2(tc2): return None - +#================================================================================ +# Selection/determination of tb for closed solid & filament +#================================================================================ +def determineClosedSolidFilament__(tb): + ##OUT - filamentBases :list of names of open geometries + ##OUT - isFilamentOnly :boolean if there is only a filament in tb + ##OUT - isOrthoProjectFirst :bolean to do orthonormal projection first + ##OUT - tb :tb of solid geometries only + ##OUT - tbFilament :tb of filament geometries only + + ## General case where only a closeSolid is in tb + ## or tb only has a filament + 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 + isOrthoProjectFirst = isFilamentOnly + + ## if tb has both a closed solid and filaments + tbFilament = Internal.copyTree(tb) + 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]