Source code for Muscat.FE.IntegrationTools

# -*- 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.Containers.Mesh import Mesh
from Muscat.Containers.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
from Muscat.FE.IntegrationRules import LagrangeIsoParamQuadrature


[docs]def IntegrateSymExpression(mesh:Mesh, fields: List[FieldBase], symexp, ef:Optional[ElementFilter]=None ) -> MuscatFloat:

    space = ConstantSpaceGlobal
    numbering = ComputeDofNumbering(mesh,space)
    testfield = FEField("Testfunc",mesh=mesh,space=space,numbering=numbering)
    Tt = GetTestField("Testfunc",1)
    integrationRule = LagrangeIsoParamQuadrature

    for f in fields:
        if isinstance(f, IPField):
            integrationRule = f.rule

    _,f = IntegrateGeneral(mesh, symexp*Tt, {}, fields, unknownFields = [ testfield ],integrationRule=integrationRule, elementFilter=ef )
    return f[0]




[docs]def IntegrateField(field, ef:Optional[ElementFilter]=None ):
    return IntegrateSymExpression(field.mesh, [field], symexp= GetField(field.name, size=1), ef=ef)



[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.Containers.Filters.FilterObjects import ElementFilter

    field = CreateFieldFromDescription(mesh,[(ElementFilter(),1)] )
    area = IntegrateField(field)
    print(f"Computed Area = {area}")
    import numpy as np
    np.testing.assert_equal(area, 2.9088166085800675)

    field = CreateFieldFromDescription(mesh,[(ElementFilter(),1)],fieldType="IP" )
    area = IntegrateField(field)
    print(f"Computed Area = {area}")
    import numpy as np
    np.testing.assert_equal(area, 2.9088166085800675)

    return "ok"



if __name__ == '__main__':
    print(CheckIntegrity(True))  # pragma: no cover