correct conflict (4aebd637) · Commits · DRTI / BasicTools

src/BasicTools/Bridges/CGNSBridge.py

+247 −181

Original line number	Diff line number	Diff line
		@@ -7,7 +7,11 @@ import numpy as np

		from BasicTools.NumpyDefs import PBasicFloatType, PBasicIndexType
		from BasicTools.Containers.UnstructuredMesh import UnstructuredMesh, ElementsContainer
		from BasicTools.Containers import UnstructuredMeshInspectionTools as UMIT
		from BasicTools.Containers import UnstructuredMeshModificationTools as UMMT
		from BasicTools.Containers import Filters as F
		import BasicTools.Containers.ElementNames as EN
		from typing import List

		CGNSNameToBasicTools = {}
		CGNSNameToBasicTools["NODE"] = EN.Point_1
		@@ -81,6 +85,7 @@ def __ReadIndexRange(pyTree):
		return np.array(res, dtype = PBasicIndexType).ravel()

		def GetRecursiveObjects(pyTree, filter, state= None):

		res = list()
		state, lres = filter(state, pyTree)
		res.extend(lres)
		@@ -113,37 +118,63 @@ def GetNodeByPath(pyTree, path):

		return GetRecursiveObjects(pyTree, fil, path.strip("/").split("/"))[1]

		def CGNSToMesh(pyTree, baseNumberOrName= 0, zoneNumberOrName = 0)-> UnstructuredMesh:

		res = UnstructuredMesh()
		def CGNSToMesh(pyTree, baseNames:List[str] = None, zoneNames:List[str] = None, FieldsAsTags:bool=False)-> UnstructuredMesh:

		if type(baseNumberOrName) is int:
		basepath = GetAllNodesByTypeSet(pyTree,"CGNSBase_t")[baseNumberOrName]
		if baseNames != None:
		basepaths0 = GetAllNodesByTypeSet(pyTree,"CGNSBase_t")
		availableBaseNames = [b.split('/')[-1] for b in basepaths0]
		basepaths = []
		for baseName in baseNames:
		if baseName in availableBaseNames:
		basepaths.append(basepaths0[availableBaseNames.index(baseName)])
		else:
		raise("baseName "+baseName+" not available in CGNS tree")
		else:
		basepath = baseNumberOrName
		basepaths = GetAllNodesByTypeSet(pyTree,"CGNSBase_t")

		base = GetNodeByPath(pyTree,basepath)[0]

		if type(zoneNumberOrName) is int:
		zonepath = GetAllNodesByTypeSet(base,"Zone_t")[zoneNumberOrName]
		else:
		zonepath = zoneNumberOrName
		res = None

		for basepath in basepaths:

		basePyTree = GetNodeByPath(pyTree, basepath)[0]

		zonePyTree = GetNodeByPath(base,"/"+ zonepath)[0]
		if zoneNames != None:

		zonepaths0 = GetAllNodesByTypeSet(basePyTree,"Zone_t")
		availableZoneNames = [b.split('/')[-1] for b in zonepaths0]
		zonepaths = []
		for zoneName in zoneNames:
		if zoneName in availableZoneNames:
		zonepaths.append(zonepaths0[availableZoneNames.index(zoneName)])
		else:
		raise("zoneName "+zoneName+" not available in base "+basepath)
		else:
		zonepaths = GetAllNodesByTypeSet(basePyTree,"Zone_t")


		for zonepath in zonepaths:
		zonePyTree = GetNodeByPath(basePyTree, zonepath)[0]
		gridCoordinatesPath = GetAllNodesByTypeSet(zonePyTree, "GridCoordinates_t")[0]

		#nodes
		if len(GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateZ')) == 0:
		node_dim = 2
		else:
		node_dim = 3

		gx = GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateX')[0][1]
		res.nodes = np.empty((gx.shape[0],3), dtype= PBasicFloatType)
		res.nodes[:,0] = gx
		res.nodes[:,1] = GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateY')[0][1]
		res.nodes[:,2] = GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateZ')[0][1]
		res.originalIDNodes = np.arange(1, res.nodes.shape[0]+1, dtype= PBasicIndexType)

		local_mesh = UnstructuredMesh()
		local_mesh.nodes = np.empty((gx.shape[0],node_dim), dtype= PBasicFloatType)
		local_mesh.nodes[:,0] = gx
		local_mesh.nodes[:,1] = GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateY')[0][1]
		if node_dim == 3:
		local_mesh.nodes[:,2] = GetNodeByPath(zonePyTree,gridCoordinatesPath+'/CoordinateZ')[0][1]
		local_mesh.originalIDNodes = np.arange(1, local_mesh.nodes.shape[0]+1, dtype= PBasicIndexType)

		#elements
		elementsPaths = GetAllNodesByTypeSet(zonePyTree, "Elements_t")

		for elementsPath in elementsPaths:
		cgnsElements = GetNodeByPath(zonePyTree,elementsPath)[0]
		cgnsUserName = "Elements_" + cgnsElements[0]
		@@ -161,46 +192,24 @@ def CGNSToMesh(pyTree, baseNumberOrName= 0, zoneNumberOrName = 0)-> Unstructured
		if CGNSToBasicToolsPermutation.get(cgnsElemType,None) is not None:
		coon = coon[CGNSToBasicToolsPermutation[cgnsElemType]]
		cpt += nbp
		elems: ElementsContainer = res.elements.GetElementsOfType(basicToolsElemType)
		elems: ElementsContainer = local_mesh.elements.GetElementsOfType(basicToolsElemType)
		elcpt = elems.AddNewElements(coon[None,:], [oid] )
		elems.tags.CreateTag(cgnsUserName,False).AddToTag(elcpt-1)
		else:
		basicToolsElemType = GetCGNSNumberToBasicTools(cgnsElemType)
		elementConnectivity = np.asarray(GetNodeByPath(cgnsElements,cgnsUserName+'/ElementConnectivity')[0][1], dtype=PBasicIndexType).reshape( (-1, EN.numberOfNodes[basicToolsElemType] ) )-1
		elems: ElementsContainer = res.elements.GetElementsOfType(basicToolsElemType)

		elems: ElementsContainer = local_mesh.elements.GetElementsOfType(basicToolsElemType)

		if CGNSToBasicToolsPermutation.get(cgnsElemType,None) is not None:
		elementConnectivity = elementConnectivity[:,CGNSToBasicToolsPermutation[cgnsElemType]]
		elcpt = elems.cpt

		nelcpt = elems.AddNewElements(elementConnectivity, originalIds)
		#elems.tags.CreateTag(cgnsUserName,False).AddToTag(np.arange(elcpt,nelcpt))

		elems.tags.CreateTag(cgnsUserName,False).AddToTag(np.arange(elcpt,nelcpt))

		datasPaths = GetAllNodesByTypeSet(zonePyTree, "FlowSolution_t")
		for dataPath in datasPaths :
		datas = GetNodeByPath(zonePyTree,dataPath)[0]
		gl = GetAllNodesByTypeSet(datas, "GridLocation_t")

		store = res.nodeFields
		if len(gl) > 0:
		if "".join(np.array(GetNodeByPath(datas,gl[0])[0][1], dtype=str) )== "CellCenter":
		store = res.elemFields
		if "".join(np.array(GetNodeByPath(datas,gl[0])[0][1], dtype=str) )== "Vertex":
		store = res.nodeFields

		fieldPaths = GetAllNodesByTypeSet(datas, "DataArray_t")
		for fieldPath in fieldPaths:
		fieldData = GetNodeByPath(datas,fieldPath)[0]
		dataName = fieldData[0]
		data = fieldData[1]
		if dataName == "OriginalIds" and store is res.nodeFields:
		res.originalIDNodes = np.asarray(data, dtype=PBasicIndexType,order='C')
		elif dataName == "OriginalIds" and store is res.elemFields:
		res.SetElementsOriginalIDs(np.asarray(data, dtype=PBasicIndexType))
		else:
		res.nodeFields[dataName] = np.asarray(data)

		#tags
		elemTags = []
		nodeTags = []
		ZoneBCPaths = GetAllNodesByTypeSet(zonePyTree, "ZoneBC_t")
		for ZoneBCPath in ZoneBCPaths:
		ZoneBC = GetNodeByPath(zonePyTree,ZoneBCPath)[0]
		@@ -216,14 +225,53 @@ def CGNSToMesh(pyTree, baseNumberOrName= 0, zoneNumberOrName = 0)-> Unstructured

		gl = GetAllNodesByTypeSet(BCNode, "GridLocation_t")
		if "".join(np.array(GetNodeByPath(BCNode,gl[0])[0][1], dtype=str) )== "CellCenter":
		#res.AddElementToTagUsingOriginalId(indices-1, BCName)
		#res.AddElementsToTag(indices-1, BCName)
		print(f"Reading Element selections is not supported. skiping {BCName}")
		#local_mesh.AddElementToTagUsingOriginalId(indices-1, BCName)
		local_mesh.AddElementsToTag(indices-1, BCName)
		elemTags.append(BCName)
		#print(indices)
		#print(f"Reading Element selections is not supported. skiping {BCName}")
		pass
		if "".join(np.array(GetNodeByPath(BCNode,gl[0])[0][1], dtype=str) )== "Vertex":
		res.nodesTags.CreateTag(BCName).SetIds(indices-1)
		local_mesh.nodesTags.CreateTag(BCName).SetIds(indices-1)
		nodeTags.append(BCName)

		#fields
		datasPaths = GetAllNodesByTypeSet(zonePyTree, "FlowSolution_t")
		for dataPath in datasPaths :
		datas = GetNodeByPath(zonePyTree,dataPath)[0]
		gl = GetAllNodesByTypeSet(datas, "GridLocation_t")

		store = local_mesh.nodeFields
		if len(gl) > 0:
		if "".join(np.array(GetNodeByPath(datas,gl[0])[0][1], dtype=str) )== "CellCenter":
		store = local_mesh.elemFields
		if "".join(np.array(GetNodeByPath(datas,gl[0])[0][1], dtype=str) )== "Vertex":
		store = local_mesh.nodeFields

		fieldPaths = GetAllNodesByTypeSet(datas, "DataArray_t")
		for fieldPath in fieldPaths:
		fieldData = GetNodeByPath(datas,fieldPath)[0]
		dataName = fieldData[0]
		data = fieldData[1]
		if dataName == "OriginalIds" and store is local_mesh.nodeFields:
		local_mesh.originalIDNodes = np.asarray(data, dtype=PBasicIndexType,order='C')
		elif dataName == "OriginalIds" and store is local_mesh.elemFields:
		local_mesh.SetElementsOriginalIDs(np.asarray(data, dtype=PBasicIndexType))
		elif store is local_mesh.nodeFields:
		if (FieldsAsTags == True) or (FieldsAsTags == False and dataName not in nodeTags):
		local_mesh.nodeFields[dataName] = np.asarray(data)
		elif store is local_mesh.elemFields:
		if (FieldsAsTags == True) or (FieldsAsTags == False and dataName not in elemTags):
		local_mesh.elemFields[dataName] = np.asarray(data)

		if res == None:
		res = local_mesh
		else:
		res.Merge(local_mesh)
		#UMMT.CleanDoubleElements(res)
		UMMT.CleanDoubleNodes(res)


		res.PrepareForOutput()
		return res

		def NewDataArray(father, name, data):
		@@ -261,107 +309,125 @@ def NewFlowSolution(father,name,gridlocation ):
		father[2].append(res)
		return res

		def MeshToCGNS( mesh: UnstructuredMesh, outputPyTree = None, baseNumberOrName= 0, zoneNumberOrName = 0) :

		if type(baseNumberOrName) is int:
		baseName = f"Base_{baseNumberOrName}"
		else:
		baseName = f"{baseNumberOrName}"

		if type(zoneNumberOrName) is int:
		zoneName = f"Zone_{baseNumberOrName}"
		else:
		zoneName = f"{baseNumberOrName}"
		def MeshToCGNS( mesh: UnstructuredMesh, outputPyTree = None, OneBasePerTopoDim:bool = False, TagsAsFields:bool=False) :

		if outputPyTree is None:
		outputPyTree=['CGNSTree', None, [], 'CGNSTree_t']
		version = ['CGNSLibraryVersion', np.array([3.4], dtype=np.float32), [], 'CGNSLibraryVersion_t']
		outputPyTree[2].append(version)
		physicalDim = mesh.GetPointsDimensionality()
		topologicalDim = mesh.GetElementsDimensionality()
		nbElemTags = len(mesh.GetNamesOfElemTags())
		nbNodesTags = len(mesh.nodesTags)
		totalNumberOfTags = nbElemTags +nbNodesTags

		base = [baseName, np.array([physicalDim, topologicalDim], dtype=np.int32), [], 'CGNSBase_t'] # name, physical dim, topological dim
		if mesh.GetPointsDimensionality() == 2:
		physicalDim = 2
		coordNames = ["X", "Y"]
		elif mesh.GetPointsDimensionality() == 3:
		physicalDim = 3
		coordNames = ["X", "Y", "Z"]
		else:
		raise('mesh.GetPointsDimensionality() should be 2 or 3')


		def ConstructZone(outputPyTree, local_mesh:UnstructuredMesh, topologicalDim:int, physicalDim:int, baseName:str, zoneName:str)->None:
		base = [baseName, np.array([topologicalDim, physicalDim], dtype=np.int32), [], 'CGNSBase_t'] # name, physical dim, topological dim
		outputPyTree[2].append(base)

		family = ['Bulk', np.array([b'B', b'u', b'l', b'k'], dtype='\|S1'), [], 'FamilyName_t']
		base[2].append(family)

		s = np.array([[mesh.GetNumberOfNodes(),mesh.GetNumberOfElements(),totalNumberOfTags]],dtype=np.int32)
		nbElemTags = len(local_mesh.GetNamesOfElemTags())
		nbNodesTags = len(local_mesh.nodesTags)
		totalNumberOfTags = nbElemTags +nbNodesTags
		totalNumberOfElem = local_mesh.GetNumberOfElements()

		s = np.array([[local_mesh.GetNumberOfNodes(),local_mesh.GetNumberOfElements(),totalNumberOfTags]],dtype=np.int32)
		grid = [zoneName, s, [['ZoneType', np.array([b'U', b'n', b's', b't', b'r', b'u', b'c', b't', b'u', b'r', b'e',b'd'], dtype='\|S1'), [], 'ZoneType_t']], 'Zone_t']
		base[2].append(grid)
		zone_family = ['FamilyName', np.array([b'B', b'u', b'l', b'k'], dtype='\|S1'), [], 'FamilyName_t']
		grid[2].append(zone_family)

		else:
		grid = GetNodeByPath(outputPyTree,f"/{baseName}/{zoneName}/")[0][1]

		#add nodes
		gridCoordinates = ['GridCoordinates', None, [], 'GridCoordinates_t']
		grid[2].append(gridCoordinates)

		if mesh.GetPointsDimensionality() == 2:
		coordNames = ["X", "Y"]
		elif mesh.GetPointsDimensionality() == 3:
		coordNames = ["X", "Y", "Z"]
		else:
		raise('mesh.GetPointsDimensionality() should be 2 or 3')
		for i,coord in enumerate(coordNames):
		da = NewDataArray(gridCoordinates, "Coordinate"+coord, np.copy(mesh.nodes[:,i]) )


		#nodes tags
		if len(mesh.nodesTags):
		zbg = ['Points_Selections', None, [], 'ZoneBC_t']
		grid[2].append(zbg)
		for tag in mesh.nodesTags:
		NewBC(zbg,tag.name,tag.GetIds()+1,pttype="IndexArray_t")
		da = NewDataArray(gridCoordinates, "Coordinate"+coord, np.copy(local_mesh.nodes[:,i]) )

		for name, data in mesh.elements.items():
		#elem tags
		nbelem = data.GetNumberOfElements()
		#add elements
		for name, data in local_mesh.elements.items():
		nbelem = data.connectivity.shape[0]
		if nbelem == 0:
		continue
		NewElements(grid, data)

		ezbg = ['Elements_Selections', None, [], 'ZoneBC_t']
		grid[2].append(ezbg)

		for name in mesh.GetNamesOfElemTags():
		ids = mesh.GetElementsInTag(name)+1
		NewBC(ezbg,name,ids, pttype="IndexArray_t", onPoints=False)

		## add nodeFields
		flowSolution = NewFlowSolution(grid, "PointData", gridlocation="Vertex")
		NewDataArray(flowSolution, "OriginalIds", mesh.originalIDNodes )
		if len(mesh.nodeFields) > 0:
		for name, pointData in mesh.nodeFields.items():
		NewDataArray(flowSolution, "OriginalIds", local_mesh.originalIDNodes)
		if len(local_mesh.nodeFields) > 0:
		for name, pointData in local_mesh.nodeFields.items():
		if pointData.dtype.char == "U":
		print(f"skipping nodeFields '{name}' because if of type: {pointData.dtype}" )
		continue
		if name=='Normals':
		assert(pointData.shape[1] == mesh.GetPointsDimensionality())
		for i,coord in enumerate(coordNames):
		NewDataArray(flowSolution, "Normals"+coord, np.copy(pointData[:,i]) )
		else:
		NewDataArray(flowSolution, name,pointData)

		## elem fields
		## add elemFields
		flowSolutionCell = NewFlowSolution(grid, "CellData", gridlocation="CellCenter")
		NewDataArray(flowSolutionCell, "OriginalIds", mesh.GetElementsOriginalIDs() )
		if len(mesh.elemFields) > 0:
		for name, cellData in mesh.elemFields.items():
		NewDataArray(flowSolutionCell, "OriginalIds", local_mesh.GetElementsOriginalIDs() )
		if len(local_mesh.elemFields) > 0:
		for name, cellData in local_mesh.elemFields.items():
		if cellData.dtype.char == "U":
		if name != 'FE Names':
		print(f"skipping elemFields '{name}' because if of type: {cellData.dtype}" )
		continue
		if name=='Normals':
		assert(cellData.shape[1] == mesh.GetPointsDimensionality())
		for i,coord in enumerate(coordNames):
		NewDataArray(flowSolutionCell, "Normals"+coord, np.copy(cellData[:,i]) )
		else:
		NewDataArray(flowSolutionCell, name, np.copy(cellData))

		#add elements tags (as field also)
		has_elem_select = False
		for name, data in local_mesh.elements.items():
		if len(data.tags.keys())>0:
		has_elem_select = True
		if has_elem_select:
		ezbg = ['Elements_Selections', None, [], 'ZoneBC_t']
		grid[2].append(ezbg)

		for name, data in local_mesh.elements.items():
		for elTag in data.tags.keys():
		ids = data.tags[elTag].GetIds()
		NewBC(ezbg,elTag,ids+1, pttype="IndexArray_t", onPoints=False)
		if TagsAsFields:
		cellData = np.zeros(totalNumberOfElem, dtype=np.int64)
		cellData[ids] = 1
		NewDataArray(flowSolutionCell, elTag+"_as_field", cellData)

		#add node tags (as field also)
		if len(local_mesh.nodesTags):
		zbg = ['Points_Selections', None, [], 'ZoneBC_t']
		grid[2].append(zbg)
		for tag in local_mesh.nodesTags:
		ids = tag.GetIds()
		NewBC(zbg,tag.name,ids+1,pttype="IndexArray_t")
		if TagsAsFields:
		pointData = np.zeros(local_mesh.GetNumberOfNodes())
		pointData[ids] = 1
		NewDataArray(flowSolution, tag.name+"_as_field", np.copy(pointData))

		el_topo_dim={1:[],2:[],3:[]}

		maxTopoDim = 0
		for name in mesh.elements.keys():
		el_topo_dim[EN.dimension[name]].append(name)
		maxTopoDim = max(maxTopoDim, EN.dimension[name])

		if OneBasePerTopoDim:
		for td in range(3):
		topologicalDim = td+1
		if len(el_topo_dim[topologicalDim])>0:
		ff = F.ElementFilter(mesh, dimensionality=topologicalDim)
		local_mesh = UMIT.ExtractElementsByElementFilter(mesh, ff)
		ConstructZone(outputPyTree, local_mesh, topologicalDim, physicalDim, "Base_"+str(topologicalDim)+"_"+str(physicalDim), "Zone")

		else:
		#ff = F.ElementFilter(mesh, dimensionality=maxTopoDim)
		#local_mesh = UMIT.ExtractElementsByElementFilter(mesh, ff)
		ConstructZone(outputPyTree, mesh, maxTopoDim, physicalDim, "Base_"+str(maxTopoDim)+"_"+str(physicalDim), "Zone")

		return outputPyTree