From: Martin Kronbichler Date: Thu, 16 Sep 2010 14:43:35 +0000 (+0000) Subject: Fix mixing interpolate_boundary_values with other types of constraints as e.g. hangin... X-Git-Tag: v8.0.0~5481 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=24217b0b1cbab235780f35747add6da16f9390d0;p=dealii.git Fix mixing interpolate_boundary_values with other types of constraints as e.g. hanging node constraints: We do not set boundary values to dofs that are already constrained. git-svn-id: https://svn.dealii.org/trunk@21992 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 2cec72f3be..ccfa873e2b 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -687,31 +687,43 @@ class VectorTools /** - * 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. - * + * Insert the (algebraic) constraints due + * to Dirichlet boundary conditions into + * a ConstraintMatrix @p + * constraints. This function identifies + * the degrees of freedom subject to + * Dirichlet boundary conditions, adds + * them to the list of constrained DoFs + * in @p constraints and sets the + * respective inhomogeneity to the value + * interpolated around the boundary. If + * this routine encounters a DoF that + * already is constrained (for instance + * by a hanging node constraint, see + * below, or any other type of + * constraint, e.g. from periodic + * boundary conditions), the old setting + * of the constraint (dofs the entry is + * constrained to, inhomogeneities) is + * kept and nothing happens. + * + * Please note that when combining + * adaptively refined meshes with hanging + * node constraints and inhomogeneous + * Dirichlet boundary conditions within + * one ConstraintMatrix object, the + * hanging node constraints should always + * be set first, and then Dirichlet + * boundary conditions should be + * interpolated. This makes sure that the + * discretization remains H1 + * conforming as is needed e.g. for the + * Laplace equation in 3D, as hanging + * nodes on the boundary should be still + * set to the weighted average of + * neighbors, and not the actual + * Dirichlet value. + * * * The parameter @p boundary_component * corresponds to the number @p * boundary_indicator of the face. 255 @@ -880,25 +892,42 @@ class VectorTools 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. + * Project a function to the boundary of + * the domain, using the given quadrature + * formula for the faces. This function + * identifies the degrees of freedom + * subject to Dirichlet boundary + * conditions, adds them to the list of + * constrained DoFs in @p constraints and + * sets the respective inhomogeneity to + * the value resulting from the + * projection operation. If this routine + * encounters a DoF that already is + * constrained (for instance by a hanging + * node constraint, see below, or any + * other type of constraint, e.g. from + * periodic boundary conditions), the old + * setting of the constraint (dofs the + * entry is constrained to, + * inhomogeneities) is kept and nothing + * happens. + * + * Please note that when combining + * adaptively refined meshes with hanging + * node constraints and inhomogeneous + * Dirichlet boundary conditions within + * one ConstraintMatrix object, the + * hanging node constraints should always + * be set first, and then Dirichlet + * boundary conditions should be + * interpolated. This makes sure that the + * discretization remains H1 + * conforming as is needed e.g. for the + * Laplace equation in 3D, as hanging + * nodes on the boundary should be still + * set to the weighted average of + * neighbors, and not the actual + * Dirichlet value. * * If @p component_mapping is empty, it * is assumed that the number of diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index f6d43aaf55..c975d44e56 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -1236,6 +1236,10 @@ VectorTools::interpolate_boundary_values (const Mapping &boundary_values, const std::vector &component_mask_) { + //ConstraintMatrix boundary_constraints(); + // interpolate_boundary_values (dof, boundary_component, boundary_function, + // boundary_constraints, component_mask); + const unsigned int dim=DH::dimension; const unsigned int spacedim=DH::space_dimension; @@ -1732,132 +1736,6 @@ 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, - outermost_cell->active_fe_index()); - 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 @@ -1869,317 +1747,20 @@ VectorTools::interpolate_boundary_values 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 boundary_values; + interpolate_boundary_values (mapping, dof, function_map, + boundary_values, component_mask_); + std::map::const_iterator boundary_value = + boundary_values.begin(); + for ( ; boundary_value !=boundary_values.end(); ++boundary_value) { - 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 + if (!constraints.is_constrained(boundary_value->first)) { - // 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)); - } + constraints.add_line (boundary_value->first); + constraints.set_inhomogeneity (boundary_value->first, + boundary_value->second); + } } - // 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 @@ -2437,26 +2017,6 @@ 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 @@ -2468,160 +2028,23 @@ VectorTools::project_boundary_values (const Mapping &mappin 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) + std::map boundary_values; + project_boundary_values (mapping, dof, boundary_functions, q, + boundary_values, component_mapping); + std::map::const_iterator boundary_value = + boundary_values.begin(); + for ( ; boundary_value !=boundary_values.end(); ++boundary_value) { -#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 + if (!constraints.is_constrained(boundary_value->first)) + { + constraints.add_line (boundary_value->first); + constraints.set_inhomogeneity (boundary_value->first, + boundary_value->second); + } } - 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 diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index fb7822eafb..f03a4e1be8 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -260,6 +260,21 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.

deal.II

    +
  1. +

    + Fixed: The methods VectorTools::interpolate_boundary_values and + VectorTools::project_boundary_values with ConstraintMatrix argument added + inhomogeneities also to DoFs that were already constrained when entering the + call without changing the list of DoFs that those were constrained to. This + lead to wrong results for inhomogeneous constraints in case e.g. + DoFTools::make_hanging_node_constraints was called before. Now, the correct + way for combining hanging node constraints and Dirichlet conditions in one + ConstraintMatrix is to first make the hanging node constraints, and then + interpolate/project the boundary values. +
    + (Luca Heltai, Martin Kronbichler 2010/09/16) +

    +

  2. Fixed: The method FEValuesViews::Vector::curl aborted the program in 2d under certain diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index d509a89f39..9aa1a95cb0 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -573,49 +573,52 @@ void StokesProblem::setup_dofs () // Now comes the implementation of // Dirichlet boundary conditions, which // should be evident after the discussion - // in the introduction. All that changed - // is that the function already appears - // in the setup functions, whereas we - // were used to see it in some assembly - // routine. Further down below where we - // set up the mesh, we will associate the - // top boundary where we impose Dirichlet - // boundary conditions with boundary - // indicator 1. We will have to pass - // this boundary indicator as second - // argument to the function below - // interpolating boundary values. There - // is one more thing, though. The - // function describing the Dirichlet - // conditions was defined for all - // components, both velocity and - // pressure. However, the Dirichlet - // conditions are to be set for the - // velocity only. To this end, we use a - // component_mask that + // in the introduction. All that changed is + // that the function already appears in the + // setup functions, whereas we were used to + // see it in some assembly routine. Further + // down below where we set up the mesh, we + // will associate the top boundary where we + // impose Dirichlet boundary conditions + // with boundary indicator 1. We will have + // to pass this boundary indicator as + // second argument to the function below + // interpolating boundary values. There is + // one more thing, though. The function + // describing the Dirichlet conditions was + // defined for all components, both + // velocity and pressure. However, the + // Dirichlet conditions are to be set for + // the velocity only. To this end, we use + // a component_mask that // filters out the pressure component, so // that the condensation is performed on - // velocity degrees of freedom - // only. Since we use adaptively refined - // grids the constraint matrix needs then - // also be filled with hanging node - // constraints needs to generated from - // the DoF handler. Note the order of the - // two functions — we first insert - // the boundary values into the - // constraint matrix, and then introduce - // the hanging node constraints. + // velocity degrees of freedom only. Since + // we use adaptively refined grids the + // constraint matrix needs to be first + // filled with hanging node constraints + // generated from the DoF handler. Note the + // order of the two functions — we + // first compute the hanging node + // constraints, and then insert the + // boundary values into the constraint + // matrix. This makes sure that we respect + // H1 conformity on boundaries + // with hanging nodes (in three space + // dimensions), where the hanging node + // needs to dominate the Dirichlet boundary + // values. { constraints.clear (); std::vector component_mask (dim+1, true); component_mask[dim] = false; + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); VectorTools::interpolate_boundary_values (dof_handler, 1, BoundaryValues(), constraints, component_mask); - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); } constraints.close (); diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h index 388c1290bc..6e454c7370 100644 --- a/deal.II/lac/include/lac/constraint_matrix.h +++ b/deal.II/lac/include/lac/constraint_matrix.h @@ -96,12 +96,20 @@ namespace internals * *

    Using the %ConstraintMatrix for Dirichlet boundary conditions

    * - * The ConstraintMatrix provides an alternative for implementinging - * Dirichlet boundary conditions (the standard way that is extensively - * discussed in the tutorial programs is to use - * MatrixTools::apply_boundary_values). The general principle of Dirichlet - * conditions are algebraic constraints of the form $x_{i} = b_i$, which - * fits into the form described above. + * The ConstraintMatrix provides an alternative for implementinging Dirichlet + * boundary conditions (the standard way that is extensively discussed in the + * tutorial programs is to use MatrixTools::apply_boundary_values). The + * general principle of Dirichlet conditions are algebraic constraints of the + * form $x_{i} = b_i$, which fits into the form described above. + * + * Please note that when combining adaptively refined meshes with hanging node + * constraints and inhomogeneous Dirichlet boundary conditions within one + * ConstraintMatrix object, the hanging node constraints should always be set + * first, and then Dirichlet boundary conditions should be interpolated. This + * makes sure that the discretization remains H1 conforming as is + * needed e.g. for the Laplace equation in 3D, as hanging nodes on the + * boundary should be still set to the weighted average of neighbors, and not + * the actual Dirichlet value. * * *

    Description of constraints

    @@ -779,6 +787,14 @@ class ConstraintMatrix : public Subscriptor const std::vector >* get_constraint_entries (unsigned int line) const; + /** + * Returns the value of the inhomogeneity + * stored in the constrained dof @p + * line. Unconstrained dofs also return a + * zero value. + */ + double get_inhomogeneity (unsigned int line) const; + /** * Print the constraint lines. Mainly * for debugging purposes. @@ -2100,6 +2116,18 @@ ConstraintMatrix::get_constraint_entries (unsigned int line) const +inline +double +ConstraintMatrix::get_inhomogeneity (unsigned int line) const +{ + if (is_constrained(line)) + return lines[lines_cache[calculate_line_index(line)]].inhomogeneity; + else + return 0; +} + + + inline unsigned int ConstraintMatrix::calculate_line_index (const unsigned int line) const {