Homec4science
Diffusion Additive Manufacturing Work (scratch_cluster_run) Browse

R11910/mechanical_model/transfer_previous_layer.py
2616fcb08aabscratch_cluster_run

Last Change
Raw File

transfer_previous_layer.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(), "muspectre/build/language_bindings/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "muspectre/build/language_bindings/libmufft/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "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
from memory_profiler import profile
import gc
class EigenStrain3:
"""
Class whose eigen_strain_func function is used to apply eigen strains
(gel expansion)
"""
def __init__(self,eigenstrain_shape, thermal_expansion,cell):
self.cell = cell
self.pixels = self.cell.pixels
self.nb_sub_grid = self.cell.nb_subdomain_grid_pts
self.sub_loc = self.cell.subdomain_locations
self.eigenstrain_shape = eigenstrain_shape
self.alpha = thermal_expansion
self.temperature_field = self.cell.get_fields().get_field('temperature').array()
self.material_field = self.cell.get_fields().get_field('material').array()
self.nb_quad_pts = self.cell.nb_quad_pts
def __call__(self, strain_field):
self.eigen_strain_func(strain_field)
def eigen_strain_func(self, strain_field):
# Looping over pixels locally and apply eigen strain based on the
# global boolean field eigen_pixels_in_structure_eigen
for i in range(self.nb_sub_grid[0]):
for j in range(self.nb_sub_grid[1]):
for k in range(self.nb_sub_grid[2]):
for q in range(self.nb_quad_pts):
if self.material_field[0, q, i,j,k]:
strain_field[:, :, q, i, j, k] -= self.alpha*self.eigenstrain_shape*self.temperature_field[0,q,i,j,k]
# def parse_args():
# """
# Function to handle arguments
# """
# parser = argparse.ArgumentParser(description='Run given layer')
# parser.add_argument("layer_no", type=int, help='layer_number')
# args = parser.parse_args()
# return args
@profile
def transfer_previous_layer(layer_no):
# args = parse_args()
# layer_no = args.layer_no
from mpi4py import MPI
comm_py = MPI.COMM_SELF
comm = muGrid.Communicator(comm_py)
#rank = comm.rank
#procs = comm.size
os.chdir('/home/ekinkubilay/Documents/AM/mechanical_model/results')
#get the geometry via fenics
mesh = Mesh(MPI.COMM_SELF)
filename = "/home/ekinkubilay/Documents/AM/Part_geometry/mesh/layer_"+str(layer_no).zfill(3)+".xdmf"
temperature_filename = "/home/ekinkubilay/Documents/AM/Thermal_model/results/temperature/temperature_layer_"+str(layer_no).zfill(3)+".h5"
f = XDMFFile(MPI.COMM_SELF, filename)
f.read(mesh)
MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m
f.close()
bounding_box = mesh.bounding_box_tree()
#nb_steps = 20
dim = 3
nb_quad = 6
thermal_expansion = 1e-5
x_max = np.max(mesh.coordinates()[:,0])
x_min = np.min(mesh.coordinates()[:,0])
y_max = np.max(mesh.coordinates()[:,1])
y_min = np.min(mesh.coordinates()[:,1])
z_max = np.max(mesh.coordinates()[:,2])
z_min = np.min(mesh.coordinates()[:,2])
xmax = np.array([np.NAN])
ymax = np.array([np.NAN])
zmax = np.array([np.NAN])
xmin = np.array([np.NAN])
ymin = np.array([np.NAN])
zmin = np.array([np.NAN])
dx, dy, dz = 0.00025, 0.00025, 0.00003
#z_lim = 0.00216
element_dimensions = [dx, dy, dz]
#nb_domain_grid_pts = [11,11,11]
domain_lengths = [x_max-x_min,y_max-y_min,z_max-z_min]
nb_domain_grid_pts = list(np.array(domain_lengths)/np.array(element_dimensions)+1)
nb_domain_grid_pts = [int(i) for i in nb_domain_grid_pts]
#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))
Youngs_modulus = 10e9
Poisson_ratio = 1/3
## 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
## macroscopic strain
strain_step = np.zeros((3,3))
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")
bulk = np.zeros(nb_domain_grid_pts)
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, z_max, nb_domain_grid_pts[2])
#assign material
for i,x in enumerate(x_dim):
for j,y in enumerate(y_dim):
for k,z in enumerate(z_dim):
if len(bounding_box.compute_collisions(Point(x+0.5*dx,y+0.5*dy,z+0.5*dz))) != 0:
bulk[i,j,k] = 1
del(bounding_box)
gc.collect()
#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.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
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)
material = msp.material.MaterialLinearElastic1_3d.make(cell, "material", Youngs_modulus , Poisson_ratio)
vacuum = msp.material.MaterialLinearElastic1_3d.make(cell, "vacuum", 0.0, Poisson_ratio)
base = msp.material.MaterialLinearElastic1_3d.make(cell, "base", 2*Youngs_modulus, Poisson_ratio)
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)
eigenstrain_shape = np.identity(3)
gradient = Stencils3D.linear_finite_elements
weights = 6 * [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)
cell.get_fields().set_nb_sub_pts("quad", nb_quad)
material_phase = cell.get_fields().register_int_field("material", 1, 'quad')
coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
temperature = cell.get_fields().register_real_field("temperature", 1 , 'quad')
vacuum_strain_field = cell.get_fields().register_real_field("vacuum_strain", (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()
# 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=["coordinates", "flux",
"grad",
"incrF",
"material",
"rhs",
"tangent", "temperature", "vacuum_strain"])
eigen_class = EigenStrain3(eigenstrain_shape, thermal_expansion, cell)
#mat_fields = material.list_fields()
#file_io_object_read.register_field_collection(field_collection=material.collection)
# sub_domain_loc = cell.subdomain_locations
# #material_phase = cell.get_fields().register_int_field("material", 1, 'quad')
# material_phase_alias = np.array(material_phase, copy=False)
# coordinates_alias = np.array(coordinates, copy=False)
# for pixel_id, pixel_coord in cell.pixels.enumerate():
# i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2]
# pix_mat = int(bulk[i, j, k])
# for q in range(nb_quad):
# # local coordinates on the processor
# a = i - sub_domain_loc[0]
# b = j - sub_domain_loc[1]
# c = k - sub_domain_loc[2]
# material_phase_alias[a, b, c] = pix_mat
# coordinates_alias[0, 0, a, b, c] = i
# coordinates_alias[1, 0, a, b, c] = j
# coordinates_alias[2, 0, a, b, c] = k
timestamp_arr = file_io_object_read.read_global_attribute('timestamp')
nb_steps = len(timestamp_arr)
nb_steps = 5
start = timer.time()
import meshio
points = msp.linear_finite_elements.get_position_3d_helper(cell, solver=solver)
nx, ny, nz = cell.nb_domain_grid_pts
xc, yc, zc = np.mgrid[:nx, :ny, :nz]
# Global node indices
def c2i(xp, yp, zp):
return xp + (nx+1) * (yp + (ny+1) * zp)
cells = cells = np.swapaxes(
[[c2i(xc, yc, zc), c2i(xc+1, yc, zc),
c2i(xc+1, yc+1, zc), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc+1, yc, zc),
c2i(xc+1, yc, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc+1, zc),
c2i(xc+1, yc+1, zc), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc+1, zc),
c2i(xc, yc+1, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc, zc+1),
c2i(xc+1, yc, zc+1), c2i(xc+1, yc+1, zc+1)],
[c2i(xc, yc, zc), c2i(xc, yc, zc+1),
c2i(xc, yc+1, zc+1), c2i(xc+1, yc+1, zc+1)]],
0, 1)
cells = cells.reshape((4, -1), order='F').T
# with meshio.xdmf.TimeSeriesWriter('layer_'+str(layer_no).zfill(3)+'.xdmf') as writer:
# writer.write_points_cells(points, {"hexahedron": cells})
file_io_object_read.read(
frame=nb_steps-1,
field_names=["coordinates", "flux",
"grad",
"material",
"temperature"])
file_io_object_read.close()
return cell
# for i in range(nb_steps):
# print('=== STEP {}/{} ==='.format(i+1, nb_steps))
# file_io_object_read.read(
# frame=i,
# field_names=["coordinates", "flux",
# "grad",
# "material",
# "temperature"])
# #strain_arr = solver.grad.field.array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
# #stress_arr = solver.flux.field.array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
# material_arr = cell.get_fields().get_field('material').array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2)
# flux_arr = cell.get_fields().get_field('flux').array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
# grad_arr = cell.get_fields().get_field('grad').array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
# coordinates_arr = cell.get_fields().get_field('coordinates').array().reshape((dim, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, 1, -1), order='F').T.swapaxes(1, 2)
# temperature_arr = cell.get_fields().get_field('temperature').array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2)
# c_data = {
# "material": np.array([material_arr]),
# "temperature": np.array([temperature_arr]),
# "stress": np.array([flux_arr]),
# "strain": np.array([grad_arr]),
# "coordinates": np.array([coordinates_arr])
# }
# #writer.write_data(timestamp_arr[i], cell_data=c_data)
# #msp.linear_finite_elements.write_3d_class(
# #"ekin_tutorial_3D_spectral_output_{}.xdmf".format(i), cell, solver, cell_data=c_data)
# #comm_py.barrier()
# #stress = cell.get_fields().get_field('flux').array()
# #strain = cell.get_fields().get_field('grad').array()
# #temperature = cell.get_fields().get_field('temperature').array()
# #np.savez('previous_fields.npz', stress=stress, strain= strain, temperature=temperature)
# file_io_object_read.close()
# return cell
c4science · Help