Skip to content

Commit

Permalink
uncomment face_nodes
Browse files Browse the repository at this point in the history
  • Loading branch information
@veenstrajelmer
veenstrajelmer committed Oct 30, 2023
2 parents ba7ea0c + f484b9f commit a631aff
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 40 deletions.
75 changes: 43 additions & 32 deletions hydrolib/core/dflowfm/net/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,9 @@ def split_by(gl: mk.GeometryList, by: float) -> list:
return lists


class Mesh2d(BaseModel): #TODO: this is an inconvenient name since meshkernel also has a Mesh2d class
class Mesh2d(
BaseModel
): # TODO: this is an inconvenient name since meshkernel also has a Mesh2d class
"""Mesh2d defines a single two dimensional grid.
Attributes:
Expand Down Expand Up @@ -137,14 +139,14 @@ def read_file(self, file_path: Path) -> None:
reader = UgridReader(file_path)
reader.read_mesh2d(self)

def _set_mesh2d(self) -> None:
mesh2d = mk.Mesh2d(
node_x=self.mesh2d_node_x,
node_y=self.mesh2d_node_y,
edge_nodes=self.mesh2d_edge_nodes.ravel(),
)
# def _set_mesh2d(self) -> None:
# mesh2d = mk.Mesh2d(
# node_x=self.mesh2d_node_x,
# node_y=self.mesh2d_node_y,
# edge_nodes=self.mesh2d_edge_nodes.ravel(),
# )

self.meshkernel.mesh2d_set(mesh2d)
# self.meshkernel.mesh2d_set(mesh2d)

def get_mesh2d(self) -> mk.Mesh2d:
"""Get the mesh2d as represented in the MeshKernel
Expand All @@ -167,24 +169,26 @@ def create_rectilinear(self, extent: tuple, dx: float, dy: float) -> None:
"""

xmin, ymin, xmax, ymax = extent

rows = int((ymax - ymin) / dy)
columns = int((xmax - xmin) / dx)

params = mk.MakeGridParameters(num_columns=columns,
num_rows=rows,
origin_x=xmin,
origin_y=ymin,
block_size_x=dx,
block_size_y=dy)

mesh2d_input = self.meshkernel #mk.MeshKernel()

params = mk.MakeGridParameters(
num_columns=columns,
num_rows=rows,
origin_x=xmin,
origin_y=ymin,
block_size_x=dx,
block_size_y=dy,
)

mesh2d_input = self.meshkernel # mk.MeshKernel()
mesh2d_input.curvilinear_compute_rectangular_grid(params)
mesh2d_input.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
mesh2d_input_m2d = mesh2d_input.mesh2d_get() #get Mesh2d object
mesh2d_input.curvilinear_convert_to_mesh2d() # convert to ugrid/mesh2d
mesh2d_input_m2d = mesh2d_input.mesh2d_get() # get Mesh2d object

# mesh2d_input_raw = mk.Mesh2d(mesh2d_input_m2d.node_x, mesh2d_input_m2d.node_y, mesh2d_input_m2d.edge_nodes)

# Process
self._process(mesh2d_input_m2d)

Expand All @@ -201,12 +205,14 @@ def create_triangular(self, geometry_list: mk.GeometryList) -> None:
self._process(self.get_mesh2d())

def _process(self, mesh2d_input) -> None:

# Add input
# self.meshkernel.mesh2d_set(mesh2d_input) #TODO: in this meshkernel function duplicates the amount of nodes. Seems not desireable, but more testbanks fail if commented.
# Get output
mesh2d_output = self.meshkernel.mesh2d_get() #better results for some testcases, comment above and: mesh2d_output = mesh2d_input

mesh2d_output = (
self.meshkernel.mesh2d_get()
) # better results for some testcases, comment above and: mesh2d_output = mesh2d_input

# Add to mesh2d variables
self.mesh2d_node_x = mesh2d_output.node_x
self.mesh2d_node_y = mesh2d_output.node_y
Expand All @@ -226,7 +232,6 @@ def _process(self, mesh2d_input) -> None:
np.ones_like(self.mesh2d_face_nodes) * np.arange(max(npf))[None, :]
) < npf[:, None]
self.mesh2d_face_nodes[idx] = mesh2d_output.face_nodes


def clip(
self,
Expand All @@ -243,7 +248,7 @@ def clip(

# Add current mesh to Mesh2d instance
# self._set_mesh2d()

# For clipping outside
if not inside:
# Check if a multipolygon was provided when clipping outside
Expand Down Expand Up @@ -639,7 +644,9 @@ def clear(self) -> None:
self.link1d2d = np.empty((0, 2), np.int32)
# The meshkernel object needs to be resetted
self.meshkernel._deallocate_state()
self.meshkernel._allocate_state(self.meshkernel.is_geographic) #TODO: .get_projection() might be better
self.meshkernel._allocate_state(
self.meshkernel.is_geographic
) # TODO: .get_projection() might be better
self.meshkernel.contacts_get()

def _process(self) -> None:
Expand Down Expand Up @@ -681,7 +688,9 @@ def _link_from_1d_to_2d(

# Computes Mesh1d-Mesh2d contacts, where each single Mesh1d node is connected to one Mesh2d face circumcenter.
# The boundary nodes of Mesh1d (those sharing only one Mesh1d edge) are not connected to any Mesh2d face.
self.meshkernel.contacts_compute_single(node_mask=node_mask, polygons=polygon, projection_factor=1.0)
self.meshkernel.contacts_compute_single(
node_mask=node_mask, polygons=polygon, projection_factor=1.0
)
self._process()

# Note that the function "contacts_compute_multiple" also computes the connections, but does not take into account
Expand Down Expand Up @@ -1114,15 +1123,17 @@ def get_node_mask(self, branchids: List[str] = None):


class Network:
def __init__(self, projection: mk.ProjectType = mk.ProjectionType.CARTESIAN) -> None:
def __init__(
self, projection: mk.ProjectType = mk.ProjectionType.CARTESIAN
) -> None:
self.meshkernel = mk.MeshKernel(projection=projection)
# Monkeypatch the meshkernel object, because the "is_geographic" is not saved
# otherwise, and needed for reinitializing the meshkernel
if projection == mk.ProjectionType.CARTESIAN:
is_geographic = False
else: #SPHERICAL or SPHERICALACCURATE
else: # SPHERICAL or SPHERICALACCURATE
is_geographic = True
#TODO: this bool does not seems to have effect, maybe discontinue?
# TODO: this bool does not seems to have effect, maybe discontinue?
self.meshkernel.is_geographic = is_geographic

self._mesh1d = Mesh1d(meshkernel=self.meshkernel)
Expand Down
17 changes: 9 additions & 8 deletions tests/dflowfm/test_net.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ def test_create_2d():
mesh2d = Mesh2d(meshkernel=MeshKernel())
mesh2d.create_rectilinear(extent=bbox, dx=0.5, dy=0.75)

#TODO: remove plotting
# TODO: remove plotting
# _, ax = plt.subplots()
# mesh2d_output.plot_edges(ax=ax)

Expand All @@ -188,7 +188,7 @@ def test_create_2d():
],
)
def test_create_clip_2d(deletemeshoption, inside, nnodes, nedgenodes, nfaces):

polygon = GeometryList(
x_coordinates=np.array([0.0, 6.0, 4.0, 2.0, 0.0]),
y_coordinates=np.array([0.0, 2.0, 7.0, 6.0, 0.0]),
Expand All @@ -202,10 +202,10 @@ def test_create_clip_2d(deletemeshoption, inside, nnodes, nedgenodes, nfaces):
mesh2d.clip(polygon, deletemeshoption=deletemeshoption, inside=inside)
mesh2d_output = mesh2d.get_mesh2d()
assert mesh2d_output.node_x.size == nnodes
assert mesh2d_output.edge_nodes.size == nedgenodes # 2x nedges
assert mesh2d_output.edge_nodes.size == nedgenodes # 2x nedges
assert mesh2d_output.face_x.size == nfaces

#TODO: remove plotting
# TODO: remove plotting
# _, ax = plt.subplots()
# mesh2d_output.plot_edges(ax=ax)
# ax.plot(polygon.x_coordinates, polygon.y_coordinates, "r-")
Expand All @@ -231,8 +231,8 @@ def test_create_refine_2d():

assert mesh2d_output.node_x.size == 114
assert mesh2d_output.edge_nodes.size == 426
#TODO: remove plotting

# TODO: remove plotting
# _, ax = plt.subplots()
# mesh2d_output = mesh2d.get_mesh2d()
# mesh2d_output.plot_edges(ax=ax)
Expand All @@ -243,6 +243,7 @@ def test_create_refine_2d():
test_input_dir / "e02/c11_korte-woerden-1d/dimr_model/dflowfm/FlowFM_net.nc",
]


@pytest.mark.parametrize(
"filepath",
cases,
Expand Down Expand Up @@ -584,7 +585,7 @@ def test_create_triangular():
)

network.mesh2d_create_triangular_within_polygon(polygon)

assert np.array_equiv(
network._mesh2d.mesh2d_node_x,
np.array([6.0, 4.0, 2.0, 0.0]),
Expand All @@ -597,7 +598,7 @@ def test_create_triangular():
network._mesh2d.mesh2d_edge_nodes,
np.array([[2, 3], [3, 0], [0, 2], [0, 1], [1, 2]]),
)


def test_add_1d2d_links():

Expand Down

0 comments on commit a631aff

Please sign in to comment.