avoid computing face centers two times

internalLoads.py

@@ -102,35 +102,6 @@ def lecture_vtp(filename):
     points=lecture_vtk(os.path.splitext(filename)[0]+'.vtk')
     return points
 
-def cellPoints(output, CoR, angles):    
-    # this function returns the cell points in the body fixed frame
-    
-    numberOfCells=output.GetNumberOfCells()
-    cellIds = vtk.vtkIdList() # cell ids store to
-    
-    faceIDs = []
-    faceCfs = []
-    for cellIndex in range(numberOfCells): # for every cell
-       output.GetCellPoints(cellIndex, cellIds) # get ids of points of the given cell
-       center = np.zeros(3)
-       points = []
-       for i in range(0, cellIds.GetNumberOfIds()): # for every points of the given cell
-           coord = output.GetPoint(cellIds.GetId(i))
-           coord = np.array([coord[0],coord[1],coord[2]])
-           coord = np.dot( np.transpose(lib.rotation_matrix(angles[0],angles[1],angles[2])),\
-                          coord-CoR)
-           center += coord
-       center = center/cellIds.GetNumberOfIds()
-       if cellIndex == 0:
-           faceCfs = center.reshape(1,3)
-       else:
-           faceCfs=np.vstack(( faceCfs, center.reshape(1,3) )) # get coordinate, type : class tuple
-       faceIDs = np.append(faceIDs, cellIndex)
-       # cellPoints[cellIndex] = {}
-       # cellPoints[cellIndex]['Points'] = points
-       # cellPoints[cellIndex]['Center'] = center/cellIds.GetNumberOfIds()
-           
-    return faceIDs, faceCfs
 
 folder = '/home/gdieci2021/Projects/VBM_results_sample/' #os.getcwd()
 pvd_name = 'patchFields/'+\
@@ -179,7 +150,7 @@ with open(dissMassFile, 'r') as f:
                 count = 0
         k += 1
 
-# sections['0']['REFPOINT'][0] = -100
+sections['0']['REFPOINT'][0] = 100
 
 print('Initialisation : reading cell centers at first time step \n')
 
@@ -211,7 +182,7 @@ angles = motionData[ini,4:7].reshape(3)
 # motion in R0 and angles in Rb
 CoR = motion + CoRInitial # CoR in R0
 
-faceIDs, faceCfs=cellPoints(output, CoR, angles) # TODO check reference frame transformation, start from fixed and transform total force
+# faceIDs, faceCfs=cellPoints(output, CoR, angles) # TODO check reference frame transformation, start from fixed and transform total force
 
 # alpha=np.array(output.GetCellData().GetArray("alpha.water")) ## les valeurs au centre des cellules
 
@@ -230,7 +201,7 @@ print('Initialisation : selection of cells within section cut')
 # faceCfs = np.loadtxt('tmp/faceCfs.txt')
 # faceSf = np.loadtxt('tmp/faceSf.txt') 
 
-faceIDs, faceWeights, faceCfs, faceSf = lib.find_and_keep_faces(output, cutPoint, faceCfs, CoR, angles) # TODO check reference frame transformation, start from fixed and transform total force
+faceIDs, faceWeights, faceCfs, faceSf = lib.find_and_keep_faces(output, cutPoint, CoR, angles) # TODO check reference frame transformation, start from fixed and transform total force
 
 
 # with open('tmp/faceIDs.txt', 'w') as f:
@@ -291,7 +262,8 @@ for dt in range(0,len(time_foamstar)):
 
     motion    = motionData[ii,1:4].reshape(3)
     # linearAcc = motionData[ii,13:16].reshape(3)
-    angles    = motionData[ii,4:7].reshape(3)
+    angles    = motionData[ii,4:7].reshape(3)   
+    
     omg       = motionData[ii,10:13].reshape(3) 
     # angAcc    = motionData[ii,16:].reshape(3)
     
@@ -395,11 +367,11 @@ ax = fig.add_subplot(111, projection='3d')
 ax.plot3D(faceCfs[:,0], faceCfs[:,1], faceCfs[:,2], 'bo')
 ax.set_box_aspect([ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in 'xyz')])
 
-# res = os.path.join(folder, 'fFluid.dat')
+# res = os.path.join(folder, 'fz.dat')
 # res = np.loadtxt(res, skiprows=4)
 
 # t = res[:, 0]
-# fz = res[:,15]
+# fz = res[:,5]
 
 # plt.figure()
 # plt.plot(time, -shearSection[:,2], 'r--')
@@ -426,22 +398,22 @@ ax.set_box_aspect([ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in
 # plt.xlabel('Time (s)')
 # plt.ylabel('Moment (N.m)')
 
-res = os.path.join(folder, 'fFluid.dat')
+res = os.path.join(folder, 'forces.dat')
 res = np.loadtxt(res, skiprows=3)
 
 t = res [:, 0]
-forces = res[:,13:16]
+forces = res[:,1:4]
 
-# for i in range(0, len(t)):
+for i in range(0, len(t)):
     
-#     angles    = motionData[i,4:7].reshape(3)
+    angles    = motionData[i,4:7].reshape(3)
 
-#     P0b = lib.rotation_matrix(angles[0],angles[1],angles[2])
+    P0b = lib.rotation_matrix(angles[0],angles[1],angles[2])
     
-#     forces[i,:] = np.dot( np.transpose(P0b) , forces[i,:] )
+    forces[i,:] = np.dot( np.transpose(P0b) , forces[i,:] )
 
 plt.figure()
-plt.plot(time, -shearSection[:,2]/2, 'r--')
+plt.plot(time, shearSection[:,2]/2, 'r--')
 plt.plot(t, forces[:,2], 'k--')
 plt.xlim(min(time), max(time))
 plt.grid()

libs.py

@@ -52,7 +52,7 @@ def polygon_area(vertices):
     for i in range(0,len(vertices[:,0])-1):
         
         if abs(np.dot(vertices[i+1,:]-vertices[i,:], n)) > 1e-9:
-            print("Coplanar error of face vertices ! err = ", abs(np.dot(vertices[i+1,:]-vertices[i,:], n)))
+            #print("Coplanar error of face vertices ! err = ", abs(np.dot(vertices[i+1,:]-vertices[i,:], n)))
             error = 1
             
     if error == 1:
@@ -210,7 +210,7 @@ def cutPlane(vertex, cutPoint):
     totalSf = polygon_area(vertex)
 
     for i in range(0, len(vertex[:,0])):
-        if vertex[i,0]>xCut: #vertex[i,2]>zCut:
+        if vertex[i,0]<xCut: #vertex[i,2]>zCut:
             vertex[i,0] = xCut #vertex[i,2] = zCut
 
         subCf += vertex[i,:]
@@ -223,7 +223,7 @@ def cutPlane(vertex, cutPoint):
         
     
 
-def find_and_keep_faces(output, cutPoint, faceCfs, CoR, angles):
+def find_and_keep_faces(output, cutPoint, CoR, angles):
 
     # zCut = cutPoint[2]  # Assuming cutPoint is a 3D vector
     xCut = cutPoint[0]  # Assuming cutPoint is a 3D vector
@@ -242,17 +242,20 @@ def find_and_keep_faces(output, cutPoint, faceCfs, CoR, angles):
         output.GetCellPoints(faceI, cellIds) # get ids of points of the given cell
         count = 0
         vertex = []
+        center = np.zeros(3)
         for i in range(0, cellIds.GetNumberOfIds()):
             coord = output.GetPoint(cellIds.GetId(i))
             coord = np.array([coord[0],coord[1],coord[2]])
             coord = np.dot( np.transpose(rotation_matrix(angles[0],angles[1],angles[2])),\
                           coord-CoR )
+            
+            center += coord
                 
             if i == 0 :
                 vertex = coord.reshape(1,3)
             else:
                 vertex = np.vstack((vertex, coord.reshape(1,3)))
-            if coord[0]>xCut: #coord[2]<zCut:
+            if coord[0]<xCut: #coord[2]<zCut:
                 count += 1
         if count > 0:
             
@@ -262,10 +265,12 @@ def find_and_keep_faces(output, cutPoint, faceCfs, CoR, angles):
     
                 faceWeights=np.append(faceWeights, 1.0)
                 
+                center = center/cellIds.GetNumberOfIds()
+                
                 if len(faceCf)>0:
-                    faceCf=np.vstack((faceCf, faceCfs[faceI]))
+                    faceCf=np.vstack((faceCf, center))
                 else:
-                    faceCf=faceCfs[faceI]
+                    faceCf=center
                     
                 if len(faceSf)>0:
                     aux = compute_normal(vertex)*polygon_area(vertex)