#280 Adding output to OpenFOAM

Open jgostick jgostick

Flags

Flags have been temporarily removed from this view while the flagging feature is refactored for better performance and user experience.

You can still use flags when viewing individual files. Flag-level thresholds will also remain on pull and merge requests in your repository provider.

More information can be found in our documentation.

Showing 2 of 4 files from the diff.

@@ -0,0 +1,446 @@
Loading
1 +
import sys
2 +
import numpy as np
3 +
from scipy.ndimage import zoom
4 +
np.set_printoptions(threshold=sys.maxsize)
5 +
6 +
7 +
class openfoam():
8 +
    """
9 +
    """
10 +
    def save(im, scale=1, zoom_factor=1, label=True):
11 +
        """
12 +
        Given a boolean numpy array where True is void and False is solid,
13 +
        creates a blockMesh file for use in OpenFOAM
14 +
        index: The boolean array
15 +
        scale: Scaling factor for vertex coords (ex. 0.001 scales to mm)
16 +
        """
17 +
18 +
        file = """
19 +
/*--------------------------------*- C++ -*----------------------------------*\
20 +
| \n =========                 |                                              |
21 +
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
22 +
|  \\    /   O peration     | Version:  v1912                                 |
23 +
|   \\  /    A nd           | Website:  www.openfoam.com                      |
24 +
|    \\/     M anipulation  |                                                 |
25 +
|*---------------------------------------------------------------------------*/
26 +
FoamFile
27 +
{
28 +
    version     2.0;
29 +
    format      ascii;
30 +
    class       dictionary;
31 +
    object      blockMeshDict;
32 +
}
33 +
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 +
35 +
scale ;
36 +
37 +
vertices
38 +
(
39 +
40 +
);
41 +
42 +
blocks
43 +
(
44 +
45 +
);
46 +
47 +
edges
48 +
(
49 +
);
50 +
51 +
boundary
52 +
(
53 +
    top
54 +
    {
55 +
        type patch;
56 +
        faces
57 +
        (
58 +
59 +
        );
60 +
    }
61 +
    back
62 +
    {
63 +
        type patch;
64 +
        faces
65 +
        (
66 +
67 +
        );
68 +
    }
69 +
    bottom
70 +
    {
71 +
        type patch;
72 +
        faces
73 +
        (
74 +
75 +
        );
76 +
    }
77 +
    front
78 +
    {
79 +
        type patch;
80 +
        faces
81 +
        (
82 +
83 +
        );
84 +
    }
85 +
    left
86 +
    {
87 +
        type patch;
88 +
        faces
89 +
        (
90 +
91 +
        );
92 +
    }
93 +
    right
94 +
    {
95 +
        type patch;
96 +
        faces
97 +
        (
98 +
99 +
        );
100 +
    }
101 +
    walls
102 +
    {
103 +
        type wall;
104 +
        faces
105 +
        (
106 +
107 +
        );
108 +
    }
109 +
);
110 +
111 +
mergePatchPairs
112 +
(
113 +
);
114 +
115 +
// ************************************************************************* //
116 +
"""
117 +
        # Applies image compression
118 +
        im = zoom(im.astype(int), zoom_factor)
119 +
120 +
        im = im.astype(bool)
121 +
122 +
        # Finds solids and unpacks coords
123 +
        x, y, z = np.where(im)
124 +
125 +
        # Setup for appending
126 +
        x = x.reshape(-1, 1) + 0.5
127 +
        y = y.reshape(-1, 1) + 0.5
128 +
        z = z.reshape(-1, 1) + 0.5
129 +
130 +
        # Combines coords to create centers matrix
131 +
        centers = np.append(x, y, axis=1)
132 +
        centers = np.append(centers, z, axis=1).reshape(-1, 3)
133 +
134 +
        # Finds boundary coords, used for labelling boundaries
135 +
        if label:
136 +
            minx = np.array([x == min(x)]).flatten()
137 +
            miny = np.array([y == min(y)]).flatten()
138 +
            minz = np.array([z == min(z)]).flatten()
139 +
140 +
            maxx = np.array([x == max(x)]).flatten()
141 +
            maxy = np.array([y == max(y)]).flatten()
142 +
            maxz = np.array([z == max(z)]).flatten()
143 +
144 +
            left_points = centers[minx] - 0.5
145 +
            bottom_points = centers[miny] - 0.5
146 +
            front_points = centers[minz] - 0.5
147 +
148 +
            right_points = centers[maxx] + 0.5
149 +
            top_points = centers[maxy] + 0.5
150 +
            back_points = centers[maxz] + 0.5
151 +
152 +
        # Matrix to convert the centers to vertices
153 +
        d = np.array([-0.5, -0.5, 0.5, 0.5, -0.5, 0.5, 0.5, 0.5, 0.5, -0.5,
154 +
                      0.5, 0.5, -0.5, -0.5, -0.5, 0.5, -0.5, -0.5, 0.5, 0.5,
155 +
                      -0.5, -0.5, 0.5, -0.5]).reshape(-1, 3)
156 +
157 +
        # Adjusting the scale
158 +
        # d = d * scale / compression
159 +
160 +
        # Matrix to face coords that have an outward facing normal
161 +
        d_left = np.array([0, 0, 0, 0, 0, 1,
162 +
                           0, 1, 1, 0, 1, 0]).reshape(-1, 3)
163 +
        d_bottom = np.array([0, 0, 0, 1, 0, 0,
164 +
                             1, 0, 1, 0, 0, 1]).reshape(-1, 3)
165 +
        d_front = np.array([0, 0, 0, 0, 1, 0,
166 +
                            1, 1, 0, 1, 0, 0]).reshape(-1, 3)
167 +
168 +
        d_right = np.array([0, 0, 0, 0, -1, 0,
169 +
                            0, -1, -1, 0, 0, -1]).reshape(-1, 3)
170 +
        d_top = np.array([0, 0, 0, 0, 0, -1,
171 +
                          -1, 0, -1, -1, 0, 0]).reshape(-1, 3)
172 +
        d_back = np.array([0, 0, 0, -1, 0, 0,
173 +
                           -1, -1, 0, 0, -1, 0]).reshape(-1, 3)
174 +
175 +
        # Adjusting the scale
176 +
        # d_left = d_left * scale / compression
177 +
        # d_bottom = d_bottom * scale / compression
178 +
        # d_front = d_front * scale / compression
179 +
180 +
        # d_right = d_right * scale / compression
181 +
        # d_top = d_top * scale / compression
182 +
        # d_back = d_back * scale / compression
183 +
184 +
        # Initialize vertices matrix
185 +
        vertices = np.zeros(len(centers) * len(d) * 3).reshape(-1, 8, 3)
186 +
187 +
        # Initialize faces matrices
188 +
        top = np.zeros((len(centers), 4), dtype=int)
189 +
        back = np.zeros((len(centers), 4), dtype=int)
190 +
        bottom = np.zeros((len(centers), 4), dtype=int)
191 +
        front = np.zeros((len(centers), 4), dtype=int)
192 +
        left = np.zeros((len(centers), 4), dtype=int)
193 +
        right = np.zeros((len(centers), 4), dtype=int)
194 +
195 +
        # Initialize blocks matrix
196 +
        blocks = np.zeros((len(centers), 8), dtype=int)
197 +
198 +
        # Combining centers and distance to create vertices, uses broadcasting
199 +
        for i in range(len(centers)):
200 +
            vertices[i] = centers[i] + d
201 +
202 +
        vertices = vertices.reshape(-1, 3)
203 +
204 +
        # Creating unique vertices matrix, as well as an index to map between
205 +
        # vertices and vertices_unique such that
206 +
        # vertices_unique[index] == vertices
207 +
        vertices_unique, index = np.unique(vertices, axis=0,
208 +
                                           return_inverse=True)
209 +
210 +
        if label:
211 +
            top_faces = boundary(top_points, d_top, vertices_unique)
212 +
            back_faces = boundary(back_points, d_back, vertices_unique)
213 +
            bottom_faces = boundary(bottom_points, d_bottom, vertices_unique)
214 +
            front_faces = boundary(front_points, d_front, vertices_unique)
215 +
            left_faces = boundary(left_points, d_left, vertices_unique)
216 +
            right_faces = boundary(right_points, d_right, vertices_unique)
217 +
218 +
            # Expanding faces (patches)
219 +
            for i in range(len(centers)):
220 +
                top[i] = np.array([2, 6, 7, 3]) + 8 * i
221 +
                back[i] = np.array([2, 3, 0, 1]) + 8 * i
222 +
                bottom[i] = np.array([4, 5, 1, 0]) + 8 * i
223 +
                front[i] = np.array([7, 6, 5, 4]) + 8 * i
224 +
                left[i] = np.array([0, 3, 7, 4]) + 8 * i
225 +
                right[i] = np.array([2, 1, 5, 6]) + 8 * i
226 +
227 +
            # Replaces duplicates vertices with their corresponding index
228 +
            # in the unique vertices list
229 +
            top = index[top]
230 +
            back = index[back]
231 +
            bottom = index[bottom]
232 +
            front = index[front]
233 +
            left = index[left]
234 +
            right = index[right]
235 +
236 +
            # Deleting all internal faces
237 +
238 +
            # Only opposing faces will overlap, makes one matrix where the face
239 +
            # definitions are sorted
240 +
            top_and_bottom = np.sort(np.vstack((top, bottom)))
241 +
            front_and_back = np.sort(np.vstack((front, back)))
242 +
            left_and_right = np.sort(np.vstack((left, right)))
243 +
244 +
            # Creates a matrix of unique face definitions as well as
245 +
            # their count
246 +
            t_and_b_unique, t_and_b_count = np.unique(top_and_bottom, axis=0,
247 +
                                                      return_counts=True)
248 +
            f_and_b_unique, f_and_b_count = np.unique(front_and_back, axis=0,
249 +
                                                      return_counts=True)
250 +
            l_and_r_unique, l_and_r_count = np.unique(left_and_right, axis=0,
251 +
                                                      return_counts=True)
252 +
253 +
            # Deletes the faces mentioned more than once
254 +
            top_and_bottom = t_and_b_unique[t_and_b_count == 1]
255 +
            front_and_back = f_and_b_unique[f_and_b_count == 1]
256 +
            left_and_right = l_and_r_unique[l_and_r_count == 1]
257 +
258 +
            # Putting ALL faces into one wall matrix
259 +
            walls = np.vstack((top_and_bottom, front_and_back, left_and_right,
260 +
                               left_faces, bottom_faces, front_faces,
261 +
                               right_faces, top_faces, back_faces))
262 +
            walls_sorted = np.sort(walls)
263 +
264 +
            # This sorts out the faces on a boundary and keeps the rest
265 +
            # Will be labelled walls
266 +
            walls_u, walls_index, walls_count = np.unique(
267 +
                walls_sorted, axis=0, return_counts=True, return_index=True
268 +
            )
269 +
            walls = walls[walls_index][walls_count == 1]
270 +
271 +
            string_top = stringify(top_faces)
272 +
            string_back = stringify(back_faces)
273 +
            string_bottom = stringify(bottom_faces)
274 +
            string_front = stringify(front_faces)
275 +
            string_left = stringify(left_faces)
276 +
            string_right = stringify(right_faces)
277 +
            string_walls = stringify(walls)
278 +
279 +
        # Expanding blocks matrix
280 +
        for i in range(len(centers)):
281 +
            blocks[i] = np.array([4, 5, 6, 7, 0, 1, 2, 3]) + 8 * i
282 +
283 +
        # Replaces duplicates vertices with their corresponding index in the
284 +
        # unique vertices list
285 +
        blocks = index[blocks]
286 +
287 +
        # OpenFOAM format
288 +
        vertices_unique = vertices_unique.reshape(-1, 3)
289 +
290 +
        # Creating the vertices string
291 +
        string_vert = str(list(vertices_unique)).replace("array([", "(")
292 +
        string_vert = string_vert.replace("])", ")\n\t\t")
293 +
        string_vert = "\t(" + string_vert.strip("[()]")
294 +
295 +
        # Creating all face strings - replaces [] with () for C++ format
296 +
297 +
        # Creating block string
298 +
        str_blocks = "[" + str(blocks).strip("[]") + "]"
299 +
        str_blocks = str_blocks.replace("[", "hex (")
300 +
        str_blocks = str_blocks.replace("]",
301 +
                                        ") (1 1 1) simpleGrading (1 1 1)\n\t")
302 +
303 +
        # Inserting vertices string
304 +
        str1 = "vertices\n(\n"
305 +
        index1 = file.find(str1) + len(str1)
306 +
        file = file[:index1] + "\t" + string_vert + file[index1:]
307 +
308 +
        # Inserting scale
309 +
        str2 = "scale "
310 +
        index2 = file.find(str2) + len(str2)
311 +
        file = file[:index2] + str(scale) + file[index2:]
312 +
313 +
        # Inserting faces
314 +
        if label:
315 +
            str3 = """
316 +
    top
317 +
    {
318 +
        type patch;
319 +
        faces
320 +
        (
321 +
"""
322 +
            index3 = file.find(str3) + len(str3)
323 +
            file = file[:index3] + str(string_top) + file[index3:]
324 +
325 +
            str4 = """
326 +
    back
327 +
    {
328 +
        type patch;
329 +
        faces
330 +
        (
331 +
"""
332 +
            index4 = file.find(str4) + len(str4)
333 +
            file = file[:index4] + str(string_back) + file[index4:]
334 +
335 +
            str5 = """
336 +
    bottom
337 +
    {
338 +
        type patch;
339 +
        faces
340 +
        (
341 +
"""
342 +
            index5 = file.find(str5) + len(str5)
343 +
            file = file[:index5] + str(string_bottom) + file[index5:]
344 +
345 +
            str6 = """
346 +
    front
347 +
    {
348 +
        type patch;
349 +
        faces
350 +
        (
351 +
"""
352 +
            index6 = file.find(str6) + len(str6)
353 +
            file = file[:index6] + str(string_front) + file[index6:]
354 +
355 +
            str7 = """
356 +
    left
357 +
    {
358 +
        type patch;
359 +
        faces
360 +
        (
361 +
"""
362 +
            index7 = file.find(str7) + len(str7)
363 +
            file = file[:index7] + str(string_left) + file[index7:]
364 +
365 +
            str8 = """
366 +
    right
367 +
    {
368 +
        type patch;
369 +
        faces
370 +
        (
371 +
"""
372 +
            index8 = file.find(str8) + len(str8)
373 +
            file = file[:index8] + str(string_right) + file[index8:]
374 +
375 +
            str9 = """
376 +
    walls
377 +
    {
378 +
        type wall;
379 +
        faces
380 +
        (
381 +
"""
382 +
383 +
            if string_walls != "()":
384 +
                index9 = file.find(str9) + len(str9)
385 +
                file = file[:index9] + str(string_walls) + file[index9:]
386 +
387 +
        # Deletes all empty labels if labeling is disabled
388 +
        if not label:
389 +
            string_start = "boundary\n("
390 +
            index_start = file.find(string_start) + len(string_start)
391 +
            string_end = "mergePatchPairs"
392 +
            index_end = file.find(string_end)
393 +
            file = file[:index_start] + "\n);\n\n" + file[index_end:]
394 +
395 +
        # Inserting blocks
396 +
        str10 = "blocks\n(\n"
397 +
        index9 = file.find(str10) + len(str10)
398 +
        file = file[:index9] + str(str_blocks) + file[index9:]
399 +
400 +
        # Gets rid of all commas in the file
401 +
        file = file.replace(",", "")
402 +
403 +
        # Returns the doc as a string
404 +
405 +
        with open("blockMeshDict", "w+") as a:
406 +
            return a.write(file)
407 +
408 +
409 +
def boundary(points, d, vert_unique):
410 +
    r""""""
411 +
    # Initialization
412 +
    boundary_points = np.zeros(len(points) * len(d) * 3).reshape(-1, 4, 3)
413 +
414 +
    # Creates all points at mesh boundary
415 +
    for i in range(len(points)):
416 +
        boundary_points[i] = points[i] + d
417 +
418 +
    boundary_points = boundary_points.reshape(-1, 3)
419 +
420 +
    # Creates tuples of each row and stores them in a numpy array
421 +
    # This allows for the vectorization of np.where
422 +
    vert_tuple = vert_unique.view(
423 +
        [('', vert_unique.dtype)] * vert_unique.shape[1]
424 +
    )
425 +
    boundary_tuple = boundary_points.view(
426 +
        [('', boundary_points.dtype)] * boundary_points.shape[1]
427 +
    )
428 +
    # Replaces vertices with corresponding index in vert_unique
429 +
    num = np.array([np.arange(0, len(vert_tuple))]).reshape(-1, 1)
430 +
    bool_mat = boundary_tuple[:, None] == vert_tuple
431 +
    faces = np.array([])
432 +
433 +
    # Replacment for np.where
434 +
    for i in bool_mat:
435 +
        faces = np.append(faces, num[i])
436 +
    faces = faces.astype(int).reshape(-1, 4)
437 +
438 +
    return faces
439 +
440 +
441 +
def stringify(thelist):
442 +
    r""""""
443 +
    string = str(thelist).replace("[", "(")
444 +
    string = string.replace("]", ")")
445 +
    string = "(" + string.strip("[()]") + ")"
446 +
    return string

@@ -290,3 +290,22 @@
Loading
290 290
        for j in range(3):
291 291
            export.vectors[i][j] = vertices[f[j], :]
292 292
    export.save(f"{filename}.stl")
293 +
294 +
295 +
def to_openfoam(im, scale=1, zoom_factor=1, label=True):
296 +
    r"""
297 +
    Save the image as an instruction file for OpenFoam to build an equivalent
298 +
    hexahedral mesh with (optionally) defined boundaries.
299 +
300 +
    Parameters
301 +
    ----------
302 +
    im : ND-array
303 +
        The image of the porous material
304 +
305 +
    filename : string or path object
306 +
        The name and location to save the file, which will have `.net` file
307 +
        extension.
308 +
309 +
    """
310 +
    from .openfoam import openfoam
311 +
    openfoam.save(im, scale, zoom_factor, label)

Learn more Showing 1 files with coverage changes found.

New file porespy/io/openfoam.py
New
Loading file...
Files Coverage
porespy -4.9% 83.3%
Project Totals (17 files) 83.3%
Loading