* Compute the volume (i.e. the dim-dimensional measure) of the
* triangulation. We compute the measure using the integral $\sum_K \int_K 1
* \; dx$ where $K$ are the cells of the given triangulation. The integral
- * is approximated via quadrature for which we need the mapping argument.
+ * is approximated via quadrature. This version of the function uses a
+ * linear mapping to compute the JxW values on each cell.
+ *
+ * If the triangulation is a dim-dimensional one embedded in a higher
+ * dimensional space of dimension spacedim, then the value returned is the
+ * dim-dimensional measure. For example, for a two-dimensional triangulation
+ * in three-dimensional space, the value returned is the area of the surface
+ * so described. (This obviously makes sense since the spacedim-dimensional
+ * measure of a dim-dimensional triangulation would always be zero if dim @<
+ * spacedim).
+ *
+ * This function also works for objects of type
+ * parallel::distributed::Triangulation, in which case the function is a
+ * collective operation.
+ *
+ * @param tria The triangulation.
+ * @return The dim-dimensional measure of the domain described by the
+ * triangulation, as discussed above.
+ */
+ template <int dim, int spacedim>
+ double
+ volume(const Triangulation<dim, spacedim> &tria);
+
+ /**
+ * Compute the volume (i.e. the dim-dimensional measure) of the
+ * triangulation. We compute the measure using the integral $\sum_K \int_K 1
+ * \; dx$ where $K$ are the cells of the given triangulation. The integral
+ * is approximated via quadrature for which we use the mapping argument.
*
* If the triangulation is a dim-dimensional one embedded in a higher
* dimensional space of dimension spacedim, then the value returned is the
* collective operation.
*
* @param tria The triangulation.
- * @param mapping An optional argument used to denote the mapping that
- * should be used when describing whether cells are bounded by straight or
- * curved faces. The default is to use a $Q_1$ mapping, which corresponds to
- * straight lines bounding the cells.
+ * @param mapping The Mapping which computes the Jacobians used to
+ * approximate the volume via quadrature. Explicitly using a higher-order
+ * Mapping (i.e., instead of using the other version of this function) will
+ * result in a more accurate approximation of the volume on Triangulations
+ * with curvature described by Manifold objects.
* @return The dim-dimensional measure of the domain described by the
* triangulation, as discussed above.
*/
template <int dim, int spacedim>
double
volume(const Triangulation<dim, spacedim> &tria,
- const Mapping<dim, spacedim> & mapping =
- (ReferenceCells::get_hypercube<dim>()
-#ifndef _MSC_VER
- .template get_default_linear_mapping<dim, spacedim>()
-#else
- .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
-#endif
- ));
+ const Mapping<dim, spacedim> & mapping);
/**
* Return an approximation of the diameter of the smallest active cell of a
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/filtered_iterator.h>
+ template <int dim, int spacedim>
+ double
+ volume(const Triangulation<dim, spacedim> &triangulation)
+ {
+ Assert(triangulation.get_reference_cells().size() == 1,
+ ExcNotImplemented());
+ const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
+ return volume(
+ triangulation,
+ reference_cell.template get_default_linear_mapping<dim, spacedim>());
+ }
+
+
+
template <int dim, int spacedim>
double
volume(const Triangulation<dim, spacedim> &triangulation,
if (const auto *p = dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
mapping_degree = p->get_degree();
else if (const auto *p =
- dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
+ dynamic_cast<const MappingFE<dim, spacedim> *>(&mapping))
mapping_degree = p->get_degree();
// then initialize an appropriate quadrature formula
- const QGauss<dim> quadrature_formula(mapping_degree + 1);
+ Assert(triangulation.get_reference_cells().size() == 1,
+ ExcNotImplemented());
+ const ReferenceCell reference_cell = triangulation.get_reference_cells()[0];
+ const Quadrature<dim> quadrature_formula =
+ reference_cell.template get_gauss_type_quadrature<dim>(mapping_degree +
+ 1);
const unsigned int n_q_points = quadrature_formula.size();
// we really want the JxW values from the FEValues object, but it
// wants a finite element. create a cheap element as a dummy
// element
- FE_Nothing<dim, spacedim> dummy_fe;
+ FE_Nothing<dim, spacedim> dummy_fe(reference_cell);
FEValues<dim, spacedim> fe_values(mapping,
dummy_fe,
quadrature_formula,
update_JxW_values);
- typename Triangulation<dim, spacedim>::active_cell_iterator
- cell = triangulation.begin_active(),
- endc = triangulation.end();
-
double local_volume = 0;
// compute the integral quantities by quadrature
- for (; cell != endc; ++cell)
+ for (const auto &cell : triangulation.active_cell_iterators())
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
diameter(
const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+ template double
+ volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+
template double
volume(const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
const Mapping<deal_II_dimension, deal_II_space_dimension> &);
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include "../tests.h"
deallog << dim << "d, "
<< "hypercube volume, " << i * 2
<< " refinements: " << GridTools::volume(tria) << std::endl;
- };
- };
+ }
+
+ Triangulation<dim> simplex_tria;
+ GridGenerator::convert_hypercube_to_simplex_mesh(tria, simplex_tria);
+ deallog << dim << "d, simplex hypercube volume: "
+ << GridTools::volume(simplex_tria) << std::endl;
+ }
// test 2: hyperball
if (dim >= 2)
Triangulation<dim> tria;
GridGenerator::hyper_ball(tria, Point<dim>(), 1);
- static const SphericalManifold<dim> boundary;
- tria.set_manifold(0, boundary);
-
for (unsigned int i = 0; i < 4; ++i)
{
tria.refine_global(1);
<< "hyperball volume, " << i
<< " refinements: " << GridTools::volume(tria) << std::endl;
}
+
+ Triangulation<dim> simplex_tria;
+ GridGenerator::convert_hypercube_to_simplex_mesh(tria, simplex_tria);
+ simplex_tria.set_all_manifold_ids(numbers::flat_manifold_id);
+ deallog << dim << "d, simplex hyperball volume: "
+ << GridTools::volume(simplex_tria) << std::endl;
+
deallog << "exact value="
<< (dim == 2 ? numbers::PI : 4. / 3. * numbers::PI) << std::endl;
}
DEAL::1d, hypercube volume, 0 refinements: 1.000
DEAL::1d, hypercube volume, 2 refinements: 1.000
+DEAL::1d, simplex hypercube volume: 1.000
DEAL::2d, hypercube volume, 0 refinements: 1.000
DEAL::2d, hypercube volume, 2 refinements: 1.000
+DEAL::2d, simplex hypercube volume: 1.000
DEAL::2d, hyperball volume, 0 refinements: 2.828
DEAL::2d, hyperball volume, 1 refinements: 3.061
DEAL::2d, hyperball volume, 2 refinements: 3.121
DEAL::2d, hyperball volume, 3 refinements: 3.137
+DEAL::2d, simplex hyperball volume: 3.137
DEAL::exact value=3.142
DEAL::3d, hypercube volume, 0 refinements: 1.000
DEAL::3d, hypercube volume, 2 refinements: 1.000
+DEAL::3d, simplex hypercube volume: 1.000
DEAL::3d, hyperball volume, 0 refinements: 3.210
DEAL::3d, hyperball volume, 1 refinements: 3.917
DEAL::3d, hyperball volume, 2 refinements: 4.119
DEAL::3d, hyperball volume, 3 refinements: 4.171
+DEAL::3d, simplex hyperball volume: 4.171
DEAL::exact value=4.189
projects / dealii.git / commitdiff