Re: Define Nodes in .comm file

  from Utilitai.partition import *
   from numpy import *
   #from math import *

   if info>1:
      print 'in CA_geeometry: '
      print mesh
      print ElemGroupStr

   PyMesh   = MAIL_PY()                 ## convert CodeAster mesh to Python
   PyMesh.FromAster(mesh);
   nonu     = PyMesh.dime_maillage[0]   ## number of nodes
   elnu     = PyMesh.dime_maillage[2]   ## number of elements
   NodeCoord    = PyMesh.cn             ## xyz coordinates of nodes (nonu*3 matrix)
   ElemConnect  = PyMesh.co             ## node numbers of elements (elnu*x matrix)
   NodeList = list(PyMesh.correspondance_noeuds)
   ElemList = list(PyMesh.correspondance_mailles)
   ElemType = list(PyMesh.tm)
   NodeGroup = PyMesh.gno               ## list of GROUP_NO groups
   ElemGroup = PyMesh.gma               ## list of GROUP_MA groups, eg groupstr

   #NodeList = list(CmeshCA.correspondance_noeuds)
   #ElemList = list(CmeshCA.correspondance_mailles)
   #ElemType = list(CmeshCA.tm)

   if info>=0:
     print 'groups by argument:%s'%(ElemGroupStr)
     print ElemGroup[ElemGroupStr]
     
     #P1 = numpy.array([0.0,2.0,3.0])   ## p1 and p2 are points defining the rotation axis
     #P2 = numpy.array([2.0,2.0,3.0])   ## define your own rotation axis in comm file
     # ideally the rotation axis is derived from the set of nodes --> to be done

   NodeCount = 0 
   TorsionDistribution=[]                   ## initialise Torsion node list group
   #NumberOfNodes=len(NodeGroup[ElemGroupStr])
   NumberOfElems=len(ElemGroup[ElemGroupStr])

   ElemNodeGroup = [0]*4*NumberOfElems
   FXX = [0]*4*NumberOfElems
   FYY = [0]*4*NumberOfElems
   ElemCount = 0
   
   print 'ElemNodeGroup: ',ElemNodeGroup,len(ElemNodeGroup)
   
   for ii in xrange(NumberOfElems):
      ElemNumber = ElemGroup[ElemGroupStr][ii]
      print 'ElemNumber: ',ii,ElemCount,ElemNumber
      ElemNode0 = ElemConnect[ElemNumber][0]
      ElemNode1 = ElemConnect[ElemNumber][1]
      ElemNode2 = ElemConnect[ElemNumber][2]
      ElemNode3 = ElemConnect[ElemNumber][3]
      print 'index: ',ElemCount,552,ii,NumberOfElems
      ElemNodeGroup[0+ElemCount]=ElemNode0
      ElemNodeGroup[1+ElemCount]=ElemNode1
      ElemNodeGroup[2+ElemCount]=ElemNode2
      ElemNodeGroup[3+ElemCount]=ElemNode3
      ElemCount+=4

   print 'ElemNodeGroup: ',ElemNodeGroup,len(ElemNodeGroup)
   #np.unique1d([1,2,6,4,2,3,2], return_index=True)
   #print 'unique: ',unique([1,2,6,4,2,3,2])
   ElemNodeGroupUnique = unique(ElemNodeGroup)
   print 'unique elemNodeGroup: ',ElemNodeGroup
   print 'unique elemNodeGroup: ',ElemNodeGroupUnique
   NumberOfNodesGroup = len(ElemNodeGroupUnique)
 
   TorsionDistribution=[]                   ## initialise Torsion node list group
   ElemCount = 0
   SumArea = 0.0
   for ii in xrange(NumberOfElems):
      ElemNumber = ElemGroup[ElemGroupStr][ii]
      print 'ElemNumber: ',ii,ElemCount,ElemNumber
      ElemNode0 = ElemConnect[ElemNumber][0]
      ElemNode1 = ElemConnect[ElemNumber][1]
      ElemNode2 = ElemConnect[ElemNumber][2]
      ElemNode3 = ElemConnect[ElemNumber][3]
      #ElemNode4 = ElemConnect[ElemNumber][4]
      #ElemNode5 = ElemConnect[ElemNumber][5]
      #ElemNode6 = ElemConnect[ElemNumber][6]
      #ElemNode7 = ElemConnect[ElemNumber][7]
      #print 'ElemNode1: ',ElemNode0,ElemNode1,ElemNode2,ElemNode3,ElemNode4,ElemNode5,ElemNode6,ElemNode7
      print 'ElemNode1: ',ElemNode0,ElemNode1,ElemNode2,ElemNode3
      
      NNxyz0 = NodeCoord[ElemNode0]
      NNxyz1 = NodeCoord[ElemNode1]
      NNxyz2 = NodeCoord[ElemNode2]
      NNxyz3 = NodeCoord[ElemNode3]
      print 'ElemCoordinates: ',NNxyz0,NNxyz1,NNxyz2,NNxyz3
      
      
			# determine area of element surface A~||(x2-x0).(x3-x1)||
      dx20 = NNxyz2[0]-NNxyz0[0]
      dy20 = NNxyz2[1]-NNxyz0[1]
      dz20 = NNxyz2[2]-NNxyz0[2]

      dx31 = NNxyz3[0]-NNxyz1[0]
      dy31 = NNxyz3[1]-NNxyz1[1]
      dz31 = NNxyz3[2]-NNxyz1[2]

      aa = dy20*dz31-dz20*dy31   ## cross product of diagonals
      bb = dz20*dx31-dx20*dz31   
      cc = dx20*dy31-dy20*dx31   
           
      AreaElem = 0.5*sqrt(aa*aa+bb*bb+cc*cc)  ## Area = 0.5*||cross(d1,d2))||
      ##AreaElem = 1.000 ## trick size of elements
      
      SumArea += AreaElem
      print 'AreaElem: ',AreaElem,SumArea
      
      #NNp1 = NodeNumber+1         # node numbers are increased by 1: Salome <--> Code Aster
      xcoord = NNxyz0[0]
      ycoord = NNxyz0[1]
      print xcoord,ycoord,AreaElem,Ftorsion
      
      Fx = -AreaElem*ycoord
      Fy = +AreaElem*xcoord
      FXX[ElemCount] = Fx 
      FYY[ElemCount] = Fy

      xcoord = NNxyz1[0]
      ycoord = NNxyz1[1]
      #print xcoord,ycoord,AreaElem,Ftorsion
      Fx = -AreaElem*ycoord
      Fy = +AreaElem*xcoord
      FXX[ElemCount+1] = Fx 
      FYY[ElemCount+1] = Fy

      xcoord = NNxyz2[0]
      ycoord = NNxyz2[1]
      #print xcoord,ycoord,AreaElem,Ftorsion
      Fx = -AreaElem*ycoord
      Fy = +AreaElem*xcoord
            
      FXX[ElemCount+2] = Fx 
      FYY[ElemCount+2] = Fy

      xcoord = NNxyz3[0]
      ycoord = NNxyz3[1]
      #print xcoord,ycoord,AreaElem,Ftorsion
      Fx = -AreaElem*ycoord
      Fy = +AreaElem*xcoord
      FXX[ElemCount+3] = Fx 
      FYY[ElemCount+3] = Fy

      #print 'eCount,Node1234,FXX/FYY: ',ElemCount+0,ElemNode0,NNxyz0[0],NNxyz0[1],FXX[ElemCount+0],FYY[ElemCount+0]
      #print 'eCount,Node1234,FXX/FYY: ',ElemCount+1,ElemNode1,NNxyz1[0],NNxyz1[1],FXX[ElemCount+1],FYY[ElemCount+1]
      #print 'eCount,Node1234,FXX/FYY: ',ElemCount+2,ElemNode2,NNxyz2[0],NNxyz2[1],FXX[ElemCount+2],FYY[ElemCount+2]
      #print 'eCount,Node1234,FXX/FYY: ',ElemCount+3,ElemNode3,NNxyz3[0],NNxyz3[1],FXX[ElemCount+3],FYY[ElemCount+3]

      ElemCount+=4
      #print 'xxxxxxxxxxxxxxxxxxxxx'

   #print 'xxxxxxxxxxxxxxxxxxxxx'
   ##for ii in xrange(len(FYY)):
   #3   print 'ii/Node/FXX/FYY: ',ii,ElemNodeGroup[ii],FXX[ii],FYY[ii]