From: kanschat Date: Wed, 31 Oct 2007 21:44:17 +0000 (+0000) Subject: make sure derivative degrees of freedom at the boundary do not cause problems X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9d0acbde636759a6dc2b236b364c9c8d351b7413;p=dealii-svn.git make sure derivative degrees of freedom at the boundary do not cause problems git-svn-id: https://svn.dealii.org/trunk@15423 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 ff2247e251..23960b172e 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -1689,7 +1690,11 @@ VectorTools::project_boundary_values (const Mapping &mapping, 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()); @@ -1715,7 +1720,7 @@ VectorTools::project_boundary_values (const Mapping &mapping, if (dim>=3) { #ifdef DEBUG -// Assert that there are no hanging nodes at the boundary +// 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) @@ -1745,11 +1750,33 @@ VectorTools::project_boundary_values (const Mapping &mapping, mass_matrix, boundary_functions, rhs, dof_to_boundary_mapping); - // same thing as above: if dim>=3 we need - // to consider constraints - Assert (dim<3, ExcNotImplemented()); - + // 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 ith 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 @@ -1764,12 +1791,15 @@ VectorTools::project_boundary_values (const Mapping &mapping, PreconditionSSOR<> prec; prec.initialize(mass_matrix, 1.2); + filtered_precondition.initialize(prec, true); // solve - cg.solve (mass_matrix, boundary_projection, rhs, prec); - + 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) + if (dof_to_boundary_mapping[i] != DoFHandler::invalid_dof_index + && ! excluded_dofs[dof_to_boundary_mapping[i]]) // this dof is on one of the // interesting boundary parts //