src/Plugins/DevelopedSurface/scripts/pp.py

   1 # Copyright (C) 2017-2022  CEA/DEN, EDF R&D
   2 #
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License, or (at your option) any later version.
   7 #
   8 # This library is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19 # Author : Anthony Geay (EDF R&D)
  20
  21 from MEDLoader import *
  22 from math import sqrt
  23 import numpy as np
  24 import scipy
  25 import scipy.sparse.linalg
  26
  27 def f0(sample,p,v,r):
  28     d=sample-p
  29     return d.magnitude()[0]-r
  30
  31 def f(sample,p,v,r):
  32     d=sample-p
  33     l=d-DataArrayDouble.Dot(d,v)[0]*v
  34     return l.magnitude()[0]-r
  35
  36 def f2(sample,zev,ff):
  37     p=zev[0,:3] ; v=zev[0,3:6] ; r=zev[0,6]
  38     return ff(sample,p,v,r)
  39
  40 def df(sample,zev,varid,ff):
  41     eps=0.0001
  42     zev=zev[:]
  43     zev2=zev[:]
  44     zev[0,varid]+=eps ; zev2[0,varid]-=eps
  45     return (f2(sample,zev,ff)-f2(sample,zev2,ff))/(2*eps)
  46
  47 #def df2(sample,p,v,r,varid):
  48 #    zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])])
  49 #    return df(sample,zev,varid)
  50
  51 def jacob(sample,p,v,r,ff):
  52     zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])])
  53     return DataArrayDouble([df(sample,zev,i,ff) for i in range(7)],1,7)
  54
  55 def jacob0(sample,p,v,r):
  56     zev=DataArrayDouble.Meld([p_s,v_s,DataArrayDouble([r_s])])
  57     return DataArrayDouble([df0(sample,zev,i) for i in range(7)],1,7)
  58
  59 mm=MEDFileMesh.New("example2.med")
  60 c=mm.getCoords()
  61 p_s=DataArrayDouble(c.accumulate(),1,3)/float(len(c))
  62 v_s=DataArrayDouble([1,0,0],1,3)
  63 o=DataArrayDouble(c.getMinMaxPerComponent(),3,2)
  64 o0,o1=o.explodeComponents()
  65 o=o1-o0
  66 o.abs()
  67 r_s=o.getMaxValue()[0]/2.
  68 #
  69 r_s=0.215598
  70 p_s=DataArrayDouble([0.,0.,0.],1,3)
  71 v_s=DataArrayDouble([0.,0.,1.],1,3)
  72 #
  73 r_s=0.2
  74 p_s=DataArrayDouble([1.,1.,1.],1,3)
  75 v_s=DataArrayDouble([1.,1.,1.],1,3)
  76 #probes=[0,979,1167,2467,2862,3706,3819]
  77 probes=[1000]+[c[:,i].getMaxValue()[1] for i in range(3)]+[c[:,i].getMinValue()[1] for i in range(3)]
  78
  79 p=p_s ; v=v_s ; r=r_s
  80 for ii in range(1):
  81     mat=DataArrayDouble.Aggregate([jacob(c[probes[0]],p,v,r,f0)]+[jacob(c[probe],p,v,r,f) for probe in probes[1:]])#+[DataArrayDouble([0,0,0,2*v[0,0],2*v[0,1],2*v[0,2],0],1,7)]
  82     y=DataArrayDouble([f0(c[probes[0]],p,v,r)]+[f(c[probe],p,v,r) for probe in probes[1:]])
  83     #y.pushBackSilent(v.magnitude()[0]-1.)
  84     delta2=y[:] ; delta2.abs()
  85     if delta2.getMaxValueInArray()<1e-5:
  86         print("finished")
  87         break;
  88     mat=scipy.matrix(mat.toNumPyArray())
  89     print ii,np.linalg.cond(mat)
  90     y=scipy.matrix(y.toNumPyArray())
  91     y=y.transpose()
  92     delta=np.linalg.solve(mat,-y)
  93     #delta=scipy.sparse.linalg.gmres(mat,-y)[0] # cg
  94     delta=DataArrayDouble(delta) ; delta.rearrange(7)
  95     #p+=delta[0,:3]
  96     #v+=delta[0,3:6]
  97     #v/=v.magnitude()[0]
  98     #r+=delta[0,6]
  99     print ii,delta.getValues()#,yy.transpose().tolist()
 100
 101     pass
 102
 103 mat1=mat[[0,1,2,6]]
 104 mat1=mat1.transpose()
 105 mat1=mat1[[0,1,2,6]]
 106 mat1=mat1.transpose()
 107 mm=np.linalg.solve(mat1,y[[0,1,2,6]])