From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Wed, 22 Jun 2016 16:12:22 +0000 (+0200)
Subject: Place new points for MappingQ along manifold
X-Git-Tag: v8.5.0-rc1~958^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=49f3a01ee6bc39bfd1983d2285caf10171713545;p=dealii.git

Place new points for MappingQ along manifold
---

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 6f82615823..0ea518fb13 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -3479,7 +3479,7 @@ namespace
       {
       case 2:
       {
-        // If two points are passed, these are the two vertices, and
+        // If two points are passed, these are the two vertices, so
         // we can only compute degree-1 intermediate points.
         for (unsigned int i=0; i<n; ++i)
           {
@@ -3500,9 +3500,8 @@ namespace
         // is n a square number
         Assert(m*m==n, ExcInternalError());
 
-        // If four points are passed, these are the two vertices, and
-        // we can only compute (degree-1)*(degree-1) intermediate
-        // points.
+        // If m points are passed, two of them are the outer vertices, and we
+        // can only compute (degree-1)*(degree-1) intermediate points.
         for (unsigned int i=0; i<m; ++i)
           {
             const double y=line_support_points.point(1+i)[0];
@@ -3522,8 +3521,42 @@ namespace
       }
 
       case 8:
-        Assert(false, ExcNotImplemented());
+      {
+        Assert(spacedim >= 3, ExcImpossibleInDim(spacedim));
+        unsigned int m=1;
+        for ( ; m < n; ++m)
+          if (m*m*m == n)
+            break;
+        // is n a cube number
+        Assert(m*m*m==n, ExcInternalError());
+
+        // If m points are passed, two of them are the outer vertices, and we can
+        // only compute (degree-1)*(degree-1)*(degree-1) intermediate points.
+        for (unsigned int k=0, c=0; k<m; ++k)
+          {
+            const double z=line_support_points.point(1+k)[0];
+            for (unsigned int i=0; i<m; ++i)
+              {
+                const double y=line_support_points.point(1+i)[0];
+                for (unsigned int j=0; j<m; ++j, ++c)
+                  {
+                    const double x=line_support_points.point(1+j)[0];
+
+                    w[0] = (1-x)*(1-y)*(1-z);
+                    w[1] =     x*(1-y)*(1-z);
+                    w[2] = (1-x)*y    *(1-z);
+                    w[3] =     x*y    *(1-z);
+                    w[4] = (1-x)*(1-y)*z    ;
+                    w[5] =     x*(1-y)*z    ;
+                    w[6] = (1-x)*y    *z    ;
+                    w[7] =     x*y    *z    ;
+                    Quadrature<spacedim> quadrature(surrounding_points, w);
+                    points[c]=manifold.get_new_point(quadrature);
+                  }
+              }
+          }
         break;
+      }
       default:
         Assert(false, ExcInternalError());
         break;
@@ -3737,9 +3770,10 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
                ExcInternalError());
 #endif
 
-      // if face at boundary, then ask boundary object to return intermediate
-      // points on it
-      if (face->at_boundary())
+      // if face at boundary or if it belongs to a manifold, then ask the
+      // boundary/manifold object to return intermediate points on it
+      if (face->at_boundary() ||
+          dynamic_cast<const Boundary<3,3> *>(&face->get_manifold()) == NULL)
         {
           get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points);
 
@@ -3877,6 +3911,13 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
         // manifold, otherwise compute them from the points around it
         if (dim != spacedim)
           add_quad_support_points(cell, a);
+        else if (dynamic_cast<const Boundary<dim,dim> *>(&cell->get_manifold()) == NULL)
+          {
+            std::vector<Point<spacedim> > quad_points (Utilities::fixed_power<dim>(polynomial_degree-1));
+            get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points);
+            for (unsigned int i=0; i<quad_points.size(); ++i)
+              a.push_back(quad_points[i]);
+          }
         else
           add_weighted_interior_points (support_point_weights_on_quad, a);
         break;
@@ -3888,7 +3929,15 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
         add_quad_support_points (cell, a);
 
         // then compute the interior points
-        add_weighted_interior_points (support_point_weights_on_hex, a);
+        if (dynamic_cast<const Boundary<dim,dim> *>(&cell->get_manifold()) == NULL)
+          {
+            std::vector<Point<spacedim> > hex_points (Utilities::fixed_power<dim>(polynomial_degree-1));
+            get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points);
+            for (unsigned int i=0; i<hex_points.size(); ++i)
+              a.push_back(hex_points[i]);
+          }
+        else
+          add_weighted_interior_points (support_point_weights_on_hex, a);
         break;
       }
 
diff --git a/tests/fe/mapping_q_manifold_02.cc b/tests/fe/mapping_q_manifold_02.cc
new file mode 100644
index 0000000000..52ca3acf7f
--- /dev/null
+++ b/tests/fe/mapping_q_manifold_02.cc
@@ -0,0 +1,242 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+// Test that high order MappingQ on a manifold puts the intermediate points
+// along the manifold
+
+#include "../tests.h"
+
+#include <fstream>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/manifold.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+using namespace dealii;
+
+double f_x(double x_m)
+{
+  double x = x_m*1000.0;
+  double y_m = 0.0;
+
+  if (x <= 9.0)
+    y_m = 0.001*(-28.0 + std::min(28.0, 2.8e1 + 6.775070969851e-3*x*x - 2.124527775800e-3*x*x*x));
+
+  else if (x > 9.0 && x <= 14.0)
+    y_m = 0.001*(-28.0 + 2.507355893131e1 + 9.754803562315e-1*x - 1.016116352781e-1*x*x + 1.889794677828e-3*x*x*x);
+
+  else if (x > 14.0 && x <= 20.0)
+    y_m = 0.001*(-28.0 + 2.579601052357e1 + 8.206693007457e-1*x - 9.055370274339e-2*x*x + 1.626510569859e-3*x*x*x);
+
+  else if (x > 20.0 && x <= 30.0)
+    y_m = 0.001*(-28.0 + 4.046435022819e1 - 1.379581654948*x + 1.945884504128e-2*x*x - 2.070318932190e-4*x*x*x);
+
+  else if (x > 30.0 && x <= 40.0)
+    y_m = 0.001*(-28.0 + 1.792461334664e1 + 8.743920332081e-1*x - 5.567361123058e-2*x*x + 6.277731764683e-4*x*x*x);
+
+  else if (x > 40.0 && x <= 54.0)
+    y_m = 0.001*(-28.0 + std::max(0.0, 5.639011190988e1 - 2.010520359035*x + 1.644919857549e-2*x*x + 2.674976141766e-5*x*x*x));
+
+  else if (x > 54.0)
+    y_m = 0.001*(-28.0);
+  return y_m;
+}
+
+template<int dim>
+class PushForward : public Function<dim>
+{
+public:
+  PushForward () : Function<dim>(dim, 0.) {} //,component(component) {}
+
+  virtual ~PushForward() {};
+
+  virtual double value (const Point<dim> &p,const unsigned int component = 0) const;
+
+private:
+  const double h = 0.028;
+
+  // data from initial block
+  const double x_max = 4.5*h;//9.0*h;
+  const double y_max = 2.036*h;
+  const double z_max = 4.5*h;
+
+  const double y_FoR = h;
+
+
+};
+
+template<int dim>
+double PushForward<dim>::value(const Point<dim> &p,const unsigned int component) const
+{
+  double result = 0;
+
+  // x component
+  if (component == 0)
+    result = p[0];
+
+  // y component
+  else if (component == 1)
+    {
+      if (p[0] <= x_max/2.0)
+        result = p[1] + (1 - (p[1] - y_FoR)/y_max)*f_x(p[0]);
+
+      else if (p[0] > x_max/2.0)
+        result = p[1] + (1 - (p[1] - y_FoR)/y_max)*f_x(x_max - p[0]);
+    }
+
+  // z component
+  else if (component == 2)
+    result = p[2];
+
+  return result;
+}
+
+template<int dim>
+class PullBack : public Function<dim>
+{
+public:
+  PullBack () : Function<dim>(dim, 0.) {} //,component(component) {}
+
+  virtual ~PullBack() {};
+
+  virtual double value (const Point<dim> &p,const unsigned int component = 0) const;
+
+private:
+  const double h = 0.028;
+  const double x_max = 4.5*h;//9.0*h;
+  const double y_max = 2.036*h;
+  const double z_max = 4.5*h;
+
+  const double y_FoR = h;
+};
+
+template<int dim>
+double PullBack<dim>::value(const Point<dim> &p,const unsigned int component) const
+{
+  double result = 0;
+
+  // x component
+  if (component == 0)
+    result = p[0];
+
+  // y component
+  else if (component == 1)
+    {
+      if (p[0] <= x_max/2.0)
+        result = (p[1] - f_x(p[0])*(1 + y_FoR/y_max))/(1.0 - f_x(p[0])/y_max);
+      else if (p[0] > x_max/2.0)
+        result = (p[1] - f_x(x_max - p[0])*(1 + y_FoR/y_max))/(1.0 - f_x(x_max - p[0])/y_max);
+    }
+
+  // z component
+  else if (component == 2)
+    result = p[2];
+
+  return result;
+}
+
+template <int dim>
+void create_tria(Triangulation<dim> &triangulation, const Manifold<dim> &manifold)
+{
+  const double h = 0.028;
+  std::vector<unsigned int> refinements(dim,1);
+  refinements[1] = 2;
+
+  Point<dim> p1, p2;
+  p2[0] = 4.5*h;//9.0*h;
+  p1[1] = h;
+  p2[1] = 2.018*h;//2.036*h;
+  if (dim == 3)
+    {
+      p1[2] = -2.25*h;
+      p2[2] = 2.25*h;
+    }
+  GridGenerator::hyper_rectangle(triangulation, p1, p2);
+
+  for (typename Triangulation<dim>::cell_iterator cell = triangulation.begin();
+       cell != triangulation.end(); ++cell)
+    cell->set_all_manifold_ids(111);
+  triangulation.set_manifold(111, manifold);
+
+  triangulation.refine_global(1);
+}
+
+template <int dim>
+void test()
+{
+  deallog << "dim: " << dim << std::endl;
+
+  PushForward<dim> push_forward;
+  PullBack<dim> pull_back;
+  FunctionManifold<dim,dim,dim> manifold(push_forward, pull_back);
+
+  Triangulation<dim> triangulation;
+  create_tria(triangulation, manifold);
+
+  for (unsigned mapping_p = 4; mapping_p <= 5; ++mapping_p)
+    {
+      deallog << "Mapping degree: " << mapping_p << std::endl;
+      MappingQ<dim> mapping(mapping_p, true);
+      std::vector<Point<dim> > points(Utilities::fixed_power<dim>(2));
+      for (unsigned int i=0, c=0; i<(dim==2?1:2); ++i)
+        for (unsigned int j=0; j<2; ++j)
+          for (unsigned int k=0; k<2; ++k, ++c)
+            {
+              points[c][0] = 0.1 + 0.8*k;
+              points[c][1] = 0.1 + 0.8*j;
+              if (dim == 3)
+                points[c][2] = 0.1 + 0.8*i;
+            }
+      FE_Nothing<dim> dummy;
+      Quadrature<dim> quad(points);
+      FEValues<dim> fe_values(mapping, dummy, quad, update_quadrature_points);
+
+      for (typename Triangulation<dim>::active_cell_iterator
+           cell=triangulation.begin_active(); cell != triangulation.end(); ++cell)
+        {
+          fe_values.reinit(cell);
+          deallog << "Points for cell with first point: " << cell->vertex(0)
+                  << std::endl;
+          for (unsigned int q=0; q<quad.size(); ++q)
+            deallog << fe_values.quadrature_point(q) << "   ";
+          deallog << std::endl;
+        }
+
+      deallog << std::endl;
+    }
+}
+
+int main ()
+{
+  deallog << std::setprecision (5);
+  deallog.attach (std::cout);
+  deallog.depth_console (0);
+  deallog.threshold_double (1e-12);
+
+  test<2>();
+  test<3>();
+}
diff --git a/tests/fe/mapping_q_manifold_02.output b/tests/fe/mapping_q_manifold_02.output
new file mode 100644
index 0000000000..93aba769c5
--- /dev/null
+++ b/tests/fe/mapping_q_manifold_02.output
@@ -0,0 +1,59 @@
+
+DEAL::dim: 2
+DEAL::Mapping degree: 4
+DEAL::Points for cell with first point: 0.0000 0.028000
+DEAL::0.0063000 0.029459   0.056700 0.0018472   0.0063000 0.040854   0.056700 0.018906   
+DEAL::Points for cell with first point: 0.063000 0.0000
+DEAL::0.069300 0.0018472   0.11970 0.029459   0.069300 0.018906   0.11970 0.040854   
+DEAL::Points for cell with first point: 0.0000 0.042252
+DEAL::0.0063000 0.043703   0.056700 0.023171   0.0063000 0.055097   0.056700 0.040229   
+DEAL::Points for cell with first point: 0.063000 0.021252
+DEAL::0.069300 0.023171   0.11970 0.043703   0.069300 0.040229   0.11970 0.055097   
+DEAL::
+DEAL::Mapping degree: 5
+DEAL::Points for cell with first point: 0.0000 0.028000
+DEAL::0.0063000 0.029303   0.056700 0.0020929   0.0063000 0.040730   0.056700 0.019101   
+DEAL::Points for cell with first point: 0.063000 0.0000
+DEAL::0.069300 0.0020929   0.11970 0.029303   0.069300 0.019101   0.11970 0.040730   
+DEAL::Points for cell with first point: 0.0000 0.042252
+DEAL::0.0063000 0.043586   0.056700 0.023353   0.0063000 0.055013   0.056700 0.040361   
+DEAL::Points for cell with first point: 0.063000 0.021252
+DEAL::0.069300 0.023353   0.11970 0.043586   0.069300 0.040361   0.11970 0.055013   
+DEAL::
+DEAL::dim: 3
+DEAL::Mapping degree: 4
+DEAL::Points for cell with first point: 0.0000 0.028000 -0.063000
+DEAL::0.0063000 0.029459 -0.056700   0.056700 0.0018472 -0.056700   0.0063000 0.040854 -0.056700   0.056700 0.018906 -0.056700   0.0063000 0.029459 -0.0063000   0.056700 0.0018472 -0.0063000   0.0063000 0.040854 -0.0063000   0.056700 0.018906 -0.0063000   
+DEAL::Points for cell with first point: 0.063000 0.0000 -0.063000
+DEAL::0.069300 0.0018472 -0.056700   0.11970 0.029459 -0.056700   0.069300 0.018906 -0.056700   0.11970 0.040854 -0.056700   0.069300 0.0018472 -0.0063000   0.11970 0.029459 -0.0063000   0.069300 0.018906 -0.0063000   0.11970 0.040854 -0.0063000   
+DEAL::Points for cell with first point: 0.0000 0.042252 -0.063000
+DEAL::0.0063000 0.043703 -0.056700   0.056700 0.023171 -0.056700   0.0063000 0.055097 -0.056700   0.056700 0.040229 -0.056700   0.0063000 0.043703 -0.0063000   0.056700 0.023171 -0.0063000   0.0063000 0.055097 -0.0063000   0.056700 0.040229 -0.0063000   
+DEAL::Points for cell with first point: 0.063000 0.021252 -0.063000
+DEAL::0.069300 0.023171 -0.056700   0.11970 0.043703 -0.056700   0.069300 0.040229 -0.056700   0.11970 0.055097 -0.056700   0.069300 0.023171 -0.0063000   0.11970 0.043703 -0.0063000   0.069300 0.040229 -0.0063000   0.11970 0.055097 -0.0063000   
+DEAL::Points for cell with first point: 0.0000 0.028000 0.0000
+DEAL::0.0063000 0.029459 0.0063000   0.056700 0.0018472 0.0063000   0.0063000 0.040854 0.0063000   0.056700 0.018906 0.0063000   0.0063000 0.029459 0.056700   0.056700 0.0018472 0.056700   0.0063000 0.040854 0.056700   0.056700 0.018906 0.056700   
+DEAL::Points for cell with first point: 0.063000 0.0000 0.0000
+DEAL::0.069300 0.0018472 0.0063000   0.11970 0.029459 0.0063000   0.069300 0.018906 0.0063000   0.11970 0.040854 0.0063000   0.069300 0.0018472 0.056700   0.11970 0.029459 0.056700   0.069300 0.018906 0.056700   0.11970 0.040854 0.056700   
+DEAL::Points for cell with first point: 0.0000 0.042252 0.0000
+DEAL::0.0063000 0.043703 0.0063000   0.056700 0.023171 0.0063000   0.0063000 0.055097 0.0063000   0.056700 0.040229 0.0063000   0.0063000 0.043703 0.056700   0.056700 0.023171 0.056700   0.0063000 0.055097 0.056700   0.056700 0.040229 0.056700   
+DEAL::Points for cell with first point: 0.063000 0.021252 0.0000
+DEAL::0.069300 0.023171 0.0063000   0.11970 0.043703 0.0063000   0.069300 0.040229 0.0063000   0.11970 0.055097 0.0063000   0.069300 0.023171 0.056700   0.11970 0.043703 0.056700   0.069300 0.040229 0.056700   0.11970 0.055097 0.056700   
+DEAL::
+DEAL::Mapping degree: 5
+DEAL::Points for cell with first point: 0.0000 0.028000 -0.063000
+DEAL::0.0063000 0.029303 -0.056700   0.056700 0.0020929 -0.056700   0.0063000 0.040730 -0.056700   0.056700 0.019101 -0.056700   0.0063000 0.029303 -0.0063000   0.056700 0.0020929 -0.0063000   0.0063000 0.040730 -0.0063000   0.056700 0.019101 -0.0063000   
+DEAL::Points for cell with first point: 0.063000 0.0000 -0.063000
+DEAL::0.069300 0.0020929 -0.056700   0.11970 0.029303 -0.056700   0.069300 0.019101 -0.056700   0.11970 0.040730 -0.056700   0.069300 0.0020929 -0.0063000   0.11970 0.029303 -0.0063000   0.069300 0.019101 -0.0063000   0.11970 0.040730 -0.0063000   
+DEAL::Points for cell with first point: 0.0000 0.042252 -0.063000
+DEAL::0.0063000 0.043586 -0.056700   0.056700 0.023353 -0.056700   0.0063000 0.055013 -0.056700   0.056700 0.040361 -0.056700   0.0063000 0.043586 -0.0063000   0.056700 0.023353 -0.0063000   0.0063000 0.055013 -0.0063000   0.056700 0.040361 -0.0063000   
+DEAL::Points for cell with first point: 0.063000 0.021252 -0.063000
+DEAL::0.069300 0.023353 -0.056700   0.11970 0.043586 -0.056700   0.069300 0.040361 -0.056700   0.11970 0.055013 -0.056700   0.069300 0.023353 -0.0063000   0.11970 0.043586 -0.0063000   0.069300 0.040361 -0.0063000   0.11970 0.055013 -0.0063000   
+DEAL::Points for cell with first point: 0.0000 0.028000 0.0000
+DEAL::0.0063000 0.029303 0.0063000   0.056700 0.0020929 0.0063000   0.0063000 0.040730 0.0063000   0.056700 0.019101 0.0063000   0.0063000 0.029303 0.056700   0.056700 0.0020929 0.056700   0.0063000 0.040730 0.056700   0.056700 0.019101 0.056700   
+DEAL::Points for cell with first point: 0.063000 0.0000 0.0000
+DEAL::0.069300 0.0020929 0.0063000   0.11970 0.029303 0.0063000   0.069300 0.019101 0.0063000   0.11970 0.040730 0.0063000   0.069300 0.0020929 0.056700   0.11970 0.029303 0.056700   0.069300 0.019101 0.056700   0.11970 0.040730 0.056700   
+DEAL::Points for cell with first point: 0.0000 0.042252 0.0000
+DEAL::0.0063000 0.043586 0.0063000   0.056700 0.023353 0.0063000   0.0063000 0.055013 0.0063000   0.056700 0.040361 0.0063000   0.0063000 0.043586 0.056700   0.056700 0.023353 0.056700   0.0063000 0.055013 0.056700   0.056700 0.040361 0.056700   
+DEAL::Points for cell with first point: 0.063000 0.021252 0.0000
+DEAL::0.069300 0.023353 0.0063000   0.11970 0.043586 0.0063000   0.069300 0.040361 0.0063000   0.11970 0.055013 0.0063000   0.069300 0.023353 0.056700   0.11970 0.043586 0.056700   0.069300 0.040361 0.056700   0.11970 0.055013 0.056700   
+DEAL::