Skip to content

Commit

Permalink
Merge pull request #5509 from anne-glerum/update_point_in_domain
Browse files Browse the repository at this point in the history
Update point_is_in_domain function
  • Loading branch information
gassmoeller authored Dec 27, 2023
2 parents 800e803 + 4b997e0 commit f6a03dc
Show file tree
Hide file tree
Showing 15 changed files with 454 additions and 35 deletions.
5 changes: 5 additions & 0 deletions doc/modules/changes/20231214_glerum
Original file line number Diff line number Diff line change
@@ -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)
14 changes: 14 additions & 0 deletions include/aspect/utilities.h
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
{

Expand Down
57 changes: 46 additions & 11 deletions source/geometry_model/box.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -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>
Expand Down
38 changes: 25 additions & 13 deletions source/geometry_model/chunk.cc
Original file line number Diff line number Diff line change
Expand Up @@ -685,21 +685,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;
}
}


Expand Down
40 changes: 38 additions & 2 deletions source/utilities.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down Expand Up @@ -214,6 +213,9 @@ namespace aspect






template <typename T>
Table<2,T>
parse_input_table (const std::string &input_string,
Expand Down Expand Up @@ -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
{

Expand Down Expand Up @@ -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); \
\
Expand Down
77 changes: 77 additions & 0 deletions tests/airy_isostasy_initial_topo.cc
Original file line number Diff line number Diff line change
@@ -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.")
}
}
5 changes: 5 additions & 0 deletions tests/airy_isostasy_initial_topo.prm
Original file line number Diff line number Diff line change
Expand Up @@ -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
23 changes: 14 additions & 9 deletions tests/airy_isostasy_initial_topo/screen-output
Original file line number Diff line number Diff line change
@@ -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)

Expand Down Expand Up @@ -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.
Expand All @@ -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.
Expand All @@ -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

Expand Down
Loading

0 comments on commit f6a03dc

Please sign in to comment.