Page MenuHomec4science
Files F88764623

adjust_previous_temperature.py
No OneTemporary
Actions

File Metadata

Created
Sun, Oct 20, 14:02

adjust_previous_temperature.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 31 09:31:24 2022
@author: ekinkubilay
"""
from fenics import *
import numpy as np
import argparse
def parse_args():
"""
Function to handle arguments
"""
parser = argparse.ArgumentParser(description='Run FEM of single scan on a simple box')
parser.add_argument("layer_no", type=int, help='layer number')
args = parser.parse_args()
return args
def main():
#get the arguments
args = parse_args()
current_layer_no = args.layer_no
#allow extrapolation for now
parameters['allow_extrapolation'] = True
#set directory names
mesh_dir = '../Part_geometry/mesh/'
temp_data_dir = 'results/temperature/'
#set file names
previous_layer_no = current_layer_no - 1
previous_mesh_name = 'layer_{}.xdmf'.format(str(previous_layer_no).zfill(3))
current_mesh_name = 'layer_{}.xdmf'.format(str(current_layer_no).zfill(3))
previous_temperature_name = 'temperature_layer_{}.h5'.format(str(previous_layer_no).zfill(3))
previous_mesh_file = mesh_dir+previous_mesh_name
current_mesh_file = mesh_dir+current_mesh_name
previous_temperature_file = temp_data_dir+previous_temperature_name
current_temperature_file = temp_data_dir+'adjusted_temperature_layer_{}.h5'.format(str(current_layer_no).zfill(3))
#load old previous layer mesh
try:
mesh = Mesh()
filename = previous_mesh_file
#filename = '../Part_geometry/mesh/layer_009.xdmf'
f = XDMFFile(MPI.comm_world, filename)
f.read(mesh)
MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m
f.close()
except RuntimeError:
return 0
z_max_old = np.max(mesh.coordinates()[:,2])
#load previous layer Temperature function
#hdf5_file = HDF5File(MPI.comm_world, 'results/trial.h5', "r")
hdf5_file = HDF5File(MPI.comm_world, previous_temperature_file, "r")
#define function space for V
V = FunctionSpace(mesh, "CG", 1)
#define temperature of old layer
T_old = Function(V)
#define the last time step
t = hdf5_file.attributes('Temperature')['count']-1
#define name of the vector
vec_name = "/Temperature/vector_%d"%t
#read the temperature of the last time step of previous layer
hdf5_file.read(T_old, vec_name)
timestamp = hdf5_file.attributes(vec_name)["timestamp"]
hdf5_file.close()
#load current layer mesh
mesh_new = Mesh()
#filename = '../Part_geometry/mesh/layer_010.xdmf'
filename = current_mesh_file
f = XDMFFile(MPI.comm_world, filename)
f.read(mesh_new)
MeshTransformation.rescale(mesh_new, 1e-3, Point(0,0,0)) #in m
f.close()
V_new = FunctionSpace(mesh_new, "CG", 1)
T_0 = Function(V_new)
y = mesh_new.coordinates()[dof_to_vertex_map(V_new)]
z_max_new = np.max(mesh_new.coordinates()[:,2])
LagrangeInterpolator.interpolate(T_0, T_old)
parameters['allow_extrapolation'] = False
array = T_0.vector()[:]
for i in range(len(array)):
if near(y[i][2], z_max_new, 0.000001):
try:
array[i] = T_old(y[i][0], y[i][1], z_max_old)
except RuntimeError:
pass
#else:
#array[i] = T_old(y[i,0], y[i,1], y[i,2])
T_0.vector().set_local(array)
hdf5_file = HDF5File(MPI.comm_world, current_temperature_file, "w")
T_0.rename("Temperature", "Temperature")
hdf5_file.write(T_0, "Temperature")
hdf5_file.close()
return 0
if __name__=='__main__':
main()

Event Timeline

c4science · Help