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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3981 - (show annotations)
Fri Sep 21 02:47:54 2012 UTC (7 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 2557 byte(s)
First pass of updating copyright notices
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 http://www.uq.edu.au
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 __all__ = ['Regularization']
24
25 import numpy as np
26 from esys.escript.linearPDEs import LinearSinglePDE
27 from esys.escript import Data, grad, inner, integrate, kronecker, vol
28
29 class Regularization(object):
30 """
31 """
32 def __init__(self, domain, m_ref=0, w0=None, w=None, location_of_set_m=Data(), tol=1e-8):
33 self.__domain=domain
34 self.__m_ref=m_ref
35 self.location_of_set_m=location_of_set_m
36 if w0 is None:
37 self._w0 = None
38 else:
39 self._w0=np.asarray(w0) / vol(self.__domain)
40 if w is None:
41 self._w = None
42 else:
43 self._w=np.asarray(w)
44
45
46 self.__projector=LinearSinglePDE(domain)
47 self.__projector.getSolverOptions().setTolerance(tol)
48 self.__projector.setValue(A=kronecker(domain), q=location_of_set_m, D=0.)
49
50 def getInner(self, f0, f1):
51 """
52 returns the inner product of two gradients.
53 """
54 return integrate(inner(grad(f0),grad(f1)))
55
56 def project(self, Y=Data(), X=Data()):
57 self.__projector.setValue(Y=Y, X=X)
58 return self.__projector.getSolution()
59
60 def getValue(self, m):
61 A=0
62 if self._w0 is not None:
63 A=(m-self.__m_ref)**2 * self._w0
64 if self._w is not None:
65 A=inner(self._w, grad(m-self.__m_ref)**2) + A
66 return integrate(A)
67
68 def getGradient(self, m):
69 if not self._w0 == None:
70 Y=2. * (m-self.__m_ref) * self._w0
71 else:
72 Y=0.
73 if not self._w == None:
74 X=2. * grad(m-self.__m_ref) * self._w
75 else:
76 X=np.zeros((self.__domain.getDim(),))
77
78 return Y, X
79
80 def getArguments(self, m):
81 return ()
82

  ViewVC Help
Powered by ViewVC 1.1.26