# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
import numpy as np
from Muscat.Types import MuscatIndex, MuscatFloat
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.Helpers.Logger import Debug
[docs]
class DofNumberingDict( ):
def __init__(self):
"""Class to generate the numbering.
The internal implementation uses Dict.
"""
super().__init__()
self.numbering = dict()
self.size = 0
self._doftopointLeft = None
self._doftopointRight= None
self._doftocellLeft = None
self._doftocellRight = None
self.almanac = dict()
self.newAlmanac = dict()
self.fromConnectivity = True
self.discontinuous = False
def __getitem__(self,key):
if key == "size":
# print("Please use the new API of DofNumbering : DofNumbering.size")
return self.size
if key == "fromConnectivity":
# print("Please use the new API of DofNumbering : DofNumbering.fromConnectivity")
return self.fromConnectivity
if key == "almanac":
# print("Please use the new API of DofNumbering : DofNumbering.almanac")
return self.almanac
if key == "doftopointLeft":
# print("Please use the new API of DofNumbering : DofNumbering.doftopointLeft")
return self.doftopointLeft
if key == "doftopointRight":
# print("Please use the new API of DofNumbering : DofNumbering.doftopointRight")
return self.doftopointRight
if key == "doftocellLeft":
# print("Please use the new API of DofNumbering : DofNumbering.doftopointLeft")
return self.doftocellLeft
if key == "doftocellRight":
# print("Please use the new API of DofNumbering : DofNumbering.doftopointRight")
return self.doftocellRight
return self.numbering[key]
[docs]
def get(self,key,default=None):
if key in self.numbering:
return self.numbering[key]
else:
return default
def __contains__(self, k):
return k in self.numbering
[docs]
def GetSize(self):
return self.size
[docs]
def ComputeNumberingFromConnectivity(self,mesh,space):
self.size = mesh.GetNumberOfNodes()
self._doftopointLeft = range(self.size)
self._doftopointRight = range(self.size)
self._doftocellLeft = []
self._doftocellRight = []
self.fromConnectivity = True
for name,data in mesh.elements.items():
self.numbering[name] = data.connectivity
almanac = {}
for i in range(self.size):
almanac[('P', i, None)] = i
self.almanac = almanac
self.mesh = mesh
return self
[docs]
def ComputeNumberingNoSharedDof(self, mesh, space, elementFilter):
self.fromConnectivity = False
self.discontinuous = True
cpt = self.size
if elementFilter is None:
elementFilter = ElementFilter()
for selection in elementFilter(mesh):
nbShapeFunctions = space[selection.elementType].GetNumberOfShapeFunctions()
if selection.elementType in self.numbering:
dofs = self.numbering[selection.elementType]
else:
dofs = np.zeros((selection.elements.GetNumberOfElements(),nbShapeFunctions ), dtype=MuscatIndex) -1
dofs[selection.indices,:] = np.arange(start=cpt,stop=selection.Size()*nbShapeFunctions, dtype=MuscatIndex).reshape(selection.Size(),nbShapeFunctions)
self.numbering[selection.elementType] = dofs
self.mesh = mesh
return self
[docs]
def ComputeNumberingGeneral(self, mesh, space, elementFilter=None):
self.fromConnectivity = False
almanac = self.almanac
if elementFilter is None:
elementFilter = ElementFilter()
cpt = self.size
Debug("bulk ")
usedDim = 0
for selection in elementFilter(mesh):
usedDim = max(usedDim,ED.dimensionality[selection.elementType] )
res = self.GetHashFor(selection.elements, space[selection.elementType], selection.indices, False)
if selection.elementType in self.numbering:
dofs = self.numbering[selection.elementType]
else:
dofs = np.zeros((selection.elements.GetNumberOfElements(),space[selection.elementType].GetNumberOfShapeFunctions()), dtype=MuscatIndex) -1
Debug(f"{selection.elementType} Done")
for i in range(len(res)):
localRes = res[i]
localDofs = dofs[:,i]
for j, elementId in enumerate(selection.indices):
d = almanac.setdefault(localRes[j],cpt)
cpt += (d == cpt)
localDofs[elementId] = d
self.numbering[selection.elementType] = dofs
Debug(f"{selection.elementType} Done Done")
Debug("bulk Done")
Debug("complementary ")
from Muscat.MeshContainers.Filters.FilterOperators import IntersectionFilter, ComplementaryFilter
outside = IntersectionFilter(filters =[ElementFilter(dimensionality=usedDim-1), ComplementaryFilter( filters = [elementFilter])] )
for selection in outside(mesh):
res = self.GetHashFor(selection.elements,space[selection.elementType],selection.indices,False)
if selection.elementType in self.numbering:
dofs = self.numbering[selection.elementType]
else:
dofs = self.numbering.setdefault(selection.elementType,np.zeros((selection.elements.GetNumberOfElements(),space[selection.elementType].GetNumberOfShapeFunctions()), dtype=MuscatIndex) -1)
for i in range(len(res)):
localRes = res[i]
localDofs = dofs[:,i]
for j,elementId in enumerate(selection.indices):
if localDofs[elementId] >= 0:
continue
localDofs[elementId] = almanac.get(localRes[j],-1)
self.numbering[selection.elementType] = dofs
self.size = cpt
#-------------------------------------------------------------------------
self.mesh = mesh
# we keep a reference to the mesh because we need it to compute the
return self
[docs]
def GetDofOfPoint(self,pid):
return self.almanac[("P",pid,None)]
@property
def doftopointLeft(self):
if self._doftopointLeft is None:
self.computeDofToPoint()
return self._doftopointLeft
@property
def doftopointRight(self):
if self._doftopointRight is None:
self.computeDofToPoint()
return self._doftopointRight
[docs]
def computeDofToPoint(self):
extractorLeftSide = np.empty(self.size,dtype=MuscatIndex)
extractorRightSide = np.empty(self.size,dtype=MuscatIndex)
cpt = 0
# if k[0] is 'P' then k[1] is the node number
for k,v in self.almanac.items():
if k[0] == 'P':
extractorLeftSide[cpt] = k[1]
extractorRightSide[cpt] = v
cpt += 1
self._doftopointLeft = np.resize(extractorLeftSide, (cpt,))
self._doftopointRight = np.resize(extractorRightSide, (cpt,))
@property
def doftocellLeft(self):
if self._doftocellLeft is None:
self.computeDofToCell()
return self._doftocellLeft
@property
def doftocellRight(self):
if self._doftocellRight is None:
self.computeDofToCell()
return self._doftocellRight
[docs]
def computeDofToCell(self):
mesh = self.mesh
extractorLeftSide = np.empty(self.size,dtype=MuscatIndex)
extractorRightSide = np.empty(self.size,dtype=MuscatIndex)
cpt = 0
# if k[0] is the elementName then k[1] is the connectivity
# we generate the same almanac with the number of each element
elemDic = {}
for name, data in mesh.elements.items():
elemDic[name] = {}
elemDic2 = elemDic[name]
sortedConnectivity = np.sort(data.connectivity,axis=1)
for i in range(data.GetNumberOfElements()):
elemDic2[tuple(sortedConnectivity[i,:])] = i
for k,v in self.almanac.items():
#if not k[0] in {'P',"F","F2","G"} :
#we need the global number of the element (not the local to the element container)
if k[0] in elemDic.keys():
localDict = elemDic[k[0]]
if k[1] in localDict.keys():
extractorLeftSide[cpt] = mesh.elements[k[0]].globaloffset + localDict[k[1]]
extractorRightSide[cpt] = v
cpt += 1
self._doftocellLeft = np.resize(extractorLeftSide, (cpt,))
self._doftocellRight = np.resize(extractorRightSide, (cpt,))
[docs]
def GetHashFor(self, data, sp, elementIds, discontinuous):
numberOfShapeFunctions = sp.GetNumberOfShapeFunctions()
res = []
name = data.elementType
elementsIdsConnectivity = data.connectivity[elementIds,:]
for j in range(numberOfShapeFunctions):
on,idxI,idxII = sp.dofAttachments[j]
if on == "P":
T = "P"
shapeFunctionConnectivity = elementsIdsConnectivity [:,idxI]
if discontinuous :
res.append( [ (T+str(elementIds),x,idxII) for i,x in zip(elementIds,shapeFunctionConnectivity) ] )
else:
res.append( [ (T,x,idxII) for x in shapeFunctionConnectivity ] )
elif on == "C":
sortedConnectivity = np.sort(elementsIdsConnectivity,axis=1)
T = name
if idxII is not None:
raise
res.append([(T,tuple(sortedConnectivity[i,:]),idxI) for i in range(len(elementIds)) ] )
elif on == "F2":
edge = ED.faces2[name][idxI]
T = edge[0]
nn = np.sort(elementsIdsConnectivity[:,edge[1]],axis=1)
if discontinuous :
res.append( [ (T+str(elementIds),tuple(x) ,0) for i,x in zip(elementIds,nn) ] )
else:
res.append( [ (T,tuple(x),0) for x in nn ] )
elif on == "F":
edge = ED.faces1[name][idxI]
T = edge[0]
nn = np.sort(elementsIdsConnectivity[:,edge[1]],axis=1)
if discontinuous :
res.append( [ (T,tuple(x),i) for i,x in zip(elementIds,nn) ] )
else:
res.append( [ (T,tuple(x),0) for x in nn ] )
elif on == "G":
"""G is for global """
key = (on,0,idxII)
res.append( [ key for x in elementIds ] )
elif on == "IP":
res.append( [ (on,tuple(localCoon),idxI) for localCoon,i in zip(elementsIdsConnectivity,elementIds) ] )
else:
print(on)
raise
return res
[docs]
def CheckIntegrity(GUI:bool=False):
import Muscat.FE.DofNumbering as DN
return DN.CheckIntegrityUsingAlgo("DictBase",GUI)
if __name__ == '__main__':
print(CheckIntegrity(True))# pragma: no cover