# -*- 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 GetCoordinate(size: MuscatIndex):
"""Function to get the coordinate of the current point. The Coordinate symbol will be populated with the proper value when processed by the integrator
Parameters
----------
size : MuscatIndex
size: MuscatIndex, dimension of the space. Usually taken equal to the problem dimension
Returns
-------
sympy.Matrix
sympy.Matrix of shape (size,1) with the special denomination "Coordinate"
"""
return GetField("Coordinate", size, sdim=size)
[docs]
def GetNormal(size: MuscatIndex):
"""
Function to get the normal of an element. The Normal symbol will be populated with the proper value when processed by the integrator
Input:
size: MuscatIndex, dimension of the normal. Usually taken equal to the problem dimension
Returns
-------
sympy.Matrix of shape (size,1) with the special denomination "Normal"
"""
return GetField("Normal", size, sdim=size)
[docs]
def GetConstant(name, size=1) -> Union[Symbol, Matrix]:
"""
Function to populate a constant field, vectorial or scalar
Input:
name: str, label of the constant field if the field is vectorial then the
Returns
-------
sympy.Symbol if size=1 or
"""
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):
"""
Factory to create a constant scalar field from multiple input types
Returns
-------
sympy.Symbol if alpha is string, float if alpha is int or float.
"""
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: str, size: MuscatIndex, sdim: MuscatIndex, extraCoordinates=[]) -> Matrix:
"""
Function to instanciate a test field. It is a simple wrapper of the function GetField with star=True
Input:
name: str, label of the field, if size > 1 then each component will get a suffix for the coordinate + and the test character
size: MuscatIndex or tuple of MuscatIndex, dimension of the field (1 for scalar, 2 for 2D, 3 for 3D), if a tuple is given a tensor field is created.
sdim : MuscatIndex, space dimension where the field will evolve
extraCoordinates: list[str], list of extra parameter of the field (time for instance)
Returns
-------
sympy.Matrix of shape (size,1) containing functions of ("x","y","z") if sdim=3, ("x","y") if sdim=2, ("x") if sdim=1, along with extraCoordinates if given
"""
return GetField(name, size, sdim=sdim, star=True, extraCoordinates=extraCoordinates)
[docs]
def GetField(name: str, size: MuscatIndex, sdim: MuscatIndex, star=False, extraCoordinates=[]) -> Matrix:
"""
Function to instanciate a scalar, a vector field or a tensor field:
Input:
name: str, label of the field, if size > 1 then each component will get a suffix for the coordinate
size: MuscatIndex or tuple of MuscatIndex, dimension of the field (1 for scalar, 2 for 2D, 3 for 3D), if a tuple is given a tensor field is created.
sdim : MuscatIndex, space dimension where the field will evolve
star: bool, whether the instanciated field is a test field (is used by the function GetTestField)
extraCoordinates: list[str], list of extra parameter of the field (time for instance)
Returns
-------
sympy.Matrix of shape (size,1) containing functions of ("x","y","z") if sdim=3, ("x","y") if sdim=2, ("x") if sdim=1, along with extraCoordinates if given
"""
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):
"""
Operator to create a matrix from a list of symbolic weak forms
Input:
arg: array or iterable or shape (n) with __len__ implemented containing weak forms (sympy expressions)
Returns
-------
Sympy Matrix of shape (n x n)
"""
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):
"""
Computing the scalar product of two vectors fields
Input:
a: vector field of shape (n,1)
b: vector field of shape (n,1)
Returns
-------
sympy.exp of the scalar product following first dimension (opposite convention than numpy)
TODO: Make implementation compliant with numpy in the case of matrices
"""
if a.shape[1]!=1 or b.shape[1]!=1:
raise NotImplementedError("Inner for matrices is not implemented yet")
return Matrix(np.inner(a.T,b.T))
[docs]
def Trace(arg):
"""
Operator to compute the trace of a sympy matrix
Input:
arg: sympy matrix
Returns
-------
sympy.expr (weak form) composed of the sum of diagonal elements
"""
if type(arg).__module__ == np.__name__:
return np.trace(arg)
else:
return sympy.trace(arg)
[docs]
def Divergence(arg, sdim=3):
"""
Operator to compute the divergence of a vector or on a matrix
Input:
arg: SympyMatrix of size (n,1) or (n,n) (n=2 or 3)
sdim: dimension of the divergence (either 2 or 3)
Returns
-------
sympy.expr if arg is a vector (shape (n,1)) or sympy.Matrix of shape (n,1) if arg is a Matrix of shape (n,n)
"""
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):
"""
Computes the gradient of a vector or scalar
Input:
arg: sympy.matrix of size (1,1) (scalar field) or (n,1) (vector field)
sdim: dimension space on which to differentiate (2 or 3 for 2D or 3D)
Returns
-------
sympy.Matrix of size (sdim,1) if arg is a scalar field else (sdim,n) if arg is a vector field. The first dim is the differentiation coordinate and the second the vector composant
"""
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):
"""
Computes the Laplacian of the scalar or vector field
Input:
arg: sympy.matrix of size (1,1) (scalar field) or (n,1) (vector field)
sdim: dimension space on which to differentiate (2 or 3 for 2D or 3D)
Returns
-------
sympy.expr if arg is a scalar field else (n,1) if arg is a vector field.
"""
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):
"""
Cross product of 3D vector field
Input:
a: sympy.matrix of size (n,1) (vector field)
b: sympy.matrix of size (n,1) (vector field)
Returns
-------
sympy.matrix of size (n,1) (vector field)
"""
shape = a.shape[0]
res = [[0]]*shape
if shape == 3:
res[0][0] = a[1,0]*b[2,0]-a[2,0]*b[1,0]
res[1][0] = -(a[0,0]*b[2,0]-a[2,0]*b[0,0])
res[2][0] = a[0,0]*b[2,0]-a[1,0]*b[0,0]
else:
raise
if type(a).__module__ == np.__name__:
return np.array(res)
else:
return Matrix(res)
[docs]
def Strain(arg, sdim=3):
"""
Computes the strain matrix of a displacement vector
Input:
arg: sympy.matrix displacement field of shape (n,1)
sdim: int dimension for the derivation (2 or 3)
Returns
-------
sympy.Matrix of shape (n,n)
"""
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,3)
u0 = GetField("u0",3,3)
pprint(Cross(u,u0))
pprint(Inner(u,u0))
ut = GetTestField("u",3,3)
f = GetField("f", 3,3)
T = GetField("T",3,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,3)
pprint(Laplacian(T))
T = GetField("T",1,3)
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