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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (hide annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 7 months ago) by sshaw
File MIME type: text/x-python
File size: 4187 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 jfenwick 3981 ##############################################################################
2 ahallam 3405 #
3 jfenwick 5593 # Copyright (c) 2009-2015 by The University of Queensland
4 jfenwick 3981 # http://www.uq.edu.au
5 ahallam 3405 #
6     # Primary Business: Queensland, Australia
7     # Licensed under the Open Software License version 3.0
8     # http://www.opensource.org/licenses/osl-3.0.php
9     #
10 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
12     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 #
14     ##############################################################################
15 sshaw 5707 from __future__ import division, print_function
16 ahallam 3405
17 jfenwick 5593 __copyright__="""Copyright (c) 2009-2015 by The University of Queensland
18 jfenwick 3981 http://www.uq.edu.au
19 ahallam 3405 Primary Business: Queensland, Australia"""
20     __license__="""Licensed under the Open Software License version 3.0
21     http://www.opensource.org/licenses/osl-3.0.php"""
22     __url__="https://launchpad.net/escript-finley"
23    
24     """
25     Author: Antony Hallam antony.hallam@uqconnect.edu.au
26     """
27    
28    
29     ############################################################FILE HEADER
30     # example10m.py
31     # Create a simple 2D mesh, which is optimised for cells close to the
32     # source. Larger elements are used to decrease computational requirements
33     # and to properly fullfil the boundary conditions.
34     #
35     #######################################################EXTERNAL MODULES
36     from esys.pycad import * #domain constructor
37     from esys.pycad.extras import layer_cake
38     from esys.pycad.gmsh import Design #Finite Element meshing package
39 sshaw 5288
40     try:
41     from esys.finley import MakeDomain
42     HAVE_FINLEY = True
43     except ImportError:
44     print("Finley module not available")
45     HAVE_FINLEY = False
46 ahallam 3405 from esys.escript import mkDir, getMPISizeWorld
47     import os
48     import subprocess as sp
49     ########################################################MPI WORLD CHECK
50     if getMPISizeWorld() > 1:
51 sshaw 4576 import sys
52     print("This example will not run in an MPI world.")
53     sys.exit(0)
54 ahallam 3405
55 sshaw 5288 if HAVE_FINLEY:
56     # make sure path exists
57     save_path= os.path.join("data","example10m")
58     mkDir(save_path)
59 ahallam 3405
60 sshaw 5288 ################################################BIG DOMAIN
61     #ESTABLISHING PARAMETERS
62     width=10000. #width of model
63     depth=10000. #depth of model
64     bele_size=500. #big element size
65     #DOMAIN CONSTRUCTION
66     p0=Point(0.0, 0.0)
67     p1=Point(width, 0.0)
68     p2=Point(width, depth)
69     p3=Point(0.0, depth)
70     # Join corners in anti-clockwise manner.
71     l01=Line(p0, p1)
72     l12=Line(p1, p2)
73     l23=Line(p2, p3)
74     l30=Line(p3, p0)
75 ahallam 3405
76 sshaw 5288 cbig=CurveLoop(l01,l12,l23,l30)
77 ahallam 3405
78 sshaw 5288 ################################################SMALL DOMAIN
79     #ESTABLISHING PARAMETERS
80     xwidth=1000.0 #x width of model
81     zdepth=1000.0 #y width of model
82     sele_size=10. #small element size
83     #TRANSFORM
84     xshift=width/2-xwidth/2
85     zshift=depth/2-zdepth/2
86     #DOMAIN CONSTRUCTION
87     p4=Point(xshift, zshift)
88     p5=Point(xwidth+xshift, zshift)
89     p6=Point(xwidth+xshift, zdepth+zshift)
90     p7=Point(xshift, zdepth+zshift)
91     # Join corners in anti-clockwise manner.
92     l45=Line(p4, p5)
93     l56=Line(p5, p6)
94     l67=Line(p6, p7)
95     l74=Line(p7, p4)
96 ahallam 3405
97 sshaw 5288 csmall=CurveLoop(l45,l56,l67,l74)
98 ahallam 3405
99 sshaw 5288 ssmall=PlaneSurface(csmall)
100     sbig=PlaneSurface(cbig,holes=[csmall])
101 ahallam 3405
102 sshaw 5288 #############################################EXPORTING MESH FOR ESCRIPT
103     # Design the geometry for the big mesh.
104     d1=Design(dim=2, element_size=bele_size, order=1)
105     d1.addItems(sbig)
106     d1.addItems(PropertySet(l01,l12,l23,l30))
107     d1.setScriptFileName(os.path.join(save_path,"example10m_big.geo"))
108     MakeDomain(d1)
109 ahallam 3405
110 sshaw 5288 # Design the geometry for the small mesh.
111     d2=Design(dim=2, element_size=sele_size, order=1)
112     d2.addItems(ssmall)
113     d2.setScriptFileName(os.path.join(save_path,"example10m_small.geo"))
114     MakeDomain(d2)
115 ahallam 3405
116 sshaw 5288 # Join the two meshes using Gmsh and then apply a 2D meshing algorithm.
117     # The small mesh must come before the big mesh in the merging call!!@!!@!
118     sp.call("gmsh -2 "+
119     os.path.join(save_path,"example10m_small.geo")+" "+
120     os.path.join(save_path,"example10m_big.geo")+" -o "+
121     os.path.join(save_path,"example10m.msh"),shell=True)

  ViewVC Help
Powered by ViewVC 1.1.26