Update point_is_in_domain function

doc/modules/changes/20231214_glerum

@@ -0,0 +1,5 @@
+ Improved: The geometry function point_is_in_domain
+ now also works for meshes that include initial topography
+ and/or mesh deformation for box and chunk geometries.
+ <br>
+ (Anne Glerum, 2023/12/14)

include/aspect/utilities.h

@@ -56,6 +56,7 @@ namespace aspect
     using namespace dealii;
     using namespace dealii::Utilities;
 
+
     /**
      * Given an array @p values, consider three cases:
      * - If it has size @p N, return the original array.
@@ -356,6 +357,19 @@ namespace aspect
 
 
 
+    /**
+     * Given a point @p point, find out if any of the MPI
+     * processes own the cell in which this point lies. If
+     * not, the point lies outside the @p triangulation.
+     */
+    template <int dim>
+    bool
+    point_is_in_triangulation(const Mapping<dim> &mapping,
+                              const parallel::distributed::Triangulation<dim> &triangulation,
+                              const Point<dim> &point,
+                              const MPI_Comm mpi_communicator);
+
+
     namespace Coordinates
     {
 

source/geometry_model/box.cc

@@ -21,7 +21,9 @@
 
 #include <aspect/geometry_model/box.h>
 #include <aspect/geometry_model/initial_topography_model/zero_topography.h>
+#include <aspect/mesh_deformation/interface.h>
 #include <aspect/simulator_signals.h>
+#include <aspect/utilities.h>
 
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria_iterator.h>
@@ -316,20 +318,53 @@ namespace aspect
     bool
     Box<dim>::point_is_in_domain(const Point<dim> &point) const
     {
-      AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
-                  this->simulator_is_past_initialization() == false,
-                  ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not."));
+      // If mesh deformation is enabled, we have to loop over all the current
+      // grid cells to see if the given point lies in the domain.
+      // If mesh deformation is not enabled, or has not happened yet,
+      // we can use the global extents of the model domain with or without
+      // initial topography instead.
+      // This function only checks if the given point lies in the domain
+      // in its current shape at the current time. It can be called before
+      // mesh deformation is applied in the first timestep (e.g., by the boundary
+      // traction plugins), and therefore there is no guarantee
+      // that the point will still lie in the domain after initial mesh deformation.
+      if (this->get_parameters().mesh_deformation_enabled &&
+          this->simulator_is_past_initialization())
+        {
+          return Utilities::point_is_in_triangulation<dim>(this->get_mapping(),
+                                                           this->get_triangulation(),
+                                                           point,
+                                                           this->get_mpi_communicator());
+        }
+      // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain.
+      else
+        {
+          // The maximal extents of the unperturbed box domain.
+          Point<dim> max_point = extents+box_origin;
 
-      AssertThrow(Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()),
-                  ExcMessage("After adding topography, this function can no longer be used "
-                             "to determine whether a point lies in the domain or not."));
+          // If mesh deformation is not enabled, but initial topography
+          // was/will be applied to the mesh, include this topography in the
+          // extent of the domain.
+          if (!Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()))
+            {
+              // Get the surface x (,y) point
+              Point<dim-1> surface_point;
+              for (unsigned int d=0; d<dim-1; ++d)
+                surface_point[d] = point[d];
+
+              // Get the surface topography at this point
+              const double topo = topo_model->value(surface_point);
+              max_point[dim-1] += topo;
+            }
 
-      for (unsigned int d = 0; d < dim; ++d)
-        if (point[d] > extents[d]+box_origin[d]+std::numeric_limits<double>::epsilon()*extents[d] ||
-            point[d] < box_origin[d]-std::numeric_limits<double>::epsilon()*extents[d])
-          return false;
+          // Check whether point lies within the min/max coordinates of the domain including initial topography.
+          for (unsigned int d = 0; d < dim; ++d)
+            if (point[d] > max_point[d]+std::numeric_limits<double>::epsilon()*extents[d] ||
+                point[d] < box_origin[d]-std::numeric_limits<double>::epsilon()*extents[d])
+              return false;
 
-      return true;
+          return true;
+        }
     }
 
     template <int dim>

source/geometry_model/chunk.cc

@@ -697,21 +697,33 @@ namespace aspect
     bool
     Chunk<dim>::point_is_in_domain(const Point<dim> &point) const
     {
-      AssertThrow(!this->get_parameters().mesh_deformation_enabled ||
-                  this->simulator_is_past_initialization() == false,
-                  ExcMessage("After displacement of the free surface, this function can no longer be used to determine whether a point lies in the domain or not."));
-
-      AssertThrow(Plugins::plugin_type_matches<const InitialTopographyModel::ZeroTopography<dim>>(this->get_initial_topography_model()),
-                  ExcMessage("After adding topography, this function can no longer be used to determine whether a point lies in the domain or not."));
-
-      const Point<dim> spherical_point = manifold->pull_back(point);
+      // If mesh deformation is enabled, we have to loop over all the current
+      // grid cells to see if the given point lies in the domain.
+      // If mesh deformation is not enabled, or has not happened yet,
+      // we can use the global extents of the model domain with or without
+      // initial topography instead.
+      // This function only checks if the given point lies in the domain
+      // in its current shape at the current time. It can be called before
+      // mesh deformation is applied in the first timestep (e.g., by the boundary
+      // traction plugins), and therefore there is no guarantee
+      // that the point will still lie in the domain after initial mesh deformation.
+      if (this->get_parameters().mesh_deformation_enabled &&
+          this->simulator_is_past_initialization())
+        {
+          return Utilities::point_is_in_triangulation(this->get_mapping(), this->get_triangulation(), point, this->get_mpi_communicator());
+        }
+      // Without mesh deformation enabled, it is much cheaper to check whether the point lies in the domain.
+      else
+        {
+          const Point<dim> spherical_point = manifold->pull_back(point);
 
-      for (unsigned int d = 0; d < dim; ++d)
-        if (spherical_point[d] > point2[d]+std::numeric_limits<double>::epsilon()*std::abs(point2[d]) ||
-            spherical_point[d] < point1[d]-std::numeric_limits<double>::epsilon()*std::abs(point2[d]))
-          return false;
+          for (unsigned int d = 0; d < dim; ++d)
+            if (spherical_point[d] > point2[d]+std::numeric_limits<double>::epsilon()*std::abs(point2[d]) ||
+                spherical_point[d] < point1[d]-std::numeric_limits<double>::epsilon()*std::abs(point2[d]))
+              return false;
 
-      return true;
+          return true;
+        }
     }
 
 

source/utilities.cc

@@ -39,8 +39,7 @@
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/signaling_nan.h>
 #include <deal.II/base/patterns.h>
-
-
+#include <deal.II/grid/grid_tools.h>
 
 #include <cerrno>
 #include <dirent.h>
@@ -214,6 +213,9 @@ namespace aspect
 
 
 
+
+
+
     template <typename T>
     Table<2,T>
     parse_input_table (const std::string &input_string,
@@ -729,6 +731,34 @@ namespace aspect
 
 
 
+    template <int dim>
+    bool
+    point_is_in_triangulation(const Mapping<dim> &mapping,
+                              const parallel::distributed::Triangulation<dim> &triangulation,
+                              const Point<dim> &point,
+                              const MPI_Comm mpi_communicator)
+    {
+      // Try to find the cell around the given point.
+      bool cell_found = false;
+      std::pair<const typename parallel::distributed::Triangulation<dim>::active_cell_iterator,
+          Point<dim>> it =
+            GridTools::find_active_cell_around_point<>(mapping, triangulation, point);
+
+      // If we found the correct cell on this MPI process, we have found the right cell.
+      if (it.first.state() == IteratorState::valid && it.first->is_locally_owned())
+        cell_found = true;
+
+      // Compute how many processes found the cell.
+      const int n_procs_cell_found = Utilities::MPI::sum(cell_found ? 1 : 0, mpi_communicator);
+      // If at least one process found the cell, the point is in the triangulation.
+      if (n_procs_cell_found > 0)
+        return true;
+      else
+        return false;
+    }
+
+
+
     namespace Coordinates
     {
 
@@ -3428,6 +3458,12 @@ namespace aspect
                                    const dealii::Point<2> &point); \
   \
   template \
+  bool point_is_in_triangulation<dim>(const Mapping<dim> &mapping, \
+                                      const parallel::distributed::Triangulation<dim> &triangulation, \
+                                      const Point<dim> &point, \
+                                      const MPI_Comm mpi_communicator); \
+  \
+  template \
   double signed_distance_to_polygon<dim>(const std::vector<Point<2>> &pointList, \
                                          const dealii::Point<2> &point); \
   \

tests/airy_isostasy_initial_topo.cc

@@ -0,0 +1,77 @@
+/*
+  Copyright (C) 2022 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#include <aspect/postprocess/interface.h>
+#include <aspect/geometry_model/interface.h>
+#include <aspect/simulator_access.h>
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    template <int dim>
+    class PointIsInDomain : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+    {
+      public:
+        virtual
+        std::pair<std::string,std::string>
+        execute (TableHandler &);
+    };
+  }
+}
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    template <int dim>
+    std::pair<std::string,std::string>
+    PointIsInDomain<dim>::execute(TableHandler &)
+    {
+      // The airy isostasy box domain is 1x1 with max 0.05 topography.
+      const Point<dim> point_in_domain(0.1, 0.1);
+      const Point<dim> point_not_in_domain(0.1, 2.0);
+      const Point<dim> point_not_in_domain_2(2.0, 0.01);
+
+      std::ostringstream screen_text;
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_in_domain) << " ";
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain) << " ";
+      screen_text << std::boolalpha << this->get_geometry_model().point_is_in_domain(point_not_in_domain_2);
+
+      return std::pair<std::string, std::string> ("Points lie in domain: ",
+                                                  screen_text.str());
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Postprocess
+  {
+    ASPECT_REGISTER_POSTPROCESSOR(PointIsInDomain,
+                                  "point is in domain",
+                                  "A postprocessor that tests whether certain points "
+                                  "lie within the domain.")
+  }
+}

tests/airy_isostasy_initial_topo.prm

@@ -30,3 +30,8 @@ subsection Boundary traction model
     set Function expression = 0;10
   end
 end
+
+# Also test the point_is_in_domain function of the geometry model.
+subsection Postprocess
+  set List of postprocessors = velocity statistics, pressure statistics, topography, point is in domain
+end

tests/airy_isostasy_initial_topo/screen-output

@@ -1,4 +1,6 @@
 
+Loading shared library <./libairy_isostasy_initial_topo.debug.so>
+
 Number of active cells: 1,024 (on 6 levels)
 Number of degrees of freedom: 17,989 (8,450+1,089+4,225+4,225)
 
@@ -38,9 +40,10 @@ Number of mesh deformation degrees of freedom: 4,162
    Solving Stokes system... 200+16 iterations.
 
    Postprocessing:
-     RMS, max velocity:    1.76e+06 m/year, 2.3e+06 m/year
-     Pressure min/avg/max: -1.132 Pa, 5.012 Pa, 10 Pa
-     Topography min/max:   0 m, 0.05 m
+     RMS, max velocity:     1.76e+06 m/year, 2.3e+06 m/year
+     Pressure min/avg/max:  -1.132 Pa, 5.012 Pa, 10 Pa
+     Topography min/max:    0 m, 0.05 m
+     Points lie in domain:  true false false
 
 *** Timestep 1:  t=1.69957e-09 years, dt=1.69957e-09 years
    Solving mesh displacement system... 8 iterations.
@@ -50,9 +53,10 @@ Number of mesh deformation degrees of freedom: 4,162
    Solving Stokes system... 200+10 iterations.
 
    Postprocessing:
-     RMS, max velocity:    1.62e+06 m/year, 2.11e+06 m/year
-     Pressure min/avg/max: -1.119 Pa, 5.004 Pa, 10.01 Pa
-     Topography min/max:   -4.848e-05 m, 0.04609 m
+     RMS, max velocity:     1.62e+06 m/year, 2.11e+06 m/year
+     Pressure min/avg/max:  -1.119 Pa, 5.004 Pa, 10.01 Pa
+     Topography min/max:    -4.848e-05 m, 0.04609 m
+     Points lie in domain:  true false false
 
 *** Timestep 2:  t=3.5e-09 years, dt=1.80043e-09 years
    Solving mesh displacement system... 8 iterations.
@@ -62,9 +66,10 @@ Number of mesh deformation degrees of freedom: 4,162
    Solving Stokes system... 200+4 iterations.
 
    Postprocessing:
-     RMS, max velocity:    1.53e+06 m/year, 1.99e+06 m/year
-     Pressure min/avg/max: -0.8445 Pa, 5.003 Pa, 10.02 Pa
-     Topography min/max:   -9.971e-05 m, 0.04229 m
+     RMS, max velocity:     1.53e+06 m/year, 1.99e+06 m/year
+     Pressure min/avg/max:  -0.8445 Pa, 5.003 Pa, 10.02 Pa
+     Topography min/max:    -9.971e-05 m, 0.04229 m
+     Points lie in domain:  true false false
 
 Termination requested by criterion: end time