From 5d509fd896f153535e5367ed9da3994bb855ce76 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 27 Mar 2001 13:31:35 +0000 Subject: [PATCH] Implement a C1 mapping. git-svn-id: https://svn.dealii.org/trunk@4300 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/mapping_c1.h | 99 ++++++++++ deal.II/deal.II/source/fe/mapping_c1.cc | 244 ++++++++++++++++++++++++ deal.II/doc/news/2001/c-3-1.html | 11 ++ tests/fe/Makefile.in | 14 ++ tests/fe/mapping_c1.cc | 115 +++++++++++ tests/fe/mapping_c1.checked | 61 ++++++ 6 files changed, 544 insertions(+) create mode 100644 deal.II/deal.II/include/fe/mapping_c1.h create mode 100644 deal.II/deal.II/source/fe/mapping_c1.cc create mode 100644 tests/fe/mapping_c1.cc create mode 100644 tests/fe/mapping_c1.checked 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 index 0000000000..36d2d05358 --- /dev/null +++ b/deal.II/deal.II/include/fe/mapping_c1.h @@ -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 + + +/** + * 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 +class MappingC1 : public MappingQ +{ + 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::cell_iterator &cell, + typename std::vector > &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::cell_iterator &cell, + typename std::vector > &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 index 0000000000..c3cf564850 --- /dev/null +++ b/deal.II/deal.II/source/fe/mapping_c1.cc @@ -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 +#include +#include +#include +#include + + + + +template +MappingC1::MappingC1 () + : + MappingQ (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 > &) 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 > &a) const +{ + const unsigned int dim = 2; + std::vector > 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::lines_per_cell; ++line_no) + { + const Triangulation::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 &boundary + = line->get_triangulation().get_boundary(line->boundary_indicator()); + + Boundary::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 straight_boundary; + straight_boundary.get_intermediate_points_on_line (line, line_points); + a.insert (a.end(), line_points.begin(), line_points.end()); + }; + }; +}; + +#endif + + + +template +void +MappingC1::add_line_support_points (const typename Triangulation::cell_iterator &cell, + typename std::vector > &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 > &) 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 > &) const +{ + const unsigned int dim = 2; + Assert (dim > 2, ExcImpossibleInDim(dim)); +}; + +#endif + + + + +template +void +MappingC1::add_quad_support_points (const typename Triangulation::cell_iterator &, + typename std::vector > &) const +{ + Assert (false, ExcNotImplemented()); +}; + + + + +// explicit instantiations +template class MappingC1; diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html index 290a397a39..5573d92097 100644 --- a/deal.II/doc/news/2001/c-3-1.html +++ b/deal.II/doc/news/2001/c-3-1.html @@ -280,6 +280,17 @@ documentation, etc.

deal.II

    +
  1. +

    + New: There is now a class MappingC1 + that implements a continous C1 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). +
    + (WB 2001/03/27) +

    +
  2. New: The static interpolate, create_right_hand_side, +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +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::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 index 0000000000..f0a193c2cc --- /dev/null +++ b/tests/fe/mapping_c1.checked @@ -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 -- 2.39.5