From e6083994f9b578539a62645a22dadc822ffabe6f Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 24 Sep 2001 13:58:28 +0000
Subject: [PATCH] Add Eulerian mapping.

git-svn-id: https://svn.dealii.org/trunk@5038 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/fe/mapping_q1_eulerian.h  | 151 ++++++++++++++++++
 .../deal.II/source/fe/mapping_q1_eulerian.cc  | 109 +++++++++++++
 2 files changed, 260 insertions(+)
 create mode 100644 deal.II/deal.II/include/fe/mapping_q1_eulerian.h
 create mode 100644 deal.II/deal.II/source/fe/mapping_q1_eulerian.cc

diff --git a/deal.II/deal.II/include/fe/mapping_q1_eulerian.h b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h
new file mode 100644
index 0000000000..c9b24d0ae6
--- /dev/null
+++ b/deal.II/deal.II/include/fe/mapping_q1_eulerian.h
@@ -0,0 +1,151 @@
+//----------------------------  mapping_q1_eulerian.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors and Michael Stadler
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  mapping_q1_eulerian.h  ---------------------------
+#ifndef __deal2__mapping_q1_eulerian_h
+#define __deal2__mapping_q1_eulerian_h
+
+#include <base/smartpointer.h>
+#include <fe/mapping_q1.h> 
+
+
+/**
+ * Eulerian mapping of general unit cells by d-linear shape
+ * functions. Each cell is thus shifted in space by values given to
+ * the mapping through a finite element field.
+ *
+ * @sect3{Usage} 
+ *
+ * The constructor of this class takes two arguments: a reference to
+ * the vector that defines the mapping from the reference
+ * configuration to the current configuration and a reference to the
+ * @ref{DoFHandler}. The vector should then represent a (flattened out
+ * version of a) vector valued field defined at nodes defined by the 
+ * the @ref{DoFHandler}, where the number of components of the vector
+ * field equals the number of space dimensions. Thus, the 
+ * @ref{DoFHandler} shall operate on a finite element that has as many 
+ * components as space dimensions. As an additional requirement, we 
+ * impose that it have as many degree of freedom per vertex as there
+ * are space dimensions; since this object only evaluates the finite 
+ * element field at the vertices, the values
+ * of all other degrees of freedom (not associated to vertices) are
+ * ignored. These requirements are met if the finite element which the
+ * given @ref{DoFHandler} operates on is constructed as a system
+ * element (@ref{FESystem}) from @p{dim} continuous @ref{FE_Q}
+ * objects.
+ *
+ * In many cases, the shift vector will also be the solution vector of
+ * the problem under investigation. If this is not the case (i.e. the
+ * number of components of the solution variable is not equal to the
+ * space dimension, e.g. for scalar problems in @p{dim>1} where the
+ * Eulerian coordinates only give a background field) or for coupled
+ * problems where more variables are computed than just the flow
+ * field), then a different @ref{DoFHandler} has to be set up on the
+ * given triangulation, and the shift vector has then to be associated
+ * to it.
+ *
+ * An example is shown below:
+ * @begin{verbatim}
+ *    FESystem<dim> fe(FE_Q<dim>(1), dim);
+ *    DoFHandler<dim> flowfield_dof_handler(triangulation);
+ *    flowfield_dof_handler.distribute_dofs(fe);
+ *    Vector<double> map_points(flowfield_dof_handler.n_dofs());
+ *    MappingQ1Eulerian<dim> mymapping(map_points, flowfield_dof_handler);
+ * @end{verbatim}
+ *
+ * Note that since the vector of shift values and the dof handler are
+ * only associated to this object at construction time, you have to
+ * make sure that whenever you use this object, the given objects
+ * still represent valid data.
+ *
+ * @author Michael Stadler, 2001
+ */
+template <int dim>
+class MappingQ1Eulerian : public MappingQ1<dim>
+{
+  public:
+
+				     /**
+				      * Constructor. It takes a
+				      * @p{Vector<double> &} as its
+				      * first argument to specify the
+				      * transformation of the whole
+				      * problem from the reference to
+				      * the current configuration.
+				      * The organization of the
+				      * elements in the @p{Vector}
+				      * must follow the concept how
+				      * deal.II stores solutions that
+				      * are associated to a
+				      * triangulation.  This is
+				      * automatically the case if the
+				      * @p{Vector} represents the
+				      * solution of the previous step
+				      * of a nonlinear problem.
+				      * Alternatively, the @p{Vector}
+				      * can be initialized by
+				      * @p{DoFObjectAccessor::set_dof_values()}.
+				      */
+    MappingQ1Eulerian ( const Vector<double>  &euler_transform_vectors,
+		        const DoFHandler<dim> &shiftmap_dof_handler);
+
+
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcWrongNoOfComponents);
+    
+				     /**
+				      * Exception.
+				      */
+    DeclException2 (ExcWrongVectorSize,
+		    int, int,
+		    << "Vector has wrong size " << arg1
+		    << ", expected size " << arg2);
+
+				     /**
+				      * Exception.
+				      */
+    DeclException0 (ExcInactiveCell);
+
+
+
+  protected:
+
+				     /**
+				      * Reference to the vector of
+				      * shifts.
+				      */
+    const Vector<double> &euler_transform_vectors;
+    
+                                     /**
+				      * Pointer to the DoFHandler to
+				      * which the mapping vector is
+				      * associated.
+				      */
+    const SmartPointer<const DoFHandler<dim> > shiftmap_dof_handler;
+    
+
+  private:    
+				     /**
+				      * Computes the support points of
+				      * the mapping. For
+				      * @p{MappingQ1Eulerian} these
+				      * are the vertices.
+				      */
+    virtual void compute_mapping_support_points(
+      const typename Triangulation<dim>::cell_iterator &cell,
+      std::vector<Point<dim> > &a) const;
+    
+};
+
+
+#endif
diff --git a/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc
new file mode 100644
index 0000000000..2df84da7fe
--- /dev/null
+++ b/deal.II/deal.II/source/fe/mapping_q1_eulerian.cc
@@ -0,0 +1,109 @@
+//----------------------------  mapping_q1_eulerian.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors and Michael Stadler
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  mapping_q1_eulerian.cc  ---------------------------
+
+
+#include <fe/mapping_q1_eulerian.h>
+#include <lac/vector.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+
+
+
+template <int dim>
+MappingQ1Eulerian<dim>::MappingQ1Eulerian ( const Vector<double>  &euler_transform_vectors,
+					    const DoFHandler<dim> &shiftmap_dof_handler)
+                     :
+		     euler_transform_vectors(euler_transform_vectors),
+		     shiftmap_dof_handler(&shiftmap_dof_handler)
+{}
+
+
+
+template <int dim>
+void
+MappingQ1Eulerian<dim>::compute_mapping_support_points(
+  const typename Triangulation<dim>::cell_iterator &cell,
+  std::vector<Point<dim> > &a) const
+{
+
+				   // 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.
+  Assert (dim == shiftmap_dof_handler->get_fe().n_dofs_per_vertex(),
+          ExcWrongNoOfComponents());
+  Assert (shiftmap_dof_handler->get_fe().n_components() == dim,
+	  ExcWrongNoOfComponents());
+
+  Assert (shiftmap_dof_handler->n_dofs() == euler_transform_vectors.size(),
+          ExcWrongVectorSize(euler_transform_vectors.size(),
+			     shiftmap_dof_handler->n_dofs()));
+
+				   // 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>::cell_iterator
+    dof_cell (const_cast<Triangulation<dim> *> (&(cell->get_triangulation())),
+	      cell->level(),
+	      cell->index(),
+	      shiftmap_dof_handler);
+
+				   // We require the cell to be
+				   // active.  This is determined by
+				   // the user when looping over all
+				   // active cells in the problem
+				   // code.
+  Assert (dof_cell->active() == true,
+          ExcInactiveCell());
+
+				   // for Q1 elements, the number of
+				   // support points should equal the
+				   // number of vertices
+  a.resize(GeometryInfo<dim>::vertices_per_cell);
+
+				   // 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<dim> 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<dim; ++j)
+	shift_vector[j] = mapping_values(vertex_mapping[i]*dim+j);
+
+				       // compute new support point by
+				       // old (reference) value and
+				       // added shift
+      a[i] = cell->vertex(vertex_mapping[i]) + shift_vector;
+    }
+}
+
+
+// explicit instantiation
+template class MappingQ1Eulerian<deal_II_dimension>;
-- 
2.39.5