# Contents of /trunk/doc/examples/cookbook/cblib1.py

Revision 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 9620 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.

 1 ######################################################## 2 # 3 # Copyright (c) 2009-2010 by University of Queensland 4 # Earth Systems Science Computational Center (ESSCC) 5 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 10 # 11 ######################################################## 12 13 __copyright__="""Copyright (c) 2009-2010 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 19 __url__= 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 from esys.weipa import saveVTK 37 38 39 # routine to find consecutive coordinates of a loop in pycad 40 def getLoopCoords(loop): 41 # return all construction points of input 42 temp = loop.getConstructionPoints() 43 #create a numpy array for xyz components or construction points 44 coords = np.zeros([len(temp),3],float) 45 #place construction points in array 46 for i in range(0,len(temp)): 47 coords[i,:]=temp[i].getCoordinates() 48 #return a numpy array 49 return coords 50 51 52 ######################################################## 53 # subroutine: cbphones 54 # Allows us to record the values of a PDE at various 55 # specified locations in the model. 56 # Arguments: 57 # domain : domain of model 58 # U : Current time state displacement solution. 59 # phones : Geophone Locations 60 # dim : model dimesions 61 # savepath: where to output the data files local is default 62 ######################################################## 63 def cbphones(domain,U,phones,dim,savepath=""): 64 #find the number of geophones 65 nphones = len(phones) 66 u_pot = np.zeros([nphones,dim],float) 67 68 for i in range(0,nphones): 69 # define the location of the phone source 70 L=Locator(domain,numpy.array(phones[i])) 71 # find potential at point source. 72 temp = L.getValue(U) 73 for j in range(0,dim): 74 u_pot[i,j]=temp[j] 75 76 # open file to save displacement at point source 77 return u_pot 78 79 ######################################################## 80 # subroutine: wavesolver2d 81 # Can solve a generic 2D wave propagation problem with a 82 # point source in a homogeneous medium. 83 # Arguments: 84 # domain : domain to solve over 85 # h : time step 86 # tend : end time 87 # lam, mu : lames constants for domain 88 # rho : density of domain 89 # U0 : magnitude of source 90 # xc : source location in domain (Vector) 91 # savepath: where to output the data files 92 ######################################################## 93 def wavesolver2d(domain,h,tend,lam,mu,rho,U0,xc,savepath,output="vtk"): 94 from esys.escript.linearPDEs import LinearPDE 95 x=domain.getX() 96 # ... open new PDE ... 97 mypde=LinearPDE(domain) 98 #mypde.setSolverMethod(LinearPDE.LUMPING) 99 mypde.setSymmetryOn() 100 kmat = kronecker(domain) 101 mypde.setValue(D=kmat*rho) 102 103 # define small radius around point xc 104 # Lsup(x) returns the maximum value of the argument x 105 src_radius = 50#2*Lsup(domain.getSize()) 106 print("src_radius = ",src_radius) 107 108 dunit=numpy.array([0.,1.]) # defines direction of point source 109 110 111 # ... set initial values .... 112 n=0 113 # initial value of displacement at point source is constant (U0=0.01) 114 # for first two time steps 115 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radius)*dunit 116 u_m1=u 117 t=0 118 119 u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2) 120 u_pc_x1 = u_pot[0,0] 121 u_pc_y1 = u_pot[0,1] 122 u_pc_x2 = u_pot[1,0] 123 u_pc_y2 = u_pot[1,1] 124 u_pc_x3 = u_pot[2,0] 125 u_pc_y3 = u_pot[2,1] 126 127 # open file to save displacement at point source 128 u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w') 129 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)) 130 131 # while t

 ViewVC Help Powered by ViewVC 1.1.26