From d37432a91ea6ca296094f808c93cadd47bd817be Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 3 Mar 2016 19:03:35 +0100 Subject: [PATCH] First shot at Jacobians --- include/deal.II/fe/mapping_manifold.h | 48 +++--- source/fe/mapping_manifold.cc | 236 ++++++++++++-------------- 2 files changed, 132 insertions(+), 152 deletions(-) diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index ccf5ea2c6d..52371c38e0 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -44,9 +44,9 @@ template class MappingQ; * reference and real cell by exploiting the geometrical information * coming from the underlying Manifold object. * - * Quadrature points computed using this mapping lye on the exact + * Quadrature points computed using this mapping lie on the exact * geometrical objects, and tangent and normal vectors computed using - * this class are normal and tangent to the underlying geometry. This + * this class are tangent and normal to the underlying geometry. This * is in constrast with the MappingQ class, which approximates the * geometry using a polynomial of some order, and then computes the * normals and tangents using the approximated surface. @@ -217,6 +217,14 @@ public: */ typename Triangulation::cell_iterator current_cell; + /** + * The actual quadrature on the reference cell. + * + * Computed once. + */ + Quadrature quad; + + /** * Values of quadrature weights for manifold quadrature formulas. * @@ -287,25 +295,25 @@ public: // */ // 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 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; + /** + * 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. diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 2832ffd07b..d6216c6779 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -51,8 +51,6 @@ std::size_t MappingManifold::InternalData::memory_consumption () const { return (Mapping::InternalDataBase::memory_consumption() ); - // MemoryConsumption::memory_consumption (shape_values) + - // MemoryConsumption::memory_consumption (shape_derivatives) + // MemoryConsumption::memory_consumption (covariant) + // MemoryConsumption::memory_consumption (contravariant) + // MemoryConsumption::memory_consumption (unit_tangentials) + @@ -76,33 +74,21 @@ initialize (const UpdateFlags update_flags, // in fill_fe_*_values() this->update_each = update_flags; + // Store the quadrature + this->quad = q; 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) compute_manifold_quadrature_weights(q); - // 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_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); @@ -405,78 +391,64 @@ namespace internal * Update the co- and contravariant matrices as well as their determinant, for the cell * described stored in the data object, but only if the update_flags of the @p data * argument indicate so. - * - * Skip the computation if possible as indicated by the first argument. */ template void - maybe_update_Jacobians (const CellSimilarity::Similarity cell_similarity, + maybe_update_Jacobians (const typename dealii::Triangulation::cell_iterator &cell, const typename dealii::QProjector::DataSetDescriptor data_set, const typename dealii::MappingManifold::InternalData &data) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_contravariant_transformation) - // if the current cell is just a - // translation of the previous one, no - // need to recompute jacobians... - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); + { + const unsigned int n_q_points = data.contravariant.size(); - std::fill(data.contravariant.begin(), data.contravariant.end(), - DerivativeForm<1,dim,spacedim>()); + 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]; + for (unsigned int point=0; point &p = data.quad.point(point+data_set); + + // Always get the maximum length from the point to the + // boundary of the reference element, to compute the + // tangent vectors from the Manifold object + for (unsigned int i=0; i ei = Point::unit_vector(i); + double ai = ei*p; + Assert(ai >=0, ExcInternalError("Was expecting a quadrature point " + "inside the unit reference element.")); + Point np(ai > .5 ? p-ai *ei : p+(1-ai)*ei); - // peel away part of sum to avoid zeroing the - // entries and adding for the first time - for (unsigned int i=0; i .5 ? -ai: 1-ai; - if (update_flags & update_covariant_transformation) - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); - for (unsigned int point=0; pointget_manifold().get_tangent_vector(p, np)/L; + } + } - if (update_flags & update_volume_elements) - if (cell_similarity != CellSimilarity::translation) - { - const unsigned int n_q_points = data.contravariant.size(); - for (unsigned int point=0; point::InternalDataBase &mapping_data, const ArrayView > &output) { - // AssertDimension (input.size(), output.size()); - // Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - // ExcInternalError()); - // const typename MappingManifold::InternalData - // &data = static_cast::InternalData &>(mapping_data); - - // switch (mapping_type) - // { - // case mapping_covariant: - // { - // Assert (data.update_each & update_contravariant_transformation, - // typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); - - // for (unsigned int i=0; i::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_covariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i > &input, const ArrayView > &output) const { -// AssertDimension (input.size(), output.size()); -// Assert (dynamic_cast(&mapping_data) != 0, -// ExcInternalError()); -// const InternalData &data = static_cast(mapping_data); + 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")); + 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