From: Martin Kronbichler Date: Fri, 13 Nov 2009 22:07:05 +0000 (+0000) Subject: Clean up parts of the code. X-Git-Tag: v8.0.0~6814 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e58b276c9da8da1efa6776fe2684dc4f0a1c9438;p=dealii.git Clean up parts of the code. git-svn-id: https://svn.dealii.org/trunk@20112 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 3674578717..266c76ee1f 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,12 +535,11 @@ 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, @@ -548,7 +547,7 @@ void VectorTools::create_right_hand_side 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); @@ -557,7 +556,7 @@ void VectorTools::create_right_hand_side 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); @@ -568,42 +567,42 @@ void VectorTools::create_right_hand_side 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 @@ -614,7 +613,7 @@ void VectorTools::create_right_hand_side { 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]; @@ -635,9 +634,9 @@ void VectorTools::create_right_hand_side weights[point]; } } - + cell->get_dof_indices (dofs); - + for (unsigned int i=0; i -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, @@ -675,14 +672,14 @@ void VectorTools::create_right_hand_side 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::FEValues 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()); @@ -693,62 +690,60 @@ void VectorTools::create_right_hand_side if (n_components==1) { std::vector 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 ()) @@ -758,7 +753,7 @@ void VectorTools::create_right_hand_side { const unsigned int component = cell->get_fe().system_to_component_index(i).first; - + cell_vector(i) += rhs_values[point](component) * fe_values.shape_value(i,point) * weights[point]; @@ -778,9 +773,9 @@ void VectorTools::create_right_hand_side weights[point]; } } - + cell->get_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); } @@ -806,17 +800,16 @@ void VectorTools::create_right_hand_side 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 > @@ -841,10 +834,9 @@ void VectorTools::create_point_source_vector 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, @@ -853,17 +845,16 @@ void VectorTools::create_point_source_vector 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 > @@ -888,10 +879,9 @@ void VectorTools::create_point_source_vector template -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, @@ -904,22 +894,21 @@ void VectorTools::create_point_source_vector 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); @@ -928,18 +917,17 @@ VectorTools::create_boundary_right_hand_side 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 () && @@ -948,19 +936,19 @@ VectorTools::create_boundary_right_hand_side 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 > 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()) @@ -977,22 +965,22 @@ VectorTools::create_boundary_right_hand_side 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; pointget_dof_indices (dofs); - + for (unsigned int i=0; i 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)); } @@ -1042,12 +1030,11 @@ 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")); @@ -1063,29 +1050,28 @@ VectorTools::create_boundary_right_hand_side 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()); @@ -1096,7 +1082,7 @@ VectorTools::create_boundary_right_hand_side if (n_components==1) { std::vector rhs_values; - + for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->at_boundary () && @@ -1111,19 +1097,19 @@ VectorTools::create_boundary_right_hand_side 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); - + 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) + + 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()) @@ -1146,22 +1132,22 @@ VectorTools::create_boundary_right_hand_side 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)); - + 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]; @@ -1183,9 +1169,9 @@ VectorTools::create_boundary_right_hand_side weights[point]; } } - + cell->get_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 &, @@ -1254,7 +1240,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); @@ -1310,7 +1296,7 @@ VectorTools::interpolate_boundary_values (const Mappingvertex(direction), function_values); - + for (unsigned int i=0; ivertex_dof_index(direction,i)] @@ -1367,7 +1353,7 @@ interpolate_boundary_values (const Mapping // immediately if (function_map.size() == 0) return; - + Assert (function_map.find(255) == function_map.end(), ExcInvalidBoundaryIndicator()); @@ -1417,7 +1403,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 @@ -1474,9 +1460,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 @@ -1487,7 +1473,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) @@ -1519,10 +1505,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); @@ -1535,7 +1521,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 @@ -1549,10 +1535,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 @@ -1660,10 +1646,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, @@ -1674,9 +1660,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 @@ -1713,7 +1699,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); } @@ -1728,7 +1714,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); } @@ -1768,7 +1754,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); @@ -1824,7 +1810,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)); } } @@ -1844,7 +1830,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, @@ -1867,7 +1853,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, @@ -1888,7 +1874,7 @@ VectorTools::interpolate_boundary_values // immediately if (function_map.size() == 0) return; - + Assert (function_map.find(255) == function_map.end(), ExcInvalidBoundaryIndicator()); @@ -1936,7 +1922,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 @@ -1987,9 +1973,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, @@ -1999,7 +1985,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) @@ -2031,10 +2017,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 @@ -2048,7 +2034,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 @@ -2060,10 +2046,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 @@ -2150,7 +2136,7 @@ VectorTools::interpolate_boundary_values - fe.get_nonzero_components(cell_i).begin()); } - + if (component_mask[component] == true) { // TODO: check whether we should clear @@ -2159,7 +2145,7 @@ VectorTools::interpolate_boundary_values constraints.set_inhomogeneity (face_dofs[i], dof_values_system[i](component)); } - } + } } else // fe has only one component, @@ -2171,9 +2157,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, @@ -2206,11 +2192,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, @@ -2218,7 +2204,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); } @@ -2227,14 +2213,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); } @@ -2293,16 +2279,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) @@ -2332,7 +2318,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) @@ -2343,7 +2329,7 @@ VectorTools::project_boundary_values (const Mapping &mappin if (level == -1) level = cell->level(); else - { + { Assert (level == cell->level(), ExcNotImplemented()); } } @@ -2351,14 +2337,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); @@ -2371,7 +2357,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); @@ -2380,7 +2366,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 @@ -2408,7 +2394,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 @@ -2491,16 +2477,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) @@ -2530,7 +2516,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) @@ -2541,7 +2527,7 @@ VectorTools::project_boundary_values (const Mapping &mappin if (level == -1) level = cell->level(); else - { + { Assert (level == cell->level(), ExcNotImplemented()); } } @@ -2549,14 +2535,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); @@ -2569,7 +2555,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); @@ -2578,7 +2564,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 @@ -2606,7 +2592,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 @@ -2664,8 +2650,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 @@ -2762,7 +2748,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], @@ -2810,14 +2796,14 @@ namespace internal dof_indices.dof_indices[1], -constraining_vector[1]/constraining_vector[2]); } - + break; } default: Assert (false, ExcNotImplemented()); } - + } @@ -2833,7 +2819,7 @@ namespace internal void compute_orthonormal_vectors (const Tensor<1,dim> &vector, Tensor<1,dim> (&orthonormals)[dim-1]) - { + { switch (dim) { case 3: @@ -2920,9 +2906,9 @@ namespace internal default: Assert (false, ExcNotImplemented()); - } + } } - + } } @@ -2940,7 +2926,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); @@ -2963,7 +2949,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; @@ -2990,7 +2976,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); @@ -3008,7 +2994,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 @@ -3142,7 +3128,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 @@ -3159,12 +3145,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); @@ -3229,11 +3215,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 @@ -3247,7 +3233,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH break; } - + // this is the case of an // edge contribution in 3d, @@ -3284,7 +3270,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 @@ -3383,7 +3369,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH Assert (std::fabs (tangent.norm()-1) < 1e-12, ExcInternalError()); - + tangential_vectors.push_back (tangent); } @@ -3406,7 +3392,7 @@ VectorTools::compute_no_normal_flux_constraints (const DH if (*t * first_tangent < 0) *t *= -1; } - + // now compute the // average tangent and // normalize it @@ -3471,7 +3457,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); @@ -3482,7 +3468,7 @@ namespace internal } difference.reinit (dof.get_tria().n_active_cells()); - + switch (norm) { case dealii::VectorTools::L2_norm: @@ -3496,7 +3482,7 @@ namespace internal default: break; } - + UpdateFlags update_flags = UpdateFlags (update_quadrature_points | update_JxW_values); switch (norm) @@ -3514,13 +3500,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 > > @@ -3543,7 +3529,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(); @@ -3555,7 +3541,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 @@ -3596,8 +3582,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); @@ -3936,7 +3922,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); @@ -4057,7 +4043,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)); @@ -4065,7 +4051,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 @@ -4079,7 +4065,7 @@ VectorTools::compute_mean_value (const Mapping &mapping, area += fe.JxW(k); }; }; - + return (mean/area); }