-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathanisotropy_mixins.py
321 lines (268 loc) · 12.3 KB
/
anisotropy_mixins.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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
"""File for anisotropic stiffness tensor mixins.
Currently, the following tensors are found here:
* Vertical transverse anisotropic stiffness tensor (for inner domain). Brute force
version.
* Simple anisotropy: Only increase stiffness in one direction.
* Vertical transverse anisotropic stiffness tensor: Non-brute force version.
Constructed generally enough to allow for some flexibility in further
developments/enhancement of it."""
import numpy as np
import porepy as pp
from utils import inner_domain_cells
from utils.stiffness_tensors import StiffnessTensorInnerVTI
class TransverselyIsotropicStiffnessTensor:
"""Note: there is a less brute-force version in the bottom of this file."""
def stiffness_tensor(self, subdomain: pp.Grid) -> pp.FourthOrderTensor:
"""Stiffness tensor [Pa].
Modified to represent a transversely isotropic stiffness tensor. It is (for now)
done in a rather brute force way. All cells corresponding to the "inner domain"
(see the function inner_domain_cells in utils) will have their stiffness tensor
values wiped. A newly created 9x9 matrix corresponding to the values of a
transversely isotropic medium is assigned to the wiped part.
Parameters:
subdomain: Subdomain where the stiffness tensor is defined.
Returns:
Cell-wise stiffness tensor in SI units.
"""
lmbda = self.solid.lame_lambda * np.ones(subdomain.num_cells)
mu = self.solid.shear_modulus * np.ones(subdomain.num_cells)
stiffness_tensor = pp.FourthOrderTensor(mu, lmbda)
width = self.params.get("inner_domain_width", 0)
inner_domain_center = self.params.get("inner_domain_center", None)
if width == 0:
return stiffness_tensor
else:
indices = inner_domain_cells(
self=self,
sd=subdomain,
width=width,
inner_domain_center=inner_domain_center,
)
anisotropic_stiffness_values = (
self.transversely_isotropic_stiffness_tensor()
)
stiffness_tensor.values[:, :, indices] = anisotropic_stiffness_values[
:, :, np.newaxis
]
return stiffness_tensor
def transversely_isotropic_stiffness_tensor(self) -> np.ndarray:
"""Matrices representing the anisotropic stiffness tensor contribution.
Stiffness tensor for transverse isotropy in z-direction is created here. The
anisotropic parameter values are found within
(self.params["anisotropy_constants"]). If no values are assigned to this
dicrionary, a zero contribution is returned. That is, the stiffness tensor is
again representing an isotropic medium.
This method is used by stiffness_tensor as a utility method.
Returns:
A 9x9 matrix with the anisotropic stiffness tensor values.
"""
lambda_mat = np.array(
[
[1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1],
]
)
lambda_parallel_mat = np.array(
[
[1, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
lambda_orthogonal_mat = np.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1],
]
)
mu_parallel_mat = np.array(
[
[2, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
mu_orthogonal_mat = np.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 1, 0],
[0, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 2],
]
)
anisotropy_constants = self.params.get("anisotropy_constants", None)
if anisotropy_constants is None:
return np.zeros((9, 9))
else:
mu_parallel = anisotropy_constants["mu_parallel"]
mu_orthogonal = anisotropy_constants["mu_orthogonal"]
volumetric_compr_lambda = self.solid.lame_lambda
lambda_parallel = anisotropy_constants["lambda_parallel"]
lambda_orthogonal = anisotropy_constants["lambda_orthogonal"]
values = (
mu_parallel * mu_parallel_mat
+ mu_orthogonal * mu_orthogonal_mat
+ lambda_orthogonal * lambda_orthogonal_mat
+ lambda_parallel * lambda_parallel_mat
+ volumetric_compr_lambda * lambda_mat
)
return values
class SimpleAnisotropy:
def stiffness_tensor(self, subdomain: pp.Grid) -> pp.FourthOrderTensor:
"""Stiffness tensor [Pa].
Parameters:
subdomain: Subdomain where the stiffness tensor is defined.
Returns:
Cell-wise stiffness tensor in SI units.
"""
lmbda = self.solid.lame_lambda * np.ones(subdomain.num_cells)
mu = self.solid.shear_modulus * np.ones(subdomain.num_cells)
stiffness_tensor = pp.FourthOrderTensor(mu, lmbda)
width = self.params.get("inner_domain_width", 0)
if width == 0:
return stiffness_tensor
else:
indices = inner_domain_cells(self=self, sd=subdomain, width=width)
anisotropic_stiffness_values = self.construct_anisotropic_contribution()
for cell_index in indices:
stiffness_tensor.values[
:, :, cell_index
] += anisotropic_stiffness_values
return stiffness_tensor
def construct_anisotropic_contribution(self) -> np.ndarray:
"""Matrices representing the anisotropic stiffness tensor contribution.
The matrices found here are a simple form of anisotropy. They increase the
stiffness in one (or more) directions parallel with the coordinate axes. Which
directions that have increased stiffness depends on the values of lmbda1, lmbda2
and lmbda3 in the params dictionary (self.params["anisotropy_constants"]). If no
values are assigned, a zero contribution is returned. That is, the stiffness
tensor is again representing an isotropic medium.
This method is used by stiffness_tensor as a utility method.
Returns:
Sum of all anisotropic contributions in the shape of a 9 by 9 matrix.
"""
lmbda1_mat = np.array(
[
[1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
lmbda2_mat = np.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
]
)
lmbda3_mat = np.array(
[
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1],
]
)
anisotropy_constants = self.params.get("anisotropy_constants", None)
if anisotropy_constants is None:
return np.zeros((9, 9))
else:
lmbda1 = anisotropy_constants["lmbda1"]
lmbda2 = anisotropy_constants["lmbda2"]
lmbda3 = anisotropy_constants["lmbda3"]
values = lmbda1_mat * lmbda1 + lmbda2_mat * lmbda2 + lmbda3_mat * lmbda3
return values
class InnerDomainVTIStiffnessTensorMixin:
"""Mixin for a stiffness tensor corresponding to an VTI inner domain, isotropic
outer domain.
Instead of using the pp.FourthOrderTensor object (found within PorePy) and then
overriding certain values, we now use another tensor object specifically tailored
for a VTI medium. That is, an isotropic medium with an inner VTI one.
"""
def stiffness_tensor(self, subdomain: pp.Grid) -> StiffnessTensorInnerVTI:
# Fetch inner domain indices such that we can distribute values of the material
# parameters in arrays according to cell numbers.
inner_domain_width = self.params.get("inner_domain_width", 0)
inner_domain_center = self.params.get("inner_domain_center", None)
inner_cell_indices = inner_domain_cells(
self=self,
sd=subdomain,
width=inner_domain_width,
inner_domain_center=inner_domain_center,
)
# Preparing basis arrays for inner and outer domains
inner = np.zeros(subdomain.num_cells)
inner[inner_cell_indices] = 1
outer = np.ones(subdomain.num_cells)
outer = outer - inner
# Standard material values: These are assigned to the outer domain
lmbda = self.solid.lame_lambda * outer
mu = self.solid.shear_modulus * outer
# Anisotropy related values: These are assigned to the inner domain
anisotropy_constants = self.params["anisotropy_constants"]
mu_parallel = anisotropy_constants["mu_parallel"] * inner
mu_orthogonal = anisotropy_constants["mu_orthogonal"] * inner
volumetric_compr_lambda = (
anisotropy_constants["volumetric_compr_lambda"] * inner
)
lambda_parallel = anisotropy_constants["lambda_parallel"] * inner
lambda_orthogonal = anisotropy_constants["lambda_orthogonal"] * inner
# Finally a call to the stiffness tensor object itself
stiffness_tensor = StiffnessTensorInnerVTI(
mu=mu,
lmbda=lmbda,
mu_parallel=mu_parallel,
mu_orthogonal=mu_orthogonal,
lambda_parallel=lambda_parallel,
lambda_orthogonal=lambda_orthogonal,
volumetric_compr_lambda=volumetric_compr_lambda,
)
return stiffness_tensor