# -*- 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
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Helpers.Logger import Debug, Info
[docs]class DofNumberingNumpy():
def __init__(self):
super().__init__()
self.numbering = dict()
self.size = 0
self._doftopointLeft = None
self._doftopointRight = None
self._doftocellLeft = None
self._doftocellRight = None
self._almanac = None
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 == "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.fromConnectivity = True
self.discontinuous = False
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
self._almanac = dict()
data = np.arange(self.size)
self._almanac[('P', 'P')] = (data, data, data)
self.totalNumberOfPoints = mesh.GetNumberOfNodes()
self.computeDofToPoint()
return self
[docs] def GetKeyFromNameAndIdxI(self, on, name, idxI):
if on == 'C':
on = ('C', name)
elif on == 'F2':
on = ('C', ED.faces2[name][idxI][0])
elif on == 'F':
on = ('C', ED.faces1[name][idxI][0])
elif on == "IP":
on = ("IP", name)
elif on == "G":
on = ('G', 'G')
elif on == 'P':
on = ('P', 'P')
else:
print(on)
raise
return on
[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):
if self._almanac is not None:
raise (Exception("Cant pass 2 time in ComputeNumberingGeneral"))
self.fromConnectivity = False
self.discontinuous = False
if elementFilter is None:
elementFilter = ElementFilter()
from Muscat.Containers.Filters.FilterTools import GetFrozenFilter
elementFilter = GetFrozenFilter(elementFilter, mesh=mesh)
self.totalNumberOfPoints = mesh.GetNumberOfNodes()
def GetNumberingColsSize(k):
if k[0] == 'C':
nn = ED.numberOfNodes[k[1]]
elif k[0] == "IP":
nn = ED.numberOfNodes[k[1]] + 1
elif k[0] == "G":
nn = 1
elif k[0] == "P":
nn = 1
# elif k == "G":
# nn = 1
# elif k == 'P':
# nn = 1
else:
print(k)
raise
nn = ED.numberOfNodes[k]
return nn
sizes = dict()
Info("Numbering counting ")
# count the total number of shape functions per element
for selection in elementFilter(mesh):
sp = space[selection.elements.elementType]
nsf = sp.GetNumberOfShapeFunctions()
for sf in range(nsf):
on, idxI, idxII = sp.dofAttachments[sf]
key = self.GetKeyFromNameAndIdxI(on, selection.elements.elementType, idxI)
sizes[key] = sizes.get(key, 0) + len(selection.indices)
Debug("Numbering memory allocation")
# allocation of matrices to store all the dofs
storage = {}
for k, v in sizes.items():
nn = GetNumberingColsSize(k)
storage[k] = np.zeros((v, nn), dtype=MuscatIndex)-1
sizes[k] = 0
Debug("Numbering generation dofs keys")
for selection in elementFilter(mesh):
Debug(f"working on {selection.elements.elementType}")
# BOO.ResetStartTime()
sp = space[selection.elements.elementType]
nsf = sp.GetNumberOfShapeFunctions()
elementIdsConnectivity = selection.elements.connectivity[selection.indices, :]
for sf in range(nsf):
nn = GetNumberingColsSize(k)
key, idxI, idxII = self.GetHashFor(selection.elements, sp, selection.indices, sf, False, elidsConnectivity=elementIdsConnectivity)
cpt = sizes[key]
storage[key][cpt:cpt+len(selection.indices), :] = idxI
sizes[key] = cpt + len(selection.indices)
cpt = self.size
tempAlmanac = dict()
Debug("Numbering generation of uniques")
# recover the unique dofs and generate the numbering
for k, v in storage.items():
unique, indices, inverse = np.unique(np.sort(v, axis=1), return_index=True, return_inverse=True, axis=0)
newDofs = np.arange(len(indices)) + cpt
tempAlmanac[k] = (unique, newDofs, inverse)
cpt += len(indices)
sizes[k] = 0
self.size = cpt
Debug("Numbering the bulk")
# push the new dofs to the almanac
maxUsedDim = 0
extractorLeftSide = np.empty(self.size, dtype=MuscatIndex)
extractorRightSide = np.empty(self.size, dtype=MuscatIndex)
extractorcpt = 0
elementcpt = 0
for elementType, data in mesh.elements.items():
ids = elementFilter.GetIdsToTreat(mesh, elementType=elementType)
if len(ids) == 0:
elementcpt += data.GetNumberOfElements()
continue
maxUsedDim = max(maxUsedDim, ED.dimensionality[elementType])
sp = space[elementType]
nsf = sp.GetNumberOfShapeFunctions()
self.numbering[elementType] = np.zeros((data.GetNumberOfElements(), nsf), dtype=MuscatIndex)-1
done = False
for sf in range(nsf):
nn = GetNumberingColsSize(k)
on, idxI, idxII = sp.dofAttachments[sf]
key = self.GetKeyFromNameAndIdxI(on, elementType, idxI)
cpt = sizes[key]
(unique, newDofs, inverse) = tempAlmanac[key]
dofs = newDofs[inverse][cpt:cpt+len(ids)]
self.numbering[elementType][ids, sf] = dofs
sizes[key] = cpt + len(ids)
# the dofs are attached to the the elements
#print(elemName,key, done)
if key[1] == elementType and not done:
done = True
extractorLeftSide[extractorcpt:extractorcpt+len(ids)] = ids
extractorLeftSide[extractorcpt:extractorcpt+len(ids)] += elementcpt
extractorRightSide[extractorcpt:extractorcpt+len(ids)] = dofs
extractorcpt += len(ids)
# raise
elementcpt += data.GetNumberOfElements()
self._doftocellLeft = np.resize(extractorLeftSide, (extractorcpt,))
self._doftocellRight = np.resize(extractorRightSide, (extractorcpt,))
from Muscat.Containers.Filters.FilterOperators import IntersectionFilter, ComplementaryFilter, UnionFilter
# we work on the elements of dim < maxuseddim not present in the original filter
complementary = ComplementaryFilter(filters=[elementFilter])
filters = [ElementFilter(dimensionality=i) for i in range(maxUsedDim)]
# filters.append(complementary)
outside = IntersectionFilter(filters=[UnionFilter(filters=filters), complementary])
# push numbering for elements touching the already numbered elements
# for example the 2D elements on the surface of 3D elements
Info("Numbering the complementary")
for selection in outside(mesh):
sp = space[selection.elementType]
nsf = sp.GetNumberOfShapeFunctions()
if selection.elementType not in self.numbering:
self.numbering[selection.elementType] = np.zeros((selection.elements.GetNumberOfElements(), nsf), dtype=MuscatIndex)-1
elementIdsConnectivity = selection.elements.connectivity[selection.indices, :]
for sf in range(nsf):
nn = GetNumberingColsSize(k)
on, idxI, idxII = sp.dofAttachments[sf]
key = self.GetKeyFromNameAndIdxI(on, selection.elementType, idxI)
if key not in tempAlmanac:
continue
(unique, newDofs, inverse) = tempAlmanac[key]
name, idxI, idxII = self.GetHashFor(selection.elements, sp, selection.indices, sf, False, elidsConnectivity=elementIdsConnectivity)
v = np.vstack((unique, idxI))
uniqueII, indices, inverse = np.unique(np.sort(v, axis=1), return_index=True, return_inverse=True, axis=0)
newnewdofs = np.hstack((newDofs, np.zeros(len(idxI), dtype=MuscatIndex)-1))[inverse][len(unique):]
self.numbering[selection.elementType][selection.indices, sf] = newnewdofs
Info("Numbering Done")
self._almanac = tempAlmanac
self.computeDofToPoint()
@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):
# print(self._almanac.keys())
self.pointDofs = np.zeros(self.totalNumberOfPoints, dtype=MuscatIndex)-1
if ('P', 'P') not in self._almanac:
return
unique, newdofs, inverse = self._almanac[('P', 'P')]
self._doftopointLeft = unique.flatten()
self._doftopointRight = newdofs.flatten()
self.pointDofs[self._doftopointLeft] = self._doftopointRight
[docs] def GetDofOfPoint(self, pointid):
return self.pointDofs[pointid]
@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):
raise
extractorLeftSide = np.empty(self.size, dtype=MuscatIndex)
extractorRightSide = np.empty(self.size, dtype=MuscatIndex)
tmpcpt = 0
for elem, data in self.mesh.elements.items():
unique, newdofs, inverse = self._almanac[elem]
extractorLeftSide[tmpcpt:tmpcpt+len(inverse)] = inverse
extractorRightSide[tmpcpt:tmpcpt+len(inverse)] = newdofs[inverse]
tmpcpt += len(newdofs)
# if k[0] is the elementname then k[1] is the connecivity
# we generate the same almanac with the number of each element
elemDic = {}
for name, data in self.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
globaloffsets = self.mesh.ComputeGlobalOffset()
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():
localdic = elemDic[k[0]]
if k[1] in localdic.keys():
globaloffset = globaloffsets[k[0]]
extractorLeftSide[tmpcpt] = globaloffset + localdic[k[1]]
extractorRightSide[tmpcpt] = v
tmpcpt += 1
self._doftocellLeft = np.resize(extractorLeftSide, (tmpcpt,))
self._doftocellRight = np.resize(extractorRightSide, (tmpcpt,))
[docs] def GetHashFor(self, data, sp, elids, sf, discontinuous=False, elidsConnectivity=None):
res = []
name = data.elementType
# Debug("2.1")
# print(type(elids))
if elidsConnectivity is None:
elidsConnectivity = data.connectivity[elids, :]
on, idxI, idxII = sp.dofAttachments[sf]
key = self.GetKeyFromNameAndIdxI(on, name, idxI)
if on == "P":
shapeFunctionConnectivity = elidsConnectivity[:, idxI]
return key, shapeFunctionConnectivity[:, None], idxII
if on == "C":
return key, elidsConnectivity, idxI
if on == "F":
edge = ED.faces1[name][idxI]
print(f"-------{name}-----------")
print(edge)
print(key)
print(elidsConnectivity)
return key, elidsConnectivity[:, edge[1]], 0
if on == "F2":
edge = ED.faces2[name][idxI]
return key, elidsConnectivity[:, edge[1]], 0
if on == "G":
"""G is for global """
return key, np.zeros((len(elids), 1), dtype=MuscatIndex), idxII
if on == "IP":
return key, np.hstack((elidsConnectivity, -idxI*np.ones(len(elids), dtype=MuscatIndex)[:, None])), idxI
return res
[docs]def CheckIntegrity(GUI: bool = False):
import Muscat.FE.DofNumbering as DN
return DN.CheckIntegrityUsingAlgo("NumpyBase", GUI)
if __name__ == '__main__':
print(CheckIntegrity(True)) # pragma: no cover