Page MenuHomec4science
Files F88373774

netcdf2h5_hea.py
No OneTemporary
Actions

File Metadata

Created
Fri, Oct 18, 11:23

netcdf2h5_hea.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 2 13:47:04 2022
@author: ekinkubilay
"""
#from netCDF4 import Dataset
"""
nc_file_name = 'spectral_3D_output.nc'
data = Dataset(nc_file_name, 'r')
print(data)
grad = data.variables['grad'][:,:,:,:,:,:]
nx = data.dimensions['nx']
print(type(grad))
hf = h5py.File('data.h5', 'w')
hf.create_dataset('dataset_1', data=grad)
hf.close()
"""
import os
import sys
import h5py
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/libmufft/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/libmugrid/python"))
sys.path.insert(0, os.path.join(os.getcwd(), "/home/kubilay/software/muspectre/build/language_bindings/python/muSpectre"))
print("current working directory: '{}'".format(os.getcwd()))
import muFFT
import muSpectre as msp
import numpy as np
import matplotlib.pyplot as plt
import muGrid
from muFFT import Stencils2D
from muFFT import Stencils3D
import cProfile, pstats
import time as timer
import fenics
from petsc4py import PETSc
import argparse
from muSpectre.gradient_integration import get_complemented_positions_class_solver
def parse_args():
"""
Function to handle arguments
"""
parser = argparse.ArgumentParser(description='Run given layer')
parser.add_argument("layer_no", type=int, help='layer_number')
parser.add_argument("factor", type=float, help='factor for material parametric analysis')
parser.add_argument("target_directory", help='target directory for output/input files')
args = parser.parse_args()
return args
def read_data_of_step(file_io_object_read, cell, solver, read_netcdf, material, read_frame=-1):
nb_quad = cell.nb_quad_pts
nb_domain_grid_pts = cell.nb_domain_grid_pts
dim = cell.dims
alpha_thermal = material.alpha_thermal
T_m = material.T_m
epsilon_sat = 0.35
sigma_for = 1.200 #GPa
if read_netcdf:
file_io_object_read.read(
frame=read_frame,
field_names=["flux", "grad", "epsilon_bar", "epsilon_plastic", "temperature", "vacuum_strain", "melt_pool", "material2"])
material_arr = cell.get_fields().get_field('material2').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)
flux_arr_2 = np.array(flux_arr, copy=True)
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)
epsilon_plastic = cell.get_globalised_current_real_field("epsilon_plastic").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
epsilon_bar = cell.get_globalised_current_real_field("epsilon_bar").array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2)
vacuum_strain = cell.get_globalised_internal_real_field("vacuum_strain").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, dim, -1), order='F').T.swapaxes(1, 2)
temperature_material = cell.get_globalised_internal_real_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)
melt_pool = cell.get_globalised_internal_int_field("melt_pool").array().reshape((1, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((1, 1, -1), order='F').T.swapaxes(1, 2)
thermal_strain = np.zeros_like(grad_arr)
thermal_strain_value = (temperature_material[:,0,0] - T_m)*alpha_thermal
for i in range(3):
thermal_strain[:,i,i] = np.array([thermal_strain_value])
elastic_strain = grad_arr - epsilon_plastic - vacuum_strain - thermal_strain
total_actual_strain = grad_arr - vacuum_strain
sigma_forest = np.zeros_like(temperature_material)
sigma_forest = sigma_for*(1-np.exp(-epsilon_bar[:,0,0]/epsilon_sat))*(1- temperature_material[:,0,0]/T_m)**0.5
stress_deviatoric = np.zeros_like(flux_arr)
trace_array = (flux_arr[:,0,0] + flux_arr[:,1,1] + flux_arr[:,2,2])/3
stress_deviatoric += flux_arr
for i in range(3):
stress_deviatoric[:,i,i] -= trace_array
sigma_vonmises = np.sqrt((3/2)*np.sum(stress_deviatoric*stress_deviatoric, axis=(1,2)))
#flux_integral = np.repeat(np.expand_dims(solver.projection.integrate(solver.flux.field).array(), 2), 6, axis=2).reshape((dim, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, 1, -1), order='F').T.swapaxes(1, 2)
#flux_integral_nonaffine = np.repeat(np.expand_dims(solver.projection.integrate_nonaffine_displacements(solver.flux.field).array(), 2), 6, axis=2).reshape((dim, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, 1, -1), order='F').T.swapaxes(1, 2)
grad_integral = np.repeat(np.expand_dims(solver.projection.integrate(solver.grad.field).array(), 2), 6, axis=2).reshape((dim, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, 1, -1), order='F').T.swapaxes(1, 2)
grad_integral_nonaffine = np.repeat(np.expand_dims(solver.projection.integrate_nonaffine_displacements(solver.grad.field).array(), 2), 6, axis=2).reshape((dim, 1, nb_quad) + tuple(nb_domain_grid_pts),order='f').reshape((dim, 1, -1), order='F').T.swapaxes(1, 2)
solver.projection.apply_projection(solver.flux.field)
c_data = {
"material": np.array([material_arr]),
"stress": np.array([flux_arr_2]),
"strain": np.array([grad_arr]),
"epsilon_plastic": np.array([epsilon_plastic]),
"epsilon_bar": np.array([epsilon_bar]),
"vacuum_strain": np.array([vacuum_strain]),
"temperature_material": np.array([temperature_material]),
"melt_pool": np.array([melt_pool]),
"epsilon_elastic": np.array([elastic_strain]),
"sigma_forest": np.array([sigma_forest]),
"sigma_vonmises": np.array([sigma_vonmises]),
"stress_deviatoric": np.array([stress_deviatoric]),
"total_actual_strain": np.array([total_actual_strain]),
"projection": np.array([flux_arr]),
#"flux_integral": np.array([flux_integral]),
"grad_integral": np.array([grad_integral]),
#"flux_integral_nonaffine": np.array([flux_integral_nonaffine]),
"grad_integral_nonaffine": np.array([grad_integral_nonaffine])
}
strain_field = cell.get_fields().get_field("grad")
material_field = cell.get_fields().get_field("material2")
strain_field_alias = np.array(strain_field, copy=False)
strain_field_alias[material_field == 0] = 0
#strain_field_alias[:,:,:,-1,:,:] = 0
#strain_field_alias[:,:,:,:,-1,:] = 0
#strain_field_alias[:,:,:,:,:,-1] = 0
#strain_field_alias[:,:,:,:,:,0] = 0
#strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5))
[x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True)
#x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
total_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')])
#file_io_object_read.read(frame=read_frame,field_names=["grad"])
strain_field = cell.get_fields().get_field("grad")
strain_field_alias = np.array(strain_field, copy=False)
vacuum_strain_field = cell.get_globalised_internal_real_field("vacuum_strain").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f')
strain_field_alias -= vacuum_strain_field
#strain_field_alias[:,:,:,-1,:,:] = 0
#strain_field_alias[:,:,:,:,-1,:] = 0
#strain_field_alias[:,:,:,:,:,-1] = 0
#strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5))
[x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True)
#x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
actual_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')])
strain_field = cell.get_fields().get_field("grad")
strain_field_alias = np.array(strain_field, copy=False)
#strain_field_alias[:,:,:,-1,:,:] = 0
#strain_field_alias[:,:,:,:,-1,:] = 0
#strain_field_alias[:,:,:,:,:,-1] = 0
#strain_field_alias[:,:,:,:,:,:] = np.roll(strain_field_alias, (nb_domain_grid_pts[0]//2, nb_domain_grid_pts[1]//2, nb_domain_grid_pts[2]//2), (3,4,5))
strain_field_alias[:,:,:,:,:] = cell.get_globalised_current_real_field("epsilon_plastic").array().reshape((dim, dim, nb_quad) + tuple(nb_domain_grid_pts),order='f')[:,:,:,:,:]
[x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver("0n", cell, solver, np.eye(3), periodically_complemented=True)
#x_displ = np.roll(x_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#y_displ = np.roll(y_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
#z_displ = np.roll(z_displ, (-(nb_domain_grid_pts[0]+1)//2, -(nb_domain_grid_pts[1]+1)//2, -(nb_domain_grid_pts[2]+1)//2), (0,1,2))
plastic_displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')])
point_data = {'actual_displacement' : actual_displacement, 'total_displacement': total_displacement, 'plastic_displacement': plastic_displacement}
file_io_object_read.read(frame=-1,field_names=["grad", "material2"])
return c_data, point_data
def main():
args = parse_args()
layer_no = args.layer_no
factor = args.factor
target_directory = args.target_directory
from mpi4py import MPI
comm_py = MPI.COMM_WORLD
comm = muGrid.Communicator(comm_py)
#rank = comm.rank
#procs = comm.size
#os.chdir('/scratch/kubilay/mechanical_model/results/mu_layer')
os.chdir(target_directory)
#get the geometry via fenics
mesh = fenics.Mesh()
filename = "/work/lammm/kubilay/mesh/layer_"+str(layer_no).zfill(3)+".xdmf"
temperature_filename = "/work/lammm/kubilay/results/temperature/temperature_layer_"+str(layer_no).zfill(3)+".h5"
f = fenics.XDMFFile(filename)
f.read(mesh)
fenics.MeshTransformation.rescale(mesh, 1e-3, fenics.Point(0,0,0)) #in m
f.close()
bounding_box = mesh.bounding_box_tree()
#nb_steps = 20
dim = 3
nb_quad = 6
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.00045, 0.0005, 0.00003
z_lim = 0.00801
#x_lim = 0.05496
#x_lim_min = -0.00004
x_lim = 0.05505
x_lim_min = 0.00005
element_dimensions = [dx, dy, dz]
#nb_domain_grid_pts = [11,11,11]
domain_lengths = [x_lim-x_lim_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]
#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*factor**(2/3) #eV
epsilon_sat = 0.35 # unitless
sigma_gb = 0.150 #GPa
sigma_s0 = 0.254*factor**(4/3) #GPa
sigma_f = 1.200 #GPa
## define the convergence tolerance for the Newton-Raphson increment
newton_tol = 1e-12
equil_tol = newton_tol
## tolerance for the solver of the linear cell
cg_tol = -1e-2
maxiter = 10000 # 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.ones(nb_domain_grid_pts)
#x_dim = np.linspace(x_min, x_lim, nb_domain_grid_pts[0])
#y_dim = np.linspace(y_min, y_max, nb_domain_grid_pts[1])
#z_dim = np.linspace(z_lim, 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(fenics.Point(x+0.5*dx,y+0.5*dy,z+0.5*dz))) != 0:
# bulk[i,j,k] = 1
#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[:,:,0] *= 2
#add base plate at the bottom and vacuum on top and sides
#increase the number of grid points and total domain length in each dimension
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, 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, 0, 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, 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]
bulk = np.insert(bulk, 0, 2, axis=2)
nb_domain_grid_pts[2] += 1
domain_lengths[2] += element_dimensions[2]
#domain discretisation TO DO:muspectre has this array already
x_dim = np.linspace(x_lim_min-element_dimensions[0], x_lim, nb_domain_grid_pts[0])+0.5*dx
y_dim = np.linspace(ymin-element_dimensions[1], ymax, nb_domain_grid_pts[1])+0.5*dy
z_dim = np.linspace(z_lim-4*element_dimensions[2], zmax, nb_domain_grid_pts[2])+0.5*dz
coords = np.zeros([nb_domain_grid_pts[0],nb_domain_grid_pts[1],nb_domain_grid_pts[2],3])
temperature_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts))
bulk_array = np.zeros_like(temperature_array)
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]
#s_list = [(0.005, 0.065)]
s_list = [(i,i+0.001) for i in np.linspace(0.005,0.050, 26)]
support_mask = np.ones_like(x_dim)
for i,x in enumerate(x_dim):
loc = x
for j,s in enumerate(s_list):
if loc>s[0] and loc<s[1]:
support_mask[i] = 0
#bulk[:,:,0] = 2
#bulk[:,-1,1] = 0
#bulk[-1,:,1] = 0
#for j,y in enumerate(y_dim):
# bulk[:,j,0] = bulk[:,j,0]*support_mask
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)
bulk_array[0,:,pixel[0], pixel[1], pixel[2]] = 1
elif bulk[tuple(pixel)]==0:
vacuum.add_pixel(pixel_index)
elif bulk[tuple(pixel)]==2:
base.add_pixel(pixel_index)
gradient = Stencils3D.linear_finite_elements_6_regular
weights = nb_quad*[1]
material.T_m = 1600
material.alpha_thermal = 22e-6
material.identify_temperature_loaded_quad_pts(cell, bulk_array.reshape(-1,1, order="F"))
## use solver
krylov_solver = msp.solvers.KrylovSolverCG(cg_tol, maxiter, msp.Verbosity.Full)
solver = msp.solvers.SolverNewtonCG(cell, krylov_solver, msp.Verbosity.Full, newton_tol, equil_tol, maxiter, gradient, weights, msp.solvers.MeanControl.stress_control)
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("material2", 1, 'quad')
material_phase_alias = np.array(material_phase, copy=False)
#coordinates = cell.get_fields().register_real_field("coordinates", 3, 'quad')
#temperature = cell.get_fields().register_real_field("temperature2", 1 , 'quad')
#vacuum_strain_field = cell.get_fields().register_real_field("vacuum_strain2", (3,3), "quad")
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()
sub_domain_loc = cell.subdomain_locations
# 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",
# "material2",
# "rhs",
# "tangent", "temperature2", "vacuum_strain2"])
#file_io_object_read.register_field_collection(field_collection=material.collection, field_names=["epsilon_bar"])
#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')
#print(timestamp_arr)
nb_steps = len(timestamp_arr)
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
T_m = 1600
alpha_thermal = 22e-6
epsilon_sat = 0.35
sigma_for = 1.200
with meshio.xdmf.TimeSeriesWriter('mu_layer_'+str(layer_no).zfill(3)+'.xdmf') as writer:
writer.write_points_cells(points, {"tetra": cells})
c_data, point_data = read_data_of_step(file_io_object_read, cell, solver, True, material, read_frame=0)
writer.write_data(timestamp_arr[0], cell_data=c_data, point_data=point_data)
if layer_no == 333:
c_data, point_data = read_data_of_step(file_io_object_read, cell, solver, True, material, read_frame=1)
writer.write_data(timestamp_arr[1], cell_data=c_data, point_data=point_data)
material.dt = 3000
res = solver.solve_load_increment(np.zeros((3,3)))
c_data, point_data = read_data_of_step(file_io_object_read, cell, solver, False, material)
writer.write_data(timestamp_arr[1]+material.dt, cell_data=c_data, point_data=point_data)
#for i in range(1):
# print('=== STEP {}/{} ==='.format(i+1, nb_steps))
# if layer_no == 333:
#flux_arr = cell.get_fields().get_field('flux').array()
#grad_arr = cell.get_fields().get_field('grad').array()
#vacuum_strain_arr = cell.get_globalised_internal_real_field("vacuum_strain").array()
#flux_arr_alias = np.array(flux_arr, copy=False)
#grad_arr_alias = np.array(grad_arr, copy=False)
#vacuum_strain_arr_alias = np.array(vacuum_strain_arr, copy=False)
#flux_arr_alias[:,:,:,:,:] = 0
#grad_arr_alias[:,:,:,:,:] = 0
#vacuum_strain_arr_alias[:,:,:,:,:] = 0
#material.update_vacuum_strain_field(cell, vacuum_strain_arr_alias.reshape(-1,1, order="F"))
# material.update_temperature_field(cell, 293*np.ones_like(bulk_array.reshape(-1,1,order="F")))
# material.dt = 3000
#material.alpha_thermal = 0.0
# res = solver.solve_load_increment(np.zeros((3,3)))
#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()
if __name__=='__main__':
main()

Event Timeline

c4science · Help