From: bangerth Date: Wed, 26 Nov 2008 01:14:45 +0000 (+0000) Subject: Reimplement parts of the VectorTools::interpolate_boundary_values (in the process... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=656816238f0f818ae550f07dc4049736a2b5db68;p=dealii-svn.git Reimplement parts of the VectorTools::interpolate_boundary_values (in the process taking care of a TODO) to make sure it doesn't use insane amounts of compute time by re-computing the same information over and over again. git-svn-id: https://svn.dealii.org/trunk@17745 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index fa07a7e949..070afd6006 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -1352,6 +1352,88 @@ interpolate_boundary_values (const Mapping &mapping, dof_values_scalar.reserve (DoFTools::max_dofs_per_face (dof)); dof_values_system.reserve (DoFTools::max_dofs_per_face (dof)); + // before we start with the loop + // over all cells create an + // hp::FEValues object that holds + // the interpolation points of all + // finite elements that may ever be + // in use + hp::FECollection finite_elements (dof.get_fe()); + hp::QCollection q_collection; + for (unsigned int f=0; f &fe = finite_elements[f]; + + // generate a quadrature rule + // on the face from the unit + // support points. this will be + // used to obtain the + // quadrature points on the + // real cell's face + // + // to do this, we check whether + // the FE has support points on + // the face at all: + if (fe.has_face_support_points()) + q_collection.push_back (Quadrature(fe.get_unit_face_support_points())); + else + { + // if not, then we should + // try a more clever + // way. the idea is that a + // finite element may not + // offer support points for + // all its shape functions, + // but maybe only some. if + // it offers support points + // for the components we + // are interested in in + // this function, then + // that's fine. if not, the + // function we call in the + // finite element will + // raise an exception. the + // support points for the + // other shape functions + // are left uninitialized + // (well, initialized by + // the default + // constructor), since we + // don't need them anyway. + // + // As a detour, we must + // make sure we only query + // face_system_to_component_index + // if the index corresponds + // to a primitive shape + // function. since we know + // that all the components + // we are interested in are + // primitive (by the above + // check), we can safely + // put such a check in + // front + std::vector > unit_support_points (fe.dofs_per_face); + + for (unsigned int i=0; i(unit_support_points)); + } + } + // now that we have a q_collection + // object with all the right + // quadrature points, create an + // hp::FEFaceValues object that we + // can use to evaluate the boundary + // values at + hp::MappingCollection mapping_collection (mapping); + hp::FEFaceValues x_fe_values (mapping_collection, finite_elements, q_collection, + update_quadrature_points); + typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); for (; cell!=endc; ++cell) @@ -1387,73 +1469,16 @@ interpolate_boundary_values (const Mapping &mapping, typename DH::face_iterator face = cell->face(face_no); const unsigned char boundary_component = face->boundary_indicator(); if (function_map.find(boundary_component) != function_map.end()) - // face is of the right component { -//TODO[?] Should work for both DoFHandlers. But probably not the most efficient -// implementation. - // next generate a quadrature rule - // on the face from the unit - // support points. this wil be used - // to obtain the quadrature points - // on the real cell's face - std::vector > - unit_support_points = fe.get_unit_face_support_points(); - - // check whether there are support - // points on the face. if not, then - // we should try a more clever - // way. the idea is that a finite - // element may not offer support - // points for all its shape - // functions, but maybe only - // some. if it offers support - // points for the components we are - // interested in in this function, - // then that's fine. if not, the - // function we call in the finite - // element will raise an - // exception. the support points - // for the other shape functions - // are left uninitialized (well, - // initialized by the default - // constructor), since we don't - // need them anyway. - // - // As a detour, we must - // make sure we only - // query - // face_system_to_component_index - // if the index - // corresponds to a - // primitive shape - // function. since we - // know that all the - // components we are - // interested in are - // primitive (by the - // above check), we can - // safely put such a - // check in front - if (unit_support_points.size() == 0) - { - unit_support_points.resize (fe.dofs_per_face); - for (unsigned int i=0; i aux_quad (unit_support_points); - FEFaceValues fe_values (mapping, fe, aux_quad, update_quadrature_points); -//TODO[?] End of inefficient code + // face is of the right component + x_fe_values.reinit(cell, face_no); + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values(); // get indices, physical location and // boundary values of dofs on this // face face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); - fe_values.reinit(cell, face_no); const std::vector > &dof_locations = fe_values.get_quadrature_points (); diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index bba9ca3ffd..bda53ba63e 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -484,6 +484,15 @@ inconvenience this causes.

deal.II

    +
  1. +

    + Fixed: The VectorTools::interpolate_boundary_values function was implemented a bit + clumsily and was using much more time than necessary. This should be fixed now. + associated vertex. +
    + (WB 2008/11/25) +

    +
  2. Fixed: The GridIn::read_msh function had a bug that made it reject