From: Wolfgang Bangerth Date: Mon, 25 May 1998 14:40:59 +0000 (+0000) Subject: Add exact mass matrix creation and test for many more internal errors. X-Git-Tag: v8.0.0~22900 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b1bd72cd5aea9099e3acbe00ab0b19a27f9e3afc;p=dealii.git Add exact mass matrix creation and test for many more internal errors. git-svn-id: https://svn.dealii.org/trunk@354 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 ad1723e1ef..f13b4a666c 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -78,6 +78,33 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, +template +void MatrixCreator::create_mass_matrix (const DoFHandler &dof, + const FiniteElement &fe, + const Boundary &boundary, + dSMatrix &matrix) { + const unsigned int total_dofs = fe.total_dofs; + + dFMatrix local_mass_matrix (total_dofs, total_dofs); + vector dofs_on_this_cell (total_dofs); + + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (dofs_on_this_cell); + fe.get_local_mass_matrix (cell, boundary, local_mass_matrix); + + for (unsigned int i=0; i void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, const FiniteElement<1> &, @@ -313,7 +340,12 @@ void MatrixTools::apply_boundary_values (const map &boundary_va dSMatrix &matrix, dVector &solution, dVector &right_hand_side) { - + // if no boundary values are to be applied + // simply return + if (boundary_values.size() == 0) + return; + + map::const_iterator dof = boundary_values.begin(), endd = boundary_values.end(); const unsigned int n_dofs = matrix.m(); @@ -419,11 +451,18 @@ template void MassMatrix::assemble (dFMatrix &cell_matrix, const FEValues &fe_values, const typename Triangulation::cell_iterator &) const { + const unsigned int total_dofs = fe_values.total_dofs, + n_q_points = fe_values.n_quadrature_points; + + Assert (cell_matrix.n() == total_dofs, + ExcWrongSize(cell_matrix.n(), total_dofs)); + Assert (cell_matrix.m() == total_dofs, + ExcWrongSize(cell_matrix.m(), total_dofs)); + Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); + const dFMatrix &values = fe_values.get_shape_values (); const vector &weights = fe_values.get_JxW_values (); - const unsigned int total_dofs = fe_values.total_dofs, - n_q_points = fe_values.n_quadrature_points; if (coefficient != 0) { @@ -455,15 +494,24 @@ void MassMatrix::assemble (dFMatrix &cell_matrix, const FEValues &fe_values, const Triangulation::cell_iterator &) const { Assert (right_hand_side != 0, ExcNoRHSSelected()); - + + const unsigned int total_dofs = fe_values.total_dofs, + n_q_points = fe_values.n_quadrature_points; + + Assert (cell_matrix.n() == total_dofs, + ExcWrongSize(cell_matrix.n(), total_dofs)); + Assert (cell_matrix.m() == total_dofs, + ExcWrongSize(cell_matrix.m(), total_dofs)); + Assert (rhs.size() == total_dofs, + ExcWrongSize(rhs.size(), total_dofs)); + Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); + Assert (rhs.all_zero(), ExcObjectNotEmpty()); + const dFMatrix &values = fe_values.get_shape_values (); const vector &weights = fe_values.get_JxW_values (); vector rhs_values (fe_values.n_quadrature_points); right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - const unsigned int total_dofs = fe_values.total_dofs, - n_q_points = fe_values.n_quadrature_points; - if (coefficient != 0) { vector coefficient_values (n_q_points); @@ -503,15 +551,18 @@ void MassMatrix::assemble (dVector &rhs, const FEValues &fe_values, const Triangulation::cell_iterator &) const { Assert (right_hand_side != 0, ExcNoRHSSelected()); - + + const unsigned int total_dofs = fe_values.total_dofs, + n_q_points = fe_values.n_quadrature_points; + + Assert (rhs.size() == total_dofs, ExcWrongSize(rhs.size(), total_dofs)); + Assert (rhs.all_zero(), ExcObjectNotEmpty()); + const dFMatrix &values = fe_values.get_shape_values (); const vector &weights = fe_values.get_JxW_values (); vector rhs_values(fe_values.n_quadrature_points); right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - const unsigned int total_dofs = fe_values.total_dofs, - n_q_points = fe_values.n_quadrature_points; - for (unsigned int point=0; point::assemble (dFMatrix &cell_matrix, const Triangulation::cell_iterator &) const { Assert (right_hand_side != 0, ExcNoRHSSelected()); + const unsigned int total_dofs = fe_values.total_dofs, + n_q_points = fe_values.n_quadrature_points; + + Assert (cell_matrix.n() == total_dofs, + ExcWrongSize(cell_matrix.n(), total_dofs)); + Assert (cell_matrix.m() == total_dofs, + ExcWrongSize(cell_matrix.m(), total_dofs)); + Assert (rhs.size() == total_dofs, + ExcWrongSize(rhs.size(), total_dofs)); + Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); + Assert (rhs.all_zero(), ExcObjectNotEmpty()); + const vector > >&gradients = fe_values.get_shape_grads (); const dFMatrix &values = fe_values.get_shape_values (); vector rhs_values(fe_values.n_quadrature_points); const vector &weights = fe_values.get_JxW_values (); right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - const unsigned int total_dofs = fe_values.total_dofs, - n_q_points = fe_values.n_quadrature_points; - if (coefficient != 0) { vector coefficient_values(n_q_points); @@ -586,12 +646,18 @@ template void LaplaceMatrix::assemble (dFMatrix &cell_matrix, const FEValues &fe_values, const Triangulation::cell_iterator &) const { - const vector > >&gradients = fe_values.get_shape_grads (); - const vector &weights = fe_values.get_JxW_values (); - const unsigned int total_dofs = fe_values.total_dofs, n_q_points = fe_values.n_quadrature_points; + Assert (cell_matrix.n() == total_dofs, + ExcWrongSize(cell_matrix.n(), total_dofs)); + Assert (cell_matrix.m() == total_dofs, + ExcWrongSize(cell_matrix.m(), total_dofs)); + Assert (cell_matrix.all_zero(), ExcObjectNotEmpty()); + + const vector > >&gradients = fe_values.get_shape_grads (); + const vector &weights = fe_values.get_JxW_values (); + if (coefficient != 0) { vector coefficient_values(n_q_points); @@ -621,15 +687,18 @@ void LaplaceMatrix::assemble (dVector &rhs, const FEValues &fe_values, const Triangulation::cell_iterator &) const { Assert (right_hand_side != 0, ExcNoRHSSelected()); - + + const unsigned int total_dofs = fe_values.total_dofs, + n_q_points = fe_values.n_quadrature_points; + + Assert (rhs.size() == total_dofs, ExcWrongSize(rhs.size(), total_dofs)); + Assert (rhs.all_zero(), ExcObjectNotEmpty()); + const dFMatrix &values = fe_values.get_shape_values (); const vector &weights = fe_values.get_JxW_values (); vector rhs_values(fe_values.n_quadrature_points); right_hand_side->value_list (fe_values.get_quadrature_points(), rhs_values); - const unsigned int total_dofs = fe_values.total_dofs, - n_q_points = fe_values.n_quadrature_points; - for (unsigned int point=0; point