From dea69b3fe288ad4a90014ee3ec3ba91959f46579 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Jan 2020 11:39:56 -0600 Subject: [PATCH 1/2] Skip cell card line if it's just a comment. --- csg2csg/MCNPInput.py | 156 +++++++++++++++++++++---------------------- 1 file changed, 78 insertions(+), 78 deletions(-) diff --git a/csg2csg/MCNPInput.py b/csg2csg/MCNPInput.py index 2a41d4f..1a10f5a 100644 --- a/csg2csg/MCNPInput.py +++ b/csg2csg/MCNPInput.py @@ -14,8 +14,8 @@ from collections import Counter from numpy import linalg as np -from numpy import dot -from numpy import cross +from numpy import dot +from numpy import cross from numpy import inf as npinf from numpy import around as nparound @@ -27,7 +27,7 @@ import re class MCNPInput(InputDeck): - """ MCNPInputDeck class - does the actuall processing + """ MCNPInputDeck class - does the actuall processing """ preserve_xsid = False @@ -74,11 +74,11 @@ def __process_importances(self): repeat = int(value.replace("r","")) - 1 # this is safe since r cannot be the first # value in the list - last_importance = importance_list[idx-1] + last_importance = importance_list[idx-1] to_insert = [last_importance]*repeat importance_list[idx] = ' '.join(str(e) for e in to_insert) # flatten the list - self.importance_list[particle] = ' '.join(str(e) for e in importance_list) + self.importance_list[particle] = ' '.join(str(e) for e in importance_list) return @@ -96,10 +96,10 @@ def __get_importances(self, start_line): # check for importance keyword if "imp" in self.file_lines[idx]: particle = self.file_lines[idx].split()[0].split(":")[1] - - # TODO mcnp allows the following forms imp:n imp:n,p etc + + # TODO mcnp allows the following forms imp:n imp:n,p etc particle = mcnpToParticle(particle) - logging.debug("%s", "found importance statement for particle " + + logging.debug("%s", "found importance statement for particle " + particleToGeneric(particle) + " on line" + str(idx)) self.importance_list[particle] = self.file_lines[idx][5:].rstrip() @@ -113,7 +113,7 @@ def __get_importances(self, start_line): continue else: # otherwise advance the line by one - idx += 1 + idx += 1 return @@ -122,7 +122,7 @@ def __get_importances(self, start_line): Supports at least two forms of reading weight window information - 1) a wwe card + 1) a wwe card """ def __get_weight_windows(self, start_line): return @@ -146,7 +146,7 @@ def __get_transform_cards(self, start_line): if "tr" in self.file_lines[line]: self.__make_transform_card(self.file_lines[line]) - + line += 1 return @@ -158,9 +158,9 @@ def __get_material_card(self, start_line): mat_num = tokens[0] mat_num = mat_num.replace("m","") # set the material number - + # rebuild the first mat string - material_string = ' '.join(tokens[1:]) + " " + material_string = ' '.join(tokens[1:]) + " " if '$' in material_string: pos = material_string.find('$') material_string = material_string[:pos] @@ -184,7 +184,7 @@ def __get_material_card(self, start_line): line = line[:pos] material_string += line else: # else we have found a new cell card - break + break # increment the line that we are looking at idx += 1 break @@ -195,7 +195,7 @@ def __get_material_card(self, start_line): # multiplier material material.material_colour = get_material_colour(len(self.material_list)) self.material_list[material.material_number] = material - + return # get the material cards definitions @@ -218,7 +218,7 @@ def __get_material_cards(self, start_line): # get the material cards definitions def __get_transform_cards(self, start_line): - idx = start_line + idx = start_line while True: #cell_line = self.file_lines[jdx] if idx == len(self.file_lines): @@ -233,13 +233,13 @@ def __get_transform_cards(self, start_line): logging.debug("%s", "trn card has continue line " + str(idx)) card_line += self.file_lines[idx] idx += 1 - idx -=1 + idx -=1 self.__make_transform_card(card_line) card_line = "" idx += 1 return - + while True: if idx == len(self.file_lines): break @@ -270,10 +270,10 @@ def __apply_universe_transformations(self): if cell.cell_universe_transformation_id is not "0": # apply the transform cell.apply_universe_transform(self.transform_list[cell.cell_universe_transformation_id]) - - # find the next free material number - def __next_free_int(self): + + # find the next free material number + def __next_free_int(self): idx = 1 while True: if str(idx) in self.material_list.keys(): @@ -282,7 +282,7 @@ def __next_free_int(self): break return str(idx) - # reorganise materials such that we get a new set of unique + # reorganise materials such that we get a new set of unique # material number/ density pairs - this to avoid mcnp's # overloadable materials def __reorganise_materials(self): @@ -361,11 +361,11 @@ def __macro_rcc_cylinder_arbitrary(self,Surface,vector): # form the gq quadratic terms gq_coeffs[0] = 1. - axis_vector[0]**2 - gq_coeffs[1] = 1. - axis_vector[1]**2 - gq_coeffs[2] = 1. - axis_vector[2]**2 + gq_coeffs[1] = 1. - axis_vector[1]**2 + gq_coeffs[2] = 1. - axis_vector[2]**2 # form the rotational terms - gq_coeffs[3] = -2.*axis_vector[0]*axis_vector[1] - gq_coeffs[4] = -2.*axis_vector[1]*axis_vector[2] + gq_coeffs[3] = -2.*axis_vector[0]*axis_vector[1] + gq_coeffs[4] = -2.*axis_vector[1]*axis_vector[2] gq_coeffs[5] = -2.*axis_vector[0]*axis_vector[2] # form the linear offset terms gq_coeffs[6] = -Surface.surface_coefficients[1]*gq_coeffs[3] \ @@ -413,7 +413,7 @@ def __macro_rcc_cylinder_arbitrary(self,Surface,vector): # plane offset 2 d2 = axis_vector[0]*(Surface.surface_coefficients[0] + vector[0]) \ + axis_vector[1]*(Surface.surface_coefficients[1] + vector[1]) \ - + axis_vector[2]*(Surface.surface_coefficients[2] + vector[2]) + + axis_vector[2]*(Surface.surface_coefficients[2] + vector[2]) self.last_free_surface_index += 1 surface_string = str(self.last_free_surface_index) + " p " @@ -461,11 +461,11 @@ def __macro_rcc_cylinder(self,Surface,vector): cylinder = " c/x " c1 = Surface.surface_coefficients[1] c2 = Surface.surface_coefficients[2] - top = Surface.surface_coefficients[3] + Surface.surface_coefficients[0] + top = Surface.surface_coefficients[3] + Surface.surface_coefficients[0] bottom = Surface.surface_coefficients[0] - # if coefficients 4 & 5 are zero then its a cz with planes at + # if coefficients 4 & 5 are zero then its a cz with planes at self.last_free_surface_index += 1 surf = MCNPSurfaceCard(str(self.last_free_surface_index) + cylinder + str(c1) + " " + @@ -539,19 +539,19 @@ def explode_macrobody(self,Surface): cell_description_outside += " : " + str(new_surf_list[4].surface_id) cell_description_outside += " : " + str(new_surf_list[5].surface_id) cell_description_outside += ")" - + cell_description = [cell_description_inside,cell_description_outside] - + elif Surface.surface_type == SurfaceCard.SurfaceType["MACRO_RCC"]: id = int(Surface.surface_id) - + vector = [Surface.surface_coefficients[3],Surface.surface_coefficients[4],Surface.surface_coefficients[5]] new_surf_list, cell_description = self.__macro_rcc_cylinder_arbitrary(Surface,vector) elif Surface.surface_type == SurfaceCard.SurfaceType["MACRO_BOX"]: id = int(Surface.surface_id) - + origin = [Surface.surface_coefficients[0], Surface.surface_coefficients[1], Surface.surface_coefficients[2]] vec1 = [Surface.surface_coefficients[3], Surface.surface_coefficients[4], Surface.surface_coefficients[5]] @@ -562,14 +562,14 @@ def explode_macrobody(self,Surface): vec2n = vec2/np.norm(vec2) vec3n = vec3/np.norm(vec3) - d1 = vec1n[0]*origin[0] + vec1n[1]*origin[1] + vec1n[2]*origin[2] + d1 = vec1n[0]*origin[0] + vec1n[1]*origin[1] + vec1n[2]*origin[2] d2 = vec1n[0]*(origin[0] + vec1[0]) + vec1n[1]*(origin[1]+vec1[1]) + vec1n[2]*(origin[2]+vec1[2]) d3 = vec2n[0]*origin[0] + vec2n[1]*origin[1] + vec2n[2]*origin[2] d4 = vec2n[0]*(origin[0] + vec2[0]) + vec2n[1]*(origin[1]+vec2[1]) + vec2n[2]*(origin[2]+vec2[2]) d5 = vec3n[0]*origin[0] + vec3n[1]*origin[1] + vec3n[2]*origin[2] d6 = vec3n[0]*(origin[0] + vec3[0]) + vec3n[1]*(origin[1]+vec3[1]) + vec3n[2]*(origin[2]+vec3[2]) - + # cannonical facet ordering is done such that the surfaces # making up the macrobody all point inwards p1 = self.__make_new_plane(vec1n,d2) @@ -589,7 +589,7 @@ def explode_macrobody(self,Surface): cell_description_inside += " -" + str(new_surf_list[2].surface_id) cell_description_inside += " -" + str(new_surf_list[3].surface_id) cell_description_inside += " -" + str(new_surf_list[4].surface_id) - cell_description_inside += " -" + str(new_surf_list[5].surface_id) + cell_description_inside += " -" + str(new_surf_list[5].surface_id) cell_description_inside += ")" cell_description_outside = "(" @@ -598,17 +598,17 @@ def explode_macrobody(self,Surface): cell_description_outside += ":" + str(new_surf_list[2].surface_id) cell_description_outside += ":" + str(new_surf_list[3].surface_id) cell_description_outside += ":" + str(new_surf_list[4].surface_id) - cell_description_outside += ":" + str(new_surf_list[5].surface_id) + cell_description_outside += ":" + str(new_surf_list[5].surface_id) cell_description_outside += ")" - + cell_description = [cell_description_inside,cell_description_outside] else: warnings.warn('Found an unsupported macrobody, files will not be correct',Warning) cell_description = ["",""] - + return cell_description, new_surf_list - # if we find a macrobody in the surface list + # if we find a macrobody in the surface list # explode it into a surface based definition def __flatten_macrobodies(self): # look through the list until we find @@ -639,14 +639,14 @@ def __flatten_macrobodies(self): while True: # cell text description is contually updated cell_text_description = cell.cell_text_description - + # if we find the surface id of the macrobdy in the text description sub = str(surf.surface_id) regex = re.compile("^-?("+str(surf.surface_id)+")(\.+[1-9])?$") matches = [m.group(0) for l in cell_text_description for m in [regex.search(l)] if m] #if str(surf.surface_id) in cell_text_description or str(surf.surface_id)+"." in cell_text_description: - - if matches: + + if matches: # loop over each component and find the macrobody for idx, surface in enumerate(cell.cell_text_description): # if it matches we have the simmple form @@ -655,9 +655,9 @@ def __flatten_macrobodies(self): cell.cell_text_description[idx] = cell_description[1] elif "-"+str(surf.surface_id) == surface: cell.cell_text_description[idx] = cell_description[0] - + # else we have the facet form - if str(surf.surface_id)+"." in surface: + if str(surf.surface_id)+"." in surface: surface_index = int(surface.split(".")[1]) # get just the mcnp surface index new_surface_id = new_surfaces[surface_index-1].surface_id # mcnp numbers them 1->n if "-" in surface: # need to take care of the -sign @@ -669,7 +669,7 @@ def __flatten_macrobodies(self): # update the text description text_string = ' '.join(cell.cell_text_description) self.cell_list[jdx].update(text_string) - + # clear up removed surfaces logging.debug("%s", "Deleting macrobody surfaces") for surf in to_remove: @@ -690,12 +690,12 @@ def __nots_remaining(self, cell): def __split_nots(self,cell): # if the cell has a not in it pos = 0 - count = 0 + count = 0 # loop over the constituents of the cell # we have to do this rookie looking stuff # because cell.OperationType is not iterable for i in cell.cell_interpreted: - # if its not an operation + # if its not an operation if not isinstance(i,cell.OperationType): if "#" in i: pos = count @@ -715,17 +715,17 @@ def __split_nots(self,cell): # build the cell into the interpreted form cell_text = [cell.OperationType(2)] + ["("] + cell_text cell_text = cell_text + [")"] - + # remove the #cell and insert the full new form cell_part = cell.cell_interpreted[:pos] - + cell_part.extend(cell_text) cell_part2 = cell.cell_interpreted[pos+1:] cell_part.extend(cell_part2) cell.cell_interpreted = cell_part return - + # loop through the cells and insert # cell definititons where needed # assuming that we have nots of the form #33 #44 and @@ -735,9 +735,9 @@ def __explode_nots(self): while self.__nots_remaining(cell): self.__split_nots(cell) continue - - # generate bounding coordinates + + # generate bounding coordinates def __generate_bounding_coordinates(self): # loop through the surfaces and generate the bounding coordinates # condisder only manifold or simple infinite surfaces, like planes @@ -755,11 +755,11 @@ def __generate_bounding_coordinates(self): self.bounding_coordinates[4] = box[4] if box[5] > self.bounding_coordinates[5]: self.bounding_coordinates[5] = box[5] - logging.debug("%s ", "bounding box of geometry is " + str(self.bounding_coordinates[0]) + " " + + logging.debug("%s ", "bounding box of geometry is " + str(self.bounding_coordinates[0]) + " " + str(self.bounding_coordinates[1]) + " " + - str(self.bounding_coordinates[2]) + " " + - str(self.bounding_coordinates[3]) + " " + - str(self.bounding_coordinates[4]) + " " + + str(self.bounding_coordinates[2]) + " " + + str(self.bounding_coordinates[3]) + " " + + str(self.bounding_coordinates[4]) + " " + str(self.bounding_coordinates[5]) + "\n") # update surfaces that need their bounding coordinates updated @@ -789,7 +789,7 @@ def __get_cell_cards(self): while True: cell_line = self.file_lines[idx] # this relies upon us checking correctly for all other - # cases where there may be blank lines due to our + # cases where there may be blank lines due to our # processing of the string if cell_line.isspace(): logging.info('%s',"found end of cell cards at line " + str(idx)) @@ -805,9 +805,9 @@ def __get_cell_cards(self): cell_comment = "" if pos_comment != -1: cell_line = cell_line[:pos_comment] - cell_comment = self.file_lines[jdx][pos_comment:] # set the comment + cell_comment = self.file_lines[jdx][pos_comment:] # set the comment self.file_lines[jdx] = cell_line # update the file data - + # mcnp continue line is indicated by 5 spaces if cell_line[0:5] == " " and not cell_line.isspace(): card_line += cell_line @@ -817,18 +817,20 @@ def __get_cell_cards(self): pass else: # else we have found a new cell card logging.debug("%s\n", "Found new cell card " + card_line) + if card_line == 'c': + continue cellcard = MCNPCellCard(card_line) # we should set the comment here self.cell_list.append(cellcard) - break + break jdx += 1 - idx = jdx + idx = jdx return idx # set the boundary conditions def __apply_boundary_conditions(self): - - # apply the importances to cells + + # apply the importances to cells if len(self.importance_list) != 0: # TODO make this loop apply to multiple particle # types but for now just do neutrons @@ -837,10 +839,10 @@ def __apply_boundary_conditions(self): for idx,value in enumerate(importances): self.cell_list[idx].cell_importance = float(value) - # loop over the cells and if the cell has + # loop over the cells and if the cell has # importance 0, all the sufaces get boundary - # condition - for cell in self.cell_list: + # condition + for cell in self.cell_list: if cell.cell_importance == 0: for surf in cell.cell_surface_list: self.get_surface_with_id(surf).boundary_condition = SurfaceCard.BoundaryCondition["VACUUM"] @@ -867,11 +869,11 @@ def __get_surface_cards(self,idx): surfacecard = MCNPSurfaceCard(surf_card) self.surface_list.append(surfacecard) # update the surface index counter - if surfacecard.surface_id > self.last_free_surface_index: + if surfacecard.surface_id > self.last_free_surface_index: self.last_free_surface_index = surfacecard.surface_id - break + break jdx += 1 - idx = jdx + idx = jdx return idx # process the mcnp input deck and read into a generic datastructure @@ -894,7 +896,7 @@ def process(self): logging.debug("%s", "Input Echo") for idx,line in enumerate(self.file_lines): logging.debug("%i %s",idx,line) - + # get the cell cards idx = self.__get_cell_cards() @@ -912,14 +914,14 @@ def process(self): # also idx will never be advanced from this point # try to get importances defined in the data block - self.__get_importances(idx) + self.__get_importances(idx) # try to get weights defined the data block self.__get_weight_windows(idx) self.__get_transform_cards(idx) self.__get_material_cards(idx) - # need to flatten first to get transformed surface in the - # correct place + # need to flatten first to get transformed surface in the + # correct place self.__flatten_macrobodies() self.__explode_nots() @@ -940,10 +942,10 @@ def process(self): self.__update_surfaces() self.split_unions() - + return - # perhaps these write functions should actually build strings + # perhaps these write functions should actually build strings # and then write at once? # write all surfaces to mcnp format def __write_mcnp_surfaces(self, filestream): @@ -973,5 +975,3 @@ def write_mcnp(self, filename, flat = True): self.__write_mcnp_surfaces(f) self.__write_mcnp_materials(f) f.close() - - From c9d247a82d95463e89ca788457e11dd0e06b9401 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Jan 2020 14:07:14 -0600 Subject: [PATCH 2/2] Adding case for comment line in cell card definitions --- csg2csg/MCNPInput.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/csg2csg/MCNPInput.py b/csg2csg/MCNPInput.py index 1a10f5a..28ee46e 100644 --- a/csg2csg/MCNPInput.py +++ b/csg2csg/MCNPInput.py @@ -796,6 +796,11 @@ def __get_cell_cards(self): idx += 1 break + if cell_line.strip() == 'c': + logging.info('%s',"Found comment only line") + idx +=1 + continue + card_line = cell_line jdx = idx + 1 # scan until we are all done @@ -817,8 +822,6 @@ def __get_cell_cards(self): pass else: # else we have found a new cell card logging.debug("%s\n", "Found new cell card " + card_line) - if card_line == 'c': - continue cellcard = MCNPCellCard(card_line) # we should set the comment here self.cell_list.append(cellcard)