--- /dev/null
+//---------------------------- mapping_c1.h ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2000, 2001 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.
+//
+//---------------------------- mapping_c1.h ---------------------------
+#ifndef __deal2__mapping_c1_h
+#define __deal2__mapping_c1_h
+
+
+#include <fe/mapping_q.h>
+
+
+/**
+ * Mapping class that uses C1 (continuous) cubic mappings of the
+ * boundary. This class is built atop of @ref{MappingQ} by simply
+ * determining the interpolation points for a cubic mapping of the
+ * boundary differently: @ref{MappingQ} chooses them such that they
+ * interpolate the boundary, while this class chooses them such that
+ * the discretized boundary is globally continuous.
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+template <int dim>
+class MappingC1 : public MappingQ<dim>
+{
+ public:
+ /**
+ * Constructor. Pass the fixed
+ * degree @p{3} down to the base
+ * class, as a cubic mapping
+ * suffices to generate a
+ * continuous mapping of the
+ * boundary.
+ */
+ MappingC1 ();
+
+ protected:
+ /**
+ * For @p{dim=2,3}. Append the
+ * support points of all shape
+ * functions located on bounding
+ * lines to the vector
+ * @p{a}. Points located on the
+ * line but on vertices are not
+ * included.
+ *
+ * Needed by the
+ * @p{compute_support_points_simple(laplace)}
+ * functions. For @p{dim=1} this
+ * function is empty.
+ *
+ * This function chooses the
+ * respective points not such
+ * that they are interpolating
+ * the boundary (as does the base
+ * class), but rather such that
+ * the resulting cubic mapping is
+ * a continuous one.
+ */
+ virtual void
+ add_line_support_points (const typename Triangulation<dim>::cell_iterator &cell,
+ typename std::vector<Point<dim> > &a) const;
+
+ /**
+ * For @p{dim=3}. Append the
+ * support points of all shape
+ * functions located on bounding
+ * faces (quads in 3d) to the
+ * vector @p{a}. Points located
+ * on the line but on vertices
+ * are not included.
+ *
+ * Needed by the
+ * @p{compute_support_points_laplace}
+ * function. For @p{dim=1} and 2
+ * this function is empty.
+ *
+ * This function chooses the
+ * respective points not such
+ * that they are interpolating
+ * the boundary (as does the base
+ * class), but rather such that
+ * the resulting cubic mapping is
+ * a continuous one.
+ */
+ virtual void
+ add_quad_support_points(const typename Triangulation<dim>::cell_iterator &cell,
+ typename std::vector<Point<dim> > &a) const;
+};
+
+
+#endif
--- /dev/null
+//---------------------------- mapping_c1.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 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.
+//
+//---------------------------- mapping_c1.cc ---------------------------
+
+
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
+#include <fe/mapping_c1.h>
+#include <cmath>
+
+
+
+
+template <int dim>
+MappingC1<dim>::MappingC1 ()
+ :
+ MappingQ<dim> (3)
+{
+ Assert (dim > 1, ExcImpossibleInDim(dim));
+};
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+MappingC1<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
+ std::vector<Point<1> > &) const
+{
+ const unsigned int dim = 1;
+ Assert (dim > 1, ExcImpossibleInDim(dim));
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void
+MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &cell,
+ std::vector<Point<2> > &a) const
+{
+ const unsigned int dim = 2;
+ std::vector<Point<dim> > line_points (2);
+
+ // loop over each of the lines,
+ // and if it is at the
+ // boundary, then first get the
+ // boundary description and
+ // second compute the points on
+ // it. if not at the boundary,
+ // get the respective points
+ // from another function
+ for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+ {
+ const Triangulation<dim>::line_iterator line = cell->line(line_no);
+
+ if (line->at_boundary())
+ {
+ // first get the tangential
+ // vectors at the two
+ // vertices of this line
+ // from the boundary
+ // description
+ const Boundary<dim> &boundary
+ = line->get_triangulation().get_boundary(line->boundary_indicator());
+
+ Boundary<dim>::FaceVertexTangents face_vertex_tangents;
+ boundary.get_tangents_at_vertices (line, face_vertex_tangents);
+
+ // then transform them into
+ // interpolation points for
+ // a cubic polynomial
+ //
+ // for this, note that if
+ // we describe the boundary
+ // curve as a polynomial in
+ // tangential coordinate
+ // @p{t=0..1} (along the
+ // line) and @p{s} in
+ // normal direction, then
+ // the cubic mapping is
+ // such that @p{s = a*t**3
+ // + b*t**2 + c*t + d}, and
+ // we want to determine the
+ // interpolation points at
+ // @p{t=1/3} and
+ // @p{t=2/3}. Since at
+ // @p{t=0,1} we want a
+ // vertex which is actually
+ // at the boundary, we know
+ // that @p{d=0} and
+ // @p{a=-b-c}. As
+ // side-conditions, we want
+ // that the derivatives at
+ // @p{t=0} and @p{t=1},
+ // i.e. at the vertices
+ // match those returned by
+ // the boundary. We then
+ // have that
+ // @p{s(1/3)=1/27(2b+8c)}
+ // and
+ // @p{s(2/3)=4/27b+10/27c}.
+ //
+ // The task is then first
+ // to determine the
+ // coefficients from the
+ // tangentials. for that,
+ // first rotate the
+ // tangentials of @p{s(t)}
+ // into the global
+ // coordinate system. they
+ // are @p{A (1,c)} and @p{A
+ // (1,b)} with @p{A} the
+ // rotation matrix, since
+ // the tangentials in the
+ // coordinate system
+ // relative to the line are
+ // @p{(1,c)} and @p{(1,b)}
+ // at the two vertices,
+ // respectively. We then
+ // have to make sure by
+ // matching @p{b,c} that
+ // these tangentials are
+ // parallel to those
+ // returned by the boundary
+ // object
+ const Tensor<1,2> coordinate_vector = line->vertex(1) - line->vertex(0);
+ const double h = std::sqrt(coordinate_vector * coordinate_vector);
+ Tensor<1,2> coordinate_axis = coordinate_vector;
+ coordinate_axis /= h;
+
+ const double alpha = std::atan2(coordinate_axis[1], coordinate_axis[0]);
+ const double b = ((face_vertex_tangents.tangents[0][0][0] * sin(alpha)
+ -face_vertex_tangents.tangents[0][0][1] * cos(alpha)) /
+ (face_vertex_tangents.tangents[0][0][0] * cos(alpha)
+ +face_vertex_tangents.tangents[0][0][1] * sin(alpha))),
+ c = ((face_vertex_tangents.tangents[1][0][0] * sin(alpha)
+ -face_vertex_tangents.tangents[1][0][1] * cos(alpha)) /
+ (face_vertex_tangents.tangents[1][0][0] * cos(alpha)
+ +face_vertex_tangents.tangents[1][0][1] * sin(alpha)));
+
+ // next evaluate the so
+ // determined cubic
+ // polynomial at the points
+ // 1/3 and 2/3, first in
+ // unit coordinates
+ const Point<2> new_unit_points[2] = { Point<2>(1./3., 1./27.*(2*b+8*c)),
+ Point<2>(2./3., 4./27.*b+10./27.*c) };
+ // then transform these
+ // points to real
+ // coordinates by rotating,
+ // scaling and shifting
+ for (unsigned int i=0; i<2; ++i)
+ {
+ Point<2> real_point (cos(alpha) * new_unit_points[i][0]
+ - sin(alpha) * new_unit_points[i][1],
+ sin(alpha) * new_unit_points[i][0]
+ + cos(alpha) * new_unit_points[i][1]);
+ real_point *= h;
+ real_point += line->vertex(0);
+ a.push_back (real_point);
+ };
+ }
+ else
+ // not at boundary
+ {
+ static const StraightBoundary<dim> straight_boundary;
+ straight_boundary.get_intermediate_points_on_line (line, line_points);
+ a.insert (a.end(), line_points.begin(), line_points.end());
+ };
+ };
+};
+
+#endif
+
+
+
+template <int dim>
+void
+MappingC1<dim>::add_line_support_points (const typename Triangulation<dim>::cell_iterator &cell,
+ typename std::vector<Point<dim> > &a) const
+{
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+MappingC1<1>::add_quad_support_points (const Triangulation<1>::cell_iterator &,
+ std::vector<Point<1> > &) const
+{
+ const unsigned int dim = 1;
+ Assert (dim > 2, ExcImpossibleInDim(dim));
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+void
+MappingC1<2>::add_quad_support_points (const Triangulation<2>::cell_iterator &,
+ std::vector<Point<2> > &) const
+{
+ const unsigned int dim = 2;
+ Assert (dim > 2, ExcImpossibleInDim(dim));
+};
+
+#endif
+
+
+
+
+template <int dim>
+void
+MappingC1<dim>::add_quad_support_points (const typename Triangulation<dim>::cell_iterator &,
+ typename std::vector<Point<dim> > &) const
+{
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+
+// explicit instantiations
+template class MappingC1<deal_II_dimension>;
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ </p>
+ New: There is now a class <code class="class">MappingC1</code>
+ that implements a continous C<sup>1</sup> mapping of the
+ boundary by using a cubic mapping with continuous derivatives
+ at cell boundaries. The class presently only implements
+ something for 2d and 1d (where it does nothing).
+ <br>
+ (WB 2001/03/27)
+ </p>
+
<li> <p>
New: The static <code class="member">interpolate</code>, <code
class="member">create_right_hand_side</code>, <code
############################################################
+mapping_c1-cc-files = mapping_c1.cc
+
+ifeq ($(debug-mode),on)
+mapping_c1-o-files = $(mapping_c1-cc-files:.cc=.go)
+else
+mapping_c1-o-files = $(mapping_c1-cc-files:.cc=.o)
+endif
+
+mapping_c1.testcase: $(mapping_c1-o-files) $(libraries)
+ @echo =====linking======= $<
+ @$(CXX) $(LDFLAGS) -g -o $@ $^
+
+############################################################
+
internals-cc-files = internals.cc
ifeq ($(debug-mode),on)
--- /dev/null
+// $Id$
+// Copyright (C) 2001 Wolfgang Bangerth
+//
+// Test the continuity of normal vectors at vertices, and thus the
+// C1-ness of the C1-mapping
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <grid/tria.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/mapping_c1.h>
+#include <fe/fe_values.h>
+
+#include <algorithm>
+#include <fstream>
+
+int main ()
+{
+ std::ofstream logfile ("mapping_c1.output");
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ // create grid of circle, somehow
+ // arbitrarily from only one cell,
+ // but since we are not interested
+ // in the quality of the mesh, this
+ // is ok
+ Point<2> center;
+ HyperBallBoundary<2> circle(center, std::sqrt(2));
+ Triangulation<2> tria;
+ GridGenerator::hyper_cube(tria, -1, 1);
+ tria.set_boundary (0, circle);
+
+
+ // refine it more or less
+ // arbitrarily
+ tria.refine_global (1);
+ if (true)
+ {
+ Triangulation<2>::active_cell_iterator cell = tria.begin_active();
+ ++cell;
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement ();
+ };
+
+ // attach a dof handler to it
+ const FE_Q<2> fe(2);
+ DoFHandler<2> dof_handler (tria);
+ dof_handler.distribute_dofs (fe);
+
+ // and loop over all exterior faces
+ // to see whether the normal
+ // vectors are indeed continuous
+ // and pointing radially outward at
+ // the vertices
+ const QTrapez<1> quadrature;
+ const MappingC1<2> c1_mapping;
+ FEFaceValues<2> c1_values (c1_mapping, fe, quadrature,
+ update_q_points | update_normal_vectors);
+
+ // to compare with, also print the
+ // normal vectors as generated by a
+ // cubic mapping
+ const MappingQ<2> q3_mapping(3);
+ FEFaceValues<2> q3_values (q3_mapping, fe, quadrature,
+ update_q_points | update_normal_vectors);
+
+ for (DoFHandler<2>::active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f)
+ if (cell->face(f)->at_boundary())
+ {
+ c1_values.reinit (cell, f);
+ q3_values.reinit (cell, f);
+
+ // there should now be two
+ // normal vectors, one for
+ // each vertex of the face
+ Assert (c1_values.get_normal_vectors().size() == 2,
+ ExcInternalError());
+
+ // check that these two
+ // normal vectors have
+ // length approximately 1
+ // and point radially
+ // outward
+ for (unsigned int i=0; i<2; ++i)
+ {
+ Point<2> radius = c1_values.quadrature_point(i);
+ radius /= std::sqrt(radius.square());
+ deallog << "Normalized radius=" << radius << endl;
+
+ deallog << "C1 normal vector " << i << ": " << c1_values.normal_vector(i) << endl;
+ deallog << "Q3 normal vector " << i << ": " << q3_values.normal_vector(i) << endl;
+ };
+
+
+ // some numerical checks for correctness
+ for (unsigned int i=0; i<2; ++i)
+ {
+ Assert (std::fabs(c1_values.normal_vector(i) *
+ c1_values.normal_vector(i) - 1) < 1e-14,
+ ExcInternalError());
+ Point<2> radius = c1_values.quadrature_point(i);
+ radius /= std::sqrt(radius.square());
+
+ Assert ((radius-c1_values.normal_vector(i)).square() < 1e-14,
+ ExcInternalError());
+ };
+ };
+};
--- /dev/null
+JobId Tue Mar 27 15:27:48 2001
+DEAL::Normalized radius=-0.707107 -0.707107
+DEAL::C1 normal vector 0: -0.707107 -0.707107
+DEAL::Q3 normal vector 0: -0.710090 -0.704111
+DEAL::Normalized radius=0.000000 -1.000000
+DEAL::C1 normal vector 1: 0.000000 -1.000000
+DEAL::Q3 normal vector 1: 0.004228 -0.999991
+DEAL::Normalized radius=-0.707107 -0.707107
+DEAL::C1 normal vector 0: -0.707107 -0.707107
+DEAL::Q3 normal vector 0: -0.704111 -0.710090
+DEAL::Normalized radius=-1.000000 0.000000
+DEAL::C1 normal vector 1: -1.000000 0.000000
+DEAL::Q3 normal vector 1: -0.999991 0.004228
+DEAL::Normalized radius=1.000000 0.000000
+DEAL::C1 normal vector 0: 1.000000 0.000000
+DEAL::Q3 normal vector 0: 0.999991 -0.004228
+DEAL::Normalized radius=0.707107 0.707107
+DEAL::C1 normal vector 1: 0.707107 0.707107
+DEAL::Q3 normal vector 1: 0.704111 0.710090
+DEAL::Normalized radius=0.000000 1.000000
+DEAL::C1 normal vector 0: 0.000000 1.000000
+DEAL::Q3 normal vector 0: -0.004228 0.999991
+DEAL::Normalized radius=0.707107 0.707107
+DEAL::C1 normal vector 1: 0.707107 0.707107
+DEAL::Q3 normal vector 1: 0.710090 0.704111
+DEAL::Normalized radius=-0.707107 0.707107
+DEAL::C1 normal vector 0: -0.707107 0.707107
+DEAL::Q3 normal vector 0: -0.710090 0.704111
+DEAL::Normalized radius=0.000000 1.000000
+DEAL::C1 normal vector 1: 0.000000 1.000000
+DEAL::Q3 normal vector 1: 0.004228 0.999991
+DEAL::Normalized radius=-1.000000 0.000000
+DEAL::C1 normal vector 0: -1.000000 0.000000
+DEAL::Q3 normal vector 0: -0.999991 -0.004228
+DEAL::Normalized radius=-0.707107 0.707107
+DEAL::C1 normal vector 1: -0.707107 0.707107
+DEAL::Q3 normal vector 1: -0.704111 0.710090
+DEAL::Normalized radius=0.000000 -1.000000
+DEAL::C1 normal vector 0: -0.000000 -1.000000
+DEAL::Q3 normal vector 0: -0.000553 -1.000000
+DEAL::Normalized radius=0.382683 -0.923880
+DEAL::C1 normal vector 1: 0.382683 -0.923880
+DEAL::Q3 normal vector 1: 0.383194 -0.923668
+DEAL::Normalized radius=0.382683 -0.923880
+DEAL::C1 normal vector 0: 0.382683 -0.923880
+DEAL::Q3 normal vector 0: 0.382173 -0.924091
+DEAL::Normalized radius=0.707107 -0.707107
+DEAL::C1 normal vector 1: 0.707107 -0.707107
+DEAL::Q3 normal vector 1: 0.707497 -0.706716
+DEAL::Normalized radius=0.707107 -0.707107
+DEAL::C1 normal vector 0: 0.707107 -0.707107
+DEAL::Q3 normal vector 0: 0.706716 -0.707497
+DEAL::Normalized radius=0.923880 -0.382683
+DEAL::C1 normal vector 1: 0.923880 -0.382683
+DEAL::Q3 normal vector 1: 0.924091 -0.382173
+DEAL::Normalized radius=0.923880 -0.382683
+DEAL::C1 normal vector 0: 0.923880 -0.382683
+DEAL::Q3 normal vector 0: 0.923668 -0.383194
+DEAL::Normalized radius=1.000000 0.000000
+DEAL::C1 normal vector 1: 1.000000 0.000000
+DEAL::Q3 normal vector 1: 1.000000 0.000553