From e2edf401c7abddaf50292bea833b782d0fd123df Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 17 Feb 2003 15:59:20 +0000 Subject: [PATCH] Fix another bug with non-primitive elements. git-svn-id: https://svn.dealii.org/trunk@7125 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 18 ++ deal.II/deal.II/source/numerics/vectors.cc | 118 +++++++- tests/bits/Makefile | 3 +- tests/bits/anna_6.cc | 326 +++++++++++++++++++++ tests/bits/anna_6.checked | 11 + 5 files changed, 469 insertions(+), 7 deletions(-) create mode 100644 tests/bits/anna_6.cc create mode 100644 tests/bits/anna_6.checked diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 19d93ac6fc..cdba28b0af 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -620,6 +620,20 @@ class VectorTools * 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. */ @@ -922,6 +936,10 @@ class VectorTools /** * Exception */ + DeclException0 (ExcFENotPrimitive); + /** + * Exception + */ DeclException0 (ExcInvalidBoundaryIndicator); /** * Exception diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 116ba292f2..1bf6f6a19d 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -865,14 +865,120 @@ interpolate_boundary_values (const Mapping &mapping, if (fe_is_system) { - function_map.find(boundary_component)->second->vector_value_list (dof_locations, - dof_values_system); + function_map.find(boundary_component)->second + ->vector_value_list (dof_locations, dof_values_system); - // enter into list + // 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; i(-1)))); + Assert (cell_i < fe.dofs_per_cell, ExcInternalError()); + + // make sure + // that if + // this is + // not a + // primitive + // shape function, + // then all + // the + // corresponding + // components + // in the + // mask are + // not set + if (!fe.is_primitive(cell_i)) + for (unsigned int c=0; c +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + + +#include +#include +#include + +#include +#include + + +template +class ImposeBC +{ + public: + ImposeBC(); + ~ImposeBC(); + void run (); + + private: + + void get_ready (); + void test_extract_boundary_DoFs (); + void test_interpolate_BC (); + + Triangulation triangulation; + + // We use a FE-System with 2 components: + // a vector-valued one for the + // field-variable u and a scalar one for + // the pressure-variable p + FESystem fe; + DoFHandler dof_handler; + + unsigned int n_u_dofs; + unsigned int n_p_dofs; + +}; + + +// some boundary function for the scalar component +template +class BoundaryFunction : public Function +{ + public: + BoundaryFunction (); + + virtual void vector_value (const Point &p, + Vector &values) const; + +}; + + +template +BoundaryFunction::BoundaryFunction () : + Function (dim+1) {}; + + + +template +inline +void BoundaryFunction::vector_value (const Point &p, + Vector &values) const +{ + + Assert (values.size() == dim+1, + ExcDimensionMismatch (values.size(), dim+1)); + + values.clear (); + values(dim) = 1.; +}; + + + + + + // Construct FE with first component: Nedelec-Element, + // second component: Q1_Element +template +ImposeBC::ImposeBC() : + fe (FE_Nedelec(1), 1, FE_Q(1), 1), + dof_handler (triangulation) + + +{}; + + + +template +ImposeBC::~ImposeBC() +{ + dof_handler.clear (); +}; + + +template +void ImposeBC::get_ready () +{ + dof_handler.distribute_dofs (fe); + std::vector dofs_per_comp; + DoFTools::count_dofs_per_component(dof_handler, dofs_per_comp); + + // For an FESystem with Nedelec-elements as + // first component and bilinear elements as + // component we have: + // dofs_per_comp[0] = dofs_per_comp[1] = # Ned-DoFs + // dofs_per_comp[2] = # Q1-DoFs + n_u_dofs = dofs_per_comp[0]; + n_p_dofs = dofs_per_comp[2]; + +}; + + +template +void ImposeBC::test_extract_boundary_DoFs () +{ + + std::map boundary_values; + std::vector bc_component_select(dim + 1); + + // extract boundary DoFs for the Nedelec-component + // and impose zero boundary condition + bc_component_select[0] = true; + bc_component_select[1] = true; + bc_component_select[2] = false; + + std::vector ned_boundary_dofs (dof_handler.n_dofs()); + std::set boundary_indicators; + boundary_indicators.insert (0); + DoFTools::extract_boundary_dofs (dof_handler, + bc_component_select, + ned_boundary_dofs, + boundary_indicators); + + + for (unsigned int i=0; i +void ImposeBC::test_interpolate_BC () +{ + + std::map boundary_values; + std::vector bc_component_select(dim + 1, false); + + + // impose inhomogeneous boundary condition + // on the scalar variable + bc_component_select.back() = true; + + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryFunction(), + boundary_values, + bc_component_select); + + + + // check + // (the pressure is assumed to be set to 1 + // on the boundary) + std::vector p_boundary_dofs (dof_handler.n_dofs()); + std::set boundary_indicators; + boundary_indicators.insert (0); + DoFTools::extract_boundary_dofs (dof_handler, + bc_component_select, + p_boundary_dofs, + boundary_indicators); + for (unsigned int i=0; i +void ImposeBC::run () +{ + GridGenerator::hyper_cube(triangulation, -1,1); + triangulation.refine_global (1); + triangulation.begin_active()->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + get_ready (); + + deallog << " Total number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl + << " Number of degrees of freedom for the field variable U: " + << n_u_dofs + << std::endl + << " Number of degrees of freedom for the pressure variable p: " + << n_p_dofs + << std::endl; + + test_extract_boundary_DoFs (); + test_interpolate_BC (); + +}; + + +int main () +{ + try + { + std::ofstream logfile("anna_6.output"); + logfile.precision (2); + deallog.attach(logfile); + deallog.depth_console(0); + + ImposeBC<2>().run (); + ImposeBC<3>().run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + return 0; +}; diff --git a/tests/bits/anna_6.checked b/tests/bits/anna_6.checked new file mode 100644 index 0000000000..e91fa2436c --- /dev/null +++ b/tests/bits/anna_6.checked @@ -0,0 +1,11 @@ + +DEAL:: Number of active cells: 7 +DEAL:: Total number of degrees of freedom: 36 +DEAL:: Number of degrees of freedom for the field variable U: 22 +DEAL:: Number of degrees of freedom for the pressure variable p: 14 +DEAL::1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +DEAL:: Number of active cells: 15 +DEAL:: Total number of degrees of freedom: 151 +DEAL:: Number of degrees of freedom for the field variable U: 105 +DEAL:: Number of degrees of freedom for the pressure variable p: 105 +DEAL::1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -- 2.39.5