Homec4science
Diffusion Additive Manufacturing Work (master) Browse

R11910/mechanical_model/transfer_previous_layer_1q.py
3e9a1d30b01fmaster

Last Change
Raw File

transfer_previous_layer_1q.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 2 13:47:04 2022
@author: ekinkubilay
"""
import os
import sys
#import h5py
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/libmufft/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/muspectre/build/language_bindings/libmugrid/python"))
print("current working directory: '{}'".format(os.getcwd()))
#import muFFT
import muSpectre as msp
import numpy as np
import muGrid
#from muFFT import Stencils2D
from muFFT import Stencils3D
#import cProfile, pstats
#import time as timer
#from fenics import Mesh, XDMFFile, MeshTransformation, Point
#from petsc4py import PETSc
#import argparse
#import gc
def transfer_previous_layer(target_directory, layer_no, bulk, nb_domain_grid_pts, domain_lengths):
#, x_max, y_max, z_max, x_min, y_min, z_min):
from mpi4py import MPI
comm_py = MPI.COMM_SELF
comm = muGrid.Communicator(comm_py)
#os.chdir('/scratch/kubilay/mechanical_model/results/mu_layer')
os.chdir(target_directory)
#nb_steps = 20
dim = 3
nb_quad = 1
dx, dy, dz = 0.0002, 0.0003125, 0.00003
# x_max -= 0.5*dx
# y_max -= 0.5*dy
# z_max -= 0.5*dz
# x_min -= 0.5*dx
# y_min -= 0.5*dy
#y_min = 0 # i don't know why this is necessary
# z_min -= 0.5*dz
# z_lim = 0.00801
#x_lim = 0.05496
#x_lim_min = -0.00004
# x_lim = 0.055
#x_lim_min = 0.00005
# z_max -= dz
# z_max -= dz
# y_max -= dy
element_dimensions = [dx, dy, dz]
#nb_domain_grid_pts = [11,11,11]
# domain_lengths = [x_lim-x_min,y_max-y_min,z_max-z_lim]
# nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions))
# nb_domain_grid_pts = [round(i) for i in nb_domain_grid_pts]
# print('transfer', nb_domain_grid_pts, domain_lengths)
#domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min]
#domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min]
#element_dimensions = np.array(domain_lengths)/np.array(nb_domain_grid_pts-np.ones(dim))
#define constant material properties
young = 162 #GPa
poisson = 0.277 #unitless
alpha_s = 0.55 #unitless
Delta_E_bs = 1.13 #eV
epsilon_sat = 0.35 # unitless
sigma_gb = 0.150 #GPa
sigma_s0 = 0.254 #GPa
sigma_f = 1.200 #GPa
## define the convergence tolerance for the Newton-Raphson increment
newton_tol = 1e-6
equil_tol = newton_tol
## tolerance for the solver of the linear cell
cg_tol = 1e-6
maxiter = 1000 # for linear cell solver
#bulk = np.load("bulk.npy")
#coords = np.load("coords.npy")
#nb_domain_grid_pts = np.load("nb_domain_grid_pts.npy")
#domain_lengths = np.load("domain_lengths.npy")
#add base plate at the bottom and vacuum on top and sides
# bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2)
# nb_domain_grid_pts[2] += 1
# domain_lengths[2] += element_dimensions[2]
# bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1)
# nb_domain_grid_pts[1] += 1
# domain_lengths[1] += element_dimensions[1]
# bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0)
# nb_domain_grid_pts[0] += 1
# domain_lengths[0] += element_dimensions[0]
"""
bulk = np.ones(nb_domain_grid_pts)
#vacuum on top layer (2)
bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
bulk = np.insert(bulk, nb_domain_grid_pts[2], 0, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
#vacuum on y plane (2)
bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1)
nb_domain_grid_pts[1] += 1
domain_lengths[1] += element_dimensions[1]
#bulk = np.insert(bulk, nb_domain_grid_pts[1], 0, axis=1)
#nb_domain_grid_pts[1] += 1
#domain_lengths[1] += element_dimensions[1]
#vacuum on x plane (2)
bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0)
nb_domain_grid_pts[0] += 1
domain_lengths[0] += element_dimensions[0]
bulk = np.insert(bulk, nb_domain_grid_pts[0], 0, axis=0)
nb_domain_grid_pts[0] += 1
domain_lengths[0] += element_dimensions[0]
#base layer on bottom layer (4)
bulk = np.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
bulk = np.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
bulk = np.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
bulk = np.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
"""
bulk = np.delete(bulk, -3, 2)
#print('bulk after', np.shape(bulk), np.shape(bulk_array))
#print('domain lengths compare', domain_lengths, domain_lengths_2)
#print('grid pts compare', nb_domain_grid_pts, nb_domain_grid_pts_2)
nb_domain_grid_pts = np.shape(bulk)
domain_lengths = nb_domain_grid_pts*np.array([dx,dy,dz])
# print(np.shape(bulk), nb_domain_grid_pts)
#print('new values', nb_domain_grid_pts_new, domain_lengths_new)
#x_dim = np.linspace(x_min, x_max, nb_domain_grid_pts[0])
#y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1])
#z_dim = np.linspace(z_min-element_dimensions[2], z_max, nb_domain_grid_pts[2])
#coords = np.zeros([nb_domain_grid_pts[0],nb_domain_grid_pts[1],nb_domain_grid_pts[2],3])
#for i in range(nb_domain_grid_pts[0]):
# coords[i,:,:,0] = x_dim[i]
#for i in range(nb_domain_grid_pts[1]):
# coords[:,i,:,1] = y_dim[i]
#for i in range(nb_domain_grid_pts[2]):
# coords[:,:,i,2] = z_dim[i]
cell = msp.cell.CellData.make(list(nb_domain_grid_pts), list(domain_lengths))
cell.nb_quad_pts = nb_quad
cell.get_fields().set_nb_sub_pts("quad", nb_quad)
#define materials
material = msp.material.MaterialHEA_3d.make(cell, "material", young, poisson, alpha_s, Delta_E_bs, epsilon_sat, sigma_gb, sigma_s0, sigma_f)
vacuum = msp.material.MaterialLinearElastic1_3d.make(cell, "vacuum", 0.0, poisson)
base = msp.material.MaterialLinearElastic1_3d.make(cell, "base", 20*young, poisson)
for pixel_index, pixel in enumerate(cell.pixels):
if bulk[tuple(pixel)]==1:
material.add_pixel(pixel_index)
elif bulk[tuple(pixel)]==0:
vacuum.add_pixel(pixel_index)
elif bulk[tuple(pixel)]==2:
base.add_pixel(pixel_index)
gradient = Stencils3D.trilinear_1q_finite_element
weights = nb_quad * [1]
## use solver
krylov_solver = msp.solvers.KrylovSolverCG(cg_tol, maxiter)
solver = msp.solvers.SolverNewtonCG(cell, krylov_solver, msp.Verbosity.Full, newton_tol, equil_tol, maxiter, gradient, weights)
solver.formulation = msp.Formulation.small_strain
solver.initialise_cell()
output_filename = "layer_"+str(layer_no).zfill(3)+".nc"
file_io_object_read = muGrid.FileIONetCDF(output_filename,
muGrid.FileIONetCDF.OpenMode.Read,
comm)
#register fields and field collection
#material_phase = cell.get_fields().register_int_field("material2", 1, 'quad')
#coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
#temperature_field = cell.get_fields().register_real_field("temperature2", 1, 'quad')
#vacuum_strain_field = cell.get_fields().register_real_field("vacuum_strain2", (3,3), "quad")
# register global fields of the cell which you want to write
# print("cell field names:\n", cell.get_field_names())
#cell_field_names = cell.get_field_names()
file_io_object_read.register_field_collection(field_collection=cell.get_fields(),field_names=cell.get_field_names())
file_io_object_read.register_field_collection(field_collection=material.collection)
# register global fields of the cell which you want to write
#print("cell field names:\n", cell.get_field_names())
#cell_field_names = cell.get_field_names()
#file_io_object_read.register_field_collection(field_collection=material.collection)
file_io_object_read.read(
frame=-1,
field_names=["flux","grad", "epsilon_bar", "epsilon_plastic","temperature", "vacuum_strain", "melt_pool"])
material.save_history_variables()
previous_stress = cell.get_fields().get_field('flux').array()
previous_strain = cell.get_fields().get_field('grad').array()
#previous_temperature = cell.get_fields().get_field('temperature2').array()
#previous_vacuum_strain = cell.get_fields().get_field('vacuum_strain2').array()
previous_epsilon_bar = np.insert(cell.get_globalised_current_real_field('epsilon_bar').array(), -1, 0, axis=-1)
#previous_epsilon_bar = previous_epsilon_bar[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
previous_melt_pool = np.insert(cell.get_globalised_internal_int_field('melt_pool').array(), -1, 0, axis=-1)
previous_melt_pool[:,:,-2] = -1
#previous_melt_pool = previous_melt_pool[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
previous_material_temperature = np.insert(cell.get_globalised_internal_real_field('temperature').array(), -1, 0, axis=-1)
#previous_material_temperature = previous_material_temperature[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
#previous_native_stress = np.insert(cell.get_globalised_internal_real_field('native_stress').array(), -1, 0, axis=-1)
#previous_native_stress = previous_native_stress[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
previous_epsilon_plastic = np.insert(cell.get_globalised_current_real_field('epsilon_plastic').array(), -1, 0, axis=-1)
#previous_epsilon_plastic = previous_epsilon_plastic[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
previous_material_vacuum_strain = np.insert(cell.get_globalised_internal_real_field('vacuum_strain').array(), -1, 0, axis=-1)
previous_material_vacuum_strain[:,:,:,:,-2] = previous_strain.reshape((dim*dim, nb_quad) + tuple(nb_domain_grid_pts),order='f')[:,:,:,:,-1]
#previous_strain_alias = np.array(previous_strain, copy=False)
#previous_strain_alias -= cell.get_globalised_internal_real_field('vacuum_strain').array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts), order='f')
#previous_material_vacuum_strain = previous_material_vacuum_strain[:,:,sub_dx:sub_dx+sub_gx,sub_dy:sub_dy+sub_gy,sub_dz:sub_dz+sub_gz]
#melt_pool_mask = np.array(previous_melt_pool, copy=True)
#melt_pool_mask[melt_pool_mask != 1] = 0
#print(np.sum(previous_material_temperature), np.sum(previous_melt_pool))
#print(np.sum(previous_melt_pool*(previous_material_temperature-1650))/np.sum(previous_melt_pool), np.sum(previous_melt_pool*previous_material_temperature)/np.sum(previous_melt_pool)-1650)
file_io_object_read.close()
return [previous_stress, previous_strain, previous_epsilon_bar, previous_melt_pool, previous_material_temperature, previous_epsilon_plastic, previous_material_vacuum_strain]
c4science · Help