Skip to content

Commit 45783fc

Browse files
vlukesvlukes
committed
add code
1 parent ec00f27 commit 45783fc

File tree

4 files changed

+375
-0
lines changed

4 files changed

+375
-0
lines changed

‎beam_shell.py

+80
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
"""
2+
Display geometry and results:
3+
4+
sfepy-view beam_shell.vtk -f mat_id:e:p0 --color-map=cool --camera-position="-0.18,0.07,0.57,0.04,-0.03,0.24,0.12,0.98,-0.19"
5+
sfepy-view results/beam_shell.vtk -f u_disp:wu_disp:f10:p0 0:vw:p0 --camera-position="-0.4,0.16,0.57,0.02,0,0.23,0.22,0.97,-0.11"
6+
"""
7+
import numpy as nm
8+
import os.path as osp
9+
from sfepy.base.base import Struct
10+
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
11+
from sfepy.discrete import Integral
12+
import sfepy.mechanics.shell10x as sh
13+
14+
wdir = osp.dirname(__file__)
15+
16+
filename_mesh, np = 'beam_shell.vtk', 13
17+
thickness = 0.002
18+
19+
20+
def pp_hook(out, pb, state, extend=False):
21+
return out
22+
23+
24+
regions = {
25+
'Omega': 'cells of group 1',
26+
'Left': ('vertices in (z < 0.0001)', 'vertex'),
27+
'Right': ('vertices in (z > 0.3999)', 'vertex'),
28+
'Top': ('vertices in (y > 0.01399)', 'vertex'),
29+
'Edge': ('r.Right *v r.Top', 'vertex'),
30+
}
31+
32+
force = 1e3
33+
pload = nm.array([[0.0, -force / np, 0.0, 0.0, 0.0, 0.0]] * np)
34+
35+
options = {
36+
'output_dir': osp.join(wdir, 'results'),
37+
'post_process_hook': pp_hook,
38+
}
39+
40+
materials = {
41+
'solid': ({
42+
'D': sh.create_elastic_tensor(young=210e9, poisson=0.3),
43+
'.drill': 1e-7,
44+
},),
45+
'load': ({'.val': pload},)
46+
}
47+
48+
fields = {
49+
'fu': ('real', 6, 'Omega', 1, 'H1', 'shell10x'),
50+
}
51+
52+
variables = {
53+
'u': ('unknown field', 'fu', 0),
54+
'v': ('test field', 'fu', 'u'),
55+
}
56+
57+
# Custom integral.
58+
aux = Integral('i', order=3)
59+
qp_coors, qp_weights = aux.get_qp('3_8')
60+
qp_coors[:, 2] = thickness * (qp_coors[:, 2] - 0.5)
61+
qp_weights *= thickness
62+
63+
integrals = {
64+
'i': ('custom', qp_coors, qp_weights),
65+
}
66+
67+
ebcs = {
68+
'Fixed': ('Left', {'u.all': 0.0}),
69+
}
70+
71+
equations = {
72+
'balance_of_forces: shell':
73+
"""dw_shell10x.i.Omega(solid.D, solid.drill, v, u)
74+
= dw_point_load.i.Edge(load.val, v)""",
75+
}
76+
77+
solvers = {
78+
'ls': ('ls.mumps', {}),
79+
'newton': ('nls.newton', {'eps_a': 1e-6}),
80+
}

‎beam_shell_foam.py

+147
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
"""
2+
Display geometry and results:
3+
4+
sfepy-view beam_shell_foam.vtk -f mat_id:e:m2:o0:p0 mat_id:e:m1:o1:p0 --color-map=cool --camera-position="-0.18,0.07,0.57,0.04,-0.03,0.24,0.12,0.98,-0.19"
5+
sfepy-view beam_shell_foam.vtk -f mat_id:e:m2:o1:p0 mat_id:e:m1:o0:p0 --color-map=cool --camera-position="-0.18,0.07,0.57,0.04,-0.03,0.24,0.12,0.98,-0.19"
6+
sfepy-view results/beam_shell_foam.vtk -f uf:wuf:f10:m2:p0 0:vw:m2:p0 us_disp:wus_disp:f10:m1:p1 --camera-position="-0.4,0.16,0.57,0.02,0,0.23,0.22,0.97,-0.11" --grid-vector1="0, 1.6, 0"
7+
"""
8+
import numpy as nm
9+
import os.path as osp
10+
from sfepy.base.base import Struct
11+
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
12+
from sfepy.discrete import Integral
13+
import sfepy.mechanics.shell10x as sh
14+
15+
wdir = osp.dirname(__file__)
16+
17+
filename_mesh, np = 'beam_shell_foam.vtk', 13
18+
thickness = 0.002
19+
20+
21+
def match_nodes(coors1, coors2):
22+
"""
23+
Match coordinates `coors1` with `coors2`.
24+
"""
25+
from sfepy.discrete.fem.mesh import find_map
26+
27+
if coors1.shape != coors2.shape:
28+
raise ValueError('incompatible shapes: %s == %s'\
29+
% (coors1.shape, coors2.shape))
30+
31+
i1, i2 = find_map(coors1, coors2, join=False)
32+
33+
if i1.shape[0] != coors1.shape[0]:
34+
print(coors1[i1])
35+
print(coors2[i2])
36+
ii = nm.setdiff1d(nm.arange(coors1.shape[0]), i1)
37+
print(coors1[ii])
38+
print(coors2[ii])
39+
raise ValueError('cannot match nodes!')
40+
41+
return i1, i2
42+
43+
44+
def pp_hook(out, pb, state, extend=False):
45+
stress = pb.evaluate('ev_cauchy_stress.if.OmegaF(foam.D, uf)',
46+
mode='el_eval')
47+
strain = pb.evaluate('ev_cauchy_strain.if.OmegaF(uf)', mode='el_eval')
48+
out['E'] = Struct(name='strain', mode='cell', data=strain, region='OmegaF')
49+
out['S'] = Struct(name='stress', mode='cell', data=stress, region='OmegaF')
50+
51+
idxs = pb.domain.regions['OmegaS'].vertices
52+
53+
v = out['us_disp'].data
54+
aux = nm.zeros((pb.domain.shape.n_nod, v.shape[1]), dtype=v.dtype)
55+
aux[idxs] = v
56+
out['us_disp'].data = aux
57+
58+
v = out['us_rot'].data
59+
aux = nm.zeros((pb.domain.shape.n_nod, v.shape[1]), dtype=v.dtype)
60+
aux[idxs] = v
61+
out['us_rot'].data = aux
62+
63+
return out
64+
65+
66+
regions = {
67+
# shell
68+
'OmegaS': 'cells of group 1',
69+
'LeftS': ('vertices in (z < 0.0001)', 'vertex', 'OmegaS'),
70+
'RightS': ('vertices in (z > 0.3999)', 'vertex', 'OmegaS'),
71+
'TopS': ('vertices in (y > 0.01399)', 'vertex', 'OmegaS'),
72+
'EdgeS': ('r.RightS *v r.TopS', 'vertex', 'OmegaS'),
73+
# solid foam
74+
'OmegaF': 'cells of group 2',
75+
'LeftF': ('vertices in (z < 0.0001)', 'facet', 'OmegaF'),
76+
'RightF': ('vertices in (z > 0.3999)', 'facet', 'OmegaF'),
77+
'SurfaceF_': ('vertices of surface', 'facet', 'OmegaF'),
78+
'SurfaceF__': ('r.SurfaceF_ -s r.LeftF', 'facet', 'OmegaF'),
79+
'SurfaceF': ('r.SurfaceF__ -s r.RightF', 'facet', 'OmegaF'),
80+
}
81+
82+
force = 1e3
83+
pload = nm.array([[0.0, -force / np, 0.0, 0.0, 0.0, 0.0]] * np)
84+
85+
options = {
86+
'output_dir': osp.join(wdir, 'results'),
87+
'post_process_hook': pp_hook,
88+
}
89+
90+
materials = {
91+
'shell': ({
92+
'D': sh.create_elastic_tensor(young=210e9, poisson=0.3),
93+
'.drill': 1e-7,
94+
},),
95+
'load': ({'.val': pload},),
96+
'foam': ({'D': stiffness_from_youngpoisson(3, 20e9, 0.25), },),
97+
}
98+
99+
fields = {
100+
'fu': ('real', 6, 'OmegaS', 1, 'H1', 'shell10x'),
101+
'displacement': ('real', 'vector', 'OmegaF', 1),
102+
}
103+
104+
variables = {
105+
'us': ('unknown field', 'fu', 0),
106+
'vs': ('test field', 'fu', 'us'),
107+
'uf': ('unknown field', 'displacement', 1),
108+
'vf': ('test field', 'displacement', 'uf'),
109+
}
110+
111+
# Custom integral.
112+
aux = Integral('i', order=3)
113+
qp_coors, qp_weights = aux.get_qp('3_8')
114+
qp_coors[:, 2] = thickness * (qp_coors[:, 2] - 0.5)
115+
qp_weights *= thickness
116+
117+
integrals = {
118+
'is': ('custom', qp_coors, qp_weights),
119+
'if': 2,
120+
}
121+
122+
ebcs = {
123+
'FixedS': ('LeftS', {'us.all': 0.0}),
124+
'FixedF': ('LeftF', {'uf.all': 0.0}),
125+
}
126+
127+
functions = {
128+
'match_nodes': (match_nodes,),
129+
}
130+
131+
lcbcs = {
132+
'match_surface': (['OmegaS', 'SurfaceF'], {'us.[0,1,2]': 'uf.[0,1,2]'},
133+
'match_nodes', 'match_dofs'),
134+
}
135+
136+
equations = {
137+
'balance_of_forces - shell':
138+
"""dw_shell10x.is.OmegaS(shell.D, shell.drill, vs, us)
139+
= dw_point_load.is.EdgeS(load.val, vs)""",
140+
'balance_of_forces - foam':
141+
"""dw_lin_elastic.if.OmegaF(foam.D, vf, uf) = 0""",
142+
}
143+
144+
solvers = {
145+
'ls': ('ls.mumps', {}),
146+
'newton': ('nls.newton', {'eps_a': 1e-6}),
147+
}

‎beam_solid.py

+69
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
"""
2+
Display geometry and results:
3+
4+
sfepy-view beam_solid.vtk -f mat_id:e:p0 --color-map=cool --camera-position="-0.18,0.07,0.57,0.04,-0.03,0.24,0.12,0.98,-0.19"
5+
sfepy-view results/beam_solid.vtk -f u:wu:f10:p0 0:vw:p0 --camera-position="-0.4,0.16,0.57,0.02,0,0.23,0.22,0.97,-0.11"
6+
"""
7+
import numpy as nm
8+
import os.path as osp
9+
from sfepy.base.base import Struct
10+
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
11+
from sfepy.discrete import Integral
12+
13+
wdir = osp.dirname(__file__)
14+
15+
filename_mesh, np = 'beam_solid.vtk', 13
16+
17+
18+
def pp_hook(out, pb, state, extend=False):
19+
return out
20+
21+
22+
regions = {
23+
'Omega': 'all',
24+
'Left': ('vertices in (z < 0.0001)', 'vertex', 'Omega'),
25+
'Right': ('vertices in (z > 0.3999)', 'vertex', 'Omega'),
26+
'Top': ('vertices in (y > 0.01499)', 'vertex', 'Omega'),
27+
'Edge': ('r.Right *v r.Top', 'vertex', 'Omega'),
28+
}
29+
30+
force = 1e3
31+
pload = nm.array([[0.0, -force / np, 0.0]] * np)
32+
33+
options = {
34+
'output_dir': osp.join(wdir, 'results'),
35+
'post_process_hook': pp_hook,
36+
}
37+
38+
materials = {
39+
'solid': ({'D': stiffness_from_youngpoisson(3, young=210e9, poisson=0.3),},),
40+
'load': ({'.val': pload},),
41+
}
42+
43+
fields = {
44+
'displacement': ('real', 'vector', 'Omega', 1),
45+
}
46+
47+
variables = {
48+
'u': ('unknown field', 'displacement', 1),
49+
'v': ('test field', 'displacement', 'u'),
50+
}
51+
52+
integrals = {
53+
'is': 2,
54+
}
55+
56+
ebcs = {
57+
'Fixed': ('Left', {'u.all': 0.0}),
58+
}
59+
60+
equations = {
61+
'balance_of_forces: solid':
62+
"""dw_lin_elastic.is.Omega(solid.D, v, u)
63+
= dw_point_load.is.Edge(load.val, v)""",
64+
}
65+
66+
solvers = {
67+
'ls': ('ls.mumps', {}),
68+
'newton': ('nls.newton', {'eps_a': 1e-6}),
69+
}

‎beam_solid_foam.py

+79
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
"""
2+
Display geometry and results:
3+
4+
sfepy-view beam_solid_foam.vtk -f mat_id:e:p0 --color-map=cool --camera-position="-0.18,0.07,0.57,0.04,-0.03,0.24,0.12,0.98,-0.19"
5+
sfepy-view results/beam_solid_foam.vtk -f u:wu:f10:p0 0:vw:p0 --camera-position="-0.4,0.16,0.57,0.02,0,0.23,0.22,0.97,-0.11"
6+
"""
7+
import numpy as nm
8+
import os.path as osp
9+
from sfepy.base.base import Struct
10+
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
11+
from sfepy.discrete import Integral
12+
13+
wdir = osp.dirname(__file__)
14+
15+
filename_mesh, np = 'beam_solid_foam.vtk', 13
16+
17+
18+
def pp_hook(out, pb, state, extend=False):
19+
stress = pb.evaluate('ev_cauchy_stress.i.Omega(solid.D, u)', mode='el_eval')
20+
strain = pb.evaluate('ev_cauchy_strain.i.Omega(u)', mode='el_eval')
21+
out['E'] = Struct(name='strain', mode='cell', data=strain)
22+
out['S'] = Struct(name='stress', mode='cell', data=stress)
23+
24+
return out
25+
26+
27+
regions = {
28+
'Omega': 'all',
29+
'Left': ('vertices in (z < 0.0001)', 'vertex', 'Omega'),
30+
'Right': ('vertices in (z > 0.3999)', 'vertex', 'Omega'),
31+
'Top': ('vertices in (y > 0.01499)', 'vertex', 'Omega'),
32+
'Edge': ('r.Right *v r.Top', 'vertex', 'Omega'),
33+
'Omega1': 'cells of group 1',
34+
'Omega2': 'cells of group 2',
35+
}
36+
37+
force = 1e3
38+
pload = nm.array([[0.0, -force / np, 0.0]] * np)
39+
40+
options = {
41+
'output_dir': osp.join(wdir, 'results'),
42+
'post_process_hook': pp_hook,
43+
}
44+
45+
materials = {
46+
'solid': ({'D': {
47+
'Omega1': stiffness_from_youngpoisson(3, 210e9, 0.3),
48+
'Omega2': stiffness_from_youngpoisson(3, 20e9, 0.25),
49+
}},),
50+
'load': ({'.val': pload},)
51+
}
52+
53+
fields = {
54+
'displacement': ('real', 'vector', 'Omega', 1),
55+
}
56+
57+
integrals = {
58+
'i': 2,
59+
}
60+
61+
variables = {
62+
'u': ('unknown field', 'displacement', 0),
63+
'v': ('test field', 'displacement', 'u'),
64+
}
65+
66+
ebcs = {
67+
'Fixed': ('Left', {'u.all': 0.0}),
68+
}
69+
70+
equations = {
71+
'balance_of_forces':
72+
"""dw_lin_elastic.i.Omega(solid.D, v, u)
73+
= dw_point_load.i.Edge(load.val, v)""",
74+
}
75+
76+
solvers = {
77+
'ls': ('ls.mumps', {}),
78+
'newton': ('nls.newton', {'eps_a': 1e-5}),
79+
}

0 commit comments

Comments
 (0)