/[escript]/trunk/doc/examples/cookbook/example11b.py
ViewVC logotype

Annotation of /trunk/doc/examples/cookbook/example11b.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3981 - (hide annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years ago) by jfenwick
File MIME type: text/x-python
File size: 3437 byte(s)
First pass of updating copyright notices
1 ahallam 3427
2 jfenwick 3981 ##############################################################################
3 ahallam 3427 #
4 jfenwick 3911 # Copyright (c) 2009-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ahallam 3427 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 ahallam 3427
16 jfenwick 3911 __copyright__="""Copyright (c) 2009-2012 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 ahallam 3427 Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23     """
24     Author: Antony Hallam antony.hallam@uqconnect.edu.au
25     """
26    
27     ############################################################FILE HEADER
28     # example10a.py
29     # Model of gravitational Potential.
30    
31     #######################################################EXTERNAL MODULES
32     # To solve the problem it is necessary to import the modules we require.
33     from esys.escript import * # This imports everything from the escript library
34     from esys.escript.unitsSI import *
35     from esys.escript.linearPDEs import LinearPDE # This defines LinearPDE as LinearPDE
36 gross 3439 from esys.weipa import saveVTK
37 ahallam 3427 from esys.finley import Rectangle # This imports the rectangle domain function from finley
38     import os, sys #This package is necessary to handle saving our data.
39     from math import pi, sqrt, sin, cos
40    
41     from esys.escript.pdetools import Projector
42    
43     import matplotlib
44     matplotlib.use('agg') #It's just here for automated testing
45    
46 gross 3439 from cblib import toRegGrid
47    
48    
49 ahallam 3427 import pylab as pl #Plotting package
50     import numpy as np
51    
52     ########################################################MPI WORLD CHECK
53     if getMPISizeWorld() > 1:
54     import sys
55 jfenwick 3892 print("This example will not run in an MPI world.")
56 ahallam 3427 sys.exit(0)
57    
58     #################################################ESTABLISHING VARIABLES
59     #Domain related.
60     mx = 1000*m #meters - model length
61     my = -250*m #meters - model depth
62     ndx = 200 # mesh steps in x direction
63     ndy = 200 # mesh steps in y direction - one dimension means one element
64     #PDE related
65     res1=500.0
66     res2=10.0
67     res3=10000.0
68     #con=1/res
69     cur=1000.
70    
71     ################################################ESTABLISHING PARAMETERS
72     #the folder to put our outputs in, leave blank "" for script path
73     save_path= os.path.join("data","example11")
74     #ensure the dir exists
75     mkDir(save_path)
76    
77     ####################################################DOMAIN CONSTRUCTION
78     domain = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy)
79     x=Solution(domain).getX()
80    
81     kro=kronecker(domain)
82     source1=[3.*mx/8.,0]; source2=[5.*mx/8.,0]
83    
84     c1=length(exp(-length(x-source1)/(10.))); c1=c1/integrate(c1)
85     c2=-length(exp(-length(x-source2)/(10.))); c2=c2/integrate(c2)
86     sourceg=cur*(c1-c2)
87    
88     res=res1*wherePositive(x[1]-my/3)+res2*whereNegative(x[1]-my/3)*wherePositive(x[1]-my*2/3)+res3*whereNegative(x[1]-my*2/3)
89     con=1/res
90     q=whereZero(x[1]-my)+whereZero(x[0])+whereZero(x[0]-mx)
91     ###############################################ESCRIPT PDE CONSTRUCTION
92    
93     mypde=LinearPDE(domain)
94     mypde.setValue(A=con*kro,Y=sourceg,q=q,r=0)
95     #mypde.setSymmetryOn()
96     sol=mypde.getSolution()
97    
98     # Save the output to file.
99     saveVTK(os.path.join(save_path,"ex11b.vtu"),\
100     source=sourceg,\
101     res_pot=sol,\
102     res=res,\
103     curden=-con*grad(sol),\
104     efield=-grad(sol))

  ViewVC Help
Powered by ViewVC 1.1.26