tnode_types.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tnode_types.cc (4058B)
---
1 /* Copyright (C) 2016, 2017, 2018 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include "node_types.hh"
21
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/error_handling.hh"
25
26 namespace pism {
27
28 /**
29 Identify node types of all the nodes (grid points).
30
31 Uses width=1 ghosts of `ice_thickness` (a "box" stencil).
32
33 A node is can be *interior*, *boundary*, or *exterior*.
34
35 An element is considered *icy* if all four of its nodes have ice thickness above
36 `thickness_threshold`.
37
38 A node is considered *interior* if all of the elements it belongs to are icy.
39
40 A node is considered *boundary* if it is not interior and at least one element it belongs to is
41 icy.
42
43 A node is considered *exterior* if it is neither interior nor boundary.
44
45 Now an element "face" (side) is a part of a boundary if and only if both of its nodes are
46 boundary nodes.
47
48 Cell layout:
49 ~~~
50 (i-1,j+1) +-------N--------+ (i+1,j+1)
51 | | |
52 | NW | NE |
53 | | |
54 W-----(i,j)------E
55 | | |
56 | SW | SE |
57 | | |
58 (i-1,j-1) +-------S--------+ (i+1,j-1)
59 ~~~
60 */
61 void compute_node_types(const IceModelVec2S &ice_thickness,
62 double thickness_threshold,
63 IceModelVec2Int &result) {
64
65 IceGrid::ConstPtr grid = ice_thickness.grid();
66
67 const IceModelVec2S &H = ice_thickness;
68 const double &H_min = thickness_threshold;
69
70 IceModelVec::AccessList list{&H, &result};
71
72 ParallelSection loop(grid->com);
73 try {
74 for (Points p(*grid); p; p.next()) {
75 const int i = p.i(), j = p.j();
76
77 // indexing shortcuts (to reduce chances of making typos below)
78 const int
79 N = j + 1,
80 E = i + 1,
81 S = j - 1,
82 W = i - 1;
83
84 // flags indicating whether the current node and its neighbors are "icy"
85 const bool
86 icy_ij = H(i, j) >= H_min,
87 icy_nw = H(W, N) >= H_min,
88 icy_n = H(i, N) >= H_min,
89 icy_ne = H(E, N) >= H_min,
90 icy_e = H(E, j) >= H_min,
91 icy_se = H(E, S) >= H_min,
92 icy_s = H(i, S) >= H_min,
93 icy_sw = H(W, S) >= H_min,
94 icy_w = H(W, j) >= H_min;
95
96 // flags indicating whether neighboring elements are "icy" (and element is icy if all its
97 // nodes are icy)
98 const bool
99 ne_element_is_icy = (icy_ij and icy_e and icy_ne and icy_n),
100 nw_element_is_icy = (icy_ij and icy_n and icy_nw and icy_w),
101 sw_element_is_icy = (icy_ij and icy_w and icy_sw and icy_s),
102 se_element_is_icy = (icy_ij and icy_s and icy_se and icy_e);
103
104 if (ne_element_is_icy and nw_element_is_icy and
105 sw_element_is_icy and se_element_is_icy) {
106 // all four elements are icy: we are at an interior node
107 result(i, j) = NODE_INTERIOR;
108 } else if (ne_element_is_icy or nw_element_is_icy or
109 sw_element_is_icy or se_element_is_icy) {
110 // at least one element is icy: we are at a boundary node
111 result(i, j) = NODE_BOUNDARY;
112 } else {
113 // all elements are ice-free: we are at an exterior node
114 result(i, j) = NODE_EXTERIOR;
115 }
116
117 }
118 } catch (...) {
119 loop.failed();
120 }
121 loop.check();
122
123 result.update_ghosts();
124 }
125
126 } // end of namespace pism