# -*- 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 Muscat.FE.SymWeakForm import GetTestSufixChar, space3D, GetField
from sympy import Expr, Function, Symbol
from Muscat.FE.WeakForms.Operators import NumericalOperators
from typing import List, Tuple, Dict
[docs]
class RPNElement:
def __init__(self, name=None, value=None, type=NumericalOperators.NONE):
self.name = name
self.value = value
self.type = type
self.der1 = -1
self.der2 = -1
def __str__(self):
return f"RPNElement({self.name}, {self.value}, {self.type},{self.der1},{self.der2}) "
def __repr__(self):
d = {
# NumericalOperators.DERIVATIVE: "d/d",
NumericalOperators.POW: "^",
NumericalOperators.PRODUCT: "*",
NumericalOperators.ADDITION: "+",
}
if self.type in d.keys():
return d[self.type]
if self.name is not None:
return str(self.name)
elif self.value is not None:
return str(self.value)
else:
return str(self.type)
[docs]
def RPNToSymWeak(rempl, exp):
from Muscat.FE.SymWeakForm import Symbol, GetField
from sympy.core.add import Add
from sympy.core.mul import Mul
from sympy.core.power import Pow
from sympy.functions.elementary.trigonometric import cos
from sympy.functions.elementary.trigonometric import sin
from sympy.functions.elementary.trigonometric import tan
from sympy.functions.elementary.exponential import log
from sympy.functions.elementary.exponential import exp as exponential
from sympy.functions.elementary.complexes import Abs
BinaryOperatorMapping = {
NumericalOperators.ADDITION: Add,
NumericalOperators.PRODUCT: Mul,
NumericalOperators.POW: Pow,
# NumericalOperators.DERIVATIVE: Derivative,
}
UnaryOperatorMapping = {
NumericalOperators.COS : cos,
NumericalOperators.SIN : sin,
NumericalOperators.TAN : tan,
NumericalOperators.LOG : log,
NumericalOperators.EXP : exponential,
NumericalOperators.ABS : Abs,
}
if rempl is None:
rempl = []
remplNames = [x[0] for x in rempl]
stack = []
for i in exp:
if i.type == NumericalOperators.COORDINATE:
stack.append(Symbol(i.name))
elif i.type == NumericalOperators.SCALAR:
if i.name in remplNames:
index = remplNames.index(i.name)
stack.append(RPNToSymWeak(rempl[0:index], rempl[index][1]))
else:
stack.append(Symbol(i.name))
elif i.type == NumericalOperators.NUMERICAL:
stack.append(i.value)
elif i.type in [
NumericalOperators.UNKNOWN_FIELD,
NumericalOperators.KNOWN_FE_FIELD,
NumericalOperators.TEST_FIELD,
]:
stack.append(GetField(i.name, size=1,sdim=3)[0])
elif i.type in BinaryOperatorMapping.keys():
b = stack.pop()
a = stack.pop()
stack.append(BinaryOperatorMapping[i.type](a, b))
elif i.type in UnaryOperatorMapping.keys():
a = stack.pop()
stack.append(UnaryOperatorMapping[i.type](a))
else:
raise Exception(i)
return stack[-1]
[docs]
def SymWeakToRPN(exp, knownFields: List[str] = [], ipFields: List[str] = [], constants: Dict = {}, useCSE=False) -> Tuple[List[RPNElement], List[RPNElement]]:
from sympy.core.add import Add
from sympy.core.mul import Mul
from sympy.core.power import Pow
from sympy.core.numbers import Integer, Float, Half, Rational, One
from sympy.core.function import Derivative
from sympy import nsimplify
from sympy.functions.elementary.trigonometric import cos
from sympy.functions.elementary.trigonometric import sin
from sympy.functions.elementary.trigonometric import tan
from sympy.functions.elementary.exponential import log
from sympy.functions.elementary.exponential import exp as exponential
from sympy.functions.elementary.complexes import Abs
BinaryOperatorMapping = {
Add: RPNElement(type=NumericalOperators.ADDITION),
Mul: RPNElement(type=NumericalOperators.PRODUCT),
Pow: RPNElement(type=NumericalOperators.POW),
}
UnaryOperatorMapping = {
cos : RPNElement(type=NumericalOperators.COS),
sin : RPNElement(type=NumericalOperators.SIN),
tan : RPNElement(type=NumericalOperators.TAN),
log : RPNElement(type=NumericalOperators.LOG),
exponential : RPNElement(type=NumericalOperators.EXP),
Abs : RPNElement(type=NumericalOperators.ABS),
}
exp = exp.expand()
try:
if exp.shape[0] == 1 and exp.shape[1] == 1:
exp = exp[0, 0]
except:
pass
exp = nsimplify(exp)
output = []
def ParseExpression(subexp: Expr, out):
operatorStack = []
if subexp.func in BinaryOperatorMapping:
op = BinaryOperatorMapping[subexp.func]
for i, arg in enumerate(subexp.args):
ParseExpression(arg, out)
if len(subexp.args) - 1 != i:
operatorStack.append(op)
elif subexp.func in UnaryOperatorMapping:
op = UnaryOperatorMapping[subexp.func]
arg = subexp.args[0]
ParseExpression(arg, out)
operatorStack.append(op)
elif isinstance(subexp, Function) and subexp.name.startswith("Coordinate"):
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.COORDINATE))
elif isinstance(subexp, Function) and subexp.name.startswith("Normal_"):
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.NORMAL))
elif isinstance(subexp, Function) and subexp.name in knownFields:
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.KNOWN_FE_FIELD))
elif isinstance(subexp, Function) and subexp.name in ipFields:
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.KNOWN_IP_FIELD))
elif isinstance(subexp, Function) and str(subexp.func)[-1] == GetTestSufixChar():
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.TEST_FIELD))
elif isinstance(subexp, Function):
out.append(RPNElement(name=str(subexp.func), type=NumericalOperators.UNKNOWN_FIELD))
elif isinstance(subexp, Symbol):
if str(subexp) in constants:
out.append(RPNElement(value=constants[str(subexp)], type=NumericalOperators.NUMERICAL))
elif str(subexp) in knownFields:
out.append(RPNElement(name=str(subexp), type=NumericalOperators.KNOWN_FE_FIELD))
elif str(subexp) in ipFields:
out.append(RPNElement(name=str(subexp), type=NumericalOperators.KNOWN_IP_FIELD))
else:
raise RuntimeError(f"ERROR:!! not possible to classify symbol {subexp}")
elif isinstance(subexp, (Integer, Float)):
out.append(RPNElement(value=subexp, type=NumericalOperators.NUMERICAL))
# derivative is stored inside the field, not as an operator
elif isinstance(subexp, One):
out.append(RPNElement(value=1.0, type=NumericalOperators.NUMERICAL))
# derivative is stored inside the field, not as an operator
elif isinstance(subexp, Derivative):
len0 = len(out)
ParseExpression(subexp.args[0], out)
assert len0 == len(out) - 1
dercpt = 0
for i in range(1, len(subexp.args)):
for derorder in range(int(subexp.args[i][1])):
if dercpt == 0:
out[-1].der1 = ["x", "y", "z"].index(str(subexp.args[i][0]))
dercpt += 1
elif dercpt == 1:
out[-1].der2 = ["x", "y", "z"].index(str(subexp.args[i][0]))
dercpt += 1
else:
raise RuntimeError("Can treat only first and second derivatives")
elif isinstance(subexp, Half):
out.append(RPNElement(value=0.5, type=NumericalOperators.NUMERICAL))
elif isinstance(subexp, Rational):
out.append(RPNElement(value=subexp.evalf(), type=NumericalOperators.NUMERICAL))
else:
raise Exception(str(subexp.func))
operatorStack.reverse()
out.extend(operatorStack)
if useCSE:
from sympy import cse
replacements, reducedExprs = cse(exp)
rep = []
for s, expr in replacements:
output = []
ParseExpression(expr, output)
rep.append((str(s), output))
output = []
ParseExpression(reducedExprs[0], output)
return rep, output
else:
ParseExpression(exp, output)
return [], output
[docs]
def SymWeakToASTWeak(swf: Expr):
pass
[docs]
def CheckIntegrity_SymWeakToASTWeakFrom(GUI: bool = False):
wfToTest = CheckIntegrityWF()
print(wfToTest[0])
from Muscat.FE.SymWeakForm import GetTestField, GetField
u = GetField("u",2,sdim=3)[0]
ut, ut2 = GetTestField("u",2,sdim=3)
A0 = Symbol("A0")
A1 = Symbol("A1")
A2 = Symbol("A2")
AL = Symbol("AL")
AL2 = Symbol("AL2")
B0 = Symbol("B0")
B1 = Symbol("B1")
B2 = Symbol("B2")
wf = B0 * A0 * u * ut + A1 * B0 * u * ut.diff("x") + ut * AL + B2 * A0 * u.diff("y") * ut + A2 * B1 * u.diff("x") * ut.diff("y") + ut.diff("z") * AL2 + ut2
constants = {"A0": 1.0, "A1": 1.1, "A2": 1.2, "B0": 2.0, "B1": 2.1, "B2": 2.2, "AL": 1.9, "AL2": 1.99}
res = SymWeakToASTWeakFormClassified(
wf,
[
str("u_0'"),
str("u_1'"),
],
["u_0"],
constants=constants,
)
for re in res:
for r in re:
print(r)
return "ok"
[docs]
def CheckIntegrityRPN(GUI: bool = False):
import numpy as np
from sympy import nsimplify
from sympy import pprint
def PrintAndTest(wf, useCSE):
#
print(
f"* using useCSE {useCSE}********************************************************************************************************************************************************************** "
)
pprint(wf)
a = SymWeakToRPN(wf, useCSE=useCSE)
print(a)
b = RPNToSymWeak(a[0], a[1])
pprint(b)
test = nsimplify(b.expand()) - nsimplify(wf.expand())
if test != 0:
print("************************ ERROR *********************** ")
print(" *********************** RPN representation *********************** ")
print(a)
print(" *********************** Original expression *********************** ")
pprint(wf)
print(" *********************** Reconstructed expression *********************** ")
pprint(b)
print(" *********************** Error *********************** ")
print((wf - b).expand())
raise Exception(f" problem treating formula")
return a
wfToTest = CheckIntegrityWF()
for wf in wfToTest:
PrintAndTest(wf, useCSE=False)
PrintAndTest(wf, useCSE=True)
return "ok"
[docs]
def CheckIntegrityWF():
from Muscat.FE.SymWeakForm import (
Symbol,
GetField,
GetTestField,
Strain,
ToVoigtEpsilon,
)
u = GetField("u",1,sdim=3)
ut = GetTestField("u",1,sdim=3)
a = Symbol("alpha")
b = Symbol("Beta")
g = Symbol("Gamma")
m = Symbol("Mu")
wfToTest = [
a,
a + b,
a + b + g,
a * b + g,
a * b + g * g,
a * b + g**0.5,
a * b / (g**0.5 + 5 + m),
a * b / (g**0.5 + 5 + m) + a * b,
u[0],
u[0].diff(Symbol("x")),
u[0].diff(Symbol("x")) * ut[0].diff(Symbol("x")),
]
from Muscat.FE.SymPhysics import MechPhysics
phys = MechPhysics(dim=3)
H00 = phys.GetHookeOperator()[
0,
0,
]
wfToTest.append(H00)
wf = phys.PostTreatmentFormulations()["squared_von_mises"][0] ** 0.5
wfToTest.append(wf)
return wfToTest
[docs]
def CheckIntegrity(GUI: bool = False):
return CheckIntegrity_SymWeakToASTWeakFrom(GUI)
return CheckIntegrityRPN()
if __name__ == "__main__":
print(CheckIntegrity(GUI=True)) # pragma: no cover