Module pf_to_dem

  # -*- encoding=utf-8 -*-

  import pickle, math, os, shutil
  from pathlib import Path
  import numpy as np
  import matplotlib.pyplot as plt
  import vtk
  from vtk.util.numpy_support import vtk_to_numpy

  # own
  import tools

  # -----------------------------------------------------------------------------#

  def index_to_str(j):
      '''
      Convert a integer into a string with the format XXX.
      '''
      if j < 10:
          return '00'+str(j)
      elif 10<=j and j < 100:
          return '0'+str(j)
      else :
          return str(j)

  # -----------------------------------------------------------------------------#

  def read_vtk(dict_user, dict_sample, j_str):
      '''
      Read the last vtk files to obtain data from MOOSE.

      Do not work calling yade.
      '''
      eta_1_map_old = dict_sample['eta_1_map'].copy()
      c_map_old = dict_sample['c_map'].copy()
      L_eta1 = []
      L_c = []
      if not dict_sample['Map_known']:
          L_limits = []
          L_XYZ = []
          L_L_i_XYZ_used = []
      else :
          L_XYZ = dict_sample['L_XYZ']

      # iterate on the proccessors used
      for i_proc in range(dict_user['n_proc']):

          # name of the file to load
          namefile = 'vtk/pf_other_'+j_str+'_'+str(i_proc)+'.vtu'

          # load a vtk file as input
          reader = vtk.vtkXMLUnstructuredGridReader()
          reader.SetFileName(namefile)
          reader.Update()

          # Grab a scalar from the vtk file
          nodes_vtk_array = reader.GetOutput().GetPoints().GetData()
          eta1_vtk_array = reader.GetOutput().GetPointData().GetArray("eta1")
          c_vtk_array = reader.GetOutput().GetPointData().GetArray("c")

          #Get the coordinates of the nodes and the scalar values
          nodes_array = vtk_to_numpy(nodes_vtk_array)
          eta1_array = vtk_to_numpy(eta1_vtk_array)
          c_array = vtk_to_numpy(c_vtk_array)

          # map is not know
          if not dict_sample['Map_known']:
              # look for limits
              x_min = None
              x_max = None
              y_min = None
              y_max = None
              # save the map
              L_i_XYZ_used = []
              # Must detect common zones between processors
              for i_XYZ in range(len(nodes_array)) :
                  XYZ = nodes_array[i_XYZ]
                  # Do not consider twice a point
                  if list(XYZ) not in L_XYZ :
                      L_XYZ.append(list(XYZ))
                      L_eta1.append(eta1_array[i_XYZ])
                      L_c.append(c_array[i_XYZ])
                      L_i_XYZ_used.append(i_XYZ)
                      # set first point
                      if x_min == None :
                          x_min = list(XYZ)[0]
                          x_max = list(XYZ)[0]
                          y_min = list(XYZ)[1]
                          y_max = list(XYZ)[1]
                      # look for limits of the processor
                      else :
                          if list(XYZ)[0] < x_min:
                              x_min = list(XYZ)[0]
                          if list(XYZ)[0] > x_max:
                              x_max = list(XYZ)[0]
                          if list(XYZ)[1] < y_min:
                              y_min = list(XYZ)[1]
                          if list(XYZ)[1] > y_max:
                              y_max = list(XYZ)[1]
              # Here the algorithm can be help as the mapping is known
              L_L_i_XYZ_used.append(L_i_XYZ_used)
              # save limits
              L_limits.append([x_min,x_max,y_min,y_max])

          # map is known
          else :
              # the last term considered is at the end of the list
              if dict_sample['L_L_i_XYZ_used'][i_proc][-1] == len(nodes_array)-1:
                  L_eta1 = list(L_eta1) + list(eta1_array)
                  L_c = list(L_c) + list(c_array)
              # the last term considered is not at the end of the list
              else :
                  L_eta1 = list(L_eta1) + list(eta1_array[dict_sample['L_L_i_XYZ_used'][i_proc][0]: dict_sample['L_L_i_XYZ_used'][i_proc][-1]+1])
                  L_c = list(L_c) + list(c_array[dict_sample['L_L_i_XYZ_used'][i_proc][0]: dict_sample['L_L_i_XYZ_used'][i_proc][-1]+1])

      if not dict_sample['Map_known']:
          # plot processors distribution
          if 'processor' in dict_user['L_figures'] and not dict_user['remesh']:
              fig, (ax1) = plt.subplots(1,1,figsize=(16,9))
              # parameters
              title_fontsize = 20
              for i_proc in range(len(L_limits)):
                  limits = L_limits[i_proc]
                  ax1.plot([limits[0],limits[1],limits[1],limits[0],limits[0]],[limits[2],limits[2],limits[3],limits[3],limits[2]], label='proc '+str(i_proc))
              ax1.legend()
              ax1.set_xlabel('X (m)')
              ax1.set_ylabel('Y (m)')
              ax1.set_title('Processor i has the priority on i+1',fontsize = title_fontsize)
              fig.suptitle('Processors ditribution',fontsize = 1.2*title_fontsize)
              fig.tight_layout()
              fig.savefig('plot/processors_distribution.png')
              plt.close(fig)
          # the map is known
          dict_sample['Map_known'] = True
          dict_sample['L_L_i_XYZ_used'] = L_L_i_XYZ_used
          dict_sample['L_XYZ'] = L_XYZ

      # rebuild map from lists
      for i_XYZ in range(len(L_XYZ)):
          # find nearest node
          L_search = list(abs(np.array(dict_sample['x_L']-L_XYZ[i_XYZ][0])))
          i_x = L_search.index(min(L_search))
          L_search = list(abs(np.array(dict_sample['y_L']-L_XYZ[i_XYZ][1])))
          i_y = L_search.index(min(L_search))
          # rewrite map
          dict_sample['eta_1_map'][-1-i_y, i_x] = L_eta1[i_XYZ]
          dict_sample['c_map'][-1-i_y, i_x] = L_c[i_XYZ]

  # -----------------------------------------------------------------------------#

  def compute_vertices(dict_user, dict_sample):
      '''
      From a phase map, compute vertices coordinates for polyhedral.
      '''
      # compute vertices for grain
      L_vertices_1 = interpolate_vertices(dict_sample['eta_1_map'], dict_sample['pos_1'], dict_user, dict_sample)

      # compute vertices for plate
      L_vertices_2 = ((min(dict_sample['x_L']), 0, 0,),
                      (max(dict_sample['x_L']), 0, 0,),
                      (max(dict_sample['x_L']), min(dict_sample['y_L']), 0,),
                      (min(dict_sample['x_L']), min(dict_sample['y_L']), 0,),
                      (min(dict_sample['x_L']), 0, 1,),
                      (max(dict_sample['x_L']), 0, 1,),
                      (max(dict_sample['x_L']), min(dict_sample['y_L']), 1,),
                      (min(dict_sample['x_L']), min(dict_sample['y_L']), 1,))

      # save data
      dict_save = {
      'L_vertices_1': L_vertices_1,
      'L_vertices_2': L_vertices_2
      }
      with open('data/planes.data', 'wb') as handle:
          pickle.dump(dict_save, handle, protocol=pickle.HIGHEST_PROTOCOL)

      # compute sample height and vertical strain
      L_x_1, L_y_1 = tools.tuplet_to_list(L_vertices_1) # from tools.py
      sample_height = max(L_y_1)
      dict_user['L_sample_height'].append(sample_height)

      # compute sphericities
      AreaSphericity, DiameterSphericity, CircleRatioSphericity, PerimeterSphericity, WidthToLengthRatioSpericity = tools.compute_sphericities(L_vertices_1) # from tools.py
      # save and plot
      dict_user['L_AreaSphericity'].append(AreaSphericity)
      dict_user['L_DiameterSphericity'].append(DiameterSphericity)
      dict_user['L_CircleRatioSphericity'].append(CircleRatioSphericity)
      dict_user['L_PerimeterSphericity'].append(PerimeterSphericity)
      dict_user['L_WidthToLengthRatioSpericity'].append(WidthToLengthRatioSpericity)

      if 'sphericities' in dict_user['L_figures']:
          fig, (ax1) = plt.subplots(1,1,figsize=(16,9))
          ax1.plot(dict_user['L_AreaSphericity'], label='Area sphericity')
          ax1.plot(dict_user['L_DiameterSphericity'], label='Diameter sphericity')
          ax1.plot(dict_user['L_CircleRatioSphericity'], label='Circle ratio sphericity')
          ax1.plot(dict_user['L_PerimeterSphericity'], label='Perimeter sphericity')
          ax1.plot(dict_user['L_WidthToLengthRatioSpericity'], label='Width/Length spericity')
          ax1.legend()
          plt.suptitle('Grain sphericities (G1)', fontsize=20)
          fig.tight_layout()
          fig.savefig('plot/sphericities.png')
          plt.close(fig)


  # -----------------------------------------------------------------------------#

  def interpolate_vertices(eta_i_map, center, dict_user, dict_sample):
      '''
      Interpolate vertices for polyhedral.
      '''
      map_phi = []
      L_phi = []
      for i_phi in range(dict_user['n_phi']):
          phi = 2*math.pi*i_phi/dict_user['n_phi']
          L_phi.append(phi)
          map_phi.append([])
      L_phi.append(2*math.pi)
      for i_x in range(len(dict_sample['x_L'])-1):
          for i_y in range(len(dict_sample['y_L'])-1):
              L_in = [] # list the nodes inside the grain
              if eta_i_map[-1-i_y    , i_x] > 0.5 :
                  L_in.append(0)
              if eta_i_map[-1-(i_y+1), i_x] > 0.5 :
                  L_in.append(1)
              if eta_i_map[-1-(i_y+1), i_x+1] > 0.5 :
                  L_in.append(2)
              if eta_i_map[-1-i_y    , i_x+1] > 0.5 :
                  L_in.append(3)
              if L_in != [] and L_in != [0,1,2,3]:
                  center_mesh = (np.array([dict_sample['x_L'][i_x], dict_sample['y_L'][i_y]])+np.array([dict_sample['x_L'][i_x+1], dict_sample['y_L'][i_y+1]]))/2
                  u = (center_mesh-np.array(center))/np.linalg.norm(center_mesh-np.array(center))
                  # compute phi
                  if u[1]>=0:
                      phi = math.acos(u[0])
                  else :
                      phi = 2*math.pi-math.acos(u[0])
                  # iterate on the lines of the mesh to find the plane intersection
                  L_p = []
                  if (0 in L_in and 1 not in L_in) or (0 not in L_in and 1 in L_in):# line 01
                      x_p = dict_sample['x_L'][i_x]
                      y_p = (0.5-eta_i_map[-1-i_y, i_x])/(eta_i_map[-1-(i_y+1), i_x]-eta_i_map[-1-i_y, i_x])*(dict_sample['y_L'][i_y+1]-dict_sample['y_L'][i_y])+dict_sample['y_L'][i_y]
                      L_p.append(np.array([x_p, y_p]))
                  if (1 in L_in and 2 not in L_in) or (1 not in L_in and 2 in L_in):# line 12
                      x_p = (0.5-eta_i_map[-1-(i_y+1), i_x])/(eta_i_map[-1-(i_y+1), i_x+1]-eta_i_map[-1-(i_y+1), i_x])*(dict_sample['x_L'][i_x+1]-dict_sample['x_L'][i_x])+dict_sample['x_L'][i_x]
                      y_p = dict_sample['y_L'][i_y+1]
                      L_p.append(np.array([x_p, y_p]))
                  if (2 in L_in and 3 not in L_in) or (2 not in L_in and 3 in L_in):# line 23
                      x_p = dict_sample['x_L'][i_x+1]
                      y_p = (0.5-eta_i_map[-1-i_y, i_x+1])/(eta_i_map[-1-(i_y+1), i_x+1]-eta_i_map[-1-i_y, i_x+1])*(dict_sample['y_L'][i_y+1]-dict_sample['y_L'][i_y])+dict_sample['y_L'][i_y]
                      L_p.append(np.array([x_p, y_p]))
                  if (3 in L_in and 0 not in L_in) or (3 not in L_in and 0 in L_in):# line 30
                      x_p = (0.5-eta_i_map[-1-i_y, i_x])/(eta_i_map[-1-i_y, i_x+1]-eta_i_map[-1-i_y, i_x])*(dict_sample['x_L'][i_x+1]-dict_sample['x_L'][i_x])+dict_sample['x_L'][i_x]
                      y_p = dict_sample['y_L'][i_y]
                      L_p.append(np.array([x_p, y_p]))
                  # compute the mean point
                  p_mean = np.array([0,0])
                  for p in L_p :
                      p_mean = p_mean + p
                  p_mean = p_mean/len(L_p)
                  # look phi in L_phi
                  i_phi = 0
                  while not (L_phi[i_phi] <= phi and phi <= L_phi[i_phi+1]) :
                      i_phi = i_phi + 1
                  # save p_mean in the map
                  map_phi[i_phi].append(p_mean)

      L_vertices = ()
      # interpolate plane (Least squares method)
      for i_phi in range(len(map_phi)):
          if map_phi[i_phi] != []:
              # mean vertices
              mean_v = np.array([0, 0])
              for v_i in map_phi[i_phi]:
                  mean_v = mean_v + v_i/len(map_phi[i_phi])
              # save
              L_vertices = L_vertices + ((mean_v[0], mean_v[1], 0,),)
              L_vertices = L_vertices + ((mean_v[0], mean_v[1], 1),)

      return L_vertices

  # -----------------------------------------------------------------------------#

  def control_force(dict_user, dict_sample):
      '''
      Control force applied in DEM with the contact volume.

      Because of the fact that convvex shape is not possible in Yade.
      '''
      # determine the volume contact in Moose (etai ad j > 0.5)
      V_contact_moose = 0
      for i_x in range(len(dict_sample['x_L'])):
          for i_y in range(len(dict_sample['y_L'])):
              if dict_sample['eta_1_map'][-1-i_y, i_x] > 0.5 and dict_sample['y_L'][i_y]<=0:
                  V_contact_moose = V_contact_moose + (dict_sample['x_L'][1]-dict_sample['x_L'][0])*(dict_sample['y_L'][1]-dict_sample['y_L'][0])*1

      # compare with the volume contact expected
      # in Yade, F=k*V and 1/k = 1/Y1 + 1/Y2
      V_expected = 2*dict_user['force_applied_target']/dict_user['E']
      # the comparison is important if multiple controls are done per one PFDEM iteration
      # if only only control is done, the adaptation of the force can be done always
      #if V_contact_moose < (1-dict_user['steady_state_detection'])*V_expected or (1+dict_user['steady_state_detection'])*V_expected < V_contact_moose:
      if True:
          # control the force applied
          dict_user['force_applied'] = dict_user['force_applied']*V_expected/V_contact_moose

      # save and plot
      dict_user['L_force_applied'].append(dict_user['force_applied'])
      if 'force_applied' in dict_user['L_figures']:
          fig, (ax1) = plt.subplots(1,1,figsize=(16,9))
          ax1.plot(dict_user['L_force_applied'])
          plt.suptitle('Force applied in Yade (DEM)', fontsize=20)
          fig.tight_layout()
          fig.savefig('plot/force_applied.png')
          plt.close(fig)