mapping_q1_eulerian.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2018 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#include <deal.II/base/std_cxx14/memory.h>

#include <deal.II/fe/mapping_q1_eulerian.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/petsc_parallel_block_vector.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/fe/fe.h>


#include <array>

DEAL_II_NAMESPACE_OPEN


template <int dim, class VectorType, int spacedim>
MappingQ1Eulerian<dim, VectorType, spacedim>::
MappingQ1Eulerian (const DoFHandler<dim,spacedim> &shiftmap_dof_handler,
                   const VectorType               &euler_transform_vectors)
  :
  MappingQGeneric<dim,spacedim>(1),
  euler_transform_vectors(&euler_transform_vectors),
  shiftmap_dof_handler(&shiftmap_dof_handler)
{}



template <int dim, class VectorType, int spacedim>
std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
MappingQ1Eulerian<dim, VectorType, spacedim>::
get_vertices
(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
{
  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
  // The assertions can not be in the constructor, since this would
  // require to call dof_handler.distribute_dofs(fe) *before* the mapping
  // object is constructed, which is not necessarily what we want.

  //TODO: Only one of these two assertions should be relevant
  AssertDimension (spacedim, shiftmap_dof_handler->get_fe().n_dofs_per_vertex());
  AssertDimension (shiftmap_dof_handler->get_fe(0).n_components(), spacedim);

  AssertDimension (shiftmap_dof_handler->n_dofs(), euler_transform_vectors->size());

  // cast the Triangulation<dim>::cell_iterator into a
  // DoFHandler<dim>::cell_iterator which is necessary for access to
  // DoFCellAccessor::get_dof_values()
  typename DoFHandler<dim,spacedim>::cell_iterator dof_cell (*cell, shiftmap_dof_handler);

  // We require the cell to be active since we can only then get nodal
  // values for the shifts
  Assert (dof_cell->active() == true, ExcInactiveCell());

  // now get the values of the shift vectors at the vertices
  Vector<double> mapping_values (shiftmap_dof_handler->get_fe().dofs_per_cell);
  dof_cell->get_dof_values (*euler_transform_vectors, mapping_values);

  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
    {
      Point<spacedim> shift_vector;

      // pick out the value of the shift vector at the present
      // vertex. since vertex dofs are always numbered first, we can
      // access them easily
      for (unsigned int j=0; j<spacedim; ++j)
        shift_vector[j] = mapping_values(i*spacedim+j);

      // compute new support point by old (reference) value and added
      // shift
      vertices[i] = cell->vertex(i) + shift_vector;
    }
  return vertices;
}



template <int dim, class VectorType, int spacedim>
std::vector<Point<spacedim> >
MappingQ1Eulerian<dim,VectorType,spacedim>::
compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
{
  const std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
  vertices = this->get_vertices(cell);

  std::vector<Point<spacedim> > a(GeometryInfo<dim>::vertices_per_cell);
  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
    a[i] = vertices[i];

  return a;
}





template <int dim, class VectorType, int spacedim>
std::unique_ptr<Mapping<dim,spacedim> >
MappingQ1Eulerian<dim, VectorType, spacedim>::clone () const
{
  return std_cxx14::make_unique<MappingQ1Eulerian<dim,VectorType,spacedim>>(*this);
}



template <int dim, class VectorType, int spacedim>
CellSimilarity::Similarity
MappingQ1Eulerian<dim,VectorType,spacedim>::
fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                const CellSimilarity::Similarity,
                const Quadrature<dim>                                     &quadrature,
                const typename Mapping<dim,spacedim>::InternalDataBase    &internal_data,
                internal::FEValuesImplementation::MappingRelatedData<dim,spacedim>      &output_data) const
{
  // call the function of the base class, but ignoring
  // any potentially detected cell similarity between
  // the current and the previous cell
  MappingQGeneric<dim,spacedim>::fill_fe_values (cell,
                                                 CellSimilarity::invalid_next_cell,
                                                 quadrature,
                                                 internal_data,
                                                 output_data);
  // also return the updated flag since any detected
  // similarity wasn't based on the mapped field, but
  // the original vertices which are meaningless
  return CellSimilarity::invalid_next_cell;
}



// explicit instantiations
#include "mapping_q1_eulerian.inst"


DEAL_II_NAMESPACE_CLOSE