7
The Linear_Interp_1D class will read a 1d solution field,
8
then find the interpolated solution value at specified
11
The driver main will read a interpolation mesh and
12
then output the interpolated solution field.
14
The x coordinate assoicated with the solution field must
17
The driver main is used via:
19
Linear_Interp.py solution_file_name interpolation_mesh_file_name interpolation_solution_file_name
26
class Linear_Interp_1D:
28
def __init__(self, solution_file_name):
29
""" Initialise Linear_Interp_1D class """
31
# open two input files
32
self.solution_file = open(solution_file_name,"r")
34
# read the solution file and x coord into lists
35
self.read_solution_file()
37
# check that the coord mesh has ascending x values
38
self.check_mesh_has_ascending_x_coord()
40
# if number of x_coord not bigger than 1 then not going to work
41
if self.num_x_coord < 2:
43
print "Error as x coord of solution has less than 2 values"
47
def read_solution_file(self):
48
""" Read the solution input file """
50
# initialise the input lists
55
self.solution_file.seek(0)
57
# This assume a particular file format
62
for line in self.solution_file:
64
# skip blank lines (ie a line that cannot be stripped)
73
# count num_x_coord found
74
self.num_x_coord = self.num_x_coord + 1
76
self.x_coord.append(float(line.split()[0]))
77
self.solution.append(float(line.split()[1]))
79
def check_mesh_has_ascending_x_coord(self):
80
""" Check that the mesh has an ascending x coordinate """
82
for node in range(2,len(self.x_coord)):
84
if self.x_coord[node] <= self.x_coord[node-1]:
86
print "Issue with solution mesh: x coordinate not ascending for nodes:",node,node-1
91
""" Find the interpolated value at coordinate x """
93
# bound the interpolation to the known x_coord min and max values
95
if x <= self.x_coord[0]:
99
elif x >= self.x_coord[-1]:
101
s = self.solution[-1]
105
# find which ele (or interval) this x resides in
106
for ele in range(1,len(self.x_coord)):
108
x_right = self.x_coord[ele]
114
# Find the x left, s left and right
115
x_left = self.x_coord[ele - 1]
117
s_right = self.solution[ele]
118
s_left = self.solution[ele-1]
121
delta_x = (x - x_left) / (x_right - x_left)
124
s = (1.0 - delta_x) * s_left + delta_x * s_right
128
if __name__ == "__main__":
132
solution_file_name = sys.argv[1]
134
interpolation_mesh_file_name = sys.argv[2]
136
interpolation_solution_file_name = sys.argv[3]
140
print "Usage: Linear_Interp_1D.py solution_file_name interpolation_mesh_file_name interpolation_solution_file_name"
144
print "Running Linear_Interp_1D.py ",solution_file_name, interpolation_mesh_file_name, interpolation_solution_file_name
147
LI = Linear_Interp_1D(solution_file_name)
149
# initialise the interpolated x coord and solution lists
153
# open the file which has the solution interpolation mesh
154
interpolation_mesh_file = open(interpolation_mesh_file_name,"r")
156
# read the interpolated solution x coord and initialise interpolated solution
157
for line in interpolation_mesh_file:
159
# skip blank lines (ie a line that cannot be stripped)
168
x_coord_interp.append(float(line.split()[0]))
170
solution_interp.append(0.0)
172
# find the interpolated solution for each of the required points
173
for x_index in range(len(x_coord_interp)):
175
solution_interp[x_index] = LI.interp(x_coord_interp[x_index])
177
# open the interpolated solution output file name
178
interpolation_solution_file = open(interpolation_solution_file_name,"w")
180
# write out the interpolated x_coord_interp solution_interp
181
for x_index in range(len(x_coord_interp)):
183
interpolation_solution_file.write(str(x_coord_interp[x_index])+" "+str(solution_interp[x_index])+"\n")
185
interpolation_solution_file.close()
added
removed
tests/darcy_impes_p1_2phase_coreyrelperm_velBCinlet_strongpressoutlet_p1satdiag_parallel/Linear_Interp_1D.py
