Page MenuHomec4science
Files F92422900

netcdf2h5.py
No OneTemporary
Actions

File Metadata

Created
Wed, Nov 20, 04:27

netcdf2h5.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(), "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 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
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()
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]):
if self.material_field[i,j,k]:
strain_field[:, :, 0, i, j, k] -= self.alpha*self.eigenstrain_shape*self.temperature_field[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
def main():
args = parse_args()
layer_no = args.layer_no
from mpi4py import MPI
comm_py = MPI.COMM_WORLD
comm = muGrid.Communicator(comm_py)
#rank = comm.rank
#procs = comm.size
os.chdir('results')
#get the geometry via fenics
mesh = fenics.Mesh()
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 = 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 = 1
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.00018
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_lim-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_lim, 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,y,z))) != 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 = 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+element_dimensions[0], nb_domain_grid_pts[0])
y_dim = np.linspace(y_min, y_max+element_dimensions[1], nb_domain_grid_pts[1])
z_dim = np.linspace(z_min-element_dimensions[2], z_lim+element_dimensions[2], 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.trilinear_1q_finite_element
weights = [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')
# 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"])
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 = 20
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 = np.swapaxes([[c2i(xc, yc, zc),c2i(xc+1, yc, zc),c2i(xc+1, yc+1, zc),c2i(xc, yc+1, zc),c2i(xc, yc, zc+1),c2i(xc+1, yc, zc+1),c2i(xc+1, yc+1, zc+1),c2i(xc, yc+1, zc+1)]],0, 1)
cells = cells.reshape((8, -1), order='F').T
with meshio.xdmf.TimeSeriesWriter('layer_'+str(layer_no).zfill(3)+'.xdmf') as writer:
writer.write_points_cells(points, {"hexahedron": cells})
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()
if __name__=='__main__':
main()

Event Timeline

c4science · Help