From bc00417764101b5874b1535e43a5be1e57c9663a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 25 Jul 2016 21:33:55 -0600 Subject: [PATCH] Work around a MS Visual Studio oddity. The compiler should just choose between two overloaded template functions. But, MS Visual Studio once again gets all bent out of shape over it. Fix this by just using a plain 'if' statement, rather than relying on template specializations and overloads. This works now because we made faces of 1d cells look very much the same as faces of regular cells, and so the dim-dimensional code actually compiles also in 1d these days. --- .../deal.II/numerics/vector_tools.templates.h | 520 ++++++++---------- 1 file changed, 240 insertions(+), 280 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index ee161501e7..95c3d50b6c 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -1690,106 +1690,15 @@ namespace VectorTools namespace { - // interpolate boundary values in - // 1d. in higher dimensions, we - // use FEValues to figure out - // what to do on faces, but in 1d - // faces are points and it is far - // easier to simply work on - // individual vertices template class M_or_MC> static inline - void do_interpolate_boundary_values - (const M_or_MC &, - const DoFHandlerType &dof, - const std::map*> &function_map, - std::map &boundary_values, - const ComponentMask &component_mask, - const dealii::internal::int2type<1>) - { - const unsigned int dim = DoFHandlerType::dimension; - const unsigned int spacedim = DoFHandlerType::space_dimension; - - Assert (component_mask.represents_n_components(dof.get_fe().n_components()), - ExcMessage ("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.")); - - // if for whatever reason we were - // passed an empty map, return - // immediately - if (function_map.size() == 0) - return; - - for (typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(); - cell != dof.end(); ++cell) - for (unsigned int direction=0; - direction::faces_per_cell; ++direction) - if (cell->at_boundary(direction) - && - (function_map.find(cell->face(direction)->boundary_id()) != function_map.end())) - { - const Function &boundary_function - = *function_map.find(cell->face(direction)->boundary_id())->second; - - // get the FE corresponding to this - // cell - const FiniteElement &fe = cell->get_fe(); - Assert (fe.n_components() == boundary_function.n_components, - ExcDimensionMismatch(fe.n_components(), - boundary_function.n_components)); - - Assert (component_mask.n_selected_components(fe.n_components()) > 0, - ComponentMask::ExcNoComponentSelected()); - - // now set the value of - // the vertex degree of - // freedom. setting - // also creates the - // entry in the map if - // it did not exist - // beforehand - // - // save some time by - // requesting values - // only once for each - // point, irrespective - // of the number of - // components of the - // function - Vector function_values (fe.n_components()); - if (fe.n_components() == 1) - function_values(0) - = boundary_function.value (cell->vertex(direction)); - else - boundary_function.vector_value (cell->vertex(direction), - function_values); - - for (unsigned int i=0; i - vertex_dof_index(direction,i, - cell->active_fe_index())] - = function_values(fe.face_system_to_component_index(i).first); - } - } - - - - // template for the case dim!=1. Since the function has a template argument - // dim_, it is clearly less specialized than the 1d function above and - // whenever possible (i.e., if dim==1), the function template above - // will be used - template class M_or_MC, int dim_> - static inline void do_interpolate_boundary_values (const M_or_MC &mapping, const DoFHandlerType &dof, const std::map*> &function_map, std::map &boundary_values, - const ComponentMask &component_mask, - const dealii::internal::int2type) + const ComponentMask &component_mask) { const unsigned int dim = DoFHandlerType::dimension; const unsigned int spacedim=DoFHandlerType::space_dimension; @@ -1810,214 +1719,267 @@ namespace VectorTools "for interior faces in your function map.")); const unsigned int n_components = DoFTools::n_components(dof); - const bool fe_is_system = (n_components != 1); - for (typename std::map*>::const_iterator i=function_map.begin(); i!=function_map.end(); ++i) Assert (n_components == i->second->n_components, ExcDimensionMismatch(n_components, i->second->n_components)); - // field to store the indices - std::vector face_dofs; - face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); - - std::vector > dof_locations; - dof_locations.reserve (DoFTools::max_dofs_per_face(dof)); - - // array to store the values of the boundary function at the boundary - // points. have two arrays for scalar and vector functions to use the - // more efficient one respectively - std::vector dof_values_scalar; - std::vector > dof_values_system; - 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 - dealii::hp::FECollection finite_elements (dof.get_fe()); - dealii::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)); - } + // interpolate boundary values in 1d. in higher dimensions, we + // use FEValues to figure out what to do on faces, but in 1d + // faces are points and it is far easier to simply work on + // individual vertices + if (dim == 1) + { + for (typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(); + cell != dof.end(); ++cell) + for (unsigned int direction=0; + direction::faces_per_cell; ++direction) + if (cell->at_boundary(direction) + && + (function_map.find(cell->face(direction)->boundary_id()) != function_map.end())) + { + const Function &boundary_function + = *function_map.find(cell->face(direction)->boundary_id())->second; + + // get the FE corresponding to this cell + const FiniteElement &fe = cell->get_fe(); + Assert (fe.n_components() == boundary_function.n_components, + ExcDimensionMismatch(fe.n_components(), + boundary_function.n_components)); + + Assert (component_mask.n_selected_components(fe.n_components()) > 0, + ComponentMask::ExcNoComponentSelected()); + + // now set the value of the vertex degree of + // freedom. setting also creates the entry in the + // map if it did not exist beforehand + // + // save some time by requesting values only once for + // each point, irrespective of the number of + // components of the function + Vector function_values (fe.n_components()); + if (fe.n_components() == 1) + function_values(0) + = boundary_function.value (cell->vertex(direction)); + else + boundary_function.vector_value (cell->vertex(direction), + function_values); + + for (unsigned int i=0; i + vertex_dof_index(direction,i, + cell->active_fe_index())] + = function_values(fe.face_system_to_component_index(i).first); + } } - // 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 - dealii::hp::MappingCollection mapping_collection (mapping); - dealii::hp::FEFaceValues x_fe_values (mapping_collection, finite_elements, q_collection, - update_quadrature_points); - - typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - if (!cell->is_artificial()) - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) + else // dim > 1 + { + const bool fe_is_system = (n_components != 1); + + // field to store the indices + std::vector face_dofs; + face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); + + std::vector > dof_locations; + dof_locations.reserve (DoFTools::max_dofs_per_face(dof)); + + // array to store the values of the boundary function at the boundary + // points. have two arrays for scalar and vector functions to use the + // more efficient one respectively + std::vector dof_values_scalar; + std::vector > dof_values_system; + 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 + dealii::hp::FECollection finite_elements (dof.get_fe()); + dealii::hp::QCollection q_collection; + for (unsigned int f=0; f &fe = cell->get_fe(); - - // we can presently deal only with primitive elements for - // boundary values. this does not preclude us using - // non-primitive elements in components that we aren't - // interested in, however. make sure that all shape functions - // that are non-zero for the components we are interested in, - // are in fact primitive - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + const FiniteElement &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 { - const ComponentMask &nonzero_component_array - = cell->get_fe().get_nonzero_components (i); - for (unsigned int c=0; cget_fe().is_primitive (i), - ExcMessage ("This function can only deal with requested boundary " - "values that correspond to primitive (scalar) base " - "elements")); + // 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)); } - - const typename DoFHandlerType::face_iterator face = cell->face(face_no); - const types::boundary_id boundary_component = face->boundary_id(); - - // see if this face is part of the boundaries for which we are - // supposed to do something, and also see if the finite element - // in use here has DoFs on the face at all - if ((function_map.find(boundary_component) != function_map.end()) - && - (cell->get_fe().dofs_per_face > 0)) + } + // 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 + dealii::hp::MappingCollection mapping_collection (mapping); + dealii::hp::FEFaceValues x_fe_values (mapping_collection, finite_elements, q_collection, + update_quadrature_points); + + typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (!cell->is_artificial()) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) { - // face is of the right component - x_fe_values.reinit(cell, face_no); - const dealii::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()); - const std::vector > &dof_locations - = fe_values.get_quadrature_points (); - - if (fe_is_system) + const FiniteElement &fe = cell->get_fe(); + + // we can presently deal only with primitive elements for + // boundary values. this does not preclude us using + // non-primitive elements in components that we aren't + // interested in, however. make sure that all shape functions + // that are non-zero for the components we are interested in, + // are in fact primitive + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) { - // resize array. avoid construction of a memory - // allocating temporary if possible - if (dof_values_system.size() < fe.dofs_per_face) - dof_values_system.resize (fe.dofs_per_face, - Vector(fe.n_components())); - else - dof_values_system.resize (fe.dofs_per_face); + const ComponentMask &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + for (unsigned int c=0; cget_fe().is_primitive (i), + ExcMessage ("This function can only deal with requested boundary " + "values that correspond to primitive (scalar) base " + "elements")); + } - function_map.find(boundary_component)->second - ->vector_value_list (dof_locations, dof_values_system); + const typename DoFHandlerType::face_iterator face = cell->face(face_no); + const types::boundary_id boundary_component = face->boundary_id(); - // enter those dofs into the list that match the - // component signature. avoid the usual complication - // that we can't just use *_system_to_component_index - // for non-primitive FEs - for (unsigned int i=0; iget_fe().dofs_per_face > 0)) + { + // face is of the right component + x_fe_values.reinit(cell, face_no); + const dealii::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()); + const std::vector > &dof_locations + = fe_values.get_quadrature_points (); + + if (fe_is_system) { - unsigned int component; - if (fe.is_primitive()) - component = fe.face_system_to_component_index(i).first; + // resize array. avoid construction of a memory + // allocating temporary if possible + if (dof_values_system.size() < fe.dofs_per_face) + dof_values_system.resize (fe.dofs_per_face, + Vector(fe.n_components())); else + dof_values_system.resize (fe.dofs_per_face); + + function_map.find(boundary_component)->second + ->vector_value_list (dof_locations, dof_values_system); + + // enter those dofs into the list that match the + // component signature. avoid the usual complication + // that we can't just use *_system_to_component_index + // for non-primitive FEs + for (unsigned int i=0; isecond - ->value_list (dof_locations, dof_values_scalar, 0); + else + // fe has only one component, so save some computations + { + // get only the one component that this function has + dof_values_scalar.resize (fe.dofs_per_face); + function_map.find(boundary_component)->second + ->value_list (dof_locations, dof_values_scalar, 0); - // enter into list + // enter into list - for (unsigned int i=0; i()); + component_mask_); } @@ -2066,8 +2027,7 @@ namespace VectorTools const ComponentMask &component_mask_) { do_interpolate_boundary_values (mapping, dof, function_map, boundary_values, - component_mask_, - dealii::internal::int2type()); + component_mask_); } -- 2.39.5