From: Martin Kronbichler Date: Thu, 5 Nov 2009 20:57:09 +0000 (+0000) Subject: A boundary mass matrix does not make much sense in 1D. X-Git-Tag: v8.0.0~6850 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=37ad314b73b34d0de42cfd1d0d1921933c3d84d6;p=dealii.git A boundary mass matrix does not make much sense in 1D. git-svn-id: https://svn.dealii.org/trunk@20042 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 266c76ee1f..3674578717 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -48,21 +48,21 @@ DEAL_II_NAMESPACE_OPEN template -void VectorTools::interpolate (const Mapping &mapping, - const DH &dof, - const Function &function, - VECTOR &vec) +void VectorTools::interpolate (const Mapping &mapping, + const DH &dof, + const Function &function, + VECTOR &vec) { const unsigned int dim=DH::dimension; Assert (dof.get_fe().n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe().n_components(), function.n_components)); - + const hp::FECollection fe (dof.get_fe()); const unsigned int n_components = fe.n_components(); const bool fe_is_system = (n_components != 1); - + typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); @@ -88,11 +88,11 @@ void VectorTools::interpolate (const Mapping Assert (unit_support_points[fe_index].size() != 0, ExcNonInterpolatingFE()); } - - // Find the support points + + // Find the support points // on a cell that - // are multiply mentioned in + // are multiply mentioned in // unit_support_points. // Mark the first representative // of each multiply mentioned @@ -104,7 +104,7 @@ void VectorTools::interpolate (const Mapping // the following vector collects all dofs i, // 0<=i std::vector > dof_to_rep_index_table(fe.size()); std::vector n_rep_points (fe.size(), 0); - + for (unsigned int fe_index=0; fe_index // support points that are placed // one after the other. for (unsigned int j=dofs_of_rep_points[fe_index].size(); j>0; --j) - if (unit_support_points[fe_index][i] + if (unit_support_points[fe_index][i] == unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]]) { dof_to_rep_index_table[fe_index].push_back(j-1); representative=false; break; } - + if (representative) { // rep_index=dofs_of_rep_points.size() @@ -151,12 +151,12 @@ void VectorTools::interpolate (const Mapping Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell, ExcInternalError()); } - + const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(), n_rep_points.end()); std::vector dofs_on_cell (fe.max_dofs_per_cell()); std::vector > rep_points (max_rep_points); - + // get space for the values of the // function at the rep support points. // @@ -175,21 +175,21 @@ void VectorTools::interpolate (const Mapping // Transformed support points are computed by // FEValues hp::MappingCollection mapping_collection (mapping); - + hp::FEValues fe_values (mapping_collection, fe, support_quadrature, update_quadrature_points); - + for (; cell!=endc; ++cell) { const unsigned int fe_index = cell->active_fe_index(); - + // for each cell: // get location of finite element // support_points fe_values.reinit(cell); const std::vector >& support_points = fe_values.get_present_fe_values().get_quadrature_points(); - + // pick out the representative // support points rep_points.resize (dofs_of_rep_points[fe_index].size()); @@ -197,7 +197,7 @@ void VectorTools::interpolate (const Mapping rep_points[j] = support_points[dofs_of_rep_points[fe_index][j]]; // get indices of the dofs on this cell - dofs_on_cell.resize (fe[fe_index].dofs_per_cell); + dofs_on_cell.resize (fe[fe_index].dofs_per_cell); cell->get_dof_indices (dofs_on_cell); @@ -236,7 +236,7 @@ void VectorTools::interpolate (const Mapping // values to the global // vector for (unsigned int i=0; i::mapping, + interpolate(StaticMappingQ1::mapping, dof, function, vec); } @@ -269,27 +269,27 @@ VectorTools::interpolate (const DoFHandler &dof_1, std::vector touch_count (dof_2.n_dofs(), 0); std::vector local_dof_indices (dof_2.get_fe().dofs_per_cell); - + typename DoFHandler::active_cell_iterator h = dof_1.begin_active(); typename DoFHandler::active_cell_iterator l = dof_2.begin_active(); const typename DoFHandler::cell_iterator endh = dof_1.end(); - + for(; h != endh; ++h, ++l) { h->get_dof_values(data_1, cell_data_1); transfer.vmult(cell_data_2, cell_data_1); l->get_dof_indices (local_dof_indices); - + // distribute cell vector - for (unsigned int j=0; j &dof_1, { Assert (touch_count[i] != 0, ExcInternalError()); - + data_2(i) /= touch_count[i]; }; } @@ -315,7 +315,7 @@ namespace internal template void - interpolate_zero_boundary_values (const dealii::DoFHandler &dof_handler, + interpolate_zero_boundary_values (const dealii::DoFHandler &dof_handler, std::map &boundary_values) { // we only need to find the @@ -337,8 +337,8 @@ namespace internal //codimension 1 void - interpolate_zero_boundary_values (const dealii::DoFHandler<1,2> &dof_handler, - std::map &boundary_values) + interpolate_zero_boundary_values (const dealii::DoFHandler<1,2> &dof_handler, + std::map &boundary_values) { // we only need to find the // left-most and right-most @@ -358,11 +358,11 @@ namespace internal #else - + template void - interpolate_zero_boundary_values (const dealii::DoFHandler &dof_handler, - std::map &boundary_values) + interpolate_zero_boundary_values (const dealii::DoFHandler &dof_handler, + std::map &boundary_values) { const FiniteElement &fe = dof_handler.get_fe(); @@ -405,7 +405,7 @@ namespace internal boundary_values[face_dof_indices[i]] = 0.; } } - + #endif } } @@ -413,15 +413,15 @@ namespace internal template -void VectorTools::project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VECTOR &vec_result, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) +void VectorTools::project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VECTOR &vec_result, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) { Assert (dof.get_fe().n_components() == function.n_components, ExcDimensionMismatch(dof.get_fe().n_components(), @@ -429,18 +429,18 @@ void VectorTools::project (const Mapping &mapping, Assert (vec_result.size() == dof.n_dofs(), ExcDimensionMismatch (vec_result.size(), dof.n_dofs())); - + // make up boundary values std::map boundary_values; - - if (enforce_zero_boundary == true) + + if (enforce_zero_boundary == true) // no need to project boundary // values, but enforce // homogeneous boundary values // anyway internal::VectorTools:: interpolate_zero_boundary_values (dof, boundary_values); - + else // no homogeneous boundary values if (project_to_boundary_first == true) @@ -500,7 +500,7 @@ void VectorTools::project (const Mapping &mapping, prec.initialize(mass_matrix, 1.2); // solve cg.solve (mass_matrix, vec, tmp, prec); - + // distribute solution constraints.distribute (vec); @@ -517,14 +517,14 @@ void VectorTools::project (const Mapping &mapping, template -void VectorTools::project (const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - VECTOR &vec, - const bool enforce_zero_boundary, - const Quadrature &q_boundary, - const bool project_to_boundary_first) +void VectorTools::project (const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + VECTOR &vec, + const bool enforce_zero_boundary, + const Quadrature &q_boundary, + const bool project_to_boundary_first) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); project(StaticMappingQ1::mapping, dof, constraints, quadrature, function, vec, @@ -535,11 +535,12 @@ void VectorTools::project (const DoFHandler &dof, template -void VectorTools::create_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) +void VectorTools::create_right_hand_side + (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector) { const FiniteElement &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, @@ -547,7 +548,7 @@ void VectorTools::create_right_hand_side (const Mapping &mappi Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); rhs_vector = 0; - + UpdateFlags update_flags = UpdateFlags(update_values | update_quadrature_points | update_JxW_values); @@ -556,7 +557,7 @@ void VectorTools::create_right_hand_side (const Mapping &mappi const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points, n_components = fe.n_components(); - + std::vector dofs (dofs_per_cell); Vector cell_vector (dofs_per_cell); @@ -567,42 +568,42 @@ void VectorTools::create_right_hand_side (const Mapping &mappi if (n_components==1) { std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) + + for (; cell!=endc; ++cell) { fe_values.reinit(cell); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; for (unsigned int point=0; pointget_dof_indices (dofs); - + for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) + + for (; cell!=endc; ++cell) { fe_values.reinit(cell); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; // Use the faster code if the // FiniteElement is primitive @@ -613,7 +614,7 @@ void VectorTools::create_right_hand_side (const Mapping &mappi { const unsigned int component = fe.system_to_component_index(i).first; - + cell_vector(i) += rhs_values[point](component) * fe_values.shape_value(i,point) * weights[point]; @@ -634,9 +635,9 @@ void VectorTools::create_right_hand_side (const Mapping &mappi weights[point]; } } - + cell->get_dof_indices (dofs); - + for (unsigned int i=0; i &mappi template -void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) +void VectorTools::create_right_hand_side + (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_right_hand_side(StaticMappingQ1::mapping, dof_handler, quadrature, - rhs_function, rhs_vector); + create_right_hand_side(StaticMappingQ1::mapping, dof_handler, + quadrature, rhs_function, rhs_vector); } template -void VectorTools::create_right_hand_side (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof_handler, - const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector) +void VectorTools::create_right_hand_side + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector) { const hp::FECollection &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, @@ -672,14 +675,14 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection x_fe_values (mapping, fe, quadrature, update_flags); const unsigned int n_components = fe.n_components(); - + std::vector dofs (fe.max_dofs_per_cell()); Vector cell_vector (fe.max_dofs_per_cell()); @@ -690,60 +693,62 @@ void VectorTools::create_right_hand_side (const hp::MappingCollection rhs_values; - - for (; cell!=endc; ++cell) + + for (; cell!=endc; ++cell) { x_fe_values.reinit(cell); - const FEValues &fe_values = x_fe_values.get_present_fe_values(); - + const FEValues &fe_values = + x_fe_values.get_present_fe_values(); + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points; rhs_values.resize (n_q_points); dofs.resize (dofs_per_cell); cell_vector.reinit (dofs_per_cell); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; for (unsigned int point=0; pointget_dof_indices (dofs); - + for (unsigned int i=0; i > rhs_values; - + for (; cell!=endc; ++cell) { x_fe_values.reinit(cell); - const FEValues &fe_values = x_fe_values.get_present_fe_values(); - + const FEValues &fe_values = + x_fe_values.get_present_fe_values(); + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points; rhs_values.resize (n_q_points, Vector(n_components)); dofs.resize (dofs_per_cell); cell_vector.reinit (dofs_per_cell); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; - + // Use the faster code if the // FiniteElement is primitive if (cell->get_fe().is_primitive ()) @@ -753,7 +758,7 @@ void VectorTools::create_right_hand_side (const hp::MappingCollectionget_fe().system_to_component_index(i).first; - + cell_vector(i) += rhs_values[point](component) * fe_values.shape_value(i,point) * weights[point]; @@ -773,9 +778,9 @@ void VectorTools::create_right_hand_side (const hp::MappingCollectionget_dof_indices (dofs); - + for (unsigned int i=0; i -void VectorTools::create_right_hand_side (const hp::DoFHandler &dof_handler, - const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector) +void VectorTools::create_right_hand_side + (const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - create_right_hand_side(hp::StaticMappingQ1::mapping_collection, + create_right_hand_side(hp::StaticMappingQ1::mapping_collection, dof_handler, quadrature, rhs_function, rhs_vector); } @@ -800,16 +806,17 @@ void VectorTools::create_right_hand_side (const hp::DoFHandler template -void VectorTools::create_point_source_vector (const Mapping &mapping, - const DoFHandler &dof_handler, - const Point &p, - Vector &rhs_vector) +void VectorTools::create_point_source_vector + (const Mapping &mapping, + const DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) { Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); Assert (dof_handler.get_fe().n_components() == 1, ExcMessage ("This function only works for scalar finite elements")); - + rhs_vector = 0; std::pair::active_cell_iterator, Point > @@ -834,9 +841,10 @@ void VectorTools::create_point_source_vector (const Mapping template -void VectorTools::create_point_source_vector (const DoFHandler &dof_handler, - const Point &p, - Vector &rhs_vector) +void VectorTools::create_point_source_vector + (const DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); create_point_source_vector(StaticMappingQ1::mapping, dof_handler, @@ -845,16 +853,17 @@ void VectorTools::create_point_source_vector (const DoFHandler template -void VectorTools::create_point_source_vector (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof_handler, - const Point &p, - Vector &rhs_vector) +void VectorTools::create_point_source_vector + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) { Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); Assert (dof_handler.get_fe().n_components() == 1, ExcMessage ("This function only works for scalar finite elements")); - + rhs_vector = 0; std::pair::active_cell_iterator, Point > @@ -879,9 +888,10 @@ void VectorTools::create_point_source_vector (const hp::MappingCollection -void VectorTools::create_point_source_vector (const hp::DoFHandler &dof_handler, - const Point &p, - Vector &rhs_vector) +void VectorTools::create_point_source_vector + (const hp::DoFHandler &dof_handler, + const Point &p, + Vector &rhs_vector) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); create_point_source_vector(hp::StaticMappingQ1::mapping_collection, @@ -894,21 +904,22 @@ void VectorTools::create_point_source_vector (const hp::DoFHandler template void -VectorTools::create_boundary_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, - const std::set &boundary_indicators) +VectorTools::create_boundary_right_hand_side + (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) { const FiniteElement &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, ExcDimensionMismatch(fe.n_components(), rhs_function.n_components)); Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - + rhs_vector = 0; - + UpdateFlags update_flags = UpdateFlags(update_values | update_quadrature_points | update_JxW_values); @@ -917,17 +928,18 @@ VectorTools::create_boundary_right_hand_side (const Mapping const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = fe_values.n_quadrature_points, n_components = fe.n_components(); - + std::vector dofs (dofs_per_cell); Vector cell_vector (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); if (n_components==1) { std::vector rhs_values(n_q_points); - + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->at_boundary () && @@ -936,19 +948,19 @@ VectorTools::create_boundary_right_hand_side (const Mapping boundary_indicators.end())) { fe_values.reinit(cell, face); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; for (unsigned int point=0; pointget_dof_indices (dofs); - + for (unsigned int i=0; i else { std::vector > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) + + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->at_boundary () && (boundary_indicators.find (cell->face(face)->boundary_indicator()) @@ -965,22 +977,22 @@ VectorTools::create_boundary_right_hand_side (const Mapping boundary_indicators.end())) { fe_values.reinit(cell, face); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; - + // Use the faster code if the // FiniteElement is primitive if (fe.is_primitive ()) - { + { for (unsigned int point=0; point weights[point]; } } - + cell->get_dof_indices (dofs); - + for (unsigned int i=0; i // Implementation for 1D template void -VectorTools::create_boundary_right_hand_side (const Mapping &, +VectorTools::create_boundary_right_hand_side (const Mapping &, const DoFHandler &, - const Quadrature &, - const Function &, - Vector &, - const std::set &) + const Quadrature &, + const Function &, + Vector &, + const std::set &) { Assert (false, ExcImpossibleInDim(dim)); } @@ -1030,11 +1042,12 @@ VectorTools::create_boundary_right_hand_side (const Mapping &, template void -VectorTools::create_boundary_right_hand_side (const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector, - const std::set &boundary_indicators) +VectorTools::create_boundary_right_hand_side + (const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -1050,28 +1063,29 @@ VectorTools::create_boundary_right_hand_side (const DoFHandler & template void -VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof_handler, - const hp::QCollection &quadrature, - const Function &rhs_function, - Vector &rhs_vector, - const std::set &boundary_indicators) +VectorTools::create_boundary_right_hand_side + (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof_handler, + const hp::QCollection &quadrature, + const Function &rhs_function, + Vector &rhs_vector, + const std::set &boundary_indicators) { const hp::FECollection &fe = dof_handler.get_fe(); Assert (fe.n_components() == rhs_function.n_components, ExcDimensionMismatch(fe.n_components(), rhs_function.n_components)); Assert (rhs_vector.size() == dof_handler.n_dofs(), ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - + rhs_vector = 0; - + UpdateFlags update_flags = UpdateFlags(update_values | update_quadrature_points | update_JxW_values); hp::FEFaceValues x_fe_values (mapping, fe, quadrature, update_flags); const unsigned int n_components = fe.n_components(); - + std::vector dofs (fe.max_dofs_per_cell()); Vector cell_vector (fe.max_dofs_per_cell()); @@ -1082,7 +1096,7 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection rhs_values; - + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->at_boundary () && @@ -1097,19 +1111,19 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &weights = fe_values.get_JxW_values (); rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; for (unsigned int point=0; pointget_dof_indices (dofs); - + for (unsigned int i=0; i > rhs_values; - - for (; cell!=endc; ++cell) + + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->at_boundary () && (boundary_indicators.find (cell->face(face)->boundary_indicator()) @@ -1132,22 +1146,22 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollection(n_components)); - + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - + cell_vector = 0; - + // Use the faster code if the // FiniteElement is primitive if (cell->get_fe().is_primitive ()) - { + { for (unsigned int point=0; pointget_fe().system_to_component_index(i).first; - + cell_vector(i) += rhs_values[point](component) * fe_values.shape_value(i,point) * weights[point]; @@ -1169,9 +1183,9 @@ VectorTools::create_boundary_right_hand_side (const hp::MappingCollectionget_dof_indices (dofs); - + for (unsigned int i=0; i void -VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &, +VectorTools::create_boundary_right_hand_side (const hp::MappingCollection &, const hp::DoFHandler &, const hp::QCollection &, const Function &, @@ -1240,7 +1254,7 @@ VectorTools::interpolate_boundary_values (const Mappingneighbor(direction).state() == IteratorState::valid) outermost_cell = outermost_cell->neighbor(direction); - + while (outermost_cell->has_children()) outermost_cell = outermost_cell->child(direction); @@ -1296,7 +1310,7 @@ VectorTools::interpolate_boundary_values (const Mappingvertex(direction), function_values); - + for (unsigned int i=0; ivertex_dof_index(direction,i)] @@ -1353,7 +1367,7 @@ interpolate_boundary_values (const Mapping // immediately if (function_map.size() == 0) return; - + Assert (function_map.find(255) == function_map.end(), ExcInvalidBoundaryIndicator()); @@ -1403,7 +1417,7 @@ interpolate_boundary_values (const Mapping 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 @@ -1460,9 +1474,9 @@ interpolate_boundary_values (const Mapping if (component_mask[fe.face_system_to_component_index(i).first] == true) unit_support_points[i] = fe.unit_face_support_point(i); - + q_collection.push_back (Quadrature(unit_support_points)); - } + } } // now that we have a q_collection // object with all the right @@ -1473,7 +1487,7 @@ interpolate_boundary_values (const Mapping 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) @@ -1505,10 +1519,10 @@ interpolate_boundary_values (const Mapping "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()) + if (function_map.find(boundary_component) != function_map.end()) { // face is of the right component x_fe_values.reinit(cell, face_no); @@ -1521,7 +1535,7 @@ interpolate_boundary_values (const Mapping 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 @@ -1535,10 +1549,10 @@ interpolate_boundary_values (const Mapping 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 @@ -1646,10 +1660,10 @@ interpolate_boundary_values (const Mapping - fe.get_nonzero_components(cell_i).begin()); } - + if (component_mask[component] == true) boundary_values[face_dofs[i]] = dof_values_system[i](component); - } + } } else // fe has only one component, @@ -1660,9 +1674,9 @@ interpolate_boundary_values (const Mapping dof_values_scalar.resize (fe.dofs_per_face); function_map.find(boundary_component)->second ->value_list (dof_locations, dof_values_scalar, 0); - + // enter into list - + for (unsigned int i=0; i void @@ -1699,7 +1713,7 @@ VectorTools::interpolate_boundary_values (const DH &dof, const std::vector &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - interpolate_boundary_values(StaticMappingQ1::mapping, + interpolate_boundary_values(StaticMappingQ1::mapping, dof, boundary_component, boundary_function, boundary_values, component_mask); } @@ -1714,7 +1728,7 @@ VectorTools::interpolate_boundary_values (const DH &dof, const std::vector &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - interpolate_boundary_values(StaticMappingQ1::mapping, + interpolate_boundary_values(StaticMappingQ1::mapping, dof, function_map, boundary_values, component_mask); } @@ -1754,7 +1768,7 @@ VectorTools::interpolate_boundary_values (const Mappingneighbor(direction).state() == IteratorState::valid) outermost_cell = outermost_cell->neighbor(direction); - + while (outermost_cell->has_children()) outermost_cell = outermost_cell->child(direction); @@ -1810,7 +1824,7 @@ VectorTools::interpolate_boundary_values (const Mappingvertex(direction), function_values); - + for (unsigned int i=0; ivertex_dof_index(direction,i); constraints.add_line (row); - constraints.set_inhomogeneity (row, + constraints.set_inhomogeneity (row, function_values(fe.face_system_to_component_index(i).first)); } } @@ -1830,7 +1844,7 @@ VectorTools::interpolate_boundary_values (const Mapping void -VectorTools::interpolate_boundary_values +VectorTools::interpolate_boundary_values (const Mapping &mapping, const DH &dof, const typename FunctionMap::type &function_map, @@ -1853,7 +1867,7 @@ VectorTools::interpolate_boundary_values template void -VectorTools::interpolate_boundary_values +VectorTools::interpolate_boundary_values (const Mapping &mapping, const DH &dof, const typename FunctionMap::type &function_map, @@ -1874,7 +1888,7 @@ VectorTools::interpolate_boundary_values // immediately if (function_map.size() == 0) return; - + Assert (function_map.find(255) == function_map.end(), ExcInvalidBoundaryIndicator()); @@ -1922,7 +1936,7 @@ VectorTools::interpolate_boundary_values 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 @@ -1973,9 +1987,9 @@ VectorTools::interpolate_boundary_values if (component_mask[fe.face_system_to_component_index(i).first] == true) unit_support_points[i] = fe.unit_face_support_point(i); - + q_collection.push_back (Quadrature(unit_support_points)); - } + } } // now that we have a q_collection object // with all the right quadrature points, @@ -1985,7 +1999,7 @@ VectorTools::interpolate_boundary_values 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) @@ -2017,10 +2031,10 @@ VectorTools::interpolate_boundary_values "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()) + if (function_map.find(boundary_component) != function_map.end()) { // face is of the right // component @@ -2034,7 +2048,7 @@ VectorTools::interpolate_boundary_values 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 @@ -2046,10 +2060,10 @@ VectorTools::interpolate_boundary_values 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 @@ -2136,7 +2150,7 @@ VectorTools::interpolate_boundary_values - fe.get_nonzero_components(cell_i).begin()); } - + if (component_mask[component] == true) { // TODO: check whether we should clear @@ -2145,7 +2159,7 @@ VectorTools::interpolate_boundary_values constraints.set_inhomogeneity (face_dofs[i], dof_values_system[i](component)); } - } + } } else // fe has only one component, @@ -2157,9 +2171,9 @@ VectorTools::interpolate_boundary_values dof_values_scalar.resize (fe.dofs_per_face); function_map.find(boundary_component)->second ->value_list (dof_locations, dof_values_scalar, 0); - + // enter into list - + for (unsigned int i=0; i void -VectorTools::interpolate_boundary_values +VectorTools::interpolate_boundary_values (const Mapping &mapping, const DH &dof, const unsigned char boundary_component, @@ -2192,11 +2206,11 @@ VectorTools::interpolate_boundary_values } #endif - + template void -VectorTools::interpolate_boundary_values +VectorTools::interpolate_boundary_values (const DH &dof, const unsigned char boundary_component, const Function &boundary_function, @@ -2204,7 +2218,7 @@ VectorTools::interpolate_boundary_values const std::vector &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); - interpolate_boundary_values(StaticMappingQ1::mapping, + interpolate_boundary_values(StaticMappingQ1::mapping, dof, boundary_component, boundary_function, constraints, component_mask); } @@ -2213,14 +2227,14 @@ VectorTools::interpolate_boundary_values template void -VectorTools::interpolate_boundary_values +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, + interpolate_boundary_values(StaticMappingQ1::mapping, dof, function_map, constraints, component_mask); } @@ -2279,16 +2293,16 @@ VectorTools::project_boundary_values (const Mapping &mappin } 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) @@ -2318,7 +2332,7 @@ VectorTools::project_boundary_values (const Mapping &mappin 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) @@ -2329,7 +2343,7 @@ VectorTools::project_boundary_values (const Mapping &mappin if (level == -1) level = cell->level(); else - { + { Assert (level == cell->level(), ExcNotImplemented()); } } @@ -2337,14 +2351,14 @@ VectorTools::project_boundary_values (const Mapping &mappin #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, + MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, mass_matrix, boundary_functions, rhs, dof_to_boundary_mapping, (const Function*) 0, component_mapping); @@ -2357,7 +2371,7 @@ VectorTools::project_boundary_values (const Mapping &mappin // 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); @@ -2366,7 +2380,7 @@ VectorTools::project_boundary_values (const Mapping &mappin for (unsigned int i=0;i max_element) max_element = mass_matrix.diag_element(i); - + for (unsigned int i=0;i &mappin mass_matrix.diag_element(i) = 1.; excluded_dofs[i] = true; } - + Vector boundary_projection (rhs.size()); // Allow for a maximum of 5*n @@ -2394,7 +2408,7 @@ VectorTools::project_boundary_values (const Mapping &mappin // solve cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition); filtered_precondition.apply_constraints(boundary_projection, true); - filtered_precondition.clear(); + filtered_precondition.clear(); // fill in boundary values for (unsigned int i=0; i::invalid_dof_index @@ -2477,16 +2491,16 @@ VectorTools::project_boundary_values (const Mapping &mappin } 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) @@ -2516,7 +2530,7 @@ VectorTools::project_boundary_values (const Mapping &mappin 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) @@ -2527,7 +2541,7 @@ VectorTools::project_boundary_values (const Mapping &mappin if (level == -1) level = cell->level(); else - { + { Assert (level == cell->level(), ExcNotImplemented()); } } @@ -2535,14 +2549,14 @@ VectorTools::project_boundary_values (const Mapping &mappin #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, + MatrixCreator::create_boundary_mass_matrix (mapping, dof, q, mass_matrix, boundary_functions, rhs, dof_to_boundary_mapping, (const Function*) 0, component_mapping); @@ -2555,7 +2569,7 @@ VectorTools::project_boundary_values (const Mapping &mappin // 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); @@ -2564,7 +2578,7 @@ VectorTools::project_boundary_values (const Mapping &mappin for (unsigned int i=0;i max_element) max_element = mass_matrix.diag_element(i); - + for (unsigned int i=0;i &mappin mass_matrix.diag_element(i) = 1.; excluded_dofs[i] = true; } - + Vector boundary_projection (rhs.size()); // Allow for a maximum of 5*n @@ -2592,7 +2606,7 @@ VectorTools::project_boundary_values (const Mapping &mappin // solve cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition); filtered_precondition.apply_constraints(boundary_projection, true); - filtered_precondition.clear(); + filtered_precondition.clear(); // fill in boundary values for (unsigned int i=0; i::invalid_dof_index @@ -2650,8 +2664,8 @@ namespace internal for (unsigned int i=0; i &other) const { for (unsigned int i=0; i &dof_indices, const Tensor<1,dim> &constraining_vector, ConstraintMatrix &constraints) - { + { // choose the DoF that has the // largest component in the // constraining_vector as the @@ -2748,7 +2762,7 @@ namespace internal (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2]))) { constraints.add_line (dof_indices.dof_indices[0]); - + if (std::fabs (constraining_vector[1]/constraining_vector[0]) > std::numeric_limits::epsilon()) constraints.add_entry (dof_indices.dof_indices[0], @@ -2796,14 +2810,14 @@ namespace internal dof_indices.dof_indices[1], -constraining_vector[1]/constraining_vector[2]); } - + break; } default: Assert (false, ExcNotImplemented()); } - + } @@ -2819,7 +2833,7 @@ namespace internal void compute_orthonormal_vectors (const Tensor<1,dim> &vector, Tensor<1,dim> (&orthonormals)[dim-1]) - { + { switch (dim) { case 3: @@ -2906,9 +2920,9 @@ namespace internal default: Assert (false, ExcNotImplemented()); - } + } } - + } } @@ -2926,7 +2940,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH ExcMessage ("This function is not useful in 1d because it amounts " "to imposing Dirichlet values on the vector-valued " "quantity.")); - + const FiniteElement &fe = dof_handler.get_fe(); std::vector face_dofs (fe.dofs_per_face); @@ -2949,7 +2963,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH // of the same cell. consequently, we also // have to store which cell a normal vector // was computed on - typedef + typedef std::multimap, std::pair, typename DH::active_cell_iterator> > DoFToNormalsMap; @@ -2976,7 +2990,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH // vectors at the locations // where they are defined: const std::vector > & - unit_support_points = fe.get_unit_face_support_points(); + unit_support_points = fe.get_unit_face_support_points(); Quadrature aux_quad (unit_support_points); FEFaceValues fe_values (mapping, fe, aux_quad, update_normal_vectors); @@ -2994,7 +3008,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH // components of vector internal::VectorTools::VectorDoFTuple vector_dofs; vector_dofs.dof_indices[0] = face_dofs[i]; - + for (unsigned int k=0; k vector_dofs.dof_indices[fe.face_system_to_component_index(k).first - first_vector_component] = face_dofs[k]; - + for (unsigned int d=0; d // dofs typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin(); - + while (p != dof_to_normals_map.end()) { // first find the range of entries in @@ -3128,7 +3142,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH // average over all normal vectors case 1: { - + // compute the average // normal vector from all // the ones that have the @@ -3145,12 +3159,12 @@ VectorTools::compute_no_normal_flux_constraints (const DH x = cell_to_normals_map.begin(); x != cell_to_normals_map.end(); ++x) normal += x->second.first; - normal /= normal.norm(); - + normal /= normal.norm(); + // then construct constraints // from this: const internal::VectorTools::VectorDoFTuple & - dof_indices = same_dof_range[0]->first; + dof_indices = same_dof_range[0]->first; internal::VectorTools::add_constraint (dof_indices, normal, constraints); @@ -3215,11 +3229,11 @@ VectorTools::compute_no_normal_flux_constraints (const DH for (unsigned int i=0; isecond.first[j]; - + Assert (std::fabs(determinant (t)) > 1e-3, ExcMessage("Found a set of normal vectors that are nearly collinear.")); } - + // so all components of // this vector dof are // constrained. enter @@ -3233,7 +3247,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH break; } - + // this is the case of an // edge contribution in 3d, @@ -3270,7 +3284,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH q != same_dof_range[1]; ++q) cell_contributions[q->second.second].push_back (q->second.first); Assert (cell_contributions.size() >= 1, ExcInternalError()); - + // now for each cell that // has contributed // determine the number @@ -3369,7 +3383,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH Assert (std::fabs (tangent.norm()-1) < 1e-12, ExcInternalError()); - + tangential_vectors.push_back (tangent); } @@ -3392,7 +3406,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH if (*t * first_tangent < 0) *t *= -1; } - + // now compute the // average tangent and // normalize it @@ -3457,7 +3471,7 @@ namespace internal // necessary, so have a read-write // version of it: double exponent = exponent_1; - + const unsigned int n_components = dof.get_fe().n_components(); const bool fe_is_system = (n_components != 1); @@ -3468,7 +3482,7 @@ namespace internal } difference.reinit (dof.get_tria().n_active_cells()); - + switch (norm) { case dealii::VectorTools::L2_norm: @@ -3482,7 +3496,7 @@ namespace internal default: break; } - + UpdateFlags update_flags = UpdateFlags (update_quadrature_points | update_JxW_values); switch (norm) @@ -3500,13 +3514,13 @@ namespace internal default: update_flags |= UpdateFlags (update_values); break; - } + } dealii::hp::FECollection fe_collection (dof.get_fe()); dealii::hp::FEValues x_fe_values(mapping, fe_collection, q, update_flags); const unsigned int max_n_q_points = q.max_n_quadrature_points (); - + std::vector< dealii::Vector > function_values (max_n_q_points, dealii::Vector(n_components)); std::vector > > @@ -3529,7 +3543,7 @@ namespace internal // scalar functions std::vector tmp_values (max_n_q_points); std::vector > tmp_gradients (max_n_q_points); - + // loop over all cells typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); @@ -3541,7 +3555,7 @@ namespace internal const dealii::FEValues &fe_values = x_fe_values.get_present_fe_values (); const unsigned int n_q_points = fe_values.n_quadrature_points; - + // resize all out scratch // arrays to the number of // quadrature points we use @@ -3582,8 +3596,8 @@ namespace internal for (unsigned int k=0;k &mapping, Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); - + const Quadrature quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); FEValues fe_values(mapping, fe, quadrature, update_values); @@ -3922,7 +3936,7 @@ VectorTools::point_difference (const Mapping &mapping, difference(0) = exact_function.value(point); else exact_function.vector_value(point, difference); - + for (unsigned int i=0; i &mapping, Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, ExcInternalError()); - + const Quadrature quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); FEValues fe_values(mapping, fe, quadrature, update_values); @@ -4043,7 +4057,7 @@ VectorTools::compute_mean_value (const Mapping &mapping, ExcDimensionMismatch (v.size(), dof.n_dofs())); Assert (component < dof.get_fe().n_components(), ExcIndexRange(component, 0, dof.get_fe().n_components())); - + FEValues fe(mapping, dof.get_fe(), quadrature, UpdateFlags(update_JxW_values | update_values)); @@ -4051,7 +4065,7 @@ VectorTools::compute_mean_value (const Mapping &mapping, typename DoFHandler::active_cell_iterator c; std::vector > values(quadrature.size(), Vector (dof.get_fe().n_components())); - + double mean = 0.; double area = 0.; // Compute mean value @@ -4065,7 +4079,7 @@ VectorTools::compute_mean_value (const Mapping &mapping, area += fe.JxW(k); }; }; - + return (mean/area); } diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 58d8142dc2..5b724676d0 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -2018,6 +2018,7 @@ void MatrixCreator::create_mass_matrix const Function * const coefficient); +#if deal_II_dimension != 1 template void MatrixCreator::create_boundary_mass_matrix @@ -2055,7 +2056,6 @@ MatrixCreator::create_boundary_mass_matrix const Function * const, std::vector); - template void MatrixCreator::create_boundary_mass_matrix (const hp::DoFHandler&, @@ -2067,6 +2067,7 @@ void MatrixCreator::create_boundary_mass_matrix const Function * const, std::vector); +#endif // non-hp versions of create_laplace_matrix @@ -2201,6 +2202,8 @@ void MatrixCreator::create_mass_matrix &rhs_vector, const Function * const coefficient); +#if deal_II_dimension == 2 + template void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, @@ -2224,7 +2227,7 @@ void MatrixCreator::create_boundary_mass_matrix * const a, std::vector); - +#endif // #if deal_II_dimension != 1 // template