Coverage for bim2sim / plugins / PluginOpenFOAM / bim2sim_openfoam / utils / openfoam_utils.py: 0%
187 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-18 09:34 +0000
1import OCC.Core.TopoDS
2from OCC.Core.Extrema import Extrema_ExtFlag_MIN
3from OCC.Core.TopOpeBRep import TopOpeBRep_ShapeIntersector
4from OCC.Core.gp import gp_Pnt
5from OCC.Core.BRepExtrema import BRepExtrema_DistShapeShape
6import bim2sim.tasks.common.inner_loop_remover as ilr
7from bim2sim.utilities.common_functions import filter_elements
8from bim2sim.utilities.pyocc_tools import PyOCCTools
9import math
10from stl import mesh
11import numpy as np
12import re
15class OpenFOAMUtils:
16 @staticmethod
17 def split_openfoam_elements(openfoam_elements: dict) -> tuple[list, list,
18 list, list, list]:
19 stl_bounds = filter_elements(openfoam_elements, 'StlBound')
20 heaters = filter_elements(openfoam_elements, 'Heater')
21 air_terminals = filter_elements(openfoam_elements, 'AirTerminal')
22 furniture = filter_elements(openfoam_elements, 'Furniture')
23 people = filter_elements(openfoam_elements, 'People')
24 return stl_bounds, heaters, air_terminals, furniture, people
26 @staticmethod
27 def get_refinement_level(dist: float, bM_size: float, mean_dist:
28 float=None) -> list:
29 """
30 Computes the refinement level based on the desired
31 refined cell size and the blockMesh cell size.
32 bM_size/2^X = min_dist
33 """
34 ref_level = math.log(bM_size / dist) / math.log(2)
35 ref_level = math.ceil(ref_level)
36 if mean_dist:
37 min_level = math.log(bM_size / mean_dist) / math.log(2)
38 min_level = int(round(min_level,0))
39 if min_level < 0:
40 min_level = 0
41 if min_level > ref_level:
42 min_level = ref_level
43 else:
44 min_level = ref_level
45 return [min_level, ref_level + 2]
47 @staticmethod
48 def get_min_refdist_between_shapes(shape1: OCC.Core.TopoDS.TopoDS_Shape,
49 shape2: OCC.Core.TopoDS.TopoDS_Shape,
50 dist_bound=0.05) -> float:
51 """
52 Computes the minimal distance between two TopoDS Shapes and returns
53 the distance divided by 3 such that in the refinement zone,
54 there will be 3 cells in between two objects. Optional argument
55 dist_bound can specify a maximal minimal distance (default: 1cm).
56 """
57 ins = TopOpeBRep_ShapeIntersector()
58 ins.InitIntersection(shape1, shape2)
59 if ins.MoreIntersection():
60 return 0
61 nb1 = PyOCCTools.get_number_of_vertices(shape1)
62 nb2 = PyOCCTools.get_number_of_vertices(shape2)
63 if nb1+nb2 > 1e3:
64 # avoid expensive distance calculations for far complex elements
65 box1 = PyOCCTools.simple_bounding_box_shape([shape1])
66 box2 = PyOCCTools.simple_bounding_box_shape([shape2])
67 box_dist = BRepExtrema_DistShapeShape(box1, box2,
68 Extrema_ExtFlag_MIN).Value()
69 if box_dist/3 > dist_bound:
70 return dist_bound
71 if (nb1+nb2) < 200:
72 dist = BRepExtrema_DistShapeShape(shape1, shape2,
73 Extrema_ExtFlag_MIN).Value()
74 else:
75 # apply point based distance computation for cases with a large
76 # amount of vertices. OpenCascade is very slow on these kind of
77 # distance computations, so an approximate solution is calculated
78 dist = PyOCCTools.calculate_point_based_distance(shape1, shape2)
79 # To ensure correctness of all boundary conditions, there must be at
80 # least 3 cells between the object and the wall.
81 min_dist = dist / 3
82 if min_dist > dist_bound:
83 min_dist = dist_bound
84 return min_dist
86 @staticmethod
87 def get_min_internal_dist(points: list[gp_Pnt]) -> float:
88 """
89 Computes the minimal internal distance of a TopoDS Shape.
90 """
91 mindist = 0.05
92 for i, p1 in enumerate(points):
93 for p2 in points[i:]:
94 dist = p1.Distance(p2)
95 if mindist > dist > 5e-04:
96 mindist = dist
97 return mindist
99 def validate_volumes_cells(self, case: str, inlet_type: str,
100 outlet_type: str) -> bool:
101 """
102 Compute the absolute difference in the volumes of the geometry of
103 the room and the volume of the discretized domain.
104 If the V directory was not created in the controlDict, it can be
105 created after meshing using (Linux):
106 postProcess -func writeCellVolumes -time 0
107 gzip -d V.gz
108 The total discretized volume can also be obtained from the checkMesh
109 functionality.
111 case must be the path to the corresponding OpenFOAM case directory.
112 inlet_type and outlet_type must be equal to the
113 corresponding sim_setting options.
114 """
115 base_dir = case + '/' + 'constant' + '/' + 'triSurface' + '/'
116 vol_geom = self.get_vol_from_geometry(base_dir, inlet_type, outlet_type)
117 vol_cells = self.get_vol_from_cells(case)
118 abs_diff = abs(vol_geom - vol_cells)
119 if abs_diff < 0.05:
120 mesh_valid = True
121 else:
122 mesh_valid = False
123 return mesh_valid
125 @staticmethod
126 def get_vol_from_geometry(base_dir: str, inlet_type: str,
127 outlet_type: str) -> float:
128 """
129 Compute the volume of the simulated room based on the STL files.
131 base_dir must be the path to the triSurface directory of the
132 corresponding OpenFOAM case.
133 inlet_type and outlet_type must be equal to the
134 corresponding sim_setting options.
135 """
136 partMesh = []
137 with open(base_dir + 'space_2RSCzLOBz4FAK$_wE8VckM.stl', 'r') as \
138 multimesh:
139 mesh_lines = multimesh.readlines()
140 start_ind = end_ind = -1
141 for line in mesh_lines:
142 if 'endsolid' in line:
143 start_ind = end_ind + 1
144 end_ind = mesh_lines.index(line)
145 with open('temp_mesh.stl', 'w') as temp_mesh:
146 temp_mesh.writelines(mesh_lines[start_ind:end_ind + 1])
147 temp_mesh.close()
148 partMesh.append(mesh.Mesh.from_file('temp_mesh.stl'))
149 partMesh.append(mesh.Mesh.from_file(base_dir + 'inlet_diffuser.stl'))
150 partMesh.append(mesh.Mesh.from_file(base_dir + 'outlet_diffuser.stl'))
151 if not inlet_type == 'SimpleStlDiffusor':
152 partMesh.append(mesh.Mesh.from_file(base_dir + 'inlet_box.stl'))
153 partMesh.append(mesh.Mesh.from_file(base_dir +
154 'inlet_source_sink.stl'))
155 if not outlet_type == 'SimpleStlDiffusor':
156 partMesh.append(mesh.Mesh.from_file(base_dir + 'outlet_box.stl'))
157 partMesh.append(mesh.Mesh.from_file(base_dir +
158 'outlet_source_sink.stl'))
159 combined = partMesh[0]
160 for part in partMesh[1:]:
161 combined = mesh.Mesh(np.concatenate((combined.data, part.data)))
162 vol_geom = float(combined.get_mass_properties()[0])
163 return vol_geom
165 @staticmethod
166 def get_vol_from_cells(case) -> float:
167 """
168 Compute the volume of the simulated room based on the FVM cells.
170 case must be the path to the corresponding OpenFOAM case directory.
171 """
172 with open(case + '/' + '0' + '/' + 'V') as f:
173 lines = f.readlines()
174 n_cells = int(lines[21])
175 total_vol = 0
176 for line in lines[23:n_cells + 23]:
177 total_vol += float(line)
178 return total_vol
180 def detriangulize(self, obj):
181 """2-step algorithm for removing triangularizaton from an obj as
182 TopoDS_Compound. 1. remove inner edges. 2. remove collinear and
183 coincident points."""
184 triang = ilr._get_triangulation(obj)
185 inner_edges, outer_edges = ilr._get_inside_outside_edges(triang,
186 must_equal=False)
187 edges = inner_edges + outer_edges
188 vertices = PyOCCTools.get_unique_vertices(outer_edges)
189 vertices = self.remove_coincident_vertices(vertices)
190 vertices = PyOCCTools.remove_collinear_vertices2(vertices)
191 return vertices, edges
193 @staticmethod
194 def get_edge_lengths(edges):
195 """
196 Calculate the lengths of edges in 3D space.
198 Parameters:
199 edges (list of tuples): A list where each element is a tuple
200 containing two 3D coordinates
201 (e.g., [((x1, y1, z1), (x2, y2, z2)), ...]).
203 Returns:
204 list: A list of lengths corresponding to the edges.
205 """
206 lengths = []
207 for (x1, y1, z1), (x2, y2, z2) in edges:
208 length = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)
209 lengths.append(length)
210 return lengths
212 @staticmethod
213 def remove_coincident_vertices(vert_list: list) -> list:
214 """Slightly modified version of the method in PyOCCTools. Remove
215 coincident vertices from list of gp_Pnt. Vertices are coincident if
216 closer than tolerance. Does not assume vertices to be sorted."""
217 tol_dist = 1e-3
218 remove_list = []
219 for i, vert in enumerate(vert_list):
220 for vert2 in vert_list[i + 1:]:
221 v = np.array(vert.Coord())
222 v2 = np.array(vert2.Coord())
223 d_b = np.linalg.norm(v - v2)
224 if d_b < tol_dist:
225 remove_list.append(vert2)
226 for v in remove_list:
227 if v in vert_list:
228 vert_list.remove(v)
229 return vert_list
231 @staticmethod
232 def string_to_dict(s):
233 # Extract the tuples from the string using regular expressions
234 matches = re.findall(r"\((\d+)\s([\d.]+)\)", s)
235 # Convert the matches to a dictionary
236 return {int(k): float(v) for k, v in matches}
238 @staticmethod
239 def dict_to_string(d):
240 # Format the dictionary back into the original string format
241 sorted_items = sorted(d.items())
242 tuples_str = " ".join(f"({k} {v})" for k, v in sorted_items)
243 return f"table ( {tuples_str} )"
245 @staticmethod
246 def float_cutoff(val: float, max_dig: int = 6) -> float:
247 """Round long float values after max_dig digits.
249 Round long float values val after max_dig digits in OpenFOAM case
250 files for higher stability in regression tests.
252 Args:
253 val: float value to be rounded.
254 max_dig: maximum number of digits in float values.
256 Returns:
257 rounded float value.
258 """
259 if val is None:
260 return val
261 return round(val, max_dig)
263 @staticmethod
264 def duplicate_table_for_restart(dict_with_string_tables: dict,
265 add_number_to_keys: int) -> dict:
266 new_dict = {}
267 for key, eq_val in dict_with_string_tables.items():
268 d = OpenFOAMUtils.string_to_dict(eq_val)
269 if d:
270 d_updated = {}
271 for d_key, d_val in d.items():
272 d_updated.update({d_key: d_val})
273 d_updated.update({d_key + add_number_to_keys: d_val})
274 s_updated = OpenFOAMUtils.dict_to_string(d_updated)
275 new_dict.update({key: s_updated})
276 else:
277 new_dict.update({key: eq_val})
278 return new_dict
280 @staticmethod
281 def prime_factors(n):
282 factors = []
283 divisor = 2
284 while n > 1:
285 while n % divisor == 0:
286 factors.append(divisor)
287 n //= divisor
288 divisor += 1
289 return factors
291 @staticmethod
292 def split_into_three_factors(n):
293 factors = OpenFOAMUtils.prime_factors(n)
294 factors.sort(reverse=True) # Start with the largest primes first
295 groups = [1, 1, 1]
297 # Distribute factors to minimize the difference among groups
298 for factor in factors:
299 groups.sort() # Always add to the smallest group
300 groups[0] *= factor
302 # Sort the final groups to ensure descending order
303 groups.sort(reverse=True)
304 return groups