const std::vector<double> & weights,
const std::vector<Tensor<1, spacedim>> &normals);
+ /**
+ * Clears weights, points and normals vectors.
+ */
+ void
+ clear();
+
/**
* Extend the given formula by an additional quadrature point.
* The point, weight and normal should be with respect to reference space,
const hp::QCollection<1> &quadratures1D,
const AdditionalData & additional_data = AdditionalData());
+ /**
+ * Clears the inside, outside and surface quadratures.
+ */
+ void
+ clear_quadratures();
+
/**
* Construct immersed quadratures rules for the incoming level set
* function over the BoundingBox.
void
generate(const Function<dim> &level_set, const BoundingBox<dim> &box);
+ /**
+ * Same as above but does not clear quadratures and appends it to the
+ * existing quadrature instead.
+ */
+ void
+ generate_append(const Function<dim> & level_set,
+ const BoundingBox<dim> &box);
+
/**
* Return the quadrature rule for the region
* $\{x \in B : \psi(x) < 0 \}$
const hp::QCollection<1> &quadratures1D,
const AdditionalData & additional_data = AdditionalData());
+ /**
+ * Clears the inside, outside and surface quadratures.
+ */
+ void
+ clear_quadratures();
+
/**
* Construct immersed quadratures rules for the incoming level set
* function on a given face of the BoundingBox.
const BoundingBox<dim> &box,
const unsigned int face_index);
+ /**
+ * Same as above but does not clear quadratures and appends it to the
+ * existing quadrature instead.
+ */
+ void
+ generate_append(const Function<dim> & level_set,
+ const BoundingBox<dim> &box,
+ const unsigned int face_index);
+
/**
* Return the quadrature rule for the region
* $\{x \in F : \psi(x) < 0 \}$
const hp::QCollection<1> &quadratures1D,
const AdditionalData & additional_data = AdditionalData());
+ /**
+ * Does nothing. Exists for compatibility reasons.
+ */
+ void
+ clear_quadratures();
+
/**
* Construct immersed quadratures rules for the incoming level set
* function on a given face of the BoundingBox.
const BoundingBox<1> &box,
const unsigned int face_index);
+ /**
+ * Same as above but does not clear quadratures and appends it to the
+ * existing quadrature instead.
+ */
+ void
+ generate_append(const Function<1> & level_set,
+ const BoundingBox<1> &box,
+ const unsigned int face_index);
+
/**
* @copydoc FaceQuadratureGenerator<dim>::get_inside_quadrature()
*/
generate(const typename Triangulation<dim>::active_cell_iterator &cell);
private:
+ /**
+ * Construct immersed quadratures for FE_Q_iso_Q1.
+ */
+ void
+ generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box);
+
/**
* Function that describes our level set function in reference space.
*/
const unsigned int face_index);
private:
+ /**
+ * Construct immersed quadratures for FE_Q_iso_Q1.
+ */
+ void
+ generate_fe_q_iso_q1(const BoundingBox<dim> &unit_box,
+ unsigned int face_index);
+
/**
* Function that describes our level set function in reference space.
*/
*/
ExtendableQuadrature(const Quadrature<dim> &quadrature);
+ /**
+ * Clears weights and points vectors.
+ */
+ void
+ clear();
+
/**
* Add a point with an associated weight to the quadrature.
*/
ExtendableQuadrature<dim> &
quadrature_by_definiteness(const Definiteness definiteness);
+ /**
+ * Clears all quadratures.
+ */
+ void
+ clear();
+
/**
* Quadrature for the region $\{x \in B : \psi_i(x) < 0 \forall i \}$ of
* the box, $B$.
virtual void
set_active_cell(
const typename Triangulation<dim>::active_cell_iterator &cell) = 0;
+
+ /**
+ * Set the dof values and the bounding box of the subcell the
+ * function should be evaluated on. Relevant for FE_Q_iso_Q1.
+ */
+ virtual void
+ set_subcell(const std::vector<unsigned int> &mask,
+ const BoundingBox<dim> & subcell_box) = 0;
+
+ /**
+ * Returns flag indicating if the finite element is FE_Q_iso_Q1.
+ */
+ virtual bool
+ is_fe_q_iso_q1() const = 0;
+
+ /**
+ * Number of subdivisions of the FE_Q_iso_Q1 element.
+ */
+ virtual unsigned int
+ n_subdivisions() const = 0;
};
} // namespace DiscreteQuadratureGeneratorImplementation
+ template <int dim, int spacedim>
+ inline void
+ ImmersedSurfaceQuadrature<dim, spacedim>::clear()
+ {
+ this->quadrature_points.clear();
+ this->weights.clear();
+ this->normals.clear();
+ }
+
+
+
template <int dim, int spacedim>
void
ImmersedSurfaceQuadrature<dim, spacedim>::push_back(
#include <deal.II/dofs/dof_accessor.h>
+#include "deal.II/fe/fe_q_iso_q1.h"
#include <deal.II/fe/fe_bernstein.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
face_index,
local_levelset_values);
+ const FiniteElement<dim> &fe =
+ level_set_description->get_fe_collection()[fe_index];
+
+ const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+ dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
+
+ const FE_Poly<dim> *fe_poly = dynamic_cast<const FE_Poly<dim> *>(&fe);
+
+ const bool is_linear = fe_q_iso_q1 != nullptr ||
+ (fe_poly != nullptr && fe_poly->get_degree() == 1);
+
+ // shortcut for linear elements
+ if (is_linear)
+ {
+ return internal::MeshClassifierImplementation::location_from_dof_signs(
+ local_levelset_values);
+ }
+
lagrange_to_bernstein_face[fe_index][face_index].solve(
local_levelset_values);
for (unsigned int i = 0; i < fe_collection.size(); i++)
{
const FiniteElement<dim> &element = fe_collection[i];
- const FE_Q<dim> *fe_q = dynamic_cast<const FE_Q<dim> *>(&element);
+ const FE_Q_Base<dim> * fe_q =
+ dynamic_cast<const FE_Q_Base<dim> *>(&element);
Assert(fe_q != nullptr, ExcNotImplemented());
const FE_Bernstein<dim> fe_bernstein(fe_q->get_degree());
#include <deal.II/base/function_tools.h>
+#include "deal.II/fe/fe_q.h"
+#include "deal.II/fe/fe_q_iso_q1.h"
+
#include <deal.II/grid/reference_cell.h>
#include <deal.II/lac/block_vector.h>
+ template <int dim>
+ inline void
+ ExtendableQuadrature<dim>::clear()
+ {
+ this->quadrature_points.clear();
+ this->weights.clear();
+ }
+
+
+
template <int dim>
void
ExtendableQuadrature<dim>::push_back(const Point<dim> &point,
+ template <int dim>
+ void
+ QPartitioning<dim>::clear()
+ {
+ positive.clear();
+ negative.clear();
+ indefinite.clear();
+ surface.clear();
+ }
+
+
+
/**
* Takes a (dim-1)-dimensional point from the cross-section (orthogonal
* to direction) of the box. Creates the two dim-dimensional points, which
void
QGeneratorBase<dim, spacedim>::clear_quadratures()
{
- q_partitioning = QPartitioning<dim>();
+ q_partitioning.clear();
}
set_active_cell(const typename Triangulation<dim>::active_cell_iterator
&cell) override;
+ /**
+ * @copydoc CellWiseFunction::set_subcell()
+ */
+ void
+ set_subcell(const std::vector<unsigned int> &mask,
+ const BoundingBox<dim> & subcell_box) override;
+
+ /**
+ * @copydoc CellWiseFunction::is_fe_q_iso_q1()
+ */
+ bool
+ is_fe_q_iso_q1() const override;
+
+ /**
+ * @copydoc CellWiseFunction::n_subdivisions()
+ */
+ unsigned int
+ n_subdivisions() const override;
+
/**
* @copydoc Function::value()
*
SmartPointer<const FiniteElement<dim>> element;
/**
- * DOF-indices of the cell in the last call to set_active_cell()
+ * DOF-indices of the cell in the last call to set_active_cell().
*/
std::vector<types::global_dof_index> local_dof_indices;
/**
* Local solution values of the cell in the last call to
- * set_active_cell()
+ * set_active_cell().
*/
std::vector<typename VectorType::value_type> local_dof_values;
+ /**
+ * Local solution values of the subcell after the last call to
+ * set_subcell().
+ */
+ std::vector<typename VectorType::value_type> local_dof_values_subcell;
+
+ /**
+ * Bounding box of the subcell after the last call to set_subcell().
+ */
+ BoundingBox<dim> subcell_box;
+
+ /**
+ * Number of subdivisions per line of the FE_Q_iso_Q1 element. Set to
+ * numbers::invalid_unsigned_int for other elements.
+ */
+ unsigned int n_subdivisions_per_line;
+
/**
* Description of the 1d polynomial basis for tensor product elements
* used for the fast path of this class using tensor product
const VectorType & dof_values)
: dof_handler(&dof_handler)
, global_dof_values(&dof_values)
+ , n_subdivisions_per_line(numbers::invalid_unsigned_int)
{
Assert(dof_handler.n_dofs() == dof_values.size(),
ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size()));
const auto dof_handler_cell =
cell->as_dof_handler_iterator(*dof_handler);
+ const FE_Q<dim> fe_q1(1);
+
// Save the element and the local dof values, since this is what we need
// to evaluate the function.
poly.clear();
element = &dof_handler_cell->get_fe();
- if (element->n_base_elements() == 1 &&
- dealii::internal::FEPointEvaluation::is_fast_path_supported(
- *element, 0))
+ const FiniteElement<dim> *fe = &dof_handler_cell->get_fe();
+
+ if (const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+ dynamic_cast<const FE_Q_iso_Q1<dim> *>(
+ &dof_handler_cell->get_fe()))
+ {
+ this->n_subdivisions_per_line = fe_q_iso_q1->get_degree();
+ fe = &fe_q1;
+ local_dof_values_subcell.resize(fe_q1.n_dofs_per_cell());
+ }
+ else
+ this->n_subdivisions_per_line = numbers::invalid_unsigned_int;
+
+ if (fe->n_base_elements() == 1 &&
+ dealii::internal::FEPointEvaluation::is_fast_path_supported(*fe,
+ 0))
{
dealii::internal::MatrixFreeFunctions::ShapeInfo<double>
shape_info;
renumber = shape_info.lexicographic_numbering;
poly =
dealii::internal::FEPointEvaluation::get_polynomial_space(
- element->base_element(0));
+ fe->base_element(0));
polynomials_are_hat_functions =
(poly.size() == 2 && poly[0].value(0.) == 1. &&
+ template <int dim, typename VectorType>
+ void
+ RefSpaceFEFieldFunction<dim, VectorType>::set_subcell(
+ const std::vector<unsigned int> &mask,
+ const BoundingBox<dim> & subcell_box)
+ {
+ for (unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
+ local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]];
+
+ this->subcell_box = subcell_box;
+ }
+
+
+
+ template <int dim, typename VectorType>
+ bool
+ RefSpaceFEFieldFunction<dim, VectorType>::is_fe_q_iso_q1() const
+ {
+ return n_subdivisions_per_line != numbers::invalid_unsigned_int;
+ }
+
+
+
+ template <int dim, typename VectorType>
+ unsigned int
+ RefSpaceFEFieldFunction<dim, VectorType>::n_subdivisions() const
+ {
+ return n_subdivisions_per_line;
+ }
+
+
+
template <int dim, typename VectorType>
bool
RefSpaceFEFieldFunction<dim, VectorType>::cell_is_set() const
// TODO: this could be extended to a component that is not zero
return dealii::internal::evaluate_tensor_product_value(
poly,
- local_dof_values,
- point,
+ this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+ local_dof_values,
+ this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point,
polynomials_are_hat_functions,
- renumber);
+ this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
}
else
{
// TODO: this could be extended to a component that is not zero
return dealii::internal::evaluate_tensor_product_value_and_gradient(
poly,
- local_dof_values,
- point,
+ this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+ local_dof_values,
+ this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) :
+ point,
polynomials_are_hat_functions,
- renumber)
+ this->is_fe_q_iso_q1() ? std::vector<unsigned int>() :
+ renumber)
.second;
}
else
{
// TODO: this could be extended to a component that is not zero
return dealii::internal::evaluate_tensor_product_hessian(
- poly, local_dof_values, point, renumber);
+ poly,
+ this->is_fe_q_iso_q1() ? local_dof_values_subcell :
+ local_dof_values,
+ this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point,
+ this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
}
else
{
return symmetrize(hessian);
}
}
+
+
+
+ template <int dim>
+ BoundingBox<dim>
+ create_subcell_box(const BoundingBox<dim> & unit_box,
+ const std::array<unsigned int, dim> &subcell_indices,
+ const unsigned int n_subdivisions)
+ {
+ std::pair<Point<dim>, Point<dim>> corners;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ AssertIndexRange(subcell_indices[d], n_subdivisions + 1);
+
+ const double split_box_side_length =
+ (1. / n_subdivisions) * unit_box.side_length(d);
+ corners.first[d] = subcell_indices[d] * split_box_side_length;
+ corners.second[d] = corners.first[d] + split_box_side_length;
+ }
+
+ return BoundingBox<dim>(corners);
+ }
+
+
+
+ template <int dim>
+ std::vector<unsigned int>
+ setup_lexicographic_mask(
+ const std::array<unsigned int, dim> &subcell_indices,
+ unsigned int dofs_per_line)
+ {
+ // This code expands the tensor product space, with `mask` increasing
+ // in size by a factor `2` for each iteration in `d`. Note how the
+ // indices in the next line are appended at the end of the vector.
+ std::vector<unsigned int> mask;
+ mask.reserve(1 << dim);
+ mask.push_back(0);
+ unsigned int stride = 1;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ const unsigned int previous_mask_size = 1U << d;
+ AssertDimension(mask.size(), previous_mask_size);
+ for (unsigned int i = 0; i < previous_mask_size; ++i)
+ {
+ mask[i] += subcell_indices[d] * stride;
+ mask.push_back(stride + mask[i]);
+ }
+ stride *= dofs_per_line;
+ }
+ return mask;
+ }
} // namespace DiscreteQuadratureGeneratorImplementation
} // namespace internal
+ template <int dim>
+ void
+ QuadratureGenerator<dim>::clear_quadratures()
+ {
+ q_generator.clear_quadratures();
+ }
+
+
+
template <int dim>
void
QuadratureGenerator<dim>::generate(const Function<dim> & level_set,
const BoundingBox<dim> &box)
+ {
+ clear_quadratures();
+ generate_append(level_set, box);
+ }
+
+
+
+ template <int dim>
+ void
+ QuadratureGenerator<dim>::generate_append(const Function<dim> & level_set,
+ const BoundingBox<dim> &box)
{
Assert(level_set.n_components == 1,
ExcMessage(
" it should have one component."));
Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume."));
- q_generator.clear_quadratures();
-
std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
level_sets.push_back(level_set);
+ template <int dim>
+ void
+ FaceQuadratureGenerator<dim>::clear_quadratures()
+ {
+ quadrature_generator.clear_quadratures();
+ surface_quadrature.clear();
+ }
+
+
+
template <int dim>
void
FaceQuadratureGenerator<dim>::generate(const Function<dim> & level_set,
const BoundingBox<dim> &box,
const unsigned int face_index)
+ {
+ clear_quadratures();
+ generate_append(level_set, box, face_index);
+ }
+
+
+
+ template <int dim>
+ void
+ FaceQuadratureGenerator<dim>::generate_append(const Function<dim> &level_set,
+ const BoundingBox<dim> &box,
+ const unsigned int face_index)
{
AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
// Reuse the lower dimensional QuadratureGenerator on the face.
const BoundingBox<dim - 1> cross_section =
box.cross_section(face_normal_direction);
- quadrature_generator.generate(face_restriction, cross_section);
+
+ quadrature_generator.generate_append(face_restriction, cross_section);
// We need the dim-dimensional normals of the zero-contour.
// Recompute these.
+ void
+ FaceQuadratureGenerator<1>::clear_quadratures()
+ {
+ // do nothing
+ }
+
+
+
void
FaceQuadratureGenerator<1>::generate(const Function<1> & level_set,
const BoundingBox<1> &box,
const unsigned int face_index)
+ {
+ clear_quadratures();
+ generate_append(level_set, box, face_index);
+ }
+
+
+
+ void
+ FaceQuadratureGenerator<1>::generate_append(const Function<1> & level_set,
+ const BoundingBox<1> &box,
+ const unsigned int face_index)
{
AssertIndexRange(face_index, GeometryInfo<1>::faces_per_cell);
+ template <int dim>
+ void
+ DiscreteQuadratureGenerator<dim>::generate_fe_q_iso_q1(
+ const BoundingBox<dim> &unit_box)
+ {
+ const unsigned int n_subdivisions =
+ reference_space_level_set->n_subdivisions();
+
+ const unsigned int dofs_per_line = n_subdivisions + 1;
+
+ this->clear_quadratures();
+
+ // The triple nested loop below expands the tensor product of the
+ // subdivisions, increasing the size of 'all_subcell_indices' by a factor
+ // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
+ // copies of the data are added at the end of the vector in each
+ // iteration over d. This code avoids a dimension-dependent
+ // number of nested loops over each direction, yet having the same
+ // effect.
+ std::vector<std::array<unsigned int, dim>> all_subcell_indices;
+ all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
+ all_subcell_indices.emplace_back();
+ for (unsigned d = 0; d < dim; ++d)
+ {
+ const unsigned int old_size = all_subcell_indices.size();
+ for (unsigned int i = 1; i < n_subdivisions; ++i)
+ {
+ for (unsigned int j = 0; j < old_size; ++j)
+ {
+ std::array<unsigned int, dim> next_entry =
+ all_subcell_indices[j];
+ next_entry[d] = i;
+ all_subcell_indices.push_back(next_entry);
+ }
+ }
+ }
+
+
+ for (const auto &subcell_indices : all_subcell_indices)
+ {
+ const std::vector<unsigned int> lexicographic_mask =
+ internal::DiscreteQuadratureGeneratorImplementation::
+ setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
+
+ const auto subcell_box =
+ internal::DiscreteQuadratureGeneratorImplementation::
+ create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
+
+ reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
+
+ QuadratureGenerator<dim>::generate_append(*reference_space_level_set,
+ subcell_box);
+ }
+ }
+
+
+
template <int dim>
void
DiscreteQuadratureGenerator<dim>::generate(
reference_space_level_set->set_active_cell(cell);
const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
- QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
+ if (reference_space_level_set->is_fe_q_iso_q1())
+ generate_fe_q_iso_q1(unit_box);
+ else
+ QuadratureGenerator<dim>::generate(*reference_space_level_set, unit_box);
}
+ template <int dim>
+ void
+ DiscreteFaceQuadratureGenerator<dim>::generate_fe_q_iso_q1(
+ const BoundingBox<dim> &unit_box,
+ unsigned int face_index)
+ {
+ const unsigned int n_subdivisions =
+ reference_space_level_set->n_subdivisions();
+
+ const unsigned int dofs_per_line = n_subdivisions + 1;
+
+ this->clear_quadratures();
+
+ const unsigned int coordinate_direction_face = face_index / 2;
+ const unsigned int coordinate_orientation_face = face_index % 2;
+
+ // The triple nested loop below expands the tensor product of the
+ // subdivisions, increasing the size of 'all_subcell_indices' by a factor
+ // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1'
+ // copies of the data are added at the end of the vector in each
+ // iteration over d. This code avoids a dimension-dependent
+ // number of nested loops over each direction, yet having the same
+ // effect. To account for the face all subcell indices in the coordinate of
+ // the face direction are filled with the appropriate index depending on the
+ // coordinate orientation. The other two indices (on the face) are
+ // set accordingly in the nested loops (with a warp around).
+ std::vector<std::array<unsigned int, dim>> all_subcell_indices;
+ all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim));
+ all_subcell_indices.emplace_back();
+ all_subcell_indices[0][coordinate_direction_face] =
+ coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1;
+ for (unsigned d = 0; d < dim - 1; ++d)
+ {
+ const unsigned int old_size = all_subcell_indices.size();
+ for (unsigned int i = 1; i < n_subdivisions; ++i)
+ {
+ for (unsigned int j = 0; j < old_size; ++j)
+ {
+ std::array<unsigned int, dim> next_entry =
+ all_subcell_indices[j];
+ next_entry[(coordinate_direction_face + d + 1) % dim] = i;
+ all_subcell_indices.push_back(next_entry);
+ }
+ }
+ }
+
+ for (const auto &subcell_indices : all_subcell_indices)
+ {
+ const std::vector<unsigned int> lexicographic_mask =
+ internal::DiscreteQuadratureGeneratorImplementation::
+ setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
+
+ const auto subcell_box =
+ internal::DiscreteQuadratureGeneratorImplementation::
+ create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
+
+ reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
+
+ FaceQuadratureGenerator<dim>::generate_append(
+ *reference_space_level_set, subcell_box, face_index);
+ }
+ }
+
+
+
template <int dim>
void
DiscreteFaceQuadratureGenerator<dim>::generate(
reference_space_level_set->set_active_cell(cell);
const BoundingBox<dim> unit_box = create_unit_bounding_box<dim>();
- FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
- unit_box,
- face_index);
+ if (reference_space_level_set->is_fe_q_iso_q1())
+ generate_fe_q_iso_q1(unit_box, face_index);
+ else
+ FaceQuadratureGenerator<dim>::generate(*reference_space_level_set,
+ unit_box,
+ face_index);
}
} // namespace NonMatching
#include "quadrature_generator.inst"
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
#include <deal.II/grid/grid_generator.h>
class Test
{
public:
- Test();
+ Test(const FE_Poly<dim> &fe);
void
run();
template <int dim>
-Test<dim>::Test()
+Test<dim>::Test(const FE_Poly<dim> &fe)
: dof_handler(triangulation)
{
- fe_collection.push_back(FE_Q<dim>(1));
- const unsigned int n_quadrature_points = 1;
+ const FE_Q_iso_Q1<dim> *fe_q_iso_q1 =
+ dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe);
+ const bool is_fe_q_iso_q1 = fe_q_iso_q1 != nullptr;
+ fe_collection.push_back(fe);
+
+ const unsigned int n_quadrature_points = is_fe_q_iso_q1 ? 1 : fe.get_degree();
q_collection1D.push_back(QGauss<1>(n_quadrature_points));
}
q_collection1D, dof_handler, level_set);
quadrature_generator.generate(triangulation.begin_active());
+ deallog << "inside" << std::endl;
print_quadrature(quadrature_generator.get_inside_quadrature());
+ deallog << "outside" << std::endl;
print_quadrature(quadrature_generator.get_outside_quadrature());
+ deallog << "surface" << std::endl;
print_surface_quadrature(quadrature_generator.get_surface_quadrature());
}
for (const auto f : cell->face_indices())
{
face_quadrature_generator.generate(cell, f);
-
+ deallog << "face_index = " << f << std::endl;
+ deallog << "inside" << std::endl;
print_quadrature(face_quadrature_generator.get_inside_quadrature());
+ deallog << "outside" << std::endl;
print_quadrature(face_quadrature_generator.get_outside_quadrature());
+ deallog << "surface" << std::endl;
print_surface_quadrature(
face_quadrature_generator.get_surface_quadrature());
}
template <int dim>
void
-run_test()
+run_test(const FE_Poly<dim> &fe)
{
deallog << "dim = " << dim << std::endl;
- Test<dim> test;
+ Test<dim> test(fe);
test.run();
deallog << std::endl;
}
{
initlog();
- run_test<1>();
- run_test<2>();
- run_test<3>();
+ std::shared_ptr<FE_Poly<1>> fe_1;
+ std::shared_ptr<FE_Poly<2>> fe_2;
+ std::shared_ptr<FE_Poly<3>> fe_3;
+ fe_1 = std::make_shared<FE_Q<1>>(1);
+ fe_2 = std::make_shared<FE_Q<2>>(1);
+ fe_3 = std::make_shared<FE_Q<3>>(1);
+ run_test<1>(*fe_1);
+ run_test<2>(*fe_2);
+ run_test<3>(*fe_3);
+ fe_1 = std::make_shared<FE_Q_iso_Q1<1>>(2);
+ fe_2 = std::make_shared<FE_Q_iso_Q1<2>>(2);
+ fe_3 = std::make_shared<FE_Q_iso_Q1<3>>(2);
+ run_test<1>(*fe_1);
+ run_test<2>(*fe_2);
+ run_test<3>(*fe_3);
}
DEAL::dim = 1
+DEAL::inside
DEAL::0.375000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.250000
+DEAL::surface
DEAL::0.750000, 1.00000, 1.00000
+DEAL::face_index = 0
+DEAL::inside
DEAL::1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
DEAL::1.00000
+DEAL::surface
DEAL::
DEAL::dim = 2
+DEAL::inside
DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
DEAL::0.500000, 1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
DEAL::0.500000, 1.00000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
DEAL::0.375000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.250000
+DEAL::surface
DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
DEAL::0.375000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.250000
+DEAL::surface
DEAL::0.750000, 1.00000, 1.00000, 0.00000
DEAL::
DEAL::dim = 3
+DEAL::inside
DEAL::0.375000, 0.500000, 0.500000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.500000, 0.500000, 0.250000
+DEAL::surface
DEAL::0.750000, 0.500000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
DEAL::0.500000, 0.500000, 1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
DEAL::0.500000, 0.500000, 1.00000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
DEAL::0.500000, 0.375000, 0.750000
+DEAL::outside
DEAL::0.500000, 0.875000, 0.250000
+DEAL::surface
DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
DEAL::0.500000, 0.375000, 0.750000
+DEAL::outside
DEAL::0.500000, 0.875000, 0.250000
+DEAL::surface
DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 4
+DEAL::inside
DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 5
+DEAL::inside
DEAL::0.375000, 0.500000, 0.750000
+DEAL::outside
DEAL::0.875000, 0.500000, 0.250000
+DEAL::surface
DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000
DEAL::
+DEAL::dim = 1
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::1.00000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::1.00000
+DEAL::surface
+DEAL::
+DEAL::dim = 2
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.750000, 0.500000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::0.250000, 0.500000
+DEAL::0.750000, 0.500000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
+DEAL::0.250000, 0.500000
+DEAL::0.625000, 0.250000
+DEAL::outside
+DEAL::0.875000, 0.250000
+DEAL::surface
+DEAL::0.750000, 1.00000, 1.00000, 0.00000
+DEAL::
+DEAL::dim = 3
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000, 0.125000
+DEAL::0.625000, 0.250000, 0.250000, 0.0625000
+DEAL::0.250000, 0.750000, 0.250000, 0.125000
+DEAL::0.625000, 0.750000, 0.250000, 0.0625000
+DEAL::0.250000, 0.250000, 0.750000, 0.125000
+DEAL::0.625000, 0.250000, 0.750000, 0.0625000
+DEAL::0.250000, 0.750000, 0.750000, 0.125000
+DEAL::0.625000, 0.750000, 0.750000, 0.0625000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.250000, 0.0625000
+DEAL::0.875000, 0.750000, 0.250000, 0.0625000
+DEAL::0.875000, 0.250000, 0.750000, 0.0625000
+DEAL::0.875000, 0.750000, 0.750000, 0.0625000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.250000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 0
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.750000, 0.750000, 0.250000
+DEAL::outside
+DEAL::surface
+DEAL::face_index = 1
+DEAL::inside
+DEAL::outside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.750000, 0.750000, 0.250000
+DEAL::surface
+DEAL::face_index = 2
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.625000, 0.125000
+DEAL::0.750000, 0.625000, 0.125000
+DEAL::outside
+DEAL::0.250000, 0.875000, 0.125000
+DEAL::0.750000, 0.875000, 0.125000
+DEAL::surface
+DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 3
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.750000, 0.250000, 0.250000
+DEAL::0.250000, 0.625000, 0.125000
+DEAL::0.750000, 0.625000, 0.125000
+DEAL::outside
+DEAL::0.250000, 0.875000, 0.125000
+DEAL::0.750000, 0.875000, 0.125000
+DEAL::surface
+DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 4
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::face_index = 5
+DEAL::inside
+DEAL::0.250000, 0.250000, 0.250000
+DEAL::0.625000, 0.250000, 0.125000
+DEAL::0.250000, 0.750000, 0.250000
+DEAL::0.625000, 0.750000, 0.125000
+DEAL::outside
+DEAL::0.875000, 0.250000, 0.125000
+DEAL::0.875000, 0.750000, 0.125000
+DEAL::surface
+DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000
+DEAL::
projects / dealii.git / commitdiff