From ccfc3cacaa3d83b35c029890f349d19018c7e1f3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 3 Sep 2015 20:57:34 -0500 Subject: [PATCH] Refactor MappingQ/Q1 into a common base class. This patch separates out the common functionality for generic polynomial mappings into a base class. It contains most of what used to be in MappingQ1 but really was not specific to p=1 mappings. --- include/deal.II/fe/mapping_q.h | 7 +- include/deal.II/fe/mapping_q1.h | 509 +----- include/deal.II/fe/mapping_q_generic.h | 604 +++++++ source/fe/CMakeLists.txt | 2 + source/fe/mapping_q1.cc | 2079 +--------------------- source/fe/mapping_q_generic.cc | 2259 ++++++++++++++++++++++++ source/fe/mapping_q_generic.inst.in | 23 + 7 files changed, 2901 insertions(+), 2582 deletions(-) create mode 100644 include/deal.II/fe/mapping_q_generic.h create mode 100644 source/fe/mapping_q_generic.cc create mode 100644 source/fe/mapping_q_generic.inst.in diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 470352e90d..7fb3cc0364 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -56,8 +57,8 @@ class MappingQ : public MappingQ1 { public: /** - * Constructor. @p p gives the degree of mapping polynomials on boundary - * cells. + * Constructor. @p polynomial_degree denotes the polynomial degree + * of the polynomials that are used to map cells boundary. * * The second argument determines whether the higher order mapping should * also be used on interior cells. If its value is false (the @@ -68,7 +69,7 @@ public: * would also like to use a higher order mapping in the interior. The * MappingQEulerian class is one such case. */ - MappingQ (const unsigned int p, + MappingQ (const unsigned int polynomial_degree, const bool use_mapping_q_on_all_cells = false); /** diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h index bf07fa14b9..bbdf4b4df8 100644 --- a/include/deal.II/fe/mapping_q1.h +++ b/include/deal.II/fe/mapping_q1.h @@ -23,7 +23,7 @@ #include #include #include -#include +#include #include @@ -37,7 +37,8 @@ DEAL_II_NAMESPACE_OPEN * Mapping of the reference to cell to a general * quadrilateral/hexahedra by $d$-linear shape functions. * - * This function maps the reference (unit) cell to a general grid cell with + * This mapping implemented by this class maps the reference (unit) cell + * to a general grid cell with * straight lines in $d$ dimensions. (Note, however, that in 3D the * faces of a general, trilinearly mapped cell may be curved, even if the * edges are not). This is the standard mapping used for polyhedral domains. It @@ -51,7 +52,7 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2000, 2001; Ralf Hartmann, 2000, 2001, 2005 */ template -class MappingQ1 : public Mapping +class MappingQ1 : public MappingQGeneric { public: /** @@ -63,12 +64,6 @@ public: virtual Mapping *clone () const; - /** - * Always returns @p true because MappingQ1 preserves vertex locations. - */ - virtual - bool preserves_vertex_locations () const; - /** * @name Mapping points between reference and real cells * @{ @@ -90,54 +85,6 @@ public: * @} */ - /** - * @name Functions to transform tensors from reference to real coordinates - * @{ - */ - - // for documentation, see the Mapping base class - virtual - void - transform (const VectorSlice > > input, - const MappingType type, - const typename Mapping::InternalDataBase &internal, - VectorSlice > > output) const; - - // for documentation, see the Mapping base class - virtual - void - transform (const VectorSlice > > input, - const MappingType type, - const typename Mapping::InternalDataBase &internal, - VectorSlice > > output) const; - - // for documentation, see the Mapping base class - virtual - void - transform (const VectorSlice > > input, - const MappingType type, - const typename Mapping::InternalDataBase &internal, - VectorSlice > > output) const; - - // for documentation, see the Mapping base class - virtual - void - transform (const VectorSlice > > input, - const MappingType type, - const typename Mapping::InternalDataBase &internal, - VectorSlice > > output) const; - - // for documentation, see the Mapping base class - virtual - void - transform (const VectorSlice > > input, - const MappingType type, - const typename Mapping::InternalDataBase &internal, - VectorSlice > > output) const; - - /** - * @} - */ /** * @name Interface with FEValues @@ -146,298 +93,13 @@ public: public: /** - * Storage for internal data of d-linear mappings. See Mapping::InternalDataBase - * for an extensive description. - * - * For the current class MappingQ1, the InternalData class store - * 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'. + * Use the InternalData class of the base class without modification + * and additions. */ - 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; - }; + typedef typename MappingQGeneric::InternalData InternalData; protected: - // 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 - typename Mapping::InternalDataBase * - get_face_data (const UpdateFlags flags, - const Quadrature& quadrature) const; - - // documentation can be found in Mapping::get_subface_data() - virtual - typename Mapping::InternalDataBase * - 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; /** * @} @@ -480,15 +142,6 @@ protected: const Point &p) const; - /** - * Declare a convenience typedef for the class that describes offsets into - * quadrature formulas projected onto faces and subfaces. - */ - typedef - typename QProjector::DataSetDescriptor - DataSetDescriptor; - - /** * Transforms a point @p p on the unit cell to the point @p p_real on the * real cell @p cell and returns @p p_real. @@ -587,154 +240,6 @@ struct StaticMappingQ1 /*@}*/ -/*----------------------------------------------------------------------*/ - -#ifndef DOXYGEN - -template -inline -const double & -MappingQ1::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 & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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> & -MappingQ1::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 -MappingQ1::preserves_vertex_locations () const -{ - return true; -} - -#endif // DOXYGEN - -/* -------------- declaration of explicit specializations ------------- */ - DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h new file mode 100644 index 0000000000..5a1e97b2da --- /dev/null +++ b/include/deal.II/fe/mapping_q_generic.h @@ -0,0 +1,604 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__mapping_q_generic_h +#define dealii__mapping_q_generic_h + + +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/*!@addtogroup mapping */ +/*@{*/ + + +/** + * A base class for all polynomial mappings. In particular, this class + * provides the basis for the MappingQ1 and MappingQ classes that + * implement (bi-, tri-)linear mappings and higher order mappings, + * respectively. + * + * + *

Implementation

+ * + * This class provides essentially the entire generic infrastructure + * for polynomial mappings. What it requires to work from derived + * classes is an implementation of the + * compute_mapping_support_points() class that provides a list of + * locations to which the support points of the mapping (e.g., the + * vertices of the cell in the lowest order case) should be mapped. + * + * + * @author Wolfgang Bangerth, 2015 + */ +template +class MappingQGeneric : 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. + */ + MappingQGeneric (const unsigned int polynomial_degree); + + /** + * Always returns @p true because MappingQ1 preserves vertex locations. + */ + virtual + bool preserves_vertex_locations () const; + + + /** + * @name Functions to transform tensors from reference to real coordinates + * @{ + */ + + // for documentation, see the Mapping base class + virtual + void + transform (const VectorSlice > > input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + VectorSlice > > output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const VectorSlice > > input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + VectorSlice > > output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const VectorSlice > > input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + VectorSlice > > output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const VectorSlice > > input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + VectorSlice > > output) const; + + // for documentation, see the Mapping base class + virtual + void + transform (const VectorSlice > > input, + const MappingType type, + const typename Mapping::InternalDataBase &internal, + VectorSlice > > 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; + }; + +protected: + + // 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 + typename Mapping::InternalDataBase * + get_face_data (const UpdateFlags flags, + const Quadrature& quadrature) const; + + // documentation can be found in Mapping::get_subface_data() + virtual + typename Mapping::InternalDataBase * + 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; + + /** + * @} + */ + + + /** + * The degree of the polynomials used as shape functions for the mapping + * of cells. + */ + const unsigned int polynomial_degree; + + /** + * An interface that derived classes have to implement and that + * computes the locations of support points for the mapping. For + * example, for MappingQ1 these are the vertices. 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. + */ + virtual + void + compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, + std::vector > &a) const = 0; +}; + + + +/*@}*/ + +/*----------------------------------------------------------------------*/ + +#ifndef DOXYGEN + +template +inline +const double & +MappingQGeneric::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 & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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> & +MappingQGeneric::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 +MappingQGeneric::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 228cd4e686..d00504a3c6 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -52,6 +52,7 @@ SET(_src mapping_cartesian.cc mapping.cc mapping_fe_field.cc + mapping_q_generic.cc mapping_q1.cc mapping_q1_eulerian.cc mapping_q.cc @@ -95,6 +96,7 @@ SET(_inst mapping_cartesian.inst.in mapping.inst.in mapping_fe_field.inst.in + mapping_q_generic.inst.in mapping_q1_eulerian.inst.in mapping_q1.inst.in mapping_q_eulerian.inst.in diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index d4f83839d3..2b522e6d0f 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -38,165 +38,12 @@ DEAL_II_NAMESPACE_OPEN -template -MappingQ1::InternalData::InternalData (const unsigned int polynomial_degree) - : - polynomial_degree (polynomial_degree), - n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) -{} - - - -template -std::size_t -MappingQ1::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 -MappingQ1::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 -MappingQ1::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); - } - } - } - } -} - - template MappingQ1::MappingQ1 () + : + MappingQGeneric (1) {} @@ -322,1030 +169,6 @@ namespace internal Assert(false, ExcInternalError()); return Point<3>(); } - - - - template - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingQ1<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::MappingQ1<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::MappingQ1<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; - } - } - } - } -} - - -namespace -{ - template - std::vector - get_dpo_vector (const unsigned int degree) - { - std::vector dpo(dim+1, 1U); - for (unsigned int i=1; i -void -MappingQ1::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; point -UpdateFlags -MappingQ1::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 MappingQ1::InternalData * -MappingQ1::get_data (const UpdateFlags update_flags, - const Quadrature &q) const -{ - InternalData *data = new InternalData(1); - data->initialize (requires_update_flags(update_flags), q, q.size()); - - return data; -} - - - -template -typename Mapping::InternalDataBase * -MappingQ1::get_face_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const -{ - InternalData *data = new InternalData(1); - data->initialize_face (requires_update_flags(update_flags), - QProjector::project_to_all_faces(quadrature), - quadrature.size()); - - return data; -} - - - -template -typename Mapping::InternalDataBase * -MappingQ1::get_subface_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const -{ - InternalData *data = new InternalData(1); - data->initialize_face (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::MappingQ1::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::MappingQ1::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::MappingQ1::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::MappingQ1::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) - { - 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::MappingQ1::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::MappingQ1::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) - { - 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::MappingQ1::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::MappingQ1::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) - { - 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::compute_mapping_support_points( -template -CellSimilarity::Similarity -MappingQ1:: -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)) - { - compute_mapping_support_points(cell, data.mapping_support_points); - data.cell_of_current_support_points = cell; - } - - internal::maybe_compute_q_points (DataSetDescriptor::cell (), - data, - output_data.quadrature_points); - internal::maybe_update_Jacobians (cell_similarity, - 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, - DataSetDescriptor::cell (), - data, - output_data.jacobian_grads); - - internal::maybe_update_jacobian_pushed_forward_grads (cell_similarity, - DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_grads); - - internal::maybe_update_jacobian_2nd_derivatives (cell_similarity, - DataSetDescriptor::cell (), - data, - output_data.jacobian_2nd_derivatives); - - internal::maybe_update_jacobian_pushed_forward_2nd_derivatives (cell_similarity, - DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_2nd_derivatives); - - internal::maybe_update_jacobian_3rd_derivatives (cell_similarity, - DataSetDescriptor::cell (), - data, - output_data.jacobian_3rd_derivatives); - - internal::maybe_update_jacobian_pushed_forward_3rd_derivatives (cell_similarity, - 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::MappingQ1 &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::MappingQ1::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 (data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d], - mapping_contravariant, - data, - 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 cell_normal; - const DerivativeForm<1,spacedim,dim> DX_t = - data.contravariant[point].transpose(); - cross_product(cell_normal,DX_t[0],DX_t[1]); - cell_normal /= cell_normal.norm(); - - // then compute the face normal from the face tangent - // and the cell normal: - cross_product (output_data.boundary_forms[point], - 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::MappingQ1 &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::MappingQ1::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 -MappingQ1:: -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)) - { - compute_mapping_support_points(cell, data.mapping_support_points); - data.cell_of_current_support_points = cell; - } - - internal::do_fill_fe_face_values (*this, - cell, face_no, numbers::invalid_unsigned_int, - 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 -MappingQ1:: -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)) - { - compute_mapping_support_points(cell, data.mapping_support_points); - data.cell_of_current_support_points = cell; - } - - internal::do_fill_fe_face_values (*this, - cell, face_no, subface_no, - 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 VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - ExcInternalError()); - const typename MappingQ1::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 VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - ExcInternalError()); - const typename MappingQ1::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 VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - ExcInternalError()); - const typename MappingQ1::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 VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) - { - AssertDimension (input.size(), output.size()); - Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - ExcInternalError()); - const typename MappingQ1::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 -MappingQ1:: -transform (const VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) const -{ - transform_fields(input, mapping_type, mapping_data, output); -} - - - -template -void -MappingQ1:: -transform (const VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > output) const -{ - transform_differential_forms(input, mapping_type, mapping_data, output); -} - - - -template -void -MappingQ1:: -transform (const VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > 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 -MappingQ1:: -transform (const VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > 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 -MappingQ1:: -transform (const VectorSlice > > input, - const MappingType mapping_type, - const typename Mapping::InternalDataBase &mapping_data, - VectorSlice > > 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()); - } - -} - - template Point diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc new file mode 100644 index 0000000000..3d8df89ddf --- /dev/null +++ b/source/fe/mapping_q_generic.cc @@ -0,0 +1,2259 @@ +// --------------------------------------------------------------------- +// +// 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 + + +DEAL_II_NAMESPACE_OPEN + + +template +MappingQGeneric::InternalData::InternalData (const unsigned int polynomial_degree) + : + polynomial_degree (polynomial_degree), + n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) +{} + + + +template +std::size_t +MappingQGeneric::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 +MappingQGeneric::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 +MappingQGeneric::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 internal +{ + namespace MappingQGeneric + { + // 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; + // fast exit if the point is not in the cell (this is the only case + // where the discriminant is negative) + if (discriminant < 0.0) + { + return Point<2>(2, 2); + } + + 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: + else if (std::abs(c/b) < 1e-12) + { + 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 + { + return Point<2>(2, 2); + } + } + } + + + + 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>(); + } + + + + template + void + compute_shape_function_values (const unsigned int n_shape_functions, + const std::vector > &unit_points, + typename dealii::MappingQGeneric<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::MappingQGeneric<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::MappingQGeneric<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; + } + } + } + } +} + + +namespace +{ + template + std::vector + get_dpo_vector (const unsigned int degree) + { + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +void +MappingQGeneric::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::MappingQGeneric::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; point +MappingQGeneric::MappingQGeneric (const unsigned int p) + : + polynomial_degree(p) +{} + + + +template +UpdateFlags +MappingQGeneric::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 MappingQGeneric::InternalData * +MappingQGeneric::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 Mapping::InternalDataBase * +MappingQGeneric::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 Mapping::InternalDataBase * +MappingQGeneric::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::MappingQGeneric::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::MappingQGeneric::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::MappingQGeneric::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::MappingQGeneric::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) + { + 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::MappingQGeneric::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::MappingQGeneric::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) + { + 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::MappingQGeneric::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::MappingQGeneric::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) + { + 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 +MappingQGeneric:: +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)) + { + compute_mapping_support_points(cell, data.mapping_support_points); + 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::MappingQGeneric &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::MappingQGeneric::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 (data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d], + mapping_contravariant, + data, + 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 cell_normal; + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); + cross_product(cell_normal,DX_t[0],DX_t[1]); + cell_normal /= cell_normal.norm(); + + // then compute the face normal from the face tangent + // and the cell normal: + cross_product (output_data.boundary_forms[point], + 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::MappingQGeneric &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::MappingQGeneric::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 +MappingQGeneric:: +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)) + { + compute_mapping_support_points(cell, data.mapping_support_points); + 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 +MappingQGeneric:: +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)) + { + compute_mapping_support_points(cell, data.mapping_support_points); + 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 VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingQGeneric::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 VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingQGeneric::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 VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingQGeneric::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 VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingQGeneric::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 +MappingQGeneric:: +transform (const VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) const +{ + transform_fields(input, mapping_type, mapping_data, output); +} + + + +template +void +MappingQGeneric:: +transform (const VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > output) const +{ + transform_differential_forms(input, mapping_type, mapping_data, output); +} + + + +template +void +MappingQGeneric:: +transform (const VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > 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 +MappingQGeneric:: +transform (const VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > 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 +MappingQGeneric:: +transform (const VectorSlice > > input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + VectorSlice > > 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()); + } + +} + + + +//--------------------------- Explicit instantiations ----------------------- +#include "mapping_q_generic.inst" + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/mapping_q_generic.inst.in b/source/fe/mapping_q_generic.inst.in new file mode 100644 index 0000000000..06156a2cce --- /dev/null +++ b/source/fe/mapping_q_generic.inst.in @@ -0,0 +1,23 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class MappingQGeneric; +#endif + } + -- 2.39.5