-
Notifications
You must be signed in to change notification settings - Fork 0
/
splinemodel.py
161 lines (128 loc) · 6.53 KB
/
splinemodel.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""Mathematical model that describes the partial derivatives of the potential ϕ
and the auxiliary quantities needed to compute them.
ϕ = ϕ(u,v,w) (3-parameter model) or ϕ = ϕ(u,v) (2-parameter model).
The derivatives ∂ϕ/∂u, ∂ϕ/∂v, ∂ϕ/∂w must be externally supplied.
In terms of these ∂ϕ/∂u, ∂ϕ/∂v, ∂ϕ/∂w, any partial derivatives introduced
by the chain rule are automatically generated.
This module does **not** implement splines; the spline implementation is
to be externally supplied by the user-defined additional stage1 interfaces.
This only handles the layer cake of auxiliary variables, to translate
from the physical fields (B, ε) to the auxiliary domain (u, v, w)
where the splines are defined.
Created on Tue Oct 24 14:07:45 2017
@author: Juha Jeronen <juha.jeronen@tut.fi>
"""
import sympy as sy
import symutil
from potentialmodelbase import PotentialModelBase
class Model(PotentialModelBase):
"""Generate mathematical expressions for the spline-based (B,ε) model."""
def __init__(self, kind="2par"): # 2par, 3par
"""Constructor.
Parameters:
kind: str
"2par" to build the 2-parameter model where ϕ = ϕ(u,v),
using only the invariants I4 and I5.
"3par" to build the 3-parameter model where ϕ = ϕ(u,v,w),
using all three invariants I4, I5 and I6.
"""
if kind not in ("2par", "3par"):
raise ValueError("Unknown kind '{invalid}'; valid: '2par', '3par'".format(invalid=kind))
self.label = kind
super().__init__()
makef = symutil.make_function
# Set up applied-function definitions for the layer cake.
I4 = makef("I4", *self.Bs) # i.e. in math notation, I4 = I4(Bx, By, Bz)
I5 = makef("I5", *(self.Bs + self.es)) # deviatoric strain!
I6 = makef("I6", *(self.Bs + self.es))
self.Is = I4, I5, I6
# u', v', w' are the raw u, v, w before normalization.
up = makef("up", I4) # applied function is a symbol; ok as dependency.
vp = makef("vp", I4, I5)
wp = makef("wp", I4, I5, I6)
self.ups = up, vp, wp
# The final u, v, w are normalized: u ∈ [0,1], v ∈ [-1,1], w ∈ [-1,1]
u = makef("u", up) # ...but here we only declare a formal dependency.
v = makef("v", vp)
w = makef("w", wp)
self.us = u, v, w
# Finally, the normalized u, v, w are the formal parameters of ϕ.
self.ϕ = makef("ϕ", u, v) if kind == "2par" else makef("ϕ", u, v, w)
def define_api(self):
"""See docstring for ``ModelBase.define_api()``."""
defs = {}
keyify = self.keyify
# Define the derivatives of ϕ(u, v, w) in terms of B and ε, while
# leaving ϕ itself unspecified (except its dependencies).
# For completeness, provide a function to evaluate ϕ(B, ε). We would
# like to say ϕ'(B, ε) = ϕ(u, v, w), and then drop the prime.
#
# But the LHS (function we export) cannot be named "ϕ", because in the
# spline model, the user-defined (additional stage1) interfaces are
# expected to provide phi(u, v, w), with which "ϕ" would conflict
# after degreeking.
#
# We want to keep that one as "phi", so that on the RHS, ϕ = ϕ(u, v, w),
# consistently with how __init__() defines self.ϕ. Hence we name our
# export as "ϕp", which stands for "phi prime".
#
# On the RHS, we put just "ϕ", thus telling stage2 that ϕ' depends
# on the user-defined function ϕ(u, v, w). Then stage2 does the rest,
# so that the public API for ϕ' indeed takes B and ε as its args.
print("model: {label} forming expression for ϕ".format(label=self.label))
sym, expr = self.dϕdq(qs=(), strip=False)
defs[keyify(sy.symbols("ϕp"))] = expr
# All 1st and 2nd derivatives of ϕ w.r.t. the independent vars B and ε.
# Formally, without inserting expressions.
defs.update(self.dϕdqs())
# Define the quantities appearing at the various layers of the ϕ cake.
print("model: {label} writing definitions".format(label=self.label))
B = sy.Matrix(self.Bs) # Magnetic flux density (as column vector)
ε = symutil.voigt_to_mat(self.εs) # Cauchy strain
e = symutil.voigt_to_mat(self.es) # Deviatoric strain
εM_expr = sy.factor(sy.S("1/3") * ε.trace()) # mean volumetric strain
defs[keyify(sy.symbols("εM"))] = εM_expr # will be inserted to e_expr; just a convenience
# e in terms of ε
val = ε - εM_expr * sy.eye(3)
assert symutil.is_symmetric(val)
for _, (r, c) in symutil.voigt_mat_idx():
defs[keyify(e[r, c])] = val[r, c]
# I4, I5, I6 in terms of (B, e)
I4, I5, I6 = self.Is
for key, val, label in ((I4, B.T * B, None),
(I5, B.T * e * B, None),
(I6, B.T * e * e * B, "3par")): # only in 3par model
if label is None or label == self.label:
assert val.shape == (1,1) # result should be scalar
expr = val[0,0] # extract scalar from matrix wrapper
defs[keyify(key)] = self.simplify(expr)
# u', v', w' in terms of (I4, I5, I6)
up, vp, wp = self.ups
for key, val, label in ((up, sy.sqrt(I4), None),
(vp, sy.S("3/2") * I5 / I4, None),
(wp, sy.sqrt(I6*I4 - I5**2) / I4, "3par")):
if label is None or label == self.label:
defs[keyify(key)] = val # no simplification possible; just save.
# u, v, w in terms of (u', v', w')
u, v, w = self.us
u0, v0, w0 = sy.symbols("u0, v0, w0")
for key, val, label in ((u, up / u0, None),
(v, vp / v0, None),
(w, wp / w0, "3par")):
if label is None or label == self.label:
defs[keyify(key)] = val
assert all(isinstance(key, (sy.Symbol, sy.Derivative)) for key in defs)
return defs
def test():
def scrub(expr):
return symutil.strip_function_arguments(symutil.derivatives_to_names_in(expr))
for kind in ("2par", "3par"):
print(kind)
m = Model(kind)
api = m.define_api()
api_humanreadable = {scrub(k): scrub(v) for k, v in api.items()}
print(api_humanreadable)
if __name__ == '__main__':
test()