--- /dev/null
+Improved: FEImmersedSurfaceValues can now be constructed using a MappingFEField.
+<br>
+(Marco Feder, 2022/10/31)
+
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ virtual void
+ fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const override;
+
/**
* @}
*/
* Constructor. Gets cell-independent data from mapping and finite element
* objects, matching the quadrature rule and update flags.
*
- * @note Currently this class is only implemented for MappingCartesian and
- * MappingQ.
+ * @note Currently this class is only implemented for MappingCartesian,
+ * MappingQ and MappingFEField.
*/
FEImmersedSurfaceValues(const Mapping<dim> & mapping,
const FiniteElement<dim> & element,
}
+
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::fill_fe_immersed_surface_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const NonMatching::ImmersedSurfaceQuadrature<dim> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const
+{
+ AssertDimension(dim, spacedim);
+ Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
+ ExcInternalError());
+ const InternalData &data = static_cast<const InternalData &>(internal_data);
+
+ const unsigned int n_q_points = quadrature.size();
+
+ update_internal_dofs(cell, data);
+
+ internal::MappingFEFieldImplementation::
+ maybe_compute_q_points<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.quadrature_points);
+
+ internal::MappingFEFieldImplementation::
+ maybe_update_Jacobians<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real);
+
+ const UpdateFlags update_flags = data.update_each;
+ const std::vector<double> &weights = quadrature.get_weights();
+
+ if (update_flags & (update_normal_vectors | update_JxW_values))
+ {
+ AssertDimension(output_data.JxW_values.size(), n_q_points);
+
+ Assert(!(update_flags & update_normal_vectors) ||
+ (output_data.normal_vectors.size() == n_q_points),
+ ExcDimensionMismatch(output_data.normal_vectors.size(),
+ n_q_points));
+
+
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ const double det = data.contravariant[point].determinant();
+
+ // check for distorted cells.
+
+ // TODO: this allows for anisotropies of up to 1e6 in 3D and
+ // 1e12 in 2D. might want to find a finer
+ // (dimension-independent) criterion
+ Assert(det > 1e-12 * Utilities::fixed_power<dim>(
+ cell->diameter() / std::sqrt(double(dim))),
+ (typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
+ cell->center(), det, point)));
+
+ // The normals are n = J^{-T} * \hat{n} before normalizing.
+ Tensor<1, spacedim> normal;
+ for (unsigned int d = 0; d < spacedim; d++)
+ normal[d] =
+ data.covariant[point][d] * quadrature.normal_vector(point);
+
+ output_data.JxW_values[point] = weights[point] * det * normal.norm();
+
+ if ((update_flags & update_normal_vectors) != 0u)
+ {
+ normal /= normal.norm();
+ output_data.normal_vectors[point] = normal;
+ }
+ }
+
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_jacobians)
+ {
+ AssertDimension(output_data.jacobians.size(), n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ output_data.jacobians[point] = data.contravariant[point];
+ }
+
+ // copy values from InternalData to vector given by reference
+ if (update_flags & update_inverse_jacobians)
+ {
+ AssertDimension(output_data.inverse_jacobians.size(), n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ output_data.inverse_jacobians[point] =
+ data.covariant[point].transpose();
+ }
+
+ // calculate derivatives of the Jacobians
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_grads<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_grads);
+
+ // calculate derivatives of the Jacobians pushed forward to real cell
+ // coordinates
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_pushed_forward_grads<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_pushed_forward_grads);
+
+ // calculate hessians of the Jacobians
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_2nd_derivatives<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_2nd_derivatives);
+
+ // calculate hessians of the Jacobians pushed forward to real cell
+ // coordinates
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_pushed_forward_2nd_derivatives<dim,
+ spacedim,
+ VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_pushed_forward_2nd_derivatives);
+
+ // calculate gradients of the hessians of the Jacobians
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_3rd_derivatives<dim, spacedim, VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_3rd_derivatives);
+
+ // calculate gradients of the hessians of the Jacobians pushed forward to
+ // real cell coordinates
+ internal::MappingFEFieldImplementation::
+ maybe_update_jacobian_pushed_forward_3rd_derivatives<dim,
+ spacedim,
+ VectorType>(
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ euler_dof_handler->get_fe(),
+ fe_mask,
+ fe_to_real,
+ output_data.jacobian_pushed_forward_3rd_derivatives);
+ }
+}
+
namespace internal
{
namespace MappingFEFieldImplementation
--- /dev/null
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/tensor_function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/mapping_collection.h>
+
+#include <deal.II/non_matching/fe_immersed_values.h>
+#include <deal.II/non_matching/immersed_surface_quadrature.h>
+
+#include <deal.II/numerics/vector_tools_interpolate.h>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+using NonMatching::FEImmersedSurfaceValues;
+
+// Set up a triangulation that we can construct an
+// FEImmersedSurfaceValues object and compute known values when using a
+// MappingFEField.
+
+
+template <int dim>
+class Test
+{
+public:
+ Test();
+
+ void
+ run();
+
+private:
+ // Construct a triangulation with a single cell
+ void
+ setup_single_cell_triangulation();
+
+ // Make the ImmersedSurfaceQuadraturewith contain a single point in the center
+ // of the cell, with weight .5 and a normal different in all components.
+ void
+ setup_single_point_quadrature();
+
+ // Test that the constant n_quadrature_points corresponds to what we sent in.
+ void
+ test_n_quadrature_points();
+
+ // Test that we can use the get_quadrature function to get the quadrature
+ // we passed to the FEImmersedSurfaceValues constructor.
+ void
+ test_get_quadrature();
+
+ // Print the quadrature point in real space to make sure that it's mapped to
+ // real space correctly.
+ void
+ test_point_mapped_correctly();
+
+ // Print the normal in real space to make sure that it's mapped to real space
+ // correctly.
+ void
+ test_normal();
+
+ // Print JxW to check that the value is correct.
+ void
+ test_JxW();
+
+ // Test that we can call the shape_surface_grad function.
+ void
+ test_shape_surface_grad();
+
+ const FE_Bernstein<dim> element;
+ const FESystem<dim> fe_system;
+ Triangulation<dim> triangulation;
+ DoFHandler<dim> dof_handler;
+ DoFHandler<dim> euler_dof_handler;
+ Vector<double> euler_vector;
+ std::unique_ptr<MappingFEField<dim>> euler_mapping;
+ NonMatching::ImmersedSurfaceQuadrature<dim> quadrature;
+};
+
+
+
+template <int dim>
+Test<dim>::Test()
+ : element(1)
+ , fe_system(element, dim)
+ , dof_handler(triangulation)
+ , euler_dof_handler(triangulation)
+{}
+
+
+
+template <int dim>
+void
+Test<dim>::run()
+{
+ setup_single_cell_triangulation();
+
+ setup_single_point_quadrature();
+
+ test_n_quadrature_points();
+ test_get_quadrature();
+ test_point_mapped_correctly();
+ test_normal();
+ test_JxW();
+ test_shape_surface_grad();
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_single_cell_triangulation()
+{
+ const Point<dim> lower_left;
+ Point<dim> upper_right;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ upper_right[d] = 2 + d;
+ }
+
+ GridGenerator::hyper_rectangle(triangulation, lower_left, upper_right);
+ dof_handler.distribute_dofs(element);
+ euler_dof_handler.distribute_dofs(fe_system);
+ euler_vector.reinit(euler_dof_handler.n_dofs());
+
+ const ComponentMask mask(dim, true);
+ VectorTools::get_position_vector(euler_dof_handler, euler_vector, mask);
+ euler_mapping = std::make_unique<MappingFEField<dim>>(euler_dof_handler,
+ euler_vector,
+ mask);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::setup_single_point_quadrature()
+{
+ Point<dim> point;
+ const double weight = .5;
+ Tensor<1, dim> normal;
+ for (unsigned int d = 0; d < dim; ++d)
+ {
+ point[d] = .5;
+ normal[d] = d + 1;
+ }
+ normal /= normal.norm();
+ quadrature.push_back(point, weight, normal);
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_n_quadrature_points()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_default);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "n_quadrature_points = " << fe_values.n_quadrature_points
+ << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_get_quadrature()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_default);
+ fe_values.reinit(triangulation.begin_active());
+
+ const NonMatching::ImmersedSurfaceQuadrature<dim> &stored_quadrature =
+ fe_values.get_quadrature();
+
+ for (unsigned int q = 0; q < stored_quadrature.size(); q++)
+ {
+ deallog << "(point, weight, normal) = ([" << stored_quadrature.point(q)
+ << "], " << stored_quadrature.weight(q) << ", ["
+ << stored_quadrature.normal_vector(q) << "])" << std::endl;
+ }
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_point_mapped_correctly()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_quadrature_points);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "point = " << fe_values.quadrature_point(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_normal()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_normal_vectors);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "normal = " << fe_values.normal_vector(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_JxW()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_JxW_values);
+ fe_values.reinit(triangulation.begin_active());
+
+ deallog << "JxW = " << fe_values.JxW(0) << std::endl;
+}
+
+
+
+template <int dim>
+void
+Test<dim>::test_shape_surface_grad()
+{
+ FEImmersedSurfaceValues<dim> fe_values(*euler_mapping,
+ element,
+ quadrature,
+ update_gradients |
+ update_normal_vectors);
+ fe_values.reinit(dof_handler.begin_active());
+
+ const unsigned int function_index = 0;
+ const unsigned int q_index = 0;
+ const unsigned int component = 0;
+
+ deallog << "shape_surface_grad = "
+ << fe_values.shape_surface_grad(function_index, q_index) << std::endl;
+
+ deallog << "shape_surface_grad_component = "
+ << fe_values.shape_surface_grad_component(function_index,
+ q_index,
+ component)
+ << std::endl;
+}
+
+
+
+template <int dim>
+void
+run_test()
+{
+ deallog << "dim = " << dim << std::endl;
+ Test<dim> test;
+ test.run();
+ deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+ run_test<1>();
+ run_test<2>();
+ run_test<3>();
+}
--- /dev/null
+
+DEAL::dim = 1
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000], 0.500000, [1.00000])
+DEAL::point = 1.00000
+DEAL::normal = 1.00000
+DEAL::JxW = 0.500000
+DEAL::shape_surface_grad = 0.00000
+DEAL::shape_surface_grad_component = 0.00000
+DEAL::
+DEAL::dim = 2
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000], 0.500000, [0.447214 0.894427])
+DEAL::point = 1.00000 1.50000
+DEAL::normal = 0.600000 0.800000
+DEAL::JxW = 1.11803
+DEAL::shape_surface_grad = -0.0800000 0.0600000
+DEAL::shape_surface_grad_component = -0.0800000 0.0600000
+DEAL::
+DEAL::dim = 3
+DEAL::n_quadrature_points = 1
+DEAL::(point, weight, normal) = ([0.500000 0.500000 0.500000], 0.500000, [0.267261 0.534522 0.801784])
+DEAL::point = 1.00000 1.50000 2.00000
+DEAL::normal = 0.445976 0.594635 0.668965
+DEAL::JxW = 3.59563
+DEAL::shape_surface_grad = -0.0593923 0.00414365 0.0359116
+DEAL::shape_surface_grad_component = -0.0593923 0.00414365 0.0359116
+DEAL::
projects / dealii.git / commitdiff