Create quarter administration

libthreedigrid/cells.f90

@@ -227,4 +227,190 @@ subroutine set_2d_computational_lines(l_u, l_v, k, m, n, mn, lg, lgrmin, area_ma
     end subroutine set_2d_computational_lines
 
 
+    subroutine set_quarter_admin(nodk, nodm, nodn, line, kcu, quarter_line, quarter_neighbour, liutot, livtot, n2dtot, n2dobc)
+
+        use parameters, only : LINE_2D_BOUNDARY_EAST, LINE_2D_BOUNDARY_WEST, LINE_2D_BOUNDARY_SOUTH, LINE_2D_BOUNDARY_NORTH
+
+        integer, intent(in) :: nodk(1:n2dtot+n2dobc)
+        integer, intent(in) :: nodm(1:n2dtot+n2dobc)
+        integer, intent(in) :: nodn(1:n2dtot+n2dobc)
+        integer, intent(in) :: line(1:liutot+livtot+n2dobc, 2)
+        integer, intent(in) :: kcu(1:liutot+livtot+n2dobc)
+        integer, intent(inout) :: quarter_line(1:4*n2dtot, 2)
+        integer, intent(inout) :: quarter_neighbour(1:4*n2dtot, 2)
+        integer, intent(in) :: liutot
+        integer, intent(in) :: livtot
+        integer, intent(in) :: n2dtot
+        integer, intent(in) :: n2dobc
+        integer :: l
+        integer :: quarter
+        integer :: nodd
+        integer :: nodu
+        integer :: kd
+        integer :: ku
+        integer :: md
+        integer :: mu
+        integer :: nd
+        integer :: nu
+        integer :: nb
+
+
+        !!!!!!!!!!!!!!!!
+        !!  3  !!  4  !!
+        !!!!!!!!!!!!!!!!
+        !!  1  !!  2  !!
+        !!!!!!!!!!!!!!!!
+        do l=1, liutot
+            nodd = line(l, 1) + 1 ! This is a 0-based python index, therefore +1
+            nodu = line(l, 2) + 1
+            kd = nodk(nodd)
+            ku = nodk(nodu)
+            md = nodm(nodd) 
+            mu = nodm(nodu)
+            nd = nodn(nodd)
+            nu = nodn(nodu)
+            if (ku==kd) then
+                quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1 ! This needs to be an index in Python again.
+                quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
+                quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
+                quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
+                quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
+                quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
+            elseif (ku == kd + 1) then
+                if (nd == 2 * nu - 1) then
+                    quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
+                endif
+                if (nd == 2 * nu) then
+                    quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
+                endif
+            else if (kd == ku + 1) then
+                if (nu == 2 * nd - 1) then
+                    quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
+                elseif (nu == 2 * nd) then
+                    quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 3), 1) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
+                endif
+            endif
+        enddo
+
+        do l=liutot + 1, liutot+livtot
+            nodd = line(l,1) + 1
+            nodu = line(l,2) + 1 
+            kd = nodk(nodd)
+            ku = nodk(nodu)
+            md = nodm(nodd) 
+            mu = nodm(nodu)
+            nd = nodn(nodd)
+            nu = nodn(nodu)
+            if (ku==kd) then
+                quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
+                quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
+                quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
+                quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
+                quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
+                quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
+            elseif (ku == kd + 1) then
+                if (md == 2 * mu - 1) then
+                    quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
+                endif
+                if (md == 2 * mu) then
+                    quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
+                endif
+            else if (kd == ku + 1) then
+                if (mu == 2 * md - 1) then
+                    quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
+                elseif (mu == 2 * md) then
+                    quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
+                    quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
+                    quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
+                    quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
+                endif
+            endif
+        enddo
+
+        do nb=1,n2dobc  !nodobc
+            l = liutot + livtot + nb
+            nodd = line(l, 1) + 1
+            nodu = line(l, 2) + 1
+            select case(kcu(l))
+            case(LINE_2D_BOUNDARY_WEST)
+                quarter_line(get_quarter_idx(nodu, 1), 1) = l - 1
+                quarter_line(get_quarter_idx(nodu, 3), 1) = l
+                quarter_neighbour(get_quarter_idx(nodu, 1), 1) = nodd - 1
+                quarter_neighbour(get_quarter_idx(nodu, 3), 1) = nodd - 1
+            case(LINE_2D_BOUNDARY_EAST)
+                quarter_line(get_quarter_idx(nodd, 2), 1) = l - 1
+                quarter_line(get_quarter_idx(nodd, 4), 1) = l - 1
+                quarter_neighbour(get_quarter_idx(nodd, 2), 1) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodd, 4), 1) = nodu - 1
+            case(LINE_2D_BOUNDARY_SOUTH)
+                quarter_line(get_quarter_idx(nodu, 1), 2) = l - 1
+                quarter_line(get_quarter_idx(nodu, 2), 2) = l - 1
+                quarter_neighbour(get_quarter_idx(nodu, 1), 2) = nodd - 1
+                quarter_neighbour(get_quarter_idx(nodu, 2), 2) = nodd - 1
+            case(LINE_2D_BOUNDARY_NORTH)
+                quarter_line(get_quarter_idx(nodd, 3), 2) = l - 1
+                quarter_line(get_quarter_idx(nodd, 4), 2) = l - 1
+                quarter_neighbour(get_quarter_idx(nodd, 3), 2) = nodu - 1
+                quarter_neighbour(get_quarter_idx(nodd, 4), 2) = nodu - 1
+            end select
+        enddo
+
+    end subroutine set_quarter_admin
+
+
+    function get_quarter_idx(nod, quadrant) result(idx)
+        !!!!!!!!!!!!!!!!
+        !!  3  !!  4  !!
+        !!!!!!!!!!!!!!!!
+        !!  1  !!  2  !!
+        !!!!!!!!!!!!!!!!
+        integer, intent(in) :: nod
+        integer, intent(in) :: quadrant
+        integer :: idx
+
+        idx = 4 * (nod - 1) + quadrant
+
+    end function get_quarter_idx
+
 end module m_cells

libthreedigrid/parameters.f90

@@ -13,5 +13,31 @@ module parameters
                           NODE_2D_BOUND = 5,&
                           NODE_2D_GW_BOUND = 6,&
                           NODE_1D_BOUND = 7
+
+    integer, parameter :: LINE_1D_EMBEDDED = 0,&
+                          LINE_1D_ISOLATED = 1,&
+                          LINE_1D_CONNECTED = 2,&
+                          LINE_1D_LONG_CRESTED = 3,&
+                          LINE_1D_SHORT_CRESTED = 4,&
+                          LINE_1D_DOUBLE_CONNECTED = 5,&
+                          LINE_1D2D_SINGLE_CONNECTED_CLOSED = 51,&
+                          LINE_1D2D_SINGLE_CONNECTED_OPEN_WATER = 52,&
+                          LINE_1D2D_DOUBLE_CONNECTED_CLOSED = 53,&
+                          LINE_1D2D_DOUBLE_CONNECTED_OPEN_WATER = 54,&
+                          LINE_1D2D_POSSIBLE_BREACH = 55,&
+                          LINE_1D2D_ACTIVE_BREACH = 56,&
+                          LINE_1D2D_GROUNDWATER = 57,&
+                          LINE_2D = 100,&
+                          LINE_2D_OBSTACLE = 101,&
+                          LINE_2D_VERTICAL = 150,&
+                          LINE_2D_GROUNDWATER = -150,&
+                          LINE_2D_BOUNDARY_WEST = 200,&
+                          LINE_2D_BOUNDARY_EAST = 300,&
+                          LINE_2D_BOUNDARY_SOUTH = 400,&
+                          LINE_2D_BOUNDARY_NORTH = 500,&
+                          LINE_2D_GROUNDWATER_BOUNDARY_WEST = 600,&
+                          LINE_2D_GROUNDWATER_BOUNDARY_EAST = 700,&
+                          LINE_2D_GROUNDWATER_BOUNDARY_SOUTH = 800,&
+                          LINE_2D_GROUNDWATER_BOUNDARY_NORTH = 900
 
 end module parameters

threedigrid_builder/application.py

@@ -93,6 +93,9 @@ def _make_gridadmin(
             node_id_counter,
             line_id_counter,
         )
+
+        grid.set_quarter_administration(quadtree)
+
         dem_average_areas = db.get_dem_average_areas()
         grid.set_dem_averaged_cells(dem_average_areas)
 
@@ -225,6 +228,7 @@ def make_gridadmin(
     progress_callback: Optional[Callable[[float, str], None]] = None,
     upgrade: bool = False,
     convert_to_geopackage: bool = False,
+    create_triangulation: bool = False,
 ):
     """Create a Grid instance from sqlite and DEM paths
 

threedigrid_builder/base/__init__.py

@@ -5,5 +5,6 @@
 from .linestrings import *  # NOQA
 from .nodes import *  # NOQA
 from .pumps import *  # NOQA
+from .quarters import *  # NOQA
 from .settings import *  # NOQA
 from .surfaces import *  # NOQA

threedigrid_builder/base/quarters.py

@@ -0,0 +1,15 @@
+from typing import Tuple
+from .array import Array
+
+__all__ = ["Quarters"]
+
+
+class Quarter:
+    id: int
+    line: Tuple[int, int]
+    neighbour_node: Tuple[int, int]
+
+
+class Quarters(Array[Quarter]):
+    """Calculation quarters."""
+    pass

threedigrid_builder/grid/grid.py

@@ -11,6 +11,7 @@
 from threedigrid_builder.base import (
     Lines,
     Nodes,
+    Quarters,
     PointsOnLine,
     Pumps,
     SurfaceMaps,
@@ -172,6 +173,7 @@ def __init__(
         breaches=None,
         meta=None,
         quadtree_stats=None,
+        quarters: Optional[Quarters] = None,
     ):
         if not isinstance(nodes, Nodes):
             raise TypeError(f"Expected Nodes instance, got {type(nodes)}")
@@ -206,8 +208,15 @@ def __init__(
             raise TypeError(f"Expected Obstacles instance, got {type(obstacles)}")
         if breaches is not None and not isinstance(breaches, PotentialBreaches):
             raise TypeError(f"Expected Breaches instance, got {type(breaches)}")
+
+        if quarters is None:
+            quarters = Quarters(id=[])
+        elif not isinstance(quarters, Quarters):
+            raise TypeError(f"Expected Quarters instance, got {type(quarters)}")
+
         self.nodes = nodes
         self.lines = lines
+        self.quarters = quarters
         self.meta = meta
         self.surfaces = surfaces
         self.surface_maps = surface_maps
@@ -605,6 +614,11 @@ def set_boundary_conditions_2d(
         self.nodes += nodes
         self.lines += lines
 
+    def set_quarter_administration(self, quadtree):
+        """Sets the quarter cell administration for lines and neighbouring calculation cells.
+        """
+        self.quarters = quadtree.get_quarters_admin(self.nodes, self.lines)
+
     def set_pumps(self, pumps):
         """Set the pumps on this grid object
 

threedigrid_builder/grid/quadtree.py

@@ -4,7 +4,7 @@
 
 import numpy as np
 
-from threedigrid_builder.base import Lines, Nodes
+from threedigrid_builder.base import Lines, Nodes, Quarters
 from threedigrid_builder.constants import LineType, NodeType
 from threedigrid_builder.exceptions import SchematisationError
 
@@ -247,3 +247,38 @@ def get_nodes_lines(self, area_mask, node_id_counter, line_id_counter):
         )
 
         return nodes, lines
+
+    def get_quarters_admin(self, nodes, lines):
+        """Compute 2D openwater quarter administration based on computed Quadtree.
+        Determine the connecting flowlines for a quarter in u and v direction.
+        Determine the adjacent calculation nodes for a quarter in u and v direction.
+
+        Return:
+          Quarters object
+        """
+
+        n_2d_nodes = np.count_nonzero(np.isin(nodes.node_type, NodeType.NODE_2D_OPEN_WATER))
+        n_2d_bnd_nodes = np.count_nonzero(np.isin(nodes.node_type, NodeType.NODE_2D_BOUNDARIES))\
+
+        quarter_line = np.full((n_2d_nodes * 4, 2), -9999, dtype=np.int32, order="F")
+        neighbour_node = np.full((n_2d_nodes * 4, 2), -9999, dtype=np.int32, order="F")
+
+        m_cells.set_quarter_admin(
+            nodes.nodk,
+            nodes.nodm,
+            nodes.nodn,
+            lines.line,
+            lines.kcu,
+            quarter_line,
+            neighbour_node,
+            self.n_lines_u,
+            self.n_lines_v,
+            n_2d_nodes,
+            n_2d_bnd_nodes,
+        )
+
+        return Quarters(
+            id=np.arange(0, 4 * n_2d_nodes),
+            line=quarter_line,
+            neighbour_node=neighbour_node
+        )

threedigrid_builder/interface/gridadmin.py

@@ -172,6 +172,7 @@ def write(self, grid: Grid):
         self.write_quadtree(grid.quadtree_stats)
         self.write_nodes(grid.nodes)
         self.write_nodes_embedded(grid.nodes_embedded)
+        self.write_quarters(grid.quarters)
         self.write_lines(grid.lines, grid.cross_sections)
         self.write_pumps(grid.pumps)
         self.write_cross_sections(grid.cross_sections)
@@ -558,6 +559,25 @@ def write_lines(self, lines: Lines, cross_sections):
         )
         self.write_dataset(group, "windshieldings", lines.windshieldings.T)
 
+    def write_quarters(self, quarters, group_name="quarters"):
+        """Write the "quarters" group in the gridadmin file
+
+        Raises a ValueError if it exists already.
+
+        Notes:
+            Some datasets were 64-bit integers, but now they are saved as 32-bit integers.
+
+        Args:
+            quarters (Quarters)
+        """
+        if quarters is None:
+            return
+        group = self._file.create_group(group_name)
+
+        self.write_dataset(group, "id", quarters.id + 1)
+        self.write_dataset(group, "quarters_line", quarters.line.T + 1)
+        self.write_dataset(group, "quarters_neighbour_node", quarters.neighbour_node.T + 1)
+
     def write_pumps(self, pumps):
         group = self._file.create_group("pumps")