src/daComposant/daNumerics/Models/TwoDimensionalInverseDistanceCS2010.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2024 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 import sys, unittest, numpy
  24
  25 # ==============================================================================
  26 class TwoDimensionalInverseDistanceCS2010:
  27     """
  28     Two-dimensional inverse distance function of parameter µ=(µ1,µ2):
  29
  30         s(x,y,µ) = 1 / sqrt( (x-µ1)² + (y-µ2)² + 0.1² )
  31
  32         avec (x,y) ∈ [0.1,0.9]² et µ=(µ1,µ2) ∈ [-1,-0.01]².
  33
  34     This is the non-linear parametric function (3.38) of the reference:
  35         Chaturantabut, S., Sorensen, D. C.,
  36         Nonlinear Model Reduction via Discrete Empirical Interpolation,
  37         SIAM Journal on Scientific Computing, 32(5), pp. 2737-2764 (2010).
  38     """
  39     def __init__(self, nx: int = 20, ny: int = 20):
  40         self.nx  = max(1, nx)
  41         self.ny  = max(1, ny)
  42         self.x   = numpy.linspace(0.1, 0.9, self.nx, dtype=float)
  43         self.y   = numpy.linspace(0.1, 0.9, self.ny, dtype=float)
  44
  45     def FunctionH(self, XX ):
  46         __mu1, __mu2 = numpy.ravel(XX)
  47         #
  48         __x, __y = numpy.meshgrid( self.x, self.y )
  49         __sxymu = 1. / numpy.sqrt( (__x - __mu1)**2 + (__y - __mu2)**2 + 0.1**2 )
  50         #
  51         return __sxymu
  52
  53     def get_x(self):
  54         return self.x, self.y
  55
  56     def get_sample_of_mu(self, ns1: int = 20, ns2: int = 20):
  57         smu1 = numpy.linspace(-1, -0.01, ns1, dtype=float)
  58         smu2 = numpy.linspace(-1, -0.01, ns2, dtype=float)
  59         smu = numpy.array([(mu1, mu2) for mu1 in smu1 for mu2 in smu2])
  60         return smu
  61
  62     def get_random_sample_of_mu(self, ns1: int = 1, ns2: int = 1):
  63         smu = []
  64         for i in range(ns1 * ns2):
  65             smu1 = numpy.random.uniform(-1, -0.01)
  66             smu2 = numpy.random.uniform(-1, -0.01)
  67             smu.append((smu1, smu2))
  68         smu = numpy.array(smu)
  69         return smu
  70
  71     def get_bounds_on_space(self):
  72         return [[min(self.x), max(self.x)], [min(self.y), max(self.y)]]
  73
  74     def get_bounds_on_parameter(self):
  75         return [[-1, -0.01]] * 2
  76
  77     OneRealisation = FunctionH
  78
  79 # ==============================================================================
  80 class LocalTest(unittest.TestCase):
  81     @classmethod
  82     def setUpClass(cls):
  83         print('\nAUTODIAGNOSTIC\n==============\n')
  84         print("    " + TwoDimensionalInverseDistanceCS2010().__doc__.strip())
  85
  86     def test001(self):
  87         numpy.random.seed(123456789)
  88         Equation = TwoDimensionalInverseDistanceCS2010()
  89         for mu in Equation.get_sample_of_mu(5, 5):
  90             solution = Equation.OneRealisation( mu )
  91             # Nappe maximale au coin (0,0)
  92             self.assertTrue(numpy.max(solution.flat) <= solution[0, 0])
  93             # Nappe minimale au coin [-1,-1]
  94             self.assertTrue(numpy.min(solution.flat) >= solution[-1, -1])
  95
  96     def tearDown(cls):
  97         print("\n    Tests OK\n")
  98
  99 # ==============================================================================
 100 if __name__ == "__main__":
 101     sys.stderr = sys.stdout
 102     unittest.main(verbosity=0)