From: Luca Heltai Date: Sat, 5 Mar 2016 13:35:38 +0000 (+0100) Subject: MappingManifold now works on faces. X-Git-Tag: v8.5.0-rc1~1131^2~16 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8b7c0e59e2e0b83334398716320540e181fa635e;p=dealii.git MappingManifold now works on faces. --- diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index 4a16e5866f..adb2641911 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -276,20 +276,20 @@ public: // */ // 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; + /** + * 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 @@ -328,10 +328,10 @@ public: */ mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant; - // /** - // * Auxiliary vectors for internal use. - // */ - // mutable std::vector > > aux; + /** + * Auxiliary vectors for internal use. + */ + mutable std::vector > > aux; // /** // * Stores the support points of the mapping shape functions on the @p @@ -354,6 +354,13 @@ public: * A Q1 Finite element, to compute weights. */ const FE_Q fe_q; + + /** + * A pointer to the Manifold in use. + * + * Updated each. + */ + mutable SmartPointer > manifold; }; diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 3a69a9ed25..a047b3af4f 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -118,63 +119,63 @@ 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); - // } - // } - // } - // } + 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); + } + } + } + } } @@ -350,6 +351,29 @@ namespace internal { namespace { + /** + * Some specialization for face Manifolds. In one dimension, there + * are no Manifolds associated to faces. + */ + template + const dealii::Manifold<1, spacedim> & + get_face_manifold(const typename dealii::Triangulation<1,spacedim>::cell_iterator &cell, + const unsigned int &) + { + return cell->get_manifold(); + } + + /** + * Some specialization for face Manifolds. + */ + template + const dealii::Manifold & + get_face_manifold(const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no) + { + return cell->face(face_no)->get_manifold(); + } + /** * Compute the locations of quadrature points on the object described by * the first argument (and the cell for which the mapping support points @@ -373,7 +397,7 @@ namespace internal for (unsigned int point=0; pointget_manifold(). + quadrature_points[point] = data.manifold-> get_new_point(Quadrature(data.vertices, data.cell_manifold_quadrature_weights[point])); } @@ -412,7 +436,7 @@ namespace internal const Point &p = data.quad.point(point+data_set); // And get its image on the manifold: - const Point P = data.cell->get_manifold(). + const Point P = data.manifold-> get_new_point(Quadrature(data.vertices, data.cell_manifold_quadrature_weights[point+data_set])); @@ -435,10 +459,10 @@ namespace internal for (unsigned int j=0; j::vertices_per_cell; ++j) weights[j] = data.fe_q.shape_value(j, np); - Point NP=data.cell->get_manifold(). + Point NP=data.manifold-> get_new_point(Quadrature(data.vertices, weights)); - Tensor<1,spacedim> T = data.cell->get_manifold().get_tangent_vector(P, NP); + Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP); for (unsigned int d=0; d::cell_iterator &cell, const unsigned int n_q_points=quadrature.size(); data.store_vertices(cell); + data.manifold = &(cell->get_manifold()); internal::maybe_compute_q_points (QProjector::DataSetDescriptor::cell (), data, @@ -1221,7 +1246,7 @@ namespace internal */ template void - do_fill_fe_face_values (const dealii::MappingManifold &mapping, + do_fill_fe_face_values (const dealii::MappingManifold &mapping, const typename dealii::Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, @@ -1230,37 +1255,43 @@ namespace internal const typename dealii::MappingManifold::InternalData &data, internal::FEValues::MappingRelatedData &output_data) { - maybe_compute_q_points (cell, - data_set, + data.store_vertices(cell); + + // This should really be get_face_manifold(cell, face_no), but + // that does not compile... At the moment this class works only + // on the cell manifold, and does not respect face manifolds. + data.manifold = &cell->get_manifold(); + + maybe_compute_q_points (data_set, data, output_data.quadrature_points); - maybe_update_Jacobians (CellSimilarity::none, - data_set, + 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_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(), @@ -1281,37 +1312,22 @@ fill_fe_face_values (const typename Triangulation::cell_iterator & const typename Mapping::InternalDataBase &internal_data, internal::FEValues::MappingRelatedData &output_data) const { -// // ensure that the following cast is really correct: -// Assert ((dynamic_cast(&internal_data) != 0), -// ExcInternalError()); -// const InternalData &data -// = static_cast(internal_data); - -// // if necessary, recompute the support points of the transformation of this cell -// // (note that we need to first check the triangulation pointer, since otherwise -// // the second test might trigger an exception if the triangulations are not the -// // same) -// if ((data.mapping_support_points.size() == 0) -// || -// (&cell->get_triangulation() != -// &data.cell_of_current_support_points->get_triangulation()) -// || -// (cell != data.cell_of_current_support_points)) -// { -// data.mapping_support_points = this->compute_mapping_support_points(cell); -// data.cell_of_current_support_points = cell; -// } - -// internal::do_fill_fe_face_values (*this, -// cell, face_no, numbers::invalid_unsigned_int, -// QProjector::DataSetDescriptor::face (face_no, -// cell->face_orientation(face_no), -// cell->face_flip(face_no), -// cell->face_rotation(face_no), -// quadrature.size()), -// quadrature, -// data, -// output_data); + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != 0), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + 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); Assert(false, ExcNotImplemented()); } @@ -1327,105 +1343,90 @@ fill_fe_subface_values (const typename Triangulation::cell_iterato const typename Mapping::InternalDataBase &internal_data, internal::FEValues::MappingRelatedData &output_data) const { -// // ensure that the following cast is really correct: -// Assert ((dynamic_cast(&internal_data) != 0), -// ExcInternalError()); -// const InternalData &data -// = static_cast(internal_data); - -// // if necessary, recompute the support points of the transformation of this cell -// // (note that we need to first check the triangulation pointer, since otherwise -// // the second test might trigger an exception if the triangulations are not the -// // same) -// if ((data.mapping_support_points.size() == 0) -// || -// (&cell->get_triangulation() != -// &data.cell_of_current_support_points->get_triangulation()) -// || -// (cell != data.cell_of_current_support_points)) -// { -// data.mapping_support_points = this->compute_mapping_support_points(cell); -// data.cell_of_current_support_points = cell; -// } - -// internal::do_fill_fe_face_values (*this, -// cell, face_no, subface_no, -// QProjector::DataSetDescriptor::subface (face_no, subface_no, -// cell->face_orientation(face_no), -// cell->face_flip(face_no), -// cell->face_rotation(face_no), -// quadrature.size(), -// cell->subface_case(face_no)), -// quadrature, -// data, -// output_data); + // ensure that the following cast is really correct: + Assert ((dynamic_cast(&internal_data) != 0), + ExcInternalError()); + const InternalData &data + = static_cast(internal_data); + + internal::do_fill_fe_face_values (*this, + cell, face_no, subface_no, + QProjector::DataSetDescriptor::subface (face_no, subface_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)), + quadrature, + data, + output_data); } namespace { - // template - // void - // transform_fields(const ArrayView > &input, - // const MappingType mapping_type, - // const typename Mapping::InternalDataBase &mapping_data, - // const ArrayView > &output) - // { - // AssertDimension (input.size(), output.size()); - // Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), - // ExcInternalError()); - // const typename MappingManifold::InternalData - // &data = static_cast::InternalData &>(mapping_data); - - // switch (mapping_type) - // { - // case mapping_contravariant: - // { - // Assert (data.update_each & update_contravariant_transformation, - // typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); - - // for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); - // Assert (data.update_each & update_volume_elements, - // typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); - // Assert (rank==1, ExcMessage("Only for rank 1")); - // if (rank!=1) - // return; - - // for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); - - // for (unsigned int i=0; i + void + transform_fields(const ArrayView > &input, + const MappingType mapping_type, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView > &output) + { + AssertDimension (input.size(), output.size()); + Assert ((dynamic_cast::InternalData *>(&mapping_data) != 0), + ExcInternalError()); + const typename MappingManifold::InternalData + &data = static_cast::InternalData &>(mapping_data); + + switch (mapping_type) + { + case mapping_contravariant: + { + Assert (data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + if (rank!=1) + return; + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_covariant_transformation")); + + for (unsigned int i=0; i @@ -1718,20 +1719,20 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { - // switch (mapping_type) - // { - // case mapping_contravariant: - // transform_fields(input, mapping_type, mapping_data, output); - // return; - - // case mapping_piola_gradient: - // case mapping_contravariant_gradient: - // case mapping_covariant_gradient: - // transform_gradients(input, mapping_type, mapping_data, output); - // return; - // default: - Assert(false, ExcNotImplemented()); - // } + 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()); + } } @@ -1793,16 +1794,16 @@ transform (const ArrayView > &input, const typename Mapping::InternalDataBase &mapping_data, const ArrayView > &output) const { -// switch (mapping_type) -// { -// case mapping_piola_hessian: -// case mapping_contravariant_hessian: -// case mapping_covariant_hessian: -// transform_hessians(input, mapping_type, mapping_data, output); -// return; -// default: - Assert(false, ExcNotImplemented()); -// } + 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()); + } } @@ -1971,127 +1972,6 @@ namespace } } - -// template -// void -// MappingManifold:: -// add_line_support_points (const typename Triangulation::cell_iterator &cell, -// std::vector > &a) const -// { -// // if we only need the midpoint, then ask for it. -// if (this->polynomial_degree==2) -// { -// for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) -// { -// const typename Triangulation::line_iterator line = -// (dim == 1 ? -// static_cast::line_iterator>(cell) : -// cell->line(line_no)); - -// const Manifold &manifold = -// ( ( line->manifold_id() == numbers::invalid_manifold_id ) && -// ( dim < spacedim ) -// ? -// cell->get_manifold() -// : -// line->get_manifold() ); -// a.push_back(manifold.get_new_point_on_line(line)); -// } -// } -// else -// // otherwise call the more complicated functions and ask for inner points -// // from the boundary description -// { -// std::vector > line_points (this->polynomial_degree-1); -// // loop over each of the lines, and if it is at the boundary, then first -// // get the boundary description and second compute the points on it -// for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) -// { -// const typename Triangulation::line_iterator -// line = (dim == 1 -// ? -// static_cast::line_iterator>(cell) -// : -// cell->line(line_no)); - -// const Manifold &manifold = -// ( ( line->manifold_id() == numbers::invalid_manifold_id ) && -// ( dim < spacedim ) -// ? -// cell->get_manifold() : -// line->get_manifold() ); - -// // get_intermediate_points_on_object (manifold, line_support_points, line, line_points); - -// if (dim==3) -// { -// // in 3D, lines might be in wrong orientation. if so, reverse -// // the vector -// if (cell->line_orientation(line_no)) -// a.insert (a.end(), line_points.begin(), line_points.end()); -// else -// a.insert (a.end(), line_points.rbegin(), line_points.rend()); -// } -// else -// // in 2D, lines always have the correct orientation. simply append -// // all points -// a.insert (a.end(), line_points.begin(), line_points.end()); -// } -// } -// } - - - -// template -// std::vector > -// MappingManifold:: -// compute_mapping_support_points(const typename Triangulation::cell_iterator &cell) const -// { -// // get the vertices first -// std::vector > a(GeometryInfo::vertices_per_cell); -// for (unsigned int i=0; i::vertices_per_cell; ++i) -// a[i] = cell->vertex(i); - -// if (this->polynomial_degree>1) -// switch (dim) -// { -// case 1: -// add_line_support_points(cell, a); -// break; -// case 2: -// // in 2d, add the points on the four bounding lines to the exterior -// // (outer) points -// add_line_support_points(cell, a); - -// // then get the support points on the quad if we are on a -// // manifold, otherwise compute them from the points around it -// if (dim != spacedim) -// add_quad_support_points(cell, a); -// else -// add_weighted_interior_points (support_point_weights_on_quad, a); -// break; - -// case 3: -// { -// // in 3d also add the points located on the boundary faces -// add_line_support_points (cell, a); -// add_quad_support_points (cell, a); - -// // then compute the interior points -// add_weighted_interior_points (support_point_weights_on_hex, a); -// break; -// } - -// default: -// Assert(false, ExcNotImplemented()); -// break; -// } - -// return a; -// } - - - //--------------------------- Explicit instantiations ----------------------- #include "mapping_manifold.inst"