# -*- 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.
#
from typing import Tuple, Optional, Union, List, Dict
import numpy as np
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
import Muscat.Containers.ElementsDescription as ED
from Muscat.Containers.Filters.FilterObjects import ElementFilter
from Muscat.Containers.Mesh import Mesh
from Muscat.Containers.MeshCreationTools import QuadToLin
from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
from Muscat.LinAlg.Transform import Transform
from Muscat.FE.Fields.FEField import FEField
# Tolerance to stop the newton, dxietaphi.squaredNorm()
D_XI_ETA_PHI_TOL_SQUARED = 1e-10
# Tolerance to detect on projection
XI_ETA_PHI_TOL_SQUARED = 1e-10
#Relative tolerance to detect inside outside (if the f in the newton is zero) with respect to the diagonal
# of the bounding box. l**2/diagonal**2 */
REL_X_Y_Z_TOL_SQUARED = 1e-10
# Maximal number of newton iters
MAX_NEWTON_ITER = 25
[docs]def GetFieldTransferOpPython(inputField: FEField, targetPoints: ArrayLike, method: Union[str, None] = None, verbose: bool = False, elementFilter: Optional[ElementFilter] = None) -> Tuple[np.ndarray, np.ndarray]:
"""Compute the transfer operator from the inputField to the target points
valueAtTargetPoints = op.dot(FEField.data)
The methods available are:
- "Interp/Nearest"
- "Nearest/Nearest"
- "Interp/Clamp"
- "Interp/Extrap"
- "Interp/ZeroFill"
Possible values (int) for the returned status are:
- 0: "Nearest"
- 1: "Interp"
- 2: "Extrap"
- 3: "Clamp"
- 4: "ZeroFill"
- 5: "MultipleExtrapPoint"
Parameters
----------
inputField : FEField
the FEField to be transferred
targetPoints : ArrayLike
Numpy array of the target points. Position to extract the values
method : Union[str,None], optional
A couple for the algorithm used when the point is inside/outside the mesh
- "Interp" -> to use the interpolation of the FEField to extract the values
- "Nearest" -> to use the closest node to extract the values
- "Clamp" -> to use the closest point to the surface if the point is outside
- "Extrap" -> to use shape function of the closest element to compute a extrapolation
- "ZeroFill" -> fill with zero in the case the point is outside
If None is provided then "Interp/Clamp" is used
verbose : bool, optional
Print a progress bar, by default False
elementFilter : Optional[ElementFilter], optional
ElementFilter to extract the information from only a subdomain, by default None
Returns
-------
Tuple [np.ndarray,np.ndarray]
a tuple with 3 object containing:
- op, sparse matrix with the operator to make the transfer
- status: vector of ints with the status transfer for each target point:
- entitiesIds: vector of ints with the node or element number of the closes entitie (in function of status)
return op, status, entitiesIds
"""
possibleMethods = ["Interp/Nearest", "Nearest/Nearest", "Interp/Clamp", "Interp/Extrap", "Interp/ZeroFill"]
possibleMethodsDict = {"Nearest": 0, "Interp": 1, "Extrap": 2, "Clamp": 3, "ZeroFill": 4}
from scipy.spatial import cKDTree
from scipy.sparse import coo_matrix
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
if method is None:
method = possibleMethods[2]
elif method not in possibleMethods:
raise (Exception(f"Method for transfer operator not know '{method}', possible options are : {possibleMethods}"))
insideMethod = possibleMethodsDict[method.split("/")[0]]
outsideMethod = possibleMethodsDict[method.split("/")[1]]
iMeshDim = inputField.mesh.GetElementsDimensionality()
originalMesh = inputField.mesh
if elementFilter == None:
elementFilter = ElementFilter(dimensionality=iMeshDim)
from Muscat.Containers.MeshModificationTools import CleanLonelyNodes
iMesh = ExtractElementsByElementFilter(originalMesh, elementFilter)
CleanLonelyNodes(iMesh)
numbering = inputField.numbering
space = inputField.space
iNodes = iMesh.nodes
nbTargetPoints = targetPoints.shape[0]
kdt = cKDTree(iNodes)
if insideMethod == 0 and outsideMethod == 0:
dist, ids = kdt.query(targetPoints)
ids = iMesh.originalIDNodes[ids]
if numbering is None or numbering.fromConnectivity:
cols = ids
else:
cols = [numbering.GetDofOfPoint(pid) for pid in ids]
row = np.arange(nbTargetPoints)
data = np.ones(nbTargetPoints)
return coo_matrix((data, (row, cols)), shape=(nbTargetPoints, iNodes.shape[0])), np.zeros(nbTargetPoints), ids
# be sure to create the spaces
for elementType, data in inputField.mesh.elements.items():
space[elementType].Create()
LagrangeSpaceGeo[elementType].Create()
facesToTreat = [ED.faces1[elementType],
ED.faces2[elementType],
ED.faces3[elementType]]
for tt in facesToTreat:
for elementType2, num in tt:
LagrangeSpaceGeo[elementType2].Create()
# we build de Dual Connectivity
from Muscat.Containers.MeshInspectionTools import ComputeNodeToElementConnectivity
dualGraph, nbUsed = ComputeNodeToElementConnectivity(iMesh)
from Muscat.Containers.MeshTools import GetElementsCenters
centers = GetElementsCenters(iMesh)
kdtCenters = cKDTree(centers)
# 30 to be sure to hold exa27 coefficients
cols = np.empty(nbTargetPoints*30, dtype=int)
rows = np.empty(nbTargetPoints*30, dtype=int)
dataS = np.empty(nbTargetPoints*30)
fillCpt = 0
cooData = [cols, rows, dataS, fillCpt]
def AddToOutput(l, col, row, dat, cooData):
fillCpt = cooData[3]
s = np.s_[fillCpt: fillCpt + l]
cooData[0][s] = col
cooData[1][s] = row
cooData[2][s] = dat
cooData[3] += l
def GetElement(iMesh, enb):
for name, data in iMesh.elements.items():
if enb < data.GetNumberOfElements():
return originalMesh.elements[name], data.originalIds[enb], data, enb
else:
enb -= data.GetNumberOfElements()
raise (Exception("Element not found"))
distTP, idsTP = kdt.query(targetPoints)
distTPcenters, idsTPcenters = kdtCenters.query(targetPoints)
if verbose:
from Muscat.Helpers.ProgressBar import PrintProgressBar
PrintProgressBar(0, nbTargetPoints, prefix=f'Building Transfer {method}:', suffix='Complete', length=50)
verboseCpt = 0
status = np.zeros(nbTargetPoints, dtype=MuscatIndex)
entities = np.zeros(nbTargetPoints, dtype=MuscatIndex)
ones = np.ones(50)
for p in range(nbTargetPoints):
if verbose:
nvc = int(p/nbTargetPoints*1000)
if verboseCpt != nvc:
PrintProgressBar(p+1, nbTargetPoints, prefix=f'Building Transfer {method}:', suffix='Complete', length=50)
verboseCpt = nvc
TP = targetPoints[p, :] # target point position
CP = idsTP[p] # closest point position
# we use the closest element (in the sens of cells center )
original_data, lenb, imesh_data, imesh_elnb = GetElement(iMesh, idsTPcenters[p])
# construct the potentialElements list (all element touching the closest element)
potentialElements = []
# Element connected to the closest point
potentialElements.extend(dualGraph[idsTP[p], 0:nbUsed[idsTP[p]]])
# Elements connected to the closest element (bases on the element center)
for elempoint in imesh_data.connectivity[imesh_elnb, :]:
potentialElements.extend(dualGraph[elempoint, 0:nbUsed[elempoint]])
potentialElements = np.unique(potentialElements)
# compute distance to elements
# for the moment we use the distance to the center, this gives a good estimate
# of the order to check the elements
dist = np.empty(len(potentialElements))
for cpt, e in enumerate(potentialElements):
#data, lenb, imesh_data, imesh_elnb = GetElement(iMesh,e)
diff = centers[e, :]-TP
dist[cpt] = diff.dot(diff)
# order the element to test, closest element first
index = np.argsort(dist)
dist = dist[index]
potentialElements = potentialElements[index]
minDistToElement = 1e10
minDistToElementProjection = 1e10
multiple_closest_elements = 0 # to clamp in the case of
#print(p, potentialElements)
for cpt, e in enumerate(potentialElements):
original_data, lenb, imesh_data, imesh_elnb = GetElement(iMesh, e)
localnumbering = numbering[original_data.elementType]
localspace = space[original_data.elementType]
posnumbering = original_data.connectivity
#posspace = LagrangeSpaceGeo[original_data.elementType]
coordAtDofs = originalMesh.nodes[posnumbering[lenb, :], :]
inside, distv, bary, baryClamped, toProjectionVector, masterDim, elemDim = ComputeInterpolationExtrapolationsBarycentricCoordinates(TP, original_data.elementType, coordAtDofs, LagrangeSpaceGeo)
# update the distance**2 with a *exact* distance
dist[cpt] = distv.dot(distv)
# compute shape function of the incomming space using the xi eta phi
shapeFunc = localspace.GetShapeFunc(bary)
shapeFuncClamped = localspace.GetShapeFunc(baryClamped)
# need to add a tolerance over the distv (real distance). this is needed because
# we can have a point that the projection is inside an element (surface or line)
# but not on the surface/line. Need to add a better test
if inside : #and normsquared(distv) <= 1e-14
sF = shapeFunc
status[p] = 1
entities[p] = e
memMasterDim = masterDim
memElemDim = elemDim
break
# store the best element (closest)
projectionDistance = normsquared(toProjectionVector)
if abs(dist[cpt] - minDistToElement) < 1e-14:
if projectionDistance < minDistToElementProjection:
minDistToElementProjection = projectionDistance
minDistToElement = dist[cpt]
memshapeFunc = shapeFunc
memshapeFuncClamped = shapeFuncClamped
memlenb = lenb
memlocalnumbering = localnumbering
memElement = e
memMasterDim = masterDim
memElemDim = elemDim
if multiple_closest_elements ==0:
multiple_closest_elements += 1
multiple_closest_elements += 1
elif dist[cpt] <= minDistToElement and inside is not None:
multiple_closest_elements = 0
minDistToElement = dist[cpt]
if inside is None:
# non is error in the inverse mapping, only clamp is available
memshapeFunc = None
else:
memshapeFunc = shapeFunc
memshapeFuncClamped = shapeFuncClamped
memlenb = lenb
memlocalnumbering = localnumbering
memElement = e
memMasterDim = masterDim
memElemDim = elemDim
#print("keep " , memElement)
else:
# we are outside
if minDistToElement == 1e10 or outsideMethod == 0 or (outsideMethod == 2 and memshapeFunc is None):
# not a single element found
col = [CP]
row = [p]
dat = [1.]
status[p] = 0
entities[p] = CP
AddToOutput(len(col), col, row, dat, cooData)
continue
if outsideMethod == 2 and memshapeFunc is not None:
entities[p] = memElement
if multiple_closest_elements+memElemDim > memMasterDim:
sF = memshapeFuncClamped
status[p] = -3
else:
sF = memshapeFunc
status[p] = -2
status[p] = memElemDim
elif outsideMethod == 3:
sF = memshapeFuncClamped
status[p] = 3
entities[p] = memElement
else:
status[p] = 4
entities[p] = memElement
continue
lenb = memlenb
localnumbering = memlocalnumbering
col = localnumbering[lenb, :]
l = len(col)
row = p*ones[0:l]
dat = sF
AddToOutput(l, col, row, dat, cooData)
if verbose:
PrintProgressBar(nbTargetPoints, nbTargetPoints, prefix=f'Building Transfer {method}:', suffix='Complete', length=50)
s = np.s_[0:cooData[3]]
return coo_matrix((cooData[2][s], (cooData[1][s], cooData[0][s])), shape=(nbTargetPoints, inputField.numbering.size)), status, entities
[docs]def ComputeInterpolationCoefficients(mask: ArrayLike, TP: ArrayLike, elementsdata, coordAtDofs: ArrayLike, posspace, localspace):
"""Compute the interpolation coefficients for the element.
Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
"""
ft = True
for cpt, (faceElementType, faceLocalConnectivity) in enumerate(elementsdata):
# work only on element touching the mask
if not np.any(mask[faceLocalConnectivity]):
continue
faceInside, onProjection, fbary, fbaryClamped, distToElementExtended, distToElement = ComputeBarycentricCoordinateOnElement(coordAtDofs[faceLocalConnectivity, :], posspace[faceElementType], TP, faceElementType)
faceDistv = posspace[faceElementType].GetShapeFunc(fbary).dot(coordAtDofs[faceLocalConnectivity, :]) - TP
faceDist = faceDistv.dot(faceDistv)
if faceInside is not None and onProjection:
if ft:
ft = False
faceDistvMem = faceDistv
faceDistMem = faceDist
fbaryMem = fbary
faceLocalConnectivityMem = faceLocalConnectivity
faceElementTypeMem = faceElementType
else:
if faceDist < faceDistMem:
faceDistMem = faceDist
faceDistvMem = faceDistv
fbaryMem = fbary
faceLocalConnectivityMem = faceLocalConnectivity
faceElementTypeMem = faceElementType
if ft:
return False, False, None, None
else:
flocalspace = posspace[faceElementTypeMem]
fshapeFunc = flocalspace.GetShapeFunc(fbaryMem)
baryClamped = fshapeFunc.dot(localspace.posN[faceLocalConnectivityMem, :])
return False, True, faceDistvMem, baryClamped
[docs]def ComputeShapeFunctionsOnElement(coordAtDofs: ArrayLike, localspace, localnumbering, point: ArrayLike, faceElementType):
inside, onProjection, xiEtaPhi, xiEtaPhiClamped, distToElementExtended, distToElement = ComputeBarycentricCoordinateOnElement(coordAtDofs, localspace, point, faceElementType)
N = localspace.GetShapeFunc(xiEtaPhi)
NClamped = localspace.GetShapeFunc(xiEtaPhiClamped)
return inside, N, NClamped
[docs]def ddf(f, xiEtaPhi, dN, GetShapeFuncDerDer, coordAtDofs, linear, x):
"""Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
"""
dNX = dN.dot(coordAtDofs)
res = dNX.dot(dNX.T)
# After some investigation the Newton is more stable using only the first part
# of the hessian, so we activate the second part only after 10 iteration.
# the number 10 is arbitrary and can be increased
if linear or x < 10:
return res
ddN = GetShapeFuncDerDer(xiEtaPhi)
# we work for every coordinate
for cpt in range(coordAtDofs.shape[1]):
# error of the component cpt
fx = f[cpt]
coordx = coordAtDofs[:, cpt]
# we build the derivative of the pos field with respect of xi chi
# d2f_i/dxidchi * pos[i,ccpt]
pp = [ddNi.dot(xi) for ddNi, xi in zip(ddN, coordx)]
p = np.add.reduce(pp)
res -= fx*p
return res
[docs]def df(f, dN, coordAtDofs):
""" Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
"""
dNX = dN.dot(coordAtDofs)
res = -f.dot(dNX.T)
return res
[docs]def vdet(A: ArrayLike) -> MuscatFloat:
"""Compute the determinant of a 3x3 matrix
Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
Parameters
----------
A : ArrayLike
a 3x3 matrix
Returns
-------
MuscatFloat
the determinant
"""
detA = A[0, 0] * (A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]) -\
A[0, 1] * (A[2, 2] * A[1, 0] - A[2, 0] * A[1, 2]) +\
A[0, 2] * (A[1, 0] * A[2, 1] - A[2, 0] * A[1, 1])
return detA
[docs]def hdinv(A: ArrayLike) -> np.ndarray:
"""Compute the inverse of a 3x3 matrix
Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
Parameters
----------
A : ArrayLike
a 3x3 matrix
Returns
-------
np.ndarray
the inverse of A
"""
invA = np.zeros_like(A)
detA = vdet(A)
invA[0, 0] = (-A[1, 2] * A[2, 1] +
A[1, 1] * A[2, 2]) / detA
invA[1, 0] = (A[1, 2] * A[2, 0] -
A[1, 0] * A[2, 2]) / detA
invA[2, 0] = (-A[1, 1] * A[2, 0] +
A[1, 0] * A[2, 1]) / detA
invA[0, 1] = (A[0, 2] * A[2, 1] -
A[0, 1] * A[2, 2]) / detA
invA[1, 1] = (-A[0, 2] * A[2, 0] +
A[0, 0] * A[2, 2]) / detA
invA[2, 1] = (A[0, 1] * A[2, 0] -
A[0, 0] * A[2, 1]) / detA
invA[0, 2] = (-A[0, 2] * A[1, 1] +
A[0, 1] * A[1, 2]) / detA
invA[1, 2] = (A[0, 2] * A[1, 0] -
A[0, 0] * A[1, 2]) / detA
invA[2, 2] = (-A[0, 1] * A[1, 0] +
A[0, 0] * A[1, 1]) / detA
return invA
[docs]def inv22(A: ArrayLike) -> np.ndarray:
"""Compute the inverse of a 2x2 matrix
Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
Parameters
----------
A : ArrayLike
a 2x2 matrix
Returns
-------
np.ndarray
the inverse of A
"""
a = A[0, 0]
b = A[0, 1]
c = A[1, 0]
d = A[1, 1]
invA = np.zeros_like(A)
invA[0, 0] = d
invA[0, 1] = -b
invA[1, 0] = -c
invA[1, 1] = a
invA /= (a*d-b*c)
return invA
[docs]def ComputeBarycentricCoordinateOnElement(coordAtDofs: ArrayLike, localspace, targetPoint: ArrayLike, elementType: ED.ElementTypeLike):
"""Newton to compute the best baricentric coordinates on an element for the target point
Warning!! This function is not intended for the final user. function used by (GetFieldTransferOp)
Parameters
----------
coordAtDofs : ArrayLike
Coordinates of the node of the element
localspace : _type_
_description_
targetPoint : ArrayLike
Target Point
elementType : str
Element type
Returns
-------
_type_
_description_
"""
linear = ED.linear[elementType]
spacedim = localspace.GetDimensionality()
# print("spacedim",spacedim)
xietaphi = np.array([0.5]*spacedim)
N = localspace.GetShapeFunc(xietaphi)
currentPoint = N.dot(coordAtDofs)
f = currentPoint - targetPoint
diag = np.max(coordAtDofs,axis=0) - np.min(coordAtDofs,axis=0)
diag2 = diag.dot(diag)
for x in range(MAX_NEWTON_ITER):
#print("----------- in iteration ",x)
dN = localspace.GetShapeFuncDer(xietaphi)
#print(" dN ", dN)
#print(" f ", f)
#print(" coordAtDofs ", coordAtDofs)
df_num = df(-f, dN, coordAtDofs)
#print(" df_num ", df_num)
H = ddf(-f, xietaphi, dN, localspace.GetShapeFuncDerDer, coordAtDofs, linear, x)
#print(" H ", H)
if spacedim == 2:
dxietaphi = inv22(H).dot(df_num)
elif spacedim == 3:
dxietaphi = hdinv(H).dot(df_num)
else:
dxietaphi = df_num/H[0, 0]
#print(" dxietaphi ", dxietaphi)
xietaphi -= dxietaphi
# if the cell is linear only one iteration is needed
N = localspace.GetShapeFunc(xietaphi)
f = N.dot(coordAtDofs) - targetPoint
if linear:
break
if normsquared(dxietaphi) < D_XI_ETA_PHI_TOL_SQUARED:
break
else:
return None, False, xietaphi, localspace.ClampParamCoordinates(xietaphi), None, None
xichietaClamped = localspace.ClampParamCoordinates(xietaphi)
distToElement = localspace.GetShapeFunc(xichietaClamped).dot(coordAtDofs) - targetPoint
# we treat very closes point as inside
inside = normsquared(distToElement) < REL_X_Y_Z_TOL_SQUARED*diag2
#inside = normsquared(xichietaClamped-xietaphi) < XI_ETA_PHI_TOL_SQUARED
#print(" xietaphi ", xietaphi)
onProjection = normsquared(xichietaClamped- xietaphi) < XI_ETA_PHI_TOL_SQUARED
return inside, onProjection, xietaphi, xichietaClamped, f, distToElement
[docs]def normsquared(v: ArrayLike) -> np.number:
"""sqared norm of a vector
Parameters
----------
v : ArrayLike
a vector
Returns
-------
np.number
the squared norm
"""
return sum([x*x for x in v])
[docs]def CheckIntegrity(GUI: bool = False):
return "OK"
if __name__ == '__main__':
print(CheckIntegrity(True)) # pragma: no cover