/[escript]/branches/symbolic_from_3470/escript/py_src/symbolic/utils.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/escript/py_src/symbolic/utils.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3864 - (show annotations)
Mon Mar 12 05:18:16 2012 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/x-python
File size: 6461 byte(s)
Symbols now allow direct operations with Data objects and grad() et al allow
specifying FunctionSpace objects directly, without having to use temporary
symbols :-)


1
2 ########################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21 __author__="Cihan Altinay"
22
23 """
24 :var __author__: name of author
25 :var __copyright__: copyrights
26 :var __license__: licence agreement
27 :var __url__: url entry point on documentation
28 :var __version__: version
29 :var __date__: date of the version
30 """
31
32 import numpy
33 import sympy
34 from symbols import Symbol
35
36 def symbols(*names, **kwargs):
37 """
38 Emulates the behaviour of sympy.symbols.
39 """
40
41 shape=kwargs.pop('shape', ())
42
43 s = names[0]
44 if not isinstance(s, list):
45 import re
46 s = re.split('\s|,', s)
47 res = []
48 for t in s:
49 # skip empty strings
50 if not t:
51 continue
52 sym = Symbol(t, shape, **kwargs)
53 res.append(sym)
54 res = tuple(res)
55 if len(res) == 0: # var('')
56 res = None
57 elif len(res) == 1: # var('x')
58 res = res[0]
59 # otherwise var('a b ...')
60 return res
61
62 def combineData(array, shape):
63 """
64 """
65
66 # array could just be a single value
67 if not hasattr(array,'__len__') and shape==():
68 return array
69
70 from esys.escript import Data
71 n=numpy.array(array) # for indexing
72
73 # find function space if any
74 dom=set()
75 fs=set()
76 for idx in numpy.ndindex(shape):
77 if isinstance(n[idx], Data):
78 fs.add(n[idx].getFunctionSpace())
79 dom.add(n[idx].getDomain())
80
81 if len(dom)>1:
82 domain=dom.pop()
83 while len(dom)>0:
84 if domain!=dom.pop():
85 raise ValueError("Mixing of domains not supported")
86
87 if len(fs)>0:
88 d=Data(0., shape, fs.pop()) # maybe interpolate instead of using first?
89 else:
90 d=numpy.zeros(shape)
91 for idx in numpy.ndindex(shape):
92 #z=numpy.zeros(shape)
93 #z[idx]=1.
94 #d+=n[idx]*z # much slower!
95 if hasattr(n[idx], "ndim") and n[idx].ndim==0:
96 d[idx]=float(n[idx])
97 else:
98 d[idx]=n[idx]
99 return d
100
101 def isSymbol(arg):
102 """
103 Returns True if the argument ``arg`` is an escript ``Symbol`` or
104 ``sympy.Basic`` object, False otherwise.
105 """
106 return isinstance(arg, Symbol) or isinstance(arg, sympy.Basic)
107
108 def removeFsFromGrad(sym):
109 """
110 Returns ``sym`` with all occurrences grad_n(a,b,c) replaced by grad_n(a,b).
111 That is, all functionspace parameters are removed.
112 """
113 from esys.escript import symfn
114 gg=sym.atoms(symfn.grad_n)
115 for g in gg:
116 if len(g.args)==3:
117 r=symfn.grad_n(*g.args[:2])
118 sym=sym.subs(g, r)
119 return sym
120
121 def getTotalDifferential(f, x, order=0):
122 """
123 Df/Dx = del_f/del_x + del_f/del_grad(x)*del_grad(x)/del_x + ...
124 \ / \ /
125 a b
126 """
127
128 from esys.escript import util
129 res=()
130 shape=util.getShape(f)
131 if not isSymbol(f):
132 res+=(numpy.zeros(shape+x.getShape()),)
133 for i in range(order):
134 x=x.grad()
135 res+=numpy.zeros(shape+x.getShape())
136
137 elif x.getRank()==0:
138 f=removeFsFromGrad(f)
139 dfdx=f.diff(x)
140 dgdx=x.grad().diff(x)
141 a=numpy.empty(shape, dtype=object)
142 if order>0:
143 b=numpy.empty(shape+dgdx.getShape(), dtype=object)
144
145 if len(shape)==0:
146 for j in numpy.ndindex(dgdx.getShape()):
147 y=dfdx
148 z=dgdx[j]
149 # expand() and coeff() are very expensive so
150 # we set the unwanted factors to zero to extract
151 # the one we need
152 for jj in numpy.ndindex(dgdx.getShape()):
153 if j==jj: continue
154 y=y.subs(dgdx[jj], 0)
155 a=y.subs(z,0) # terms in x and constants
156 if order>0:
157 b[j]=y.subs(z,1)-a
158 else:
159 for i in numpy.ndindex(shape):
160 for j in numpy.ndindex(dgdx.getShape()):
161 y=dfdx[i]
162 z=dgdx[j]
163 for jj in numpy.ndindex(dgdx.getShape()):
164 if j==jj: continue
165 y=y.subs(dgdx[jj], 0)
166 a[i]=y.subs(z,0) # terms in x and constants
167 if order>0:
168 b[i+j]=y.subs(z,1)-a[i]
169 res+=(Symbol(a, dim=f.getDim(), subs=f.getDataSubstitutions()),)
170 if order>0:
171 res+=(Symbol(b, dim=f.getDim(), subs=f.getDataSubstitutions()),)
172
173 elif x.getRank()==1:
174 f=removeFsFromGrad(f)
175 dfdx=f.diff(x)
176 dgdx=x.grad().diff(x).transpose(2)
177 a=numpy.empty(shape+x.getShape(), dtype=object)
178 if order>0:
179 b=numpy.empty(shape+x.grad().getShape(), dtype=object)
180
181 if len(shape)==0:
182 raise NotImplementedError('f scalar, x vector')
183 else:
184 for i in numpy.ndindex(shape):
185 for k,l in numpy.ndindex(x.grad().getShape()):
186 if dgdx[k,k,l]==0:
187 a[i+(k,)]=0
188 if order>0:
189 b[i+(k,l)]=0
190 else:
191 y=dfdx[i+(k,)]
192 z=dgdx[k,k,l]
193 for kk,ll in numpy.ndindex(x.grad().getShape()):
194 if k==kk and l==ll: continue
195 y=y.subs(dgdx[kk,kk,ll], 0)
196 a[i+(k,)]=y.subs(z,0) # terms in x and constants
197 if order>0:
198 b[i+(k,l)]=y.subs(z,1)-a[i+(k,)]
199
200 res+=(Symbol(a, dim=f.getDim(), subs=f.getDataSubstitutions()),)
201 if order>0:
202 res+=(Symbol(b, dim=f.getDim(), subs=f.getDataSubstitutions()),)
203
204 if len(res)==1:
205 return res[0]
206 else:
207 return res
208

  ViewVC Help
Powered by ViewVC 1.1.26