From b744529db97e62ed40a95e76dc106679a985f3f4 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 1 Oct 2015 10:16:51 -0500 Subject: [PATCH] Move the remaining pieces of transform_r_to_u() from MappingQ to MappingQGeneric. --- include/deal.II/fe/mapping.h | 2 +- include/deal.II/fe/mapping_q1_eulerian.h | 3 +- include/deal.II/fe/mapping_q_eulerian.h | 10 +++++ source/fe/mapping_q.cc | 50 +--------------------- source/fe/mapping_q_eulerian.cc | 21 +++++++++- source/fe/mapping_q_generic.cc | 53 ++++++++++++++++++++---- 6 files changed, 78 insertions(+), 61 deletions(-) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index be9671888c..3a99fbc542 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -369,7 +369,7 @@ public: * equals the given point @p p. If this is the case then this function * throws an exception of type Mapping::ExcTransformationFailed . Whether * the given point @p p lies outside the cell can therefore be determined by - * checking whether the return reference coordinates lie inside or outside + * checking whether the returned reference coordinates lie inside or outside * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell()) or * whether the exception mentioned above has been thrown. * diff --git a/include/deal.II/fe/mapping_q1_eulerian.h b/include/deal.II/fe/mapping_q1_eulerian.h index bd2e3078a5..fc17a1734a 100644 --- a/include/deal.II/fe/mapping_q1_eulerian.h +++ b/include/deal.II/fe/mapping_q1_eulerian.h @@ -108,8 +108,7 @@ public: */ virtual std_cxx11::array, GeometryInfo::vertices_per_cell> - get_vertices - (const typename Triangulation::cell_iterator &cell) const; + get_vertices (const typename Triangulation::cell_iterator &cell) const; /** * Return a pointer to a copy of the present object. The caller of this copy diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index c81e62a688..f4bd91316b 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -116,6 +116,16 @@ public: const VECTOR &euler_vector, const DoFHandler &euler_dof_handler) DEAL_II_DEPRECATED; + /** + * Return the mapped vertices of the cell. For the current class, this function does + * not use the support points from the geometry of the current cell but + * instead evaluates an externally given displacement field in addition to + * the geometry of the cell. + */ + virtual + std_cxx11::array, GeometryInfo::vertices_per_cell> + get_vertices (const typename Triangulation::cell_iterator &cell) const; + /** * Return a pointer to a copy of the present object. The caller of this copy * then assumes ownership of it. diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 400d2998d6..e2162e64f6 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -445,60 +445,14 @@ MappingQ:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // first a Newton iteration based on a Q1 mapping to get a good starting - // point, the idea being that this is cheaper than trying to start with the - // real mapping and likely also more robust. - // - // that said, this doesn't always work: there are cases where the point is - // outside the cell and the inverse mapping doesn't converge. in that case, - // use the center point of the cell as a starting point if we are to go on - // using the full mapping, or just propagate up the exception if we had no - // intention of continuing with the full mapping - Point initial_p_unit; - try - { - initial_p_unit = q1_mapping->transform_real_to_unit_cell(cell, p); - } - catch (const typename Mapping::ExcTransformationFailed &) - { - // mirror the conditions of the code below to determine if we need to - // use an arbitrary starting point or if we just need to rethrow the - // exception - if (cell->has_boundary_lines() - || - use_mapping_q_on_all_cells - || - (dim!=spacedim) ) - { - for (unsigned int d=0; dhas_boundary_lines() || use_mapping_q_on_all_cells || (dim!=spacedim) ) - { - // use the full mapping. in case the function above should have given us - // something back that lies outside the unit cell (that might happen - // because we may have given a point 'p' that lies inside the cell with - // the higher order mapping, but outside the Q1-mapped reference cell), - // then project it back into the reference cell in hopes that this gives - // a better starting point to the following iteration - initial_p_unit = GeometryInfo::project_to_unit_cell(initial_p_unit); - - return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit); - } + return MappingQGeneric::transform_real_to_unit_cell(cell, p); else - return initial_p_unit; + return q1_mapping->transform_real_to_unit_cell(cell, p); } diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index c7a91efb09..51b8654498 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -126,13 +126,32 @@ SupportQuadrature (const unsigned int map_degree) // .... COMPUTE MAPPING SUPPORT POINTS +template +std_cxx11::array, GeometryInfo::vertices_per_cell> +MappingQEulerian:: +get_vertices +(const typename Triangulation::cell_iterator &cell) const +{ + // get the vertices as the first 2^dim mapping support points + std::vector > a; + compute_mapping_support_points(cell, a); + + std_cxx11::array, GeometryInfo::vertices_per_cell> vertex_locations; + std::copy (a.begin(), + a.begin()+GeometryInfo::vertices_per_cell, + vertex_locations.begin()); + + return vertex_locations; +} + + + template void MappingQEulerian:: compute_mapping_support_points (const typename Triangulation::cell_iterator &cell, std::vector > &a) const { - // first, basic assertion with respect to vector size, const types::global_dof_index n_dofs = euler_dof_handler->n_dofs(); diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 034f2b4edb..946a655788 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -1859,15 +1860,49 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it } - // Find the initial value for the Newton iteration by a normal - // projection to the least square plane determined by the vertices - // of the cell - std::vector > a; - compute_mapping_support_points (cell,a); - Assert(a.size() == GeometryInfo::vertices_per_cell, - ExcInternalError()); - const Point initial_p_unit = - internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + Point initial_p_unit; + if (polynomial_degree == 1) + { + // Find the initial value for the Newton iteration by a normal + // projection to the least square plane determined by the vertices + // of the cell + std::vector > a; + compute_mapping_support_points (cell,a); + Assert(a.size() == GeometryInfo::vertices_per_cell, + ExcInternalError()); + initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + } + else + { + try + { + // Find the initial value for the Newton iteration by a normal + // projection to the least square plane determined by the vertices + // of the cell + // + // we do this by first getting all support points, then + // throwing away all but the vertices, and finally calling + // the same function as above + std::vector > a; + compute_mapping_support_points (cell,a); + a.resize(GeometryInfo::vertices_per_cell); + initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + for (unsigned int d=0; d::project_to_unit_cell(initial_p_unit); + } // perform the Newton iteration and return the result. note that // this statement may throw an exception, which we simply pass up to -- 2.39.5