-namespace
-{
- // Construct a quadrature formula containing the vertices of the reference
- // cell in dimension dim (with invalid weights)
- template <int dim>
- Quadrature<dim> &
- get_vertex_quadrature()
- {
- static Quadrature<dim> quad;
- if (quad.size() == 0)
- {
- std::vector<Point<dim>> points(GeometryInfo<dim>::vertices_per_cell);
- for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
- points[i] = GeometryInfo<dim>::unit_cell_vertex(i);
- quad = Quadrature<dim>(points);
- }
- return quad;
- }
-} // namespace
-
-
-
template <int dim, int spacedim, typename VectorType>
MappingFEField<dim, spacedim, VectorType, void>::MappingFEField(
const DoFHandler<dim, spacedim> &euler_dof_handler,
true))
, fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
, fe_values(this->euler_dof_handler->get_fe(),
- get_vertex_quadrature<dim>(),
+ ReferenceCell::get_nodal_type_quadrature<dim>(
+ this->euler_dof_handler->get_fe().reference_cell_type()),
update_values)
{
unsigned int size = 0;
true))
, fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
, fe_values(this->euler_dof_handler->get_fe(),
- get_vertex_quadrature<dim>(),
+ ReferenceCell::get_nodal_type_quadrature<dim>(
+ this->euler_dof_handler->get_fe().reference_cell_type()),
update_values)
{
unsigned int size = 0;
true))
, fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
, fe_values(this->euler_dof_handler->get_fe(),
- get_vertex_quadrature<dim>(),
+ ReferenceCell::get_nodal_type_quadrature<dim>(
+ this->euler_dof_handler->get_fe().reference_cell_type()),
update_values)
{
unsigned int size = 0;
, fe_mask(mapping.fe_mask)
, fe_to_real(mapping.fe_to_real)
, fe_values(mapping.euler_dof_handler->get_fe(),
- get_vertex_quadrature<dim>(),
+ ReferenceCell::get_nodal_type_quadrature<dim>(
+ this->euler_dof_handler->get_fe().reference_cell_type()),
update_values)
{}
*cell, euler_dof_handler);
Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell());
- AssertDimension(GeometryInfo<dim>::vertices_per_cell,
- fe_values.n_quadrature_points);
+ AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points);
AssertDimension(fe_to_real.size(),
euler_dof_handler->get_fe().n_components());
if (uses_level_dofs)
boost::container::small_vector<Point<spacedim>,
GeometryInfo<dim>::vertices_per_cell>
- vertices(GeometryInfo<dim>::vertices_per_cell);
+ vertices(cell->n_vertices());
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int comp = fe_to_real
typename VectorType::value_type value =
internal::ElementAccess<VectorType>::get(vector, dof_indices[i]);
if (euler_dof_handler->get_fe().is_primitive(i))
- for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
+ for (const unsigned int v : cell->vertex_indices())
vertices[v][comp] += fe_values.shape_value(i, v) * value;
else
Assert(false, ExcNotImplemented());
data.aux.resize(
dim - 1, std::vector<Tensor<1, spacedim>>(n_original_q_points));
+
+ // TODO: only a single reference cell type possible...
+ const auto reference_cell_type =
+ this->euler_dof_handler->get_fe().reference_cell_type();
+ const auto n_faces =
+ ReferenceCell::internal::Info::get_cell(reference_cell_type)
+ .n_faces();
+
// Compute tangentials to the unit cell.
- for (const unsigned int i : GeometryInfo<dim>::face_indices())
+ for (unsigned int i = 0; i < n_faces; ++i)
{
data.unit_tangentials[i].resize(n_original_q_points);
std::fill(data.unit_tangentials[i].begin(),
data.unit_tangentials[i].end(),
- GeometryInfo<dim>::unit_tangential_vectors[i][0]);
+ ReferenceCell::unit_tangential_vectors<dim>(
+ reference_cell_type, i, 0));
if (dim > 2)
{
- data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .resize(n_original_q_points);
- std::fill(
- data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .begin(),
- data.unit_tangentials[GeometryInfo<dim>::faces_per_cell + i]
- .end(),
- GeometryInfo<dim>::unit_tangential_vectors[i][1]);
+ data.unit_tangentials[n_faces + i].resize(
+ n_original_q_points);
+ std::fill(data.unit_tangentials[n_faces + i].begin(),
+ data.unit_tangentials[n_faces + i].end(),
+ ReferenceCell::unit_tangential_vectors<dim>(
+ reference_cell_type, i, 1));
}
}
}
// 0.
for (unsigned int d = 0; d != dim - 1; ++d)
{
- Assert(face_no + GeometryInfo<dim>::faces_per_cell * d <
+ Assert(face_no + cell->n_faces() * d <
data.unit_tangentials.size(),
ExcInternalError());
Assert(
data.aux[d].size() <=
- data
- .unit_tangentials[face_no +
- GeometryInfo<dim>::faces_per_cell * d]
- .size(),
+ data.unit_tangentials[face_no + cell->n_faces() * d].size(),
ExcInternalError());
mapping.transform(
make_array_view(
- data
- .unit_tangentials[face_no +
- GeometryInfo<dim>::faces_per_cell * d]),
+ data.unit_tangentials[face_no + cell->n_faces() * d]),
mapping_contravariant,
data,
make_array_view(data.aux[d]));
if (subface_no != numbers::invalid_unsigned_int)
{
+ // TODO
const double area_ratio =
GeometryInfo<dim>::subface_ratio(
cell->subface_case(face_no), subface_no);
initial_p_unit[d] = 0.5;
}
+ // TODO
initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
// for (unsigned int d=0; d<dim; ++d)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test MappingFEField for simplex meshes.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Solution : public Function<dim>
+{
+public:
+ Solution()
+ : Function<dim>(dim)
+ {}
+
+ double
+ value(const Point<dim> &point, const unsigned int compontent) const
+ {
+ return point[compontent];
+ }
+};
+
+void
+test()
+{
+ const int dim = 2;
+
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+ Simplex::FE_P<dim> fe(2);
+ FESystem<dim> euler_fe(fe, dim);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ DoFHandler<dim> euler_dof_handler(tria);
+ euler_dof_handler.distribute_dofs(euler_fe);
+
+ Vector<double> euler_vector(euler_dof_handler.n_dofs());
+
+ // TODO: not working (missing mapping)
+ // VectorTools::get_position_vector(euler_dof_handler, euler_vector);
+
+ MappingFE<dim> mapping_interpolation(Simplex::FE_P<dim>(1));
+ VectorTools::interpolate(mapping_interpolation,
+ euler_dof_handler,
+ Solution<dim>(),
+ euler_vector);
+
+ MappingFEField<dim> mapping(euler_dof_handler, euler_vector);
+
+ Simplex::QGauss<dim> quadrature_formula(1);
+
+ FEValues<dim> fe_values(mapping,
+ fe,
+ quadrature_formula,
+ update_values | update_gradients);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+ }
+
+ {
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+
+ Vector<double> solution(dof_handler.n_dofs());
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(mapping, 2);
+
+#if false
+ std::ofstream output("test.vtk");
+ data_out.write_vtk(output);
+#else
+ data_out.write_vtk(deallog.get_file_stream());
+#endif
+ }
+}
+
+int
+main()
+{
+ initlog();
+
+ test();
+}
\ No newline at end of file
projects / dealii.git / commitdiff