/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4154 by jfenwick, Tue Jan 22 09:30:23 2013 UTC revision 4237 by jfenwick, Thu Feb 21 07:27:36 2013 UTC
# Line 441  class Locator(object): Line 441  class Locator(object):
441        if isinstance(data, escore.Data):        if isinstance(data, escore.Data):
442           if data.getFunctionSpace()!=self.getFunctionSpace():           if data.getFunctionSpace()!=self.getFunctionSpace():
443             raise TypeError("setValue: FunctionSpace of Locator and Data object must match.")             raise TypeError("setValue: FunctionSpace of Locator and Data object must match.")
444           data.expand()             data.expand()  
445           id=self.getId()           id=self.getId()
446           if isinstance(id, list):           if isinstance(id, list):
447            for i in id:            for i in id:
# Line 459  def getInfLocator(arg): Line 459  def getInfLocator(arg):
459      if not isinstance(arg, escore.Data):      if not isinstance(arg, escore.Data):
460         raise TypeError("getInfLocator: Unknown argument type.")         raise TypeError("getInfLocator: Unknown argument type.")
461      a_inf=util.inf(arg)      a_inf=util.inf(arg)
462      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint()     # This gives us the location but not coords
463      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
464      x_min=x.getTupleForGlobalDataPoint(*loc)      x_min=x.getTupleForGlobalDataPoint(*loc)
465      return Locator(arg.getFunctionSpace(),x_min)      return Locator(arg.getFunctionSpace(),x_min)
# Line 471  def getSupLocator(arg): Line 471  def getSupLocator(arg):
471      if not isinstance(arg, escore.Data):      if not isinstance(arg, escore.Data):
472         raise TypeError("getInfLocator: Unknown argument type.")         raise TypeError("getInfLocator: Unknown argument type.")
473      a_inf=util.sup(arg)      a_inf=util.sup(arg)
474      loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords      loc=util.length(arg-a_inf).minGlobalDataPoint()     # This gives us the location but not coords
475      x=arg.getFunctionSpace().getX()      x=arg.getFunctionSpace().getX()
476      x_min=x.getTupleForGlobalDataPoint(*loc)      x_min=x.getTupleForGlobalDataPoint(*loc)
477      return Locator(arg.getFunctionSpace(),x_min)      return Locator(arg.getFunctionSpace(),x_min)
478                
479    
480  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
481     """     """
# Line 746  def NewtonGMRES(defect, x, iter_max=100, Line 746  def NewtonGMRES(defect, x, iter_max=100,
746     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
747              if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)              if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
748              #              #
749          #   adjust subtol_              #   adjust subtol_
750          #              #
751              if iter > 1:              if iter > 1:
752                 rat=fnrm/fnrmo                 rat=fnrm/fnrmo
753                 subtol_old=subtol                 subtol_old=subtol
754                 subtol=gamma*rat**2                 subtol=gamma*rat**2
755                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
756                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
757          #              #
758          # calculate newton increment xc              # calculate newton increment xc
759              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
760              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
761              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
# Line 1394  class ArithmeticTuple(object): Line 1394  class ArithmeticTuple(object):
1394             if l!=len(self):             if l!=len(self):
1395                 raise ValueError("length of arguments don't match.")                 raise ValueError("length of arguments don't match.")
1396             for i in range(l):             for i in range(l):
1397          if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):                  if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1398              out.append(escore.Data())                      out.append(escore.Data())
1399          else:                  else:
1400              out.append(self[i]*other[i])                      out.append(self[i]*other[i])
1401         except TypeError:         except TypeError:
1402          for i in range(len(self)):               for i in range(len(self)):  
1403          if self.__isEmpty(self[i]) or self.__isEmpty(other):                  if self.__isEmpty(self[i]) or self.__isEmpty(other):
1404              out.append(escore.Data())                      out.append(escore.Data())
1405          else:                  else:
1406              out.append(self[i]*other)                      out.append(self[i]*other)
1407         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1408    
1409     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1417  class ArithmeticTuple(object): Line 1417  class ArithmeticTuple(object):
1417        """        """
1418        out=[]        out=[]
1419        try:        try:
1420        l=len(other)            l=len(other)
1421        if l!=len(self):            if l!=len(self):
1422            raise ValueError("length of arguments don't match.")                raise ValueError("length of arguments don't match.")
1423        for i in range(l):            for i in range(l):
1424          if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):                  if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1425              out.append(escore.Data())                      out.append(escore.Data())
1426          else:                  else:
1427              out.append(other[i]*self[i])                      out.append(other[i]*self[i])
1428        except TypeError:        except TypeError:
1429        for i in range(len(self)):              for i in range(len(self)):  
1430          if self.__isEmpty(self[i]) or self.__isEmpty(other):                  if self.__isEmpty(self[i]) or self.__isEmpty(other):
1431              out.append(escore.Data())                      out.append(escore.Data())
1432          else:                  else:
1433              out.append(other*self[i])                      out.append(other*self[i])
1434        return ArithmeticTuple(*tuple(out))        return ArithmeticTuple(*tuple(out))
1435    
1436     def __div__(self,other):     def __div__(self,other):
# Line 1455  class ArithmeticTuple(object): Line 1455  class ArithmeticTuple(object):
1455        """        """
1456        out=[]        out=[]
1457        try:        try:
1458        l=len(other)            l=len(other)
1459        if l!=len(self):            if l!=len(self):
1460            raise ValueError("length of arguments don't match.")                raise ValueError("length of arguments don't match.")
1461                    
1462        for i in range(l):            for i in range(l):
1463          if self.__isEmpty(self[i]):                  if self.__isEmpty(self[i]):
1464              raise ZeroDivisionError("in component %s"%i)                      raise ZeroDivisionError("in component %s"%i)
1465          else:                  else:
1466              if self.__isEmpty(other[i]):                      if self.__isEmpty(other[i]):
1467              out.append(escore.Data())                          out.append(escore.Data())
1468              else:                      else:
1469              out.append(other[i]/self[i])                          out.append(other[i]/self[i])
1470        except TypeError:        except TypeError:
1471        for i in range(len(self)):            for i in range(len(self)):
1472          if self.__isEmpty(self[i]):                  if self.__isEmpty(self[i]):
1473              raise ZeroDivisionError("in component %s"%i)                      raise ZeroDivisionError("in component %s"%i)
1474          else:                  else:
1475              if self.__isEmpty(other):                      if self.__isEmpty(other):
1476              out.append(escore.Data())                          out.append(escore.Data())
1477              else:                      else:
1478              out.append(other/self[i])                          out.append(other/self[i])
1479        return ArithmeticTuple(*tuple(out))        return ArithmeticTuple(*tuple(out))
1480    
1481     def __iadd__(self,other):     def __iadd__(self,other):
# Line 1486  class ArithmeticTuple(object): Line 1486  class ArithmeticTuple(object):
1486        :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1487        """        """
1488        if len(self) != len(other):        if len(self) != len(other):
1489        raise ValueError("tuple lengths must match.")            raise ValueError("tuple lengths must match.")
1490        for i in range(len(self)):        for i in range(len(self)):
1491        if self.__isEmpty(self.__items[i]):            if self.__isEmpty(self.__items[i]):
1492            self.__items[i]=other[i]                self.__items[i]=other[i]
1493        else:            else:
1494            self.__items[i]+=other[i]                self.__items[i]+=other[i]
1495                            
1496        return self        return self
1497    
1498     def __add__(self,other):     def __add__(self,other):
# Line 1504  class ArithmeticTuple(object): Line 1504  class ArithmeticTuple(object):
1504        """        """
1505        out=[]        out=[]
1506        try:        try:
1507        l=len(other)            l=len(other)
1508        if l!=len(self):            if l!=len(self):
1509            raise ValueError("length of arguments don't match.")                raise ValueError("length of arguments don't match.")
1510        for i in range(l):            for i in range(l):
1511          if self.__isEmpty(self[i]):                  if self.__isEmpty(self[i]):
1512              out.append(other[i])                      out.append(other[i])
1513          elif self.__isEmpty(other[i]):                  elif self.__isEmpty(other[i]):
1514              out.append(self[i])                      out.append(self[i])
1515          else:                  else:
1516              out.append(self[i]+other[i])                      out.append(self[i]+other[i])
1517        except TypeError:        except TypeError:
1518          for i in range(len(self)):                  for i in range(len(self)):    
1519          if self.__isEmpty(self[i]):                  if self.__isEmpty(self[i]):
1520              out.append(other)                      out.append(other)
1521          elif self.__isEmpty(other):                  elif self.__isEmpty(other):
1522              out.append(self[i])                      out.append(self[i])
1523          else:                  else:
1524              out.append(self[i]+other)                      out.append(self[i]+other)
1525        return ArithmeticTuple(*tuple(out))        return ArithmeticTuple(*tuple(out))
1526    
1527     def __sub__(self,other):     def __sub__(self,other):
# Line 1533  class ArithmeticTuple(object): Line 1533  class ArithmeticTuple(object):
1533        """        """
1534        out=[]        out=[]
1535        try:        try:
1536        l=len(other)            l=len(other)
1537        if l!=len(self):            if l!=len(self):
1538            raise ValueError("length of arguments don't match.")                raise ValueError("length of arguments don't match.")
1539        for i in range(l):            for i in range(l):
1540          if self.__isEmpty(other[i]):                  if self.__isEmpty(other[i]):
1541              out.append(self[i])                      out.append(self[i])
1542          elif self.__isEmpty(self[i]):                  elif self.__isEmpty(self[i]):
1543              out.append(-other[i])                      out.append(-other[i])
1544          else:                  else:
1545              out.append(self[i]-other[i])                      out.append(self[i]-other[i])
1546        except TypeError:        except TypeError:
1547          for i in range(len(self)):                  for i in range(len(self)):    
1548          if  self.__isEmpty(other):                  if  self.__isEmpty(other):
1549              out.append(self[i])                      out.append(self[i])
1550          elif self.__isEmpty(self[i]):                  elif self.__isEmpty(self[i]):
1551              out.append(-other)                      out.append(-other)
1552          else:                  else:
1553              out.append(self[i]-other)                      out.append(self[i]-other)
1554                                    
1555        return ArithmeticTuple(*tuple(out))        return ArithmeticTuple(*tuple(out))
1556    
1557     def __isub__(self,other):     def __isub__(self,other):
# Line 1562  class ArithmeticTuple(object): Line 1562  class ArithmeticTuple(object):
1562        :type other: ``ArithmeticTuple``        :type other: ``ArithmeticTuple``
1563        """        """
1564        if len(self) != len(other):        if len(self) != len(other):
1565        raise ValueError("tuple length must match.")            raise ValueError("tuple length must match.")
1566        for i in range(len(self)):        for i in range(len(self)):
1567        if not self.__isEmpty(other[i]):            if not self.__isEmpty(other[i]):
1568            if self.__isEmpty(self.__items[i]):                if self.__isEmpty(self.__items[i]):
1569            self.__items[i]=-other[i]                    self.__items[i]=-other[i]
1570            else:                else:
1571            self.__items[i]=other[i]                    self.__items[i]=other[i]
1572        return self        return self
1573    
1574     def __neg__(self):     def __neg__(self):
# Line 1577  class ArithmeticTuple(object): Line 1577  class ArithmeticTuple(object):
1577        """        """
1578        out=[]        out=[]
1579        for i in range(len(self)):        for i in range(len(self)):
1580        if self.__isEmpty(self[i]):            if self.__isEmpty(self[i]):
1581            out.append(escore.Data())                out.append(escore.Data())
1582        else:            else:
1583            out.append(-self[i])                out.append(-self[i])
1584                    
1585        return ArithmeticTuple(*tuple(out))        return ArithmeticTuple(*tuple(out))
1586     def __isEmpty(self, d):     def __isEmpty(self, d):
1587      if isinstance(d, escore.Data):      if isinstance(d, escore.Data):
1588      return d.isEmpty()          return d.isEmpty()
1589      else:      else:
1590      return False          return False
1591    
1592    
1593  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
# Line 1802  class HomogeneousSaddlePointProblem(obje Line 1802  class HomogeneousSaddlePointProblem(obje
1802            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1803            if f<0: raise ValueError("negative pressure norm.")            if f<0: raise ValueError("negative pressure norm.")
1804            return math.sqrt(f)            return math.sqrt(f)
1805                    
1806        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1807           """           """
1808           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.

Legend:
Removed from v.4154  
changed lines
  Added in v.4237

  ViewVC Help
Powered by ViewVC 1.1.26