# -*- 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 Optional
import numpy as np
from scipy.sparse import coo_matrix, csr_matrix
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.DofNumbering import ComputeDofNumbering
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.Types import MuscatIndex, MuscatFloat
from Muscat.FE.IntegrationRules import LagrangeIsoParamQuadrature
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceGeo
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Mesh import Mesh
[docs]
def PrepareFEComputation(mesh, elementFilter=None, numberOfComponents=None, space=LagrangeSpaceGeo, integrationRule=LagrangeIsoParamQuadrature):
"""
Prepares a finite element computation with Lagrange isoparametric finite
elements by computing FESpace, numberings, offset and NGauss
Parameters
----------
mesh : Mesh
mesh on which the function is applied
elementFilter : ElementFilter, optional
Filter on which the function is applied
numberOfComponents : int
the number of components
Returns
-------
FESpace
numberings
offset : list of int
NGauss : int
"""
dim = mesh.GetPointsDimensionality()
if elementFilter == None:
elementFilter = ElementFilter()
elementFilter.SetDimensionality(dim)
if numberOfComponents == None:
numberOfComponents = dim
NGauss = 0
for selection in elementFilter(mesh):
elementaryQuadrature = integrationRule[selection.elements.elementType]
NGauss += len(selection.indices)*elementaryQuadrature.GetNumberOfPoints()
numbering = ComputeDofNumbering(mesh, space, fromConnectivity=True)
numberings = [numbering]*numberOfComponents
offset = []
totaldofs = 0
for n in numberings:
offset.append(totaldofs)
totaldofs += n.size
return space, numberings, offset, NGauss
[docs]
def ComputeL2ScalarProductMatrix(mesh, numberOfComponents, elementFilter=None, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the L2 scalar product used to compute the correlations
between the primal solution snapshots. The numberOfComponents
depends on the solution type: 3 for solid mechanics in 3D, or 1 for
thermal in any dimension. (Lagrange isoparametric finite elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
numberOfComponents : int
the number of components
elementFilter : ElementFilter, optional
Filter on which the function is applied
Returns
-------
scipy.sparse.csr_matrix
the sparse matrix of the L2 scalar product
"""
from Muscat.FE.Integration import IntegrateGeneral
dim = mesh.GetPointsDimensionality()
if elementFilter == None:
elementFilter = ElementFilter()
elementFilter.SetDimensionality(dim)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, elementFilter, numberOfComponents, integrationRule=integrationRule)
from Muscat.FE.SymWeakForm import GetField
from Muscat.FE.SymWeakForm import GetTestField
T = GetField("T", numberOfComponents, sdim= mesh.GetPointsDimensionality())
Tt = GetTestField("T", numberOfComponents, sdim= mesh.GetPointsDimensionality())
wf = T.T*Tt
constants = {}
fields = []
if numberOfComponents == 1:
unknownFields = [FEField("T", mesh=mesh, space=spaces, numbering=numberings[0])]
else:
unknownFields = [FEField("T_"+str(i), mesh=mesh, space=spaces, numbering=numberings[i]) for i in range(numberOfComponents)]
K, _ = IntegrateGeneral(mesh=mesh, wform=wf, constants=constants, fields=fields,
unknownFields=unknownFields, elementFilter=elementFilter)
return K
[docs]
def ComputeH10ScalarProductMatrix(mesh, numberOfComponents, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the H10 scalar product used to compute the correlations
between the primal solution snapshots. The numberOfComponents
depends on the solution type: 3 for solid mechanics in 3D, or 1 for
thermal in any dimension (Lagrange isoparametric finite elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
numberOfComponents : int
the number of components of the primal variable snapshots
Returns
-------
scipy.sparse.csr_matrix
the sparse matrix of the H10 scalar product
"""
nbNodes = mesh.GetNumberOfNodes()
dim = mesh.GetPointsDimensionality()
ff = ElementFilter()
ff.SetDimensionality(dim)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ff, numberOfComponents, integrationRule=integrationRule)
ev = []
ei = []
ej = []
for selection in ff(mesh):
elementaryQuadrature = integrationRule[selection.elementType]
w = elementaryQuadrature.weights
lenNumbering = len(numberings[0][selection.elementType][0, :])
ones = np.ones(lenNumbering, dtype=int)
space_ipvalues = spaces[selection.elementType].GetSpaceEvaluatedAt(elementaryQuadrature.points, elementaryQuadrature.weights)
for el in selection.indices:
xcoor = mesh.nodes[selection.elements.connectivity[el], :]
leftNumberings = [numberings[j][selection.elementType][el, :]+offset[j] for j in range(numberOfComponents)]
for ip in range(elementaryQuadrature.GetNumberOfPoints()):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
BxByBzI = Jinv(space_ipvalues.valdphidxi[ip])
for j in range(numberOfComponents):
ev.extend((w[ip]*Jdet)*np.tensordot(BxByBzI, BxByBzI, axes=(0, 0)).ravel())
for i in leftNumberings[j]:
ei.extend(i*ones)
ej.extend(leftNumberings[j].ravel())
mat = coo_matrix((ev, (ei, ej)), shape=(numberOfComponents*nbNodes, numberOfComponents*nbNodes)).tocsr()
return mat
[docs]
def ComputeInterpolationMatrix_FE_GaussPoint(mesh, feSpace, integrationRule, feNumbering=None, ipNumbering=None, elementFilter=None):
from Muscat.FE.Integration import IntegrateGeneral
if elementFilter is None:
elementFilter = ElementFilter()
dim = mesh.GetPointsDimensionality()
elementFilter.SetDimensionality(dim)
else:
pass
if feNumbering is None:
numberingLeft = ComputeDofNumbering(mesh, space=feSpace, elementFilter=elementFilter)
else:
numberingLeft = feNumbering
leftField = FEField(name="P1", numbering=numberingLeft, mesh=mesh, space=feSpace)
from Muscat.FE.Spaces.IPSpaces import GenerateSpaceForIntegrationPointInterpolation
gaussSpace = GenerateSpaceForIntegrationPointInterpolation(integrationRule)
if ipNumbering is None:
numberingRight = ComputeDofNumbering(mesh, space=gaussSpace, elementFilter=elementFilter)
else:
numberingRight = ipNumbering
rightField = FEField(name="Gauss'", numbering=numberingRight, mesh=mesh, space=gaussSpace)
from Muscat.FE.SymWeakForm import GetField, GetTestField
LF = GetField("P1", 1, sdim=mesh.GetPointsDimensionality())
RF = GetTestField("Gauss", 1, sdim=mesh.GetPointsDimensionality())
symForm = LF.T*RF
interpMatrixMatrix, _ = IntegrateGeneral(mesh=mesh, constants={}, fields=[], wform=symForm, unknownFields=[leftField], testFields=[
rightField], onlyEvaluation=True, integrationRule=integrationRule, elementFilter=elementFilter)
return interpMatrixMatrix
[docs]
def ComputeJdetAtIntegPoint(mesh, elementSets=None, relativeDimension=0, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes determinant of the Jacobian of the transformation of the
transformation between the reference element and the mesh element, at
the integration points. (Lagrange isoparametric finite elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
elementSets : list if strings
sets of elements on which the det-Jacobians are computed
relativeDimension : int (0, -1 or -2)
difference between the dimension of the elements on which the function
is computed and the dimensionality of the mesh
Returns
-------
np.ndarray
of size (NGauss,)
"""
ff = ElementFilter()
dimension = mesh.GetPointsDimensionality() + relativeDimension
ff.SetDimensionality(dimension)
if elementSets:
assert type(elementSets) == list, "elementSets should be a list of elementSets"
ff.SetETags(elementSets)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ff, integrationRule=integrationRule)
jDet = np.zeros(NGauss)
countTotal = 0
for name, data, ids in ff(mesh):
elementaryQuadrature = integrationRule[name]
p, w = elementaryQuadrature.points, elementaryQuadrature.weights
NGaussperEl = elementaryQuadrature.GetNumberOfPoints()
space_ipvalues = spaces[name].GetSpaceEvaluatedAt(p, w)
for el in ids:
xcoor = mesh.nodes[data.connectivity[el], :]
for ip in range(NGaussperEl):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
jDet[countTotal] = Jdet
countTotal += 1
return jDet
[docs]
def ComputePhiAtIntegPoint(mesh, elementSets=None, relativeDimension=0, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the value of the finite element shape functions at the integration
points. (Lagrange isoparametric finite elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
elementSets : list if strings
sets of elements on which the det-Jacobian are computed
relativeDimension : int (0, -1 or -2)
difference between the dimension of the elements on which the function
is computed and the dimensionality of the mesh
Returns
-------
np.ndarray
of size (NGauss,)
coo_matrix
of size (NGauss, nbNodes)
"""
ff = ElementFilter()
dimension = mesh.GetPointsDimensionality() + relativeDimension
ff.SetDimensionality(dimension)
if elementSets:
assert type(elementSets) == list, "elementSets should be a list of elementSets"
ff.SetETags(elementSets)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ff, integrationRule=integrationRule)
nbNodes = mesh.GetNumberOfNodes()
phiAtIntegPointIndices = []
phiAtIntegPointValues = []
row = []
integrationWeights = np.zeros(NGauss)
countTotal = 0
for selection in ff(mesh):
elementaryQuadrature = integrationRule[selection.elementType]
p, w = elementaryQuadrature.points, elementaryQuadrature.weights
NGaussPerEl = elementaryQuadrature.GetNumberOfPoints()
lenNumbering = len(numberings[0][selection.elementType][0, :])
ones = np.ones(lenNumbering, dtype=int)
for el in selection.indices:
xcoor = mesh.nodes[selection.elements.connectivity[el], :]
space_ipvalues = spaces[selection.elementType].GetSpaceEvaluatedAt(p, w)
for ip in range(NGaussPerEl):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
integrationWeights[countTotal] = w[ip]*Jdet
leftNumbering = numberings[0][selection.elementType][el, :]
left = space_ipvalues.valN[ip]
phiAtIntegPointIndices.extend(leftNumbering)
phiAtIntegPointValues.extend(left)
row.extend(countTotal*ones)
countTotal += 1
phiAtIntegPointMatrix = coo_matrix((phiAtIntegPointValues, (row, phiAtIntegPointIndices)), shape=(NGauss, nbNodes))
return integrationWeights, phiAtIntegPointMatrix
[docs]
def ComputePhiAtIntegPointFromElFilter(mesh, elFilter, space=LagrangeSpaceGeo, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the shape functions on the integration
points and the integration weights associated with the integration scheme
Parameters
----------
mesh : Mesh
elFilter : elementFilter
Returns
-------
np.ndarray
integrationWeights
coo_matrix of size (NGauss, nbNodes)
phiAtIntegPointMatrix
"""
numbering = ComputeDofNumbering(mesh, space, fromConnectivity=True)
nbNodes = mesh.GetNumberOfNodes()
phiAtIntegPointIndices = []
phiAtIntegPointValues = []
row = []
integrationWeights = []
countTotal = 0
for name, data, ids in elFilter(mesh):
elementaryQuadrature = integrationRule[name]
p, w = elementaryQuadrature.points, elementaryQuadrature.weights
NGaussPerEl = len(w)
lenNumbering = len(numbering[name][0, :])
ones = np.ones(lenNumbering, dtype=int)
for el in ids:
xcoor = mesh.nodes[data.connectivity[el], :]
space_ipvalues = space[name].GetSpaceEvaluatedAt(p, w)
for ip in range(NGaussPerEl):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
integrationWeights.append(w[ip]*Jdet)
leftNumbering = numbering[name][el, :]
left = space_ipvalues.valN[ip]
phiAtIntegPointIndices.extend(leftNumbering)
phiAtIntegPointValues.extend(left)
row.extend(countTotal*ones)
countTotal += 1
NGauss = len(integrationWeights)
integrationWeights = np.array(integrationWeights)
phiAtIntegPointMatrix = coo_matrix((phiAtIntegPointValues, (row, phiAtIntegPointIndices)), shape=(NGauss, nbNodes))
return integrationWeights, phiAtIntegPointMatrix
[docs]
def ComputeGradPhiAtIntegPoint(mesh, elementSets=None, relativeDimension=0, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the components of the gradient of the shape functions on the
integration points and the integration weights associated with the
integration scheme
Parameters
----------
mesh : Mesh
mesh on which the function is applied
elementSets : list of strings
element sets, support for the shape functions
relativeDimension : int
0, -1, or -2: the dimension of the element set relative to the dimension of the mesh
Returns
-------
np.ndarray
integrationWeights
list (length dimensionality of the mesh) coo_matrix of size (NGauss, nbNodes)
gradPhiAtIntegPoint
"""
ff = ElementFilter()
dimension = mesh.GetPointsDimensionality() + relativeDimension
ff.SetDimensionality(dimension)
if elementSets:
assert type(elementSets) == list, "elementSets should be a list of elementSets"
ff.SetETags(elementSets)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ff, integrationRule=integrationRule)
nbNodes = mesh.GetNumberOfNodes()
meshDimension = mesh.GetPointsDimensionality()
GradPhiAtIntegPointIndices = []
GradPhiAtIntegPointValues = [[] for i in range(meshDimension)]
row = []
integrationWeights = np.zeros(NGauss)
countTotal = 0
for selection in ff(mesh):
elementaryRule = integrationRule[selection.elementType]
p = elementaryRule.points
w = elementaryRule.weights
NGaussPerEl = len(w)
lenNumbering = len(numberings[0][selection.elementType][0, :])
ones = np.ones(lenNumbering, dtype=int)
for el in selection.indices:
xcoor = mesh.nodes[selection.elements.connectivity[el], :]
space_ipvalues = spaces[selection.elementType].GetSpaceEvaluatedAt(p, w)
for ip in range(NGaussPerEl):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
BxByBzI = Jinv(space_ipvalues.valdphidxi[ip])
integrationWeights[countTotal] = w[ip]*Jdet
leftNumbering = numberings[0][selection.elementType][el, :]
GradPhiAtIntegPointIndices.extend(leftNumbering)
for k in range(meshDimension):
GradPhiAtIntegPointValues[k].extend(BxByBzI[k])
row.extend(countTotal*ones)
countTotal += 1
GradPhiAtIntegPointMatrix = [coo_matrix((GradPhiAtIntegPointValues[k], (row, GradPhiAtIntegPointIndices)), shape=(NGauss, nbNodes)) for k in range(meshDimension)]
return integrationWeights, GradPhiAtIntegPointMatrix
[docs]
def ComputeNormalsAtPoints(mesh: Mesh, ef: Optional[ElementFilter] = None) -> np.ndarray:
"""Compute the normal at the nodes selected by the element filter
Parameters
----------
mesh : Mesh
mesh to use, must have 2D element for a 3D mesh, or 1D elements for a 2D mesh
ef : ElementFilter, optional
Selection of elements to compute the normal on , by default None
Returns
-------
np.ndarray
a (number of nodes, dimensionality) array with the normal at the nodes (zero filled)
"""
pDim = mesh.GetPointsDimensionality()
edim = mesh.GetElementsDimensionality()#
if ef is None:
ef = ElementFilter(dimensionality = pDim-1)
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.FE.Fields.FieldTools import VectorToFEFieldsData
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.IntegrationRules import GetIntegrationRuleByName
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ef, integrationRule=GetIntegrationRuleByName("NodalEvaluationIsoParam"))
N = [FEField(name="N_"+str(x), mesh=mesh,space=spaces, numbering= numberings[0]) for x in range(pDim)]
from Muscat.FE.SymWeakForm import GetNormal, GetField, GetTestField
normal = GetNormal(pDim)
NTsym = GetTestField("N", pDim, sdim=pDim)
normalAtNodes = normal.T*NTsym
_, F = IntegrateGeneral(mesh=mesh, wform=normalAtNodes, constants={}, fields=[],
unknownFields=N,
elementFilter=ef)
res = np.zeros( (mesh.GetNumberOfNodes(), pDim) , dtype=MuscatFloat)
F2 = np.reshape(F.view(), (F.size//pDim,pDim) , 'F')
norm = np.linalg.norm(F2,axis=1)[:, None]
norm[norm==0] = 1
F2 /= norm
res[numberings[0].doftopointLeft,:] = F2[numberings[0].doftopointRight,:]
return res
[docs]
def ComputeNormalsAtIntegPoint(mesh, elementSets, integrationRule=LagrangeIsoParamQuadrature):
"""
Computes the normals at the elements from the sets elementSets in the
ambiant space (i.e. if mesh is of dimensionality 3, elementSets should
contain 2D elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
elementSets : list of strings
sets of elements on which the function is computed
Returns
-------
np.ndarray
of size (dimensionality, NGauss)
"""
ff = ElementFilter()
dimension = mesh.GetPointsDimensionality() - 1
ff.SetDimensionality(dimension)
if elementSets:
assert type(elementSets) == list, "elementSets should be a list of elementSets"
for elementSet in elementSets:
if elementSet:
ff.AddETag(elementSet)
spaces, numberings, offset, NGauss = PrepareFEComputation(mesh, ff, integrationRule=integrationRule)
normalsAtIntegPoint = np.empty((mesh.GetPointsDimensionality(), NGauss))
countTotal = 0
for selection in ff(mesh):
elementaryQuadrature = integrationRule[selection.elementType]
p, w = elementaryQuadrature.points, elementaryQuadrature.weights
NGaussperEl = elementaryQuadrature.GetNumberOfPoints()
for el in selection.indices:
xcoor = mesh.nodes[selection.elements.connectivity[el], :]
space_ipvalues = spaces[selection.elementType].GetSpaceEvaluatedAt(p, w)
for ip in range(NGaussperEl):
Jacobian, Jdet, Jinv = space_ipvalues.GetJacobianDeterminantAndInverseAtIP(ip, xcoor)
normal = space_ipvalues.GetNormal(Jacobian)
normalsAtIntegPoint[:, countTotal] = normal
countTotal += 1
return normalsAtIntegPoint
[docs]
def CellDataToIntegrationPointsData(mesh, scalarFields, elementSet=None, relativeDimension=0, integrationRule=LagrangeIsoParamQuadrature):
"""
Change the representation of scalarFields from data constant by cell
(elements) to data at integration points. (Lagrange isoparametric finite
elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
scalarFields : np.ndarray of size (nbe of fields, nbe of elements) or dict
with "nbe of fields" as keys and np.ndarray of size (nbe of elements,) or floats
as values fields whose representation is changed by the function
elementSet : elementSet defining the elements on which the function is
computed. If None, takes all the elements of considered dimension
relativeDimension : int (0, -1 or -2)
difference between the dimension of the elements on which the function
is computed and the dimensionality of the mesh
Returns
-------
np.ndarray
of size (nbe of fields, numberOfIntegrationPoints)
"""
dimension = mesh.GetPointsDimensionality() + relativeDimension
ff = ElementFilter()
ff.SetDimensionality(dimension)
if elementSet != None:
ff.AddETag(elementSet)
_, _, _, NGauss = PrepareFEComputation(mesh, ff, dimension, integrationRule=integrationRule)
if type(scalarFields) == dict:
keymap = list(scalarFields.keys())
elif type(scalarFields) == np.ndarray:
keymap = np.arange(scalarFields.shape[0])
else:
raise RuntimeError("scalarFields is neither an np.ndarray nor a dict")
numberOfFields = len(keymap)
typeField = []
for f in range(numberOfFields):
if type(scalarFields[keymap[f]]) == float:
typeField.append("scalar")
else:
typeField.append("vector")
integrationPointsData = np.empty((numberOfFields, NGauss))
countEl = 0
countIp = 0
for selection in ff(mesh):
elementaryQuadrature = integrationRule[selection.elementType]
numberOfIPperEl = elementaryQuadrature.GetNumberOfPoints()
for el in selection.indices:
for f in range(numberOfFields):
if typeField[f] == "vector":
integrationPointsData[f, countIp:countIp+numberOfIPperEl] = scalarFields[keymap[f]][countEl]
else:
integrationPointsData[f, countIp:countIp+numberOfIPperEl] = scalarFields[keymap[f]]
countEl += 1
countIp += numberOfIPperEl
return integrationPointsData
[docs]
def IntegrationPointsToCellData(mesh, scalarFields, integrationRule=LagrangeIsoParamQuadrature):
"""
Change the representation of scalarFields from data at integration points
to data constant by cell (elements) - taking the mean of values at the
integration points in each elements. (Lagrange isoparametric finite
elements)
Parameters
----------
mesh : Mesh
mesh on which the function is applied
scalarFields : np.ndarray of size (nbe of fields, nbe of elements) or dict
with "nbe of fields" as keys and np.ndarray of size (nbe of elements,) as
values fields whose representation in changed by the function
Returns
-------
np.ndarray
of size (nbe of elements,)
"""
if type(scalarFields) == dict:
keymap = list(scalarFields.keys())
elif type(scalarFields) == np.ndarray:
keymap = np.arange(scalarFields.shape[0])
else:
raise RuntimeError("scalarFields is neither an np.ndarray nor a dict")
numberOfFields = len(keymap)
from Muscat.FE.Fields import IPField as IPF
cellData = []
ef = ElementFilter(dimensionality=mesh.GetPointsDimensionality())
for f in range(numberOfFields):
iPField = IPF.IPField("", mesh, rule=integrationRule)
iPField.Allocate()
iPField.SetDataFromNumpy(scalarFields[keymap[f]], ef)
cellData.append(iPField.GetCellRepresentation())
return np.array(cellData)
[docs]
def CheckIntegrity(GUI: bool = False):
from Muscat.FE.SymPhysics import MechPhysics
mechPhysics = MechPhysics()
wform = mechPhysics.GetBulkFormulation(1.0, 0.3)
res, rhs = GetElementaryMatrixForFormulation(ED.Hexahedron_8, wform, unknownNames=mechPhysics.GetPrimalNames())
import Muscat.TestData as MuscatTestData
from Muscat.IO import GeofReader as GR
mesh = GR.ReadGeof(MuscatTestData.GetTestDataPath()+"cube2.geof")
ComputeL2ScalarProductMatrix(mesh, 1)
ComputeL2ScalarProductMatrix(mesh, 3)
ComputeH10ScalarProductMatrix(mesh, 1)
ComputeH10ScalarProductMatrix(mesh, 3)
ComputePhiAtIntegPoint(mesh)
ComputeGradPhiAtIntegPoint(mesh)
ComputeIntegrationPointsTags(mesh, 3)
ComputeNormalsAtIntegPoint(mesh, ["x0"])
length = len(mesh.elements[ED.Quadrangle_4].tags["x0"].GetIds())
fields = {'a': np.ones(length), 'b': 1.}
CellDataToIntegrationPointsData(mesh, fields, "x0", relativeDimension=-1)
res = CellDataToIntegrationPointsData(mesh, np.ones((2, 216)))
nGauss = res.shape[1]
IntegrationPointsToCellData(mesh, np.ones((2, nGauss)))
from Muscat.FE.IntegrationRules import LagrangeP1Quadrature
mesh.elements[ED.Hexahedron_8].tags.CreateTag("Transfer").SetIds([0, 1])
elementFilter = ElementFilter(eTag="Transfer")
data = ComputeInterpolationMatrix_FE_GaussPoint(mesh=mesh, feSpace=LagrangeSpaceP1, integrationRule=LagrangeP1Quadrature, elementFilter=elementFilter)
print(np.dot(data.toarray(), np.arange(12)))
from Muscat.IO.StlReader import ReadStl
mesh = ReadStl(MuscatTestData.GetTestDataPath()+"stlOpenSphere.stl")
normals = ComputeNormalsAtPoints(mesh, ElementFilter(dimensionality=2) )
if GUI:
from Muscat.Actions.OpenInParaView import OpenInParaView
mesh.nodeFields["CalculatedNormals"] = normals
OpenInParaView(mesh)
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(GUI=True))
print("Done")