# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE', which is part of this source code package.
#
from __future__ import annotations
from Muscat.FE.Spaces.FESpaces import ConstantSpaceGlobal
from Muscat.FE.DofNumbering import ComputeDofNumbering
from Muscat.FE.Fields.FEField import FEField
from Muscat.FE.Fields.IPField import IPField
from Muscat.FE.SymWeakForm import GetField
from Muscat.FE.SymWeakForm import GetTestField
from Muscat.FE.Integration import IntegrateGeneral
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
[docs]def IPFieldIntegration(ipfield:IPField, ff:ElementFilter)->float:
"""Compte integral of a field over a region
Parameters
----------
ipfield : IPField
Field to be integrated
ff : ElementFilter
integrand (region)
Returns
-------
float
field integrated over region
"""
mesh = ipfield.mesh
space = ConstantSpaceGlobal
ipfield_name = ipfield.name
numbering = ComputeDofNumbering(mesh,space)
testfield = FEField("Testfunc",mesh=mesh,space=space,numbering=numbering)
Tt = GetTestField("Testfunc",1,sdim=mesh.GetPointsDimensionality())
rho = GetField(ipfield_name,1,sdim=mesh.GetPointsDimensionality())
wf = rho.T*Tt
constants = {}
fields = [ipfield]
unknownFields = [ testfield ]
_,F = IntegrateGeneral(mesh=mesh,
wform=wf,
constants=constants,
fields=fields,
unknownFields=unknownFields,
integrationRule=ipfield.rule,
elementFilter=ff)
res = F[0]
return res
[docs]def CheckIntegrity():
import Muscat.MeshContainers.ElementsDescription as ED
from Muscat.MeshTools.MeshCreationTools import CreateCube
mesh = CreateCube([2.,3.,4.],[-1.0,-1.0,-1.0],[2./10, 2./10,2./10])
from Muscat.FE.Fields.IPField import IPField
from Muscat.FE.IntegrationRules import GetIntegrationRuleByName
LP1 = GetIntegrationRuleByName("LagrangeP1Quadrature")
rho_field = IPField("rho",mesh=mesh,rule=LP1)
factor = 1
rho_field.Allocate(factor)
print("field for Hexaedron 6 elements, 8 integration points")
print(rho_field.GetDataFor(ED.Hexahedron_8))
from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter
ff = ElementFilter(dimensionality=3)
print("volume is " , IPFieldIntegration(rho_field,ff ) )
ff = ElementFilter(dimensionality=2)
print("surface is " , IPFieldIntegration(rho_field,ff ) )
return "ok"
if __name__ == "__main__":# pragma: no cover
print(CheckIntegrity())