From: Martin Kronbichler Date: Mon, 7 Sep 2009 19:52:00 +0000 (+0000) Subject: Found a smarter way to set the diagonal. X-Git-Tag: v8.0.0~7164 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d0d29cce460286b177094160f991ff86b8a302ef;p=dealii.git Found a smarter way to set the diagonal. git-svn-id: https://svn.dealii.org/trunk@19409 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc index b146710be2..9cee8280a9 100644 --- a/deal.II/examples/step-37/step-37.cc +++ b/deal.II/examples/step-37/step-37.cc @@ -142,42 +142,42 @@ void Coefficient::value_list (const std::vector > &points, template class MatrixFree : public Subscriptor { -public: - MatrixFree (); - - void reinit (const unsigned int n_dofs, - const unsigned int n_cells, - const FullMatrix &cell_matrix, - const unsigned int n_points_per_cell); - void clear(); - - unsigned int m () const; - unsigned int n () const; - ConstraintMatrix & get_constraints (); - - void set_local_dof_indices (const unsigned int cell_no, - const std::vector &local_dof_indices); - void set_derivative_data (const unsigned int cell_no, - const unsigned int quad_point, - const Transformation &trans_in); - - template - void vmult (Vector &dst, - const Vector &src) const; - template - void Tvmult (Vector &dst, - const Vector &src) const; - template - void vmult_add (Vector &dst, - const Vector &src) const; - template - void Tvmult_add (Vector &dst, - const Vector &src) const; - - number el (const unsigned int row, const unsigned int col) const; - void calculate_diagonal () const; - - std::size_t memory_consumption () const; + public: + MatrixFree (); + + void reinit (const unsigned int n_dofs, + const unsigned int n_cells, + const FullMatrix &cell_matrix, + const unsigned int n_points_per_cell); + void clear(); + + unsigned int m () const; + unsigned int n () const; + ConstraintMatrix & get_constraints (); + + void set_local_dof_indices (const unsigned int cell_no, + const std::vector &local_dof_indices); + void set_derivative_data (const unsigned int cell_no, + const unsigned int quad_point, + const Transformation &trans_in); + + template + void vmult (Vector &dst, + const Vector &src) const; + template + void Tvmult (Vector &dst, + const Vector &src) const; + template + void vmult_add (Vector &dst, + const Vector &src) const; + template + void Tvmult_add (Vector &dst, + const Vector &src) const; + + number el (const unsigned int row, const unsigned int col) const; + void calculate_diagonal () const; + + std::size_t memory_consumption () const; // The private member variables of the // MatrixFree class are a @@ -195,28 +195,28 @@ public: // matrix for handling boundary conditions // as well as a few other variables that // store matrix properties. -private: - template - void vmult_on_subrange (const unsigned int first_cell, - const unsigned int last_cell, - Vector &dst, - const Vector &src) const; - - FullMatrix small_matrix; - Table<2,unsigned int> indices_local_to_global; - Table<2,Transformation> derivatives; - - ConstraintMatrix constraints; - - mutable Vector diagonal_values; - mutable bool diagonal_is_calculated; - - struct MatrixSizes - { - unsigned int n_dofs, n_cells; - unsigned int m, n; - unsigned int n_points, n_comp; - } matrix_sizes; + private: + template + void vmult_on_subrange (const unsigned int first_cell, + const unsigned int last_cell, + Vector &dst, + const Vector &src) const; + + FullMatrix small_matrix; + Table<2,unsigned int> indices_local_to_global; + Table<2,Transformation> derivatives; + + ConstraintMatrix constraints; + + mutable Vector diagonal_values; + mutable bool diagonal_is_calculated; + + struct MatrixSizes + { + unsigned int n_dofs, n_cells; + unsigned int m, n; + unsigned int n_points, n_comp; + } matrix_sizes; }; @@ -226,8 +226,8 @@ private: // nothing. template MatrixFree::MatrixFree () -: - Subscriptor() + : + Subscriptor() {} @@ -671,20 +671,22 @@ MatrixFree::vmult_add (Vector &dst, constraints.condense (dst); // One thing to be cautious about: The - // deal.II classes expect that the matrix - // still contains a diagonal entry for - // constrained dofs (otherwise, the matrix - // would be singular, which is not what we - // want). Since the condense - // command of the constraint matrix sets - // those constrained elements to zero, we - // have to circumvent that problem by using - // the diagonal element which we have - // access to together with the solution - // function. + // deal.II classes expect that the + // matrix still contains a diagonal + // entry for constrained dofs + // (otherwise, the matrix would be + // singular, which is not what we + // want). Since the + // condense command of the + // constraint matrix sets those + // constrained elements to zero, we + // have to circumvent that problem by + // setting the diagonal to some + // non-zero value. We simply set it to + // one. for (unsigned int i=0; i::calculate_diagonal() const diag_value += calculation[q] * small_matrix(dof,q); diagonal_values(indices_local_to_global(cell,dof)) += diag_value; } + constraints.condense (diagonal_values); + for (unsigned int i=0; i::transform (number* result) const if (dim == 2) { const number temp = result[0]; - result[0] = transformation[0] * temp + transformation[1]*result[1]; - result[1] = transformation[1] * temp + transformation[2]*result[1]; + result[0] = transformation[0] * temp + transformation[1] * result[1]; + result[1] = transformation[1] * temp + transformation[2] * result[1]; } else if (dim == 3) { const number temp1 = result[0]; const number temp2 = result[1]; - result[0] = transformation[0] * temp1 + transformation[1]*temp2 + + result[0] = transformation[0] * temp1 + transformation[1] * temp2 + transformation[2] * result[2]; - result[1] = transformation[1] * temp1 + transformation[3]*temp2 + + result[1] = transformation[1] * temp1 + transformation[3] * temp2 + transformation[4] * result[2]; - result[2] = transformation[2] * temp1 + transformation[4]*temp2 + + result[2] = transformation[2] * temp1 + transformation[4] * temp2 + transformation[5] * result[2]; } else @@ -1413,7 +1420,7 @@ void LaplaceProblem::run () int main () { deallog.depth_console (0); - LaplaceProblem<2> laplace_problem (2); + LaplaceProblem<3> laplace_problem (2); laplace_problem.run (); return 0;