# -*- 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 Union
from numbers import Number
import sympy
from sympy import Symbol, Function
from sympy.matrices import Matrix
import numpy as np
from Muscat.Types import MuscatIndex
from Muscat.Helpers.Logger import Warning
testcharacter = "'"
space3D = Matrix([Symbol('x'), Symbol('y'), Symbol('z')])
space2D = Matrix([Symbol('x'), Symbol('y')])
space1D = Matrix([Symbol('x'),])
[docs]def GetTestSufixChar() -> str:
return testcharacter
[docs]def GetNormal(size: MuscatIndex):
return GetField("Normal", size)
[docs]def GetConstant(name, size=1) -> Union[Symbol, Matrix]:
if size == 1:
return Symbol(name)
else:
res = []
for i in range(size):
res.append(Symbol(name+"_"+str(i)))
return Matrix([res]).T
[docs]def GetScalarField(alpha):
if alpha is None:
a = 1.
elif isinstance(alpha, str):
a = GetConstant(alpha)
elif isinstance(alpha, (float, int)):
a = float(alpha)
else:
a = alpha
return a
[docs]def GetTestField(name, size, sdim=3, extraCoordinates=[]):
return GetField(name, size, star=True, sdim=sdim, extraCoordinates=extraCoordinates)
[docs]def GetField(name: str, size: MuscatIndex, star=False, sdim=3, extraCoordinates=[]) -> Matrix:
res = []
suffix = ""
if star:
suffix = testcharacter
s = space3D[0:sdim]
s.extend(extraCoordinates)
if isinstance(size,(MuscatIndex,int)):
if size == 1:
if len(s) == 0:
res.append(Function(name+suffix)())
else:
res.append(Function(name+suffix)(*s))
else:
for i in range(size):
res.append(Function(name+"_"+str(i)+suffix)(*s))
res = [res]
return (Matrix(res)).T
else:
res = [[None]*size[1] for x in range(size[0])]
for i in range(size[0]):
for j in range(size[1]):
t = Function(name+"_"+str(i)+"_"+str(j)+suffix)(*s)
res[i][j] = t
return Matrix(res).T
[docs]def ArrayToMatrix(array, normalize=False) -> Matrix:
"""Convert a list of object into a sympy.Matrix (1D).
str objects are converted into a symbols.
If `normalize` is `True`, the array is rescaled by its L2-norm,
if all values are numerical; a warning is issued otherwise.
Parameters
----------
array : arraylike
a list of objects
normalize : bool, optional
if true and all the values in the array are numerical
then the vector is rescaled (norm(array) ==1), by default False
Returns
-------
Matrix
a sympy matrix with str converted to symbols
"""
a = list(map(GetScalarField, array))
if normalize:
is_number = lambda x: isinstance(x, Number)
if all(map(is_number, a)):
a = np.array(a) / np.linalg.norm(a)
else:
Warning( f"Could not ensure that array {array} has unit norm.")
return Matrix(np.atleast_2d(a)).T
########################## Mathematical operators ##############################
[docs]def Diag(arg):
shape = len(arg)
res = [[0]*shape for i in range(shape)]
for i, v in enumerate(arg):
res[i][i] = v
if type(arg).__module__ == np.__name__ or type(arg) == list:
return np.array(res)
else:
return Matrix(res)
[docs]def Inner(a, b):
return a@b
[docs]def Trace(arg):
if type(arg).__module__ == np.__name__:
return np.trace(arg)
else:
return sympy.trace(arg)
[docs]def Divergence(arg, sdim=3):
if arg.shape[1]==1:
return Trace(Gradient(arg, sdim=sdim))
else:
return Matrix([Divergence(arg[:,[i]],sdim=sdim) for i in range(arg.shape[1])])
[docs]def Gradient(arg, sdim=3):
shape = arg.shape[0]
res = [[0]*shape for i in range(sdim)]
for s in range(shape):
for d in range(sdim):
res[d][s] = arg[s].diff(space3D[d])
if type(arg).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def Laplacian(arg,sdim=3):
if arg.shape[1]==1:
return Divergence(Gradient(arg,sdim))
else:
return Matrix([Laplacian(arg[:,[i]],sdim=sdim) for i in range(arg.shape[1])])
[docs]def Cross(a, b):
shape = a.shape[0]
res = [[0]*shape]
if shape == 3:
res[0][0] = a[1]*b[2]-a[2]*b[1]
res[0][1] = -(a[0]*b[2]-a[2]*b[0])
res[0][2] = a[0]*b[2]-a[1]*b[0]
else:
raise
if type(a).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def Strain(arg, sdim=3):
G = Gradient(arg, sdim)
return (G+G.T)/2
[docs]def StrainAxyCol(arg, radius):
# strain definition for axisymmetric mechanical problem
u = arg[0]
v = arg[1]
r = space2D[0]
h = space2D[1]
dudr = u.diff(r)
dudh = u.diff(h)
dvdr = v.diff(r)
dvdh = v.diff(h)
# E_r, E_z, E_theta, E_rz
res = [dudr, dvdh, u/radius, dudh+dvdr]
if type(arg).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def ToVoigtEpsilon(arg):
""" https://en.wikipedia.org/wiki/Voigt_notation
"""
if arg.shape[0] == 3:
res = [arg[0, 0], arg[1, 1], arg[2, 2], 2*arg[1, 2], 2*arg[0, 2], 2*arg[0, 1], ]
elif arg.shape[0] == 2:
res = [arg[0, 0], arg[1, 1], 2*arg[0, 1]]
elif arg.shape[0] == 1:
res = [arg[0, 0]]
else:
raise ()
if type(arg).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def ToVoigtSigma(arg):
""" https://en.wikipedia.org/wiki/Voigt_notation
"""
if arg.shape[0] == 3:
res = [arg[0, 0], arg[1, 1], arg[2, 2], arg[1, 2], arg[0, 2], arg[0, 1], ]
elif arg.shape[0] == 2:
res = [arg[0, 0], arg[1, 1], arg[0, 1]]
elif arg.shape[0] == 1:
res = [arg[0, 0]]
else:
raise ()
if type(arg).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def FromVoigtSigma(arg):
""" https://en.wikipedia.org/wiki/Voigt_notation
"""
if arg.shape[0] == 6:
res = [[arg[0], arg[5], arg[4]],
[arg[5], arg[1], arg[3]],
[arg[4], arg[3], arg[2]],]
elif arg.shape[0] == 3:
res = [[arg[0], arg[2]],
[arg[2], arg[1]]]
elif arg.shape[0] == 1:
res = [arg[0]]
else:
raise ()
if type(arg).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]def CheckIntegrity(GUI: bool = False):
from sympy import pprint
# init_session()
u = GetField("u", 3)
u0 = GetField("u0", 3)
ut = GetTestField("u", 3)
f = GetField("f", 3)
T = GetField("T",3)
alpha = Symbol("alpha")
globalconstant = GetField("g", 1, sdim=0)
print(globalconstant)
print(u)
print(u.shape)
print(u.diff(Symbol("x")))
print(ut.diff(Symbol("x")))
print("-----------------")
pprint(u, use_unicode=GUI)
pprint(Gradient(u), use_unicode=GUI)
pprint(Strain(u), use_unicode=GUI)
pprint(u[0].diff(space3D[1]), use_unicode=GUI)
from Muscat.FE.MaterialHelp import HookeIso
K = HookeIso(1, 0.3)
pprint(K, use_unicode=GUI)
ener = ToVoigtEpsilon(Strain(u+u0)).T*K*ToVoigtEpsilon(Strain(ut)) + f.T*ut*alpha
pprint(ener, use_unicode=GUI)
pprint(Gradient(u))
pprint(Trace(Gradient(u)))
pprint(Divergence(u)) # divergence of vector
pprint(Divergence(FromVoigtSigma(K*ToVoigtEpsilon(Strain(ut))))) # divergence of matrix
T = GetField("T",3)
pprint(Laplacian(T))
T = GetField("T",1)
pprint(Laplacian(T))
print(type(Gradient(u)))
print(type(u))
print(type(u).__module__)
print(type(sympy))
print(sympy.__name__)
import logging
from Muscat.Helpers.Logger import CatchLog
# Test ArrayToMatrix
m = ArrayToMatrix([0, 0, 1])
assert isinstance(m, Matrix)
assert m == Matrix([[0.], [0.], [1.]])
m = ArrayToMatrix([3, -4, 0])
assert m == Matrix([[3.], [-4.], [0.]])
m = ArrayToMatrix([3, -4, 0], normalize=True)
assert m == Matrix([[0.6], [-0.8], [0.]])
m = ArrayToMatrix([3, "foo", 0])
assert m == Matrix([[3.], [Symbol("foo")], [0.]])
with CatchLog() as e:
m = ArrayToMatrix([3, "foo", 0], normalize=True)
assert e.messages['warning'][0].startswith("Could not ensure that")
assert m == Matrix([[3.], [Symbol("foo")], [0.]])
return "OK"
if __name__ == '__main__':
print(CheckIntegrity(GUI=True)) # pragma: no cover