# -*- 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, List
from Muscat.Types import MuscatFloat
from Muscat.MeshContainers.Mesh import Mesh
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.FE.Fields.FieldTools import IntegrationPointWrapper
from Muscat.FE.Fields.FieldBase import FieldBase
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Fields.IPField import IPField
from Muscat.FE.Spaces.FESpaces import ConstantSpaceGlobal
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.FE.SymWeakForm import GetTestField, GetField, GetNormal, GetCoordinate
from Muscat.FE.IntegrationRules import LagrangeIsoParamQuadrature
[docs]
def IntegrateSymExpression(mesh:Mesh, symexp, fields: List[FieldBase]=[], elementFilter:Optional[ElementFilter]=None ) -> MuscatFloat:
space = ConstantSpaceGlobal
numbering = ComputeDofNumbering(mesh,space)
testfield = FEField("Testfunc",mesh=mesh,space=space,numbering=numbering)
Tt = GetTestField("Testfunc",1, sdim=mesh.GetPointsDimensionality())
integrationRule = LagrangeIsoParamQuadrature
for f in fields:
if isinstance(f, IPField):
integrationRule = f.rule
_,f = IntegrateGeneral(mesh, symexp*Tt, {}, fields, unknownFields = [ testfield ],integrationRule=integrationRule, elementFilter=elementFilter )
return f[0]
[docs]
def IntegrateField(field, elementFilter:Optional[ElementFilter]=None ):
return IntegrateSymExpression(field.mesh, symexp= GetField(field.name, size=1,sdim=field.mesh.GetPointsDimensionality()), fields= [field], elementFilter=elementFilter)
[docs]
def CheckIntegrity(GUI= False):
from Muscat.IO.StlReader import ReadStl
from Muscat.TestData import GetTestDataPath
mesh = ReadStl(GetTestDataPath() + "stlsphere.stl")
from Muscat.FE.Fields.FieldTools import CreateFieldFromDescription
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
import numpy as np
field = CreateFieldFromDescription(mesh,[(ElementFilter(),1.)] )
field.name = "rho"
area = IntegrateField(field)
np.testing.assert_allclose(area, 2.9088166085800675,2e-14)
field = CreateFieldFromDescription(mesh, [(ElementFilter(),1)], fieldType="IP" )
field.name = "rho"
area = IntegrateField(field)
np.testing.assert_allclose(area, 2.908816608580072,2e-14)
return CheckIntegrityRPNIntegrator(GUI)
[docs]
def CheckIntegrityRPNIntegrator(GUI= False):
from Muscat.FE.Integration import UseKokkos
if UseKokkos == False:
print("This test is only valid for the RPN Kokkos integrator")
return "skip"
from Muscat.IO.StlReader import ReadStl
from Muscat.TestData import GetTestDataPath
mesh = ReadStl(GetTestDataPath() + "stlsphere.stl")
import numpy as np
from sympy import sin, pi, cos, tan, Abs, sqrt
coords = GetCoordinate(3)
symnormals = GetNormal(3)
exprToTests = [ (coords[0], lambda pos: pos[:,0]),
(coords[1], lambda pos: pos[:,1]),
(coords[2], lambda pos: pos[:,2]),
( cos(0), lambda pos: np.ones(pos.shape[0])*np.cos(0)),
( cos(0.8)+1, lambda pos: np.ones(pos.shape[0])*np.cos(0.8)+1),
( sin(0), lambda pos: np.ones(pos.shape[0])*np.sin(0)),
( cos(coords[1])+1.2, lambda pos: np.cos(pos[:,1])+1.2),
( tan(coords[2]+0.1), lambda pos: np.tan(pos[:,2]+0.1)),
( Abs(coords[2]-0.2) , lambda pos: np.abs(pos[:,2]-0.2)),
( sqrt(coords[2]**2+0.2) , lambda pos: np.sqrt(pos[:,2]**2+0.2)),
( Abs(coords[2]-0.2)/sqrt(coords[2]**2+0.2)+0.1 , lambda pos: np.abs(pos[:,2]-0.2)/np.sqrt(pos[:,2]**2+0.2)+0.1 ),
( symnormals[0] , lambda pos: numnormals[:,0] ),
]
from Muscat.FE.IntegrationRules import elementCenterQuadrature
from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP0
from Muscat.MeshTools.MeshTools import GetElementsCenters
from Muscat.MeshTools.MeshInspectionTools import ComputeNormalsAtElements
centers = GetElementsCenters(mesh)
numnormals = ComputeNormalsAtElements(mesh)
for symExpr, lambdaFunction in exprToTests:
print(f"-------- working on {symExpr} ------------")
testfield = FEField("Testfunc",mesh=mesh, space= LagrangeSpaceP0 )
Tt = GetTestField("Testfunc",1,sdim=3)
_,f_I = IntegrateGeneral(mesh, symExpr*Tt, constants = {}, fields = [], unknownFields = [ testfield ], integrationRule=elementCenterQuadrature, elementFilter=ElementFilter(), onlyEvaluation=True )
f_II = lambdaFunction(centers)
np.testing.assert_allclose(f_I, f_II)
return "ok"
if __name__ == '__main__':
print(CheckIntegrity(True)) # pragma: no cover