]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove a few 1D specializations for mapping. 572/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 19 Feb 2015 09:50:55 +0000 (10:50 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 19 Feb 2015 10:18:19 +0000 (11:18 +0100)
The old implementation of MappingQEulerian only moved the vertices
in 1D, not interior points for higher order mappings. While it might
be debatable how useful non-linear mappings are in 1D, this was not
documented and lead to surprising behavior. By making the class more
general we can actually make our code base simpler because we can remove
a few template specialization and use the same code as in higher
dimensions.

doc/news/changes.h
include/deal.II/fe/mapping_q.h
include/deal.II/grid/tria_boundary.h
source/fe/mapping_q.cc
source/grid/tria_boundary.cc
tests/bits/mapping_q_eulerian_06.cc [new file with mode: 0644]
tests/bits/mapping_q_eulerian_06.output [new file with mode: 0644]

index 23518e1e709fadbe0af1d438b5101ec13d434a45..77ca21aa83cfd8c203cc94958128e10587124c2a 100644 (file)
@@ -350,6 +350,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: MappingQEulerian would previously not move interior points
+  in 1D for higher order mappings. This has been fixed by removing a few
+  specializations of MappingQ for 1D that are no longer necessary.
+  <br>
+  (Martin Kronbichler, 2015/02/19)
+  </li>
+
   <li> Fixed: The implementation of the class GrowingVectorMemory has been
   moved from source/lac/vector_memory.cc to the new file
   include/deal.II/lac/vector_memory.templates.h. This allows users to
index 87a8e628eeab15c897613a2e79acdca9dccc9719..9f8957068b84aec830d92caed8e857521f46f966 100644 (file)
@@ -487,13 +487,6 @@ private:
 
 #ifndef DOXYGEN
 
-template<> MappingQ<1>::MappingQ (const unsigned int,
-                                  const bool);
-template<> MappingQ<1>::~MappingQ ();
-
-template<>
-void MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
-                                          MappingQ1<1>::InternalData   &data) const;
 template <>
 void MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const;
 
index c5d59a4ddf4aaf45ca2af3921d0fa20129faac62..145d69461d2fd88cd7625bbc693f9f2623b3794b 100644 (file)
@@ -465,12 +465,6 @@ Point<3>
 StraightBoundary<3,3>::
 get_new_point_on_quad (const Triangulation<3,3>::quad_iterator &quad) const;
 
-template <>
-void
-StraightBoundary<1,1>::
-get_intermediate_points_on_line (const Triangulation<1,1>::line_iterator &,
-                                 std::vector<Point<1> > &) const;
-
 template <>
 void
 StraightBoundary<3,3>::
index 43bc3f7814f39ed689147375dc094394bea97478..79d53e46d9d5b31ac10d2ddef8e46b0287eb4434 100644 (file)
@@ -57,44 +57,6 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 
 
 
-// in 1d, it is irrelevant which polynomial degree to use, since all
-// cells are scaled linearly. Unless codimension is equal to two
-template<>
-MappingQ<1>::MappingQ (const unsigned int,
-                       const bool /*use_mapping_q_on_all_cells*/)
-  :
-  degree(1),
-  n_inner(0),
-  n_outer(0),
-  tensor_pols(0),
-  n_shape_functions(2),
-  renumber(0),
-  use_mapping_q_on_all_cells (false),
-  feq(degree),
-  line_support_points(degree+1)
-{}
-
-
-template<>
-MappingQ<1>::MappingQ (const MappingQ<1> &m):
-  MappingQ1<1> (),
-  degree(1),
-  n_inner(0),
-  n_outer(0),
-  tensor_pols(0),
-  n_shape_functions(2),
-  renumber(0),
-  use_mapping_q_on_all_cells (m.use_mapping_q_on_all_cells),
-  feq(degree),
-  line_support_points(degree+1)
-{}
-
-template<>
-MappingQ<1>::~MappingQ ()
-{}
-
-
-
 namespace
 {
   template <int dim>
@@ -179,16 +141,6 @@ MappingQ<dim,spacedim>::~MappingQ ()
 
 
 
-template<>
-void
-MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
-                                     MappingQ1<1>::InternalData   &data) const
-{
-  MappingQ1<1>::compute_shapes_virtual(unit_points, data);
-}
-
-
-
 template<int dim, int spacedim>
 void
 MappingQ<dim,spacedim>::compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
index e32916554fcee4103bb6c7941e43a2f0095de74a..e68d058cd586370bbf70b6ccb37e2818ea977a3c 100644 (file)
@@ -335,42 +335,6 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
 
 
 
-template <>
-void
-StraightBoundary<1>::
-get_intermediate_points_on_line (const Triangulation<1>::line_iterator &,
-                                 std::vector<Point<1> > &) const
-{
-  Assert(false, ExcImpossibleInDim(1));
-}
-
-template <>
-void
-StraightBoundary<1, 2>::
-get_intermediate_points_on_line (const Triangulation<1, 2>::line_iterator &line,
-                                 std::vector<Point<2> > &points) const
-{
-  const unsigned int spacedim = 2;
-  const unsigned int n=points.size();
-  Assert(n>0, ExcInternalError());
-
-  // Use interior points of QGaussLobatto quadrature formula support points
-  // for consistency with MappingQ
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-  const Point<spacedim> vertices[2] = { line->vertex(0),
-                                        line->vertex(1)
-                                      };
-
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[i+1][0];
-      points[i] = (1-x)*vertices[0] + x*vertices[1];
-    }
-}
-
-
-
-
 template <int dim, int spacedim>
 void
 StraightBoundary<dim, spacedim>::
diff --git a/tests/bits/mapping_q_eulerian_06.cc b/tests/bits/mapping_q_eulerian_06.cc
new file mode 100644 (file)
index 0000000..4d627e1
--- /dev/null
@@ -0,0 +1,91 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// computes points in real space for 1D Eulerian mapping where middle points
+// are moved
+
+#include "../tests.h"
+#include <fstream>
+#include <deal.II/base/logstream.h>
+
+// all include files you need here
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <fstream>
+#include <string>
+
+std::ofstream logfile("output");
+
+void test(unsigned int degree)
+{
+  const unsigned int dim = 1;
+  const unsigned int spacedim = 1;
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::hyper_cube(tria, 0, 1);
+  FE_Q<dim, spacedim> fe(degree);
+
+  DoFHandler<dim, spacedim> shift_dh(tria);
+
+  shift_dh.distribute_dofs(fe);
+
+  // Shift just the interior points but not the boundary points
+  Vector<double> shift(shift_dh.n_dofs());
+  for (unsigned int i=2; i<=degree; ++i)
+    shift(i) = 0.1;
+
+  QGauss<dim> quad(degree+1);
+  MappingQEulerian<dim,Vector<double>,spacedim> mapping(degree,shift, shift_dh);
+
+  Triangulation<dim,spacedim>::active_cell_iterator cell=tria.begin_active(),
+                                                    endc=tria.end() ;
+  Point<spacedim> real;
+  Point<dim> unit;
+  double eps = 1e-10;
+  for (; cell!=endc; ++cell)
+    {
+      deallog<<cell<< std::endl;
+      for (unsigned int q=0; q<quad.size(); ++q)
+        {
+          real = mapping.transform_unit_to_real_cell(cell, quad.point(q));
+          unit = mapping.transform_real_to_unit_cell(cell, real);
+          deallog<<quad.point(q)<< " -> " << real << std::endl;
+          if ( (unit-quad.point(q)).norm()>eps)
+            deallog<<quad.point(q)<< " != " << unit << std::endl;
+        }
+    }
+
+}
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test(1);
+  test(2);
+  test(3);
+  return 0;
+}
+
diff --git a/tests/bits/mapping_q_eulerian_06.output b/tests/bits/mapping_q_eulerian_06.output
new file mode 100644 (file)
index 0000000..780c1cb
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL::0.0
+DEAL::0.211325 -> 0.211325
+DEAL::0.788675 -> 0.788675
+DEAL::0.0
+DEAL::0.112702 -> 0.152702
+DEAL::0.500000 -> 0.600000
+DEAL::0.887298 -> 0.927298
+DEAL::0.0
+DEAL::0.0694318 -> 0.0985068
+DEAL::0.330009 -> 0.429506
+DEAL::0.669991 -> 0.769487
+DEAL::0.930568 -> 0.959643

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.