]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Allow construction of FEImmersedSurfaceValues with a MappingFEField 14380/head
authorMarco Feder <marco.feder@sissa.it>
Sat, 29 Oct 2022 16:51:30 +0000 (18:51 +0200)
committerMarco Feder <marco.feder@sissa.it>
Tue, 1 Nov 2022 08:42:35 +0000 (09:42 +0100)
doc/news/changes/minor/20221031Feder [new file with mode: 0644]
include/deal.II/fe/mapping_fe_field.h
include/deal.II/non_matching/fe_immersed_values.h
source/fe/mapping_fe_field.cc
tests/non_matching/fe_immersed_surface_values_02.cc [new file with mode: 0644]
tests/non_matching/fe_immersed_surface_values_02.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20221031Feder b/doc/news/changes/minor/20221031Feder
new file mode 100644 (file)
index 0000000..a01350b
--- /dev/null
@@ -0,0 +1,4 @@
+Improved: FEImmersedSurfaceValues can now be constructed using a MappingFEField.
+<br>
+(Marco Feder, 2022/10/31)
+
index 9562157930e624f5dc2e7e80ddd0b112e37c4929..676369469ddef1c5b90afa229b984a784d1a6cce 100644 (file)
@@ -534,6 +534,14 @@ protected:
     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;
+
   /**
    * @}
    */
index 8ee37ac395ee5cadd9df343ee1e255db4752018d..48897206785265a035f39b3d08d70617a9242cfa 100644 (file)
@@ -60,8 +60,8 @@ namespace NonMatching
      * 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,
index 749783d9dcf5e1adffab24d8588eabb586cd2f66..58be345c395c169dbf618e91775744c0b549a2da 100644 (file)
@@ -1771,6 +1771,170 @@ MappingFEField<dim, spacedim, VectorType>::fill_fe_subface_values(
 }
 
 
+
+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
diff --git a/tests/non_matching/fe_immersed_surface_values_02.cc b/tests/non_matching/fe_immersed_surface_values_02.cc
new file mode 100644 (file)
index 0000000..3affc77
--- /dev/null
@@ -0,0 +1,292 @@
+#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>();
+}
diff --git a/tests/non_matching/fe_immersed_surface_values_02.output b/tests/non_matching/fe_immersed_surface_values_02.output
new file mode 100644 (file)
index 0000000..4d4f0d8
--- /dev/null
@@ -0,0 +1,28 @@
+
+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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.