From: bangerth Date: Thu, 28 Dec 2006 16:40:36 +0000 (+0000) Subject: Check in forgotten function. hp/matrices now works and produces exactly the same... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=381e6c2a70a5ebca4f0e8266816ed602208b82c6;p=dealii-svn.git Check in forgotten function. hp/matrices now works and produces exactly the same output as deal.II/matrices. git-svn-id: https://svn.dealii.org/trunk@14282 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index 1a38085eb7..6959df271a 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -1000,7 +1000,6 @@ create_boundary_mass_matrix_1 (const Mapping &mapping, Assert (matrix.n() == matrix.m(), ExcInternalError()); Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof.get_fe() == fe, ExcInternalError()); Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError()); Assert (n_components == boundary_functions.begin()->second->n_components, @@ -1437,10 +1436,9 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, const Function * const coefficient, const IteratorRange > range, Threads::ThreadMutex &mutex) -{ -/* - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); +{ + const hp::FECollection &fe_collection = dof.get_fe(); + const unsigned int n_components = fe_collection.n_components(); const bool fe_is_system = (n_components != 1); Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions), @@ -1448,7 +1446,6 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, Assert (matrix.n() == matrix.m(), ExcInternalError()); Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); Assert (boundary_functions.size() != 0, ExcInternalError()); - Assert (dof.get_fe() == fe, ExcInternalError()); Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError()); Assert (n_components == boundary_functions.begin()->second->n_components, @@ -1469,35 +1466,34 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, }; #endif - const unsigned int dofs_per_cell = fe.dofs_per_cell, - dofs_per_face = fe.dofs_per_face; + const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(), + max_dofs_per_face = fe_collection.max_dofs_per_face(); - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_vector(dofs_per_cell); + FullMatrix cell_matrix(max_dofs_per_cell, max_dofs_per_cell); + Vector cell_vector(max_dofs_per_cell); UpdateFlags update_flags = UpdateFlags (update_values | update_JxW_values | update_q_points); - FEFaceValues fe_values (mapping, fe, q, update_flags); + hp::FEFaceValues x_fe_values (mapping, fe_collection, q, update_flags); // two variables for the coefficient, // one for the two cases indicated in // the name - std::vector coefficient_values (fe_values.n_quadrature_points); - std::vector > coefficient_vector_values (fe_values.n_quadrature_points, - Vector(n_components)); + std::vector coefficient_values; + std::vector > coefficient_vector_values; - std::vector rhs_values_scalar (fe_values.n_quadrature_points); - std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_components)); + std::vector rhs_values_scalar; + std::vector > rhs_values_system; - std::vector dofs (dofs_per_cell); - std::vector dofs_on_face_vector (dofs_per_face); + std::vector dofs (max_dofs_per_cell); + std::vector dofs_on_face_vector (max_dofs_per_face); // for each dof on the cell, have a // flag whether it is on the face - std::vector dof_is_on_face(dofs_per_cell); + std::vector dof_is_on_face(max_dofs_per_cell); + typename hp::DoFHandler::active_cell_iterator cell = range.first; for (; cell!=range.second; ++cell) @@ -1507,14 +1503,24 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, if (boundary_functions.find(cell->face(face)->boundary_indicator()) != boundary_functions.end()) { - cell_matrix = 0; - cell_vector = 0; + x_fe_values.reinit (cell, face); + + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values (); + + const FiniteElement &fe = cell->get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int dofs_per_face = fe.dofs_per_face; - fe_values.reinit (cell, face); + cell_matrix.reinit (dofs_per_cell, dofs_per_cell); + cell_vector.reinit (dofs_per_cell); + cell_matrix = 0; + cell_vector = 0; if (fe_is_system) // FE has several components { + rhs_values_system.resize (fe_values.n_quadrature_points, + Vector(n_components)); boundary_functions.find(cell->face(face)->boundary_indicator()) ->second->vector_value_list (fe_values.get_quadrature_points(), rhs_values_system); @@ -1523,6 +1529,7 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, { if (coefficient->n_components==1) { + coefficient_values.resize (fe_values.n_quadrature_points); coefficient->value_list (fe_values.get_quadrature_points(), coefficient_values); for (unsigned int point=0; point &mapping, } else { + coefficient_vector_values.resize (fe_values.n_quadrature_points, + Vector(n_components)); coefficient->vector_value_list (fe_values.get_quadrature_points(), coefficient_vector_values); for (unsigned int point=0; point &mapping, else // FE is a scalar one { + rhs_values_scalar.resize (fe_values.n_quadrature_points); boundary_functions.find(cell->face(face)->boundary_indicator()) ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); if (coefficient != 0) { + coefficient_values.resize (fe_values.n_quadrature_points); coefficient->value_list (fe_values.get_quadrature_points(), coefficient_values); for (unsigned int point=0; point &mapping, // inefficient, so we copy the dofs // into a set, which enables binary // searches. + dofs.resize (dofs_per_cell); + dofs_on_face_vector.resize (dofs_per_face); + dof_is_on_face.resize (dofs_per_cell); + cell->get_dof_indices (dofs); - cell->face(face)->get_dof_indices (dofs_on_face_vector); + cell->face(face)->get_dof_indices (dofs_on_face_vector, + cell->active_fe_index()); // check for each of the // dofs on this cell @@ -1781,9 +1797,8 @@ create_boundary_mass_matrix_1 (const hp::MappingCollection &mapping, // smallness Assert (std::fabs(cell_vector(j)) <= 1e-10 * max_diag_entry, ExcInternalError()); - }; - }; -*/ + } + } }