From 3bc00a998be7fcab6fef139ba7d05c2b81c17c56 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 10 Mar 2009 16:55:43 +0000 Subject: [PATCH] Add functions interpolate_boundary_values and project_boundary_values for a ConstraintMatrix argument by using inhomogeneous constraints instead of std::map. There is a lot of duplicated code in vectors.templates.h right now, need to fix that later. git-svn-id: https://svn.dealii.org/trunk@18474 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 264 ++++++- .../include/numerics/vectors.templates.h | 721 +++++++++++++++++- deal.II/deal.II/source/numerics/vectors.cc | 40 + 3 files changed, 986 insertions(+), 39 deletions(-) diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 92c14cad83..e63d0b2e5f 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -682,53 +682,179 @@ class VectorTools const typename FunctionMap::type &function_map, std::map &boundary_values, const std::vector &component_mask = std::vector()); + + + /** + * Insert the (algebraic) constraints + * due to Dirichlet boundary conditions + * to the ConstraintMatrix. This + * function makes up the list of + * degrees of freedom subject to + * Dirichlet boundary conditions and + * the values to be assigned to them, + * by interpolation around the + * boundary. If the ConstraintMatrix @p + * constraints contained values or + * other constraints before, the new + * ones are added, or the old ones + * overwritten if a node of the + * boundary part to be used was already + * in the list of constraints. This is + * handled by using inhomogeneous + * constraints. Please note that when + * combining adaptive meshes and this + * kind of constraints, the Dirichlet + * conditions should be set first, and + * then completed by hanging node + * constraints, in order to make sure + * that the discretization remains + * consistent. + * + * The parameter @p boundary_component + * corresponds to the number @p + * boundary_indicator of the face. 255 + * is an illegal value, since it is + * reserved for interior faces. + * + * The flags in the last parameter, @p + * component_mask denote which + * components of the finite element + * space shall be interpolated. If it + * is left as specified by the default + * value (i.e. an empty array), all + * components are interpolated. If it + * is different from the default value, + * it is assumed that the number of + * entries equals the number of + * components in the boundary functions + * and the finite element, and those + * components in the given boundary + * function will be used for which the + * respective flag was set in the + * component mask. + * + * It is assumed that the number of + * components of the function in @p + * boundary_function matches that of + * the finite element used by @p dof. + * + * If the finite element used has shape + * functions that are non-zero in more + * than one component (in deal.II + * speak: they are non-primitive), then + * these components can presently not + * be used for interpolating boundary + * values. Thus, the elements in the + * component mask corresponding to the + * components of these non-primitive + * shape functions must be @p false. + * + * See the general doc for more + * information. + */ + template + static + void + interpolate_boundary_values (const Mapping &mapping, + const DH &dof, + const typename FunctionMap::type &function_map, + ConstraintMatrix &constraints, + const std::vector &component_mask = std::vector()); + /** + * @deprecated This function is there + * mainly for backward compatibility. + * + * Same function as above, but taking + * only one pair of boundary indicator + * and corresponding boundary + * function. Calls the other function + * with remapped arguments. + * + */ + template + static + void + interpolate_boundary_values (const Mapping &mapping, + const DH &dof, + const unsigned char boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const std::vector &component_mask = std::vector()); + + /** + * Calls the other + * interpolate_boundary_values() + * function, see above, with + * mapping=MappingQ1@(). + */ + template + static + void + interpolate_boundary_values (const DH &dof, + const unsigned char boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const std::vector &component_mask = std::vector()); + + /** + * Calls the other + * interpolate_boundary_values() + * function, see above, with + * mapping=MappingQ1@(). + */ + template + static + void + interpolate_boundary_values (const DH &dof, + const typename FunctionMap::type &function_map, + ConstraintMatrix &constraints, + const std::vector &component_mask = std::vector()); + + /** * Project a function to the boundary - * of the domain, using the given quadrature - * formula for the faces. If the - * @p boundary_values contained values - * before, the new ones are added, or - * the old one overwritten if a node - * of the boundary part to be projected - * on already was in the variable. - * - * If @p component_mapping is - * empty, it is assumed that the - * number of components of @p - * boundary_function matches that - * of the finite element used by - * @p dof. + * of the domain, using the given + * quadrature formula for the faces. If + * the @p boundary_values contained + * values before, the new ones are + * added, or the old one overwritten if + * a node of the boundary part to be + * projected on already was in the + * variable. + * + * If @p component_mapping is empty, it + * is assumed that the number of + * components of @p boundary_function + * matches that of the finite element + * used by @p dof. * * In 1d, projection equals * interpolation. Therefore, * interpolate_boundary_values is * called. * - * @arg @p boundary_values: the - * result of this function, a map - * containing all indices of - * degrees of freedom at the - * boundary (as covered by the + * @arg @p boundary_values: the result + * of this function, a map containing + * all indices of degrees of freedom at + * the boundary (as covered by the * boundary parts in @p - * boundary_functions) and the - * computed dof value for this - * degree of freedom. - * - * @arg @p component_mapping: if - * the components in @p - * boundary_functions and @p dof - * do not coincide, this vector - * allows them to be remapped. If - * the vector is not empty, it - * has to have one entry for each - * component in @p dof. This - * entry is the component number - * in @p boundary_functions that - * should be used for this - * component in @p dof. By - * default, no remapping is + * boundary_functions) and the computed + * dof value for this degree of + * freedom. + * + * @arg @p component_mapping: if the + * components in @p boundary_functions + * and @p dof do not coincide, this + * vector allows them to be + * remapped. If the vector is not + * empty, it has to have one entry for + * each component in @p dof. This entry + * is the component number in @p + * boundary_functions that should be + * used for this component in @p + * dof. By default, no remapping is * applied. */ template @@ -751,6 +877,72 @@ class VectorTools std::map &boundary_values, std::vector component_mapping = std::vector()); + /** + * Project a function to the boundary + * of the domain, using the given + * quadrature formula for the faces. If + * the ConstraintMatrix @p constraints + * contained values or other + * constraints before, the new ones are + * added, or the old ones overwritten + * if a node of the boundary part to be + * used was already in the list of + * constraints. This is handled by + * using inhomogeneous + * constraints. Please note that when + * combining adaptive meshes and this + * kind of constraints, the Dirichlet + * conditions should be set first, and + * then completed by hanging node + * constraints, in order to make sure + * that the discretization remains + * consistent. + * + * If @p component_mapping is empty, it + * is assumed that the number of + * components of @p boundary_function + * matches that of the finite element + * used by @p dof. + * + * In 1d, projection equals + * interpolation. Therefore, + * interpolate_boundary_values is + * called. + * + * @arg @p component_mapping: if the + * components in @p boundary_functions + * and @p dof do not coincide, this + * vector allows them to be + * remapped. If the vector is not + * empty, it has to have one entry for + * each component in @p dof. This entry + * is the component number in @p + * boundary_functions that should be + * used for this component in @p + * dof. By default, no remapping is + * applied. + */ + template + static void project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping = std::vector()); + + /** + * Calls the project_boundary_values() + * function, see above, with + * mapping=MappingQ1@(). + */ + template + static void project_boundary_values (const DoFHandler &dof, + const typename FunctionMap::type &boundary_function, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping = std::vector()); + + /** * Compute the constraints that * correspond to boundary conditions of diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 1677cbd99f..350ad1504f 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -1213,6 +1213,7 @@ VectorTools::create_boundary_right_hand_side (const hp::DoFHandler +// ----------- interpolate_boundary_values for std::map -------------------- #if deal_II_dimension == 1 @@ -1223,10 +1224,10 @@ VectorTools::create_boundary_right_hand_side (const hp::DoFHandler template void -VectorTools::interpolate_boundary_values (const Mapping &, +VectorTools::interpolate_boundary_values (const Mapping &, const DH &dof, - const unsigned char boundary_component, - const Function &boundary_function, + const unsigned char boundary_component, + const Function &boundary_function, std::map &boundary_values, const std::vector &component_mask_) { @@ -1720,6 +1721,516 @@ VectorTools::interpolate_boundary_values (const DH &dof, } + + +// ----------- interpolate_boundary_values for ConstraintMatrix -------------- + +// TODO (M.K.): There is a lot of duplicated code with the above +// function. We should unify all these interpolate_boundary_values functions +// in one way or the other. + +#if deal_II_dimension == 1 + +//TODO[?] Actually the Mapping object should be a MappingCollection object for the +// hp::DoFHandler. + +//template class DH, int spacedim> + +template +void +VectorTools::interpolate_boundary_values (const Mapping &, + const DH &dof, + const unsigned char boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const std::vector &component_mask_) +{ + const unsigned int dim=DH::dimension; + const unsigned int spacedim=DH::space_dimension; + + Assert (boundary_component != 255, + ExcInvalidBoundaryIndicator()); + Assert ((component_mask_.size() == 0) || + (component_mask_.size() == 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.")); + + // check whether boundary values at + // the left or right boundary of + // the line are + // requested. direction denotes + // the neighboring direction in + // which we seek the boundary, + // i.e. 0 is left boundary and 1 is + // right. + const unsigned int direction = boundary_component; + Assert (direction < 2, ExcInvalidBoundaryIndicator()); + + // first find the outermost active + // cell by first traversing the coarse + // grid to its end and then going + // to the children + typename DH::cell_iterator outermost_cell = dof.begin(0); + while (outermost_cell->neighbor(direction).state() == IteratorState::valid) + outermost_cell = outermost_cell->neighbor(direction); + + while (outermost_cell->has_children()) + outermost_cell = outermost_cell->child(direction); + + // get the FE corresponding to this + // cell + const FiniteElement &fe = outermost_cell->get_fe(); + Assert (fe.n_components() == boundary_function.n_components, + ExcDimensionMismatch(fe.n_components(), boundary_function.n_components)); + + // set the component mask to either + // the original value or a vector + // of trues + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (fe.n_components(), true) : + component_mask_); + Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, + ExcNoComponentSelected()); + + // now set the value of the + // outermost 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 (outermost_cell->vertex(direction)); + else + boundary_function.vector_value (outermost_cell->vertex(direction), + function_values); + + for (unsigned int i=0; ivertex_dof_index(direction,i); + constraints.add_line (row); + constraints.set_inhomogeneity (row, + function_values(fe.face_system_to_component_index(i).first)); + } +} + + +//TODO[?] Actually the Mapping object should be a MappingCollection object for the +// hp::DoFHandler. + +// Implementation for 1D +template +void +VectorTools::interpolate_boundary_values + (const Mapping &mapping, + const DH &dof, + const typename FunctionMap::type &function_map, + ConstraintMatrix &constraints, + const std::vector &component_mask) +{ + for (typename FunctionMap::type::const_iterator i=function_map.begin(); + i!=function_map.end(); ++i) + interpolate_boundary_values (mapping, dof, i->first, *i->second, + constraints, component_mask); +} + + +//TODO[?] Actually the Mapping object should be a MappingCollection object for the +// hp::DoFHandler. + + +#else + + +template +void +VectorTools::interpolate_boundary_values + (const Mapping &mapping, + const DH &dof, + const typename FunctionMap::type &function_map, + ConstraintMatrix &constraints, + const std::vector &component_mask_) +{ + const unsigned int dim=DH::dimension; + + Assert ((component_mask_.size() == 0) || + (component_mask_.size() == 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; + + Assert (function_map.find(255) == function_map.end(), + ExcInvalidBoundaryIndicator()); + + const unsigned int n_components = DoFTools::n_components(dof); + const bool fe_is_system = (n_components != 1); + + for (typename FunctionMap::type::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)); + + // set the component mask to either + // the original value or a vector + // of trues + const std::vector component_mask ((component_mask_.size() == 0) ? + std::vector (n_components, true) : + component_mask_); + Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, + ExcNoComponentSelected()); + + // 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 to 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 + 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) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + 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) + { + const std::vector &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")); + } + + 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 + 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()); + const std::vector > &dof_locations + = fe_values.get_quadrature_points (); + + if (fe_is_system) + { + // 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); + + // enter into list + + for (unsigned int i=0; i +void +VectorTools::interpolate_boundary_values + (const Mapping &mapping, + const DH &dof, + const unsigned char boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const std::vector &component_mask) +{ + typename FunctionMap::type function_map; + function_map[boundary_component] = &boundary_function; + interpolate_boundary_values (mapping, dof, function_map, constraints, + component_mask); +} + +#endif + + +template +void +VectorTools::interpolate_boundary_values + (const DH &dof, + const unsigned char boundary_component, + const Function &boundary_function, + ConstraintMatrix &constraints, + const std::vector &component_mask) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + interpolate_boundary_values(StaticMappingQ1::mapping, + dof, boundary_component, + boundary_function, constraints, component_mask); +} + + + +template +void +VectorTools::interpolate_boundary_values + (const DH &dof, + const typename FunctionMap::type &function_map, + ConstraintMatrix &constraints, + const std::vector &component_mask) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + interpolate_boundary_values(StaticMappingQ1::mapping, + dof, function_map, + constraints, component_mask); +} + + + + +// -------- implementation for project_boundary_values with std::map -------- + #if deal_II_dimension == 1 // Implementation for 1D @@ -1916,6 +2427,210 @@ VectorTools::project_boundary_values (const DoFHandler &dof, +// ----- implementation for project_boundary_values with ConstraintMatrix ----- + +#if deal_II_dimension == 1 + +// Implementation for 1D +template +void +VectorTools::project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &, + ConstraintMatrix &constraints, + std::vector component_mapping) +{ + Assert (component_mapping.size() == 0, ExcNotImplemented()); + // projection in 1d is equivalent + // to interpolation + interpolate_boundary_values (mapping, dof, boundary_functions, + constraints, std::vector()); +} + +#else + + +template +void +VectorTools::project_boundary_values (const Mapping &mapping, + const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping) +{ +//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity +// structures, matrices and right hand sides or distribution of +// solution vectors is performed. This is ok for dim<3 because then +// there are no constrained nodes on the boundary, but is not +// acceptable for higher dimensions. Fix this. + + if (component_mapping.size() == 0) + { + AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components); + // I still do not see why i + // should create another copy + // here + component_mapping.resize(dof.get_fe().n_components()); + for (unsigned int i= 0 ;i < component_mapping.size() ; ++i) + component_mapping[i] = i; + } + else + AssertDimension (dof.get_fe().n_components(), component_mapping.size()); + + std::vector dof_to_boundary_mapping; + std::set selected_boundary_components; + for (typename FunctionMap::type::const_iterator i=boundary_functions.begin(); + i!=boundary_functions.end(); ++i) + selected_boundary_components.insert (i->first); + + DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components, + dof_to_boundary_mapping); + + // Done if no degrees of freedom on + // the boundary + if (dof.n_boundary_dofs (boundary_functions) == 0) + return; + // set up sparsity structure + SparsityPattern sparsity(dof.n_boundary_dofs (boundary_functions), + dof.max_couplings_between_boundary_dofs()); + DoFTools::make_boundary_sparsity_pattern (dof, + boundary_functions, + dof_to_boundary_mapping, + sparsity); + + // note: for three or more dimensions, there + // may be constrained nodes on the boundary + // in this case the boundary mass matrix has + // to be condensed and the solution is to + // be distributed afterwards, which is not + // yet implemented. The reason for this is + // that we cannot simply use the condense + // family of functions, since the matrices + // and vectors do not use the global + // numbering but rather the boundary + // numbering, i.e. the condense function + // needs to use another indirection. There + // should be not many technical problems, + // but it needs to be implemented + if (dim>=3) + { +#ifdef DEBUG +// Assert that there are no hanging nodes at the boundary + int level = -1; + for (typename DoFHandler::active_cell_iterator cell = dof.begin_active(); + cell != dof.end(); ++cell) + for (unsigned int f=0;f::faces_per_cell;++f) + { + if (cell->at_boundary(f)) + { + if (level == -1) + level = cell->level(); + else + { + Assert (level == cell->level(), ExcNotImplemented()); + } + } + } +#endif + } + sparsity.compress(); + + + // make mass matrix and right hand side + SparseMatrix mass_matrix(sparsity); + Vector rhs(sparsity.n_rows()); + + + MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, + mass_matrix, boundary_functions, + rhs, dof_to_boundary_mapping, (const Function*) 0, + component_mapping); + + // For certain weird elements, + // there might be degrees of + // freedom on the boundary, but + // their shape functions do not + // have support there. Let's + // eliminate them here. + +//TODO: Maybe we should figure out if the element really needs this + + FilteredMatrix > filtered_mass_matrix(mass_matrix, true); + FilteredMatrix > filtered_precondition; + std::vector excluded_dofs(mass_matrix.m(), false); + + double max_element = 0.; + for (unsigned int i=0;i max_element) + max_element = mass_matrix.diag_element(i); + + for (unsigned int i=0;i boundary_projection (rhs.size()); + + // Allow for a maximum of 5*n + // steps to reduce the residual by + // 10^-12. n steps may not be + // sufficient, since roundoff + // errors may accumulate for badly + // conditioned matrices + ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false); + GrowingVectorMemory<> memory; + SolverCG<> cg(control,memory); + + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); + filtered_precondition.initialize(prec, true); + // solve + cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition); + filtered_precondition.apply_constraints(boundary_projection, true); + filtered_precondition.clear(); + // fill in boundary values + for (unsigned int i=0; i::invalid_dof_index + && ! excluded_dofs[dof_to_boundary_mapping[i]]) + // this dof is on one of the + // interesting boundary parts + // + // remember: i is the global dof + // number, dof_to_boundary_mapping[i] + // is the number on the boundary and + // thus in the solution vector + { + // TODO: check whether we should clear + // the entries in this line first. + constraints.add_line (i); + constraints.set_inhomogeneity (i, boundary_projection(dof_to_boundary_mapping[i])); + } +} + +#endif + +template +void +VectorTools::project_boundary_values (const DoFHandler &dof, + const typename FunctionMap::type &boundary_functions, + const Quadrature &q, + ConstraintMatrix &constraints, + std::vector component_mapping) +{ + Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); + project_boundary_values(StaticMappingQ1::mapping, dof, boundary_functions, q, + constraints, component_mapping); +} + + + + namespace internal { namespace VectorTools diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 542ffbc45d..0f1ffa8e85 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -137,6 +137,30 @@ void VectorTools::interpolate_boundary_values ( std::map &, const std::vector &); +template +void VectorTools::interpolate_boundary_values ( + const DoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const hp::DoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const MGDoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + template void VectorTools::project_boundary_values (const Mapping &, @@ -153,6 +177,22 @@ void VectorTools::project_boundary_values std::map&, std::vector); +template +void VectorTools::project_boundary_values +(const Mapping &, + const DoFHandler &, + const FunctionMap::type &, + const Quadrature&, + ConstraintMatrix&, std::vector); + +template +void VectorTools::project_boundary_values +(const DoFHandler &, + const FunctionMap::type &, + const Quadrature&, + ConstraintMatrix&, std::vector); + + #if deal_II_dimension != 1 template void -- 2.39.5