Source code for Muscat.FE.WeakForms.NumericalWeakForm

# -*- 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.
#


import numpy as np
from sympy.core.containers import Tuple
from sympy import Symbol, Function

from Muscat.FE.SymWeakForm import GetNormal


UseCpp = False
# Include the correct implementation depending on the UseCpp value
try:
    UseCpp = True
    # CPP integration
    from Muscat.FE.WeakForms.NativeNumericalWeakForm import (
        PyWeakTerm,
        PyWeakMonom,
        PyWeakForm,
    )
except:
    UseCpp = False
    print("Warning NativeNumericalWeakForm (cpp) not available, using python variant")
    # PYTHON integration
    from Muscat.FE.WeakForms.PythonNumericalWeakFrom import (
        PyWeakTerm,
        PyWeakMonom,
        PyWeakForm,
    )




def SymWeakMonomToNumWeakMono(exp):
    from sympy.core.mul import Mul
    from sympy.core.power import Pow

    if exp.func == Mul:
        res = PyWeakMonom()
        for arg in exp.args:
            if arg.is_Number:
                res.prefactor = float(arg)
                continue

            if isinstance(arg, Pow):
                term = ConverTermToProd(arg.args[0])
                if term is None:
                    print(type(arg.args[0]))
                    print(arg.args[0])
                    raise (
                        Exception("Unable to treat term " + str(arg.args[0])))
                for i in range(arg.args[1]):
                    res.AddProd(term)
                continue

            term = ConverTermToProd(arg)
            if term is not None:
                # res.prod.append(term)
                res.AddProd(term)
                continue

            print(type(arg))
            print(arg)

            raise
        return res
    else:
        # only one term no product
        res = PyWeakMonom()
        term = ConverTermToProd(exp)
        if term is not None:
            res.AddProd(term)
        return res



__cached_normal = GetNormal(3)




def ConverTermToProd(arg):
    if isinstance(arg, Symbol):
        t = PyWeakTerm()
        t.constant = True
        t.derDegree = 0
        t.fieldName = str(arg)
        return t

    from sympy.core.function import Derivative

    if type(arg) == Derivative:
        t = PyWeakTerm()

        t.fieldName = str(arg.args[0].func)
        # Python 3
        t.derDegree = 1
        if Tuple == type(arg.args[1]):
            # sympy 1.2
            t.derCoordName = str(arg.args[1][0])
        else:
            # sympy 1.1.1
            t.derCoordName = str(arg.args[1])

        sn = []
        for i in range(0, len(arg.args[0].args)):
            sn.append(str(arg.args[0].args[i]))

        t.derCoordIndex_ = sn.index(t.derCoordName)
        return t

    if isinstance(arg, Function):
        t = PyWeakTerm()
        # print(str(arg.func)+"* * * * * * ")
        # print(arg.func)
        # print (N)
        # print ([arg == nc for nc in N])
        if np.any([arg == nc for nc in __cached_normal]):
            t.normal = True
            t.derDegree = int(str(arg.func).split("_")[1])
        else:
            t.derDegree = 0

        t.fieldName = str(arg.func)
        return t

    raise Exception(f"Error!! arg of type ({type(arg)}) not supported: {arg}")





def SymWeakToNumWeak(exp):
    from sympy.core.add import Add
    from sympy.core.mul import Mul

    exp = exp.expand()
    res = PyWeakForm()
    try:
        if exp.shape[0] == 1 and exp.shape[1] == 1:
            exp = exp[0, 0]
    except:
        pass

    if exp.func == Mul:
        # res.AddTermform.append(SymWeakMonomToNumWeakMono(exp))
        res.AddTerm(SymWeakMonomToNumWeakMono(exp))

    elif exp.func == Add:
        for monom in exp.args:
            # res.form.append(SymWeakMonomToNumWeakMono(monom))
            res.AddTerm(SymWeakMonomToNumWeakMono(monom))

    else:
        mono = PyWeakMonom()
        mono.AddProd(ConverTermToProd(exp))
        res.AddTerm(mono)
        # raise(Exception("error treating formulation term"))

    return res





def CheckIntegrity(GUI: bool = False):
    return "ok"



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