From: wolf Date: Thu, 6 Apr 2000 15:24:25 +0000 (+0000) Subject: Fix a bug in interpolate_boundary_values for 1d when discontinuous X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e7b2a72aee4c714728dddf939e9651a26aa0cf48;p=dealii-svn.git Fix a bug in interpolate_boundary_values for 1d when discontinuous elements are involved. git-svn-id: https://svn.dealii.org/trunk@2686 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 4771820ac0..fbc735a53d 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -242,7 +242,7 @@ void VectorTools::interpolate (const DoFHandler &dof, cell->get_dof_indices (dofs_on_cell); -if (fe_is_system) + if (fe_is_system) { // get function values at // these points. Here: get @@ -408,7 +408,7 @@ void VectorTools::project (const DoFHandler &dof, }; -// set up mass matrix and right hand side + // set up mass matrix and right hand side vec.reinit (dof.n_dofs()); SparsityPattern sparsity(dof.n_dofs(), dof.n_dofs(), @@ -501,8 +501,6 @@ VectorTools::interpolate_boundary_values (const DoFHandler<1> &dof, const FiniteElement<1> &fe = dof.get_fe(); Assert (fe.n_components() == boundary_function.n_components, - ExcComponentMismatch()); - Assert (fe.dofs_per_vertex == fe.n_components(), ExcComponentMismatch()); // set the component mask to either @@ -514,56 +512,51 @@ VectorTools::interpolate_boundary_values (const DoFHandler<1> &dof, Assert (count(component_mask.begin(), component_mask.end(), true) > 0, ExcComponentMismatch()); - // check whether boundary values at the - // left boundary of the line are requested - if (boundary_component == 0) - { - // first find the leftmost active - // cell by first traversing the coarse - // grid to its left end and then going - // to the children - DoFHandler<1>::cell_iterator leftmost_cell = dof.begin(0); - while (leftmost_cell->neighbor(0).state() == valid) - leftmost_cell = leftmost_cell->neighbor(0); - - while (leftmost_cell->has_children()) - leftmost_cell = leftmost_cell->child(0); - - // now set the value of the - // leftmost degree of - // freedom. setting also - // created the entry in the map - // if it did not exist - // beforehand - for (unsigned int i=0; ivertex_dof_index(0,i)] - = boundary_function.value(leftmost_cell->vertex(0), i); - }; - - // same for the right boundary of - // the line are requested - if (boundary_component == 1) - { - // first find the leftmost active - // cell by first traversing the coarse - // grid to its left end and then going - // to the children - DoFHandler<1>::cell_iterator rightmost_cell = dof.last(0); - while (rightmost_cell->neighbor(1).state() == valid) - rightmost_cell = rightmost_cell->neighbor(1); - - while (rightmost_cell->has_children()) - rightmost_cell = rightmost_cell->child(1); - - // now set the value of the rightmost - // degree of freedom - for (unsigned int i=0; ivertex_dof_index(1,i)] - = boundary_function.value(rightmost_cell->vertex(1), i); - }; + // 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 + DoFHandler<1>::cell_iterator outermost_cell = dof.begin(0); + while (outermost_cell->neighbor(direction).state() == valid) + outermost_cell = outermost_cell->neighbor(direction); + + while (outermost_cell->has_children()) + outermost_cell = outermost_cell->child(direction); + + // 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)] + = function_values(fe.face_system_to_component_index(i).first); }; @@ -640,8 +633,8 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // get only the one component that // this function has boundary_function.value_list (dof_locations, - dof_values_scalar, - 0); + dof_values_scalar, + 0); // enter into list @@ -693,12 +686,12 @@ VectorTools::project_boundary_values (const DoFHandler &dof, Assert (false, ExcNotImplemented()); -// make mass matrix and right hand side + // make mass matrix and right hand side SparseMatrix mass_matrix(sparsity); Vector rhs(sparsity.n_rows()); -MatrixTools::create_boundary_mass_matrix (dof, q, + MatrixTools::create_boundary_mass_matrix (dof, q, mass_matrix, boundary_functions, rhs, dof_to_boundary_mapping); @@ -707,7 +700,7 @@ MatrixTools::create_boundary_mass_matrix (dof, q, Assert (dim<3, ExcNotImplemented()); -Vector boundary_projection (rhs.size()); + Vector boundary_projection (rhs.size()); SolverControl control(1000, 1e-16); PrimitiveVectorMemory<> memory;