← Back to branch summary

~chaffra/ufl/main-old

« back to all changes in this revision

Viewing changes to ufl/algorithms/pdiffs.py

  • Committer: Chaffra Affouda
  • Date: 2012-06-26 02:44:28 UTC
  • mfrom: (1170.1.257 scratch)
  • Revision ID: chaffra@gmail.com-20120626024428-g3py2piveuv0ssjg
merge

Show diffs side-by-side

added added

removed removed

Lines of Context:
ufl/algorithms/pdiffs.py
39
39
        #self._spatial_dim = spatial_dim
40
40
    def __init__(self):
41
41
        MultiFunction.__init__(self)
42
 
    
 
42
 
43
43
    # TODO: Make sure we have implemented partial derivatives of all operators.
44
44
    #        At least non-compound ones should be covered, but compound ones
45
45
    #        may be a good idea in future versions.
46
 
    
 
46
 
47
47
    def expr(self, o):
48
48
        error("No partial derivative defined for %s" % type(o))
49
 
    
 
49
 
50
50
    # --- Basic algebra operators
51
 
    
 
51
 
52
52
    def index_sum(self, f):
53
53
        "d/dx sum_j x = TODO"
54
54
        TODO
55
 
    
 
55
 
56
56
    def sum(self, f):
57
57
        "d/dx_i sum_j x_j = 1"
58
58
        #_1 = IntValue(1, o.free_indices(), o.index_dimensions())
59
59
        _1 = IntValue(1) # TODO: Handle non-scalars
60
60
        return (_1,)*len(f.operands())
61
 
    
 
61
 
62
62
    def product(self, f):
63
63
        a, b = f.operands() # TODO: Assuming binary operator for now
64
64
        da = b # TODO: Is this right even for non-scalar b?
65
65
        db = a
66
66
        return (da, db)
67
 
    
 
67
 
68
68
    def division(self, f):
69
69
        """f = x/y
70
70
        d/dx x/y = 1/y
75
75
        ufl_assert(y.shape() == (), "Expecting scalars in division.")
76
76
        d = 1 / y
77
77
        return (d, -f*d)
78
 
    
 
78
 
79
79
    def power(self, f):
80
80
        """f = x**y
81
81
        d/dx x**y = y*x**(y-1) = y*f/x
84
84
        dx = y*f/x
85
85
        dy = ln(x)*f
86
86
        return (dx, dy)
87
 
    
 
87
 
88
88
    def abs(self, f):
89
89
        """f = |x|
90
90
        d/dx |x| = { +1, if x > 0
92
92
                   {  0, if x == 0 (not strictly correct, but better than leaving it undefined?)
93
93
        """
94
94
        x, = f.operands()
95
 
        dx = sign(x) 
 
95
        dx = sign(x)
96
96
        return (dx,)
97
 
    
 
97
 
98
98
    # --- Mathfunctions
99
 
    
 
99
 
100
100
    def sqrt(self, f):
101
101
        "d/dx sqrt(x) = 1 / (2*sqrt(x))"
102
102
        return (0.5/f,)
103
 
    
 
103
 
104
104
    def exp(self, f):
105
105
        "d/dx exp(x) = exp(x)"
106
106
        return (f,)
107
 
    
 
107
 
108
108
    def ln(self, f):
109
109
        "d/dx ln x = 1 / x"
110
110
        x, = f.operands()
111
111
        return (1/x,)
112
 
    
 
112
 
113
113
    def cos(self, f):
114
114
        "d/dx cos x = -sin(x)"
115
115
        x, = f.operands()
116
116
        return (-sin(x),)
117
 
    
 
117
 
118
118
    def sin(self, f):
119
119
        "d/dx sin x = cos(x)"
120
120
        x, = f.operands()
129
129
        "d/dx acos x = -1/sqrt(1 - x^2)"
130
130
        x, = f.operands()
131
131
        return (-1.0/sqrt(1.0 - x**2),)
132
 
    
 
132
 
133
133
    def asin(self, f):
134
134
        "d/dx asin x = 1/sqrt(1 - x^2)"
135
135
        x, = f.operands()
147
147
 
148
148
    def bessel_function(self, nu, x):
149
149
        return NotImplemented
150
 
    
 
150
 
151
151
    # --- Shape and indexing manipulators
152
152
 
153
153
    def indexed(self, f): # TODO: Is this right? Fix for non-scalars too.
156
156
        ufl_assert(s == (), "TODO: Assuming a scalar expression.")
157
157
        _1 = IntValue(1) # TODO: Non-scalars
158
158
        return (_1, None)
159
 
    
 
159
 
160
160
    def list_tensor(self, f): # TODO: Is this right? Fix for higher order tensors too.
161
161
        "d/dx_i [x_0, ..., x_n-1] = e_i (unit vector)"
162
162
        ops = f.operands()
164
164
        s = ops[0].shape()
165
165
        ufl_assert(s == (), "TODO: Assuming a vector, i.e. scalar operands.")
166
166
        return unit_vectors(n) # TODO: Non-scalars
167
 
    
 
167
 
168
168
    def component_tensor(self, f):
169
169
        x, i = f.operands()
170
170
        s = f.shape()
172
172
        n, = s
173
173
        d = ListTensor([1]*n) # TODO: Non-scalars
174
174
        return (d, None)
175
 
    
 
175
 
176
176
    # --- Restrictions
177
 
    
 
177
 
178
178
    def positive_restricted(self, f):
179
179
        _1 = IntValue(1)
180
 
        return (_1,) # or _1('+')? TODO: is this right? 
 
180
        return (_1,) # or _1('+')? TODO: is this right?
181
181
        # Note that _1('+') would become 0 with the current implementation
182
 
    
 
182
 
183
183
    def negative_restricted(self, f):
184
184
        _1 = IntValue(1)
185
185
        return (_1,) # or _1('-')? TODO: is this right?
186
 
    
 
186
 
187
187
    # --- Conditionals
188
 
    
 
188
 
189
189
    def condition(self, f):
190
190
        return (None, None)
191
 
    
 
191
 
192
192
    def conditional(self, f): # TODO: Is this right? What about non-scalars?
193
193
        c, a, b = f.operands()
194
194
        s = f.shape()
198
198
        da = conditional(c, _1, _0)
199
199
        db = conditional(c, _0, _1)
200
200
        return (None, da, db)
201
 
    
 
201
 
202
202
    # --- Derivatives
203
 
    
 
203
 
204
204
    def spatial_derivative(self, f):
205
205
        error("Partial derivative of spatial_derivative not implemented, "\
206
206
              "when is this called? apply_ad should make sure it isn't called.")
207
207
        x, i = f.operands()
208
208
        return (None, None)
209
 
    
 
209
 
210
210
    def variable_derivative(self, f):
211
211
        error("Partial derivative of variable_derivative not implemented, "\
212
212
              "when is this called? apply_ad should make sure it isn't called.")
213
213
        x, v = f.operands()
214
214
        return (None, None)
215
 
    
 
215
 
216
216
    def coefficient_derivative(self, f):
217
217
        error("Partial derivative of coefficient_derivative not implemented, "\
218
218
              "when is this called? apply_ad should make sure it isn't called.")