src/TorseurCIH/Test.py

   1 # Copyright (C) 2021  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
  20 from medcoupling import *
  21 import math
  22
  23 def ReturnInertia(p,OM,area_field):
  24     base_X = DataArrayDouble(len(OM),3) ; base_X[:]=p
  25     dist_to_base_X = (OM-DataArrayDouble.Dot(OM,base_X)*base_X).magnitude()
  26     inertia = (dist_to_base_X*dist_to_base_X*area_field).accumulate()[0]
  27     return inertia
  28
  29 def fffff(initialVect,normalFace,OM,area_field, posToIterate):
  30     li=[]
  31     for zePos in posToIterate:
  32         p = initialVect.deepCopy()
  33         MEDCouplingPointSet.Rotate3DAlg([0,0,0],normalFace.getValues(),zePos/float(180)*math.pi,p)
  34         inertia = ReturnInertia(p,OM,area_field)
  35         li.append((zePos,p.deepCopy(),inertia))
  36     return max(li,key=lambda x:x[2])
  37
  38 def fff(initialVect,normalFace,OM,area_field):
  39     pos = fffff(initialVect,normalFace,OM,area_field,[i*float(10) for i in range(18)])[0]
  40     for expo in range(5):
  41         pos,p,v = fffff(initialVect,normalFace,OM,area_field,[pos+i*(10**-expo) for i in range(-9,10)])
  42     return pos,p,v
  43
  44 """
  45 EDF24091 : probleme de calcul dans TorseurCIH
  46 """
  47
  48 fname = "Case_22_09_21.med"
  49 fieldName = "LIN_____SIEF_NOEU"
  50 mm=MEDFileMesh.New(fname)
  51 m=mm[0]
  52 f1ts = MEDFileField1TS(fname,fieldName)
  53 f = f1ts.field(mm)
  54 m = f.getMesh()
  55 area_field = m.getMeasureField(True)
  56 area = area_field.accumulate()[0] # 1
  57 centerOfMassField = m.computeCellCenterOfMass()
  58 centerOfMass = DataArrayDouble([elt/area for elt in (centerOfMassField*area_field.getArray()).accumulate()],1,3) # 2
  59 m.unPolyze()
  60 tri = MEDCoupling1SGTUMesh(m)
  61 assert(tri.getCellModelEnum()==NORM_TRI3)
  62 #
  63 fCell = f.nodeToCellDiscretization()
  64 (fCell.getArray()[:,[0,1,2]]*area).accumulate()
  65 ids = area_field.getArray().findIdsLowerThan(1e-7).buildComplement(m.getNumberOfCells())
  66 #fCell[ids]
  67 eqn = m[ids].computePlaneEquationOf3DFaces()[:,:3]
  68 eqn /= eqn.magnitude()
  69 area_vector = eqn*area_field.getArray()[ids]
  70 matrix = fCell[ids].getArray()
  71 #
  72 F_x = matrix[:,0]*eqn[:,0] + matrix[:,3]*eqn[:,1] + matrix[:,4]*eqn[:,2]
  73 F_y = matrix[:,3]*eqn[:,0] + matrix[:,1]*eqn[:,1] + matrix[:,5]*eqn[:,2]
  74 F_z = matrix[:,4]*eqn[:,0] + matrix[:,5]*eqn[:,1] + matrix[:,2]*eqn[:,2]
  75 #
  76 F = DataArrayDouble.Meld([F_x,F_y,F_z])
  77 F[:] *= area_vector
  78 #
  79 ZeForce = DataArrayDouble(F.accumulate(),1,3)
  80 normalFace = DataArrayDouble(eqn.accumulate(),1,3)
  81 normalFace /= normalFace.magnitude()[0]
  82 ForceNormale = DataArrayDouble.Dot(ZeForce,normalFace)[0]*normalFace # 3
  83 TangentForce = ZeForce-ForceNormale # 4
  84 #
  85 bary = m[ids].computeCellCenterOfMass()
  86 OM = bary-centerOfMass
  87 momentum = DataArrayDouble(DataArrayDouble.CrossProduct(OM,F).accumulate(),1,3) # 5
  88 # Inertie
  89 InertiaNormale = (DataArrayDouble.Dot(OM,OM)*area_field.getArray()[ids]).accumulate()[0] # 6_A
  90 base = DataArrayDouble(DataArrayDouble.GiveBaseForPlane(normalFace),3,3)
  91 angle, tangentPrincipal, inertiaPrincipal = fff(base[0],normalFace,OM,area_field.getArray()[ids])
  92 tangentOther = DataArrayDouble.CrossProduct(normalFace,tangentPrincipal)
  93 inertiaOther = ReturnInertia(tangentOther,OM,area_field.getArray()[ids])
  94 """
  95 base_X = DataArrayDouble(len(ids),3) ; base_X[:]=base[0]
  96 base_Y = DataArrayDouble(len(ids),3) ; base_Y[:]=base[1]
  97 dist_to_base_X = (OM-DataArrayDouble.Dot(OM,base_X)*base_X).magnitude()
  98 dist_to_base_Y = (OM-DataArrayDouble.Dot(OM,base_Y)*base_Y).magnitude()
  99 inertia_mat_0 = (dist_to_base_Y*dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
 100 inertia_mat_1 = (dist_to_base_X*dist_to_base_X*area_field.getArray()[ids]).accumulate()[0]
 101 inertia_mat_01 = -(dist_to_base_X*dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
 102 from numpy import linalg as LA
 103 import numpy as np
 104 mat = np.matrix([[inertia_mat_0, inertia_mat_01], [inertia_mat_01, inertia_mat_1]])
 105 v,w = LA.eig(mat)
 106 pos_of_max = max([(i,elt) for (i,elt) in enumerate(v)],key=lambda x: x[1])[0]
 107 u0 = DataArrayDouble(np.array(w[:,pos_of_max])) ; u0.rearrange(2)
 108 v0 = DataArrayDouble(np.array(w[:,1-pos_of_max])) ; v0.rearrange(2)
 109 #
 110 I_new_base_0 = v[pos_of_max] # 6_B
 111 new_base_0 = u0[0,0]*base[0]+u0[0,1]*base[1] # 6_B
 112 #new_base_1 = v0[0,0]*base[0]+v0[0,1]*base[1]
 113 new_base_1 = DataArrayDouble.CrossProduct(normalFace,new_base_0)
 114 new_base_Y = DataArrayDouble(len(ids),3) ; new_base_Y[:]=new_base_1
 115 new_dist_to_base_Y = (OM-DataArrayDouble.Dot(OM,new_base_Y)*new_base_Y).magnitude()
 116 I_new_base_1 = (new_dist_to_base_Y*new_dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
 117 """
 118 """
 119 new_base_X = DataArrayDouble(len(ids),3) ; new_base_X[:]=new_base_0
 120 new_dist_to_base_X = (OM-DataArrayDouble.Dot(OM,new_base_X)*new_base_X).magnitude()
 121 I_new_base_0 = (new_dist_to_base_X*new_dist_to_base_X*area_field.getArray()[ids]).accumulate()[0]
 122 tmp=m[ids] ; tmp.zipCoords()
 123 f=MEDCouplingFieldDouble(ON_CELLS)
 124 f.setMesh(tmp)
 125 f.setArray(new_dist_to_base_X)
 126 f.setName("dist")
 127 f.writeVTK("distEig.vtu")"""
 128 #mat*w[:,0]
 129
 130 print("{:g}".format(ForceNormale[0,0]))