FireDynamics / fdsreader

Python reader for FDS data
GNU General Public License v3.0
44 stars 18 forks source link

Get obstruction mask for all meshes #43

Closed fliszer closed 1 year ago

fliszer commented 2 years ago

Hi,

I'm wondering if it is possible to get obst mask for silce in multimesh case? I also get error for single mesh (in multimesh case). Below I tried to use your code but without success:

def main():
    sim = fds.Simulation("my_path_to_sim_dir")

    # Get the first mesh defined in fds file
    mesh = sim.meshes[0]
    # Get the slice
    slc = sim.slices[1]

    # Get subslice that cuts through our mesh
    subslice = slc[mesh]

    # Mask the data
    mask = mesh.get_obstruction_mask_slice(subslice) # <---------- here I get error - details below

    # Timestep
    t = -1
    # Fill value for mask
    fill = 0
    sslc_data = np.where(mask[t], subslice.data[t], fill)

    # Plot the slice
    plt.imshow(sslc_data.T, origin="lower")
    plt.colorbar()
    plt.show()

And the error details:

Traceback (most recent call last):
  File "C:/Users/fliszer/PycharmProjects/PMV-PPD/fds_read_tests.py", line 40, in <module>
    main()
  File "C:/Users/fliszer/PycharmProjects/PMV-PPD/fds_read_tests.py", line 24, in main
    mask = mesh.get_obstruction_mask_slice(subslice)
  File "C:\Users\fliszer\PycharmProjects\PMV-PPD\venv\lib\site-packages\fdsreader\fds_classes\mesh.py", line 80, in get_obstruction_mask_slice
    return np.squeeze(self.get_obstruction_mask(subslice.times, cell_centered=cell_centered)[mask_indices])
  File "C:\Users\fliszer\PycharmProjects\PMV-PPD\venv\lib\site-packages\fdsreader\fds_classes\mesh.py", line 60, in get_obstruction_mask
    for t, _ in enumerate(subobst.visible_times(times)):
  File "C:\Users\fliszer\PycharmProjects\PMV-PPD\venv\lib\site-packages\fdsreader\bndf\obstruction.py", line 292, in visible_times
    for time in self.times:
  File "C:\Users\fliszer\PycharmProjects\PMV-PPD\venv\lib\site-packages\fdsreader\bndf\obstruction.py", line 283, in times
    return next(iter(self._boundary_data.values())).times
StopIteration
JanVogelsang commented 2 years ago

Hi, thanks for posting your bug report and feature request.

In version 1.8.2 the bug you are experiencing should be fixed, it was an edge case I didn't consider that is now covered.

Regarding your request to get a slice obstruction mask for multimesh cases, I added the option to get the global slice (using the to_global method") with applied obstruction masked using the "masked=True" parameter and optional fill value parameter (e.g. "fill=50"). Adding a method to get the mask as an independent array would need substantial work, and I currently don't have the time for that. Tell me if the solution I provided is sufficient for your needs.

fliszer commented 1 year ago

Hi,

thanks for feedback but I got another error according to arrays sizes:

Traceback (most recent call last):
  File "C:/Users/fliszer/PycharmProjects/PMV-PPD/fds_read_tests.py", line 22, in <module>
    main()
  File "C:/Users/fliszer/PycharmProjects/PMV-PPD/fds_read_tests.py", line 13, in main
    data = slc.to_global(masked=True, fill=100)
  File "C:\Users\fliszer\PycharmProjects\PMV-PPD\venv\lib\site-packages\fdsreader\slcf\slice.py", line 595, in to_global
    slc_data = np.where(mask, slc_data, fill)
  File "<__array_function__ internals>", line 6, in where
ValueError: operands could not be broadcast together with shapes (201,84,1,79) (201,83,1,79) ()

I did some debugging and I found out that variables in slice.py file in to_global method not equals:

print(mesh.coordinates[dim][-1])
print(global_max[dim])

so this and previous condition (if) not works (lines 559, 569):

# Add border points back again if needed
if not self.cell_centered and mesh.coordinates[dim][-1] == global_max[dim]:
    slc_data = np.concatenate((slc_data, temp_data), axis=axis + 1)

In my case I have below meshes"

# ---- Mesh ----
&MESH ID='MESH1', COLOR='YELLOW', IJK=83,59,78, XB=159.6,184.5, -62.4,-44.7, -0.3,23.1 /
&MESH ID='MESH2', COLOR='YELLOW', IJK=83,59,78, XB=159.6,184.5, -44.7,-27, -0.3,23.1 /
&MESH ID='MESH3', COLOR='YELLOW', IJK=83,59,78, XB=159.6,184.5, -27,-9.3, -0.3,23.1 /
&MESH ID='MESH4', COLOR='YELLOW', IJK=83,53,78, XB=159.6,184.5, -9.3,6.6, -0.3,23.1 /
&MESH ID='MESH5', COLOR='YELLOW', IJK=83,59,78, XB=184.5,209.4, -62.4,-44.7, -0.3,23.1 /
&MESH ID='MESH6', COLOR='YELLOW', IJK=83,59,78, XB=184.5,209.4, -44.7,-27, -0.3,23.1 /
&MESH ID='MESH7', COLOR='YELLOW', IJK=83,59,78, XB=184.5,209.4, -27,-9.3, -0.3,23.1 /
&MESH ID='MESH8', COLOR='YELLOW', IJK=83,53,78, XB=184.5,209.4, -9.3,6.6, -0.3,23.1 /
&MESH ID='MESH9', COLOR='YELLOW', IJK=83,59,78, XB=209.4,234.3, -62.4,-44.7, -0.3,23.1 /
&MESH ID='MESH10', COLOR='YELLOW', IJK=83,59,78, XB=209.4,234.3, -44.7,-27, -0.3,23.1 /
&MESH ID='MESH11', COLOR='YELLOW', IJK=83,59,78, XB=209.4,234.3, -27,-9.3, -0.3,23.1 /
&MESH ID='MESH12', COLOR='YELLOW', IJK=83,53,78, XB=209.4,234.3, -9.3,6.6, -0.3,23.1 /
&MESH ID='MESH13', COLOR='YELLOW', IJK=75,59,78, XB=234.3,256.8, -62.4,-44.7, -0.3,23.1 /
&MESH ID='MESH14', COLOR='YELLOW', IJK=75,59,78, XB=234.3,256.8, -44.7,-27, -0.3,23.1 /
&MESH ID='MESH15', COLOR='YELLOW', IJK=75,59,78, XB=234.3,256.8, -27,-9.3, -0.3,23.1 /
&MESH ID='MESH16', COLOR='YELLOW', IJK=75,53,78, XB=234.3,256.8, -9.3,6.6, -0.3,23.1 /

so you can see that first mesh max x coordinate is lower than global max x coordinate. I think this is the clue of the problem. When I removed if conditions code works well but I know it's not a solution :)