Package Owntools
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
This file contains the different functions used in the simulation.
Expand source code
# -*- coding: utf-8 -*-
"""
@author: Alexandre Sac--Morane
alexandre.sac-morane@uclouvain.be
This file contains the different functions used in the simulation.
"""
#-------------------------------------------------------------------------------
#Librairy
#-------------------------------------------------------------------------------
import os
import math
import numpy as np
from pathlib import Path
from scipy.ndimage import binary_dilation
#-------------------------------------------------------------------------------
def index_to_str(j):
'''
An integer is converted to a float with 3 components
Input :
a number (a float)
Output :
a string with 3 components (a string)
'''
if j < 10:
j_str = '00'+str(j)
elif 10 <= j and j < 100:
j_str = '0'+str(j)
else :
j_str = str(j)
return j_str
#-------------------------------------------------------------------------------
def Sort_Files(dict_algorithm):
'''
Sort files generated by MOOSE to different directories
Input :
an algorithm dictionnary (a dict)
Output :
Nothing but files are sorted
'''
os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_out.e','Output/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_out.e')
os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'.i','Input/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'.i')
j = 0
j_str = index_to_str(j)
filepath = Path(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu')
while filepath.exists():
for i_proc in range(dict_algorithm['np_proc']):
os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'_'+str(i_proc)+'.vtu','Output/Ite_'+str(dict_algorithm['i_PFDEM'])+'/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'_'+str(i_proc)+'.vtu')
os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu','Output/Ite_'+str(dict_algorithm['i_PFDEM'])+'/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu')
j = j + 1
j_str = index_to_str(j)
filepath = Path(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu')
return index_to_str(j-1)
#-------------------------------------------------------------------------------
def Extract_solute_at_p(dict_sample,ij_p):
'''
Extract the value of the solute concentration at a given point.
Input :
a sample dictionnary (a dict)
a coordinate of the point (a tuple of int)
Output :
the value of the solute concentration (a float)
'''
return dict_sample['solute_M'][-1-ij_p[0]][ij_p[1]]
#-------------------------------------------------------------------------------
def Cosine_Profile(R,r,w):
'''
Compute the phase field variable at some point.
A cosine profile is assumed (see https://mooseframework.inl.gov/source/ics/SmoothCircleIC.html).
Input :
the radius R of the grain in the direction (a float)
the distance r between the current point and the center (a float)
the width w of the interface (a float)
Output :
the value of the phase field variable (a float)
'''
#inside the grain
if r<R-w/2:
return 1
#outside the grain
elif r>R+w/2:
return 0
#inside the interface
else :
return 0.5*(1 + math.cos(math.pi*(r-R+w/2)/w))
#-------------------------------------------------------------------------------
def error_on_ymax_f(dy,overlap_L,k_L,Force_target) :
"""
Compute the function f to control the upper wall. It is the difference between the force applied and the target value.
Input :
an increment of the upper wall position (a float)
a list of overlap for contact between grain and upper wall (a list)
a list of spring for contact between grain and upper wall (a list)
a confinement force (a float)
Output :
the difference between the force applied and the confinement (a float)
"""
f = Force_target
for i in range(len(overlap_L)):
f = f - k_L[i]*(max(overlap_L[i]-dy,0))**(3/2)
return f
#-------------------------------------------------------------------------------
def error_on_ymax_df(dy,overlap_L,k_L) :
"""
Compute the derivative function df to control the upper wall (error_on_ymax_f()).
Input :
an increment of the upper wall position (a float)
a list of overlap for contact between grain and upper wall (a list)
a list of spring for contact between grain and upper wall (a list)
Output :
the derivative of error_on_ymax_f() (a float)
"""
df = 0
for i in range(len(overlap_L)):
df = df + 3/2*k_L[i]*(max(overlap_L[i]-dy,0))**(1/2)
return df
#-------------------------------------------------------------------------------
def Control_y_max_NR(dict_sample,dict_sollicitation):
"""
Control the upper wall to apply force.
A Newton-Raphson method is applied to verify the confinement.
Input :
a sample dictionnary (a dict)
a sollicitations dictionnary (a dict)
Output :
Nothing, but the sample dictionnary is updated, concerning the upper wall position and force applied (two floats)
"""
#-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
#load data needed
Force_target = dict_sollicitation['Vertical_Confinement_Force']
#-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
F = 0
overlap_L = []
k_L = []
for contact in dict_sample['L_contact_gw']:
if contact.nature == 'gwy_max':
F = F + contact.Fwg_n
overlap_L.append(contact.overlap)
k_L.append(contact.k)
#compute force applied, save contact overlap and spring
if overlap_L != []:
i_NR = 0
dy = 0
ite_criteria = True
if -0.01*Force_target<error_on_ymax_f(dy,overlap_L,k_L,Force_target) and error_on_ymax_f(dy,overlap_L,k_L,Force_target)<0.01*Force_target:
ite_criteria = False
while ite_criteria :
i_NR = i_NR + 1
dy = dy - error_on_ymax_f(dy,overlap_L,k_L,Force_target)/error_on_ymax_df(dy,overlap_L,k_L)
if i_NR > 100:
ite_criteria = False
if -0.01*Force_target<error_on_ymax_f(dy,overlap_L,k_L,Force_target) and error_on_ymax_f(dy,overlap_L,k_L,Force_target)<0.01*Force_target:
ite_criteria = False
dict_sample['y_box_max'] = dict_sample['y_box_max'] + dy
else :
#if there is no contact with the upper wall, the wall is reset
dict_sample['y_box_max'] = Reset_y_max(dict_sample['L_g'],Force_target)
for contact in dict_sample['L_contact_gw']:
if contact.nature == 'gwy_max':
#reactualisation
contact.limit = dict_sample['y_box_max']
#Update dict
dict_sollicitation['Force_on_upper_wall'] = F
#-------------------------------------------------------------------------------
def Reset_y_max(L_g,Force):
"""
The upper wall is located as a single contact verify the target value.
Input :
the list of temporary grains (a list)
the confinement force (a float)
Output :
the upper wall position (a float)
"""
print('Reset of the y_max on the upper grain')
y_max = None
id_grain_max = None
for id_grain in range(len(L_g)):
grain = L_g[id_grain]
y_max_grain = max(grain.l_border_y)
if y_max != None and y_max_grain > y_max:
y_max = y_max_grain
id_grain_max = id_grain
elif y_max == None:
y_max = y_max_grain
id_grain_max = id_grain
k = 5*4/3*L_g[id_grain_max].y/(1-L_g[id_grain_max].nu*L_g[id_grain_max].nu)*math.sqrt(L_g[id_grain_max].r_mean)
y_max = y_max - (Force/k)**(2/3)
return y_max
#-------------------------------------------------------------------------------
def Interpolate_solute_out_grains(dict_algorithm, dict_sample) :
"""
For all node of the mesh, if there is solute in one grain it is moved to the nearest node out of the grain.
A node in multiple grain (contact area) is not considered in a grain a solute can exist.
Input :
an algorithm dictionnary (a dict)
a sample dictionnary (a dict)
Output :
Nothing, but the solute map is updated (a ny x nx numpy array)
"""
#Create maps of different available nodes
node_available_M = np.array(np.ones((len(dict_sample['y_L']),len(dict_sample['x_L'])),dtype = bool))
for l in range(len(dict_sample['y_L'])):
for c in range(len(dict_sample['x_L'])):
n_grain = 0
for etai in dict_sample['L_etai']:
if etai.etai_M[l][c] > 0.5 :
n_grain = n_grain + 1
if n_grain == 1:
node_available_M[l][c] = False
#dilatation
dilated_M = binary_dilation(node_available_M, dict_algorithm['struct_element'])
#iteration on solute map to see if some solute is in not available position
for l in range(len(dict_sample['y_L'])):
for c in range(len(dict_sample['x_L'])):
if dict_sample['solute_M'][l][c] > 0 and not dilated_M[l][c]:
solute_moved = False
size_window = 1
while not solute_moved :
i_window = 0
while not solute_moved and i_window <= size_window:
n_node_available = 0
#Look to move horizontaly and vertically
if i_window == 0 :
top_available = False
down_available = False
left_available = False
right_available = False
#to the top
if l - size_window > 0:
top_available = dilated_M[l-size_window][c]
if dilated_M[l-size_window][c] :
n_node_available = n_node_available + 1
#to the down
if l + size_window < len(dict_sample['y_L']):
down_available = dilated_M[l+size_window][c]
if dilated_M[l+size_window][c] :
n_node_available = n_node_available + 1
#to the left
if c - size_window > 0:
left_available = dilated_M[l][c-size_window]
if dilated_M[l][c-size_window] :
n_node_available = n_node_available + 1
#to the right
if c + size_window < len(dict_sample['x_L']):
right_available = dilated_M[l][c+size_window]
if dilated_M[l][c+size_window] :
n_node_available = n_node_available + 1
#move solute if et least one node is available
if n_node_available != 0 :
#to the top
if top_available:
dict_sample['solute_M'][l-size_window][c] = dict_sample['solute_M'][l-size_window][c] + dict_sample['solute_M'][l][c]/n_node_available
#to the down
if down_available:
dict_sample['solute_M'][l+size_window][c] = dict_sample['solute_M'][l+size_window][c] + dict_sample['solute_M'][l][c]/n_node_available
#to the left
if left_available:
dict_sample['solute_M'][l][c-size_window] = dict_sample['solute_M'][l][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available
#to the right
if right_available:
dict_sample['solute_M'][l][c+size_window] = dict_sample['solute_M'][l][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available
dict_sample['solute_M'][l][c] = 0
solute_moved = True
#Look to move diagonally
else :
top_min_available = False
top_max_available = False
down_min_available = False
down_max_available = False
left_min_available = False
left_max_available = False
right_min_available = False
right_max_available = False
#to the top
if l - size_window > 0:
if c - i_window > 0 :
top_min_available = dilated_M[l-size_window][c-i_window]
if dilated_M[l-size_window][c-i_window] :
n_node_available = n_node_available + 1
if c + i_window < len(dict_sample['x_L']):
top_max_available = dilated_M[l-size_window][c+i_window]
if dilated_M[l-size_window][c+i_window] :
n_node_available = n_node_available + 1
#to the down
if l + size_window < len(dict_sample['y_L']):
if c - i_window > 0 :
down_min_available = dilated_M[l+size_window][c-i_window]
if dilated_M[l+size_window][c-i_window] :
n_node_available = n_node_available + 1
if c + i_window < len(dict_sample['x_L']):
down_max_available = dilated_M[l+size_window][c+i_window]
if dilated_M[l+size_window][c+i_window] :
n_node_available = n_node_available + 1
#to the left
if c - size_window > 0:
if l - i_window > 0 :
left_min_available = dilated_M[l-i_window][c-size_window]
if dilated_M[l-i_window][c-size_window] :
n_node_available = n_node_available + 1
if l + i_window < len(dict_sample['y_L']):
left_max_available = dilated_M[l+i_window][c-size_window]
if dilated_M[l+i_window][c-size_window] :
n_node_available = n_node_available + 1
#to the right
if c + size_window < len(dict_sample['x_L']):
if l - i_window > 0 :
right_min_available = dilated_M[l-i_window][c+size_window]
if dilated_M[l-i_window][c+size_window] :
n_node_available = n_node_available + 1
if l + i_window < len(dict_sample['y_L']):
right_max_available = dilated_M[l+i_window][c+size_window]
if dilated_M[l+i_window][c+size_window] :
n_node_available = n_node_available + 1
#move solute if et least one node is available
if n_node_available != 0 :
#to the top
if top_min_available:
dict_sample['solute_M'][l-size_window][c-i_window] = dict_sample['solute_M'][l-size_window][c-i_window] + dict_sample['solute_M'][l][c]/n_node_available
if top_max_available:
dict_sample['solute_M'][l-size_window][c+i_window] = dict_sample['solute_M'][l-size_window][c+i_window] + dict_sample['solute_M'][l][c]/n_node_available
#to the down
if down_min_available:
dict_sample['solute_M'][l+size_window][c-i_window] = dict_sample['solute_M'][l+size_window][c-i_window] + dict_sample['solute_M'][l][c]/n_node_available
if down_max_available:
dict_sample['solute_M'][l+size_window][c+i_window] = dict_sample['solute_M'][l+size_window][c+i_window] + dict_sample['solute_M'][l][c]/n_node_available
#to the left
if left_min_available:
dict_sample['solute_M'][l-i_window][c-size_window] = dict_sample['solute_M'][l-i_window][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available
if left_max_available:
dict_sample['solute_M'][l+i_window][c-size_window] = dict_sample['solute_M'][l+i_window][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available
#to the right
if right_min_available:
dict_sample['solute_M'][l-i_window][c+size_window] = dict_sample['solute_M'][l-i_window][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available
if right_max_available:
dict_sample['solute_M'][l+i_window][c+size_window] = dict_sample['solute_M'][l+i_window][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available
dict_sample['solute_M'][l][c] = 0
solute_moved = True
i_window = i_window + 1
size_window = size_window + 1Sub-modules
Owntools.Compute-
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
Functions to compute parameters or variables. Owntools.PFtoDEM_Multi-
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
Functions to convert PF data into DEM data. Owntools.Plot-
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
Functions to plot data. Owntools.Save-
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
Functions to save data. Owntools.Write-
@author: Alexandre Sac–Morane alexandre.sac-morane@uclouvain.be
Functions to write .txt or .i files.
Functions
def Control_y_max_NR(dict_sample, dict_sollicitation)-
Control the upper wall to apply force.
A Newton-Raphson method is applied to verify the confinement. Input : a sample dictionnary (a dict) a sollicitations dictionnary (a dict) Output : Nothing, but the sample dictionnary is updated, concerning the upper wall position and force applied (two floats)
Expand source code
def Control_y_max_NR(dict_sample,dict_sollicitation): """ Control the upper wall to apply force. A Newton-Raphson method is applied to verify the confinement. Input : a sample dictionnary (a dict) a sollicitations dictionnary (a dict) Output : Nothing, but the sample dictionnary is updated, concerning the upper wall position and force applied (two floats) """ #-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. #load data needed Force_target = dict_sollicitation['Vertical_Confinement_Force'] #-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-. F = 0 overlap_L = [] k_L = [] for contact in dict_sample['L_contact_gw']: if contact.nature == 'gwy_max': F = F + contact.Fwg_n overlap_L.append(contact.overlap) k_L.append(contact.k) #compute force applied, save contact overlap and spring if overlap_L != []: i_NR = 0 dy = 0 ite_criteria = True if -0.01*Force_target<error_on_ymax_f(dy,overlap_L,k_L,Force_target) and error_on_ymax_f(dy,overlap_L,k_L,Force_target)<0.01*Force_target: ite_criteria = False while ite_criteria : i_NR = i_NR + 1 dy = dy - error_on_ymax_f(dy,overlap_L,k_L,Force_target)/error_on_ymax_df(dy,overlap_L,k_L) if i_NR > 100: ite_criteria = False if -0.01*Force_target<error_on_ymax_f(dy,overlap_L,k_L,Force_target) and error_on_ymax_f(dy,overlap_L,k_L,Force_target)<0.01*Force_target: ite_criteria = False dict_sample['y_box_max'] = dict_sample['y_box_max'] + dy else : #if there is no contact with the upper wall, the wall is reset dict_sample['y_box_max'] = Reset_y_max(dict_sample['L_g'],Force_target) for contact in dict_sample['L_contact_gw']: if contact.nature == 'gwy_max': #reactualisation contact.limit = dict_sample['y_box_max'] #Update dict dict_sollicitation['Force_on_upper_wall'] = F def Cosine_Profile(R, r, w)-
Compute the phase field variable at some point.
A cosine profile is assumed (see https://mooseframework.inl.gov/source/ics/SmoothCircleIC.html).
Input : the radius R of the grain in the direction (a float) the distance r between the current point and the center (a float) the width w of the interface (a float) Output : the value of the phase field variable (a float)
Expand source code
def Cosine_Profile(R,r,w): ''' Compute the phase field variable at some point. A cosine profile is assumed (see https://mooseframework.inl.gov/source/ics/SmoothCircleIC.html). Input : the radius R of the grain in the direction (a float) the distance r between the current point and the center (a float) the width w of the interface (a float) Output : the value of the phase field variable (a float) ''' #inside the grain if r<R-w/2: return 1 #outside the grain elif r>R+w/2: return 0 #inside the interface else : return 0.5*(1 + math.cos(math.pi*(r-R+w/2)/w)) def Extract_solute_at_p(dict_sample, ij_p)-
Extract the value of the solute concentration at a given point.
Input : a sample dictionnary (a dict) a coordinate of the point (a tuple of int) Output : the value of the solute concentration (a float)
Expand source code
def Extract_solute_at_p(dict_sample,ij_p): ''' Extract the value of the solute concentration at a given point. Input : a sample dictionnary (a dict) a coordinate of the point (a tuple of int) Output : the value of the solute concentration (a float) ''' return dict_sample['solute_M'][-1-ij_p[0]][ij_p[1]] def Interpolate_solute_out_grains(dict_algorithm, dict_sample)-
For all node of the mesh, if there is solute in one grain it is moved to the nearest node out of the grain.
A node in multiple grain (contact area) is not considered in a grain a solute can exist.
Input : an algorithm dictionnary (a dict) a sample dictionnary (a dict) Output : Nothing, but the solute map is updated (a ny x nx numpy array)Expand source code
def Interpolate_solute_out_grains(dict_algorithm, dict_sample) : """ For all node of the mesh, if there is solute in one grain it is moved to the nearest node out of the grain. A node in multiple grain (contact area) is not considered in a grain a solute can exist. Input : an algorithm dictionnary (a dict) a sample dictionnary (a dict) Output : Nothing, but the solute map is updated (a ny x nx numpy array) """ #Create maps of different available nodes node_available_M = np.array(np.ones((len(dict_sample['y_L']),len(dict_sample['x_L'])),dtype = bool)) for l in range(len(dict_sample['y_L'])): for c in range(len(dict_sample['x_L'])): n_grain = 0 for etai in dict_sample['L_etai']: if etai.etai_M[l][c] > 0.5 : n_grain = n_grain + 1 if n_grain == 1: node_available_M[l][c] = False #dilatation dilated_M = binary_dilation(node_available_M, dict_algorithm['struct_element']) #iteration on solute map to see if some solute is in not available position for l in range(len(dict_sample['y_L'])): for c in range(len(dict_sample['x_L'])): if dict_sample['solute_M'][l][c] > 0 and not dilated_M[l][c]: solute_moved = False size_window = 1 while not solute_moved : i_window = 0 while not solute_moved and i_window <= size_window: n_node_available = 0 #Look to move horizontaly and vertically if i_window == 0 : top_available = False down_available = False left_available = False right_available = False #to the top if l - size_window > 0: top_available = dilated_M[l-size_window][c] if dilated_M[l-size_window][c] : n_node_available = n_node_available + 1 #to the down if l + size_window < len(dict_sample['y_L']): down_available = dilated_M[l+size_window][c] if dilated_M[l+size_window][c] : n_node_available = n_node_available + 1 #to the left if c - size_window > 0: left_available = dilated_M[l][c-size_window] if dilated_M[l][c-size_window] : n_node_available = n_node_available + 1 #to the right if c + size_window < len(dict_sample['x_L']): right_available = dilated_M[l][c+size_window] if dilated_M[l][c+size_window] : n_node_available = n_node_available + 1 #move solute if et least one node is available if n_node_available != 0 : #to the top if top_available: dict_sample['solute_M'][l-size_window][c] = dict_sample['solute_M'][l-size_window][c] + dict_sample['solute_M'][l][c]/n_node_available #to the down if down_available: dict_sample['solute_M'][l+size_window][c] = dict_sample['solute_M'][l+size_window][c] + dict_sample['solute_M'][l][c]/n_node_available #to the left if left_available: dict_sample['solute_M'][l][c-size_window] = dict_sample['solute_M'][l][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available #to the right if right_available: dict_sample['solute_M'][l][c+size_window] = dict_sample['solute_M'][l][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available dict_sample['solute_M'][l][c] = 0 solute_moved = True #Look to move diagonally else : top_min_available = False top_max_available = False down_min_available = False down_max_available = False left_min_available = False left_max_available = False right_min_available = False right_max_available = False #to the top if l - size_window > 0: if c - i_window > 0 : top_min_available = dilated_M[l-size_window][c-i_window] if dilated_M[l-size_window][c-i_window] : n_node_available = n_node_available + 1 if c + i_window < len(dict_sample['x_L']): top_max_available = dilated_M[l-size_window][c+i_window] if dilated_M[l-size_window][c+i_window] : n_node_available = n_node_available + 1 #to the down if l + size_window < len(dict_sample['y_L']): if c - i_window > 0 : down_min_available = dilated_M[l+size_window][c-i_window] if dilated_M[l+size_window][c-i_window] : n_node_available = n_node_available + 1 if c + i_window < len(dict_sample['x_L']): down_max_available = dilated_M[l+size_window][c+i_window] if dilated_M[l+size_window][c+i_window] : n_node_available = n_node_available + 1 #to the left if c - size_window > 0: if l - i_window > 0 : left_min_available = dilated_M[l-i_window][c-size_window] if dilated_M[l-i_window][c-size_window] : n_node_available = n_node_available + 1 if l + i_window < len(dict_sample['y_L']): left_max_available = dilated_M[l+i_window][c-size_window] if dilated_M[l+i_window][c-size_window] : n_node_available = n_node_available + 1 #to the right if c + size_window < len(dict_sample['x_L']): if l - i_window > 0 : right_min_available = dilated_M[l-i_window][c+size_window] if dilated_M[l-i_window][c+size_window] : n_node_available = n_node_available + 1 if l + i_window < len(dict_sample['y_L']): right_max_available = dilated_M[l+i_window][c+size_window] if dilated_M[l+i_window][c+size_window] : n_node_available = n_node_available + 1 #move solute if et least one node is available if n_node_available != 0 : #to the top if top_min_available: dict_sample['solute_M'][l-size_window][c-i_window] = dict_sample['solute_M'][l-size_window][c-i_window] + dict_sample['solute_M'][l][c]/n_node_available if top_max_available: dict_sample['solute_M'][l-size_window][c+i_window] = dict_sample['solute_M'][l-size_window][c+i_window] + dict_sample['solute_M'][l][c]/n_node_available #to the down if down_min_available: dict_sample['solute_M'][l+size_window][c-i_window] = dict_sample['solute_M'][l+size_window][c-i_window] + dict_sample['solute_M'][l][c]/n_node_available if down_max_available: dict_sample['solute_M'][l+size_window][c+i_window] = dict_sample['solute_M'][l+size_window][c+i_window] + dict_sample['solute_M'][l][c]/n_node_available #to the left if left_min_available: dict_sample['solute_M'][l-i_window][c-size_window] = dict_sample['solute_M'][l-i_window][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available if left_max_available: dict_sample['solute_M'][l+i_window][c-size_window] = dict_sample['solute_M'][l+i_window][c-size_window] + dict_sample['solute_M'][l][c]/n_node_available #to the right if right_min_available: dict_sample['solute_M'][l-i_window][c+size_window] = dict_sample['solute_M'][l-i_window][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available if right_max_available: dict_sample['solute_M'][l+i_window][c+size_window] = dict_sample['solute_M'][l+i_window][c+size_window] + dict_sample['solute_M'][l][c]/n_node_available dict_sample['solute_M'][l][c] = 0 solute_moved = True i_window = i_window + 1 size_window = size_window + 1 def Reset_y_max(L_g, Force)-
The upper wall is located as a single contact verify the target value.
Input : the list of temporary grains (a list) the confinement force (a float) Output : the upper wall position (a float)Expand source code
def Reset_y_max(L_g,Force): """ The upper wall is located as a single contact verify the target value. Input : the list of temporary grains (a list) the confinement force (a float) Output : the upper wall position (a float) """ print('Reset of the y_max on the upper grain') y_max = None id_grain_max = None for id_grain in range(len(L_g)): grain = L_g[id_grain] y_max_grain = max(grain.l_border_y) if y_max != None and y_max_grain > y_max: y_max = y_max_grain id_grain_max = id_grain elif y_max == None: y_max = y_max_grain id_grain_max = id_grain k = 5*4/3*L_g[id_grain_max].y/(1-L_g[id_grain_max].nu*L_g[id_grain_max].nu)*math.sqrt(L_g[id_grain_max].r_mean) y_max = y_max - (Force/k)**(2/3) return y_max def Sort_Files(dict_algorithm)-
Sort files generated by MOOSE to different directories
Input : an algorithm dictionnary (a dict) Output : Nothing but files are sorted
Expand source code
def Sort_Files(dict_algorithm): ''' Sort files generated by MOOSE to different directories Input : an algorithm dictionnary (a dict) Output : Nothing but files are sorted ''' os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_out.e','Output/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_out.e') os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'.i','Input/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'.i') j = 0 j_str = index_to_str(j) filepath = Path(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu') while filepath.exists(): for i_proc in range(dict_algorithm['np_proc']): os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'_'+str(i_proc)+'.vtu','Output/Ite_'+str(dict_algorithm['i_PFDEM'])+'/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'_'+str(i_proc)+'.vtu') os.rename(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu','Output/Ite_'+str(dict_algorithm['i_PFDEM'])+'/'+dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu') j = j + 1 j_str = index_to_str(j) filepath = Path(dict_algorithm['namefile']+'_'+str(dict_algorithm['i_PFDEM'])+'_other_'+j_str+'.pvtu') return index_to_str(j-1) def error_on_ymax_df(dy, overlap_L, k_L)-
Compute the derivative function df to control the upper wall (error_on_ymax_f()).
Input : an increment of the upper wall position (a float) a list of overlap for contact between grain and upper wall (a list) a list of spring for contact between grain and upper wall (a list) Output : the derivative of error_on_ymax_f() (a float)Expand source code
def error_on_ymax_df(dy,overlap_L,k_L) : """ Compute the derivative function df to control the upper wall (error_on_ymax_f()). Input : an increment of the upper wall position (a float) a list of overlap for contact between grain and upper wall (a list) a list of spring for contact between grain and upper wall (a list) Output : the derivative of error_on_ymax_f() (a float) """ df = 0 for i in range(len(overlap_L)): df = df + 3/2*k_L[i]*(max(overlap_L[i]-dy,0))**(1/2) return df def error_on_ymax_f(dy, overlap_L, k_L, Force_target)-
Compute the function f to control the upper wall. It is the difference between the force applied and the target value.
Input : an increment of the upper wall position (a float) a list of overlap for contact between grain and upper wall (a list) a list of spring for contact between grain and upper wall (a list) a confinement force (a float) Output : the difference between the force applied and the confinement (a float)Expand source code
def error_on_ymax_f(dy,overlap_L,k_L,Force_target) : """ Compute the function f to control the upper wall. It is the difference between the force applied and the target value. Input : an increment of the upper wall position (a float) a list of overlap for contact between grain and upper wall (a list) a list of spring for contact between grain and upper wall (a list) a confinement force (a float) Output : the difference between the force applied and the confinement (a float) """ f = Force_target for i in range(len(overlap_L)): f = f - k_L[i]*(max(overlap_L[i]-dy,0))**(3/2) return f def index_to_str(j)-
An integer is converted to a float with 3 components
Input : a number (a float) Output : a string with 3 components (a string)
Expand source code
def index_to_str(j): ''' An integer is converted to a float with 3 components Input : a number (a float) Output : a string with 3 components (a string) ''' if j < 10: j_str = '00'+str(j) elif 10 <= j and j < 100: j_str = '0'+str(j) else : j_str = str(j) return j_str