pycad/py_src/extras.py

# -*- coding: utf-8 -*-

########################################################
#
# Copyright (c) 2003-2012 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) 2003-2012 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"

"""
Layer Cake

This script uses the pycad module to build a rectangular layered
system of prisms for modelling.

:var __author__: name of author
:var __copyright__: copyrights
:var __license__: licence agreement
:var __url__: url entry point on documentation
:var __version__: version
:var __date__: date of the version
"""

__author__="Antony Hallam antony.hallam@uqconnect.edu.au"

from esys.pycad import * #domain constructor
from esys.pycad.gmsh import Design #Finite Element meshing package

def buildFreeSurface(xwidth,ywidth):
    '''
    Build a free surface to start the layer cake model.
        This surface is planar.

    Parameters:
        xwidth :: width in x direction in meters.\
        ywidth :: width in y direction in meters.\
    '''

    # Layer Corners
    corner_points=[]
    corner_points.append(Point(0.0,    0.0,      0.0))
    corner_points.append(Point(xwidth, 0.0,      0.0))
    corner_points.append(Point(xwidth, ywidth,   0.0))
    corner_points.append(Point(0.0,    ywidth,   0.0))
    corner_points.append(corner_points[0]) #repeated point for line looping

    # Edges of Free Surface 
    hor_lines=[]
    for i in range(0,4): # loop through four sides
        hor_lines.append(Line(corner_points[i],corner_points[i+1]))

    # Create Free Surface
    free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))

    # Return Surface and primative arrays.
    return free_surf,hor_lines,corner_points

def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
    '''
    Builds a boxlike volume and returns primatives for layered model
    construction.

    Parameters:
        xwidth :: width in x direction in meters.\
        ywidth :: width in y direction in meters.\
        depth  :: depth to bottom of layer in meters.\
        lay_surf :: surface at top of layer (see buildFreeSurf)\
        hor_lines :: lines of lay_surf\
        corner_points :: points of hor_lines\
    '''
    # Layer Corners
    corner_points.append(Point(0.0,    0.0,    depth))
    corner_points.append(Point(xwidth, 0.0,    depth))
    corner_points.append(Point(xwidth, ywidth, depth))
    corner_points.append(Point(0.0,    ywidth, depth))
    corner_points.append(corner_points[5]) #repeated point for line looping

    # Build the bottom surface edges.
    for i in range(0,4): # loop through four edges
        hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))

    # Join corners vertically.
    ver_lines=[]
    for i in range(0,4): # loop through four corners
        ver_lines.append(Line(corner_points[i],corner_points[i+5]))
    ver_lines.append(ver_lines[0]) #repeated edge for surface looping   

    # Build surface array.
    lay_surf=[-lay_surf] #Negative of top surface
    # Bottom Surface
    lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
    for i in range(0,4): # loop through four sides
        lay_surf.append(PlaneSurface(CurveLoop(\
                        -ver_lines[i],-hor_lines[i+4],\
                        ver_lines[i+1],hor_lines[i] )))

    # Build Layer Volume
    lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
    
    # Return new volume, and primatives for next volume layer.    
    return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]


def layer_cake(domain,xwidth,ywidth,depths):
    '''
    Builds a horizontally layered box like model. All layers are 
    tagged as 'intface_i' where i is the python style integer denoting
    that layer. For example, the free surface is tagged 'interface_0'.
    Volumes are similarly tagged as 'volume_i'.

    Parameters:
       domain :: output of Pycad.Design - It needs to be dim 3.\
       xwidth :: width in x direction in meters.\
       ywidth :: width in y direction in meters.\
       depth  :: depth to bottom of layer in meters.\

    One may save the domain using:
        # Output settings.\    
        domain.setScriptFileName(os.path.join(save_path,fname+".geo"))\
        domain.setMeshFileName(os.path.join(save_path,fname+".msh"))\
        findomain=fin.MakeDomain(domain) #  make the finley domain:\
    
        Create a file that can be read back in to python with ReadMesh.\
        findomain.write(os.path.join(save_path,fname+".fly"))\             
    '''
    
    #get number of layers
    if not hasattr(depths,"__len__"): depths = [ depths, ]
    ndepths=len(depths)

    if domain.getDim() != 3:
        raise TypeError("domain must be of dimension order 3.")        

    # Build the First Surface and add it to the domain
    fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
    domain.addItems(PropertySet('intface_%d'%(0),fsuf))        

    # Build each layer sequentially
    # Set up temporary variables.
    tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
    # Loop through all depths.
    for i in range(0,ndepths):
        # Build new layer.
        tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
                                     tsuf,tsurl,tsurp)
        # Add the new interface to the domain.
        domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))        
        # Add the new volume/layer to the domain.
        domain.addItems(PropertySet('volume_%d'%i,tvol))
    
    return domain

1	# -- coding: utf-8 --
2
3	########################################################
4	#
5	# Copyright (c) 2003-2012 by University of Queensland
6	# Earth Systems Science Computational Center (ESSCC)
7	# http://www.uq.edu.au/esscc
8	#
9	# Primary Business: Queensland, Australia
10	# Licensed under the Open Software License version 3.0
11	# http://www.opensource.org/licenses/osl-3.0.php
12	#
13	########################################################
14
15	__copyright__="""Copyright (c) 2003-2012 by University of Queensland
16	Earth Systems Science Computational Center (ESSCC)
17	http://www.uq.edu.au/esscc
18	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	Layer Cake
25
26	This script uses the pycad module to build a rectangular layered
27	system of prisms for modelling.
28
29	:var __author__: name of author
30	:var __copyright__: copyrights
31	:var __license__: licence agreement
32	:var __url__: url entry point on documentation
33	:var __version__: version
34	:var __date__: date of the version
35	"""
36
37	__author__="Antony Hallam antony.hallam@uqconnect.edu.au"
38
39	from esys.pycad import * #domain constructor
40	from esys.pycad.gmsh import Design #Finite Element meshing package
41
42	def buildFreeSurface(xwidth,ywidth):
43	'''
44	Build a free surface to start the layer cake model.
45	This surface is planar.
46
47	Parameters:
48	xwidth :: width in x direction in meters.\
49	ywidth :: width in y direction in meters.\
50	'''
51
52	# Layer Corners
53	corner_points=[]
54	corner_points.append(Point(0.0, 0.0, 0.0))
55	corner_points.append(Point(xwidth, 0.0, 0.0))
56	corner_points.append(Point(xwidth, ywidth, 0.0))
57	corner_points.append(Point(0.0, ywidth, 0.0))
58	corner_points.append(corner_points[0]) #repeated point for line looping
59
60	# Edges of Free Surface
61	hor_lines=[]
62	for i in range(0,4): # loop through four sides
63	hor_lines.append(Line(corner_points[i],corner_points[i+1]))
64
65	# Create Free Surface
66	free_surf = PlaneSurface(CurveLoop(*tuple(hor_lines[0:4])))
67
68	# Return Surface and primative arrays.
69	return free_surf,hor_lines,corner_points
70
71	def buildLayer(xwidth,ywidth,depth,lay_surf,hor_lines,corner_points):
72	'''
73	Builds a boxlike volume and returns primatives for layered model
74	construction.
75
76	Parameters:
77	xwidth :: width in x direction in meters.\
78	ywidth :: width in y direction in meters.\
79	depth :: depth to bottom of layer in meters.\
80	lay_surf :: surface at top of layer (see buildFreeSurf)\
81	hor_lines :: lines of lay_surf\
82	corner_points :: points of hor_lines\
83	'''
84	# Layer Corners
85	corner_points.append(Point(0.0, 0.0, depth))
86	corner_points.append(Point(xwidth, 0.0, depth))
87	corner_points.append(Point(xwidth, ywidth, depth))
88	corner_points.append(Point(0.0, ywidth, depth))
89	corner_points.append(corner_points[5]) #repeated point for line looping
90
91	# Build the bottom surface edges.
92	for i in range(0,4): # loop through four edges
93	hor_lines.append(Line(corner_points[5+i],corner_points[6+i]))
94
95	# Join corners vertically.
96	ver_lines=[]
97	for i in range(0,4): # loop through four corners
98	ver_lines.append(Line(corner_points[i],corner_points[i+5]))
99	ver_lines.append(ver_lines[0]) #repeated edge for surface looping
100
101	# Build surface array.
102	lay_surf=[-lay_surf] #Negative of top surface
103	# Bottom Surface
104	lay_surf.append(PlaneSurface(CurveLoop(*tuple(hor_lines[4:8]))))
105	for i in range(0,4): # loop through four sides
106	lay_surf.append(PlaneSurface(CurveLoop(\
107	-ver_lines[i],-hor_lines[i+4],\
108	ver_lines[i+1],hor_lines[i] )))
109
110	# Build Layer Volume
111	lay_vol=Volume(-SurfaceLoop(*tuple(lay_surf)))
112
113	# Return new volume, and primatives for next volume layer.
114	return lay_vol,-lay_surf[1],hor_lines[4:8],corner_points[5:10]
115
116
117	def layer_cake(domain,xwidth,ywidth,depths):
118	'''
119	Builds a horizontally layered box like model. All layers are
120	tagged as 'intface_i' where i is the python style integer denoting
121	that layer. For example, the free surface is tagged 'interface_0'.
122	Volumes are similarly tagged as 'volume_i'.
123
124	Parameters:
125	domain :: output of Pycad.Design - It needs to be dim 3.\
126	xwidth :: width in x direction in meters.\
127	ywidth :: width in y direction in meters.\
128	depth :: depth to bottom of layer in meters.\
129
130	One may save the domain using:
131	# Output settings.\
132	domain.setScriptFileName(os.path.join(save_path,fname+".geo"))\
133	domain.setMeshFileName(os.path.join(save_path,fname+".msh"))\
134	findomain=fin.MakeDomain(domain) # make the finley domain:\
135
136	Create a file that can be read back in to python with ReadMesh.\
137	findomain.write(os.path.join(save_path,fname+".fly"))\
138	'''
139
140	#get number of layers
141	if not hasattr(depths,"__len__"): depths = [ depths, ]
142	ndepths=len(depths)
143
144	if domain.getDim() != 3:
145	raise TypeError("domain must be of dimension order 3.")
146
147	# Build the First Surface and add it to the domain
148	fsuf,fsurl,fsurp=buildFreeSurface(xwidth,ywidth)
149	domain.addItems(PropertySet('intface_%d'%(0),fsuf))
150
151	# Build each layer sequentially
152	# Set up temporary variables.
153	tsuf=fsuf; tsurl=fsurl; tsurp=fsurp
154	# Loop through all depths.
155	for i in range(0,ndepths):
156	# Build new layer.
157	tvol,tsuf,tsurl,tsurp=buildLayer(xwidth,ywidth,depths[i],\
158	tsuf,tsurl,tsurp)
159	# Add the new interface to the domain.
160	domain.addItems(PropertySet('intface_%d'%(i+1),tsuf))
161	# Add the new volume/layer to the domain.
162	domain.addItems(PropertySet('volume_%d'%i,tvol))
163
164	return domain
165