/[escript]/trunk/downunder/py_src/regularizations.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4387 - (hide annotations)
Thu May 2 00:46:49 2013 UTC (6 years, 6 months ago) by jfenwick
File MIME type: text/x-python
File size: 21675 byte(s)
The future is now
1 jfenwick 4387 from __future__ import print_function
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
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 caltinay 3946
16 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 caltinay 3946 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 caltinay 3947 __all__ = ['Regularization']
24    
25 jfenwick 4387
26 jfenwick 4238 from .costfunctions import CostFunction
27 gross 4074
28 caltinay 3946 import numpy as np
29 caltinay 4211 from esys.escript import Function, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
30 gross 4074 from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
31     from esys.escript.pdetools import ArithmeticTuple
32 gross 4373 from .coordinates import makeTranformation
33 caltinay 3946
34 jfenwick 4238
35 gross 4074 class Regularization(CostFunction):
36 caltinay 3946 """
37 caltinay 4097 The regularization term for the level set function ``m`` within the cost
38     function J for an inversion:
39    
40 caltinay 4285 *J(m)=1/2 * sum_k integrate( mu[k] * ( w0[k] * m_k**2 * w1[k,i] * m_{k,i}**2) + sum_l<k mu_c[l,k] wc[l,k] * | curl(m_k) x curl(m_l) |^2*
41 caltinay 4097
42 gross 4099 where w0[k], w1[k,i] and wc[k,l] are non-negative weighting factors and
43     mu[k] and mu_c[l,k] are trade-off factors which may be altered
44     during the inversion. The weighting factors are normalized such that their
45 caltinay 4097 integrals over the domain are constant:
46    
47 gross 4099 *integrate(w0[k] + inner(w1[k,:],1/L[:]**2))=scale[k]* volume(domain)*
48     *integrate(wc[l,k]*1/L**4)=scale_c[k]* volume(domain) *
49 caltinay 4097
50 caltinay 3946 """
51 caltinay 4097 def __init__(self, domain, numLevelSets=1,
52 gross 4099 w0=None, w1=None, wc=None,
53 caltinay 4097 location_of_set_m=Data(),
54 gross 4099 useDiagonalHessianApproximation=False, tol=1e-8,
55 gross 4373 coordinates=None,
56 gross 4099 scale=None, scale_c=None):
57 caltinay 4007 """
58 caltinay 4095 initialization.
59 caltinay 4097
60     :param domain: domain
61 gross 4074 :type domain: `Domain`
62     :param numLevelSets: number of level sets
63     :type numLevelSets: ``int``
64 gross 4099 :param w0: weighting factor for the m**2 term. If not set zero is assumed.
65     :type w0: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
66 caltinay 4097 (``numLevelSets`` ,) if ``numLevelSets`` > 1
67 gross 4099 :param w1: weighting factor for the grad(m_i) terms. If not set zero is assumed
68     :type w1: ``Vector`` if ``numLevelSets`` == 1 or `Data` object of shape
69 caltinay 4097 (``numLevelSets`` , DIM) if ``numLevelSets`` > 1
70 gross 4099 :param wc: weighting factor for the cross gradient terms. If not set
71 caltinay 4097 zero is assumed. Used for the case if ``numLevelSets`` > 1
72 gross 4099 only. Only values ``wc[l,k]`` in the lower triangle (l<k)
73 caltinay 4097 are used.
74 gross 4099 :type wc: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
75 caltinay 4097 :param location_of_set_m: marks location of zero values of the level
76     set function ``m`` by a positive entry.
77     :type location_of_set_m: ``Scalar`` if ``numLevelSets`` == 1 or `Data`
78     object of shape (``numLevelSets`` ,) if ``numLevelSets`` > 1
79 gross 4099 :param useDiagonalHessianApproximation: if True cross gradient terms
80 caltinay 4097 between level set components are ignored when calculating
81     approximations of the inverse of the Hessian Operator.
82     This can speed-up the calculation of the inverse but may
83     lead to an increase of the number of iteration steps in the
84     inversion.
85 gross 4074 :type useDiagonalHessianApproximation: ``bool``
86 caltinay 4097 :param tol: tolerance when solving the PDE for the inverse of the
87     Hessian Operator
88 gross 4074 :type tol: positive ``float``
89 caltinay 4213
90 gross 4373 :param coordinates: defines coordinate system to be used
91     :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
92 caltinay 4213 :param scale: weighting factor for level set function variation terms. If not set one is used.
93 gross 4099 :type scale: ``Scalar`` if ``numLevelSets`` == 1 or `Data` object of shape
94     (``numLevelSets`` ,) if ``numLevelSets`` > 1
95     :param scale_c: scale for the cross gradient terms. If not set
96     one is assumed. Used for the case if ``numLevelSets`` > 1
97     only. Only values ``scale_c[l,k]`` in the lower triangle (l<k)
98     are used.
99     :type scale_c: `Data` object of shape (``numLevelSets`` , ``numLevelSets``)
100 caltinay 4095
101 caltinay 4007 """
102 gross 4099 if w0 == None and w1==None:
103     raise ValueError("Values for w0 or for w1 must be given.")
104 caltinay 4132 if wc == None and numLevelSets>1:
105     raise ValueError("Values for wc must be given.")
106 caltinay 4097
107 caltinay 3946 self.__domain=domain
108 gross 4074 DIM=self.__domain.getDim()
109     self.__numLevelSets=numLevelSets
110 gross 4373 self.__trafo=makeTranformation(domain, coordinates)
111 caltinay 4213
112 caltinay 4095 self.__pde=LinearPDE(self.__domain, numEquations=self.__numLevelSets)
113 gross 4074 self.__pde.getSolverOptions().setTolerance(tol)
114     self.__pde.setSymmetryOn()
115     try:
116 caltinay 4097 self.__pde.setValue(q=location_of_set_m)
117     except IllegalCoefficientValue:
118     raise ValueError("Unable to set location of fixed level set function.")
119 caltinay 4213
120     # =========== check the shape of the scales: ========================
121 gross 4099 if scale is None:
122 caltinay 4132 if numLevelSets == 1 :
123     scale = 1.
124 gross 4074 else:
125 caltinay 4132 scale = np.ones((numLevelSets,))
126     else:
127     scale=np.asarray(scale)
128     if numLevelSets == 1 :
129     if scale.shape == ():
130     if not scale > 0 :
131     raise ValueError("Value for scale must be positive.")
132     else:
133     raise ValueError("Unexpected shape %s for scale."%scale.shape)
134     else:
135     if scale.shape is (numLevelSets,):
136     if not min(scale) > 0:
137     raise ValueError("All value for scale must be positive.")
138     else:
139     raise ValueError("Unexpected shape %s for scale."%scale.shape)
140 caltinay 4213
141 gross 4099 if scale_c is None or numLevelSets < 2:
142 caltinay 4132 scale_c = np.ones((numLevelSets,numLevelSets))
143     else:
144     scale_c=np.asarray(scale_c)
145     if scale_c.shape == (numLevelSets,numLevelSets):
146 caltinay 4263 if not all( [ [ scale_c[l,k] > 0. for l in range(k) ] for k in range(1,numLevelSets) ]):
147 caltinay 4132 raise ValueError("All values in the lower triangle of scale_c must be positive.")
148 gross 4099 else:
149 caltinay 4132 raise ValueError("Unexpected shape %s for scale."%scale_c.shape)
150 caltinay 4213 # ===== check the shape of the weights: =============================
151 gross 4099 if w0 is not None:
152 caltinay 4132 w0 = interpolate(w0,self.__pde.getFunctionSpaceForCoefficient('D'))
153     s0=w0.getShape()
154     if numLevelSets == 1 :
155     if not s0 == () :
156     raise ValueError("Unexpected shape %s for weight w0."%s0)
157 gross 4099 else:
158 caltinay 4132 if not s0 == (numLevelSets,):
159 caltinay 4213 raise ValueError("Unexpected shape %s for weight w0."%s0)
160 gross 4373 if not self.__trafo.isCartesian():
161     w0*=self.__trafo.getVolumeFactor()
162 caltinay 4132 if not w1 is None:
163     w1 = interpolate(w1,self.__pde.getFunctionSpaceForCoefficient('A'))
164     s1=w1.getShape()
165 gross 4373 if numLevelSets == 1 :
166 caltinay 4132 if not s1 == (DIM,) :
167     raise ValueError("Unexpected shape %s for weight w1."%s1)
168 gross 4099 else:
169 caltinay 4132 if not s1 == (numLevelSets,DIM):
170 caltinay 4213 raise ValueError("Unexpected shape %s for weight w1."%s1)
171 gross 4373 if not self.__trafo.isCartesian():
172     f= self.__trafo.getScalingFactors()**2*self.__trafo.getVolumeFactor()
173     if numLevelSets == 1 :
174     w1*=f
175     else:
176     for i in range(numLevelSets): w1[i,:]*=f
177    
178 caltinay 4213 if numLevelSets == 1:
179 gross 4099 wc=None
180     else:
181     wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
182     sc=wc.getShape()
183 gross 4124 if not sc == (numLevelSets, numLevelSets):
184     raise ValueError("Unexpected shape %s for weight wc."%(sc,))
185 gross 4373 if not self.__trafo.isCartesian():
186     raise ValueError("Non-cartesian coordinates for cross gradient term is not supported yet.")
187 caltinay 4213 # ============= now we rescale weights: =============================
188 gross 4099 L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
189     L4=1/np.sum(1/L2s)**2
190 caltinay 4213 if numLevelSets == 1:
191 gross 4099 A=0
192     if w0 is not None:
193 caltinay 4132 A = integrate(w0)
194     if w1 is not None:
195     A += integrate(inner(w1, 1/L2s))
196     if A > 0:
197     f = scale/A
198     if w0 is not None:
199     w0*=f
200     if w1 is not None:
201     w1*=f
202 gross 4099 else:
203 caltinay 4213 raise ValueError("Non-positive weighting factor detected.")
204 gross 4099 else:
205 gross 4124
206 caltinay 4263 for k in range(numLevelSets):
207 caltinay 4132 A=0
208 gross 4099 if w0 is not None:
209 caltinay 4132 A = integrate(w0[k])
210     if w1 is not None:
211     A += integrate(inner(w1[k,:], 1/L2s))
212     if A > 0:
213     f = scale[k]/A
214     if w0 is not None:
215     w0[k]*=f
216     if w1 is not None:
217     w1[k,:]*=f
218 caltinay 4095 else:
219 caltinay 4213 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
220    
221 caltinay 4132 # and now the cross-gradient:
222     if wc is not None:
223 caltinay 4263 for l in range(k):
224 caltinay 4132 A = integrate(wc[l,k])/L4
225     if A > 0:
226     f = scale_c[l,k]/A
227     wc[l,k]*=f
228 gross 4124 # else:
229 caltinay 4213 # raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
230    
231 gross 4099 self.__w0=w0
232     self.__w1=w1
233     self.__wc=wc
234 caltinay 4097
235 caltinay 4213 self.__pde_is_set=False
236 gross 4099 if self.__numLevelSets > 1:
237     self.__useDiagonalHessianApproximation=useDiagonalHessianApproximation
238     else:
239     self.__useDiagonalHessianApproximation=True
240     self._update_Hessian=True
241 caltinay 4097
242 gross 4099 self.__num_tradeoff_factors=numLevelSets+((numLevelSets-1)*numLevelSets)/2
243     self.setTradeOffFactors()
244 gross 4102 self.__vol_d=vol(self.__domain)
245 caltinay 4213
246 gross 4074 def getDomain(self):
247 caltinay 3946 """
248 caltinay 4097 returns the domain of the regularization term
249    
250 gross 4074 :rtype: ``Domain``
251 caltinay 3946 """
252 gross 4074 return self.__domain
253 gross 4373
254     def getCoordinateTransformation(self):
255     """
256     returns the coordinate transformation being used
257 caltinay 4097
258 gross 4373 :rtype: ``CoordinateTransformation``
259     """
260     return self.__trafo
261    
262 gross 4074 def getNumLevelSets(self):
263     """
264 caltinay 4097 returns the number of level set functions
265    
266 gross 4074 :rtype: ``int``
267     """
268     return self.__numLevelSets
269 caltinay 3946
270 gross 4074 def getPDE(self):
271 caltinay 4007 """
272 gross 4074 returns the linear PDE to be solved for the Hessian Operator inverse
273 caltinay 4097
274     :rtype: `LinearPDE`
275 caltinay 4007 """
276 gross 4074 return self.__pde
277 caltinay 4213
278 gross 4074 def getDualProduct(self, m, r):
279     """
280 caltinay 4097 returns the dual product of a gradient represented by X=r[1] and Y=r[0]
281     with a level set function m:
282    
283 gross 4074 *Y_i*m_i + X_ij*m_{i,j}*
284 caltinay 4097
285     :type m: `Data`
286     :type r: `ArithmeticTuple`
287     :rtype: ``float``
288 gross 4074 """
289     A=0
290     if not r[0].isEmpty(): A+=integrate(inner(r[0], m))
291     if not r[1].isEmpty(): A+=integrate(inner(r[1], grad(m)))
292 caltinay 4095 return A
293 gross 4099 def getNumTradeOffFactors(self):
294     """
295     returns the number of trade-off factors being used.
296 caltinay 3946
297 gross 4099 :rtype: ``int``
298     """
299     return self.__num_tradeoff_factors
300    
301     def setTradeOffFactors(self, mu=None):
302     """
303     sets the trade-off factors for the level-set variation and the cross-gradient
304 caltinay 4213
305     :param mu: new values for the trade-off factors where values mu[:numLevelSets] are the
306     trade-off factors for the level-set variation and the remaining values for
307 gross 4099 the cross-gradient part with mu_c[l,k]=mu[numLevelSets+l+((k-1)*k)/2] (l<k).
308     If no values for mu is given ones are used. Values must be positive.
309     :type mu: ``list`` of ``float`` or ```numpy.array```
310     """
311     numLS=self.getNumLevelSets()
312     numTF=self.getNumTradeOffFactors()
313     if mu is None:
314 caltinay 4132 mu = np.ones((numTF,))
315     else:
316     mu = np.asarray(mu)
317 gross 4099
318 caltinay 4132 if mu.shape == (numTF,):
319     self.setTradeOffFactorsForVariation(mu[:numLS])
320     mu_c2=np.zeros((numLS,numLS))
321 caltinay 4263 for k in range(numLS):
322     for l in range(k):
323 caltinay 4213 mu_c2[l,k] = mu[numLS+l+((k-1)*k)/2]
324 caltinay 4132 self.setTradeOffFactorsForCrossGradient(mu_c2)
325     elif mu.shape == () and numLS ==1:
326     self.setTradeOffFactorsForVariation(mu)
327     else:
328 caltinay 4213 raise ValueError("Unexpected shape %s for mu."%(mu.shape,))
329    
330 gross 4099 def setTradeOffFactorsForVariation(self, mu=None):
331     """
332     sets the trade-off factors for the level-set variation part
333 caltinay 4213
334 gross 4099 :param mu: new values for the trade-off factors. Values must be positive.
335 caltinay 4132 :type mu: ``float``, ``list`` of ``float`` or ```numpy.array```
336 gross 4099 """
337     numLS=self.getNumLevelSets()
338     if mu is None:
339 caltinay 4132 if numLS == 1:
340     mu = 1.
341     else:
342     mu = np.ones((numLS,))
343 gross 4099
344 caltinay 4132 mu=np.asarray(mu)
345     if numLS == 1:
346 caltinay 4213 if mu.shape == (1,): mu=mu[0]
347 caltinay 4132 if mu.shape == ():
348     if mu > 0:
349     self.__mu= mu
350     self._new_mu=True
351     else:
352 caltinay 4213 raise ValueError("Value for trade-off factor must be positive.")
353 gross 4099 else:
354 caltinay 4132 raise ValueError("Unexpected shape %s for mu."%mu.shape)
355     else:
356     if mu.shape == (numLS,):
357     if min(mu) > 0:
358     self.__mu= mu
359     self._new_mu=True
360     else:
361     raise ValueError("All value for mu must be positive.")
362     else:
363 caltinay 4213 raise ValueError("Unexpected shape %s for trade-off factor."%mu.shape)
364 gross 4099
365     def setTradeOffFactorsForCrossGradient(self, mu_c=None):
366     """
367 caltinay 4132 sets the trade-off factors for the cross-gradient terms
368 caltinay 4213
369 caltinay 4132 :param mu_c: new values for the trade-off factors for the cross-gradient
370     terms. Values must be positive. If no value is given ones
371     are used. Onky value mu_c[l,k] for l<k are used.
372     :type mu_c: ``float``, ``list`` of ``float`` or ``numpy.array``
373 gross 4099 """
374     numLS=self.getNumLevelSets()
375     if mu_c is None or numLS < 2:
376 caltinay 4132 self.__mu_c = np.ones((numLS,numLS))
377     if isinstance(mu_c, float) or isinstance(mu_c, int):
378     self.__mu_c = np.zeros((numLS,numLS))
379 caltinay 4213 self.__mu_c[:,:]=mu_c
380 caltinay 4132 else:
381     mu_c=np.asarray(mu_c)
382     if mu_c.shape == (numLS,numLS):
383 caltinay 4263 if not all( [ [ mu_c[l,k] > 0. for l in range(k) ] for k in range(1,numLS) ]):
384 caltinay 4132 raise ValueError("All trade-off factors in the lower triangle of mu_c must be positive.")
385     else:
386     self.__mu_c = mu_c
387     self._new_mu=True
388 gross 4099 else:
389 caltinay 4132 raise ValueError("Unexpected shape %s for mu."%mu_c.shape)
390 caltinay 4213
391 gross 4099 def getArguments(self, m):
392     """
393     """
394     return ( grad(m),)
395    
396 gross 4074 def getValue(self, m, grad_m):
397 caltinay 4007 """
398 caltinay 4097 returns the value of the cost function J with respect to m.
399    
400 gross 4074 :rtype: ``float``
401 caltinay 4007 """
402 gross 4099 mu=self.__mu
403     mu_c=self.__mu_c
404 gross 4074 DIM=self.getDomain().getDim()
405     numLS=self.getNumLevelSets()
406 caltinay 4097
407 caltinay 3946 A=0
408 gross 4099 if self.__w0 is not None:
409     A+=inner(integrate(m**2 * self.__w0), mu)
410 caltinay 4097
411 gross 4099 if self.__w1 is not None:
412 caltinay 4095 if numLS == 1:
413 gross 4099 A+=integrate(inner(grad_m**2, self.__w1))*mu
414 caltinay 4095 else:
415 caltinay 4263 for k in range(numLS):
416 jfenwick 4387 print("C = ",integrate(inner(grad_m[k,:]**2,self.__w1[k,:])))
417 gross 4099 A+=mu[k]*integrate(inner(grad_m[k,:]**2,self.__w1[k,:]))
418 caltinay 4213
419 gross 4099 if numLS > 1:
420 caltinay 4263 for k in range(numLS):
421 caltinay 4097 gk=grad_m[k,:]
422     len_gk=length(gk)
423 caltinay 4263 for l in range(k):
424 caltinay 4097 gl=grad_m[l,:]
425 jfenwick 4387 print("CC =",integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 ))
426 gross 4099 A+= mu_c[l,k] * integrate( self.__wc[l,k] * ( len_gk * length(gl) )**2 - inner(gk, gl)**2 )
427 jfenwick 4387 print("A =",A)
428 gross 4074 return A/2
429 caltinay 4097
430 gross 4074 def getGradient(self, m, grad_m):
431     """
432 caltinay 4097 returns the gradient of the cost function J with respect to m.
433     The function returns Y_k=dPsi/dm_k and X_kj=dPsi/dm_kj
434 gross 4074 """
435 caltinay 4097
436 gross 4099 mu=self.__mu
437     mu_c=self.__mu_c
438 gross 4074 DIM=self.getDomain().getDim()
439     numLS=self.getNumLevelSets()
440 caltinay 4097
441 caltinay 4211 grad_m=grad(m, Function(m.getDomain()))
442 gross 4099 if self.__w0 is not None:
443     Y = m * self.__w0 * mu
444 caltinay 4095 else:
445 caltinay 4132 if numLS == 1:
446     Y = Scalar(0, grad_m.getFunctionSpace())
447     else:
448 gross 4122 Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
449 caltinay 3946
450 gross 4099 if self.__w1 is not None:
451 caltinay 4132 X=grad_m*self.__w1
452 caltinay 4095 if numLS == 1:
453 gross 4099 X=grad_m* self.__w1*mu
454 caltinay 4095 else:
455 caltinay 4263 for k in range(numLS):
456 gross 4099 X[k,:]*=mu[k]
457 caltinay 4095 else:
458 gross 4122 X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
459 caltinay 4213
460 gross 4122 # cross gradient terms:
461 caltinay 4213 if numLS > 1:
462 caltinay 4263 for k in range(numLS):
463 caltinay 4132 grad_m_k=grad_m[k,:]
464     l2_grad_m_k = length(grad_m_k)**2
465 caltinay 4263 for l in range(k):
466 caltinay 4132 grad_m_l=grad_m[l,:]
467     l2_grad_m_l = length(grad_m_l)**2
468     grad_m_lk = inner(grad_m_l, grad_m_k)
469     f= mu_c[l,k]* self.__wc[l,k]
470     X[l,:] += f * ( l2_grad_m_l * grad_m_l - grad_m_lk * grad_m_k )
471     X[k,:] += f * ( l2_grad_m_k * grad_m_k - grad_m_lk * grad_m_l )
472 caltinay 4213
473 gross 4074 return ArithmeticTuple(Y, X)
474 caltinay 4097
475 gross 4074 def getInverseHessianApproximation(self, m, r, grad_m):
476     """
477     """
478     if self._new_mu or self._update_Hessian:
479 caltinay 4097
480 caltinay 4095 self._new_mu=False
481     self._update_Hessian=False
482 gross 4099 mu=self.__mu
483     mu_c=self.__mu_c
484 caltinay 4097
485 gross 4074 DIM=self.getDomain().getDim()
486     numLS=self.getNumLevelSets()
487 gross 4099 if self.__w0 is not None:
488 caltinay 4095 if numLS == 1:
489 gross 4099 D=self.__w0 * mu
490 caltinay 4097 else:
491 gross 4074 D=self.getPDE().createCoefficient("D")
492 caltinay 4263 for k in range(numLS): D[k,k]=self.__w0[k] * mu[k]
493 caltinay 4095 self.getPDE().setValue(D=D)
494 caltinay 3946
495 caltinay 4095 A=self.getPDE().createCoefficient("A")
496 gross 4099 if self.__w1 is not None:
497 caltinay 4095 if numLS == 1:
498 caltinay 4263 for i in range(DIM): A[i,i]=self.__w1[i] * mu
499 caltinay 4095 else:
500 caltinay 4263 for k in range(numLS):
501     for i in range(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
502 caltinay 4213
503 gross 4099 if numLS > 1:
504 caltinay 4263 for k in range(numLS):
505 caltinay 4132 grad_m_k=grad_m[k,:]
506     l2_grad_m_k = length(grad_m_k)**2
507     o_kk=outer(grad_m_k, grad_m_k)
508 caltinay 4263 for l in range(k):
509 caltinay 4132 grad_m_l=grad_m[l,:]
510     l2_grad_m_l = length(grad_m_l)**2
511     i_lk = inner(grad_m_l, grad_m_k)
512     o_lk = outer(grad_m_l, grad_m_k)
513     o_kl = outer(grad_m_k, grad_m_l)
514     o_ll=outer(grad_m_l, grad_m_l)
515     f= mu_c[l,k]* self.__wc[l,k]
516 caltinay 4213
517 caltinay 4132 A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
518     A[l,:,k,:] += f * ( 2 * o_lk - o_kl - i_lk * kronecker(DIM) )
519     A[k,:,l,:] += f * ( 2 * o_kl - o_lk - i_lk * kronecker(DIM) )
520     A[k,:,k,:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
521 gross 4074 self.getPDE().setValue(A=A)
522 gross 4102 #self.getPDE().resetRightHandSideCoefficients()
523     #self.getPDE().setValue(X=r[1])
524     #print "X only: ",self.getPDE().getSolution()
525     #self.getPDE().resetRightHandSideCoefficients()
526     #self.getPDE().setValue(Y=r[0])
527     #print "Y only: ",self.getPDE().getSolution()
528 caltinay 4097
529 gross 4079 self.getPDE().resetRightHandSideCoefficients()
530 gross 4074 self.getPDE().setValue(X=r[1], Y=r[0])
531     return self.getPDE().getSolution()
532 caltinay 4097
533 gross 4074 def updateHessian(self):
534     """
535     notify the class to recalculate the Hessian operator
536     """
537 caltinay 4095 if not self.__useDiagonalHessianApproximation:
538     self._update_Hessian=True
539 caltinay 4213
540 gross 4100 def getNorm(self, m):
541     """
542     returns the norm of ``m``
543    
544     :param m: level set function
545     :type m: `Data`
546     :rtype: ``float``
547     """
548 caltinay 4132 return sqrt(integrate(length(m)**2)/self.__vol_d)
549 caltinay 4213

  ViewVC Help
Powered by ViewVC 1.1.26