<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
};
};
-// 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
//---------------------------------------------------------------------------
// $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
internal::SubfaceCase<1>
CellAccessor<1>::subface_case(const unsigned int) const
{
- Assert(false, ExcImpossibleInDim(1));
return internal::SubfaceCase<1>::case_none;
}
internal::SubfaceCase<1>
CellAccessor<1,2>::subface_case(const unsigned int) const
{
- Assert(false, ExcImpossibleInDim(1));
return internal::SubfaceCase<1>::case_none;
}
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 <>
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);
}
// $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
GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
const unsigned int)
{
- Assert(false, ExcImpossibleInDim(1));
return 1;
}
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);
}
}
-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
-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());
-}
{
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,
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
{
{
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(
--- /dev/null
+//-------------------------------------------------------
+// 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>();
+}
+
--- /dev/null
+
+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::
projects / dealii-svn.git / commitdiff