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 fname = "slice.med"
  45 mm=MEDFileMesh.New(fname)
  46 m=mm[0]
  47 f1ts = MEDFileField1TS(fname,"RESUME__SIEF_NOEU")
  48 f = f1ts.field(mm)
  49 m = f.getMesh()
  50 area_field = m.getMeasureField(True)
  51 area = area_field.accumulate()[0] # 1
  52 centerOfMassField = m.computeCellCenterOfMass()
  53 centerOfMass = DataArrayDouble([elt/area for elt in (centerOfMassField*area_field.getArray()).accumulate()],1,3) # 2
  54 m.unPolyze()
  55 tri = MEDCoupling1SGTUMesh(m)
  56 assert(tri.getCellModelEnum()==NORM_TRI3)
  57 #
  58 fCell = f.nodeToCellDiscretization()
  59 (fCell.getArray()[:,[0,1,2]]*area).accumulate()
  60 ids = area_field.getArray().findIdsLowerThan(1e-7).buildComplement(m.getNumberOfCells())
  61 #fCell[ids]
  62 eqn = m[ids].computePlaneEquationOf3DFaces()[:,:3]
  63 eqn /= eqn.magnitude()
  64 area_vector = eqn*area_field.getArray()[ids]
  65 matrix = fCell[ids].getArray()
  66 #
  67 F_x = matrix[:,0]*eqn[:,0] + matrix[:,3]*eqn[:,1] + matrix[:,4]*eqn[:,2]
  68 F_y = matrix[:,3]*eqn[:,0] + matrix[:,1]*eqn[:,1] + matrix[:,5]*eqn[:,2]
  69 F_z = matrix[:,4]*eqn[:,0] + matrix[:,5]*eqn[:,1] + matrix[:,2]*eqn[:,2]
  70 #
  71 F = DataArrayDouble.Meld([F_x,F_y,F_z])
  72 #
  73 ZeForce = DataArrayDouble(F.accumulate(),1,3)
  74 normalFace = DataArrayDouble(eqn.accumulate(),1,3)
  75 normalFace /= normalFace.magnitude()[0]
  76 ForceNormale = DataArrayDouble.Dot(ZeForce,normalFace)[0]*normalFace # 3
  77 TangentForce = ZeForce-ForceNormale # 4
  78 #
  79 bary = m[ids].computeCellCenterOfMass()
  80 OM = bary-centerOfMass
  81 momentum = DataArrayDouble(DataArrayDouble.CrossProduct(OM,F).accumulate(),1,3) # 5
  82 # Inertie
  83 InertiaNormale = (DataArrayDouble.Dot(OM,OM)*area_field.getArray()[ids]).accumulate()[0] # 6_A
  84 base = DataArrayDouble(DataArrayDouble.GiveBaseForPlane(normalFace),3,3)
  85 angle, tangentPrincipal, inertiaPrincipal = fff(base[0],normalFace,OM,area_field.getArray()[ids])
  86 tangentOther = DataArrayDouble.CrossProduct(normalFace,tangentPrincipal)
  87 inertiaOther = ReturnInertia(tangentOther,OM,area_field.getArray()[ids])
  88 """
  89 base_X = DataArrayDouble(len(ids),3) ; base_X[:]=base[0]
  90 base_Y = DataArrayDouble(len(ids),3) ; base_Y[:]=base[1]
  91 dist_to_base_X = (OM-DataArrayDouble.Dot(OM,base_X)*base_X).magnitude()
  92 dist_to_base_Y = (OM-DataArrayDouble.Dot(OM,base_Y)*base_Y).magnitude()
  93 inertia_mat_0 = (dist_to_base_Y*dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
  94 inertia_mat_1 = (dist_to_base_X*dist_to_base_X*area_field.getArray()[ids]).accumulate()[0]
  95 inertia_mat_01 = -(dist_to_base_X*dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
  96 from numpy import linalg as LA
  97 import numpy as np
  98 mat = np.matrix([[inertia_mat_0, inertia_mat_01], [inertia_mat_01, inertia_mat_1]])
  99 v,w = LA.eig(mat)
 100 pos_of_max = max([(i,elt) for (i,elt) in enumerate(v)],key=lambda x: x[1])[0]
 101 u0 = DataArrayDouble(np.array(w[:,pos_of_max])) ; u0.rearrange(2)
 102 v0 = DataArrayDouble(np.array(w[:,1-pos_of_max])) ; v0.rearrange(2)
 103 #
 104 I_new_base_0 = v[pos_of_max] # 6_B
 105 new_base_0 = u0[0,0]*base[0]+u0[0,1]*base[1] # 6_B
 106 #new_base_1 = v0[0,0]*base[0]+v0[0,1]*base[1]
 107 new_base_1 = DataArrayDouble.CrossProduct(normalFace,new_base_0)
 108 new_base_Y = DataArrayDouble(len(ids),3) ; new_base_Y[:]=new_base_1
 109 new_dist_to_base_Y = (OM-DataArrayDouble.Dot(OM,new_base_Y)*new_base_Y).magnitude()
 110 I_new_base_1 = (new_dist_to_base_Y*new_dist_to_base_Y*area_field.getArray()[ids]).accumulate()[0]
 111 """
 112 """
 113 new_base_X = DataArrayDouble(len(ids),3) ; new_base_X[:]=new_base_0
 114 new_dist_to_base_X = (OM-DataArrayDouble.Dot(OM,new_base_X)*new_base_X).magnitude()
 115 I_new_base_0 = (new_dist_to_base_X*new_dist_to_base_X*area_field.getArray()[ids]).accumulate()[0]
 116 tmp=m[ids] ; tmp.zipCoords()
 117 f=MEDCouplingFieldDouble(ON_CELLS)
 118 f.setMesh(tmp)
 119 f.setArray(new_dist_to_base_X)
 120 f.setName("dist")
 121 f.writeVTK("distEig.vtu")"""
 122 #mat*w[:,0]