diff --git a/mechanical_model/mechanical_model.log b/mechanical_model/mechanical_model.log index 4e77700..3fc59cf 100644 --- a/mechanical_model/mechanical_model.log +++ b/mechanical_model/mechanical_model.log @@ -1,52 +1,17 @@ -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -=== STEP 20/3308 === -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -=== STEP 1/20 === -=== STEP 2/20 === -=== STEP 3/20 === -=== STEP 4/20 === -=== STEP 5/20 === -=== STEP 6/20 === -=== STEP 7/20 === -=== STEP 8/20 === -=== STEP 9/20 === -=== STEP 10/20 === -=== STEP 11/20 === -=== STEP 12/20 === -=== STEP 13/20 === -=== STEP 14/20 === -=== STEP 15/20 === -=== STEP 16/20 === -=== STEP 17/20 === -=== STEP 18/20 === -=== STEP 19/20 === -=== STEP 20/20 === -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -=== STEP 20/2221 === -current working directory: '/home/ekinkubilay/Documents/AM/mechanical_model' -=== STEP 1/20 === -=== STEP 2/20 === -=== STEP 3/20 === -=== STEP 4/20 === -=== STEP 5/20 === -=== STEP 6/20 === -=== STEP 7/20 === -=== STEP 8/20 === -=== STEP 9/20 === -=== STEP 10/20 === -=== STEP 11/20 === -=== STEP 12/20 === -=== STEP 13/20 === -=== STEP 14/20 === -=== STEP 15/20 === -=== STEP 16/20 === -=== STEP 17/20 === -=== STEP 18/20 === -=== STEP 19/20 === -=== STEP 20/20 === +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' +current working directory: '/scratch/kubilay/additive-manufacturing-work/additive-manufacturing-work/mechanical_model' diff --git a/mechanical_model/netcdf2h5.py b/mechanical_model/netcdf2h5.py index 08a7ab8..e51e8af 100644 --- a/mechanical_model/netcdf2h5.py +++ b/mechanical_model/netcdf2h5.py @@ -1,396 +1,395 @@ #!/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")) sys.path.insert(0, os.path.join(os.getcwd(), "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 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 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('/home/ekinkubilay/Documents/AM/mechanical_model/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" + filename = "/scratch/kubilay/layer_"+str(layer_no).zfill(3)+".xdmf" + temperature_filename = "/scratch/kubilay/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 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 + dx, dy, dz = 0.00020, 0.00020, 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(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 = 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') # 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 = 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, {"tetra": 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]) } [x_0, y_0, z_0], [x_displ, y_displ, z_displ] = get_complemented_positions_class_solver( "0d", cell, solver, np.eye(3), periodically_complemented=True) points = np.transpose([x_0.ravel(order='F'), y_0.ravel(order='F'), z_0.ravel(order='F')]) displacement = np.transpose([x_displ.ravel(order='F'), y_displ.ravel(order='F'), z_displ.ravel(order='F')]) point_data = {'displacement' : displacement} writer.write_data(timestamp_arr[i], cell_data=c_data, point_data=point_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() diff --git a/mechanical_model/residual_stress_analysis.py b/mechanical_model/residual_stress_analysis.py index df30dd1..0e98182 100644 --- a/mechanical_model/residual_stress_analysis.py +++ b/mechanical_model/residual_stress_analysis.py @@ -1,574 +1,578 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Jul 25 14:38:47 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")) +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 petsc4py import PETSc import argparse 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 update_temperature_field(cell, temperature_array): sub_domain_loc = cell.subdomain_locations #temperature = np.zeros((cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2])) temperature = np.zeros((1,cell.nb_quad_pts, cell.nb_subdomain_grid_pts[0],cell.nb_subdomain_grid_pts[1],cell.nb_subdomain_grid_pts[2])) for pixel_id, pixel_coord in cell.pixels.enumerate(): i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2] a = i - sub_domain_loc[0] b = j - sub_domain_loc[1] c = k - sub_domain_loc[2] for q in range(cell.nb_quad_pts): temperature[0,q,a,b,c] = temperature_array[0,q,i,j,k] return temperature 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("mu_rank", type=int, help='number of processes for muSpectre') args = parser.parse_args() return args def main(): args = parse_args() layer_no = args.layer_no - os.chdir('/home/ekinkubilay/Documents/AM/mechanical_model/results') + os.chdir('/scratch/kubilay/additive-manufacturing-work/mechanical_model/results') #comm_py is WORLD communicator (with even size), it is split into two #set muspectre communicator from mpi4py import MPI comm_py = MPI.COMM_WORLD rank_py = comm_py.rank size_py = comm_py.size #size must be even, color is tags for the split ranks mu_rank = args.mu_rank color = np.sign(int(rank_py/(size_py-mu_rank))) #split communicator comm = comm_py.Split(color, rank_py) rank = comm.rank size = comm.size - #print(rank_py, size_py, rank, size, color) + print(rank_py, size_py, rank, size, color) #define empty variable for broadcast buffer nb_steps = None - dx, dy, dz = 0.00025, 0.00025, 0.00018 + dx, dy, dz = 0.00020, 0.00020, 0.00018 #z_lim = 0.00216 element_dimensions = [dx, dy, dz] #define global simulation parameters dim = 3 nb_quad = 6 thermal_expansion = 1e-5 #nb_domain_grid_pts = [11,11,11] #define constant material properties 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)) #no applied external strain maxiter = 1000 # for linear cell solver if color == 0: #local size and rank on the split communicator, this communicator is called comm rank = comm.rank procs = comm.size other_size = size_py - size start_time = MPI.Wtime() send_ranks = np.array_split(np.arange(0,other_size, 1)+size, size)[rank] #print(send_ranks) - + from dolfin import Mesh, XDMFFile, MeshTransformation, Point, HDF5File, FunctionSpace, Function #get the geometry and temperature via fenics mesh = Mesh(comm) #read the mesh and the temperature function at initial timestep #using Fenics within the FEM communicator - 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" + + filename = "/scratch/kubilay/layer_"+str(layer_no).zfill(3)+".xdmf" + temperature_filename = "/scratch/kubilay/temperature_layer_"+str(layer_no).zfill(3)+".h5" f = XDMFFile(comm, filename) f.read(mesh) MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m f.close() del(f) gc.collect() bounding_box = mesh.bounding_box_tree() hdf5_file = HDF5File(comm, temperature_filename, "r") V = FunctionSpace(mesh, "CG", 1) temperature_field_store = Function(V) no = 0 vec_name = "/Temperature/vector_%d"%no timestamp = hdf5_file.attributes(vec_name)['timestamp'] nb_steps = hdf5_file.attributes('Temperature')['count']-1 hdf5_file.read(temperature_field_store, vec_name) timestamp_arr = [0.0]*nb_steps # timestamp_arr[0] = timestamp #obtain local domain dimensions xmax = np.max(mesh.coordinates()[:,0]) xmin = np.min(mesh.coordinates()[:,0]) ymax = np.max(mesh.coordinates()[:,1]) ymin = np.min(mesh.coordinates()[:,1]) zmax = np.max(mesh.coordinates()[:,2]) zmin = np.min(mesh.coordinates()[:,2]) #communicate max/min domain values within FEM communicator xmax = comm.allreduce(xmax, op=MPI.MAX) ymax = comm.allreduce(ymax, op=MPI.MAX) zmax = comm.allreduce(zmax, op=MPI.MAX) xmin = comm.allreduce(xmin, op=MPI.MIN) ymin = comm.allreduce(ymin, op=MPI.MIN) zmin = comm.allreduce(zmin, op=MPI.MIN) #calculate domain lengths and element dimensions domain_lengths = [xmax-xmin,ymax-ymin,zmax-zmin] 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] #element_dimensions = np.array(domain_lengths)/np.array(nb_domain_grid_pts-np.ones(dim)) #comm_py.send(domain_lengths, rank_py+2, tag=0) #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)) #array to store material (and not vacuum or base plate) bulk = np.zeros(nb_domain_grid_pts) #domain discretisation TO DO:muspectre has this array already x_dim = np.linspace(xmin, xmax, nb_domain_grid_pts[0]) y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1]) z_dim = np.linspace(zmin, zmax, nb_domain_grid_pts[2]) #obtain the material (and not internal vacuum) by observing if the point #is in or out the submesh. This creates comm.size number of a partially #filled array where it is non-zero only on the sub-mesh present at that rank #obtain the temperature by searching every point at every sub-mesh domain #pass if fenics raises an error (due to point not being present on that domain) 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() + #gather the material across all FEM ranks and broadcast to all FEM ranks #bulk and temperature_array are global (size corresponds to complete domain) bulk = comm.allreduce(bulk, op=MPI.SUM) bulk[bulk > 1] = 1 #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, 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] #domain discretisation TO DO:muspectre has this array already x_dim = np.linspace(xmin, xmax, nb_domain_grid_pts[0]) y_dim = np.linspace(ymin, ymax, nb_domain_grid_pts[1]) z_dim = np.linspace(zmin-element_dimensions[2], zmax, nb_domain_grid_pts[2]) 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)) 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] # print("element_dimensions", element_dimensions) # print("element_dimensions", nb_domain_grid_pts) # print(x_dim) #store the temperature on the array with the right index if point is found #this creates comm.size number of temperature_arrays each non-zero only on #the sub-mesh present in that rank #mesh and temperature are divided differently accross communicators!! TO LOOK! # for i,x in enumerate(x_dim): # for j,y in enumerate(y_dim): # for k,z in enumerate(z_dim): # if bulk[i,j,k] == 1: # for q in range(nb_quad): # try: # temperature_array[0,q,i,j,k] = temperature_field_store(Point(x,y,z)) # except RuntimeError: # pass #send global arrays (material/temperature/coordinates) from FEM communicator #to muSpectre communicator for proc_rank in send_ranks: comm_py.send(rank, proc_rank, tag=0) comm_py.send(bulk, proc_rank, tag=1) comm_py.send(nb_domain_grid_pts, proc_rank, tag=3) comm_py.send(domain_lengths, proc_rank, tag=4) comm_py.send(coords, proc_rank, tag=5) comm_py.send(timestamp_arr, proc_rank, tag=6) comm_py.send(temperature_array, proc_rank, tag=7) end_time = MPI.Wtime() print("FENICS INIT \n Local rank = {}, took {} seconds".format(rank, end_time-start_time)) #for muSepctre processes if color == 1: #define muSpectre communicator, this communicator is called comm_mu comm_mu = muGrid.Communicator(comm) #local size and rank on the split communicator rank = comm.rank procs = comm.size start_time = MPI.Wtime() #recieve from FEM communicator material arrays and system size parameters recieve_rank = comm_py.recv(source = MPI.ANY_SOURCE, tag=0) bulk = comm_py.recv(source = recieve_rank, tag=1) nb_domain_grid_pts = comm_py.recv(source = recieve_rank, tag=3) domain_lengths = comm_py.recv(source = recieve_rank, tag=4) coords = comm_py.recv(source = recieve_rank, tag=5) #recieve empty timestamp array from FEM communicator timestamp_arr = comm_py.recv(source = recieve_rank, tag=6) temperature_array = comm_py.recv(source= recieve_rank, tag=7) #create cell using muSpectre communicator cell = msp.cell.CellData.make_parallel(list(nb_domain_grid_pts), list(domain_lengths), comm_mu) cell.nb_quad_pts = nb_quad #define materials 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) #assign materials 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) #define solver parameters 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.Silent, newton_tol, equil_tol, maxiter, gradient, weights) solver.formulation = msp.Formulation.small_strain solver.initialise_cell() #define file input/output parameters for muSpectre output # output_filename = "layer_"+str(layer_no-1).zfill(3)+".nc" # if comm_mu.rank == 0 and os.path.exists(output_filename): # os.remove(output_filename) output_filename = "layer_"+str(layer_no).zfill(3)+".nc" if comm_mu.rank == 0 and os.path.exists(output_filename): os.remove(output_filename) + + comm_mu.barrier() file_io_object_write = muGrid.FileIONetCDF(output_filename, muGrid.FileIONetCDF.OpenMode.Write, comm_mu) #register fields and field collection 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_field = cell.get_fields().register_real_field("temperature", 1, 'quad') file_io_object_write.register_field_collection(field_collection=cell.get_fields(),field_names=cell.get_field_names()) #initialize global time attribute (to be updated once filled) file_io_object_write.write_global_attribute('timestamp', timestamp_arr) sub_domain_loc = cell.subdomain_locations #create pointer-like copies for fields material_phase_alias = np.array(material_phase, copy=False) coordinates_alias = np.array(coordinates, copy=False) temperature_alias = np.array(temperature_field, copy=False) #assign values to the above created aliases for pixel_id, pixel_coord in cell.pixels.enumerate(): i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2] 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[0, q, a, b, c] = int(bulk[i, j, k]) coordinates_alias[0, q, a, b, c] = coords[i,j,k][0] coordinates_alias[1, q, a, b, c] = coords[i,j,k][1] coordinates_alias[2, q, a, b, c] = coords[i,j,k][2] #use the temperature array to change the values using temperature #array pointer (alias) temperature_alias[:,:,:,:,:] = update_temperature_field(cell, temperature_array)[:,:,:,:,:] if layer_no != 60: strain_field = cell.get_fields().get_field("grad") strain_field_alias = np.array(strain_field, copy=False) stress_field = cell.get_fields().get_field("flux") stress_field_alias = np.array(stress_field, copy=False) - if rank == 0: from transfer_previous_layer import transfer_previous_layer previous_cell = transfer_previous_layer(layer_no-1) previous_stress = previous_cell.get_fields().get_field('flux').array() previous_strain = previous_cell.get_fields().get_field('grad').array() previous_temperature = previous_cell.get_fields().get_field('temperature').array() #npzfile = np.load("previous_fields.npz") #previous_stress = npzfile["stress"] #previous_strain = npzfile["strain"] #previous_temperature = npzfile["temperature"] #if comm_mu.rank == 0 : # os.remove("previous_fields.npz") - + #print(np.shape(strain_field_alias), np.shape(previous_strain)) #assign values to the above created aliases for pixel_id, pixel_coord in cell.pixels.enumerate(): i, j, k = pixel_coord[0], pixel_coord[1], pixel_coord[2] 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] strain_field_alias[:, :, q, a, b, c] = previous_strain[:, :, q, i, j, k] stress_field_alias[:, :, q, a, b, c] = previous_stress[:, :, q, i, j, k] temperature_alias[0,q,a,b,c] = previous_temperature[0,q,i,j,k] comm.bcast(strain_field_alias, root=0) comm.bcast(stress_field_alias, root=0) comm.bcast(temperature_alias, root=0) - + #define eigen_class required by solve_increment eigen_class = EigenStrain3(eigenstrain_shape, thermal_expansion, cell) end_time = MPI.Wtime() print("MUSPECTRE INIT \n Local rank = {}, took {} seconds".format(rank, end_time-start_time)) #make nb_Steps global by broadcasting nb_steps = comm_py.bcast(nb_steps, root = 0) start = timer.time() #seperate operations for FEM and muSpectre communicators depending on ranks if color == 1: #temperature_array = np.zeros(nb_domain_grid_pts) start_time = MPI.Wtime() #for all steps for s in range(nb_steps): print('=== STEP {}/{} ==='.format(s+1, nb_steps)) #recieve temp_array from constant sender request #temperature_array = np.zeros(nb_domain_grid_pts) print("MUSPECTRE Step {} --- Recieve request created: {}".format(s, MPI.Wtime())) req = comm_py.Irecv(temperature_array, source = recieve_rank, tag=s) print("MUSPECTRE Step {} --- Recieve request done: {}".format(s, MPI.Wtime())) #print(np.max(cell.get_fields().get_field("grad").array())) #solve machanics problem res = solver.solve_load_increment(strain_step, eigen_class.eigen_strain_func) comm_mu.barrier() #write all fields in .nc format TO DO:not all fields file_io_object_write.append_frame().write(["coordinates", "flux", "grad", "material", "temperature"]) #make sure temperature array is recieved req.wait() print("MUSPECTRE Step {} --- Recieve request Wait done: {}".format(s, MPI.Wtime())) #use the temperature array to change the values using temperature #array pointer (alias) temperature_alias[:,:,:] = update_temperature_field(cell, temperature_array)[:,:,:] if (s+1)%20 == 0: print('=== STEP {}/{} ==='.format(s+1, nb_steps)) break print("MUSPECTRE Step {} --- Step done: {}".format(s, MPI.Wtime())) end_time = MPI.Wtime() print("MUSPECTRE Step {} --- RUN \n Local rank = {}, took {} seconds".format(s, rank, end_time-start_time)) if color == 0: start_time = MPI.Wtime() quad_coords = [np.array([0.75,0.5,0.25]), np.array([0.75,0.25,0.5]), np.array([0.5,0.75,0.25]), np.array([0.25,0.75,0.5]), np.array([0.5,0.25,0.75]), np.array([0.25,0.5,0.75])] quad_coords = quad_coords*np.array([dx,dy,dz]) for s in range(nb_steps): #start constant sending request print("FENICS Step {} --- Send Requests created: {}".format(s, MPI.Wtime())) reqs = [None] * len(send_ranks) for p, proc_rank in enumerate(send_ranks): reqs[p] = comm_py.Isend([temperature_array, MPI.DOUBLE], proc_rank, tag=s) print("FENICS Step {} --- Send Requests done: {}".format(s, MPI.Wtime())) #re-read the temperature for the next step no += 1 vec_name = "/Temperature/vector_%d"%no timestamp = hdf5_file.attributes(vec_name)['timestamp'] timestamp_arr[s] = timestamp hdf5_file.read(temperature_field_store, vec_name) #store it on a dummy array since we don't want to change #temperature array before the sending/recieving is done #which will be ensured by MPI.wait dummy_array = np.zeros((1,nb_quad)+tuple(nb_domain_grid_pts)) #store the temperature on the array with the right index if point is found #this creates comm.size number of temperature_arrays each non-zero only on #the sub-mesh present in that rank #mesh and temperature are divided differently accross communicators!! TO LOOK! for i,x in enumerate(x_dim): for j,y in enumerate(y_dim): for k,z in enumerate(z_dim): if bulk[i,j,k] == 1: for q in range(nb_quad): try: #dummy_array[0,q,i,j,k] = temperature_field_store(Point(x,y,z)) #if z > 0.00162: print(x,y,z,x+quad_coords[q,0],y+quad_coords[q,1],z+quad_coords[q,2]) dummy_array[0,q,i,j,k] = temperature_field_store(Point(x+quad_coords[q,0],y+quad_coords[q,1],z+quad_coords[q,2])) except RuntimeError: pass - #make sure temperature sending is complete and create the new #temperature array for p, proc_rank in enumerate(send_ranks): reqs[p].wait() print("FENICS Step {} --- Send Requests Wait done: {}".format(s, MPI.Wtime())) temperature_array = comm.allreduce(dummy_array, op=MPI.SUM) if (s+1)%20 == 0: print('=== STEP {}/{} ==='.format(s+1, nb_steps)) break print("FENICS Step {} --- Step done: {}".format(s, MPI.Wtime())) end_time = MPI.Wtime() print("FENICS Step {} --- RUN \n Local rank = {}, took {} seconds".format(s, rank, end_time-start_time)) #make timestamp array global by broadcasting timestamp_arr = comm_py.bcast(timestamp_arr, root = 0) #if on the muSpectre communicator, update the global timestamp if color == 1: file_io_object_write.update_global_attribute('timestamp', 'timestamp', timestamp_arr) file_io_object_write.close() print("Total run time:", timer.time()-start) if __name__=='__main__': main() diff --git a/mechanical_model/run_mechanical_model_fidis.sh b/mechanical_model/run_mechanical_model_fidis.sh index 3bcb151..2876fee 100644 --- a/mechanical_model/run_mechanical_model_fidis.sh +++ b/mechanical_model/run_mechanical_model_fidis.sh @@ -1,30 +1,32 @@ #!/bin/bash +#SBATCH --partition parallel #SBATCH --nodes 1 -#SBATCH --cpus-per-task 20 +#SBATCH --cpus-per-task 24 #SBATCH --ntasks-per-node 1 -#SBATCH --mem 48G -#SBATCH --time 15:59:00 +#SBATCH --mem 60G +#SBATCH --time 00:59:00 #SBATCH --error=slurm-%j.stderr #SBATCH --output=slurm-%j.stdout #SBATCH --job-name=mech_model source /home/kubilay/anaconda3/bin/activate conda activate muspectre_env2 export LD_LIBRARY_PATH=/home/kubilay/pnetcdf/lib/:$LD_LIBRARY_PATH echo "*** STARTING JOB ***" LOG_FILE="mechanical_model.log" +mpirun -np 24 python3 residual_stress_analysis.py 60 23 > $LOG_FILE +#srun python3 netcdf2h5.py 60 >> $LOG_FILE +mpirun -np 24 python3 residual_stress_analysis.py 61 23 >> $LOG_FILE +#srun -n 1 --exclusive python3 netcdf2h5.py 61 >> $LOG_FILE -mpirun -np 20 python3 residual_stress_analysis.py 60 19 > $LOG_FILE -python3 netcdf2h5.py 60 >> $LOG_FILE -mpirun -np 20 python3 residual_stress_analysis.py 61 19 >> $LOG_FILE -python3 netcdf2h5.py 61 >> $LOG_FILE +conda deactivate echo "*** JOB FINISHED ***" diff --git a/mechanical_model/transfer_previous_layer.py b/mechanical_model/transfer_previous_layer.py index ec7e59d..0f696e0 100644 --- a/mechanical_model/transfer_previous_layer.py +++ b/mechanical_model/transfer_previous_layer.py @@ -1,372 +1,371 @@ #!/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 +import gc 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') + os.chdir('/scratch/kubilay/additive-manufacturing-work/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" + filename = "/scratch/kubilay/layer_"+str(layer_no).zfill(3)+".xdmf" + temperature_filename = "/scratch/kubilay/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.00018 + dx, dy, dz = 0.00020, 0.00020, 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_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') # 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 = 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 - \ No newline at end of file +