22 #include "pism/util/petscwrappers/DM.hh"
23 #include "pism/util/petscwrappers/Vec.hh"
24 #include "pism/util/connected_components.hh"
26 #include "pism/frontretreat/util/IcebergRemoverFEM.hh"
28 #include "pism/util/fem/Element.hh"
29 #include "pism/util/array/CellType.hh"
30 #include "pism/util/Mask.hh"
31 #include "pism/util/interpolation.hh"
38 m_mask(grid,
"temporary_mask") {
71 mask_grounded_ice = 1,
72 mask_floating_ice = 2;
88 for (
auto p =
m_grid->points(); p; p.next()) {
89 const int i = p.i(), j = p.j();
124 if (
m_grid->rank() == 0) {
127 true, mask_grounded_ice);
142 auto da =
m_grid->get_dm(1, 0);
143 PetscErrorCode ierr = DMDAGetLocalInfo(*da, &info);
145 throw std::runtime_error(
"Failed to get DMDA info");
153 double mask_iceberg[] = {1.0, 1.0, 1.0, 1.0};
154 double mask_grounded[] = {-1.0, -1.0, -1.0, -1.0};
157 for (
int j = info.gys; j < info.gys + info.gym - 1; j++) {
158 for (
int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
191 list.add(ice_thickness);
193 for (
auto p =
m_grid->points(); p; p.next()) {
194 const int i = p.i(), j = p.j();
197 ice_thickness(i,j) = 0.0;
const std::shared_ptr< const Grid > m_grid
grid used by this component
void failed()
Indicates a failure of a parallel section.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void set_interpolation_type(InterpolationType type)
void get_from_proc0(petsc::Vec &onp0)
Gets a local Array2 from processor 0.
void set(double c)
Result: v[j] <- c for all j.
void put_on_proc0(petsc::Vec &onp0) const
Puts a local array::Scalar on processor 0.
void update_ghosts()
Updates ghost points.
unsigned int stencil_width() const
Get the stencil width of the current Array. Returns 0 if ghosts are not available.
IcebergRemoverFEM(std::shared_ptr< const Grid > g)
void update_impl(const array::Scalar &bc_mask, array::CellType1 &cell_type, array::Scalar &ice_thickness)
std::shared_ptr< petsc::Vec > m_mask_p0
array::Scalar m_iceberg_mask
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
void nodal_values(const array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
Q1 element with sides parallel to X and Y axes.
The 1-point Gaussian quadrature on the square [-1,1]*[-1,1].
Wrapper around VecGetArray and VecRestoreArray.
void label_connected_components(double *image, int nrows, int ncols, bool identify_icebergs, int mask_grounded, int first_label)
bool icy(int M)
Ice-filled cell (grounded or floating).
bool grounded(int M)
Grounded cell (grounded ice or ice-free).