From efa8ef37d91f4cd26303deadb7fe9e6d4d6283f7 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 4 Jul 2012 08:29:19 +0000 Subject: [PATCH] Handle the fact that Mapping::transform_real_to_unit_cell can now also throw exceptions. git-svn-id: https://svn.dealii.org/trunk@25672 0785d39b-7218-0410-832d-ea1e28bc413d --- .../numerics/fe_field_function.templates.h | 111 ++++++++++++------ 1 file changed, 75 insertions(+), 36 deletions(-) diff --git a/deal.II/include/deal.II/numerics/fe_field_function.templates.h b/deal.II/include/deal.II/numerics/fe_field_function.templates.h index c343871594..5af4c2960a 100644 --- a/deal.II/include/deal.II/numerics/fe_field_function.templates.h +++ b/deal.II/include/deal.II/numerics/fe_field_function.templates.h @@ -215,7 +215,7 @@ namespace Functions unsigned int nq = qpoints[i].size(); // Construct a quadrature formula std::vector< double > ww(nq, 1./((double) nq)); - + quadrature_collection.push_back (Quadrature (qpoints[i], ww)); } // Get a function value object @@ -274,7 +274,7 @@ namespace Functions unsigned int nq = qpoints[i].size(); // Construct a quadrature formula std::vector< double > ww(nq, 1./((double) nq)); - + quadrature_collection.push_back (Quadrature (qpoints[i], ww)); } // Get a function value object @@ -333,7 +333,7 @@ namespace Functions unsigned int nq = qpoints[i].size(); // Construct a quadrature formula std::vector< double > ww(nq, 1./((double) nq)); - + quadrature_collection.push_back (Quadrature (qpoints[i], ww)); } // Get a function value object @@ -397,25 +397,51 @@ namespace Functions // Current quadrature point typename DH::active_cell_iterator cell = cell_hint.get(); if (cell == dh->end()) - cell = dh->begin_active(); - Point qp = mapping.transform_real_to_unit_cell(cell, points[0]); + cell = dh->begin_active(); - // Check if we already have a - // valid cell for the first point - if (!GeometryInfo::is_inside_unit_cell(qp)) - { - const std::pair > - my_pair = GridTools::find_active_cell_around_point - (mapping, *dh, points[0]); - cell = my_pair.first; - qp = my_pair.second; - point_flags[0] = true; - } + { + // see if the point is + // inside the + // cell. there are two + // ways that + // transform_real_to_unit_cell + // can indicate that a + // point is outside: by + // returning + // coordinates that lie + // outside the + // reference cell, or + // by throwing an + // exception. handle + // both + bool point_is_inside; + Point qp; + try + { + qp = mapping.transform_real_to_unit_cell(cell, points[0]); + point_is_inside = GeometryInfo::is_inside_unit_cell(qp); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + point_is_inside = false; + } + + if (point_is_inside) + { + const std::pair > + my_pair = GridTools::find_active_cell_around_point + (mapping, *dh, points[0]); + cell = my_pair.first; + qp = my_pair.second; + point_flags[0] = true; + } + + // Put in the first point. + cells.push_back(cell); + qpoints.push_back(std::vector >(1, qp)); + maps.push_back(std::vector (1, 0)); + } - // Put in the first point. - cells.push_back(cell); - qpoints.push_back(std::vector >(1, qp)); - maps.push_back(std::vector (1, 0)); // Check if we need to do anything else if (points.size() > 1) @@ -439,22 +465,35 @@ namespace Functions // If we found one in this cell, keep looking in the same cell for (unsigned int p=first_outside; p qpoint = mapping.transform_real_to_unit_cell(cells[c], points[p]); - if (GeometryInfo::is_inside_unit_cell(qpoint)) - { - point_flags[p] = true; - qpoints[c].push_back(qpoint); - maps[c].push_back(p); - } - else - { - // Set things up for next round - if (left_over == false) - first_outside = p; - left_over = true; - } - } + if (point_flags[p] == false) + { + // same logic as above + bool point_is_inside; + Point qpoint; + try + { + qpoint = mapping.transform_real_to_unit_cell(cells[c], points[p]); + point_is_inside = GeometryInfo::is_inside_unit_cell(qpoint); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + point_is_inside = false; + } + + if (point_is_inside) + { + point_flags[p] = true; + qpoints[c].push_back(qpoint); + maps[c].push_back(p); + } + else + { + // Set things up for next round + if (left_over == false) + first_outside = p; + left_over = true; + } + } // If we got here and there is // no left over, we are // done. Else we need to find -- 2.39.5