Commit 79e985a8 authored by Felipe Bordeu's avatar Felipe Bordeu
Browse files

(SymPhysics.py) Clean refactor

parent cb1bc6b6

src/BasicTools/FE/SymPhysics.py

+307 −225
Original line number Diff line number Diff line
@@ -12,12 +12,16 @@ from sympy.matrices import Matrix 
from BasicTools.NumpyDefs import PBasicFloatType

from BasicTools.Containers.Filters import ElementFilter

from BasicTools.Helpers.BaseOutputObject import BaseOutputObject as BOO

from BasicTools.FE.SymWeakForm import Gradient,Divergence, GetField,GetTestField, GetScalarField, Inner

from BasicTools.FE.SymWeakForm import Gradient, Divergence, GetField, GetTestField, GetScalarField, Inner, Trace

import BasicTools.FE.SymWeakForm as swf

import BasicTools.Helpers.ParserHelper as PH

from BasicTools.Helpers.BaseOutputObject import froze_it





@froze_it

class Physics(BOO):

    """Basic class to hold the information about symbolic terms"""



    def __init__(self):

        self.integrationRule = None

        self.spaces = [None]
@@ -26,6 +30,10 @@ class Physics(BOO): 
        self.numberings = None

        self.spaceDimension = 3

        self.extraRHSTerms = []

        self.coeffs = {}



    def GetCoeff(self, var_name):

        return self.coeffs.get(var_name, GetScalarField(var_name))



    def AddToRHSTerm(self, nf, val):

        """ nf  : a NodalFilter
@@ -50,10 +58,10 @@ class Physics(BOO): 
        return u.T*ut*a



    def SetSpaceToLagrange(self, P=None, isoParam=None):

        if P is None and isoParam is None:

        if P is None and isoParam is None:  # pragma: no cover

            raise (ValueError("Please set the type of integration P=1,2 or isoParam=True"))



        if P is not None and isoParam is not None:

        if P is not None and isoParam is not None:  # pragma: no cover

            raise (ValueError("Please set the type of integration P=1,2 or isoParam=True"))



        if isoParam is None or isoParam == False:
@@ -65,8 +73,8 @@ class Physics(BOO): 
                from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceP2

                space = LagrangeSpaceP2

                self.integrationRule = "LagrangeP2"

            else:

                raise(ValueError("I don't understand"))

            else:  # pragma: no cover

                raise (ValueError(f" P cant be : {P} , must be 1, 2 or None"))

        else:

            from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceGeo

            space = LagrangeSpaceGeo
@@ -85,7 +93,6 @@ class Physics(BOO): 


    def AddLFormulation(self, zoneOrElementFilter, data):





        if type(zoneOrElementFilter) == str:

            ef = ElementFilter(tag=zoneOrElementFilter)

        else:
@@ -118,14 +125,17 @@ class Physics(BOO): 
                    break

            else:

                if fromConnectivity:

                    self.numberings[d] = ComputeDofNumbering(mesh,self.spaces[d],fromConnectivity = True ,dofs=self.numberings[d])

                    self.numberings[d] = ComputeDofNumbering(

                        mesh, self.spaces[d], fromConnectivity=True, dofs=self.numberings[d])

                else:

                    self.numberings[d] = ComputeDofNumbering(mesh,self.spaces[d],fromConnectivity =False,elementFilter=allFilters,dofs=self.numberings[d])

                    self.numberings[d] = ComputeDofNumbering(

                        mesh, self.spaces[d], fromConnectivity=False, elementFilter=allFilters, dofs=self.numberings[d])



    def ComputeDofNumberingFromConnectivity(self, mesh):

        self.ComputeDofNumbering(mesh, fromConnectivity=True)





@froze_it

class MecaPhysics(Physics):

    """Weak forms for mechanical problems
@@ -137,26 +147,25 @@ class MecaPhysics(Physics): 
        The type of model used, by default "isotropic"

        option are : "isotropic", "orthotropic", "anisotropic"

    """



    def __init__(self, dim=3, elasticModel="isotropic"):

        super(MecaPhysics,self).__init__()

        self.dim = dim

        super().__init__()

        self.UnFrozen()

        self.spaceDimension = dim



        self.mecaPrimalName = ("u",self.dim)

        self.mecaPrimalName = ("u", self.spaceDimension)

        self.SetMecaPrimalName(self.mecaPrimalName[0])



        self.mecaSpace = None



        self.coeffs = defaultdict(lambda : 0.)

        self.coeffs["young"] = 1.

        self.coeffs["poisson"] = 0.3

        self.coeffs["density"] = 1.



        if elasticModel not in ["isotropic", "orthotropic", "anisotropic"]:  # pragma: no cover

            raise Exception(f'elasticModel ({elasticModel}) not available : options are "isotropic", "orthotropic", "anisotropic" ')



            raise Exception(

                f'elasticModel ({elasticModel}) not available : options are "isotropic", "orthotropic", "anisotropic" ')

        self.elasticModel = elasticModel



        self.materialOrientations = np.eye(dim, dtype = PBasicFloatType)

        self.HookeLocalOperator = None



        self.materialOrientations = np.eye(self.spaceDimension, dtype=PBasicFloatType)

        self.planeStress = True



    def SetYoung(self, val):
@@ -166,7 +175,7 @@ class MecaPhysics(Physics): 
        self.coeffs["poisson"] = PH.ReadFloat(val)



    def SetMecaPrimalName(self, name):

        self.mecaPrimalName = (name,self.dim)

        self.mecaPrimalName = (name, self.spaceDimension)

        self.primalUnknown = GetField(self.mecaPrimalName[0], self.mecaPrimalName[1])

        self.primalTest = GetTestField(self.mecaPrimalName[0], self.mecaPrimalName[1])
@@ -190,28 +199,38 @@ class MecaPhysics(Physics): 
        from BasicTools.FE.MaterialHelp import HookeLaw



        if young is None:

            young = self.coeffs["young"]

            young = self.GetCoeff("young")

        if poisson is None:

            poisson = self.coeffs["poisson"]

            poisson = self.GetCoeff("poisson")



        young = GetScalarField(young)*GetScalarField(factor)



        op = HookeLaw()

        op.Read({"E": young, "nu": poisson})

        return op.HookeIso(dim=self.dim,planeStress=self.planeStress)



        return op.HookeIso(dim=self.spaceDimension, planeStress=self.planeStress)



    def GetHookeOperatorOrthotropic(self, factor=None):



        if self.spaceDimension == 2:

            hookOrtho = [["C11", "C12",    0],

                         ["C12", "C22",    0],

                         [0,     0, "C33"]]

        elif self.spaceDimension == 3:

            hookOrtho = [["C11", "C12", "C13",     0,     0,    0],

                         ["C12", "C22", "C23",     0,     0,    0],

                         ["C13", "C23", "C33",     0,     0,    0],

                         [0,     0,     0, "C44",     0,    0],

                         [0,     0,     0,     0, "C55",    0],

                         [0,     0,     0,     0,     0, "C66"]]

        else:  # pragma: no cover

            raise Exception("Orthotropic material no available for dimension : {self.spaceDimension}")



        return  np.array([[ GetScalarField(self.coeffs[c]) for c in line] for line in hookOrtho]) *GetScalarField(factor)

        return np.array([[self.GetCoeff(c) for c in line] for line in hookOrtho]) * GetScalarField(factor)



    def GetHookeOperatorAnisotropic(self, factor=None):

        if self.spaceDimension != 3:  # pragma: no cover

            raise Exception("Anisotropic material no available for dimension : {self.spaceDimension}")



        hookAnisotropic = [["C1111", "C1122", "C1133", "C1123", "C1131", "C1112"],

                           ["C2211", "C2222", "C2233", "C2223", "C2231", "C2212"],

                           ["C3311", "C3322", "C3333", "C3323", "C3331", "C3312"],
@@ -219,8 +238,7 @@ class MecaPhysics(Physics): 
                           ["C3111", "C3122", "C3133", "C3123", "C3131", "C3112"],

                           ["C1211", "C1222", "C1233", "C1223", "C1231", "C1212"]]



        return  np.array([[ GetScalarField(self.coeffs[c]) for c in line] for line in hookAnisotropic]) *GetScalarField(factor)



        return np.array([[self.GetCoeff(c) for c in line] for line in hookAnisotropic]) * GetScalarField(factor)



    def GetStressVoigt(self, utGlobal, HookeLocalOperator=None):

        from BasicTools.FE.SymWeakForm import ToVoigtEpsilon, Strain
@@ -228,7 +246,7 @@ class MecaPhysics(Physics): 
            HookeLocalOperator = self.HookeLocalOperator



        uLocal = Inner(self.materialOrientations, utGlobal)

        return swf.Inner(ToVoigtEpsilon(Strain(uLocal,self.dim)).T,HookeLocalOperator)

        return swf.Inner(ToVoigtEpsilon(Strain(uLocal, self.spaceDimension)).T, HookeLocalOperator)



    def GetBulkFormulation(self, young=None, poisson=None, alpha=None):

        from BasicTools.FE.SymWeakForm import ToVoigtEpsilon, Strain
@@ -240,7 +258,7 @@ class MecaPhysics(Physics): 
        HookeLocalOperator = self.GetHookeOperator(young, poisson, alpha)

        stress = self.GetStressVoigt(uGlobal, HookeLocalOperator)



        return  stress*ToVoigtEpsilon(Strain(utLocal,self.dim))

        return stress*ToVoigtEpsilon(Strain(utLocal, self.spaceDimension))



    def GetPressureFormulation(self, pressure):

        ut = self.primalTest
@@ -248,7 +266,7 @@ class MecaPhysics(Physics): 
        p = GetScalarField(pressure)



        from BasicTools.FE.SymWeakForm import GetNormal

        Normal = GetNormal(self.dim)

        Normal = GetNormal(self.spaceDimension)



        return p*Normal.T*ut
@@ -269,7 +287,7 @@ class MecaPhysics(Physics): 
    def GetAccelerationFormulation(self, direction, density=None):



        if density is None:

            density = self.coeffs["density"]

            density = self.GetCoeff("density")



        ut = self.primalTest

        density = GetScalarField(density)
@@ -280,11 +298,11 @@ class MecaPhysics(Physics): 


        return density*direction.T*ut



    def GetCentrifugalTerm(self, axe=[0,0,1], pointOnAxe=[0,0,0], angularSpeed=1, density=1):

    def GetCentrifugalTerm(self, axe=[0, 0, 1], pointOnAxe=[0, 0, 0], angularSpeed=1, density=None):

        ut = self.primalTest



        if density is None:

            density = self.coeffs["density"]

            density = self.GetCoeff("density")



        density = GetScalarField(density)
@@ -296,7 +314,7 @@ class MecaPhysics(Physics): 
            pointOnAxe = [GetScalarField(d) for d in pointOnAxe]

            pointOnAxe = Matrix([pointOnAxe]).T



        pos = GetField("Pos",self.dim) - pointOnAxe

        pos = GetField("Pos", self.spaceDimension) - pointOnAxe

        r = pos - pos.dot(axe)*pos



        return -density*angularSpeed**2*r.T*ut
@@ -309,14 +327,15 @@ class MecaPhysics(Physics): 
        utLocal = Inner(self.materialOrientations, uGlobal)



        nodalEnergyT = GetTestField("elastic_energy", 1)

        symEnergy = 0.5*wf.ToVoigtEpsilon(wf.Strain(utLocal)).T*self.HookeLocalOperator*wf.ToVoigtEpsilon(wf.Strain(utLocal))*nodalEnergyT



        symEnergy = 0.5*wf.ToVoigtEpsilon(wf.Strain(utLocal)).T*self.HookeLocalOperator * \

            wf.ToVoigtEpsilon(wf.Strain(utLocal))*nodalEnergyT



        trStrainT = GetTestField("tr_strain_", 1)

        symTrStrain = wf.Trace(wf.Strain(uGlobal))*trStrainT



        trStressT = GetTestField("tr_stress_", 1)

        symTrStress = wf.Trace(wf.FromVoigtSigma(wf.ToVoigtEpsilon(wf.Strain(utLocal)).T*self.HookeLocalOperator))*trStressT

        symTrStress = wf.Trace(wf.FromVoigtSigma(wf.ToVoigtEpsilon(

            wf.Strain(utLocal)).T*self.HookeLocalOperator))*trStressT



        postQuantities = {"elastic_energy": symEnergy,

                          "tr_strain_": symTrStrain,
@@ -327,7 +346,7 @@ class MecaPhysics(Physics): 


class MecaPhysicsAxi(MecaPhysics):

    def __init__(self):

        super(MecaPhysics,self).__init__(dim=2)

        super().__init__(dim=2)



    def GetFieldR(self):

        return GetScalarField("r")
@@ -339,22 +358,23 @@ class MecaPhysicsAxi(MecaPhysics): 
        ut = self.primalTest



        if young is None:

            young = self.young

            young = self.GetCoeff("young")

        if poisson is None:

            poisson = self.poisson

            poisson = self.GetCoeff("poisson")



        young = GetScalarField(young)*GetScalarField(alpha)



        op = HookeLaw()

        op.Read({"E": young, "nu": poisson})

        self.HookeLocalOperator = op.HookeIso(dim=self.dim,planeStress=self.planeStress, axisymetric=self.axisymetric)

        self.HookeLocalOperator = op.HookeIso(dim=self.spaceDimension, planeStress=self.planeStress, axisymetric=True)

        print(self.HookeLocalOperator)



        r = self.GetFieldR()



        from BasicTools.FE.SymWeakForm import StrainAxyCol

        epsilon_u = StrainAxyCol(u, r)

        epsilon_ut = StrainAxyCol(ut, r)

        return 2*np.pi*epsilon_u*self.HookeLocalOperator*epsilon_ut*r

        return 2*np.pi*epsilon_u.T*self.HookeLocalOperator*epsilon_ut*r



    def GetPressureFormulation(self, pressure):

        return super().GetPressureFormulation(pressure)*self.GetFieldR()
@@ -370,20 +390,23 @@ class MecaPhysicsAxi(MecaPhysics): 
        import BasicTools.FE.SymWeakForm as wf

        symDep = self.primalUnknown



        pir2 = 2*np.pi*self.GetFieldR()

        r = self.GetFieldR()

        pir2 = 2*np.pi*r



        nodalEnergyT = GetTestField("strain_energy", 1)

        symEnergy = pir2*0.5*wf.ToVoigtEpsilon(wf.Strain(symDep)).T*self.HookeLocalOperator*wf.ToVoigtEpsilon(wf.Strain(symDep))*nodalEnergyT



        symEnergy = pir2*0.5*wf.StrainAxyCol(symDep, r).T*self.HookeLocalOperator * \

            wf.StrainAxyCol(symDep, r)*nodalEnergyT



        from sympy import prod



        trStrainT = GetTestField("tr(strain)", 1)

        strain = wf.StrainAxyCol(symDep)

        symTrStrain = prod(strain[0:3]).T*trStrainT

        strain = wf.StrainAxyCol(symDep, r)

        symTrStrain = prod(strain[0:3])*trStrainT



        trStressT = GetTestField("tr(stress)", 1)

        stress = strain*self.HookeLocalOperator

        symTrStress = prod(stress[0:3]).T*trStressT

        stress = strain.T*self.HookeLocalOperator

        symTrStress = prod(stress[0:3])*trStressT



        postQuantities = {"strain_energy": symEnergy,

                          "tr(strain)": symTrStrain,
@@ -391,9 +414,12 @@ class MecaPhysicsAxi(MecaPhysics): 


        return postQuantities





@froze_it

class BasicPhysics(Physics):

    def __init__(self):

        super(BasicPhysics,self).__init__()

        super().__init__()

        self.UnFrozen()

        self.PrimalNameTrial = ("u", 1)

        self.PrimalNameTest = ("u", 1)

        self.Space = None
@@ -407,14 +433,12 @@ class BasicPhysics(Physics): 
        self.primalUnknown = GetField(self.PrimalNameTrial[0], self.PrimalNameTrial[1])

        self.primalTest = GetTestField(self.PrimalNameTest[0], self.PrimalNameTest[1])





    def GetPrimalNames(self):

        return [self.PrimalNameTrial[0]]



    def GetPrimalDims(self):

        return [self.PrimalNameTrial[1]]





    def GetBulkFormulation_dudi_dtdj(self, u=0, t=0, i=0, j=0, alpha=1.):



        a = GetScalarField(alpha)
@@ -447,9 +471,12 @@ class BasicPhysics(Physics): 


        return f*tt





@froze_it

class ThermalPhysics(Physics):

    def __init__(self):

        super(ThermalPhysics,self).__init__()

        super().__init__()

        self.UnFrozen()

        self.thermalPrimalName = ("t", 1)

        self.SetPrimalNames(self.thermalPrimalName)

        self.thermalSpace = None
@@ -501,28 +528,27 @@ class ThermalPhysics(Physics): 
        Tinf = GetScalarField(Tinf)

        return beta*(t - Tinf)*tt





@froze_it

class StokesPhysics(Physics):

    def __init__(self):

        super(StokesPhysics,self).__init__()

        self.velocityPrimalName = ("v",3)

        self.pressurePrimalName = ("p",1)

        self.SetPrimalNames(self.velocityPrimalName[0])



        super().__init__()

        self.UnFrozen()

        self.SetPrimalNames("v", "p")



    def GetPrimalNames(self):

        res = [self.velocityPrimalName[0] + "_" + str(c) for c in range(self.velocityPrimalName[1])]

        res.append(self.pressurePrimalName)

        return res



    def SetMecaPrimalName(self, vName, pName):

        self.velocityPrimalName = (vName,self.dim)

    def SetPrimalNames(self, vName, pName):

        self.velocityPrimalName = (vName, self.spaceDimension)

        self.pressurePrimalName = (pName, 1)

        self.primalUnknownV = GetField(self.velocityPrimalName[0], self.velocityPrimalName[1])

        self.primalUnknownP = GetField(self.pressurePrimalName[0], self.pressurePrimalName[1])

        self.primalTestV = GetTestField(self.velocityPrimalName[0], self.velocityPrimalName[1])

        self.primalTestP = GetTestField(self.pressurePrimalName[0], self.pressurePrimalName[1])





    def SetSpaceToLagrange(self, P=None, isoParam=None):

        from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceP1

        from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceP2
@@ -531,22 +557,24 @@ class StokesPhysics(Physics): 
        self.spaces.append(LagrangeSpaceP1)

        self.integrationRule = "LagrangeP2"





    def GetBulkFormulation(self, mu=1.):



        mu = GetScalarField(mu)

        v = self.primalUnknownV

        vt = self.primalTestV

        p  = self.primalUnknownP

        pt = self.primalTestP



        res = Gradient(v,self.spaceDimension).T*mu*Gradient(vt,self.spaceDimension)  -  Divergence(vt,self.spaceDimension)*p + pt*Divergence(v,self.spaceDimension)

        p = self.primalUnknownP[0, 0]

        pt = self.primalTestP[0, 0]

        a = Trace(Gradient(v, self.spaceDimension)*mu*Gradient(vt, self.spaceDimension))

        b = Divergence(vt, self.spaceDimension)*p

        c = pt*Divergence(v, self.spaceDimension)

        return a - b + c



        return res



@froze_it

class ThermoMecaPhysics(Physics):

    def __init__(self):

        super(ThermoMecaPhysics,self).__init__()

        super().__init__()

        self.UnFrozen()

        self.mecaPhys = MecaPhysics()

        self.thermalPhys = ThermalPhysics()
@@ -567,20 +595,49 @@ class ThermoMecaPhysics(Physics): 




def CheckIntegrity(GUI=False):



    bp = BasicPhysics()

    bp.Reset()

    bp.GetBulkMassFormulation()

    bp.SetSpaceToLagrange(1)

    bp.SetSpaceToLagrange(2)

    bp.SetSpaceToLagrange(isoParam=True)

    bp.AddBFormulation("tagname", bp.GetBulkMassFormulation())



    from BasicTools.Containers.Filters import ElementFilter

    bp.AddBFormulation(ElementFilter(), bp.GetBulkMassFormulation())

    bp.AddLFormulation('ElementTagName', bp.GetBulkMassFormulation())

    bp.AddLFormulation(ElementFilter(), bp.GetBulkMassFormulation())

    print(bp.GetNumberOfUnkownFields())



    assert (bp.ExpandNames(("t", 1)) == "t")



    bp.spaceDimension = 3

    bp.SetPrimalName("U", "V", 3, 3)

    print(bp.primalUnknown)

    print(bp.primalTest)

    print(bp.GetBulkFormulation_dudi_dtdj(u=0, i=1, t=1, j=2))

    print(bp.GetBulkFormulation_dudi_dtdj())

    print(bp.GetBulkLaplacian())

    print(bp.GetFlux())



    ###############################################################

    res = ThermoMecaPhysics()

    print(res.GetBulkFormulation())

    print(res.GetPrimalNames())



    print(BasicPhysics().GetBulkFormulation_dudi_dtdj())

    t = BasicPhysics()

    t.spaceDimension = 3

    t.SetPrimalName("U","V",3,3)

    print(t.primalUnknown)

    print(t.primalTest)

    print(t.GetBulkFormulation_dudi_dtdj(u=0,i=1,t=1,j=2) )



    ###############################################################

    M2D = MecaPhysics(dim=2)

    print(M2D.GetHookeOperator())

    M2D.SetYoung(1)

    M2D.SetPoisson(0)

    np.testing.assert_allclose(np.array(M2D.GetHookeOperator()), [[1, 0, 0], [0, 1, 0], [0, 0, 0.5]])

    ###############################################################



    M2D = MecaPhysics(dim=2, elasticModel="orthotropic")

    print(np.array(M2D.GetHookeOperator()))



    ###############################################################



    M3DI = MecaPhysics(dim=3)

    from itertools import product

    print(M3DI.GetHookeOperator())
@@ -598,9 +655,34 @@ def CheckIntegrity(GUI=False): 
    print(M3DA.GetPressureFormulation(1))

    print(M3DA.GetAccelerationFormulation([1, 0, 0]))

    print(M3DA.GetDistributedForceFormulation([1, 0, 0]))





    print(M3DA.GetForceFormulation([1, 0, 0]))

    print(M3DA.GetCentrifugalTerm())



    print(M3DA.PostTreatmentFormulations())



    ###############################################################

    MPA = MecaPhysicsAxi()

    MPA.GetBulkFormulation()

    MPA.GetPressureFormulation(1)

    MPA.GetForceFormulation([1, 0], 1)

    MPA.GetAccelerationFormulation([1, 0], 1)

    MPA.PostTreatmentFormulations()

    ###############################################################



    TP = ThermalPhysics()

    TP.SetThermalPrimalName("u")

    TP.GetBulkFormulation([1, 0, 0])

    TP.GetMassOperator()

    TP.GetNormalFlux()

    TP.GetRobinFlux(beta=1, Tinf=24)

    ###############################################################



    SP = StokesPhysics()

    SP.GetPrimalNames()

    SP.SetSpaceToLagrange()

    SP.GetBulkFormulation()

    return "ok"





if __name__ == '__main__':

    print(CheckIntegrity(GUI=True))  # pragma: no cover