]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make FEFaceValues and FESubFaceValues work for faces in 1d, i.e. on points.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 12 Aug 2011 22:44:18 +0000 (22:44 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 12 Aug 2011 22:44:18 +0000 (22:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@24055 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/fe_system.h
deal.II/include/deal.II/grid/tria_accessor.templates.h
deal.II/source/base/geometry_info.cc
deal.II/source/base/quadrature.cc
deal.II/source/fe/fe_system.cc
deal.II/source/fe/mapping_q1.cc
tests/fe/fe_face_values_1d.cc [new file with mode: 0644]
tests/fe/fe_face_values_1d/cmp/generic [new file with mode: 0644]

index 6ef4cab8f353b255c2982e4ee48b46ce7fb0e926..f2eddb807efac5ad3ce4807f7dd71218cd08fc6a 100644 (file)
@@ -91,6 +91,14 @@ changed to use std_cxx1x::_1, std_cxx1x::_2, etc from now on.
 <h3>General</h3>
 
 <ol>
+<li> Extended: Many operations on objects of type Point<0>, Quadrature<0>, etc
+(including creation) were previously forbidden since such objects do not make
+much sense. However, this prevented a lot of code that could otherwise work
+in a dimension independent way, from working in 1d, e.g. integration on
+faces. Many of these places have now been cleaned up and work.
+<br>
+(Wolfgang Bangerth, 2011/08/12)
+
 <li> Extended: The classes Tensor, SymmetricTensor and Point now have an
 additional template argument for the number type. It is now possible to base
 these classes on any abstract data type that implements basic arithmetic
index e3e9ebc0b7c39d844721061ddeb9d0bdb471cee0..839d11996f1af47ff0980c2c5da8e09b7feff351 100644 (file)
@@ -1267,24 +1267,6 @@ class FESystem : public FiniteElement<dim,spacedim>
     };
 };
 
-// Force explicit implementation of the following templates
-template <>
-Mapping<1,1>::InternalDataBase*
-FESystem<1,1>::get_face_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const;
-
-template <>
-Mapping<1,1>::InternalDataBase*
-FESystem<1,1>::get_subface_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const;
-
-template <>
-Mapping<1,2>::InternalDataBase*
-FESystem<1,2>::get_face_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const;
-
-template <>
-Mapping<1,2>::InternalDataBase*
-FESystem<1,2>::get_subface_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const;
-
-
 
 DEAL_II_NAMESPACE_CLOSE
 
index f3dd876f1e8c167c5d5fed77a55ad796954d46ee..7ac0529a8852497b490c428138a7b316160a42e5 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -2600,7 +2600,6 @@ inline
 internal::SubfaceCase<1>
 CellAccessor<1>::subface_case(const unsigned int) const
 {
-  Assert(false, ExcImpossibleInDim(1));
   return internal::SubfaceCase<1>::case_none;
 }
 
@@ -2609,7 +2608,6 @@ inline
 internal::SubfaceCase<1>
 CellAccessor<1,2>::subface_case(const unsigned int) const
 {
-  Assert(false, ExcImpossibleInDim(1));
   return internal::SubfaceCase<1>::case_none;
 }
 
@@ -2622,7 +2620,9 @@ CellAccessor<2>::subface_case(const unsigned int face_no) const
   Assert(active(), TriaAccessorExceptions::ExcCellNotActive());
   Assert(face_no<GeometryInfo<2>::faces_per_cell,
         ExcIndexRange(face_no,0,GeometryInfo<2>::faces_per_cell));
-  return (face(face_no)->has_children()) ? internal::SubfaceCase<2>::case_x : internal::SubfaceCase<2>::case_none;
+  return ((face(face_no)->has_children()) ?
+         internal::SubfaceCase<2>::case_x :
+         internal::SubfaceCase<2>::case_none);
 }
 
 template <>
@@ -2633,7 +2633,9 @@ CellAccessor<2,3>::subface_case(const unsigned int face_no) const
   Assert(active(), TriaAccessorExceptions::ExcCellNotActive());
   Assert(face_no<GeometryInfo<2>::faces_per_cell,
         ExcIndexRange(face_no,0,GeometryInfo<2>::faces_per_cell));
-  return (face(face_no)->has_children()) ? internal::SubfaceCase<2>::case_x : internal::SubfaceCase<2>::case_none;
+  return ((face(face_no)->has_children()) ?
+         internal::SubfaceCase<2>::case_x :
+         internal::SubfaceCase<2>::case_none);
 }
 
 
index e6c9c3e4d63f3d783b3fd9641b854f9060104e57..3a567a9612654b3b1f33ba695b7d0193444e1d4a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -280,7 +280,6 @@ double
 GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
                               const unsigned int)
 {
-  Assert(false, ExcImpossibleInDim(1));
   return 1;
 }
 
index 1925c7e146d1991963cc050d751f65adab643759..57a07f443a80bee353833ed20bc8a573bab0964f 100644 (file)
@@ -1148,16 +1148,22 @@ face (const unsigned int face_no,
 template <>
 QProjector<1>::DataSetDescriptor
 QProjector<1>::DataSetDescriptor::
-subface (const unsigned int,
-         const unsigned int,
+subface (const unsigned int face_no,
+         const unsigned int subface_no,
          const bool,
          const bool,
          const bool,
-         const unsigned int,
+         const unsigned int n_quadrature_points,
         const internal::SubfaceCase<1>)
 {
-  Assert (false, ExcInternalError());
-  return deal_II_numbers::invalid_unsigned_int;
+  Assert (face_no < GeometryInfo<1>::faces_per_cell,
+          ExcInternalError());
+  Assert (subface_no < GeometryInfo<1>::max_children_per_face,
+          ExcInternalError());
+
+  return ((face_no * GeometryInfo<1>::max_children_per_face +
+          subface_no)
+         * n_quadrature_points);
 }
 
 
index 3b1d0c996d6f2d15ccd04cbc65cd485f0a4efcf5..1f063c86a5293454febb2f0f6876778440d1ec51 100644 (file)
@@ -864,41 +864,6 @@ FESystem<dim,spacedim>::get_subface_data (
 }
 
 
-template <>
-Mapping<1,1>::InternalDataBase *
-FESystem<1,1>::get_face_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const
-{
-  Assert (false, ExcNotImplemented ());
-  return 0;
-}
-
-
-template <>
-Mapping<1,1>::InternalDataBase *
-FESystem<1,1>::get_subface_data (const UpdateFlags, const Mapping<1,1>&, const Quadrature<0> &) const
-{
-  Assert (false, ExcNotImplemented ());
-  return 0;
-}
-
-
-template <>
-Mapping<1,2>::InternalDataBase *
-FESystem<1,2>::get_face_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const
-{
-  Assert (false, ExcNotImplemented ());
-  return 0;
-}
-
-
-template <>
-Mapping<1,2>::InternalDataBase *
-FESystem<1,2>::get_subface_data (const UpdateFlags, const Mapping<1,2>&, const Quadrature<0> &) const
-{
-  Assert (false, ExcNotImplemented ());
-  return 0;
-}
-
 
 template <int dim, int spacedim>
 void
index ad814bd73be1124418ffe54ab29aa81e8e3c0e01..799a5f4e14ee08c6b19fc4de5554f90046890f35 100644 (file)
@@ -911,73 +911,6 @@ MappingQ1<dim,spacedim>::fill_fe_values (
 
 
 
-template <>
-void
-MappingQ1<1,1>::fill_fe_face_values (
-  const Triangulation<1,1>::cell_iterator &,
-  const unsigned,
-  const Quadrature<0>&,
-  Mapping<1,1>::InternalDataBase&,
-  std::vector<Point<1> >&,
-  std::vector<double>&,
-  std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <>
-void
-MappingQ1<1,2>::fill_fe_face_values (
-  const Triangulation<1,2>::cell_iterator &,
-  const unsigned,
-  const Quadrature<0>&,
-  Mapping<1,2>::InternalDataBase&,
-  std::vector<Point<2> >&,
-  std::vector<double>&,
-  std::vector<Tensor<1,2> >&,
-  std::vector<Point<2> >&) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <>
-void
-MappingQ1<1,1>::fill_fe_subface_values (
-  const Triangulation<1,1>::cell_iterator &,
-  const unsigned,
-  const unsigned,
-  const Quadrature<0>&,
-  Mapping<1,1>::InternalDataBase&,
-  std::vector<Point<1> >&,
-  std::vector<double>&,
-  std::vector<Tensor<1,1> >&,
-  std::vector<Point<1> >&) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <>
-void
-MappingQ1<1,2>::fill_fe_subface_values (
-  const Triangulation<1,2>::cell_iterator &,
-  const unsigned,
-  const unsigned,
-  const Quadrature<0>&,
-  Mapping<1,2>::InternalDataBase&,
-  std::vector<Point<2> >&,
-  std::vector<double>&,
-  std::vector<Tensor<1,2> >&,
-  std::vector<Point<2> >&) const
-{
-  Assert(false, ExcNotImplemented());
-}
 
 
 
@@ -985,24 +918,6 @@ namespace internal
 {
   namespace
   {
-    template <int spacedim>
-    void
-    compute_fill_face (const dealii::MappingQ1<1,spacedim> &,
-                      const typename dealii::Triangulation<1,spacedim>::cell_iterator &,
-                      const unsigned int,
-                      const unsigned int,
-                      const unsigned int,
-                      const std::vector<double> &,
-                      typename dealii::MappingQ1<1,spacedim>::InternalData &,
-                      std::vector<double> &,
-                      std::vector<Tensor<1,spacedim> > &,
-                      std::vector<Point<spacedim> > &)
-    {
-      Assert(false, ExcNotImplemented());
-    }
-
-
-
     template <int dim, int spacedim>
     void
     compute_fill_face (const dealii::MappingQ1<dim,spacedim> &mapping,
@@ -1053,12 +968,32 @@ namespace internal
          if (dim == spacedim)
            {
              for (unsigned int i=0; i<n_q_points; ++i)
-               if (dim == 2)
-                 cross_product (boundary_forms[i], data.aux[0][i]);
-               else if (dim == 3)
-                 cross_product (boundary_forms[i], data.aux[0][i], data.aux[1][i]);
-               else
-                 Assert(false, ExcNotImplemented());
+               switch (dim)
+                 {
+                   case 1:
+                                                          // in 1d, we don't
+                                                          // have access to
+                                                          // any of the
+                                                          // data.aux fields,
+                                                          // but we can still
+                                                          // compute the
+                                                          // boundary form
+                                                          // simply by simply
+                                                          // looking at the
+                                                          // number of the
+                                                          // face
+                         boundary_forms[i][0] = (face_no == 0 ?
+                                                 -1 : +1);
+                         break;
+                   case 2:
+                         cross_product (boundary_forms[i], data.aux[0][i]);
+                         break;
+                   case 3:
+                         cross_product (boundary_forms[i], data.aux[0][i], data.aux[1][i]);
+                         break;
+                   default:
+                         Assert(false, ExcNotImplemented());
+                 }
            }
          else
            {
@@ -1140,7 +1075,8 @@ namespace internal
              {
                if (update_flags & update_JxW_values)
                  {
-                   JxW_values[i] = boundary_forms[i].norm() * weights[i];
+                   JxW_values[i] = boundary_forms[i].norm() * weights[i];
+
                    if (subface_no!=deal_II_numbers::invalid_unsigned_int)
                      {
                        const double area_ratio=GeometryInfo<dim>::subface_ratio(
diff --git a/tests/fe/fe_face_values_1d.cc b/tests/fe/fe_face_values_1d.cc
new file mode 100644 (file)
index 0000000..df73cbf
--- /dev/null
@@ -0,0 +1,189 @@
+//-------------------------------------------------------
+//    rt_11.cc,v 1.3 2003/06/09 16:00:38 wolf Exp
+//    Version: 
+//
+//    Copyright (C) 2011 by the deal.II authors
+//
+//    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.
+//
+//-------------------------------------------------------
+
+// A variant of the mapping.cc test: shows the shape functions on the faces of
+// 1d cells
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/mapping_cartesian.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe.h>
+#include <vector>
+#include <fstream>
+#include <iomanip>
+#include <string>
+#include <sstream>
+
+#define PRECISION 4
+
+
+
+template<int dim>
+inline void
+plot_faces(Mapping<dim> &mapping,
+          FiniteElement<dim> &fe,
+          typename DoFHandler<dim>::cell_iterator &cell)
+{
+                                  // create a QGauss<0>(4), which should
+                                  // still only have 1 quadrature point
+  QGauss<dim-1> q(4);
+  Assert (q.size() == 1, ExcInternalError());
+
+  FEFaceValues<dim> fe_values(mapping, fe, q,
+                             UpdateFlags(update_q_points
+                                         | update_JxW_values
+                                         | update_values
+                                         | update_gradients
+                                         | update_hessians
+                                         | update_normal_vectors));
+
+  for (unsigned int face_nr=0;
+       face_nr < GeometryInfo<dim>::faces_per_cell;
+       ++ face_nr)
+    {
+      deallog << "Face=" << face_nr << std::endl;
+      fe_values.reinit(cell, face_nr);
+
+                                      // print some data on the location
+      deallog << "x=" << fe_values.quadrature_point(0) << std::endl;
+      deallog << "n=" << fe_values.normal_vector(0) << std::endl;
+      deallog << "JxW=" << fe_values.JxW(0) << std::endl;
+
+                                      // now print some data on the shape
+                                      // functions
+      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+       deallog << "shape_function " << i << ":" << std::endl
+               << "  phi=" << fe_values.shape_value(i,0) << std::endl
+               << "  grad phi=" << fe_values.shape_grad(i,0) << std::endl
+               << "  grad^2 phi=" << fe_values.shape_hessian(i,0) << std::endl;
+
+      deallog << std::endl;
+    }
+}
+
+
+
+template<int dim>
+inline void
+plot_subfaces(Mapping<dim> &mapping,
+             FiniteElement<dim> &fe,
+             typename DoFHandler<dim>::cell_iterator &cell)
+{
+
+                                  // create a QGauss<0>(4), which should
+                                  // still only have 1 quadrature point
+  QGauss<dim-1> q(4);
+  Assert (q.size() == 1, ExcInternalError());
+
+  FESubfaceValues<dim> fe_values(mapping, fe, q,
+                                UpdateFlags(update_q_points
+                                            | update_JxW_values
+                                            | update_values
+                                            | update_gradients
+                                            | update_hessians
+                                            | update_normal_vectors));
+  for (unsigned int face_nr=0;
+       face_nr < GeometryInfo<dim>::faces_per_cell;
+       ++ face_nr)
+    for (unsigned int sub_nr=0;
+        sub_nr < GeometryInfo<dim>::max_children_per_face;
+        ++ sub_nr)
+      {
+       deallog << "Face=" << face_nr
+               << ", subface=" << sub_nr
+               << std::endl;
+       
+       fe_values.reinit(cell, face_nr, sub_nr);
+       
+                                        // print some data on the location
+       deallog << "x=" << fe_values.quadrature_point(0) << std::endl;
+       deallog << "n=" << fe_values.normal_vector(0) << std::endl;
+       deallog << "JxW=" << fe_values.JxW(0) << std::endl;
+
+                                        // now print some data on the shape
+                                        // functions
+       for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+         deallog << "shape_function " << i << ":" << std::endl
+                 << "  phi=" << fe_values.shape_value(i,0) << std::endl
+                 << "  grad phi=" << fe_values.shape_grad(i,0) << std::endl
+                 << "  grad^2 phi=" << fe_values.shape_hessian(i,0) << std::endl;
+
+       deallog << std::endl;
+      }
+}
+
+  
+template<int dim>
+void mapping_test()
+{
+  deallog << "dim=" << dim << std::endl;
+  
+  std::vector<Mapping<dim> *> mapping_ptr;
+  std::vector<std::string> mapping_strings;
+
+  MappingQ1<dim> mapping;
+  std::string mapping_name = "MappingQ1";
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+
+  FE_Q<dim> fe_q4(2);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe_q4);      
+
+  typename DoFHandler<dim>::cell_iterator cell = dof.begin_active();
+  {
+    std::ostringstream ost;
+    ost << "MappingFace" << dim << "d-"
+       << mapping_name;
+    deallog << ost.str() << std::endl;     
+    plot_faces(mapping, fe_q4, cell);
+  }
+
+  {
+    std::ostringstream ost;
+    ost << "MappingSubface" << dim << "d-"
+       << mapping_name;
+    deallog << ost.str() << std::endl;     
+    plot_subfaces(mapping, fe_q4, cell);
+  }
+}
+
+
+
+int main()
+{
+  std::ofstream logfile ("fe_face_values_1d/output");
+  deallog << std::setprecision(PRECISION);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+                                  // ----------------------- 
+                                  // Tests for dim=1
+                                  // -----------------------
+  mapping_test<1>();
+}
+
diff --git a/tests/fe/fe_face_values_1d/cmp/generic b/tests/fe/fe_face_values_1d/cmp/generic
new file mode 100644 (file)
index 0000000..47f9094
--- /dev/null
@@ -0,0 +1,72 @@
+
+DEAL::dim=1
+DEAL::MappingFace1d-MappingQ1
+DEAL::Face=0
+DEAL::x=0.000
+DEAL::n=-1.000
+DEAL::JxW=1.000
+DEAL::shape_function 0:
+DEAL::  phi=1.000
+DEAL::  grad phi=-3.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 1:
+DEAL::  phi=0
+DEAL::  grad phi=-1.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 2:
+DEAL::  phi=0
+DEAL::  grad phi=4.000
+DEAL::  grad^2 phi=-8.000
+DEAL::
+DEAL::Face=1
+DEAL::x=1.000
+DEAL::n=1.000
+DEAL::JxW=1.000
+DEAL::shape_function 0:
+DEAL::  phi=0
+DEAL::  grad phi=1.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 1:
+DEAL::  phi=1.000
+DEAL::  grad phi=3.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 2:
+DEAL::  phi=0
+DEAL::  grad phi=-4.000
+DEAL::  grad^2 phi=-8.000
+DEAL::
+DEAL::MappingSubface1d-MappingQ1
+DEAL::Face=0, subface=0
+DEAL::x=0.000
+DEAL::n=-1.000
+DEAL::JxW=1.000
+DEAL::shape_function 0:
+DEAL::  phi=1.000
+DEAL::  grad phi=-3.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 1:
+DEAL::  phi=0
+DEAL::  grad phi=-1.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 2:
+DEAL::  phi=0
+DEAL::  grad phi=4.000
+DEAL::  grad^2 phi=-8.000
+DEAL::
+DEAL::Face=1, subface=0
+DEAL::x=1.000
+DEAL::n=1.000
+DEAL::JxW=1.000
+DEAL::shape_function 0:
+DEAL::  phi=0
+DEAL::  grad phi=1.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 1:
+DEAL::  phi=1.000
+DEAL::  grad phi=3.000
+DEAL::  grad^2 phi=4.000
+DEAL::shape_function 2:
+DEAL::  phi=0
+DEAL::  grad phi=-4.000
+DEAL::  grad^2 phi=-8.000
+DEAL::

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.