]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a C1 mapping.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 27 Mar 2001 13:31:35 +0000 (13:31 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 27 Mar 2001 13:31:35 +0000 (13:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@4300 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/mapping_c1.h [new file with mode: 0644]
deal.II/deal.II/source/fe/mapping_c1.cc [new file with mode: 0644]
deal.II/doc/news/2001/c-3-1.html
tests/fe/Makefile.in
tests/fe/mapping_c1.cc [new file with mode: 0644]
tests/fe/mapping_c1.checked [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/fe/mapping_c1.h b/deal.II/deal.II/include/fe/mapping_c1.h
new file mode 100644 (file)
index 0000000..36d2d05
--- /dev/null
@@ -0,0 +1,99 @@
+//----------------------------  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
diff --git a/deal.II/deal.II/source/fe/mapping_c1.cc b/deal.II/deal.II/source/fe/mapping_c1.cc
new file mode 100644 (file)
index 0000000..c3cf564
--- /dev/null
@@ -0,0 +1,244 @@
+//----------------------------  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>;
index 290a397a391e04d0c1868e894ad44f05e866b03c..5573d920975d3c14c682507701ac4525b76a1fd1 100644 (file)
@@ -280,6 +280,17 @@ documentation, etc</a>.
 <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
index 4eeb6576bdb2ee085004858c939a23f182eb8f1c..216def6d13f874cde3c7fee2d67f74586463bb7d 100644 (file)
@@ -164,6 +164,20 @@ mapping.testcase: $(mapping-o-files) $(libraries)
 
 ############################################################
 
+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)
diff --git a/tests/fe/mapping_c1.cc b/tests/fe/mapping_c1.cc
new file mode 100644 (file)
index 0000000..ed7e726
--- /dev/null
@@ -0,0 +1,115 @@
+// $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());
+           };
+       };  
+};
diff --git a/tests/fe/mapping_c1.checked b/tests/fe/mapping_c1.checked
new file mode 100644 (file)
index 0000000..f0a193c
--- /dev/null
@@ -0,0 +1,61 @@
+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

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.