Skip to content

Commit 2de3c8f

Browse files
vlukesvlukes
committed
init commit
0 parents  commit 2de3c8f

6 files changed

+4579
-0
lines changed

‎README.md

+9
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
Two-scale numerical simulation of a large deforming fluid-saturated porous structure
2+
====================================================================================
3+
4+
5+
* [Example Page](http://sfepy.org/sfepy_examples/example_largedef_porous/)
6+
7+
* [The Paper](https://arxiv.org/abs/2012.03730)
8+
9+
* [SfePy](https://sfepy.org)

‎largedef_porous_mac.py

+359
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,359 @@
1+
# This example implements homogenization of porous structures undergoing finite strains.
2+
#
3+
# largedef_porous_mac.py - problem at (global) macroscopic level
4+
# largedef_porous_mic.py - local subproblems, homogenized coefficients
5+
#
6+
# The mathematical model and numerical results are described in:
7+
#
8+
# LUKEŠ V., ROHAN E.
9+
# Homogenization of large deforming fluid-saturated porous structures
10+
# https://arxiv.org/abs/2012.03730
11+
#
12+
# Run simulation:
13+
#
14+
# ./simple.py example_largedef_porous-1/largedef_porous_mac.py
15+
#
16+
# The results are stored in `example_largedef_porous-1/results` directory.
17+
#
18+
19+
import numpy as nm
20+
import six
21+
import os.path as osp
22+
from sfepy import data_dir
23+
from sfepy.base.base import Struct, output, debug
24+
from sfepy.terms.terms_hyperelastic_ul import HyperElasticULFamilyData
25+
from sfepy.homogenization.micmac import get_homog_coefs_nonlinear
26+
import sfepy.linalg as la
27+
from sfepy.solvers.ts import TimeStepper
28+
from sfepy.discrete.state import State
29+
30+
wdir = osp.dirname(__file__)
31+
32+
hyperelastic_data = {
33+
'update_materials': True,
34+
'state': {'u': None, 'du': None,
35+
'p': None, 'dp': None},
36+
'mapping0': None,
37+
'coors0': None,
38+
'macro_data': None,
39+
}
40+
41+
def post_process(out, pb, state, extend=False):
42+
ts = hyperelastic_data['ts']
43+
44+
if isinstance(state, dict):
45+
pass
46+
else:
47+
stress = pb.evaluate('ev_volume_integrate_mat.i.Omega(solid.S, u)',
48+
mode='el_avg')
49+
50+
out['cauchy_stress'] = Struct(name='output_data',
51+
mode='cell',
52+
data=stress)
53+
54+
ret_stress = pb.evaluate('ev_volume_integrate_mat.i.Omega(solid.Q, u)',
55+
mode='el_avg')
56+
57+
out['retardation_stress'] = Struct(name='output_data',
58+
mode='cell',
59+
data=ret_stress)
60+
61+
strain = pb.evaluate('ev_volume_integrate_mat.i.Omega(solid.E, u)',
62+
mode='el_avg')
63+
64+
out['green_strain'] = Struct(name='output_data',
65+
mode='cell',
66+
data=strain)
67+
68+
he_state = hyperelastic_data['state']
69+
for ch in pb.conf.chs:
70+
plab = 'p%d' % ch
71+
out['p0_%d' % ch] = Struct(name='output_data',
72+
mode='vertex',
73+
data=he_state[plab][:, nm.newaxis])
74+
75+
dvel = pb.evaluate('ev_diffusion_velocity.i.Omega(solid.C%d, %s)' % (ch, plab),
76+
mode='el_avg')
77+
out['w%d' % ch] = Struct(name='output_data',
78+
mode='cell',
79+
data=dvel)
80+
81+
out['u0'] = Struct(name='output_data',
82+
mode='vertex',
83+
data=he_state['u'])
84+
85+
return out
86+
87+
88+
def homog_macro_map(ccoors, macro, nel):
89+
nqpe = ccoors.shape[0] // nel
90+
macro_ = {k: nm.sum(v.reshape((nel, nqpe) + v.shape[1:]), axis=1) / nqpe
91+
for k, v in macro.items()}
92+
macro_['recovery_idxs'] = []
93+
ccoors_ = nm.sum(ccoors.reshape((nel, nqpe) + ccoors.shape[1:]), axis=1) / nqpe
94+
95+
return ccoors_, macro_
96+
97+
98+
def homog_macro_remap(homcf, ncoor):
99+
nqpe = ncoor // homcf['Volume_total'].shape[0]
100+
homcf_ = {k: nm.repeat(v, nqpe, axis=0) for k, v in homcf.items()
101+
if not k == 'Volume_total'}
102+
103+
return homcf_
104+
105+
106+
def get_homog_mat(ts, coors, mode, term=None, problem=None, **kwargs):
107+
hyperela = hyperelastic_data
108+
ts = hyperela['ts']
109+
110+
output('get_homog_mat: mode=%s, update=%s'\
111+
% (mode, hyperela['update_materials']))
112+
113+
if not mode == 'qp':
114+
return
115+
116+
if not hyperela['update_materials']:
117+
out = hyperela['homog_mat']
118+
return {k: nm.array(v) for k, v in six.iteritems(out)}
119+
120+
dim = problem.domain.mesh.dim
121+
nqp = coors.shape[0]
122+
123+
state_u = problem.equations.variables['u']
124+
if len(state_u.field.mappings0) == 0:
125+
state_u.field.get_mapping(term.region, term.integral,
126+
term.integration)
127+
state_u.field.save_mappings()
128+
129+
state_u.field.clear_mappings()
130+
state_u.set_data(hyperela['state']['u'].ravel()) # + state_u.data[-1][state_u.indx]
131+
132+
mtx_f = problem.evaluate('ev_def_grad.i.Omega(u)',
133+
mode='qp').reshape(-1, dim, dim)
134+
135+
# relative deformation gradient
136+
if hasattr(problem, 'mtx_f_prev'):
137+
rel_mtx_f = la.dot_sequences(mtx_f, nm.linalg.inv(problem.mtx_f_prev),
138+
'AB')
139+
else:
140+
rel_mtx_f = mtx_f
141+
142+
problem.mtx_f_prev = mtx_f.copy()
143+
144+
macro_data = {
145+
'mtx_e_rel': rel_mtx_f - nm.eye(dim), # relative macro strain
146+
}
147+
148+
for ch in problem.conf.chs:
149+
plab = 'p%d' % ch
150+
state_p = problem.equations.variables[plab]
151+
state_p.set_data(hyperela['state'][plab])
152+
macro_data['p%d_0' % ch] = \
153+
problem.evaluate('ev_volume_integrate.i.Omega(p%d)' % ch,
154+
mode='qp').reshape(-1, 1, 1)
155+
macro_data['gp%d_0' % ch] = \
156+
problem.evaluate('ev_grad.i.Omega(p%d)' % ch,
157+
mode='qp').reshape(-1, dim, 1)
158+
159+
state_p.set_data(hyperela['state']['d' + plab])
160+
macro_data['dp%d_0' % ch] = \
161+
problem.evaluate('ev_volume_integrate.i.Omega(p%d)' % ch,
162+
mode='qp').reshape(-1, 1, 1)
163+
macro_data['gdp%d_0' % ch] = \
164+
problem.evaluate('ev_grad.i.Omega(p%d)' % ch,
165+
mode='qp').reshape(-1, dim, 1)
166+
167+
nel = term.region.entities[-1].shape[0]
168+
ccoors0, macro_data0 = homog_macro_map(coors, macro_data, nel)
169+
macro_data0['macro_ccoor'] = ccoors0
170+
out0 = get_homog_coefs_nonlinear(ts, ccoors0, mode, macro_data0,
171+
term=term, problem=problem,
172+
iteration=ts.step, **kwargs)
173+
out0['C1'] += nm.eye(2) * 1e-12 # ! auxiliary permeability
174+
out0['C2'] += nm.eye(2) * 1e-12 # ! auxiliary permeability
175+
out = homog_macro_remap(out0, nqp)
176+
177+
# Green strain
178+
out['E'] = 0.5 * (la.dot_sequences(mtx_f, mtx_f, 'ATB') - nm.eye(dim))
179+
180+
for ch in problem.conf.chs:
181+
out['B%d' % ch] = out['B%d' % ch].reshape((nqp, dim, dim))
182+
out['Q'] = out['Q'].reshape((nqp, dim, dim))
183+
184+
hyperela['time'] = ts.step
185+
hyperela['homog_mat'] = \
186+
{k: nm.array(v) for k, v in six.iteritems(out)}
187+
hyperela['update_materials'] = False
188+
hyperela['macro_data'] = macro_data
189+
190+
return out
191+
192+
193+
def incremental_algorithm(pb):
194+
hyperela = hyperelastic_data
195+
chs = pb.conf.chs
196+
ts = pb.conf.ts
197+
198+
hyperela['ts'] = ts
199+
hyperela['ofn_trunk'] = pb.ofn_trunk + '_%03d'
200+
pb.domain.mesh.coors_act = pb.domain.mesh.coors.copy()
201+
202+
pbvars = pb.get_variables()
203+
he_state = hyperela['state']
204+
205+
out = []
206+
out_data ={}
207+
208+
coors0 = pbvars['u'].field.get_coor()
209+
210+
he_state['coors0'] = coors0.copy()
211+
he_state['u'] = nm.zeros_like(coors0)
212+
he_state['du'] = nm.zeros_like(coors0)
213+
214+
for ch in chs:
215+
plab = 'p%d' % ch
216+
press0 = pbvars[plab].field.get_coor()[:, 0].squeeze()
217+
he_state[plab] = nm.zeros_like(press0)
218+
he_state['d' + plab] = nm.zeros_like(press0)
219+
220+
for step, time in ts:
221+
print('>>> step %d (%e):' % (step, time))
222+
hyperela['update_materials'] = True
223+
224+
pb.ofn_trunk = hyperela['ofn_trunk'] % step
225+
226+
yield pb, out
227+
228+
state = out[-1][1]
229+
result = state.get_parts()
230+
du = result['u']
231+
232+
he_state['u'] += du.reshape(he_state['du'].shape)
233+
he_state['du'][:] = du.reshape(he_state['du'].shape)
234+
pb.set_mesh_coors(he_state['u'] + he_state['coors0'],
235+
update_fields=True, actual=True, clear_all=False)
236+
237+
for ch in chs:
238+
plab = 'p%d' % ch
239+
dp = result[plab]
240+
he_state[plab] += dp
241+
he_state['d' + plab][:] = dp
242+
243+
out_data = post_process(out_data, pb, state, extend=False)
244+
filename = pb.get_output_name()
245+
pb.save_state(filename, out=out_data)
246+
247+
yield None
248+
249+
print('<<< step %d finished' % step)
250+
251+
252+
def move(ts, coor, problem=None, ramp=0.4, **kwargs):
253+
ts = problem.conf.ts
254+
nrs = round(ts.n_step * ramp)
255+
switch = 1 if (ts.step <= nrs) and (ts.step > 0) else 0
256+
displ = nm.ones((coor.shape[0],)) * problem.conf.move_val / nrs * switch
257+
258+
return displ
259+
260+
261+
def define():
262+
chs = [1, 2]
263+
264+
ts = TimeStepper(0, 0.15, n_step=30)
265+
266+
options = {
267+
'output_dir': osp.join(wdir, 'results'),
268+
'micro_filename': osp.join(wdir, 'largedef_porous_mic.py'),
269+
'parametric_hook': 'incremental_algorithm',
270+
}
271+
272+
materials = {
273+
'solid': 'get_homog',
274+
}
275+
276+
fields = {
277+
'displacement': ('real', 'vector', 'Omega', 1),
278+
'pressure': ('real', 'scalar', 'Omega', 1),
279+
}
280+
281+
variables = {
282+
'u': ('unknown field', 'displacement', 0),
283+
'v': ('test field', 'displacement', 'u'),
284+
'p1': ('unknown field', 'pressure', 1),
285+
'q1': ('test field', 'pressure', 'p1'),
286+
'p2': ('unknown field', 'pressure', 2),
287+
'q2': ('test field', 'pressure', 'p2'),
288+
}
289+
290+
filename_mesh = osp.join(wdir, 'macro_mesh_3x2.vtk')
291+
mesh_d, move_val = 0.24, -0.04
292+
293+
regions = {
294+
'Omega': 'all',
295+
'Left': ('vertices in (x < 0.0001)', 'facet'),
296+
'Right': ('vertices in (x > %e)' % (mesh_d * 0.999), 'facet'),
297+
'Recovery': ('cell 1', 'cell'),
298+
}
299+
300+
ebcs = {
301+
'left_fix_all': ('Left', {'u.all': 0.0}),
302+
'right_fix_x': ('Right', {'u.0': 0.0}),
303+
'right_move_x': ('Right', {'u.1': 'move'}),
304+
}
305+
306+
micro_args = {
307+
'eps0': mesh_d / 3,
308+
'dt': ts.dt,
309+
}
310+
311+
functions = {
312+
'move': (move,),
313+
'get_homog': (lambda ts, coors, mode, **kwargs:
314+
get_homog_mat(ts, coors, mode, define_args=micro_args, **kwargs),),
315+
}
316+
317+
integrals = {
318+
'i': 3,
319+
}
320+
321+
equations = {
322+
# eq. (60)
323+
'balance_of_forces': """
324+
dw_nonsym_elastic.i.Omega(solid.A, v, u)
325+
- dw_biot.i.Omega(solid.B1, v, p1)
326+
- dw_biot.i.Omega(solid.B2, v, p2)
327+
=
328+
- dw_lin_prestress.i.Omega(solid.S, v)
329+
- dw_lin_prestress.i.Omega(solid.Q, v)""",
330+
# eq. (61), alpha = 1
331+
'mass_conservation_1': """
332+
- %e * dw_biot.i.Omega(solid.B1, u, q1)
333+
- dw_volume_dot.i.Omega(solid.G11, q1, p1)
334+
- dw_volume_dot.i.Omega(solid.G12, q1, p2)
335+
- dw_diffusion.i.Omega(solid.C1, q1, p1)
336+
=
337+
dw_volume_lvf.i.Omega(solid.Z1, q1)
338+
""" % (1 / ts.dt),
339+
# eq. (61), alpha = 2
340+
'mass_conservation_2': """
341+
- %e * dw_biot.i.Omega(solid.B2, u, q2)
342+
- dw_volume_dot.i.Omega(solid.G21, q2, p1)
343+
- dw_volume_dot.i.Omega(solid.G22, q2, p2)
344+
- dw_diffusion.i.Omega(solid.C2, q2, p2)
345+
=
346+
dw_volume_lvf.i.Omega(solid.Z2, q2)
347+
""" % (1. / ts.dt),
348+
}
349+
350+
solvers = {
351+
'ls': ('ls.scipy_direct', {}),
352+
'newton': ('nls.newton', {
353+
'eps_a': 1e-3,
354+
'eps_r': 1e-3,
355+
'i_max': 1,
356+
}),
357+
}
358+
359+
return locals()

0 commit comments

Comments
 (0)