mapping_eulerian.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by Luca Heltai
//
// This file is part of the FSI-suite platform, based on the deal.II library.
//
// The FSI-suite platform 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 3.0 of the License,
// or (at your option) any later version. The full text of the license can be
// found in the file LICENSE at the top level of the FSI-suite platform
// distribution.
//
// ---------------------------------------------------------------------
 
#ifndef parsed_tools_mapping_fe_field_h
#define parsed_tools_mapping_fe_field_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/algorithms/general_data_storage.h>
 
#include <deal.II/base/function_parser.h>
#include <deal.II/base/parameter_acceptor.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/mapping_q_eulerian.h>
 
#include <deal.II/numerics/vector_tools.h>
 
namespace ParsedTools
{
  template <int dim, int spacedim = dim>
  class MappingEulerian : public dealii::ParameterAcceptor
  {
  public:
    MappingEulerian(
      const dealii::DoFHandler<dim, spacedim> &dh,
      const std::string                       &section_name              = "",
      const std::string           &initial_configuration_or_displacement = "",
      const bool                   use_displacement = false,
      const dealii::ComponentMask &mask             = dealii::ComponentMask());
 
    template <typename VectorType>
    void
    initialize(VectorType &configuration_or_displacement);
 
 
    template <typename VectorType>
    void
    initialize(VectorType &configuration_or_displacement,
               VectorType &locally_relevant_configuration_or_displacement);
 
    operator const dealii::Mapping<dim, spacedim> &() const;
 
    const dealii::Mapping<dim, spacedim> &
    operator()() const;
 
  private:
    const dealii::SmartPointer<const dealii::DoFHandler<dim, spacedim>>
      dof_handler;
 
    const dealii::ComponentMask mask;
 
    bool use_displacement;
 
    std::string initial_configuration_or_displacement_expression;
 
    std::unique_ptr<dealii::Mapping<dim, spacedim>> mapping;
  };
 
 
// Template implementation
#ifndef DOXYGEN
  template <int dim, int spacedim>
  template <typename VectorType>
  void
  MappingEulerian<dim, spacedim>::initialize(
    VectorType &configuration_or_displacement)
  {
    initialize(configuration_or_displacement, configuration_or_displacement);
  }
 
 
 
  template <int dim, int spacedim>
  template <typename VectorType>
  void
  MappingEulerian<dim, spacedim>::initialize(
    VectorType &configuration_or_displacement,
    VectorType &locally_relevant_configuration_or_displacement)
  {
    // Interpolate the user data
    if (initial_configuration_or_displacement_expression != "")
      {
        dealii::FunctionParser<spacedim> initial_configuration(
          initial_configuration_or_displacement_expression);
 
        AssertThrow(initial_configuration.n_components ==
                      dof_handler->get_fe().n_components(),
                    dealii::ExcMessage("The number of components in the "
                                       "configuration expression must be equal "
                                       "to the number of components in the "
                                       "fe."));
 
 
        dealii::VectorTools::interpolate(*dof_handler,
                                         initial_configuration,
                                         configuration_or_displacement,
                                         mask);
      }
    else
      {
        if (use_displacement)
          configuration_or_displacement = 0;
        else
          {
            if (dof_handler->get_triangulation()
                  .all_reference_cells_are_hyper_cube())
              dealii::VectorTools::get_position_vector(
                *dof_handler, configuration_or_displacement, mask);
            else
              {
                // Use interpolation to get the position vector.
                const std::string id[3] = {"x", "y", "z"};
                std::string       id_expression;
                std::string       sep = "";
                unsigned int      j   = 0;
                for (unsigned int i = 0;
                     i < dof_handler->get_fe().n_components();
                     ++i)
                  {
                    if (mask[i] && j < spacedim)
                      id_expression += sep + id[j++];
                    else
                      id_expression += sep + "0";
                    sep = "; ";
                  }
 
                dealii::FunctionParser<spacedim> initial_configuration(
                  id_expression);
 
                dealii::VectorTools::interpolate(*dof_handler,
                                                 initial_configuration,
                                                 configuration_or_displacement,
                                                 mask);
              }
          }
      }
    // Copy to the locally relevant vector
    locally_relevant_configuration_or_displacement =
      configuration_or_displacement;
 
    if (use_displacement)
      // Finallly initialize the mapping with our own configuration
      // vector.
      {
        AssertThrow(
          dof_handler->get_triangulation().all_reference_cells_are_hyper_cube(),
          dealii::ExcMessage("The displacement mapping is only "
                             "supported for hypercube grids."));
 
        AssertThrow(mask == dealii::ComponentMask(),
                    dealii::ExcMessage(
                      "Using the displacement is only possible with default "
                      "component mask."));
 
        const auto degree = dof_handler->get_fe().degree;
        mapping.reset(new dealii::MappingQEulerian<dim, VectorType, spacedim>(
          degree,
          *dof_handler,
          locally_relevant_configuration_or_displacement));
      }
    else
      {
        // Finallly initialize the mapping with our own configuration
        // vector.
        mapping.reset(new dealii::MappingFEField<dim, spacedim, VectorType>(
          *dof_handler, locally_relevant_configuration_or_displacement, mask));
      }
  }
#endif
 
} // namespace ParsedTools
#endif