Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
mpim-sw
cdo
Commits
257074b2
Commit
257074b2
authored
Oct 07, 2020
by
Uwe Schulzweida
Browse files
verifygrid: changed datatype to Point3D.
parent
d0a60fa2
Pipeline
#4532
passed with stages
in 16 minutes and 32 seconds
Changes
3
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
src/Samplegridicon.cc
View file @
257074b2
...
...
@@ -316,8 +316,8 @@ compute_child_from_bounds(CellIndex *cellindex2, long ncells2, double *grid_cent
size_t
*
nbr_addr
=
&
knnWeights
.
m_addr
[
0
];
int
ncorner
=
3
;
double
center_point_xyz
[
3
]
;
double
cell_corners_xyz
[
12
]
;
Point3D
center_point_xyz
;
Varray
<
Point3D
>
cell_corners_xyz
(
4
)
;
Point
center_point_plane_projection
;
Varray
<
Point
>
cell_corners_plane_projection
(
4
);
...
...
@@ -328,11 +328,13 @@ compute_child_from_bounds(CellIndex *cellindex2, long ncells2, double *grid_cent
for
(
int
k
=
0
;
k
<
MAX_CHILDS
;
++
k
)
child2
[
cell_no2
*
MAX_CHILDS
+
k
]
=
-
1
;
set_cell_corners_xyz
(
ncorner
,
&
grid_corner_lon2
[
cell_no2
*
ncorner
],
&
grid_corner_lat2
[
cell_no2
*
ncorner
],
cell_corners_xyz
);
cell_corners_xyz
[
ncorner
]
=
cell_corners_xyz
[
0
];
const
auto
coordinate_to_ignore
=
find_coordinate_to_ignore
(
cell_corners_xyz
);
bool
invert_result
=
false
;
if
(
cell_corners_xyz
[
coordinate_to_ignore
-
1
]
<
0
)
invert_result
=
true
;
const
auto
cval
=
(
coordinate_to_ignore
==
1
)
?
cell_corners_xyz
[
0
].
X
:
((
coordinate_to_ignore
==
2
)
?
cell_corners_xyz
[
0
].
Y
:
cell_corners_xyz
[
0
].
Z
);
if
(
cval
<
0.0
)
invert_result
=
true
;
set_cell_corners_plane_projection
(
coordinate_to_ignore
,
ncorner
,
cell_corners_xyz
,
cell_corners_plane_projection
);
...
...
@@ -345,12 +347,16 @@ compute_child_from_bounds(CellIndex *cellindex2, long ncells2, double *grid_cent
gridSearchPoint
(
gps
,
grid_center_lon2
[
cell_no2
],
grid_center_lat2
[
cell_no2
],
knnWeights
);
int
k
=
0
;
double
center_coordinates
[
3
];
for
(
int
i
=
0
;
i
<
MAX_SEARCH
;
++
i
)
{
const
auto
cell_no1
=
nbr_addr
[
i
];
if
(
cell_no1
<
SIZE_MAX
)
{
gcLLtoXYZ
(
grid_center_lon1
[
cell_no1
],
grid_center_lat1
[
cell_no1
],
center_point_xyz
);
gcLLtoXYZ
(
grid_center_lon1
[
cell_no1
],
grid_center_lat1
[
cell_no1
],
center_coordinates
);
center_point_xyz
.
X
=
center_coordinates
[
0
];
center_point_xyz
.
Y
=
center_coordinates
[
1
];
center_point_xyz
.
Z
=
center_coordinates
[
2
];
set_center_point_plane_projection
(
coordinate_to_ignore
,
center_point_xyz
,
center_point_plane_projection
);
...
...
src/Verifygrid.cc
View file @
257074b2
...
...
@@ -119,13 +119,13 @@ determinant(const double (&matrix)[3][3])
}
static
void
find_unit_normal
(
const
double
a
[
3
],
const
double
b
[
3
],
const
double
c
[
3
],
double
(
&
unit_normal
)
[
3
])
find_unit_normal
(
const
Point3D
&
a
,
const
Point3D
&
b
,
const
Point3D
&
c
,
Point3D
&
unit_normal
)
{
// Calculates the unit normal for a plane defined on three points a, b, c in Euclidean space.
const
double
matrix_for_x
[
3
][
3
]
=
{
{
1
,
a
[
1
],
a
[
2
]
},
{
1
,
b
[
1
],
b
[
2
]
},
{
1
,
c
[
1
],
c
[
2
]
}
};
const
double
matrix_for_y
[
3
][
3
]
=
{
{
a
[
0
]
,
1
,
a
[
2
]
},
{
b
[
0
]
,
1
,
b
[
2
]
},
{
c
[
0
]
,
1
,
c
[
2
]
}
};
const
double
matrix_for_z
[
3
][
3
]
=
{
{
a
[
0
],
a
[
1
]
,
1
},
{
b
[
0
],
b
[
1
]
,
1
},
{
c
[
0
],
c
[
1
]
,
1
}
};
const
double
matrix_for_x
[
3
][
3
]
=
{
{
1
,
a
.
Y
,
a
.
Z
},
{
1
,
b
.
Y
,
b
.
Z
},
{
1
,
c
.
Y
,
c
.
Z
}
};
const
double
matrix_for_y
[
3
][
3
]
=
{
{
a
.
X
,
1
,
a
.
Z
},
{
b
.
X
,
1
,
b
.
Z
},
{
c
.
X
,
1
,
c
.
Z
}
};
const
double
matrix_for_z
[
3
][
3
]
=
{
{
a
.
X
,
a
.
Y
,
1
},
{
b
.
X
,
b
.
Y
,
1
},
{
c
.
X
,
c
.
Y
,
1
}
};
const
auto
x
=
determinant
(
matrix_for_x
);
const
auto
y
=
determinant
(
matrix_for_y
);
...
...
@@ -133,28 +133,28 @@ find_unit_normal(const double a[3], const double b[3], const double c[3], double
const
auto
magnitude
=
std
::
sqrt
(
x
*
x
+
y
*
y
+
z
*
z
);
unit_normal
[
0
]
=
x
/
magnitude
;
unit_normal
[
1
]
=
y
/
magnitude
;
unit_normal
[
2
]
=
z
/
magnitude
;
unit_normal
.
X
=
x
/
magnitude
;
unit_normal
.
Y
=
y
/
magnitude
;
unit_normal
.
Z
=
z
/
magnitude
;
}
int
find_coordinate_to_ignore
(
const
double
*
cell_corners_xyz
)
find_coordinate_to_ignore
(
const
Varray
<
Point3D
>
&
cell_corners_xyz
)
{
// Takes the first three corners/vertices of the cell and calculates the unit normal via determinants.
const
auto
pcorner_coordinates1
=
&
cell_corners_xyz
[
0
];
const
auto
pcorner_coordinates2
=
&
cell_corners_xyz
[
3
];
const
auto
pcorner_coordinates3
=
&
cell_corners_xyz
[
6
];
const
auto
&
pcorner_coordinates1
=
cell_corners_xyz
[
0
];
const
auto
&
pcorner_coordinates2
=
cell_corners_xyz
[
1
];
const
auto
&
pcorner_coordinates3
=
cell_corners_xyz
[
2
];
double
surface_normal_of_the_cell
[
3
]
;
Point3D
surface_normal_of_the_cell
;
find_unit_normal
(
pcorner_coordinates1
,
pcorner_coordinates2
,
pcorner_coordinates3
,
surface_normal_of_the_cell
);
// The surface normal is used to choose the coordinate to ignore.
const
auto
abs_x
=
std
::
fabs
(
surface_normal_of_the_cell
[
0
]
);
const
auto
abs_y
=
std
::
fabs
(
surface_normal_of_the_cell
[
1
]
);
const
auto
abs_z
=
std
::
fabs
(
surface_normal_of_the_cell
[
2
]
);
const
auto
abs_x
=
std
::
fabs
(
surface_normal_of_the_cell
.
X
);
const
auto
abs_y
=
std
::
fabs
(
surface_normal_of_the_cell
.
Y
);
const
auto
abs_z
=
std
::
fabs
(
surface_normal_of_the_cell
.
Z
);
int
coordinate_to_ignore
=
3
;
...
...
@@ -208,7 +208,7 @@ winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners,
const
auto
&
point2
=
cell_corners
[
k
+
1
];
if
(
is_point_left_of_edge
(
point1
,
point2
,
point
)
>
0
)
winding_number
++
;
//
printf(" 1: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
//printf(" 1: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
}
}
else
...
...
@@ -219,7 +219,7 @@ winding_numbers_algorithm(const Varray<Point> &cell_corners, int number_corners,
const
auto
&
point2
=
cell_corners
[
k
+
1
];
if
(
is_point_left_of_edge
(
point1
,
point2
,
point
)
<
0
)
winding_number
--
;
//
printf(" 2: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
//printf(" 2: %d %d %g %g %g %g %g %g\n", k, winding_number, point1.x, point1.y, point2.x, point2.y, point.x, point.y);
}
}
}
...
...
@@ -378,7 +378,7 @@ get_no_unique_center_points(size_t gridsize, size_t nx, const Varray<double> &ce
return
no_unique_center_points
;
}
void
set_cell_corners_xyz
(
int
ncorner
,
const
double
*
cell_corners_lon
,
const
double
*
cell_corners_lat
,
double
*
cell_corners_xyz
)
void
set_cell_corners_xyz
(
int
ncorner
,
const
double
*
cell_corners_lon
,
const
double
*
cell_corners_lat
,
Varray
<
Point3D
>
&
cell_corners_xyz
)
{
double
corner_coordinates
[
3
];
for
(
int
corner_no
=
0
;
corner_no
<
ncorner
;
corner_no
++
)
...
...
@@ -387,29 +387,25 @@ void set_cell_corners_xyz(int ncorner, const double *cell_corners_lon, const dou
gcLLtoXYZ
(
cell_corners_lon
[
corner_no
],
cell_corners_lat
[
corner_no
],
corner_coordinates
);
// The components of the result vector are appended to the list of cell corner coordinates.
const
auto
off
=
corner_no
*
3
;
cell_corners_xyz
[
off
+
0
]
=
corner_coordinates
[
0
];
cell_corners_xyz
[
off
+
1
]
=
corner_coordinates
[
1
];
cell_corners_xyz
[
off
+
2
]
=
corner_coordinates
[
2
];
cell_corners_xyz
[
corner_no
].
X
=
corner_coordinates
[
0
];
cell_corners_xyz
[
corner_no
].
Y
=
corner_coordinates
[
1
];
cell_corners_xyz
[
corner_no
].
Z
=
corner_coordinates
[
2
];
}
cell_corners_xyz
[
ncorner
*
3
+
0
]
=
cell_corners_xyz
[
0
];
cell_corners_xyz
[
ncorner
*
3
+
1
]
=
cell_corners_xyz
[
1
];
cell_corners_xyz
[
ncorner
*
3
+
2
]
=
cell_corners_xyz
[
2
];
}
void
set_center_point_plane_projection
(
int
coordinate_to_ignore
,
const
double
(
&
center_point_xyz
)[
3
]
,
Point
&
center_point_plane_projection
)
void
set_center_point_plane_projection
(
int
coordinate_to_ignore
,
const
Point3D
&
center_point_xyz
,
Point
&
center_point_plane_projection
)
{
// clang-format off
switch
(
coordinate_to_ignore
)
{
case
1
:
center_point_plane_projection
=
Point
{
center_point_xyz
[
1
]
,
center_point_xyz
[
2
]
};
break
;
case
2
:
center_point_plane_projection
=
Point
{
center_point_xyz
[
2
]
,
center_point_xyz
[
0
]
};
break
;
case
3
:
center_point_plane_projection
=
Point
{
center_point_xyz
[
0
]
,
center_point_xyz
[
1
]
};
break
;
case
1
:
center_point_plane_projection
=
Point
{
center_point_xyz
.
Y
,
center_point_xyz
.
Z
};
break
;
case
2
:
center_point_plane_projection
=
Point
{
center_point_xyz
.
Z
,
center_point_xyz
.
X
};
break
;
case
3
:
center_point_plane_projection
=
Point
{
center_point_xyz
.
X
,
center_point_xyz
.
Y
};
break
;
}
// clang-format on
}
void
set_cell_corners_plane_projection
(
int
coordinate_to_ignore
,
int
ncorner
,
const
double
*
cell_corners_xyz
,
void
set_cell_corners_plane_projection
(
int
coordinate_to_ignore
,
int
ncorner
,
const
Varray
<
Point3D
>
&
cell_corners_xyz
,
Varray
<
Point
>
&
cell_corners_plane_projection
)
{
switch
(
coordinate_to_ignore
)
...
...
@@ -417,27 +413,75 @@ void set_cell_corners_plane_projection(int coordinate_to_ignore, int ncorner, co
case
1
:
for
(
int
corner_no
=
0
;
corner_no
<=
ncorner
;
corner_no
++
)
{
cell_corners_plane_projection
[
corner_no
].
x
=
cell_corners_xyz
[
corner_no
*
3
+
1
];
cell_corners_plane_projection
[
corner_no
].
y
=
cell_corners_xyz
[
corner_no
*
3
+
2
];
cell_corners_plane_projection
[
corner_no
]
=
Point
{
cell_corners_xyz
[
corner_no
].
Y
,
cell_corners_xyz
[
corner_no
].
Z
};
}
break
;
case
2
:
for
(
int
corner_no
=
0
;
corner_no
<=
ncorner
;
corner_no
++
)
{
cell_corners_plane_projection
[
corner_no
].
x
=
cell_corners_xyz
[
corner_no
*
3
+
2
];
cell_corners_plane_projection
[
corner_no
].
y
=
cell_corners_xyz
[
corner_no
*
3
+
0
];
cell_corners_plane_projection
[
corner_no
]
=
Point
{
cell_corners_xyz
[
corner_no
].
Z
,
cell_corners_xyz
[
corner_no
].
X
};
}
break
;
case
3
:
for
(
int
corner_no
=
0
;
corner_no
<=
ncorner
;
corner_no
++
)
{
cell_corners_plane_projection
[
corner_no
].
x
=
cell_corners_xyz
[
corner_no
*
3
+
0
];
cell_corners_plane_projection
[
corner_no
].
y
=
cell_corners_xyz
[
corner_no
*
3
+
1
];
cell_corners_plane_projection
[
corner_no
]
=
Point
{
cell_corners_xyz
[
corner_no
].
X
,
cell_corners_xyz
[
corner_no
].
Y
};
}
break
;
}
}
int
get_actual_number_of_corners
(
int
ncorner
,
const
Varray
<
Point3D
>
&
cell_corners_xyz_open_cell
)
{
auto
actual_number_of_corners
=
ncorner
;
for
(
int
corner_no
=
ncorner
-
1
;
corner_no
>
0
;
corner_no
--
)
{
if
(
IS_EQUAL
(
cell_corners_xyz_open_cell
[
corner_no
].
X
,
cell_corners_xyz_open_cell
[
corner_no
-
1
].
X
)
&&
IS_EQUAL
(
cell_corners_xyz_open_cell
[
corner_no
].
Y
,
cell_corners_xyz_open_cell
[
corner_no
-
1
].
Y
)
&&
IS_EQUAL
(
cell_corners_xyz_open_cell
[
corner_no
].
Z
,
cell_corners_xyz_open_cell
[
corner_no
-
1
].
Z
))
actual_number_of_corners
--
;
else
break
;
}
return
actual_number_of_corners
;
}
int
get_no_duplicates
(
int
actual_number_of_corners
,
const
Varray
<
Point3D
>
&
cell_corners_xyz_open_cell
,
std
::
vector
<
bool
>
&
marked_duplicate_indices
)
{
for
(
int
i
=
0
;
i
<
actual_number_of_corners
;
i
++
)
marked_duplicate_indices
[
i
]
=
false
;
int
no_duplicates
=
0
;
for
(
int
i
=
0
;
i
<
actual_number_of_corners
;
i
++
)
for
(
int
j
=
i
+
1
;
j
<
actual_number_of_corners
;
j
++
)
if
(
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
].
X
-
cell_corners_xyz_open_cell
[
j
].
X
)
<
eps
&&
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
].
Y
-
cell_corners_xyz_open_cell
[
j
].
Y
)
<
eps
&&
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
].
Z
-
cell_corners_xyz_open_cell
[
j
].
Z
)
<
eps
)
{
no_duplicates
+=
1
;
marked_duplicate_indices
[
j
]
=
true
;
}
return
no_duplicates
;
}
void
copy_unique_corners
(
int
actual_number_of_corners
,
const
Varray
<
Point3D
>
&
cell_corners_xyz_open_cell
,
std
::
vector
<
bool
>
&
marked_duplicate_indices
,
Varray
<
Point3D
>
&
cell_corners_xyz_without_duplicates
)
{
int
unique_corner_number
=
0
;
for
(
int
corner_no
=
0
;
corner_no
<
actual_number_of_corners
;
corner_no
++
)
{
if
(
marked_duplicate_indices
[
corner_no
]
==
false
)
{
cell_corners_xyz_without_duplicates
[
unique_corner_number
]
=
cell_corners_xyz_open_cell
[
corner_no
];
unique_corner_number
++
;
}
}
}
static
void
verify_grid
(
int
gridtype
,
size_t
gridsize
,
size_t
nx
,
int
gridno
,
int
ngrids
,
int
ncorner
,
double
*
grid_center_lon
,
double
*
grid_center_lat
,
double
*
grid_corner_lon
,
double
*
grid_corner_lat
)
...
...
@@ -458,10 +502,9 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
The results of the tests are printed on stdout.
*/
double
center_point_xyz
[
3
]
;
Varray
<
double
>
cell_corners_xyz_open_cell
(
3
*
ncorner
);
Point3D
center_point_xyz
;
Varray
<
Point3D
>
cell_corners_xyz_open_cell
(
ncorner
);
double
corner_coordinates
[
3
];
Point
center_point_plane_projection
;
size_t
no_of_cells_with_duplicates
=
0
;
...
...
@@ -486,8 +529,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
// used only actual_number_of_corners
std
::
vector
<
bool
>
marked_duplicate_indices
(
ncorner
);
Varray
<
double
>
cell_corners_xyz_without_duplicates
(
3
*
ncorner
);
Varray
<
double
>
cell_corners_xyz
(
3
*
(
ncorner
+
1
));
Varray
<
Point3D
>
cell_corners_xyz
(
ncorner
+
1
);
Varray
<
Point
>
cell_corners_plane_projection
(
ncorner
+
1
);
Varray
<
size_t
>
cells_with_duplicates
;
cells_with_duplicates
.
reserve
(
gridsize
);
...
...
@@ -501,8 +543,13 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
for
(
size_t
cell_no
=
0
;
cell_no
<
gridsize
;
cell_no
++
)
{
// Conversion of center point spherical coordinates to Cartesian coordinates.
gcLLtoXYZdeg
(
grid_center_lon
[
cell_no
],
grid_center_lat
[
cell_no
],
center_point_xyz
);
double
center_coordinates
[
3
];
gcLLtoXYZdeg
(
grid_center_lon
[
cell_no
],
grid_center_lat
[
cell_no
],
center_coordinates
);
center_point_xyz
.
X
=
center_coordinates
[
0
];
center_point_xyz
.
Y
=
center_coordinates
[
1
];
center_point_xyz
.
Z
=
center_coordinates
[
2
];
double
corner_coordinates
[
3
];
for
(
int
corner_no
=
0
;
corner_no
<
ncorner
;
corner_no
++
)
{
// Conversion of corner spherical coordinates to Cartesian coordinates.
...
...
@@ -510,10 +557,9 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
corner_coordinates
);
// The components of the result vector are appended to the list of cell corner coordinates.
const
int
off
=
corner_no
*
3
;
cell_corners_xyz_open_cell
[
off
+
0
]
=
corner_coordinates
[
0
];
cell_corners_xyz_open_cell
[
off
+
1
]
=
corner_coordinates
[
1
];
cell_corners_xyz_open_cell
[
off
+
2
]
=
corner_coordinates
[
2
];
cell_corners_xyz_open_cell
[
corner_no
].
X
=
corner_coordinates
[
0
];
cell_corners_xyz_open_cell
[
corner_no
].
Y
=
corner_coordinates
[
1
];
cell_corners_xyz_open_cell
[
corner_no
].
Z
=
corner_coordinates
[
2
];
}
/*
...
...
@@ -522,24 +568,11 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
following identifies the surplus corners and gives the correct length of the cell.
*/
auto
actual_number_of_corners
=
ncorner
;
for
(
int
corner_no
=
ncorner
-
1
;
corner_no
>
0
;
corner_no
--
)
{
const
auto
off
=
corner_no
*
3
;
const
auto
off2
=
(
corner_no
-
1
)
*
3
;
if
(
IS_EQUAL
(
cell_corners_xyz_open_cell
[
off
+
0
],
cell_corners_xyz_open_cell
[
off2
+
0
])
&&
IS_EQUAL
(
cell_corners_xyz_open_cell
[
off
+
1
],
cell_corners_xyz_open_cell
[
off2
+
1
])
&&
IS_EQUAL
(
cell_corners_xyz_open_cell
[
off
+
2
],
cell_corners_xyz_open_cell
[
off2
+
2
]))
actual_number_of_corners
--
;
else
break
;
}
auto
actual_number_of_corners
=
get_actual_number_of_corners
(
ncorner
,
cell_corners_xyz_open_cell
);
no_cells_with_a_specific_no_of_corners
[
actual_number_of_corners
-
1
]
+=
1
;
/* If there are less than three corners in the cell, it is unusable and considered degenerate. No area can be
* computed. */
// If there are less than three corners in the cell, it is unusable and considered degenerate. No area can be computed.
if
(
actual_number_of_corners
<
3
)
{
...
...
@@ -555,48 +588,33 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
no_usable_cells
++
;
/* Checks if there are any duplicate vertices in the list of corners. Note that the last (additional) corner has
* not been set yet. */
// Checks if there are any duplicate vertices in the list of corners. Note that the last (additional) corner has not been set yet.
for
(
int
i
=
0
;
i
<
actual_number_of_corners
;
i
++
)
marked_duplicate_indices
[
i
]
=
false
;
auto
no_duplicates
=
get_no_duplicates
(
actual_number_of_corners
,
cell_corners_xyz_open_cell
,
marked_duplicate_indices
)
;
int
no_duplicates
=
0
;
for
(
int
i
=
0
;
i
<
actual_number_of_corners
*
3
;
i
=
i
+
3
)
for
(
int
j
=
i
+
3
;
j
<
actual_number_of_corners
*
3
;
j
=
j
+
3
)
if
(
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
+
0
]
-
cell_corners_xyz_open_cell
[
j
+
0
])
<
eps
&&
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
+
1
]
-
cell_corners_xyz_open_cell
[
j
+
1
])
<
eps
&&
std
::
fabs
(
cell_corners_xyz_open_cell
[
i
+
2
]
-
cell_corners_xyz_open_cell
[
j
+
2
])
<
eps
)
if
(
Options
::
cdoVerbose
)
{
for
(
int
i
=
0
;
i
<
actual_number_of_corners
;
i
++
)
{
if
(
Options
::
cdoVerbose
)
if
(
marked_duplicate_indices
[
i
]
)
{
fprintf
(
stdout
,
"Duplicate vertex [lon=%.5g lat=%.5g] was found in cell no %zu"
,
grid_corner_lon
[
cell_no
*
ncorner
+
i
/
3
],
grid_corner_lat
[
cell_no
*
ncorner
+
i
/
3
],
cell_no
+
1
);
grid_corner_lon
[
cell_no
*
ncorner
+
i
],
grid_corner_lat
[
cell_no
*
ncorner
+
i
],
cell_no
+
1
);
if
(
nx
)
printIndex2D
(
cell_no
,
nx
);
fprintf
(
stdout
,
"
\n
"
);
}
no_duplicates
+=
1
;
marked_duplicate_indices
[
j
/
3
]
=
true
;
}
}
// Writes the unique corner vertices in a new array.
int
unique_corner
_number
=
0
;
copy_
unique_corner
s
(
actual_number_of_corners
,
cell_corners_xyz_open_cell
,
marked_duplicate_indices
,
cell_corners_xyz
)
;
for
(
int
corner_no
=
0
;
corner_no
<
actual_number_of_corners
;
corner_no
++
)
{
if
(
marked_duplicate_indices
[
corner_no
]
==
false
)
{
const
auto
off
=
corner_no
*
3
;
const
auto
off2
=
unique_corner_number
*
3
;
cell_corners_xyz_without_duplicates
[
off2
+
0
]
=
cell_corners_xyz_open_cell
[
off
+
0
];
cell_corners_xyz_without_duplicates
[
off2
+
1
]
=
cell_corners_xyz_open_cell
[
off
+
1
];
cell_corners_xyz_without_duplicates
[
off2
+
2
]
=
cell_corners_xyz_open_cell
[
off
+
2
];
unique_corner_number
+=
1
;
}
}
actual_number_of_corners
-=
no_duplicates
;
actual_number_of_corners
=
actual_number_of_corners
-
no_duplicates
;
// We are creating a closed polygon/cell by setting the additional last corner to be the same as the first one.
cell_corners_xyz
[
actual_number_of_corners
]
=
cell_corners_xyz
[
0
];
if
(
no_duplicates
!=
0
)
{
...
...
@@ -618,21 +636,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
continue
;
}
// We are creating a closed polygon/cell by setting the additional last corner to be the same as the first one.
for
(
int
corner_no
=
0
;
corner_no
<
actual_number_of_corners
;
corner_no
++
)
{
const
auto
off
=
corner_no
*
3
;
cell_corners_xyz
[
off
+
0
]
=
cell_corners_xyz_without_duplicates
[
off
+
0
];
cell_corners_xyz
[
off
+
1
]
=
cell_corners_xyz_without_duplicates
[
off
+
1
];
cell_corners_xyz
[
off
+
2
]
=
cell_corners_xyz_without_duplicates
[
off
+
2
];
}
cell_corners_xyz
[
actual_number_of_corners
*
3
+
0
]
=
cell_corners_xyz
[
0
];
cell_corners_xyz
[
actual_number_of_corners
*
3
+
1
]
=
cell_corners_xyz
[
1
];
cell_corners_xyz
[
actual_number_of_corners
*
3
+
2
]
=
cell_corners_xyz
[
2
];
const
auto
coordinate_to_ignore
=
find_coordinate_to_ignore
(
cell_corners_xyz
.
data
());
const
auto
coordinate_to_ignore
=
find_coordinate_to_ignore
(
cell_corners_xyz
);
/* The remaining two-dimensional coordinates are extracted into one array for all the cell's corners and into one
* array for the center point. */
...
...
@@ -643,10 +647,11 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
vice versa. */
bool
invert_result
=
false
;
if
(
cell_corners_xyz
[
coordinate_to_ignore
-
1
]
<
0
)
invert_result
=
true
;
const
auto
cval
=
(
coordinate_to_ignore
==
1
)
?
cell_corners_xyz
[
0
].
X
:
((
coordinate_to_ignore
==
2
)
?
cell_corners_xyz
[
0
].
Y
:
cell_corners_xyz
[
0
].
Z
);
if
(
cval
<
0.0
)
invert_result
=
true
;
set_center_point_plane_projection
(
coordinate_to_ignore
,
center_point_xyz
,
center_point_plane_projection
);
set_cell_corners_plane_projection
(
coordinate_to_ignore
,
actual_number_of_corners
,
cell_corners_xyz
.
data
()
,
cell_corners_plane_projection
);
set_cell_corners_plane_projection
(
coordinate_to_ignore
,
actual_number_of_corners
,
cell_corners_xyz
,
cell_corners_plane_projection
);
// Checking for convexity of the cell.
...
...
@@ -702,8 +707,7 @@ verify_grid(int gridtype, size_t gridsize, size_t nx, int gridno, int ngrids, in
center_point_plane_projection.x, center_point_plane_projection.y);
*/
// The winding numbers algorithm is used to test whether the presumed center point is within the bounds of the cell.
auto
winding_number
=
winding_numbers_algorithm
(
cell_corners_plane_projection
,
actual_number_of_corners
+
1
,
center_point_plane_projection
);
auto
winding_number
=
winding_numbers_algorithm
(
cell_corners_plane_projection
,
actual_number_of_corners
+
1
,
center_point_plane_projection
);
// if ( winding_number == 0 ) printf("%d,", cell_no+1);
if
(
winding_number
==
0
)
no_of_cells_with_center_points_out_of_bounds
+=
1
;
...
...
src/verifygrid.h
View file @
257074b2
...
...
@@ -17,12 +17,13 @@
#ifndef VERIFYGRID_H
#define VERIFYGRID_H
struct
Point
{
double
x
,
y
;
};
struct
Point
{
double
x
,
y
;
};
struct
Point3D
{
double
X
,
Y
,
Z
;
};
void
set_cell_corners_xyz
(
int
ncorner
,
const
double
*
cell_corners_lon
,
const
double
*
cell_corners_lat
,
double
*
cell_corners_xyz
);
void
set_center_point_plane_projection
(
int
coordinate_to_ignore
,
const
double
(
&
center_point_xyz
)[
3
]
,
Point
&
center_point_plane_projection
);
void
set_cell_corners_plane_projection
(
int
coordinate_to_ignore
,
int
ncorner
,
const
double
*
cell_corners_xyz
,
Varray
<
Point
>
&
cell_corners_plane_projection
);
int
find_coordinate_to_ignore
(
const
double
*
cell_corners_xyz
);
void
set_cell_corners_xyz
(
int
ncorner
,
const
double
*
cell_corners_lon
,
const
double
*
cell_corners_lat
,
Varray
<
Point3D
>
&
cell_corners_xyz
);
void
set_center_point_plane_projection
(
int
coordinate_to_ignore
,
const
Point3D
&
center_point_xyz
,
Point
&
center_point_plane_projection
);
void
set_cell_corners_plane_projection
(
int
coordinate_to_ignore
,
int
ncorner
,
const
Varray
<
Point3D
>
&
cell_corners_xyz
,
Varray
<
Point
>
&
cell_corners_plane_projection
);
int
find_coordinate_to_ignore
(
const
Varray
<
Point3D
>
&
cell_corners_xyz
);
double
calculate_the_polygon_area
(
const
Varray
<
Point
>
&
cell_corners
,
int
number_corners
);
bool
are_polygon_vertices_arranged_in_clockwise_order
(
double
cell_area
);
int
winding_numbers_algorithm
(
const
Varray
<
Point
>
&
cell_corners
,
int
number_corners
,
const
Point
&
point
);
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment