-
Notifications
You must be signed in to change notification settings - Fork 0
/
util.py
106 lines (90 loc) · 2.77 KB
/
util.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
"""
util.py
Some shared utility functions.
"""
import numpy as np
#### SOME THRESHOLD DISTRIBUTION GENERATORS ####
def retconst(const):
"""Generator for the constant const.
"""
def f(s=None):
if s is None:
return const
else:
return np.ones(s)*const
return f
def uni_dist(lo,hi):
"""Generator for uniform distribution between lo and hi.
"""
def f(s=None):
if s is None:
return lo + (hi-lo)*np.random.rand()
elif len(np.atleast_1d(s))==1:
return lo + (hi-lo)*np.random.rand(s)
else:
return lo + (hi-lo)*np.random.rand(*s)
return f
def weib_dist(l,a):
"""Weibull distribution with factor l and shape a
"""
def f(s=None):
if s is None:
return np.random.weibull(a)*l
elif len(np.atleast_1d(s))==1:
return np.random.weibull(a)*l
else:
return np.random.weibull(a)*l
return f
def norm_dist(mu,sig):
"""Normal distribution, positive only
"""
def f(s=None):
if s is None:
return np.abs((np.random.randn()+mu)*sig)
elif len(np.atleast_1d(s))==1:
return np.abs((np.random.randn(s)+mu)*sig)
else:
return np.abs((np.random.randn(*s)+mu)*sig)
return f
def strict_dist(lo=0,hi=1):
"""Generator for strictly increasing thresholds lo to hi, excluding lo.
"""
def f(s=None):
assert s is not None, 'No way to set one strict threshold'
if len(np.atleast_1d(s))==1:
return np.linspace(lo,hi,s+1)[1:]
else:
return np.repeat(np.linspace(lo,hi,s[1]+1)[1:][None,:],s[0],0)
return f
#### SOME STATE VARIABLE FUNCTIONS ####
def strain_thresh(x, ssaModel):
"""
Compute strain from grounding line.
"""
# Interpolate to locations x
U = ssaModel.U(x)
# Analytical form for strain in 1D
strains = np.log(U/ssaModel.U0)
return strains
def strain_ddt(x, ssaModel):
strainFunc = project(grad(ssaModel.U)[0,0], ssaModel.Q_cg)
return strainFunc(x)
def compute_stress(ssaModel):
def epsilon(u):
return 0.5*(nabla_grad(u)+nabla_grad(u).T)
def sigma(u):
return (2*nabla_grad(u))
def epsII(u):
eps = epsilon(u)
epsII_UFL = sqrt(eps**2 + Constant(1e-16)**2)
return epsII_UFL
def eta(u):
return Constant(ssaModel.B)*epsII(u)**(1./3-1.0)
tau11 = project(eta(ssaModel.U)*grad(ssaModel.U)[0,0], ssaModel.Q_cg)
return tau11
def strain_ddt_criticalstress(stress_c, stressFunc=compute_stress):
def strain_ddt(x, ssaModel):
stress = stressFunc(ssaModel)(x)
strainFunc = project(grad(ssaModel.U)[0,0], ssaModel.Q_cg)
return strainFunc(x)*(stress>stress_c)
return strain_ddt