From b7fc02857288dd20e44121e1a5d5aa5e26ccaf7f Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 12 Aug 2011 22:44:18 +0000 Subject: [PATCH] Make FEFaceValues and FESubFaceValues work for faces in 1d, i.e. on points. git-svn-id: https://svn.dealii.org/trunk@24055 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 8 + deal.II/include/deal.II/fe/fe_system.h | 18 -- .../deal.II/grid/tria_accessor.templates.h | 12 +- deal.II/source/base/geometry_info.cc | 3 +- deal.II/source/base/quadrature.cc | 16 +- deal.II/source/fe/fe_system.cc | 35 ---- deal.II/source/fe/mapping_q1.cc | 120 +++-------- tests/fe/fe_face_values_1d.cc | 189 ++++++++++++++++++ tests/fe/fe_face_values_1d/cmp/generic | 72 +++++++ 9 files changed, 316 insertions(+), 157 deletions(-) create mode 100644 tests/fe/fe_face_values_1d.cc create mode 100644 tests/fe/fe_face_values_1d/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 6ef4cab8f3..f2eddb807e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -91,6 +91,14 @@ changed to use std_cxx1x::_1, std_cxx1x::_2, etc from now on.

General

    +
  1. 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. +
    +(Wolfgang Bangerth, 2011/08/12) +
  2. 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 diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index e3e9ebc0b7..839d11996f 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -1267,24 +1267,6 @@ class FESystem : public FiniteElement }; }; -// 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 diff --git a/deal.II/include/deal.II/grid/tria_accessor.templates.h b/deal.II/include/deal.II/grid/tria_accessor.templates.h index f3dd876f1e..7ac0529a88 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -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::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::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); } diff --git a/deal.II/source/base/geometry_info.cc b/deal.II/source/base/geometry_info.cc index e6c9c3e4d6..3a567a9612 100644 --- a/deal.II/source/base/geometry_info.cc +++ b/deal.II/source/base/geometry_info.cc @@ -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; } diff --git a/deal.II/source/base/quadrature.cc b/deal.II/source/base/quadrature.cc index 1925c7e146..57a07f443a 100644 --- a/deal.II/source/base/quadrature.cc +++ b/deal.II/source/base/quadrature.cc @@ -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); } diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 3b1d0c996d..1f063c86a5 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -864,41 +864,6 @@ FESystem::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 void diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index ad814bd73b..799a5f4e14 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -911,73 +911,6 @@ MappingQ1::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 >&, - std::vector&, - std::vector >&, - std::vector >&) 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 >&, - std::vector&, - std::vector >&, - std::vector >&) 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 >&, - std::vector&, - std::vector >&, - std::vector >&) 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 >&, - std::vector&, - std::vector >&, - std::vector >&) const -{ - Assert(false, ExcNotImplemented()); -} @@ -985,24 +918,6 @@ namespace internal { namespace { - template - 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 &, - typename dealii::MappingQ1<1,spacedim>::InternalData &, - std::vector &, - std::vector > &, - std::vector > &) - { - Assert(false, ExcNotImplemented()); - } - - - template void compute_fill_face (const dealii::MappingQ1 &mapping, @@ -1053,12 +968,32 @@ namespace internal if (dim == spacedim) { for (unsigned int i=0; i::subface_ratio( diff --git a/tests/fe/fe_face_values_1d.cc b/tests/fe/fe_face_values_1d.cc new file mode 100644 index 0000000000..df73cbf35a --- /dev/null +++ b/tests/fe/fe_face_values_1d.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define PRECISION 4 + + + +template +inline void +plot_faces(Mapping &mapping, + FiniteElement &fe, + typename DoFHandler::cell_iterator &cell) +{ + // create a QGauss<0>(4), which should + // still only have 1 quadrature point + QGauss q(4); + Assert (q.size() == 1, ExcInternalError()); + + FEFaceValues 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::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 tria; + GridGenerator::hyper_cube (tria); + + FE_Q fe_q4(2); + DoFHandler dof(tria); + dof.distribute_dofs(fe_q4); + + typename DoFHandler::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 index 0000000000..47f9094c22 --- /dev/null +++ b/tests/fe/fe_face_values_1d/cmp/generic @@ -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:: -- 2.39.5