test/python/drucker_prager.py


##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################

from __future__ import division, print_function

__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"

# run by scons build/posix/modellib/test/python/drucker_prager.passed from check out cd.

import os
from esys.modellib.mechanics import DruckerPrager
from esys.escript.modelframe import Link,Simulation
from esys.modellib.geometry import RectangularDomain, VectorConstrainerOverBox, UpdateGeometry
from esys.modellib.input import Sequencer, InterpolateOverBox
from esys.modellib.visualization import WriteVTK

try:
   WORKDIR=os.environ['MODELLIB_WORKDIR']
except KeyError:
   WORKDIR='.'


dbg=True

dom=RectangularDomain(debug=dbg)
dom.dim=3
dom.l=[0.5,1.,1.]
dom.n=[30,6,6]
dom.order=1
dom.integrationOrder=-1


sq=Sequencer(debug=dbg)
sq.t=0
sq.t_end=0.8
sq.dt_max=100.

iob=InterpolateOverBox(debug=dbg)
iob.domain=Link(dom,"domain")
iob.value_left_bottom_front=[-1.,0.,0.]
iob.value_right_bottom_front=[0.,0.,0.]
iob.value_left_bottom_back=[-1.,0.,0.]
iob.value_right_bottom_back=[0.,0.,0.]
iob.value_left_top_front=[-1.,0.,0.]
iob.value_right_top_front=[0.,0.,0.]
iob.value_left_top_back=[-1.,0.,0.]
iob.value_right_top_back=[0.,0.,0.]

# iob.value_left_bottom_front=[-1.,0.]
# iob.value_right_bottom_front=[0.,0.]
# iob.value_left_bottom_back=[-1.,0.]
# iob.value_right_bottom_back=[0.,0.]
# iob.value_left_top_front=[-1.,0.]
# iob.value_right_top_front=[0.,0.]
# iob.value_left_top_back=[-1.,0.]
# iob.value_right_top_back=[0.,0.]

m=DruckerPrager(debug=dbg)
m.domain=Link(dom,"domain")

cv=VectorConstrainerOverBox(debug=dbg)
cv.domain=Link(dom,"domain")
cv.value=Link(iob,"out")

cv.left=[True, False, False] 
cv.right= [True, True, True]
cv.bottom= [False, False, False]
cv.top= [False, False, False]
cv.front= [False, False, False]
cv.back= [False, False, False]

m.velocity=Link(iob,"out")
m.prescribed_velocity=Link(cv,"value_of_constraint")
m.location_prescribed_velocity=Link(cv,"location_of_constraint")
m.rel_tol=0.01

m.expansion_coefficient= 0.
m.bulk_modulus=1000.
m.shear_modulus=1.
m.plastic_stress=0.
m.friction_parameter=0.
m.dilatancy_parameter=0.
m.shear_length=m.shear_modulus*0.9


ug=UpdateGeometry(debug=dbg)
ug.domain=Link(dom,"domain")
ug.displacement=Link(m,"displacement")

vis=WriteVTK(debug=dbg)
vis.t=Link(sq)
vis.data0=Link(m,"plastic_stress")
vis.data1=Link(m,"velocity")
vis.data2=Link(m,"stress")
vis.dt=0.01
vis.filename=os.path.join(WORKDIR,"temp.vtu")

s=Simulation([sq,cv,m,vis],debug=True)
s.writeXML()
s.run()
1	ksteube	1809
2	jfenwick	3981	##############################################################################
3	ksteube	1312	#
4	jfenwick	6651	# Copyright (c) 2003-2018 by The University of Queensland
5	jfenwick	3981	# http://www.uq.edu.au
6	ksteube	1312	#
7	ksteube	1809	# Primary Business: Queensland, Australia
8	jfenwick	6112	# Licensed under the Apache License, version 2.0
9			# http://www.apache.org/licenses/LICENSE-2.0
10	ksteube	1312	#
11	jfenwick	3981	# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12	jfenwick	4657	# Development 2012-2013 by School of Earth Sciences
13			# Development from 2014 by Centre for Geoscience Computing (GeoComp)
14	jfenwick	3981	#
15			##############################################################################
16	gross	821
17	sshaw	5706	from __future__ import division, print_function
18
19	jfenwick	6651	__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20	jfenwick	3981	http://www.uq.edu.au
21	ksteube	1809	Primary Business: Queensland, Australia"""
22	jfenwick	6112	__license__="""Licensed under the Apache License, version 2.0
23			http://www.apache.org/licenses/LICENSE-2.0"""
24	jfenwick	2344	__url__="https://launchpad.net/escript-finley"
25	gross	821
26			# run by scons build/posix/modellib/test/python/drucker_prager.passed from check out cd.
27
28	gross	829	import os
29	gross	821	from esys.modellib.mechanics import DruckerPrager
30			from esys.escript.modelframe import Link,Simulation
31	gross	911	from esys.modellib.geometry import RectangularDomain, VectorConstrainerOverBox, UpdateGeometry
32	gross	821	from esys.modellib.input import Sequencer, InterpolateOverBox
33	gross	829	from esys.modellib.visualization import WriteVTK
34	gross	821
35	gross	829	try:
36			WORKDIR=os.environ['MODELLIB_WORKDIR']
37			except KeyError:
38			WORKDIR='.'
39
40
41	gross	918	dbg=True
42	gross	821
43	gross	918	dom=RectangularDomain(debug=dbg)
44	gross	872	dom.dim=3
45			dom.l=[0.5,1.,1.]
46			dom.n=[30,6,6]
47			dom.order=1
48			dom.integrationOrder=-1
49	gross	821
50
51	gross	918	sq=Sequencer(debug=dbg)
52	gross	821	sq.t=0
53	gross	872	sq.t_end=0.8
54			sq.dt_max=100.
55	gross	821
56	gross	918	iob=InterpolateOverBox(debug=dbg)
57	gross	821	iob.domain=Link(dom,"domain")
58	gross	872	iob.value_left_bottom_front=[-1.,0.,0.]
59	gross	821	iob.value_right_bottom_front=[0.,0.,0.]
60	gross	872	iob.value_left_bottom_back=[-1.,0.,0.]
61	gross	821	iob.value_right_bottom_back=[0.,0.,0.]
62	gross	872	iob.value_left_top_front=[-1.,0.,0.]
63	gross	821	iob.value_right_top_front=[0.,0.,0.]
64	gross	872	iob.value_left_top_back=[-1.,0.,0.]
65	gross	821	iob.value_right_top_back=[0.,0.,0.]
66
67	gross	872	# iob.value_left_bottom_front=[-1.,0.]
68			# iob.value_right_bottom_front=[0.,0.]
69			# iob.value_left_bottom_back=[-1.,0.]
70			# iob.value_right_bottom_back=[0.,0.]
71			# iob.value_left_top_front=[-1.,0.]
72			# iob.value_right_top_front=[0.,0.]
73			# iob.value_left_top_back=[-1.,0.]
74			# iob.value_right_top_back=[0.,0.]
75	gross	829
76	gross	918	m=DruckerPrager(debug=dbg)
77	gross	821	m.domain=Link(dom,"domain")
78
79	gross	918	cv=VectorConstrainerOverBox(debug=dbg)
80	gross	821	cv.domain=Link(dom,"domain")
81			cv.value=Link(iob,"out")
82	gross	825
83	gross	821	cv.left=[True, False, False]
84	gross	872	cv.right= [True, True, True]
85			cv.bottom= [False, False, False]
86	gross	821	cv.top= [False, False, False]
87	gross	872	cv.front= [False, False, False]
88	gross	821	cv.back= [False, False, False]
89
90	gross	829	m.velocity=Link(iob,"out")
91	gross	821	m.prescribed_velocity=Link(cv,"value_of_constraint")
92			m.location_prescribed_velocity=Link(cv,"location_of_constraint")
93	gross	841	m.rel_tol=0.01
94	gross	821
95	gross	825	m.expansion_coefficient= 0.
96	gross	834	m.bulk_modulus=1000.
97	gross	825	m.shear_modulus=1.
98			m.plastic_stress=0.
99			m.friction_parameter=0.
100			m.dilatancy_parameter=0.
101	gross	872	m.shear_length=m.shear_modulus*0.9
102	gross	821
103	gross	825
104	gross	918	ug=UpdateGeometry(debug=dbg)
105	gross	823	ug.domain=Link(dom,"domain")
106			ug.displacement=Link(m,"displacement")
107
108	gross	918	vis=WriteVTK(debug=dbg)
109	gross	829	vis.t=Link(sq)
110	gross	912	vis.data0=Link(m,"plastic_stress")
111			vis.data1=Link(m,"velocity")
112			vis.data2=Link(m,"stress")
113	gross	872	vis.dt=0.01
114	caltinay	3344	vis.filename=os.path.join(WORKDIR,"temp.vtu")
115	gross	829
116	gross	941	s=Simulation([sq,cv,m,vis],debug=True)
117	gross	917	s.writeXML()
118	gross	821	s.run()