1 # Copyright (C) 2017-2020 CEA/DEN, EDF R&D
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.
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.
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
17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 # Author : Anthony Geay (EDF R&D)
21 from MEDLoader import *
25 import scipy.sparse.linalg
29 return d.magnitude()[0]-r
33 l=d-DataArrayDouble.Dot(d,v)[0]*v
34 return l.magnitude()[0]-r
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)
40 def df(sample,zev,varid,ff):
44 zev[0,varid]+=eps ; zev2[0,varid]-=eps
45 return (f2(sample,zev,ff)-f2(sample,zev2,ff))/(2*eps)
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)
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)
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)
59 mm=MEDFileMesh.New("example2.med")
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()
67 r_s=o.getMaxValue()[0]/2.
70 p_s=DataArrayDouble([0.,0.,0.],1,3)
71 v_s=DataArrayDouble([0.,0.,1.],1,3)
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)]
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:
88 mat=scipy.matrix(mat.toNumPyArray())
89 print ii,np.linalg.cond(mat)
90 y=scipy.matrix(y.toNumPyArray())
92 delta=np.linalg.solve(mat,-y)
93 #delta=scipy.sparse.linalg.gmres(mat,-y)[0] # cg
94 delta=DataArrayDouble(delta) ; delta.rearrange(7)
99 print ii,delta.getValues()#,yy.transpose().tolist()
104 mat1=mat1.transpose()
106 mat1=mat1.transpose()
107 mm=np.linalg.solve(mat1,y[[0,1,2,6]])