From: Luca Heltai Date: Mon, 29 Feb 2016 15:51:53 +0000 (+0100) Subject: Copy of mapping generic. X-Git-Tag: v8.5.0-rc1~1131^2~37 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b1de4eef72359bf6f2969cd1adca1482624a5ed3;p=dealii.git Copy of mapping generic. --- diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h new file mode 100644 index 0000000000..2b072b9615 --- /dev/null +++ b/include/deal.II/fe/mapping_manifold.h @@ -0,0 +1,782 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii__mapping_manifold_h +#define dealii__mapping_manifold_h + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +template class MappingQ; + + +/*!@addtogroup mapping */ +/*@{*/ + + +/** + * This class implements the functionality for polynomial mappings $Q_p$ of + * polynomial degree $p$ that will be used on all cells of the mesh. The + * MappingQ1 and MappingQ classes specialize this behavior slightly. + * + * The class is poorly named. It should really have been called MappingQ + * because it consistently uses $Q_p$ mappings on all cells of a + * triangulation. However, the name MappingQ was already taken when we rewrote + * the entire class hierarchy for mappings. One might argue that one should + * always use MappingQGeneric over the existing class MappingQ (which, unless + * explicitly specified during the construction of the object, only uses + * mappings of degree $p$ on cells at the boundary of the domain). On + * the other hand, there are good reasons to use MappingQ in many situations: + * in many situations, curved domains are only provided with information about + * how exactly edges at the boundary are shaped, but we do not know anything + * about internal edges. Thus, in the absence of other information, we can + * only assume that internal edges are straight lines, and in that case + * internal cells may as well be treated is bilinear quadrilaterals or + * trilinear hexahedra. (An example of how such meshes look is shown in step-1 + * already, but it is also discussed in the "Results" section of step-6.) + * Because bi-/trilinear mappings are significantly cheaper to compute than + * higher order mappings, it is advantageous in such situations to use the + * higher order mapping only on cells at the boundary of the domain -- i.e., + * the behavior of MappingQ. Of course, MappingQGeneric also uses bilinear + * mappings for interior cells as long as it has no knowledge about curvature + * of interior edges, but it implements this the expensive way: as a general + * $Q_p$ mapping where the mapping support points just happen to be + * arranged along linear or bilinear edges or faces. + * + * There are a number of special cases worth considering: + * - If you really want to use a higher order mapping for all cells, + * you can do this using the current class, but this only makes sense if you + * can actually provide information about how interior edges and faces of the + * mesh should be curved. This is typically done by associating a Manifold + * with interior cells and edges. A simple example of this is discussed in the + * "Results" section of step-6; a full discussion of manifolds is provided in + * step-53. + * - If you are working on meshes that describe a (curved) manifold + * embedded in higher space dimensions, i.e., if dim!=spacedim, then every + * cell is at the boundary of the domain you will likely already have attached + * a manifold object to all cells that can then also be used by the mapping + * classes for higher order mappings. + * + * + * @author Wolfgang Bangerth, 2015 + */ +template +class MappingManifold : public Mapping +{ +public: + /** + * Constructor. @p polynomial_degree denotes the polynomial degree of the + * polynomials that are used to map cells from the reference to the real + * cell. + */ + MappingManifold (const unsigned int polynomial_degree); + + /** + * Copy constructor. + */ + MappingManifold (const MappingManifold &mapping); + + // for documentation, see the Mapping base class + virtual + Mapping *clone () const; + + /** + * Return the degree of the mapping, i.e. the value which was passed to the + * constructor. + */ + unsigned int get_degree () const; + + /** + * Always returns @p true because the default implementation of functions in + * this class preserves vertex locations. + */ + virtual + bool preserves_vertex_locations () const; + + /** + * @name Mapping points between reference and real cells + * @{ + */ + + // for documentation, see the Mapping base class + virtual + Point + transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const; + + // for documentation, see the Mapping base class + virtual + Point + transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const; + + /** + * @} + */ + + /** + * @name Functions to transform tensors from reference to real coordinates + * @{ + */ + + // for documentation, see the Mapping base class + virtual + void + transform (const ArrayView > &input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + const ArrayView > &output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const ArrayView > &input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + const ArrayView > &output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const ArrayView > &input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + const ArrayView > &output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const ArrayView > &input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + const ArrayView > &output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const ArrayView > &input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + const ArrayView > &output) const; + + /** + * @} + */ + + /** + * @name Interface with FEValues + * @{ + */ + +public: + /** + * Storage for internal data of polynomial mappings. See + * Mapping::InternalDataBase for an extensive description. + * + * For the current class, the InternalData class stores data that is + * computed once when the object is created (in get_data()) as well as data + * the class wants to store from between the call to fill_fe_values(), + * fill_fe_face_values(), or fill_fe_subface_values() until possible later + * calls from the finite element to functions such as transform(). The + * latter class of member variables are marked as 'mutable'. + */ + class InternalData : public Mapping::InternalDataBase + { + public: + /** + * Constructor. The argument denotes the polynomial degree of the mapping + * to which this object will correspond. + */ + InternalData(const unsigned int polynomial_degree); + + /** + * Initialize the object's member variables related to cell data based on + * the given arguments. + * + * The function also calls compute_shape_function_values() to actually set + * the member variables related to the values and derivatives of the + * mapping shape functions. + */ + void + initialize (const UpdateFlags update_flags, + const Quadrature &quadrature, + const unsigned int n_original_q_points); + + /** + * Initialize the object's member variables related to cell and face data + * based on the given arguments. In order to initialize cell data, this + * function calls initialize(). + */ + void + initialize_face (const UpdateFlags update_flags, + const Quadrature &quadrature, + const unsigned int n_original_q_points); + + /** + * Compute the values and/or derivatives of the shape functions used for + * the mapping. + * + * Which values, derivatives, or higher order derivatives are computed is + * determined by which of the member arrays have nonzero sizes. They are + * typically set to their appropriate sizes by the initialize() and + * initialize_face() functions, which indeed call this function + * internally. However, it is possible (and at times useful) to do the + * resizing by hand and then call this function directly. An example is in + * a Newton iteration where we update the location of a quadrature point + * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to re- + * compute the mapping and its derivatives at this location, but have + * already sized all internal arrays correctly. + */ + void compute_shape_function_values (const std::vector > &unit_points); + + + /** + * Shape function at quadrature point. Shape functions are in tensor + * product order, so vertices must be reordered to obtain transformation. + */ + const double &shape (const unsigned int qpoint, + const unsigned int shape_nr) const; + + /** + * Shape function at quadrature point. See above. + */ + double &shape (const unsigned int qpoint, + const unsigned int shape_nr); + + /** + * Gradient of shape function in quadrature point. See above. + */ + const Tensor<1,dim> &derivative (const unsigned int qpoint, + const unsigned int shape_nr) const; + + /** + * Gradient of shape function in quadrature point. See above. + */ + Tensor<1,dim> &derivative (const unsigned int qpoint, + const unsigned int shape_nr); + + /** + * Second derivative of shape function in quadrature point. See above. + */ + const Tensor<2,dim> &second_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const; + + /** + * Second derivative of shape function in quadrature point. See above. + */ + Tensor<2,dim> &second_derivative (const unsigned int qpoint, + const unsigned int shape_nr); + + /** + * third derivative of shape function in quadrature point. See above. + */ + const Tensor<3,dim> &third_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const; + + /** + * third derivative of shape function in quadrature point. See above. + */ + Tensor<3,dim> &third_derivative (const unsigned int qpoint, + const unsigned int shape_nr); + + /** + * fourth derivative of shape function in quadrature point. See above. + */ + const Tensor<4,dim> &fourth_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const; + + /** + * fourth derivative of shape function in quadrature point. See above. + */ + Tensor<4,dim> &fourth_derivative (const unsigned int qpoint, + const unsigned int shape_nr); + + /** + * Return an estimate (in bytes) or the memory consumption of this object. + */ + virtual std::size_t memory_consumption () const; + + /** + * Values of shape functions. Access by function @p shape. + * + * Computed once. + */ + std::vector shape_values; + + /** + * Values of shape function derivatives. Access by function @p derivative. + * + * Computed once. + */ + std::vector > shape_derivatives; + + /** + * Values of shape function second derivatives. Access by function @p + * second_derivative. + * + * Computed once. + */ + std::vector > shape_second_derivatives; + + /** + * Values of shape function third derivatives. Access by function @p + * second_derivative. + * + * Computed once. + */ + std::vector > shape_third_derivatives; + + /** + * Values of shape function fourth derivatives. Access by function @p + * second_derivative. + * + * Computed once. + */ + std::vector > shape_fourth_derivatives; + + /** + * Unit tangential vectors. Used for the computation of boundary forms and + * normal vectors. + * + * This vector has (dim-1)GeometryInfo::faces_per_cell entries. The first + * GeometryInfo::faces_per_cell contain the vectors in the first + * tangential direction for each face; the second set of + * GeometryInfo::faces_per_cell entries contain the vectors in the second + * tangential direction (only in 3d, since there we have 2 tangential + * directions per face), etc. + * + * Filled once. + */ + std::vector > > unit_tangentials; + + /** + * The polynomial degree of the mapping. Since the objects here are also + * used (with minor adjustments) by MappingQ, we need to store this. + */ + unsigned int polynomial_degree; + + /** + * Number of shape functions. If this is a Q1 mapping, then it is simply + * the number of vertices per cell. However, since also derived classes + * use this class (e.g. the Mapping_Q() class), the number of shape + * functions may also be different. + * + * In general, it is $(p+1)^\text{dim}$, where $p$ is the polynomial + * degree of the mapping. + */ + const unsigned int n_shape_functions; + + /** + * Tensors of covariant transformation at each of the quadrature points. + * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * + * Jacobian, is the first fundamental form of the map; if dim=spacedim + * then it reduces to the transpose of the inverse of the Jacobian matrix, + * which itself is stored in the @p contravariant field of this structure. + * + * Computed on each cell. + */ + mutable std::vector > covariant; + + /** + * Tensors of contravariant transformation at each of the quadrature + * points. The contravariant matrix is the Jacobian of the transformation, + * i.e. $J_{ij}=dx_i/d\hat x_j$. + * + * Computed on each cell. + */ + mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant; + + /** + * Auxiliary vectors for internal use. + */ + mutable std::vector > > aux; + + /** + * Stores the support points of the mapping shape functions on the @p + * cell_of_current_support_points. + */ + mutable std::vector > mapping_support_points; + + /** + * Stores the cell of which the @p mapping_support_points are stored. + */ + mutable typename Triangulation::cell_iterator cell_of_current_support_points; + + /** + * The determinant of the Jacobian in each quadrature point. Filled if + * #update_volume_elements. + */ + mutable std::vector volume_elements; + }; + + + // documentation can be found in Mapping::requires_update_flags() + virtual + UpdateFlags + requires_update_flags (const UpdateFlags update_flags) const; + + // documentation can be found in Mapping::get_data() + virtual + InternalData * + get_data (const UpdateFlags, + const Quadrature &quadrature) const; + + // documentation can be found in Mapping::get_face_data() + virtual + InternalData * + get_face_data (const UpdateFlags flags, + const Quadrature& quadrature) const; + + // documentation can be found in Mapping::get_subface_data() + virtual + InternalData * + get_subface_data (const UpdateFlags flags, + const Quadrature& quadrature) const; + + // documentation can be found in Mapping::fill_fe_values() + virtual + CellSimilarity::Similarity + fill_fe_values (const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + dealii::internal::FEValues::MappingRelatedData &output_data) const; + + // documentation can be found in Mapping::fill_fe_face_values() + virtual void + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + dealii::internal::FEValues::MappingRelatedData &output_data) const; + + // documentation can be found in Mapping::fill_fe_subface_values() + virtual void + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + dealii::internal::FEValues::MappingRelatedData &output_data) const; + + /** + * @} + */ + +protected: + + /** + * The degree of the polynomials used as shape functions for the mapping of + * cells. + */ + const unsigned int polynomial_degree; + + /* + * The default line support points. These are used when computing + * the location in real space of the support points on lines and + * quads, which are asked to the Manifold class. + * + * The number of quadrature points depends on the degree of this + * class, and it matches the number of degrees of freedom of an + * FE_Q<1>(this->degree). + */ + QGaussLobatto<1> line_support_points; + + /** + * An FE_Q object which is only needed in 3D, since it knows how to reorder + * shape functions/DoFs on non-standard faces. This is used to reorder + * support points in the same way. + */ + const std_cxx11::unique_ptr > fe_q; + + /** + * A table of weights by which we multiply the locations of the support + * points on the perimeter of a quad to get the location of interior support + * points. + * + * Sizes: support_point_weights_on_quad.size()= number of inner + * unit_support_points support_point_weights_on_quad[i].size()= number of + * outer unit_support_points, i.e. unit_support_points on the boundary of + * the quad + * + * For the definition of this vector see equation (8) of the `mapping' + * report. + */ + Table<2,double> support_point_weights_on_quad; + + /** + * A table of weights by which we multiply the locations of the support + * points on the perimeter of a hex to get the location of interior support + * points. + * + * For the definition of this vector see equation (8) of the `mapping' + * report. + */ + Table<2,double> support_point_weights_on_hex; + + /** + * Return the locations of support points for the mapping. For example, for + * $Q_1$ mappings these are the vertices, and for higher order polynomial + * mappings they are the vertices plus interior points on edges, faces, and + * the cell interior that are placed in consultation with the Manifold + * description of the domain and its boundary. However, other classes may + * override this function differently. In particular, the MappingQ1Eulerian + * class does exactly this by not computing the support points from the + * geometry of the current cell but instead evaluating an externally given + * displacement field in addition to the geometry of the cell. + * + * The default implementation of this function is appropriate for most + * cases. It takes the locations of support points on the boundary of the + * cell from the underlying manifold. Interior support points (ie. support + * points in quads for 2d, in hexes for 3d) are then computed using the + * solution of a Laplace equation with the position of the outer support + * points as boundary values, in order to make the transformation as smooth + * as possible. + * + * The function works its way from the vertices (which it takes from the + * given cell) via the support points on the line (for which it calls the + * add_line_support_points() function) and the support points on the quad + * faces (in 3d, for which it calls the add_quad_support_points() function). + * It then adds interior support points that are either computed by + * interpolation from the surrounding points using weights computed by + * solving a Laplace equation, or if dim > + compute_mapping_support_points (const typename Triangulation::cell_iterator &cell) const; + + /** + * Transforms the point @p p on the real cell to the corresponding point on + * the unit cell @p cell by a Newton iteration. + */ + Point + transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit) const; + + /** + * For dim=2,3. Append the support points of all shape functions + * located on bounding lines of the given cell to the vector @p a. Points + * located on the vertices of a line are not included. + * + * Needed by the @p compute_support_points() function. For dim=1 + * this function is empty. The function uses the underlying manifold object + * of the line (or, if none is set, of the cell) for the location of the + * requested points. + * + * This function is made virtual in order to allow derived classes to choose + * shape function support points differently than the present class, which + * chooses the points as interpolation points on the boundary. + */ + virtual + void + add_line_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; + + /** + * For dim=3. Append the support points of all shape functions + * located on bounding faces (quads in 3d) of the given cell to the vector + * @p a. Points located on the vertices or lines of a quad are not included. + * + * Needed by the @p compute_support_points() function. For dim=1 + * and dim=2 this function is empty. The function uses the + * underlying manifold object of the quad (or, if none is set, of the cell) + * for the location of the requested points. + * + * This function is made virtual in order to allow derived classes to choose + * shape function support points differently than the present class, which + * chooses the points as interpolation points on the boundary. + */ + virtual + void + add_quad_support_points(const typename Triangulation::cell_iterator &cell, + std::vector > &a) const; + + /** + * Make MappingQ a friend since it needs to call the fill_fe_values() + * functions on its MappingManifold(1) sub-object. + */ + template friend class MappingQ; +}; + + + +/*@}*/ + +/*----------------------------------------------------------------------*/ + +#ifndef DOXYGEN + +template +inline +const double & +MappingManifold::InternalData::shape (const unsigned int qpoint, + const unsigned int shape_nr) const +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_values.size())); + return shape_values [qpoint*n_shape_functions + shape_nr]; +} + + + +template +inline +double & +MappingManifold::InternalData::shape (const unsigned int qpoint, + const unsigned int shape_nr) +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_values.size())); + return shape_values [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +const Tensor<1,dim> & +MappingManifold::InternalData::derivative (const unsigned int qpoint, + const unsigned int shape_nr) const +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_derivatives.size())); + return shape_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + + +template +inline +Tensor<1,dim> & +MappingManifold::InternalData::derivative (const unsigned int qpoint, + const unsigned int shape_nr) +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_derivatives.size())); + return shape_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +const Tensor<2,dim> & +MappingManifold::InternalData::second_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_second_derivatives.size())); + return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +Tensor<2,dim> & +MappingManifold::InternalData::second_derivative (const unsigned int qpoint, + const unsigned int shape_nr) +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_second_derivatives.size())); + return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; +} + +template +inline +const Tensor<3,dim> & +MappingManifold::InternalData::third_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_third_derivatives.size())); + return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +Tensor<3,dim> & +MappingManifold::InternalData::third_derivative (const unsigned int qpoint, + const unsigned int shape_nr) +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_third_derivatives.size())); + return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +const Tensor<4,dim> & +MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, + const unsigned int shape_nr) const +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_fourth_derivatives.size())); + return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + +template +inline +Tensor<4,dim> & +MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, + const unsigned int shape_nr) +{ + Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), + ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, + shape_fourth_derivatives.size())); + return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; +} + + + +template +inline +bool +MappingManifold::preserves_vertex_locations () const +{ + return true; +} + +#endif // DOXYGEN + +/* -------------- declaration of explicit specializations ------------- */ + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index d00504a3c6..daf002449b 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -57,6 +57,7 @@ SET(_src mapping_q1_eulerian.cc mapping_q.cc mapping_q_eulerian.cc + mapping_manifold.cc ) SET(_inst diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc new file mode 100644 index 0000000000..e777c0e301 --- /dev/null +++ b/source/fe/mapping_manifold.cc @@ -0,0 +1,3894 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2000 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + namespace MappingQ1 + { + namespace + { + + // These are left as templates on the spatial dimension (even though dim + // == spacedim must be true for them to make sense) because templates are + // expanded before the compiler eliminates code due to the 'if (dim == + // spacedim)' statement (see the body of the general + // transform_real_to_unit_cell). + template + Point<1> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<1>::vertices_per_cell> &vertices, + const Point &p) + { + Assert(spacedim == 1, ExcInternalError()); + return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0))); + } + + + + template + Point<2> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<2>::vertices_per_cell> &vertices, + const Point &p) + { + Assert(spacedim == 2, ExcInternalError()); + const double x = p(0); + const double y = p(1); + + const double x0 = vertices[0](0); + const double x1 = vertices[1](0); + const double x2 = vertices[2](0); + const double x3 = vertices[3](0); + + const double y0 = vertices[0](1); + const double y1 = vertices[1](1); + const double y2 = vertices[2](1); + const double y3 = vertices[3](1); + + const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3); + const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1 + - (x - x1)*y2 + (x - x0)*y3; + const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1; + + const double discriminant = b*b - 4*a*c; + // exit if the point is not in the cell (this is the only case where the + // discriminant is negative) + if (discriminant < 0.0) + { + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + } + + double eta1; + double eta2; + // special case #1: if a is zero, then use the linear formula + if (a == 0.0 && b != 0.0) + { + eta1 = -c/b; + eta2 = -c/b; + } + // special case #2: if c is very small or the square root of the + // discriminant is nearly b. + else if (std::abs(c) < 1e-12*std::abs(b) + || std::abs(std::sqrt(discriminant) - b) <= 1e-14*std::abs(b)) + { + eta1 = (-b - std::sqrt(discriminant)) / (2*a); + eta2 = (-b + std::sqrt(discriminant)) / (2*a); + } + // finally, use the numerically stable version of the quadratic formula: + else + { + eta1 = 2*c / (-b - std::sqrt(discriminant)); + eta2 = 2*c / (-b + std::sqrt(discriminant)); + } + // pick the one closer to the center of the cell. + const double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2; + + /* + * There are two ways to compute xi from eta, but either one may have a + * zero denominator. + */ + const double subexpr0 = -eta*x2 + x0*(eta - 1); + const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0; + const double max_x = std::max(std::max(std::abs(x0), std::abs(x1)), + std::max(std::abs(x2), std::abs(x3))); + + if (std::abs(xi_denominator0) > 1e-10*max_x) + { + const double xi = (x + subexpr0)/xi_denominator0; + return Point<2>(xi, eta); + } + else + { + const double max_y = std::max(std::max(std::abs(y0), std::abs(y1)), + std::max(std::abs(y2), std::abs(y3))); + const double subexpr1 = -eta*y2 + y0*(eta - 1); + const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1; + if (std::abs(xi_denominator1) > 1e-10*max_y) + { + const double xi = (subexpr1 + y)/xi_denominator1; + return Point<2>(xi, eta); + } + else // give up and try Newton iteration + { + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + } + } + // bogus return to placate compiler. It should not be possible to get + // here. + Assert(false, ExcInternalError()); + return Point<2>(std::numeric_limits::quiet_NaN(), + std::numeric_limits::quiet_NaN()); + } + + + + template + Point<3> + transform_real_to_unit_cell + (const std_cxx11::array, GeometryInfo<3>::vertices_per_cell> &/*vertices*/, + const Point &/*p*/) + { + // It should not be possible to get here + Assert(false, ExcInternalError()); + return Point<3>(); + } + + + + /** + * Compute an initial guess to pass to the Newton method in + * transform_real_to_unit_cell. For the initial guess we proceed in the + * following way: + *
    + *
  • find the least square dim-dimensional plane approximating the cell + * vertices, i.e. we find an affine map A x_hat + b from the reference cell + * to the real space. + *
  • Solve the equation A x_hat + b = p for x_hat + *
  • This x_hat is the initial solution used for the Newton Method. + *
+ * + * @note if dim + struct TransformR2UInitialGuess + { + static const double KA[GeometryInfo::vertices_per_cell][dim]; + static const double Kb[GeometryInfo::vertices_per_cell]; + }; + + + /* + Octave code: + M=[0 1; 1 1]; + K1 = transpose(M) * inverse (M*transpose(M)); + printf ("{%f, %f},\n", K1' ); + */ + template <> + const double + TransformR2UInitialGuess<1>:: + KA[GeometryInfo<1>::vertices_per_cell][1] = + { + {-1.000000}, + {1.000000} + }; + + template <> + const double + TransformR2UInitialGuess<1>:: + Kb[GeometryInfo<1>::vertices_per_cell] = {1.000000, 0.000000}; + + + /* + Octave code: + M=[0 1 0 1;0 0 1 1;1 1 1 1]; + K2 = transpose(M) * inverse (M*transpose(M)); + printf ("{%f, %f, %f},\n", K2' ); + */ + template <> + const double + TransformR2UInitialGuess<2>:: + KA[GeometryInfo<2>::vertices_per_cell][2] = + { + {-0.500000, -0.500000}, + { 0.500000, -0.500000}, + {-0.500000, 0.500000}, + { 0.500000, 0.500000} + }; + + /* + Octave code: + M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1]; + K3 = transpose(M) * inverse (M*transpose(M)) + printf ("{%f, %f, %f, %f},\n", K3' ); + */ + template <> + const double + TransformR2UInitialGuess<2>:: + Kb[GeometryInfo<2>::vertices_per_cell] = + {0.750000,0.250000,0.250000,-0.250000 }; + + + template <> + const double + TransformR2UInitialGuess<3>:: + KA[GeometryInfo<3>::vertices_per_cell][3] = + { + {-0.250000, -0.250000, -0.250000}, + { 0.250000, -0.250000, -0.250000}, + {-0.250000, 0.250000, -0.250000}, + { 0.250000, 0.250000, -0.250000}, + {-0.250000, -0.250000, 0.250000}, + { 0.250000, -0.250000, 0.250000}, + {-0.250000, 0.250000, 0.250000}, + { 0.250000, 0.250000, 0.250000} + + }; + + + template <> + const double + TransformR2UInitialGuess<3>:: + Kb[GeometryInfo<3>::vertices_per_cell] = + {0.500000,0.250000,0.250000,0.000000,0.250000,0.000000,0.000000,-0.250000}; + + template + Point + transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, + const Point &p) + { + Point p_unit; + + dealii::FullMatrix KA(GeometryInfo::vertices_per_cell, dim); + dealii::Vector Kb(GeometryInfo::vertices_per_cell); + + KA.fill( (double *)(TransformR2UInitialGuess::KA) ); + for (unsigned int i=0; i::vertices_per_cell; ++i) + Kb(i) = TransformR2UInitialGuess::Kb[i]; + + FullMatrix Y(spacedim, GeometryInfo::vertices_per_cell); + for (unsigned int v=0; v::vertices_per_cell; v++) + for (unsigned int i=0; i A(spacedim,dim); + Y.mmult(A,KA); // A = Y*KA + dealii::Vector b(spacedim); + Y.vmult(b,Kb); // b = Y*Kb + + for (unsigned int i=0; i dest(dim); + + FullMatrix A_1(dim,spacedim); + if (dim + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingManifold<1,spacedim>::InternalData &data) + { + (void)n_shape_functions; + const unsigned int n_points=unit_points.size(); + for (unsigned int k = 0 ; k < n_points ; ++k) + { + double x = unit_points[k](0); + + if (data.shape_values.size()!=0) + { + Assert(data.shape_values.size()==n_shape_functions*n_points, + ExcInternalError()); + data.shape(k,0) = 1.-x; + data.shape(k,1) = x; + } + if (data.shape_derivatives.size()!=0) + { + Assert(data.shape_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.derivative(k,0)[0] = -1.; + data.derivative(k,1)[0] = 1.; + } + if (data.shape_second_derivatives.size()!=0) + { + // the following may or may not + // work if dim != spacedim + Assert (spacedim == 1, ExcNotImplemented()); + + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.second_derivative(k,0)[0][0] = 0; + data.second_derivative(k,1)[0][0] = 0; + } + if (data.shape_third_derivatives.size()!=0) + { + // if lower order derivative don't work, neither should this + Assert (spacedim == 1, ExcNotImplemented()); + + Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + + Tensor<3,1> zero; + data.third_derivative(k,0) = zero; + data.third_derivative(k,1) = zero; + } + if (data.shape_fourth_derivatives.size()!=0) + { + // if lower order derivative don't work, neither should this + Assert (spacedim == 1, ExcNotImplemented()); + + Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + + Tensor<4,1> zero; + data.fourth_derivative(k,0) = zero; + data.fourth_derivative(k,1) = zero; + } + } + } + + + template + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingManifold<2,spacedim>::InternalData &data) + { + (void)n_shape_functions; + const unsigned int n_points=unit_points.size(); + for (unsigned int k = 0 ; k < n_points ; ++k) + { + double x = unit_points[k](0); + double y = unit_points[k](1); + + if (data.shape_values.size()!=0) + { + Assert(data.shape_values.size()==n_shape_functions*n_points, + ExcInternalError()); + data.shape(k,0) = (1.-x)*(1.-y); + data.shape(k,1) = x*(1.-y); + data.shape(k,2) = (1.-x)*y; + data.shape(k,3) = x*y; + } + if (data.shape_derivatives.size()!=0) + { + Assert(data.shape_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.derivative(k,0)[0] = (y-1.); + data.derivative(k,1)[0] = (1.-y); + data.derivative(k,2)[0] = -y; + data.derivative(k,3)[0] = y; + data.derivative(k,0)[1] = (x-1.); + data.derivative(k,1)[1] = -x; + data.derivative(k,2)[1] = (1.-x); + data.derivative(k,3)[1] = x; + } + if (data.shape_second_derivatives.size()!=0) + { + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.second_derivative(k,0)[0][0] = 0; + data.second_derivative(k,1)[0][0] = 0; + data.second_derivative(k,2)[0][0] = 0; + data.second_derivative(k,3)[0][0] = 0; + data.second_derivative(k,0)[0][1] = 1.; + data.second_derivative(k,1)[0][1] = -1.; + data.second_derivative(k,2)[0][1] = -1.; + data.second_derivative(k,3)[0][1] = 1.; + data.second_derivative(k,0)[1][0] = 1.; + data.second_derivative(k,1)[1][0] = -1.; + data.second_derivative(k,2)[1][0] = -1.; + data.second_derivative(k,3)[1][0] = 1.; + data.second_derivative(k,0)[1][1] = 0; + data.second_derivative(k,1)[1][1] = 0; + data.second_derivative(k,2)[1][1] = 0; + data.second_derivative(k,3)[1][1] = 0; + } + if (data.shape_third_derivatives.size()!=0) + { + Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + + Tensor<3,2> zero; + for (unsigned int i=0; i<4; ++i) + data.third_derivative(k,i) = zero; + } + if (data.shape_fourth_derivatives.size()!=0) + { + Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + Tensor<4,2> zero; + for (unsigned int i=0; i<4; ++i) + data.fourth_derivative(k,i) = zero; + } + } + } + + + + template + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingManifold<3,spacedim>::InternalData &data) + { + (void)n_shape_functions; + const unsigned int n_points=unit_points.size(); + for (unsigned int k = 0 ; k < n_points ; ++k) + { + double x = unit_points[k](0); + double y = unit_points[k](1); + double z = unit_points[k](2); + + if (data.shape_values.size()!=0) + { + Assert(data.shape_values.size()==n_shape_functions*n_points, + ExcInternalError()); + data.shape(k,0) = (1.-x)*(1.-y)*(1.-z); + data.shape(k,1) = x*(1.-y)*(1.-z); + data.shape(k,2) = (1.-x)*y*(1.-z); + data.shape(k,3) = x*y*(1.-z); + data.shape(k,4) = (1.-x)*(1.-y)*z; + data.shape(k,5) = x*(1.-y)*z; + data.shape(k,6) = (1.-x)*y*z; + data.shape(k,7) = x*y*z; + } + if (data.shape_derivatives.size()!=0) + { + Assert(data.shape_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.derivative(k,0)[0] = (y-1.)*(1.-z); + data.derivative(k,1)[0] = (1.-y)*(1.-z); + data.derivative(k,2)[0] = -y*(1.-z); + data.derivative(k,3)[0] = y*(1.-z); + data.derivative(k,4)[0] = (y-1.)*z; + data.derivative(k,5)[0] = (1.-y)*z; + data.derivative(k,6)[0] = -y*z; + data.derivative(k,7)[0] = y*z; + data.derivative(k,0)[1] = (x-1.)*(1.-z); + data.derivative(k,1)[1] = -x*(1.-z); + data.derivative(k,2)[1] = (1.-x)*(1.-z); + data.derivative(k,3)[1] = x*(1.-z); + data.derivative(k,4)[1] = (x-1.)*z; + data.derivative(k,5)[1] = -x*z; + data.derivative(k,6)[1] = (1.-x)*z; + data.derivative(k,7)[1] = x*z; + data.derivative(k,0)[2] = (x-1)*(1.-y); + data.derivative(k,1)[2] = x*(y-1.); + data.derivative(k,2)[2] = (x-1.)*y; + data.derivative(k,3)[2] = -x*y; + data.derivative(k,4)[2] = (1.-x)*(1.-y); + data.derivative(k,5)[2] = x*(1.-y); + data.derivative(k,6)[2] = (1.-x)*y; + data.derivative(k,7)[2] = x*y; + } + if (data.shape_second_derivatives.size()!=0) + { + // the following may or may not + // work if dim != spacedim + Assert (spacedim == 3, ExcNotImplemented()); + + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + data.second_derivative(k,0)[0][0] = 0; + data.second_derivative(k,1)[0][0] = 0; + data.second_derivative(k,2)[0][0] = 0; + data.second_derivative(k,3)[0][0] = 0; + data.second_derivative(k,4)[0][0] = 0; + data.second_derivative(k,5)[0][0] = 0; + data.second_derivative(k,6)[0][0] = 0; + data.second_derivative(k,7)[0][0] = 0; + data.second_derivative(k,0)[1][1] = 0; + data.second_derivative(k,1)[1][1] = 0; + data.second_derivative(k,2)[1][1] = 0; + data.second_derivative(k,3)[1][1] = 0; + data.second_derivative(k,4)[1][1] = 0; + data.second_derivative(k,5)[1][1] = 0; + data.second_derivative(k,6)[1][1] = 0; + data.second_derivative(k,7)[1][1] = 0; + data.second_derivative(k,0)[2][2] = 0; + data.second_derivative(k,1)[2][2] = 0; + data.second_derivative(k,2)[2][2] = 0; + data.second_derivative(k,3)[2][2] = 0; + data.second_derivative(k,4)[2][2] = 0; + data.second_derivative(k,5)[2][2] = 0; + data.second_derivative(k,6)[2][2] = 0; + data.second_derivative(k,7)[2][2] = 0; + + data.second_derivative(k,0)[0][1] = (1.-z); + data.second_derivative(k,1)[0][1] = -(1.-z); + data.second_derivative(k,2)[0][1] = -(1.-z); + data.second_derivative(k,3)[0][1] = (1.-z); + data.second_derivative(k,4)[0][1] = z; + data.second_derivative(k,5)[0][1] = -z; + data.second_derivative(k,6)[0][1] = -z; + data.second_derivative(k,7)[0][1] = z; + data.second_derivative(k,0)[1][0] = (1.-z); + data.second_derivative(k,1)[1][0] = -(1.-z); + data.second_derivative(k,2)[1][0] = -(1.-z); + data.second_derivative(k,3)[1][0] = (1.-z); + data.second_derivative(k,4)[1][0] = z; + data.second_derivative(k,5)[1][0] = -z; + data.second_derivative(k,6)[1][0] = -z; + data.second_derivative(k,7)[1][0] = z; + + data.second_derivative(k,0)[0][2] = (1.-y); + data.second_derivative(k,1)[0][2] = -(1.-y); + data.second_derivative(k,2)[0][2] = y; + data.second_derivative(k,3)[0][2] = -y; + data.second_derivative(k,4)[0][2] = -(1.-y); + data.second_derivative(k,5)[0][2] = (1.-y); + data.second_derivative(k,6)[0][2] = -y; + data.second_derivative(k,7)[0][2] = y; + data.second_derivative(k,0)[2][0] = (1.-y); + data.second_derivative(k,1)[2][0] = -(1.-y); + data.second_derivative(k,2)[2][0] = y; + data.second_derivative(k,3)[2][0] = -y; + data.second_derivative(k,4)[2][0] = -(1.-y); + data.second_derivative(k,5)[2][0] = (1.-y); + data.second_derivative(k,6)[2][0] = -y; + data.second_derivative(k,7)[2][0] = y; + + data.second_derivative(k,0)[1][2] = (1.-x); + data.second_derivative(k,1)[1][2] = x; + data.second_derivative(k,2)[1][2] = -(1.-x); + data.second_derivative(k,3)[1][2] = -x; + data.second_derivative(k,4)[1][2] = -(1.-x); + data.second_derivative(k,5)[1][2] = -x; + data.second_derivative(k,6)[1][2] = (1.-x); + data.second_derivative(k,7)[1][2] = x; + data.second_derivative(k,0)[2][1] = (1.-x); + data.second_derivative(k,1)[2][1] = x; + data.second_derivative(k,2)[2][1] = -(1.-x); + data.second_derivative(k,3)[2][1] = -x; + data.second_derivative(k,4)[2][1] = -(1.-x); + data.second_derivative(k,5)[2][1] = -x; + data.second_derivative(k,6)[2][1] = (1.-x); + data.second_derivative(k,7)[2][1] = x; + } + if (data.shape_third_derivatives.size()!=0) + { + // if lower order derivative don't work, neither should this + Assert (spacedim == 3, ExcNotImplemented()); + + Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + + for (unsigned int i=0; i<3; ++i) + for (unsigned int j=0; j<3; ++j) + for (unsigned int l=0; l<3; ++l) + if ((i==j)||(j==l)||(l==i)) + { + for (unsigned int m=0; m<8; ++m) + data.third_derivative(k,m)[i][j][l] = 0; + } + else + { + data.third_derivative(k,0)[i][j][l] = -1.; + data.third_derivative(k,1)[i][j][l] = 1.; + data.third_derivative(k,2)[i][j][l] = 1.; + data.third_derivative(k,3)[i][j][l] = -1.; + data.third_derivative(k,4)[i][j][l] = 1.; + data.third_derivative(k,5)[i][j][l] = -1.; + data.third_derivative(k,6)[i][j][l] = -1.; + data.third_derivative(k,7)[i][j][l] = 1.; + } + + } + if (data.shape_fourth_derivatives.size()!=0) + { + // if lower order derivative don't work, neither should this + Assert (spacedim == 3, ExcNotImplemented()); + + Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + Tensor<4,3> zero; + for (unsigned int i=0; i<8; ++i) + data.fourth_derivative(k,i) = zero; + } + } + } + } + } +} + + + + + +template +MappingManifold::InternalData::InternalData (const unsigned int polynomial_degree) + : + polynomial_degree (polynomial_degree), + n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) +{} + + + +template +std::size_t +MappingManifold::InternalData::memory_consumption () const +{ + return (Mapping::InternalDataBase::memory_consumption() + + MemoryConsumption::memory_consumption (shape_values) + + MemoryConsumption::memory_consumption (shape_derivatives) + + MemoryConsumption::memory_consumption (covariant) + + MemoryConsumption::memory_consumption (contravariant) + + MemoryConsumption::memory_consumption (unit_tangentials) + + MemoryConsumption::memory_consumption (aux) + + MemoryConsumption::memory_consumption (mapping_support_points) + + MemoryConsumption::memory_consumption (cell_of_current_support_points) + + MemoryConsumption::memory_consumption (volume_elements) + + MemoryConsumption::memory_consumption (polynomial_degree) + + MemoryConsumption::memory_consumption (n_shape_functions)); +} + + +template +void +MappingManifold::InternalData:: +initialize (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points) +{ + // store the flags in the internal data object so we can access them + // in fill_fe_*_values() + this->update_each = update_flags; + + const unsigned int n_q_points = q.size(); + + // see if we need the (transformation) shape function values + // and/or gradients and resize the necessary arrays + if (this->update_each & update_quadrature_points) + shape_values.resize(n_shape_functions * n_q_points); + + if (this->update_each & (update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values + | update_boundary_forms + | update_normal_vectors + | update_jacobians + | update_jacobian_grads + | update_inverse_jacobians + | update_jacobian_pushed_forward_grads + | update_jacobian_2nd_derivatives + | update_jacobian_pushed_forward_2nd_derivatives + | update_jacobian_3rd_derivatives + | update_jacobian_pushed_forward_3rd_derivatives)) + shape_derivatives.resize(n_shape_functions * n_q_points); + + if (this->update_each & update_covariant_transformation) + covariant.resize(n_original_q_points); + + if (this->update_each & update_contravariant_transformation) + contravariant.resize(n_original_q_points); + + if (this->update_each & update_volume_elements) + volume_elements.resize(n_original_q_points); + + if (this->update_each & + (update_jacobian_grads | update_jacobian_pushed_forward_grads) ) + shape_second_derivatives.resize(n_shape_functions * n_q_points); + + if (this->update_each & + (update_jacobian_2nd_derivatives | update_jacobian_pushed_forward_2nd_derivatives) ) + shape_third_derivatives.resize(n_shape_functions * n_q_points); + + if (this->update_each & + (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) ) + shape_fourth_derivatives.resize(n_shape_functions * n_q_points); + + // now also fill the various fields with their correct values + compute_shape_function_values (q.get_points()); +} + + + +template +void +MappingManifold::InternalData:: +initialize_face (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points) +{ + initialize (update_flags, q, n_original_q_points); + + if (dim > 1) + { + if (this->update_each & update_boundary_forms) + { + aux.resize (dim-1, std::vector > (n_original_q_points)); + + // Compute tangentials to the + // unit cell. + const unsigned int nfaces = GeometryInfo::faces_per_cell; + unit_tangentials.resize (nfaces*(dim-1), + std::vector > (n_original_q_points)); + if (dim==2) + { + // ensure a counterclockwise + // orientation of tangentials + static const int tangential_orientation[4]= {-1,1,1,-1}; + for (unsigned int i=0; i tang; + tang[1-i/2]=tangential_orientation[i]; + std::fill (unit_tangentials[i].begin(), + unit_tangentials[i].end(), tang); + } + } + else if (dim==3) + { + for (unsigned int i=0; i tang1, tang2; + + const unsigned int nd= + GeometryInfo::unit_normal_direction[i]; + + // first tangential + // vector in direction + // of the (nd+1)%3 axis + // and inverted in case + // of unit inward normal + tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; + // second tangential + // vector in direction + // of the (nd+2)%3 axis + tang2[(nd+2)%dim]=1.; + + // same unit tangents + // for all quadrature + // points on this face + std::fill (unit_tangentials[i].begin(), + unit_tangentials[i].end(), tang1); + std::fill (unit_tangentials[nfaces+i].begin(), + unit_tangentials[nfaces+i].end(), tang2); + } + } + } + } +} + + + +namespace +{ + template + std::vector + get_dpo_vector (const unsigned int degree) + { + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +void +MappingManifold::InternalData:: +compute_shape_function_values (const std::vector > &unit_points) +{ + // if the polynomial degree is one, then we can simplify code a bit + // by using hard-coded shape functions. + if ((polynomial_degree == 1) + && + (dim == spacedim)) + internal::MappingQ1::compute_shape_function_values (n_shape_functions, + unit_points, *this); + else + // otherwise ask an object that describes the polynomial space + { + const unsigned int n_points=unit_points.size(); + + // Construct the tensor product polynomials used as shape functions for the + // Qp mapping of cells at the boundary. + const QGaussLobatto<1> line_support_points (polynomial_degree + 1); + const TensorProductPolynomials + tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); + Assert (n_shape_functions==tensor_pols.n(), + ExcInternalError()); + + // then also construct the mapping from lexicographic to the Qp shape function numbering + const std::vector + renumber (FETools:: + lexicographic_to_hierarchic_numbering ( + FiniteElementData (get_dpo_vector(polynomial_degree), 1, + polynomial_degree))); + + std::vector values; + std::vector > grads; + if (shape_values.size()!=0) + { + Assert(shape_values.size()==n_shape_functions*n_points, + ExcInternalError()); + values.resize(n_shape_functions); + } + if (shape_derivatives.size()!=0) + { + Assert(shape_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grads.resize(n_shape_functions); + } + + std::vector > grad2; + if (shape_second_derivatives.size()!=0) + { + Assert(shape_second_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad2.resize(n_shape_functions); + } + + std::vector > grad3; + if (shape_third_derivatives.size()!=0) + { + Assert(shape_third_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad3.resize(n_shape_functions); + } + + std::vector > grad4; + if (shape_fourth_derivatives.size()!=0) + { + Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points, + ExcInternalError()); + grad4.resize(n_shape_functions); + } + + + if (shape_values.size()!=0 || + shape_derivatives.size()!=0 || + shape_second_derivatives.size()!=0 || + shape_third_derivatives.size()!=0 || + shape_fourth_derivatives.size()!=0 ) + for (unsigned int point=0; pointsupport_point_weights_on_quad(hex) arrays. + * + * Called by the compute_support_point_weights_on_quad(hex) functions if the + * data is not yet hardcoded. + * + * For the definition of the support_point_weights_on_quad(hex) please + * refer to equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_laplace_vector(const unsigned int polynomial_degree) + { + Table<2,double> lvs; + + Assert(lvs.n_rows()==0, ExcInternalError()); + Assert(dim==2 || dim==3, ExcNotImplemented()); + + // for degree==1, we shouldn't have to compute any support points, since all + // of them are on the vertices + Assert(polynomial_degree>1, ExcInternalError()); + + const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); + const unsigned int n_outer = (dim==1) ? 2 : + ((dim==2) ? + 4+4*(polynomial_degree-1) : + 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); + + + // compute the shape gradients at the quadrature points on the unit cell + const QGauss quadrature(polynomial_degree+1); + const unsigned int n_q_points=quadrature.size(); + + typename MappingManifold::InternalData quadrature_data(polynomial_degree); + quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * + n_q_points); + quadrature_data.compute_shape_function_values(quadrature.get_points()); + + // Compute the stiffness matrix of the inner dofs + FullMatrix S(n_inner); + for (unsigned int point=0; point T(n_inner, n_outer); + for (unsigned int point=0; point S_1(n_inner); + S_1.invert(S); + + FullMatrix S_1_T(n_inner, n_outer); + + // S:=S_1*T + S_1.mmult(S_1_T,T); + + // Resize and initialize the lvs + lvs.reinit (n_inner, n_outer); + for (unsigned int i=0; iMappingQ for dim= 2 and 3. + * + * For degree<4 this function sets the @p support_point_weights_on_quad to + * the hardcoded data. For degree>=4 and MappingQ<2> this vector is + * computed. + * + * For the definition of the @p support_point_weights_on_quad please refer to + * equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_support_point_weights_on_quad(const unsigned int polynomial_degree) + { + Table<2,double> loqvs; + + // in 1d, there are no quads, so return an empty object + if (dim == 1) + return loqvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return loqvs; + + const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1); + const unsigned int n_outer_2d=4+4*(polynomial_degree-1); + + // first check whether we have precomputed the values for some polynomial + // degree; the sizes of arrays is n_inner_2d*n_outer_2d + if (polynomial_degree == 2) + { + // (checked these values against the output of compute_laplace_vector + // again, and found they're indeed right -- just in case someone wonders + // where they come from -- WB) + static const double loqv2[1*8] + = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; + Assert (sizeof(loqv2)/sizeof(loqv2[0]) == + n_inner_2d * n_outer_2d, + ExcInternalError()); + + // copy and return + loqvs.reinit(n_inner_2d, n_outer_2d); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } + + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_pointMappingQ<3>. + * + * For degree==2 this function sets the @p support_point_weights_on_hex to + * the hardcoded data. For degree>2 this vector is computed. + * + * For the definition of the @p support_point_weights_on_hex please refer to + * equation (8) of the `mapping' report. + */ + template + Table<2,double> + compute_support_point_weights_on_hex(const unsigned int polynomial_degree) + { + Table<2,double> lohvs; + + // in 1d and 2d, there are no hexes, so return an empty object + if (dim < 3) + return lohvs; + + // we are asked to compute weights for interior support points, but + // there are no interior points if degree==1 + if (polynomial_degree == 1) + return lohvs; + + const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); + const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1); + + // first check whether we have precomputed the values for some polynomial + // degree; the sizes of arrays is n_inner_2d*n_outer_2d + if (polynomial_degree == 2) + { + static const double lohv2[26] + = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., + 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., + 7/192., 7/192., 7/192., 7/192., + 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. + }; + + // copy and return + lohvs.reinit(n_inner, n_outer); + for (unsigned int unit_point=0; unit_point(polynomial_degree); + } + + // the sum of weights of the points at the outer rim should be one. check + // this + for (unsigned int unit_point=0; unit_point +MappingManifold::MappingManifold (const unsigned int p) + : + polynomial_degree(p), + line_support_points(this->polynomial_degree+1), + fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : 0), + support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), + support_point_weights_on_hex (compute_support_point_weights_on_hex(this->polynomial_degree)) +{ + Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings " + "with a polynomial degree greater or equal to one.")); +} + + + +template +MappingManifold::MappingManifold (const MappingManifold &mapping) + : + polynomial_degree(mapping.polynomial_degree), + line_support_points(mapping.line_support_points), + fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : 0), + support_point_weights_on_quad (mapping.support_point_weights_on_quad), + support_point_weights_on_hex (mapping.support_point_weights_on_hex) +{} + + + + +template +Mapping * +MappingManifold::clone () const +{ + return new MappingManifold(*this); +} + + + + +template +unsigned int +MappingManifold::get_degree() const +{ + return polynomial_degree; +} + + + +template +Point +MappingManifold:: +transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + // set up the polynomial space + const QGaussLobatto<1> line_support_points (polynomial_degree + 1); + const TensorProductPolynomials + tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); + Assert (tensor_pols.n() == Utilities::fixed_power(polynomial_degree+1), + ExcInternalError()); + + // then also construct the mapping from lexicographic to the Qp shape function numbering + const std::vector + renumber (FETools:: + lexicographic_to_hierarchic_numbering ( + FiniteElementData (get_dpo_vector(polynomial_degree), 1, + polynomial_degree))); + + const std::vector > support_points + = this->compute_mapping_support_points(cell); + + Point mapped_point; + for (unsigned int i=0; i when +// seeing which of the overloaded versions of +// do_transform_real_to_unit_cell_internal() to call. This leads to bad +// error messages and, generally, nothing very good. Avoid this by ensuring +// that this class exists, but does not have an inner InternalData +// type, thereby ruling out the codim-1 version of the function +// below when doing overload resolution. +template <> +class MappingManifold<3,4> +{}; + +namespace +{ + /** + * Using the relative weights of the shape functions evaluated at + * one point on the reference cell (and stored in data.shape_values + * and accessed via data.shape(0,i)) and the locations of mapping + * support points (stored in data.mapping_support_points), compute + * the mapped location of that point in real space. + */ + template + Point + compute_mapped_location_of_point (const typename MappingManifold::InternalData &data) + { + AssertDimension (data.shape_values.size(), + data.mapping_support_points.size()); + + // use now the InternalData to compute the point in real space. + Point p_real; + for (unsigned int i=0; i + Point + do_transform_real_to_unit_cell_internal + (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename MappingManifold::InternalData &mdata) + { + const unsigned int spacedim = dim; + + const unsigned int n_shapes=mdata.shape_values.size(); + (void)n_shapes; + Assert(n_shapes!=0, ExcInternalError()); + AssertDimension (mdata.shape_derivatives.size(), n_shapes); + + std::vector > &points=mdata.mapping_support_points; + AssertDimension (points.size(), n_shapes); + + + // Newton iteration to solve + // f(x)=p(x)-p=0 + // where we are looking for 'x' and p(x) is the forward transformation + // from unit to real cell. We solve this using a Newton iteration + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + // The start value is set to be the linear approximation to the cell + + // The shape values and derivatives of the mapping at this point are + // previously computed. + + Point p_unit = initial_p_unit; + + mdata.compute_shape_function_values(std::vector > (1, p_unit)); + + Point p_real = compute_mapped_location_of_point(mdata); + Tensor<1,spacedim> f = p_real-p; + + // early out if we already have our point + if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) + return p_unit; + + // we need to compare the position of the computed p(x) against the given + // point 'p'. We will terminate the iteration and return 'x' if they are + // less than eps apart. The question is how to choose eps -- or, put maybe + // more generally: in which norm we want these 'p' and 'p(x)' to be eps + // apart. + // + // the question is difficult since we may have to deal with very elongated + // cells where we may achieve 1e-12*h for the distance of these two points + // in the 'long' direction, but achieving this tolerance in the 'short' + // direction of the cell may not be possible + // + // what we do instead is then to terminate iterations if + // \| p(x) - p \|_A < eps + // where the A-norm is somehow induced by the transformation of the cell. + // in particular, we want to measure distances relative to the sizes of + // the cell in its principal directions. + // + // to define what exactly A should be, note that to first order we have + // the following (assuming that x* is the solution of the problem, i.e., + // p(x*)=p): + // p(x) - p = p(x) - p(x*) + // = -grad p(x) * (x*-x) + higher order terms + // This suggest to measure with a norm that corresponds to + // A = {[grad p(x]^T [grad p(x)]}^{-1} + // because then + // \| p(x) - p \|_A \approx \| x - x* \| + // Consequently, we will try to enforce that + // \| p(x) - p \|_A = \| f \| <= eps + // + // Note that using this norm is a bit dangerous since the norm changes + // in every iteration (A isn't fixed by depends on xk). However, if the + // cell is not too deformed (it may be stretched, but not twisted) then + // the mapping is almost linear and A is indeed constant or nearly so. + const double eps = 1.e-11; + const unsigned int newton_iteration_limit = 20; + + unsigned int newton_iteration = 0; + double last_f_weighted_norm; + do + { +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << "Newton iteration " << newton_iteration << std::endl; +#endif + + // f'(x) + Tensor<2,spacedim> df; + for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); + const Point &point=points[k]; + + for (unsigned int i=0; i df_inverse = invert(df); + const Tensor<1,spacedim> delta = df_inverse * static_cast&>(f); + +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << " delta=" << delta << std::endl; +#endif + + // do a line search + double step_length = 1; + do + { + // update of p_unit. The spacedim-th component of transformed point + // is simply ignored in codimension one case. When this component is + // not zero, then we are projecting the point to the surface or + // curve identified by the cell. + Point p_unit_trial = p_unit; + for (unsigned int i=0; i > (1, p_unit_trial)); + + // f(x) + Point p_real_trial = compute_mapped_location_of_point(mdata); + const Tensor<1,spacedim> f_trial = p_real_trial-p; + +#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL + std::cout << " step_length=" << step_length << std::endl + << " ||f || =" << f.norm() << std::endl + << " ||f*|| =" << f_trial.norm() << std::endl + << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; +#endif + + // see if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again + // + // strictly speaking, we should probably use the same norm as we use + // for the outer algorithm. in practice, line search is just a + // crutch to find a "reasonable" step length, and so using the l2 + // norm is probably just fine + if (f_trial.norm() < f.norm()) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + f = f_trial; + break; + } + else if (step_length > 0.05) + step_length /= 2; + else + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + } + while (true); + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) + AssertThrow (false, + (typename Mapping::ExcTransformationFailed())); + last_f_weighted_norm = (df_inverse * f).norm(); + } + while (last_f_weighted_norm > eps); + + return p_unit; + } + + + + /** + * Implementation of transform_real_to_unit_cell for dim==spacedim-1 + */ + template + Point + do_transform_real_to_unit_cell_internal_codim1 + (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + typename MappingManifold::InternalData &mdata) + { + const unsigned int spacedim = dim+1; + + const unsigned int n_shapes=mdata.shape_values.size(); + (void)n_shapes; + Assert(n_shapes!=0, ExcInternalError()); + Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); + Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); + + std::vector > &points=mdata.mapping_support_points; + Assert(points.size()==n_shapes, ExcInternalError()); + + Point p_minus_F; + + Tensor<1,spacedim> DF[dim]; + Tensor<1,spacedim> D2F[dim][dim]; + + Point p_unit = initial_p_unit; + Point f; + Tensor<2,dim> df; + + // Evaluate first and second derivatives + mdata.compute_shape_function_values(std::vector > (1, p_unit)); + + for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); + const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); + const Point &point_k = points[k]; + + for (unsigned int j=0; j(mdata); + + + for (unsigned int j=0; jdiameter(); + const unsigned int loop_limit = 10; + + unsigned int loop=0; + + while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(f); + p_unit -= d; + + for (unsigned int j=0; j > (1, p_unit)); + + for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); + const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); + const Point &point_k = points[k]; + + for (unsigned int j=0; j(mdata); + + for (unsigned int j=0; j::ExcTransformationFailed())); + + return p_unit; + } + + +} + + + +// visual studio freaks out when trying to determine if +// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good +// candidate. So instead of letting the compiler pick the correct overload, we +// use template specialization to make sure we pick up the right function to +// call: + +template +Point +MappingManifold:: +transform_real_to_unit_cell_internal +(const typename Triangulation::cell_iterator &, + const Point &, + const Point &) const +{ + // default implementation (should never be called) + Assert(false, ExcInternalError()); + return Point(); +} + +template<> +Point<1> +MappingManifold<1,1>:: +transform_real_to_unit_cell_internal +(const Triangulation<1,1>::cell_iterator &cell, + const Point<1> &p, + const Point<1> &initial_p_unit) const +{ + const int dim = 1; + const int spacedim = 1; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); + + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata); +} + +template<> +Point<2> +MappingManifold<2, 2>:: +transform_real_to_unit_cell_internal +(const Triangulation<2, 2>::cell_iterator &cell, + const Point<2> &p, + const Point<2> &initial_p_unit) const +{ + const int dim = 2; + const int spacedim = 2; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); + + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata); +} + +template<> +Point<3> +MappingManifold<3, 3>:: +transform_real_to_unit_cell_internal +(const Triangulation<3, 3>::cell_iterator &cell, + const Point<3> &p, + const Point<3> &initial_p_unit) const +{ + const int dim = 3; + const int spacedim = 3; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); + + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata); +} + +template<> +Point<1> +MappingManifold<1, 2>:: +transform_real_to_unit_cell_internal +(const Triangulation<1, 2>::cell_iterator &cell, + const Point<2> &p, + const Point<1> &initial_p_unit) const +{ + const int dim = 1; + const int spacedim = 2; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); + + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata); +} + +template<> +Point<2> +MappingManifold<2, 3>:: +transform_real_to_unit_cell_internal +(const Triangulation<2, 3>::cell_iterator &cell, + const Point<3> &p, + const Point<2> &initial_p_unit) const +{ + const int dim = 2; + const int spacedim = 3; + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_quadrature_points | update_jacobians; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std_cxx11::unique_ptr mdata (get_data(update_flags, + point_quadrature)); + + mdata->mapping_support_points = this->compute_mapping_support_points (cell); + + // dispatch to the various specializations for spacedim=dim, + // spacedim=dim+1, etc + return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata); +} + +template<> +Point<1> +MappingManifold<1, 3>:: +transform_real_to_unit_cell_internal +(const Triangulation<1, 3>::cell_iterator &, + const Point<3> &, + const Point<1> &) const +{ + Assert (false, ExcNotImplemented()); + return Point<1>(); +} + + + +template +Point +MappingManifold:: +transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + // Use an exact formula if one is available. this is only the case + // for Q1 mappings in 1d, and in 2d if dim==spacedim + if ((polynomial_degree == 1) && + ((dim == 1) + || + ((dim == 2) && (dim == spacedim)))) + { + // The dimension-dependent algorithms are much faster (about 25-45x in + // 2D) but fail most of the time when the given point (p) is not in the + // cell. The dimension-independent Newton algorithm given below is + // slower, but more robust (though it still sometimes fails). Therefore + // this function implements the following strategy based on the + // p's dimension: + // + // * In 1D this mapping is linear, so the mapping is always invertible + // (and the exact formula is known) as long as the cell has non-zero + // length. + // * In 2D the exact (quadratic) formula is called first. If either the + // exact formula does not succeed (negative discriminant in the + // quadratic formula) or succeeds but finds a solution outside of the + // unit cell, then the Newton solver is called. The rationale for the + // second choice is that the exact formula may provide two different + // answers when mapping a point outside of the real cell, but the + // Newton solver (if it converges) will only return one answer. + // Otherwise the exact formula successfully found a point in the unit + // cell and that value is returned. + // * In 3D there is no (known to the authors) exact formula, so the Newton + // algorithm is used. + const std_cxx11::array, GeometryInfo::vertices_per_cell> + vertices = this->get_vertices(cell); + try + { + switch (dim) + { + case 1: + { + // formula not subject to any issues in 1d + if (spacedim == 1) + return internal::MappingQ1::transform_real_to_unit_cell(vertices, p); + else + { + const std::vector > a (vertices.begin(), + vertices.end()); + return internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + } + } + + case 2: + { + const Point point + = internal::MappingQ1::transform_real_to_unit_cell(vertices, p); + + // formula not guaranteed to work for points outside of + // the cell. only take the computed point if it lies + // inside the reference cell + const double eps = 1e-15; + if (-eps <= point(1) && point(1) <= 1 + eps && + -eps <= point(0) && point(0) <= 1 + eps) + { + return point; + } + else + break; + } + + default: + { + // we should get here, based on the if-condition at the top + Assert(false, ExcInternalError()); + } + } + } + catch (const typename Mapping::ExcTransformationFailed &) + { + // simply fall through and continue on to the standard Newton code + } + } + else + { + // we can't use an explicit formula, + } + + + Point initial_p_unit; + if (polynomial_degree == 1) + { + // Find the initial value for the Newton iteration by a normal + // projection to the least square plane determined by the vertices + // of the cell + const std::vector > a + = this->compute_mapping_support_points (cell); + Assert(a.size() == GeometryInfo::vertices_per_cell, + ExcInternalError()); + initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + } + else + { + try + { + // Find the initial value for the Newton iteration by a normal + // projection to the least square plane determined by the vertices + // of the cell + // + // we do this by first getting all support points, then + // throwing away all but the vertices, and finally calling + // the same function as above + std::vector > a + = this->compute_mapping_support_points (cell); + a.resize(GeometryInfo::vertices_per_cell); + initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + for (unsigned int d=0; d::project_to_unit_cell(initial_p_unit); + } + + // perform the Newton iteration and return the result. note that + // this statement may throw an exception, which we simply pass up to + // the caller + return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit); +} + + + +template +UpdateFlags +MappingManifold::requires_update_flags (const UpdateFlags in) const +{ + // add flags if the respective quantities are necessary to compute + // what we need. note that some flags appear in both the conditions + // and in subsequent set operations. this leads to some circular + // logic. the only way to treat this is to iterate. since there are + // 5 if-clauses in the loop, it will take at most 5 iterations to + // converge. do them: + UpdateFlags out = in; + for (unsigned int i=0; i<5; ++i) + { + // The following is a little incorrect: + // If not applied on a face, + // update_boundary_forms does not + // make sense. On the other hand, + // it is necessary on a + // face. Currently, + // update_boundary_forms is simply + // ignored for the interior of a + // cell. + if (out & (update_JxW_values + | update_normal_vectors)) + out |= update_boundary_forms; + + if (out & (update_covariant_transformation + | update_JxW_values + | update_jacobians + | update_jacobian_grads + | update_boundary_forms + | update_normal_vectors)) + out |= update_contravariant_transformation; + + if (out & (update_inverse_jacobians + | update_jacobian_pushed_forward_grads + | update_jacobian_pushed_forward_2nd_derivatives + | update_jacobian_pushed_forward_3rd_derivatives) ) + out |= update_covariant_transformation; + + // The contravariant transformation + // used in the Piola transformation, which + // requires the determinant of the + // Jacobi matrix of the transformation. + // Because we have no way of knowing here whether the finite + // elements wants to use the contravariant of the Piola + // transforms, we add the JxW values to the list of flags to be + // updated for each cell. + if (out & update_contravariant_transformation) + out |= update_JxW_values; + + if (out & update_normal_vectors) + out |= update_JxW_values; + } + + return out; +} + + + +template +typename MappingManifold::InternalData * +MappingManifold::get_data (const UpdateFlags update_flags, + const Quadrature &q) const +{ + InternalData *data = new InternalData(polynomial_degree); + data->initialize (this->requires_update_flags(update_flags), q, q.size()); + + return data; +} + + + +template +typename MappingManifold::InternalData * +MappingManifold::get_face_data (const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + InternalData *data = new InternalData(polynomial_degree); + data->initialize_face (this->requires_update_flags(update_flags), + QProjector::project_to_all_faces(quadrature), + quadrature.size()); + + return data; +} + + + +template +typename MappingManifold::InternalData * +MappingManifold::get_subface_data (const UpdateFlags update_flags, + const Quadrature& quadrature) const +{ + InternalData *data = new InternalData(polynomial_degree); + data->initialize_face (this->requires_update_flags(update_flags), + QProjector::project_to_all_subfaces(quadrature), + quadrature.size()); + + return data; +} + + + +namespace internal +{ + namespace + { + /** + * Compute the locations of quadrature points on the object described by + * the first argument (and the cell for which the mapping support points + * have already been set), but only if the update_flags of the @p data + * argument indicate so. + */ + template + void + maybe_compute_q_points (const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_quadrature_points) + { + for (unsigned int point=0; point result = (shape[0] * + data.mapping_support_points[0]); + for (unsigned int k=1; k + void + maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_contravariant_transformation) + // if the current cell is just a + // translation of the previous one, no + // need to recompute jacobians... + if (cell_similarity != CellSimilarity::translation) + { + const unsigned int n_q_points = data.contravariant.size(); + + std::fill(data.contravariant.begin(), data.contravariant.end(), + DerivativeForm<1,dim,spacedim>()); + + Assert (data.n_shape_functions > 0, ExcInternalError()); + const Tensor<1,spacedim> *supp_pts = + &data.mapping_support_points[0]; + + for (unsigned int point=0; point *data_derv = + &data.derivative(point+data_set, 0); + + double result [spacedim][dim]; + + // peel away part of sum to avoid zeroing the + // entries and adding for the first time + for (unsigned int i=0; i + void + maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_pushed_forward_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_grads) + { + const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_pushed_forward_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *third = + &data.third_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim][dim]; + for (unsigned int i=0; i + void + maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingManifold::InternalData &data, + std::vector > &jacobian_pushed_forward_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); + + if (cell_similarity != CellSimilarity::translation) + { + double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point=0; point *fourth = + &data.fourth_derivative(point+data_set, 0); + double result [spacedim][dim][dim][dim][dim]; + for (unsigned int i=0; i +CellSimilarity::Similarity +MappingManifold:: +fill_fe_values (const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValues::MappingRelatedData &output_data) const +{ + // ensure that the following static_cast is really correct: + Assert (dynamic_cast(&internal_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + const unsigned int n_q_points=quadrature.size(); + + // if necessary, recompute the support points of the transformation of this cell + // (note that we need to first check the triangulation pointer, since otherwise + // the second test might trigger an exception if the triangulations are not the + // same) + if ((data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &data.cell_of_current_support_points->get_triangulation()) + || + (cell != data.cell_of_current_support_points)) + { + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; + } + + internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), + data, + output_data.quadrature_points); + internal::maybe_update_Jacobians (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data); + + const UpdateFlags update_flags = data.update_each; + const std::vector &weights=quadrature.get_weights(); + + // Multiply quadrature weights by absolute value of Jacobian determinants or + // the area element g=sqrt(DX^t DX) in case of codim > 0 + + if (update_flags & (update_normal_vectors + | update_JxW_values)) + { + AssertDimension (output_data.JxW_values.size(), n_q_points); + + Assert( !(update_flags & update_normal_vectors ) || + (output_data.normal_vectors.size() == n_q_points), + ExcDimensionMismatch(output_data.normal_vectors.size(), n_q_points)); + + + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point=0; point 1e-12*Utilities::fixed_power(cell->diameter()/ + std::sqrt(double(dim))), + (typename Mapping::ExcDistortedMappedCell(cell->center(), det, point))); + + output_data.JxW_values[point] = weights[point] * det; + } + // if dim==spacedim, then there is no cell normal to + // compute. since this is for FEValues (and not FEFaceValues), + // there are also no face normals to compute + else //codim>0 case + { + Tensor<1, spacedim> DX_t [dim]; + for (unsigned int i=0; i G; //First fundamental form + for (unsigned int i=0; idirection_flag() == false) + output_data.normal_vectors[point] *= -1.; + } + + } + } //codim>0 case + + } + } + + + + // copy values from InternalData to vector given by reference + if (update_flags & update_jacobians) + { + AssertDimension (output_data.jacobians.size(), n_q_points); + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point=0; point (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_grads); + + internal::maybe_update_jacobian_pushed_forward_grads (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_grads); + + internal::maybe_update_jacobian_2nd_derivatives (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_2nd_derivatives); + + internal::maybe_update_jacobian_pushed_forward_2nd_derivatives (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + + internal::maybe_update_jacobian_3rd_derivatives (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_3rd_derivatives); + + internal::maybe_update_jacobian_pushed_forward_3rd_derivatives (cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_3rd_derivatives); + + return cell_similarity; +} + + + + + + +namespace internal +{ + namespace + { + /** + * Depending on what information is called for in the update flags of the + * @p data object, compute the various pieces of information that is required + * by the fill_fe_face_values() and fill_fe_subface_values() functions. + * This function simply unifies the work that would be done by + * those two functions. + * + * The resulting data is put into the @p output_data argument. + */ + template + void + maybe_compute_face_data (const dealii::MappingManifold &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points, + const std::vector &weights, + const typename dealii::MappingManifold::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_boundary_forms) + { + AssertDimension (output_data.boundary_forms.size(), n_q_points); + if (update_flags & update_normal_vectors) + AssertDimension (output_data.normal_vectors.size(), n_q_points); + if (update_flags & update_JxW_values) + AssertDimension (output_data.JxW_values.size(), n_q_points); + + // map the unit tangentials to the real cell. checking for d!=dim-1 + // eliminates compiler warnings regarding unsigned int expressions < + // 0. + for (unsigned int d=0; d!=dim-1; ++d) + { + Assert (face_no+GeometryInfo::faces_per_cell*d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert (data.aux[d].size() <= + data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), + ExcInternalError()); + + mapping.transform (make_array_view(data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } + + // if dim==spacedim, we can use the unit tangentials to compute the + // boundary form by simply taking the cross product + if (dim == spacedim) + { + for (unsigned int i=0; i DX_t = + data.contravariant[point].transpose(); + + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); + + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } + } + } + + if (update_flags & (update_normal_vectors + | update_JxW_values)) + for (unsigned int i=0; i::subface_ratio( + cell->subface_case(face_no), subface_no); + output_data.JxW_values[i] *= area_ratio; + } + } + + if (update_flags & update_normal_vectors) + output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); + } + + if (update_flags & update_jacobians) + for (unsigned int point=0; point + void + do_fill_fe_face_values (const dealii::MappingManifold &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename QProjector::DataSetDescriptor data_set, + const Quadrature &quadrature, + const typename dealii::MappingManifold::InternalData &data, + internal::FEValues::MappingRelatedData &output_data) + { + maybe_compute_q_points (data_set, + data, + output_data.quadrature_points); + maybe_update_Jacobians (CellSimilarity::none, + data_set, + data); + maybe_update_jacobian_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_grads); + maybe_update_jacobian_pushed_forward_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_grads); + maybe_update_jacobian_2nd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_2nd_derivatives); + maybe_update_jacobian_pushed_forward_2nd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + maybe_update_jacobian_3rd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_3rd_derivatives); + maybe_update_jacobian_pushed_forward_3rd_derivatives (CellSimilarity::none, + data_set, + data, + output_data.jacobian_pushed_forward_3rd_derivatives); + + maybe_compute_face_data (mapping, + cell, face_no, subface_no, quadrature.size(), + quadrature.get_weights(), data, + output_data); + } + } +} + + + +template +void +MappingManifold:: +fill_fe_face_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValues::MappingRelatedData &output_data) const +{ + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != 0), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + // if necessary, recompute the support points of the transformation of this cell + // (note that we need to first check the triangulation pointer, since otherwise + // the second test might trigger an exception if the triangulations are not the + // same) + if ((data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &data.cell_of_current_support_points->get_triangulation()) + || + (cell != data.cell_of_current_support_points)) + { + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; + } + + internal::do_fill_fe_face_values (*this, + cell, face_no, numbers::invalid_unsigned_int, + QProjector::DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size()), + quadrature, + data, + output_data); +} + + + +template +void +MappingManifold:: +fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValues::MappingRelatedData &output_data) const +{ + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != 0), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + // if necessary, recompute the support points of the transformation of this cell + // (note that we need to first check the triangulation pointer, since otherwise + // the second test might trigger an exception if the triangulations are not the + // same) + if ((data.mapping_support_points.size() == 0) + || + (&cell->get_triangulation() != + &data.cell_of_current_support_points->get_triangulation()) + || + (cell != data.cell_of_current_support_points)) + { + data.mapping_support_points = this->compute_mapping_support_points(cell); + data.cell_of_current_support_points = cell; + } + + internal::do_fill_fe_face_values (*this, + cell, face_no, subface_no, + QProjector::DataSetDescriptor::subface (face_no, subface_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)), + quadrature, + data, + output_data); +} + + + +namespace +{ + template + void + transform_fields(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_contravariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + if (rank!=1) + return; + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i + void + transform_gradients(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_contravariant_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); + + for (unsigned int i=0; i A = + apply_transformation(data.contravariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } + + return; + } + + case mapping_covariant_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (rank==2, ExcMessage("Only for rank 2")); + + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], transpose(input[i]) ); + output[i] = apply_transformation(data.covariant[i], A.transpose() ); + } + + return; + } + + case mapping_piola_gradient: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==2, ExcMessage("Only for rank 2")); + + for (unsigned int i=0; i A = + apply_transformation(data.covariant[i], input[i] ); + Tensor<2,spacedim> T = + apply_transformation(data.contravariant[i], A.transpose() ); + + output[i] = transpose(T); + output[i] /= data.volume_elements[i]; + } + + return; + } + + default: + Assert(false, ExcNotImplemented()); + } + } + + + + + template + void + transform_hessians(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_contravariant_hessian: + { + Assert (data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int q=0; q::ExcAccessToUninitializedField("update_covariant_transformation")); + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + + for (unsigned int q=0; q + void + transform_differential_forms(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_covariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i +void +MappingManifold:: +transform (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) const +{ + transform_fields(input, mapping_type, mapping_data, output); +} + + + +template +void +MappingManifold:: +transform (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) const +{ + transform_differential_forms(input, mapping_type, mapping_data, output); +} + + + +template +void +MappingManifold:: +transform (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) const +{ + switch (mapping_type) + { + case mapping_contravariant: + transform_fields(input, mapping_type, mapping_data, output); + return; + + case mapping_piola_gradient: + case mapping_contravariant_gradient: + case mapping_covariant_gradient: + transform_gradients(input, mapping_type, mapping_data, output); + return; + default: + Assert(false, ExcNotImplemented()); + } +} + + + +template +void +MappingManifold:: +transform (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) const +{ + + AssertDimension (input.size(), output.size()); + Assert (dynamic_cast(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_type) + { + case mapping_covariant_gradient: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int q=0; q +void +MappingManifold:: +transform (const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) const +{ + switch (mapping_type) + { + case mapping_piola_hessian: + case mapping_contravariant_hessian: + case mapping_covariant_hessian: + transform_hessians(input, mapping_type, mapping_data, output); + return; + default: + Assert(false, ExcNotImplemented()); + } +} + + + +namespace +{ + /** + * Ask the manifold descriptor to return intermediate points on lines or + * faces. The function needs to return one or multiple points (depending on + * the number of elements in the output vector @p points that lie inside a + * line, quad or hex). Whether it is a line, quad or hex doesn't really + * matter to this function but it can be inferred from the number of input + * points in the @p surrounding_points vector. + */ + template + void + get_intermediate_points (const Manifold &manifold, + const QGaussLobatto<1> &line_support_points, + const std::vector > &surrounding_points, + std::vector > &points) + { + Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required")); + const unsigned int n=points.size(); + Assert(n>0, ExcMessage("You can't ask for 0 intermediate points.")); + std::vector w(surrounding_points.size()); + + switch (surrounding_points.size()) + { + case 2: + { + // If two points are passed, these are the two vertices, and + // we can only compute degree-1 intermediate points. + for (unsigned int i=0; i quadrature(surrounding_points, w); + points[i] = manifold.get_new_point(quadrature); + } + break; + } + + case 4: + { + Assert(spacedim >= 2, ExcImpossibleInDim(spacedim)); + const unsigned m= + static_cast(std::sqrt(static_cast(n))); + // is n a square number + Assert(m*m==n, ExcInternalError()); + + // If four points are passed, these are the two vertices, and + // we can only compute (degree-1)*(degree-1) intermediate + // points. + for (unsigned int i=0; i quadrature(surrounding_points, w); + points[i*m+j]=manifold.get_new_point(quadrature); + } + } + break; + } + + case 8: + Assert(false, ExcNotImplemented()); + break; + default: + Assert(false, ExcInternalError()); + break; + } + } + + + + + /** + * Ask the manifold descriptor to return intermediate points on the object + * pointed to by the TriaIterator @p iter. This function tries to be + * backward compatible with respect to the differences between + * Boundary and Manifold, querying the first + * whenever the passed @p manifold can be upgraded to a + * Boundary. + */ + template + void get_intermediate_points_on_object(const Manifold &manifold, + const QGaussLobatto<1> &line_support_points, + const TriaIterator &iter, + std::vector > &points) + { + const unsigned int structdim = TriaIterator::AccessorType::structure_dimension; + + // Try backward compatibility option. + if (const Boundary *boundary + = dynamic_cast *>(&manifold)) + // This is actually a boundary. Call old methods. + { + switch (structdim) + { + case 1: + { + const typename Triangulation::line_iterator line = iter; + boundary->get_intermediate_points_on_line(line, points); + return; + } + case 2: + { + const typename Triangulation::quad_iterator quad = iter; + boundary->get_intermediate_points_on_quad(quad, points); + return; + } + default: + Assert(false, ExcInternalError()); + return; + } + } + else + { + std::vector > sp(GeometryInfo::vertices_per_cell); + for (unsigned int i=0; ivertex(i); + get_intermediate_points(manifold, line_support_points, sp, points); + } + } + + + /** + * Take a support_point_weights_on_hex(quad) and apply it to the vector + * @p a to compute the inner support points as a linear combination of the + * exterior points. + * + * The vector @p a initially contains the locations of the @p n_outer + * points, the @p n_inner computed inner points are appended. + * + * See equation (7) of the `mapping' report. + */ + template + void add_weighted_interior_points(const Table<2,double> &lvs, + std::vector > &a) + { + const unsigned int n_inner_apply=lvs.n_rows(); + const unsigned int n_outer_apply=lvs.n_cols(); + Assert(a.size()==n_outer_apply, + ExcDimensionMismatch(a.size(), n_outer_apply)); + + // compute each inner point as linear combination of the outer points. the + // weights are given by the lvs entries, the outer points are the first + // (existing) elements of a + for (unsigned int unit_point=0; unit_point p; + for (unsigned int k=0; k +void +MappingManifold:: +add_line_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const +{ + // if we only need the midpoint, then ask for it. + if (this->polynomial_degree==2) + { + for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) + { + const typename Triangulation::line_iterator line = + (dim == 1 ? + static_cast::line_iterator>(cell) : + cell->line(line_no)); + + const Manifold &manifold = + ( ( line->manifold_id() == numbers::invalid_manifold_id ) && + ( dim < spacedim ) + ? + cell->get_manifold() + : + line->get_manifold() ); + a.push_back(manifold.get_new_point_on_line(line)); + } + } + else + // otherwise call the more complicated functions and ask for inner points + // from the boundary description + { + std::vector > line_points (this->polynomial_degree-1); + // 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 + for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) + { + const typename Triangulation::line_iterator + line = (dim == 1 + ? + static_cast::line_iterator>(cell) + : + cell->line(line_no)); + + const Manifold &manifold = + ( ( line->manifold_id() == numbers::invalid_manifold_id ) && + ( dim < spacedim ) + ? + cell->get_manifold() : + line->get_manifold() ); + + get_intermediate_points_on_object (manifold, line_support_points, line, line_points); + + if (dim==3) + { + // in 3D, lines might be in wrong orientation. if so, reverse + // the vector + if (cell->line_orientation(line_no)) + a.insert (a.end(), line_points.begin(), line_points.end()); + else + a.insert (a.end(), line_points.rbegin(), line_points.rend()); + } + else + // in 2D, lines always have the correct orientation. simply append + // all points + a.insert (a.end(), line_points.begin(), line_points.end()); + } + } +} + + + +template <> +void +MappingManifold<3,3>:: +add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, + std::vector > &a) const +{ + const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, + vertices_per_face = GeometryInfo<3>::vertices_per_face, + lines_per_face = GeometryInfo<3>::lines_per_face, + vertices_per_cell = GeometryInfo<3>::vertices_per_cell; + + static const StraightBoundary<3> straight_boundary; + // used if face quad at boundary or entirely in the interior of the domain + std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + // used if only one line of face quad is at boundary + std::vector > b(4*polynomial_degree); + + // Used by the new Manifold interface. This vector collects the + // vertices used to compute the intermediate points. + std::vector > vertices(4); + + // loop over all faces and collect points on them + for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); + + // select the correct mappings for the present face + const bool face_orientation = cell->face_orientation(face_no), + face_flip = cell->face_flip (face_no), + face_rotation = cell->face_rotation (face_no); + +#ifdef DEBUG + // some sanity checks up front + for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( + GeometryInfo<3>::face_to_cell_vertices(face_no, i, + face_orientation, + face_flip, + face_rotation)), + ExcInternalError()); + + // indices of the lines that bound a face are given by GeometryInfo<3>:: + // face_to_cell_lines + for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( + face_no, i, face_orientation, face_flip, face_rotation)), + ExcInternalError()); +#endif + + // if face at boundary, then ask boundary object to return intermediate + // points on it + if (face->at_boundary()) + { + get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points); + + // in 3D, the orientation, flip and rotation of the face might not + // match what we expect here, namely the standard orientation. thus + // reorder points accordingly. since a Mapping uses the same shape + // function as an FE_Q, we can ask a FE_Q to do the reordering for us. + for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, + face_orientation, + face_flip, + face_rotation)]); + } + else + { + // face is not at boundary, but maybe some of its lines are. count + // them + unsigned int lines_at_boundary=0; + for (unsigned int i=0; iline(i)->at_boundary()) + ++lines_at_boundary; + + Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); + + // if at least one of the lines bounding this quad is at the + // boundary, then collect points separately + if (lines_at_boundary>0) + { + // call of function add_weighted_interior_points increases size of b + // about 1. There resize b for the case the mentioned function + // was already called. + b.resize(4*polynomial_degree); + + // b is of size 4*degree, make sure that this is the right size + Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1), + ExcDimensionMismatch(b.size(), + vertices_per_face+lines_per_face*(polynomial_degree-1))); + + // sort the points into b. We used access from the cell (not + // from the face) to fill b, so we can assume a standard face + // orientation. Doing so, the calculated points will be in + // standard orientation as well. + for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; + + for (unsigned int i=0; i::face_to_cell_lines( + face_no, i)*(polynomial_degree-1)+j]; + + // Now b includes the support points on the quad and we can + // apply the laplace vector + add_weighted_interior_points (support_point_weights_on_quad, b); + AssertDimension (b.size(), + 4*this->polynomial_degree + + (this->polynomial_degree-1)*(this->polynomial_degree-1)); + + for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i) + a.push_back(b[4*polynomial_degree+i]); + } + else + { + // face is entirely in the interior. get intermediate + // points from the relevant manifold object. + vertices.resize(4); + for (unsigned int i=0; i<4; ++i) + vertices[i] = face->vertex(i); + get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points); + // in 3D, the orientation, flip and rotation of the face might + // not match what we expect here, namely the standard + // orientation. thus reorder points accordingly. since a Mapping + // uses the same shape function as an FE_Q, we can ask a FE_Q to + // do the reordering for us. + for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, + face_orientation, + face_flip, + face_rotation)]); + } + } + } +} + + + +template <> +void +MappingManifold<2,3>:: +add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, + std::vector > &a) const +{ + std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + get_intermediate_points_on_object (cell->get_manifold(), line_support_points, + cell, quad_points); + for (unsigned int i=0; i +void +MappingManifold:: +add_quad_support_points(const typename Triangulation::cell_iterator &, + std::vector > &) const +{ + Assert (false, ExcInternalError()); +} + + + +template +std::vector > +MappingManifold:: +compute_mapping_support_points(const typename Triangulation::cell_iterator &cell) const +{ + // get the vertices first + std::vector > a(GeometryInfo::vertices_per_cell); + for (unsigned int i=0; i::vertices_per_cell; ++i) + a[i] = cell->vertex(i); + + if (this->polynomial_degree>1) + switch (dim) + { + case 1: + add_line_support_points(cell, a); + break; + case 2: + // in 2d, add the points on the four bounding lines to the exterior + // (outer) points + add_line_support_points(cell, a); + + // then get the support points on the quad if we are on a + // manifold, otherwise compute them from the points around it + if (dim != spacedim) + add_quad_support_points(cell, a); + else + add_weighted_interior_points (support_point_weights_on_quad, a); + break; + + case 3: + { + // in 3d also add the points located on the boundary faces + add_line_support_points (cell, a); + add_quad_support_points (cell, a); + + // then compute the interior points + add_weighted_interior_points (support_point_weights_on_hex, a); + break; + } + + default: + Assert(false, ExcNotImplemented()); + break; + } + + return a; +} + + + +//--------------------------- Explicit instantiations ----------------------- +#include "mapping_q_generic.inst" + + +DEAL_II_NAMESPACE_CLOSE