From c852c94c18c36dbf7f3ad4dac26c52d8d4ff4ace Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 28 Mar 2016 12:52:55 +0200 Subject: [PATCH] Fixed JxW values on faces. --- source/fe/mapping_manifold.cc | 667 +--------------------------- tests/fe/mapping_manifold_07.cc | 54 ++- tests/fe/mapping_manifold_07.output | 214 +-------- 3 files changed, 51 insertions(+), 884 deletions(-) diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 0357c2e8b4..a921dce6d5 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -96,21 +96,6 @@ initialize (const UpdateFlags update_flags, 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()); } @@ -483,440 +468,11 @@ namespace internal } } } - - /** - * Update the Hessian of the transformation from unit to real cell, the - * Jacobian gradients. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) - { - const unsigned int n_q_points = jacobian_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_grads (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_pushed_forward_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) - { - const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_2nd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_pushed_forward_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *third = - &data.third_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim][dim]; - for (unsigned int i=0; i - void - maybe_update_jacobian_pushed_forward_3rd_derivatives (const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingManifold::InternalData &data, - std::vector > &jacobian_pushed_forward_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); - - if (cell_similarity != CellSimilarity::translation) - { - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point=0; point *fourth = - &data.fourth_derivative(point+data_set, 0); - double result [spacedim][dim][dim][dim][dim]; - for (unsigned int i=0; i CellSimilarity::Similarity MappingManifold:: @@ -1049,36 +605,6 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, output_data.inverse_jacobians[point] = data.covariant[point].transpose(); } -// internal::maybe_update_jacobian_grads (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; } @@ -1262,31 +788,6 @@ namespace internal maybe_update_Jacobians (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, @@ -1686,7 +1187,7 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - // transform_fields(input, mapping_type, mapping_data, output); + transform_fields(input, mapping_type, mapping_data, output); } @@ -1799,172 +1300,6 @@ transform (const ArrayView > &input, } } - - -namespace -{ - /** - * Ask the manifold descriptor to return intermediate points on lines or - * faces. The function needs to return one or multiple points (depending on - * the number of elements in the output vector @p points that lie inside a - * line, quad or hex). Whether it is a line, quad or hex doesn't really - * matter to this function but it can be inferred from the number of input - * points in the @p surrounding_points vector. - */ - template - void - get_intermediate_points (const Manifold &manifold, - const QGaussLobatto<1> &line_support_points, - const std::vector > &surrounding_points, - std::vector > &points) - { - Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required")); - const unsigned int n=points.size(); - Assert(n>0, ExcMessage("You can't ask for 0 intermediate points.")); - std::vector w(surrounding_points.size()); - - switch (surrounding_points.size()) - { - case 2: - { - // If two points are passed, these are the two vertices, and - // we can only compute degree-1 intermediate points. - for (unsigned int i=0; i quadrature(surrounding_points, w); - points[i] = manifold.get_new_point(quadrature); - } - break; - } - - case 4: - { - Assert(spacedim >= 2, ExcImpossibleInDim(spacedim)); - const unsigned m= - static_cast(std::sqrt(static_cast(n))); - // is n a square number - Assert(m*m==n, ExcInternalError()); - - // If four points are passed, these are the two vertices, and - // we can only compute (degree-1)*(degree-1) intermediate - // points. - for (unsigned int i=0; i quadrature(surrounding_points, w); - points[i*m+j]=manifold.get_new_point(quadrature); - } - } - break; - } - - case 8: - Assert(false, ExcNotImplemented()); - break; - default: - Assert(false, ExcInternalError()); - break; - } - } - - - - - /** - * Ask the manifold descriptor to return intermediate points on the object - * pointed to by the TriaIterator @p iter. This function tries to be - * backward compatible with respect to the differences between - * Boundary and Manifold, querying the first - * whenever the passed @p manifold can be upgraded to a - * Boundary. - */ - template - void get_intermediate_points_on_object(const Manifold &manifold, - const QGaussLobatto<1> &line_support_points, - const TriaIterator &iter, - std::vector > &points) - { - const unsigned int structdim = TriaIterator::AccessorType::structure_dimension; - - // Try backward compatibility option. - if (const Boundary *boundary - = dynamic_cast *>(&manifold)) - // This is actually a boundary. Call old methods. - { - switch (structdim) - { - case 1: - { - const typename Triangulation::line_iterator line = iter; - boundary->get_intermediate_points_on_line(line, points); - return; - } - case 2: - { - const typename Triangulation::quad_iterator quad = iter; - boundary->get_intermediate_points_on_quad(quad, points); - return; - } - default: - Assert(false, ExcInternalError()); - return; - } - } - else - { - std::vector > sp(GeometryInfo::vertices_per_cell); - for (unsigned int i=0; ivertex(i); - get_intermediate_points(manifold, line_support_points, sp, points); - } - } - - - /** - * Take a support_point_weights_on_hex(quad) and apply it to the vector - * @p a to compute the inner support points as a linear combination of the - * exterior points. - * - * The vector @p a initially contains the locations of the @p n_outer - * points, the @p n_inner computed inner points are appended. - * - * See equation (7) of the `mapping' report. - */ - template - void add_weighted_interior_points(const Table<2,double> &lvs, - std::vector > &a) - { - const unsigned int n_inner_apply=lvs.n_rows(); - const unsigned int n_outer_apply=lvs.n_cols(); - Assert(a.size()==n_outer_apply, - ExcDimensionMismatch(a.size(), n_outer_apply)); - - // compute each inner point as linear combination of the outer points. the - // weights are given by the lvs entries, the outer points are the first - // (existing) elements of a - for (unsigned int unit_point=0; unit_point p; - for (unsigned int k=0; k map_manifold; - FE_Q fe(1); - const QGauss quad(5); - - FEFaceValues fe_v(map_manifold, fe, quad, - update_quadrature_points); - - out << "set size ratio -1" << std::endl - << "set terminal aqua " << dim+spacedim << std::endl - << (spacedim == 3 ? "s" : "") << "plot '-' w l " << std::endl; - - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary()) - { - fe_v.reinit(cell, f); - const std::vector > &qps = fe_v.get_quadrature_points(); - for (unsigned int i=0; i map_manifold; + FE_Q fe(1); + const QGauss quad(3); + + FEFaceValues fe_v(map_manifold, + fe, quad, + update_JxW_values); + double area = 0; + + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) { - out << qps[i] << std::endl; - if (std::abs(qps[i].distance(center) -radius) > 1e-10) - out << "# Error! This should be on the sphere, but it's not!" << std::endl; + fe_v.reinit(cell, f); + for (unsigned int i=0; i