examples/cookbook/cblib1.py

########################################################
#
# Copyright (c) 2009 by University of Queensland
# Earth Systems Science Computational Center (ESSCC)
# http://www.uq.edu.au/esscc
#
# Primary Business: Queensland, Australia
# Licensed under the Open Software License version 3.0
# http://www.opensource.org/licenses/osl-3.0.php
#
########################################################

__copyright__="""Copyright (c) 2009 by University of Queensland
Earth Systems Science Computational Center (ESSCC)
http://www.uq.edu.au/esscc
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Open Software License version 3.0
http://www.opensource.org/licenses/osl-3.0.php"""
__url__="https://launchpad.net/escript-finley"

"""
A collection of routines to use in cookbook examples.
Author: Antony Hallam antony.hallam@uqconnect.edu.au
"""

from esys.pycad import CurveLoop

# numpy for array handling
import numpy as np

# tools for dealing with PDEs - contains locator
from esys.escript.pdetools import Locator, Projector

from esys.escript import *
from esys.escript.linearPDEs import LinearPDE

# routine to find consecutive coordinates of a loop in pycad
def getLoopCoords(loop):
    # return all construction points of input
    temp = loop.getConstructionPoints()
    #create a numpy array for xyz components or construction points
    coords = np.zeros([len(temp),3],float)
    #place construction points in array
    for i in range(0,len(temp)):
        coords[i,:]=temp[i].getCoordinates()
    #return a numpy array
    return coords
    

########################################################
# subroutine: cbphones
# Allows us to record the values of a PDE at various 
# specified locations in the model.
# Arguments:
#   domain  : domain of model
#   U       : Current time state displacement solution.
#   phones  : Geophone Locations
#   dim     : model dimesions
#   savepath: where to output the data files local is default
########################################################
def cbphones(domain,U,phones,dim,savepath=""):
   #find the number of geophones
   nphones = len(phones)
   u_pot = np.zeros([nphones,dim],float)
   
   for i in range(0,nphones):
     # define the location of the phone source 
     L=Locator(domain,numpy.array(phones[i]))
     # find potential at point source.
     temp = L.getValue(U)
     for j in range(0,dim):
       u_pot[i,j]=temp[j]

   # open file to save displacement at point source
   return u_pot

########################################################
# subroutine: wavesolver2d
# Can solve a generic 2D wave propagation problem with a
# point source in a homogeneous medium.
# Arguments:
#   domain  : domain to solve over
#   h       : time step
#   tend    : end time
#   lam, mu : lames constants for domain
#   rho : density of domain
#   U0  : magnitude of source
#   xc  : source location in domain (Vector)
#   savepath: where to output the data files
########################################################
def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
   from esys.escript.linearPDEs import LinearPDE
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   #mypde.setSolverMethod(LinearPDE.LUMPING)
   mypde.setSymmetryOn()
   kmat = kronecker(domain)
   mypde.setValue(D=kmat*rho)

   # define small radius around point xc
   # Lsup(x) returns the maximum value of the argument x
   src_radius = 50#2*Lsup(domain.getSize())
   print "src_radius = ",src_radius

   dunit=numpy.array([0.,1.]) # defines direction of point source

   
   # ... set initial values ....
   n=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
   u_m1=u
   t=0

   u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
   u_pc_x1 = u_pot[0,0]
   u_pc_y1 = u_pot[0,1]
   u_pc_x2 = u_pot[1,0]
   u_pc_y2 = u_pot[1,1]
   u_pc_x3 = u_pot[2,0]
   u_pc_y3 = u_pot[2,1]

   # open file to save displacement at point source
   u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
   u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
 
   while t<tend:
     # ... get current stress ....
     
##OLD WAY
     g=grad(u)
     stress=lam*trace(g)*kmat+mu*(g+transpose(g))
     ### ... get new acceleration ....
     #mypde.setValue(X=-stress)          
     #a=mypde.getSolution()
     ### ... get new displacement ...
     #u_p1=2*u-u_m1+h*h*a
###NEW WAY
     mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))
     u_p1 = mypde.getSolution()
     # ... shift displacements ....
     u_m1=u
     u=u_p1
     #stress = 
     t+=h
     n+=1
     print n,"-th time step t ",t
     u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)

#     print "u at point charge=",u_pc
     u_pc_x1 = u_pot[0,0]
     u_pc_y1 = u_pot[0,1]
     u_pc_x2 = u_pot[1,0]
     u_pc_y2 = u_pot[1,1]
     u_pc_x3 = u_pot[2,0]
     u_pc_y3 = u_pot[2,1]
           
     # save displacements at point source to file for t > 0
     u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
 
     # ... save current acceleration in units of gravity and displacements 
     #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
     #displacement = length(u), tensor = stress, Ux = u[0] )
     saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)

   u_pc_data.close()
   

########################################################
# subroutine: wavesolver2d
# Can solve a generic 2D wave propagation problem with a
# point source in a homogeneous medium with friction.
# Arguments:
#   domain  : domain to solve over
#   h       : time step
#   tend    : end time
#   lam, mu : lames constants for domain
#   rho : density of domain
#   U0  : magnitude of source
#   xc  : source location in domain (Vector)
#   savepath: where to output the data files
########################################################
def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   #mypde.setSolverMethod(LinearPDE.LUMPING)
   mypde.setSymmetryOn()
   kmat = kronecker(domain)
   mypde.setValue(D=kmat)
   b=0.9

   # define small radius around point xc
   # Lsup(x) returns the maximum value of the argument x
   src_radius = 50#2*Lsup(domain.getSize())
   print "src_radius = ",src_radius

   dunit=numpy.array([0.,1.]) # defines direction of point source

   
   # ... set initial values ....
   n=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit
   u_m1=u
   t=0

   u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
   u_pc_x1 = u_pot[0,0]
   u_pc_y1 = u_pot[0,1]
   u_pc_x2 = u_pot[1,0]
   u_pc_y2 = u_pot[1,1]
   u_pc_x3 = u_pot[2,0]
   u_pc_y3 = u_pot[2,1]

   # open file to save displacement at point source
   u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
   u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
 
   while t<tend:
     # ... get current stress ....
     
##OLD WAY
     g=grad(u)
     stress=lam*trace(g)*kmat+mu*(g+transpose(g))
     ### ... get new acceleration ....
     #mypde.setValue(X=-stress)          
     #a=mypde.getSolution()
     ### ... get new displacement ...
     #u_p1=2*u-u_m1+h*h*a
###NEW WAY
     y = ((rho/(-rho-b*h))*(u_m1-2*u))+(((b*h)/(-rho-(b*h)))*-u)
     mypde.setValue(X=-stress*((h*h)/(-rho-h*b)),Y=y)
     u_p1 = mypde.getSolution()
     # ... shift displacements ....
     u_m1=u
     u=u_p1
     #stress = 
     t+=h
     n+=1
     print n,"-th time step t ",t
     u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)

#     print "u at point charge=",u_pc
     u_pc_x1 = u_pot[0,0]
     u_pc_y1 = u_pot[0,1]
     u_pc_x2 = u_pot[1,0]
     u_pc_y2 = u_pot[1,1]
     u_pc_x3 = u_pot[2,0]
     u_pc_y3 = u_pot[2,1]
           
     # save displacements at point source to file for t > 0
     u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
 
     # ... save current acceleration in units of gravity and displacements 
     #saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
     #displacement = length(u), tensor = stress, Ux = u[0] )
     saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)

   u_pc_data.close()

# joel wrote this to create directories for you
# It only needs to be this complicated to accomodate our
# automated tests
import os
def needdirs(dirlist):
    for name in dirlist:
    if name == '':
        continue
    if not os.path.exists(name):
           try:
        os.makedirs(name)
           except OSError:
        if not os.path.exists(save_path):
            raise
1	########################################################
2	#
3	# Copyright (c) 2009 by University of Queensland
4	# Earth Systems Science Computational Center (ESSCC)
5	# http://www.uq.edu.au/esscc
6	#
7	# Primary Business: Queensland, Australia
8	# Licensed under the Open Software License version 3.0
9	# http://www.opensource.org/licenses/osl-3.0.php
10	#
11	########################################################
12
13	__copyright__="""Copyright (c) 2009 by University of Queensland
14	Earth Systems Science Computational Center (ESSCC)
15	http://www.uq.edu.au/esscc
16	Primary Business: Queensland, Australia"""
17	__license__="""Licensed under the Open Software License version 3.0
18	http://www.opensource.org/licenses/osl-3.0.php"""
19	__url__="https://launchpad.net/escript-finley"
20
21	"""
22	A collection of routines to use in cookbook examples.
23	Author: Antony Hallam antony.hallam@uqconnect.edu.au
24	"""
25
26	from esys.pycad import CurveLoop
27
28	# numpy for array handling
29	import numpy as np
30
31	# tools for dealing with PDEs - contains locator
32	from esys.escript.pdetools import Locator, Projector
33
34	from esys.escript import *
35	from esys.escript.linearPDEs import LinearPDE
36
37	# routine to find consecutive coordinates of a loop in pycad
38	def getLoopCoords(loop):
39	# return all construction points of input
40	temp = loop.getConstructionPoints()
41	#create a numpy array for xyz components or construction points
42	coords = np.zeros([len(temp),3],float)
43	#place construction points in array
44	for i in range(0,len(temp)):
45	coords[i,:]=temp[i].getCoordinates()
46	#return a numpy array
47	return coords
48
49
50	########################################################
51	# subroutine: cbphones
52	# Allows us to record the values of a PDE at various
53	# specified locations in the model.
54	# Arguments:
55	# domain : domain of model
56	# U : Current time state displacement solution.
57	# phones : Geophone Locations
58	# dim : model dimesions
59	# savepath: where to output the data files local is default
60	########################################################
61	def cbphones(domain,U,phones,dim,savepath=""):
62	#find the number of geophones
63	nphones = len(phones)
64	u_pot = np.zeros([nphones,dim],float)
65
66	for i in range(0,nphones):
67	# define the location of the phone source
68	L=Locator(domain,numpy.array(phones[i]))
69	# find potential at point source.
70	temp = L.getValue(U)
71	for j in range(0,dim):
72	u_pot[i,j]=temp[j]
73
74	# open file to save displacement at point source
75	return u_pot
76
77	########################################################
78	# subroutine: wavesolver2d
79	# Can solve a generic 2D wave propagation problem with a
80	# point source in a homogeneous medium.
81	# Arguments:
82	# domain : domain to solve over
83	# h : time step
84	# tend : end time
85	# lam, mu : lames constants for domain
86	# rho : density of domain
87	# U0 : magnitude of source
88	# xc : source location in domain (Vector)
89	# savepath: where to output the data files
90	########################################################
91	def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath):
92	from esys.escript.linearPDEs import LinearPDE
93	x=domain.getX()
94	# ... open new PDE ...
95	mypde=LinearPDE(domain)
96	#mypde.setSolverMethod(LinearPDE.LUMPING)
97	mypde.setSymmetryOn()
98	kmat = kronecker(domain)
99	mypde.setValue(D=kmat*rho)
100
101	# define small radius around point xc
102	# Lsup(x) returns the maximum value of the argument x
103	src_radius = 50#2*Lsup(domain.getSize())
104	print "src_radius = ",src_radius
105
106	dunit=numpy.array([0.,1.]) # defines direction of point source
107
108
109	# ... set initial values ....
110	n=0
111	# initial value of displacement at point source is constant (U0=0.01)
112	# for first two time steps
113	u=U0(cos(length(x-xc)3.1415/src_radius)+1)whereNegative(length(x-xc)-src_radius)dunit
114	u_m1=u
115	t=0
116
117	u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
118	u_pc_x1 = u_pot[0,0]
119	u_pc_y1 = u_pot[0,1]
120	u_pc_x2 = u_pot[1,0]
121	u_pc_y2 = u_pot[1,1]
122	u_pc_x3 = u_pot[2,0]
123	u_pc_y3 = u_pot[2,1]
124
125	# open file to save displacement at point source
126	u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
127	u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
128
129	while t<tend:
130	# ... get current stress ....
131
132	##OLD WAY
133	g=grad(u)
134	stress=lamtrace(g)kmat+mu*(g+transpose(g))
135	### ... get new acceleration ....
136	#mypde.setValue(X=-stress)
137	#a=mypde.getSolution()
138	### ... get new displacement ...
139	#u_p1=2u-u_m1+hh*a
140	###NEW WAY
141	mypde.setValue(X=-stress(hh),Y=(rho2u-rho*u_m1))
142	u_p1 = mypde.getSolution()
143	# ... shift displacements ....
144	u_m1=u
145	u=u_p1
146	#stress =
147	t+=h
148	n+=1
149	print n,"-th time step t ",t
150	u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
151
152	# print "u at point charge=",u_pc
153	u_pc_x1 = u_pot[0,0]
154	u_pc_y1 = u_pot[0,1]
155	u_pc_x2 = u_pot[1,0]
156	u_pc_y2 = u_pot[1,1]
157	u_pc_x3 = u_pot[2,0]
158	u_pc_y3 = u_pot[2,1]
159
160	# save displacements at point source to file for t > 0
161	u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
162
163	# ... save current acceleration in units of gravity and displacements
164	#saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
165	#displacement = length(u), tensor = stress, Ux = u[0] )
166	saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
167
168	u_pc_data.close()
169
170
171	########################################################
172	# subroutine: wavesolver2d
173	# Can solve a generic 2D wave propagation problem with a
174	# point source in a homogeneous medium with friction.
175	# Arguments:
176	# domain : domain to solve over
177	# h : time step
178	# tend : end time
179	# lam, mu : lames constants for domain
180	# rho : density of domain
181	# U0 : magnitude of source
182	# xc : source location in domain (Vector)
183	# savepath: where to output the data files
184	########################################################
185	def wavesolver2df(domain,h,tend,lam,mu,rho,U0,xc,savepath):
186	x=domain.getX()
187	# ... open new PDE ...
188	mypde=LinearPDE(domain)
189	#mypde.setSolverMethod(LinearPDE.LUMPING)
190	mypde.setSymmetryOn()
191	kmat = kronecker(domain)
192	mypde.setValue(D=kmat)
193	b=0.9
194
195	# define small radius around point xc
196	# Lsup(x) returns the maximum value of the argument x
197	src_radius = 50#2*Lsup(domain.getSize())
198	print "src_radius = ",src_radius
199
200	dunit=numpy.array([0.,1.]) # defines direction of point source
201
202
203	# ... set initial values ....
204	n=0
205	# initial value of displacement at point source is constant (U0=0.01)
206	# for first two time steps
207	u=U0(cos(length(x-xc)3.1415/src_radius)+1)whereNegative(length(x-xc)-src_radius)dunit
208	u_m1=u
209	t=0
210
211	u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
212	u_pc_x1 = u_pot[0,0]
213	u_pc_y1 = u_pot[0,1]
214	u_pc_x2 = u_pot[1,0]
215	u_pc_y2 = u_pot[1,1]
216	u_pc_x3 = u_pot[2,0]
217	u_pc_y3 = u_pot[2,1]
218
219	# open file to save displacement at point source
220	u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
221	u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
222
223	while t<tend:
224	# ... get current stress ....
225
226	##OLD WAY
227	g=grad(u)
228	stress=lamtrace(g)kmat+mu*(g+transpose(g))
229	### ... get new acceleration ....
230	#mypde.setValue(X=-stress)
231	#a=mypde.getSolution()
232	### ... get new displacement ...
233	#u_p1=2u-u_m1+hh*a
234	###NEW WAY
235	y = ((rho/(-rho-bh))(u_m1-2u))+(((bh)/(-rho-(bh)))-u)
236	mypde.setValue(X=-stress((hh)/(-rho-h*b)),Y=y)
237	u_p1 = mypde.getSolution()
238	# ... shift displacements ....
239	u_m1=u
240	u=u_p1
241	#stress =
242	t+=h
243	n+=1
244	print n,"-th time step t ",t
245	u_pot = cbphones(domain,u,[[300.,200.],[500.,200.],[750.,200.]],2)
246
247	# print "u at point charge=",u_pc
248	u_pc_x1 = u_pot[0,0]
249	u_pc_y1 = u_pot[0,1]
250	u_pc_x2 = u_pot[1,0]
251	u_pc_y2 = u_pot[1,1]
252	u_pc_x3 = u_pot[2,0]
253	u_pc_y3 = u_pot[2,1]
254
255	# save displacements at point source to file for t > 0
256	u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
257
258	# ... save current acceleration in units of gravity and displacements
259	#saveVTK(os.path.join(savepath,"usoln.%i.vtu"%n),acceleration=length(a)/9.81,
260	#displacement = length(u), tensor = stress, Ux = u[0] )
261	saveVTK(os.path.join(savepath,"tonysol.%i.vtu"%n),output1 = length(u),tensor=stress)
262
263	u_pc_data.close()
264
265	# joel wrote this to create directories for you
266	# It only needs to be this complicated to accomodate our
267	# automated tests
268	import os
269	def needdirs(dirlist):
270	for name in dirlist:
271	if name == '':
272	continue
273	if not os.path.exists(name):
274	try:
275	os.makedirs(name)
276	except OSError:
277	if not os.path.exists(save_path):
278	raise