From c83e668daa7b62c98ee47cbf8d4b995a2b2a316c Mon Sep 17 00:00:00 2001 From: heister Date: Tue, 5 Feb 2013 16:53:54 +0000 Subject: [PATCH] make indent git-svn-id: https://svn.dealii.org/trunk@28237 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-36/step-36.cc | 2 +- deal.II/examples/step-49/step-49.cc | 68 +-- deal.II/include/deal.II/base/utilities.h | 2 +- deal.II/include/deal.II/fe/fe_q_dg0.h | 2 +- deal.II/include/deal.II/grid/grid_generator.h | 44 +- deal.II/include/deal.II/lac/block_indices.h | 2 +- .../include/deal.II/lac/block_vector_base.h | 4 +- .../deal.II/lac/constraint_matrix.templates.h | 26 +- deal.II/include/deal.II/lac/slepc_solver.h | 36 +- .../lac/sparse_decomposition.templates.h | 14 +- deal.II/include/deal.II/lac/sparse_matrix.h | 24 +- .../deal.II/lac/sparse_mic.templates.h | 22 +- .../deal.II/matrix_free/dof_info.templates.h | 12 +- .../matrix_free/mapping_info.templates.h | 14 +- .../include/deal.II/matrix_free/matrix_free.h | 4 +- .../deal.II/numerics/vector_tools.templates.h | 2 +- deal.II/source/base/mpi.cc | 134 ++--- deal.II/source/base/utilities.cc | 4 +- deal.II/source/dofs/dof_handler_policy.cc | 4 +- deal.II/source/dofs/dof_tools.cc | 12 +- deal.II/source/fe/fe_q_dg0.cc | 488 +++++++++--------- deal.II/source/grid/grid_generator.cc | 156 +++--- deal.II/source/lac/petsc_matrix_base.cc | 2 +- deal.II/source/lac/petsc_parallel_vector.cc | 44 +- deal.II/source/lac/petsc_solver.cc | 6 +- deal.II/source/lac/petsc_vector_base.cc | 2 +- deal.II/source/lac/slepc_solver.cc | 114 ++-- deal.II/source/lac/trilinos_sparse_matrix.cc | 20 +- deal.II/source/multigrid/mg_tools.cc | 4 +- deal.II/source/numerics/matrix_tools.cc | 26 +- 30 files changed, 647 insertions(+), 647 deletions(-) diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc index bce17c5bb5..743049b17e 100644 --- a/deal.II/examples/step-36/step-36.cc +++ b/deal.II/examples/step-36/step-36.cc @@ -409,7 +409,7 @@ namespace Step36 const unsigned int n_iterations = solve (); std::cout << " Solver converged in " << n_iterations - << " iterations." << std::endl; + << " iterations." << std::endl; output_results (); diff --git a/deal.II/examples/step-49/step-49.cc b/deal.II/examples/step-49/step-49.cc index bdfb6261e9..3e564a9f3a 100644 --- a/deal.II/examples/step-49/step-49.cc +++ b/deal.II/examples/step-49/step-49.cc @@ -27,38 +27,38 @@ using namespace dealii; // @sect3{Generate some output for a given mesh} template -void mesh_info(Triangulation & tria, const char * filename) +void mesh_info(Triangulation &tria, const char *filename) { - // some general info: + // some general info: std::cout << "Mesh info:" << std::endl - << " dimension: " << dim << std::endl - << " no. of cells: " << tria.n_active_cells() << std::endl; + << " dimension: " << dim << std::endl + << " no. of cells: " << tria.n_active_cells() << std::endl; // loop over all the cells and find how often each boundary indicator is used: { std::map boundary_count; Triangulation<2>::active_cell_iterator - cell = tria.begin_active(), - endc = tria.end(); + cell = tria.begin_active(), + endc = tria.end(); for (; cell!=endc; ++cell) { - for (unsigned int face=0; face::faces_per_cell; ++face) - { - if (cell->face(face)->at_boundary()) - boundary_count[cell->face(face)->boundary_indicator()]++; - } + for (unsigned int face=0; face::faces_per_cell; ++face) + { + if (cell->face(face)->at_boundary()) + boundary_count[cell->face(face)->boundary_indicator()]++; + } } - + std::cout << " boundary indicators: "; for (std::map::iterator it=boundary_count.begin(); - it!=boundary_count.end(); - ++it) + it!=boundary_count.end(); + ++it) { - std::cout << it->first << "(" << it->second << ") "; + std::cout << it->first << "(" << it->second << ") "; } std::cout << std::endl; } - + // Now we want to write a graphical representation of the mesh to an output // file. std::ofstream out (filename); @@ -79,9 +79,9 @@ void grid_1 () gridin.attach_triangulation(triangulation); std::ifstream f("untitled.msh"); gridin.read_msh(f); - + mesh_info(triangulation, "grid-1.eps"); - + } @@ -93,16 +93,16 @@ void grid_2 () Triangulation<2> tria1, tria2; - + GridGenerator::hyper_cube_with_cylindrical_hole (tria1, 0.25, 1.0); std::vector< unsigned int > repetitions(2); repetitions[0]=3; repetitions[1]=2; GridGenerator::subdivided_hyper_rectangle (tria2, repetitions, - Point<2>(1.0,-1.0), Point<2>(4.0,1.0)); + Point<2>(1.0,-1.0), Point<2>(4.0,1.0)); + + GridGenerator::merge_triangulations (tria1, tria2, triangulation); - GridGenerator::merge_triangulations (tria1, tria2, triangulation); - mesh_info(triangulation, "grid-2.eps"); } @@ -113,27 +113,27 @@ void grid_3 () GridGenerator::hyper_cube_with_cylindrical_hole (triangulation, 0.25, 1.0); Triangulation<2>::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); + cell = triangulation.begin_active(), + endc = triangulation.end(); for (; cell!=endc; ++cell) { - for (unsigned int i=0;i::vertices_per_cell;++i) - { - Point<2> & v = cell->vertex(i); - if (std::abs(v(1)-1.0)<1e-5) - v(1)+=0.5; - } - + for (unsigned int i=0; i::vertices_per_cell; ++i) + { + Point<2> &v = cell->vertex(i); + if (std::abs(v(1)-1.0)<1e-5) + v(1)+=0.5; + } + } - + const HyperBallBoundary<2> boundary_description(Point<2>(0,0),0.25); triangulation.set_boundary (1, boundary_description); triangulation.refine_global(2); mesh_info(triangulation, "grid-3.eps"); - // remove boundary object from Triangulation again: - triangulation.set_boundary (1); + // remove boundary object from Triangulation again: + triangulation.set_boundary (1); } diff --git a/deal.II/include/deal.II/base/utilities.h b/deal.II/include/deal.II/base/utilities.h index c2bf55e65a..cf321eff61 100644 --- a/deal.II/include/deal.II/base/utilities.h +++ b/deal.II/include/deal.II/base/utilities.h @@ -725,7 +725,7 @@ namespace Utilities // we need operator- for the iterator type anyway, do the // test as follows, rather than via 'last >= first' Assert (last - first >= 0, - ExcMessage ("The given iterators do not satisfy the proper ordering.")); + ExcMessage ("The given iterators do not satisfy the proper ordering.")); unsigned int len = last-first; diff --git a/deal.II/include/deal.II/fe/fe_q_dg0.h b/deal.II/include/deal.II/fe/fe_q_dg0.h index abca28a15f..189ad4ab14 100644 --- a/deal.II/include/deal.II/fe/fe_q_dg0.h +++ b/deal.II/include/deal.II/fe/fe_q_dg0.h @@ -531,7 +531,7 @@ private: * constructor. */ void initialize_embedding (); - + /** * Initialize the restriction * matrices. Called from the diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h index e7b3ee3def..39f53c565b 100644 --- a/deal.II/include/deal.II/grid/grid_generator.h +++ b/deal.II/include/deal.II/grid/grid_generator.h @@ -297,8 +297,8 @@ public: static void parallelogram(Triangulation &tria, - const Point (&corners)[dim], - const bool colorize=false); + const Point (&corners)[dim], + const bool colorize=false); /** @@ -308,8 +308,8 @@ public: static void parallelogram(Triangulation &tria, - const Tensor<2,dim> &corners, - const bool colorize=false) DEAL_II_DEPRECATED; + const Tensor<2,dim> &corners, + const bool colorize=false) DEAL_II_DEPRECATED; /** * A parallelepiped. The first corner point is the origin. The @@ -329,12 +329,12 @@ public: * function. */ template - static - void - parallelepiped (Triangulation &tria, - const Point (&corners) [dim], - const bool colorize = false); - + static + void + parallelepiped (Triangulation &tria, + const Point (&corners) [dim], + const bool colorize = false); + /** * A subdivided parallelepiped. The first corner point is the * origin. The dim adjacent points are vectors describing @@ -348,12 +348,12 @@ public: * function. */ template - static - void - subdivided_parallelepiped (Triangulation &tria, - const unsigned int n_subdivisions, - const Point (&corners) [dim], - const bool colorize = false); + static + void + subdivided_parallelepiped (Triangulation &tria, + const unsigned int n_subdivisions, + const Point (&corners) [dim], + const bool colorize = false); /** * A subdivided parallelepiped, ie. the same as above, but where the @@ -364,12 +364,12 @@ public: * function. */ template - static - void - subdivided_parallelepiped (Triangulation &tria, - const unsigned int ( n_subdivisions) [dim], - const Point (&corners) [dim], - const bool colorize = false); + static + void + subdivided_parallelepiped (Triangulation &tria, + const unsigned int ( n_subdivisions) [dim], + const Point (&corners) [dim], + const bool colorize = false); /** diff --git a/deal.II/include/deal.II/lac/block_indices.h b/deal.II/include/deal.II/lac/block_indices.h index 031b9fd5d7..a05d6492f5 100644 --- a/deal.II/include/deal.II/lac/block_indices.h +++ b/deal.II/include/deal.II/lac/block_indices.h @@ -408,7 +408,7 @@ BlockIndices::global_to_local (const unsigned int i) const --block; return std::pair(block, - i-start_indices[block]); + i-start_indices[block]); } diff --git a/deal.II/include/deal.II/lac/block_vector_base.h b/deal.II/include/deal.II/lac/block_vector_base.h index 2de0135288..3d50e964fd 100644 --- a/deal.II/include/deal.II/lac/block_vector_base.h +++ b/deal.II/include/deal.II/lac/block_vector_base.h @@ -1880,8 +1880,8 @@ BlockVectorBase::in_local_range (const types::global_dof_index global_index) const { const std::pair local_index - = block_indices.global_to_local (global_index); - + = block_indices.global_to_local (global_index); + return components[local_index.first].in_local_range (global_index); } diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h index 9d4207339b..c95e72b0a0 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h @@ -148,7 +148,7 @@ ConstraintMatrix::condense (VectorType &vec) const // parallel vectors, this can only work if we can put a compress() // in between, but we don't want to call compress() twice per entry for (std::vector::const_iterator - constraint_line = lines.begin(); + constraint_line = lines.begin(); constraint_line!=lines.end(); ++constraint_line) { // in case the constraint is @@ -160,7 +160,7 @@ ConstraintMatrix::condense (VectorType &vec) const const typename VectorType::value_type old_value = vec(constraint_line->line); for (unsigned int q=0; q!=constraint_line->entries.size(); ++q) - if (vec.in_local_range(constraint_line->entries[q].first) == true) + if (vec.in_local_range(constraint_line->entries[q].first) == true) vec(constraint_line->entries[q].first) += (static_cast (old_value) * @@ -170,9 +170,9 @@ ConstraintMatrix::condense (VectorType &vec) const vec.compress(); for (std::vector::const_iterator - constraint_line = lines.begin(); + constraint_line = lines.begin(); constraint_line!=lines.end(); ++constraint_line) - if (vec.in_local_range(constraint_line->line) == true) + if (vec.in_local_range(constraint_line->line) == true) vec(constraint_line->line) = 0.; } @@ -261,12 +261,12 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, // copy entries if column will not // be condensed away, distribute // otherwise - for (typename SparseMatrix::const_iterator - p = uncondensed.begin(row); - p != uncondensed.end(row); ++p) + for (typename SparseMatrix::const_iterator + p = uncondensed.begin(row); + p != uncondensed.end(row); ++p) if (new_line[p->column()] != -1) condensed.add (new_line[row], - new_line[p->column()], + new_line[p->column()], p->value()); else { @@ -299,9 +299,9 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, else // line must be distributed { - for (typename SparseMatrix::const_iterator - p = uncondensed.begin(row); - p != uncondensed.end(row); ++p) + for (typename SparseMatrix::const_iterator + p = uncondensed.begin(row); + p != uncondensed.end(row); ++p) // for each column: distribute if (new_line[p->column()] != -1) // column is not constrained @@ -332,7 +332,7 @@ ConstraintMatrix::condense (const SparseMatrix &uncondensed, if (use_vectors == true) for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) condensed_vector (new_line[next_constraint->entries[q].first]) - -= p->value() * + -= p->value() * next_constraint->entries[q].second * c->inhomogeneity; } @@ -1790,7 +1790,7 @@ namespace internals &sparsity_struct[row_start[row]]; typename SparseMatrix::iterator val_ptr = (sparsity.row_length(row) == 0 ? - sparse_matrix->end() : + sparse_matrix->end() : sparse_matrix->begin(row)); const bool optimize_diagonal = sparsity.optimize_diagonal(); unsigned int counter = optimize_diagonal; diff --git a/deal.II/include/deal.II/lac/slepc_solver.h b/deal.II/include/deal.II/lac/slepc_solver.h index 68f72fbac1..8277f35b0f 100644 --- a/deal.II/include/deal.II/lac/slepc_solver.h +++ b/deal.II/include/deal.II/lac/slepc_solver.h @@ -311,19 +311,19 @@ namespace SLEPcWrappers /** * A function that can be used in SLEPc as a callback to check on - * convergence. + * convergence. * * @note This function is redundant. */ - static + static int convergence_test (EPS eps, - PetscScalar kr, - PetscScalar ki, - PetscReal residual_norm, - PetscReal *estimated_error, + PetscScalar kr, + PetscScalar ki, + PetscReal residual_norm, + PetscReal *estimated_error, void *solver_control); - + /** * Objects of this type are explicitly created, but are destroyed * when the surrounding solver object goes out of scope, or when @@ -501,7 +501,7 @@ namespace SLEPcWrappers * solver. Usage: Largest values of spectrum only, all problem * types, complex. * - * @ingroup SLEPcWrappers + * @ingroup SLEPcWrappers * * @author Toby D. Young 2010 */ @@ -543,7 +543,7 @@ namespace SLEPcWrappers * An implementation of the solver interface using the SLEPc * Davidson solver. Usage (incomplete/untested): All problem types. * - * @ingroup SLEPcWrappers + * @ingroup SLEPcWrappers * * @author Toby D. Young 2010 */ @@ -564,8 +564,8 @@ namespace SLEPcWrappers * change that. */ SolverGeneralizedDavidson (SolverControl &cn, - const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, - const AdditionalData &data = AdditionalData()); + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const AdditionalData &data = AdditionalData()); protected: @@ -585,7 +585,7 @@ namespace SLEPcWrappers * An implementation of the solver interface using the SLEPc * Jacobi-Davidson solver. Usage: All problem types. * - * @ingroup SLEPcWrappers + * @ingroup SLEPcWrappers * * @author Toby D. Young 2013 */ @@ -606,8 +606,8 @@ namespace SLEPcWrappers * change that. */ SolverJacobiDavidson (SolverControl &cn, - const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, - const AdditionalData &data = AdditionalData()); + const MPI_Comm &mpi_communicator = PETSC_COMM_SELF, + const AdditionalData &data = AdditionalData()); protected: @@ -638,12 +638,12 @@ namespace SLEPcWrappers { unsigned int n_converged = 0; - // Set the matrices of the problem + // Set the matrices of the problem set_matrices (A); - // and solve + // and solve solve (n_eigenvectors, &n_converged); - + if (n_converged > n_eigenvectors) n_converged = n_eigenvectors; AssertThrow (n_converged == n_eigenvectors, @@ -667,7 +667,7 @@ namespace SLEPcWrappers { unsigned int n_converged = 0; - // Set the matrices of the problem + // Set the matrices of the problem set_matrices (A, B); // and solve diff --git a/deal.II/include/deal.II/lac/sparse_decomposition.templates.h b/deal.II/include/deal.II/lac/sparse_decomposition.templates.h index 6527a49ce9..0f4ca7017f 100644 --- a/deal.II/include/deal.II/lac/sparse_decomposition.templates.h +++ b/deal.II/include/deal.II/lac/sparse_decomposition.templates.h @@ -117,9 +117,9 @@ void SparseLUDecomposition::initialize ( // and recreate own_sparsity = new SparsityPattern(matrix_sparsity, - matrix_sparsity.max_entries_per_row() - +2*data.extra_off_diagonals, - data.extra_off_diagonals); + matrix_sparsity.max_entries_per_row() + +2*data.extra_off_diagonals, + data.extra_off_diagonals); own_sparsity->compress(); sparsity_pattern_to_use = own_sparsity; } @@ -256,16 +256,16 @@ SparseLUDecomposition::strengthen_diagonal_impl () // non-diagonal element in this row Assert (this->cols->optimize_diagonal(), ExcNotImplemented()); typename SparseMatrix::iterator - diagonal_element = this->begin(row); + diagonal_element = this->begin(row); number rowsum = 0; for (typename SparseMatrix::iterator - p = diagonal_element + 1; - p != this->end(row); ++p) + p = diagonal_element + 1; + p != this->end(row); ++p) rowsum += std::fabs(p->value()); diagonal_element->value() += this->get_strengthen_diagonal (rowsum, row) * - rowsum; + rowsum; } } diff --git a/deal.II/include/deal.II/lac/sparse_matrix.h b/deal.II/include/deal.II/lac/sparse_matrix.h index 9c8527715a..0a2d7e08ed 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.h @@ -308,7 +308,7 @@ namespace SparseMatrixIterators * of the current kind. */ typedef - const Accessor & value_type; + const Accessor &value_type; /** * Constructor. Create an iterator into the matrix @p matrix for the given @@ -1556,7 +1556,7 @@ private: * doing here. */ number nth_entry_in_row (const unsigned int row, - const unsigned int index) const; + const unsigned int index) const; // make all other sparse matrices friends template friend class SparseMatrix; @@ -1893,7 +1893,7 @@ template inline number SparseMatrix::nth_entry_in_row (const unsigned int row, - const unsigned int index) const + const unsigned int index) const { Assert(rowrows, ExcIndexRange(row,0,cols->rows)); Assert(indexrow_length(row), @@ -2029,7 +2029,7 @@ namespace SparseMatrixIterators Accessor::Reference::operator number() const { return accessor->matrix->nth_entry_in_row(accessor->a_row, - accessor->a_index); + accessor->a_index); } @@ -2271,17 +2271,17 @@ namespace SparseMatrixIterators const unsigned int this_position = (*this != (*this)->get_matrix().end() - ? - sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index() - : - sparsity.get_rowstart_indices()[sparsity.n_rows()]); + ? + sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index() + : + sparsity.get_rowstart_indices()[sparsity.n_rows()]); const unsigned int other_position = (other != (*this)->get_matrix().end() - ? - sparsity.get_rowstart_indices()[other->row()] + other->index() - : - sparsity.get_rowstart_indices()[sparsity.n_rows()]); + ? + sparsity.get_rowstart_indices()[other->row()] + other->index() + : + sparsity.get_rowstart_indices()[sparsity.n_rows()]); return (this_position - other_position); } diff --git a/deal.II/include/deal.II/lac/sparse_mic.templates.h b/deal.II/include/deal.II/lac/sparse_mic.templates.h index d4882c07fc..fc322ef102 100644 --- a/deal.II/include/deal.II/lac/sparse_mic.templates.h +++ b/deal.II/include/deal.II/lac/sparse_mic.templates.h @@ -143,9 +143,9 @@ void SparseMIC::decompose (const SparseMatrix &matrix, // work on the lower left part of the matrix. we know // it's symmetric, so we can work with this alone for (typename SparseMatrix::const_iterator - p = matrix.begin(row)+1; - (p != matrix.end(row)) && (p->column() < p->row()); - ++p) + p = matrix.begin(row)+1; + (p != matrix.end(row)) && (p->column() < p->row()); + ++p) temp1 += p->value() / diag[p->column()] * inner_sums[p->column()]; Assert(temp-temp1 > 0, ExcStrengthenDiagonalTooSmall()); @@ -166,8 +166,8 @@ SparseMIC::get_rowsum (const unsigned int row) const const unsigned int *const first_after_diagonal = this->prebuilt_lower_bound[row]; number rowsum = 0; for (typename SparseMatrix::const_iterator - p = this->begin(row)+1; - p != this->end(row); ++p) + p = this->begin(row)+1; + p != this->end(row); ++p) if (p->column() > p->row()) rowsum += p->value(); @@ -204,9 +204,9 @@ SparseMIC::vmult (Vector &dst, // get start of this row. skip // the diagonal element for (typename SparseMatrix::const_iterator - p = this->begin(row)+1; - (p != this->end(row)) && (p->column() < row); - ++p) + p = this->begin(row)+1; + (p != this->end(row)) && (p->column() < row); + ++p) dst(row) -= p->value() * dst(p->column()); dst(row) *= inv_diag[row]; @@ -221,9 +221,9 @@ SparseMIC::vmult (Vector &dst, { // get end of this row for (typename SparseMatrix::const_iterator - p = this->begin(row)+1; - p != this->end(row); - ++p) + p = this->begin(row)+1; + p != this->end(row); + ++p) if (p->column() > static_cast(row)) dst(row) -= p->value() * dst(p->column()); diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index bf4003fedd..1673b078eb 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -1767,8 +1767,8 @@ not_connect: ++cell) { const unsigned int - *it = begin_indices (renumbering[cell]), - *end_cell = end_indices (renumbering[cell]); + *it = begin_indices (renumbering[cell]), + *end_cell = end_indices (renumbering[cell]); for ( ; it != end_cell; ++it) if (row_lengths[*it]>0) connectivity_dof.add(*it, block); @@ -1780,8 +1780,8 @@ not_connect: else { const unsigned int - *it = begin_indices (block), - *end_cell = end_indices (block); + *it = begin_indices (block), + *end_cell = end_indices (block); for ( ; it != end_cell; ++it) if (row_lengths[*it]>0) connectivity_dof.add(*it, block); @@ -1811,8 +1811,8 @@ not_connect: { // apply renumbering when we do blocking const unsigned int - *it = begin_indices (renumbering[cell]), - *end_cell = end_indices (renumbering[cell]); + *it = begin_indices (renumbering[cell]), + *end_cell = end_indices (renumbering[cell]); for ( ; it != end_cell; ++it) if (row_lengths[*it] > 0) { diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h index 59ed7db300..4d964d1289 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h @@ -270,9 +270,9 @@ end_set: FPArrayComparator comparator(jacobian_size); typedef Tensor<1,VectorizedArray::n_array_elements,Number> VEC_ARRAY; std::map, unsigned int, - FPArrayComparator > cartesians(comparator); + FPArrayComparator > cartesians(comparator); std::map, unsigned int, - FPArrayComparator > affines(comparator); + FPArrayComparator > affines(comparator); // loop over all cells for (unsigned int cell=0; cell, - unsigned int, FPArrayComparator >::iterator - it = cartesians.begin(); it != cartesians.end(); ++it) + unsigned int, FPArrayComparator >::iterator + it = cartesians.begin(); it != cartesians.end(); ++it) { VectorizedArray det = make_vectorized_array(1.); for (unsigned int d=0; d, - unsigned int, FPArrayComparator >::iterator - it = affines.begin(); it != affines.end(); ++it) + unsigned int, FPArrayComparator >::iterator + it = affines.begin(); it != affines.end(); ++it) { Tensor<2,dim,VectorizedArray > jac; for (unsigned int d=0; d::constraint_pool_begin (const unsigned int row) const { AssertIndexRange (row, constraint_pool_row_index.size()-1); return constraint_pool_data.empty() ? 0 : - &constraint_pool_data[0] + constraint_pool_row_index[row]; + &constraint_pool_data[0] + constraint_pool_row_index[row]; } @@ -1080,7 +1080,7 @@ MatrixFree::constraint_pool_end (const unsigned int row) const { AssertIndexRange (row, constraint_pool_row_index.size()-1); return constraint_pool_data.empty() ? 0 : - &constraint_pool_data[0] + constraint_pool_row_index[row+1]; + &constraint_pool_data[0] + constraint_pool_row_index[row+1]; } diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index b90e36daa3..8a81c86a50 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -2147,7 +2147,7 @@ namespace VectorTools // set up sparsity structure CompressedSparsityPattern c_sparsity(dof.n_boundary_dofs (boundary_functions), - dof.n_boundary_dofs (boundary_functions)); + dof.n_boundary_dofs (boundary_functions)); DoFTools::make_boundary_sparsity_pattern (dof, boundary_functions, dof_to_boundary_mapping, diff --git a/deal.II/source/base/mpi.cc b/deal.II/source/base/mpi.cc index 5737110d99..ed1d846780 100644 --- a/deal.II/source/base/mpi.cc +++ b/deal.II/source/base/mpi.cc @@ -90,7 +90,7 @@ namespace Utilities { int n_jobs=1; (void) MPI_Comm_size (mpi_communicator, &n_jobs); - + return n_jobs; } @@ -99,26 +99,26 @@ namespace Utilities { int rank=0; (void) MPI_Comm_rank (mpi_communicator, &rank); - + return rank; } - - + + MPI_Comm duplicate_communicator (const MPI_Comm &mpi_communicator) { MPI_Comm new_communicator; MPI_Comm_dup (mpi_communicator, &new_communicator); return new_communicator; } - - + + std::vector compute_point_to_point_communication_pattern (const MPI_Comm &mpi_comm, const std::vector &destinations) { unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm); unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_comm); - + for (unsigned int i=0; i(in_lhs_); MinMaxAvg *inout_rhs = static_cast(inout_rhs_); - + Assert(*len==1, ExcInternalError()); - + inout_rhs->sum += in_lhs->sum; if (inout_rhs->min>in_lhs->min) { @@ -204,7 +204,7 @@ namespace Utilities if (inout_rhs->min_index > in_lhs->min_index) inout_rhs->min_index = in_lhs->min_index; } - + if (inout_rhs->max < in_lhs->max) { inout_rhs->max = in_lhs->max; @@ -218,24 +218,24 @@ namespace Utilities } } } - - - + + + MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &mpi_communicator) { MinMaxAvg result; - + const unsigned int my_id = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); const unsigned int numproc = dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); - + MPI_Op op; int ierr = MPI_Op_create((MPI_User_function *)&max_reduce, true, &op); AssertThrow(ierr == MPI_SUCCESS, ExcInternalError()); - + MinMaxAvg in; in.sum = in.min = in.max = my_value; in.min_index = in.max_index = my_id; @@ -264,34 +264,34 @@ namespace Utilities } #else - + unsigned int n_mpi_processes (const MPI_Comm &) { return 1; } - - - + + + unsigned int this_mpi_process (const MPI_Comm &) { return 0; } - - + + MPI_Comm duplicate_communicator (const MPI_Comm &mpi_communicator) { return mpi_communicator; } - - - - + + + + MinMaxAvg min_max_avg(const double my_value, const MPI_Comm &) { MinMaxAvg result; - + result.sum = my_value; result.avg = my_value; result.min = my_value; @@ -333,7 +333,7 @@ namespace Utilities MPI_Initialized(&MPI_has_been_started); AssertThrow (MPI_has_been_started == 0, ExcMessage ("MPI error. You can only start MPI once!")); - + int mpi_err; mpi_err = MPI_Init (&argc, &argv); AssertThrow (mpi_err == 0, @@ -346,15 +346,15 @@ namespace Utilities (void)argv; (void)owns_mpi; #endif - + constructor_has_already_run = true; } - - + + MPI_InitFinalize::~MPI_InitFinalize() { #ifdef DEAL_II_COMPILER_SUPPORTS_MPI - + // make memory pool release all MPI-based vectors that are no // longer used at this point. this is relevant because the // static object destructors run for these vectors at the end of @@ -362,11 +362,11 @@ namespace Utilities // to errors # if defined(DEAL_II_USE_TRILINOS) GrowingVectorMemory - ::release_unused_memory (); + ::release_unused_memory (); GrowingVectorMemory - ::release_unused_memory (); + ::release_unused_memory (); # endif - + // Same for PETSc. only do this if PETSc hasn't been // terminated yet since PETSc deletes all vectors that // have been allocated but not freed at the time of @@ -375,31 +375,31 @@ namespace Utilities // yield errors of double deallocations #ifdef DEAL_II_USE_PETSC if ((PetscInitializeCalled == PETSC_TRUE) - && - (PetscFinalizeCalled == PETSC_FALSE)) - { - GrowingVectorMemory - ::release_unused_memory (); - GrowingVectorMemory - ::release_unused_memory (); - GrowingVectorMemory - ::release_unused_memory (); - GrowingVectorMemory - ::release_unused_memory (); - + && + (PetscFinalizeCalled == PETSC_FALSE)) + { + GrowingVectorMemory + ::release_unused_memory (); + GrowingVectorMemory + ::release_unused_memory (); + GrowingVectorMemory + ::release_unused_memory (); + GrowingVectorMemory + ::release_unused_memory (); + # ifdef DEAL_II_USE_SLEPC - // and now end SLEPc (with PETSc) - SlepcFinalize(); + // and now end SLEPc (with PETSc) + SlepcFinalize(); # else - // or just end PETSc. - PetscFinalize(); + // or just end PETSc. + PetscFinalize(); # endif - } + } #else - + int mpi_err = 0; - + int MPI_has_been_started = 0; MPI_Initialized(&MPI_has_been_started); if (Utilities::System::job_supports_mpi() == true && owns_mpi == true && @@ -415,17 +415,17 @@ namespace Utilities else mpi_err = MPI_Finalize(); } - - + + AssertThrow (mpi_err == 0, ExcMessage ("An error occurred while calling MPI_Finalize()")); #endif #endif } - - + + } // end of namespace MPI - + } // end of namespace Utilities DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/base/utilities.cc b/deal.II/source/base/utilities.cc index b96b284b27..9122243fca 100644 --- a/deal.II/source/base/utilities.cc +++ b/deal.II/source/base/utilities.cc @@ -82,10 +82,10 @@ namespace Utilities (digits==3 && i>=1000) || (digits==4 && i>=10000)|| (digits==5 && i>=100000)|| - (digits==6 && i>=1000000)|| + (digits==6 && i>=1000000)|| (digits==7 && i>=10000000)|| (digits==8 && i>=100000000)|| - (digits==9 && i>=1000000000)|| + (digits==9 && i>=1000000000)|| (i>=1000000000)), ExcInvalidNumber2StringConversersion(i, digits)); diff --git a/deal.II/source/dofs/dof_handler_policy.cc b/deal.II/source/dofs/dof_handler_policy.cc index 4ef5c615d0..7de2d8a957 100644 --- a/deal.II/source/dofs/dof_handler_policy.cc +++ b/deal.II/source/dofs/dof_handler_policy.cc @@ -2135,8 +2135,8 @@ namespace internal Assert (tr != 0, ExcInternalError()); AssertThrow( - (tr->settings & parallel::distributed::Triangulation< dim, spacedim >::construct_multigrid_hierarchy), - ExcMessage("Multigrid DoFs can only be distributed on a parallel Triangulation if the flag construct_multigrid_hierarchy is set in the constructor.")); + (tr->settings & parallel::distributed::Triangulation< dim, spacedim >::construct_multigrid_hierarchy), + ExcMessage("Multigrid DoFs can only be distributed on a parallel Triangulation if the flag construct_multigrid_hierarchy is set in the constructor.")); const unsigned int diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index 02e722a86f..af93bf0d45 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -3013,10 +3013,10 @@ namespace DoFTools // parallel::distributed::Triangulation in serial { typedef parallel::distributed::Triangulation PTRIA; - const PTRIA* ptria_p = dynamic_cast (&dof_handler.get_tria()); + const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); + ExcMessage ("This function can not be used with distributed triangulations." + "See the documentation for more information.")); } #endif @@ -3107,10 +3107,10 @@ namespace DoFTools // parallel::distributed::Triangulation in serial { typedef typename parallel::distributed::Triangulation PTRIA; - const PTRIA* ptria_p = dynamic_cast (&dof_handler.get_tria()); + const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); + ExcMessage ("This function can not be used with distributed triangulations." + "See the documentation for more information.")); } #endif diff --git a/deal.II/source/fe/fe_q_dg0.cc b/deal.II/source/fe/fe_q_dg0.cc index f6bae5220c..50e01f5834 100644 --- a/deal.II/source/fe/fe_q_dg0.cc +++ b/deal.II/source/fe/fe_q_dg0.cc @@ -71,7 +71,7 @@ namespace FE_Q_DG0_Helper for (unsigned int d=0; d > points1d (dofs1d); for (unsigned int i=0; i(unit_support_points[j](0)); - for (unsigned int d=1; d(unit_support_points[j](0)); + for (unsigned int d=1; d::FE_Q_DG0 (const unsigned int degree) FiniteElementData::L2), get_riaf_vector(degree), std::vector(1ll , std::vector(1,true))), - face_index_map(FE_Q_DG0_Helper::invert_numbering( - face_lexicographic_to_hierarchic_numbering (degree))) + face_index_map(FE_Q_DG0_Helper::invert_numbering( + face_lexicographic_to_hierarchic_numbering (degree))) { Assert (degree > 0, ExcMessage ("This element can only be used for polynomial degrees " @@ -529,7 +529,7 @@ FE_Q_DG0::FE_Q_DG0 (const unsigned int degree) this->reinit_restriction_and_prolongation_matrices(); initialize_embedding(); initialize_restriction(); - + initialize_quad_dof_index_permutation(); } @@ -546,8 +546,8 @@ FE_Q_DG0::FE_Q_DG0 (const Quadrature<1> &points) get_riaf_vector(points.size()-1), std::vector(1, std::vector(1,true))), - face_index_map(FE_Q_DG0_Helper::invert_numbering( - face_lexicographic_to_hierarchic_numbering (points.size()-1))) + face_index_map(FE_Q_DG0_Helper::invert_numbering( + face_lexicographic_to_hierarchic_numbering (points.size()-1))) { const int degree = points.size()-1; @@ -579,7 +579,7 @@ FE_Q_DG0::FE_Q_DG0 (const Quadrature<1> &points) initialize_embedding(); initialize_restriction(); - initialize_quad_dof_index_permutation(); + initialize_quad_dof_index_permutation(); } @@ -672,12 +672,12 @@ FE_Q_DG0::interpolate(std::vector &local_dofs, { Assert (values.size() == this->unit_support_points.size(), ExcDimensionMismatch(values.size(), - this->unit_support_points.size())); + this->unit_support_points.size())); Assert (local_dofs.size() == this->dofs_per_cell, ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (this->n_components() == 1, ExcDimensionMismatch(this->n_components(), 1)); - + std::copy(values.begin(), values.end(), local_dofs.begin()); //We don't need the discontinuous function for local interpolation local_dofs[local_dofs.size()-1]=0.; @@ -691,25 +691,25 @@ FE_Q_DG0::interpolate(std::vector &local_dofs, { Assert (values.size() == this->unit_support_points.size(), ExcDimensionMismatch(values.size(), - this->unit_support_points.size())); + this->unit_support_points.size())); Assert (local_dofs.size() == this->dofs_per_cell, ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (values[0].size() >= offset+this->n_components(), ExcDimensionMismatch(values[0].size(),offset+this->n_components())); - + for (unsigned int i=0; idofs_per_cell-1; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[i](offset+index.first); - } + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[i](offset+index.first); + } //We don't need the discontinuous function for local interpolation local_dofs[local_dofs.size()-1]=0.; } - - - - + + + + template void FE_Q_DG0::interpolate( @@ -718,18 +718,18 @@ FE_Q_DG0::interpolate( { Assert (values[0].size() == this->unit_support_points.size(), ExcDimensionMismatch(values.size(), - this->unit_support_points.size())); + this->unit_support_points.size())); Assert (local_dofs.size() == this->dofs_per_cell, ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (values.size() == this->n_components(), ExcDimensionMismatch(values.size(), this->n_components())); - + for (unsigned int i=0; idofs_per_cell-1; ++i) - { - const std::pair index - = this->system_to_component_index(i); - local_dofs[i] = values[index.first][i]; - } + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[index.first][i]; + } //We don't need the discontinuous function for local interpolation local_dofs[local_dofs.size()-1]=0.; } @@ -782,16 +782,16 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, const Point p = this->unit_support_points[j]; for (unsigned int i=0; idofs_per_cell-1; ++i) cell_interpolation(j,i) = this->poly_space.compute_value (i, p); - + for (unsigned int i=0; idofs_per_cell-1, - this->dofs_per_cell-1)=1.; - source_interpolation(this->dofs_per_cell-1, - source_fe.dofs_per_cell-1)=1.; + //the discontinuous node is not transformed + cell_interpolation(this->dofs_per_cell-1, + this->dofs_per_cell-1)=1.; + source_interpolation(this->dofs_per_cell-1, + source_fe.dofs_per_cell-1)=1.; // then compute the // interpolation matrix @@ -1179,7 +1179,7 @@ hp_line_dof_identities (const FiniteElement &fe_other) const // FE_Q_DG0s or if the other one is an // FE_Nothing if (const FE_Q_DG0 *fe_q_dg0_other = dynamic_cast*>(&fe_other)) + FE_Q_DG0 *>(&fe_other)) { // dofs are located along lines, so two // dofs are identical if they are @@ -1211,7 +1211,7 @@ FE_Q_DG0*>(&fe_other)) for (unsigned int j=0; junit_support_points[index_map_inverse[i+1]][0]- -fe_q_dg0_other->unit_support_points[index_map_inverse_other[j+1]][0]) + fe_q_dg0_other->unit_support_points[index_map_inverse_other[j+1]][0]) < 1e-14) identities.push_back (std::make_pair(i,j)); @@ -1244,7 +1244,7 @@ hp_quad_dof_identities (const FiniteElement &fe_other) cons // FE_Q_DG0s or if the other one is an // FE_Nothing if (const FE_Q_DG0 *fe_q_dg0_other = dynamic_cast*>(&fe_other)) + FE_Q_DG0 *>(&fe_other)) { // this works exactly like the line // case above, except that now we have @@ -1271,12 +1271,12 @@ FE_Q_DG0*>(&fe_other)) for (unsigned int j2=0; j2unit_support_points[index_map_inverse[i1+1]][0]- -fe_q_dg0_other->unit_support_points[index_map_inverse_other[j1+1]][0]) + fe_q_dg0_other->unit_support_points[index_map_inverse_other[j1+1]][0]) < 1e-14) && (std::fabs(this->unit_support_points[index_map_inverse[i2+1]][0]- -fe_q_dg0_other->unit_support_points[index_map_inverse_other[j2+1]][0]) + fe_q_dg0_other->unit_support_points[index_map_inverse_other[j2+1]][0]) < 1e-14)) identities.push_back (std::make_pair(i1*(p-1)+i2, j1*(q-1)+j2)); @@ -1594,8 +1594,8 @@ std::vector FE_Q_DG0::get_riaf_vector(const unsigned int deg) { std::vector riaf - (FiniteElementData (get_dpo_vector(deg),1,deg).dofs_per_cell, - false); + (FiniteElementData (get_dpo_vector(deg),1,deg).dofs_per_cell, + false); riaf[riaf.size()-1]=true; return riaf; } @@ -1681,154 +1681,154 @@ FE_Q_DG0::initialize_embedding () // cell matrix. The value eps is used a threshold to decide when certain // evaluations of the Lagrange polynomials are zero or one. const double eps = 1e-15*this->degree*dim; - - #ifdef DEBUG + +#ifdef DEBUG // in DEBUG mode, check that the evaluation of support points (except for // the discontinuous node) in the current numbering gives the identity // operation for (unsigned int i=0; idofs_per_cell-1; ++i) - { - Assert (std::fabs (1.-this->poly_space.compute_value - (i, this->unit_support_points[i])) < eps, - ExcInternalError()); - for (unsigned int j=0; jdofs_per_cell-1; ++j) - if (j!=i) - Assert (std::fabs (this->poly_space.compute_value - (i, this->unit_support_points[j])) < eps, - ExcInternalError()); - } - #endif - + { + Assert (std::fabs (1.-this->poly_space.compute_value + (i, this->unit_support_points[i])) < eps, + ExcInternalError()); + for (unsigned int j=0; jdofs_per_cell-1; ++j) + if (j!=i) + Assert (std::fabs (this->poly_space.compute_value + (i, this->unit_support_points[j])) < eps, + ExcInternalError()); + } +#endif + // to efficiently evaluate the polynomial at the subcell, make use of the // tensor product structure of this element and only evaluate 1D // information from the polynomial. This makes the cost of this function // almost negligible also for high order elements const unsigned int dofs1d = this->degree+1; std::vector > subcell_evaluations - (dim, Table<2,double>(dofs1d, dofs1d)); + (dim, Table<2,double>(dofs1d, dofs1d)); const std::vector &index_map_inverse = this->poly_space.get_numbering_inverse(); - + // recreate 1D polynomials, the last entry is not used std::vector > poly_space1d = FE_Q_DG0_Helper::generate_poly_space1d (this->unit_support_points, - index_map_inverse, dofs1d); - + index_map_inverse, dofs1d); + // helper value: step size how to walk through diagonal and how many // points we have left apart from the first dimension unsigned int step_size_diag = 0; { unsigned int factor = 1; for (unsigned int d=0; d::isotropic_refinement; ++ref) for (unsigned int child=0; child::n_children(RefinementCase(ref+1)); ++child) - { - // go through the points in diagonal to capture variation in all - // directions simultaneously - for (unsigned int j=0; j p_subcell = this->unit_support_points[diag_comp]; - const Point p_cell = + { + // go through the points in diagonal to capture variation in all + // directions simultaneously + for (unsigned int j=0; j p_subcell = this->unit_support_points[diag_comp]; + const Point p_cell = GeometryInfo::child_to_cell_coordinates - (p_subcell, child, RefinementCase(ref+1)); - for (unsigned int i=0; i (j_indices); - for (unsigned int j=0; jdofs_per_cell-1; j+=dofs1d) - { - unsigned int i_indices[dim]; - FE_Q_DG0_Helper::zero_indices (i_indices); - for (unsigned int i=0; idofs_per_cell-1; i+=dofs1d) - { - double val_extra_dim = 1.; - for (unsigned int d=1; dprolongation[ref][child](j,i) = - // this->poly_space.compute_value (i, p_cell); - for (unsigned int jj=0; jjprolongation[ref][child](j_ind,index_map_inverse[i+ii]) - = val_extra_dim * subcell_evaluations[0](jj,ii); - } - - // update indices that denote the tensor product position. a bit - // fuzzy and therefore not done for innermost x_0 direction - FE_Q_DG0_Helper::increment_indices (i_indices, dofs1d); - } - Assert (i_indices[dim-1] == 1, ExcInternalError()); - FE_Q_DG0_Helper::increment_indices (j_indices, dofs1d); - } - - // the discontinuous node is simply mapped on the discontinuous node - // on the child cell - for (unsigned int j=0; jdofs_per_cell-1; j++) - { - this->prolongation[ref][child](j,this->dofs_per_cell-1) = 0.; - this->prolongation[ref][child](this->dofs_per_cell-1,j) = 0.; - } - this->prolongation[ref][child](this->dofs_per_cell-1, - this->dofs_per_cell-1) = 1.; - - // make sure that the row sum is one. this must be so since for this - // element, the continuous shape functions add up to one and the - // discontinuous node is mapped to the discontinuous node on the - // child cell - #ifdef DEBUG - for (unsigned int row=0; rowdofs_per_cell; ++row) - { - double sum = 0; - for (unsigned int col=0; coldofs_per_cell; ++col) - sum += this->prolongation[ref][child](row,col); - Assert (std::fabs(sum-1.) < eps, ExcInternalError()); - } - #endif - } + (p_subcell, child, RefinementCase(ref+1)); + for (unsigned int i=0; i (j_indices); + for (unsigned int j=0; jdofs_per_cell-1; j+=dofs1d) + { + unsigned int i_indices[dim]; + FE_Q_DG0_Helper::zero_indices (i_indices); + for (unsigned int i=0; idofs_per_cell-1; i+=dofs1d) + { + double val_extra_dim = 1.; + for (unsigned int d=1; dprolongation[ref][child](j,i) = + // this->poly_space.compute_value (i, p_cell); + for (unsigned int jj=0; jjprolongation[ref][child](j_ind,index_map_inverse[i+ii]) + = val_extra_dim * subcell_evaluations[0](jj,ii); + } + + // update indices that denote the tensor product position. a bit + // fuzzy and therefore not done for innermost x_0 direction + FE_Q_DG0_Helper::increment_indices (i_indices, dofs1d); + } + Assert (i_indices[dim-1] == 1, ExcInternalError()); + FE_Q_DG0_Helper::increment_indices (j_indices, dofs1d); + } + + // the discontinuous node is simply mapped on the discontinuous node + // on the child cell + for (unsigned int j=0; jdofs_per_cell-1; j++) + { + this->prolongation[ref][child](j,this->dofs_per_cell-1) = 0.; + this->prolongation[ref][child](this->dofs_per_cell-1,j) = 0.; + } + this->prolongation[ref][child](this->dofs_per_cell-1, + this->dofs_per_cell-1) = 1.; + + // make sure that the row sum is one. this must be so since for this + // element, the continuous shape functions add up to one and the + // discontinuous node is mapped to the discontinuous node on the + // child cell +#ifdef DEBUG + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[ref][child](row,col); + Assert (std::fabs(sum-1.) < eps, ExcInternalError()); + } +#endif + } } @@ -1862,11 +1862,11 @@ FE_Q_DG0::initialize_restriction () // about this. // for the discontinuous node we just take the mean of all the child cells' // contributions. - + const double eps = 1e-15*this->degree*dim; const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); - + this->poly_space.get_numbering_inverse(); + // recreate 1D polynomials for faster evaluation of polynomial, the last // entry is not used const unsigned int dofs1d = this->degree+1; @@ -1874,90 +1874,90 @@ FE_Q_DG0::initialize_restriction () FE_Q_DG0_Helper::generate_poly_space1d (this->unit_support_points, index_map_inverse, dofs1d); std::vector > evaluations1d (dofs1d); - + for (unsigned int i=0; idofs_per_cell-1; ++i) - { - unsigned int mother_dof = index_map_inverse[i]; - const Point p_cell = this->unit_support_points[mother_dof]; - - // then find the children on which the interpolation point is located - for (unsigned int ref=RefinementCase::cut_x; - ref<=RefinementCase::isotropic_refinement; ++ref) - for (unsigned int child=0; - child::n_children(RefinementCase(ref)); - ++child) + { + unsigned int mother_dof = index_map_inverse[i]; + const Point p_cell = this->unit_support_points[mother_dof]; + + // then find the children on which the interpolation point is located + for (unsigned int ref=RefinementCase::cut_x; + ref<=RefinementCase::isotropic_refinement; ++ref) + for (unsigned int child=0; + child::n_children(RefinementCase(ref)); + ++child) + { + // check whether this interpolation point is inside this child + // cell + const Point p_subcell = + GeometryInfo::cell_to_child_coordinates + (p_cell, child, RefinementCase(ref)); + if (GeometryInfo::is_inside_unit_cell (p_subcell)) { - // check whether this interpolation point is inside this child - // cell - const Point p_subcell = - GeometryInfo::cell_to_child_coordinates - (p_cell, child, RefinementCase(ref)); - if (GeometryInfo::is_inside_unit_cell (p_subcell)) - { - // same logic as in initialize_embedding to evaluate the - // polynomial faster than from the tensor product: since we - // evaluate all polynomials, it is much faster to just - // compute the 1D values for all polynomials before and then - // get the dim-data. - for (unsigned int j=0; j (j_indices); - double sum_check = 0; - for (unsigned int j = 0; jdofs_per_cell-1; j += dofs1d) + // same logic as in initialize_embedding to evaluate the + // polynomial faster than from the tensor product: since we + // evaluate all polynomials, it is much faster to just + // compute the 1D values for all polynomials before and then + // get the dim-data. + for (unsigned int j=0; j (j_indices); + double sum_check = 0; + for (unsigned int j = 0; jdofs_per_cell-1; j += dofs1d) { double val_extra_dim = 1.; for (unsigned int d=1; drestriction[ref-1][child] - (mother_dof,child_dof)=1.; - else if (std::fabs(val) > eps) - this->restriction[ref-1][child] - (mother_dof,child_dof)=val; - sum_check += val; - } + { + + // find the child shape function(s) corresponding + // to this point. Usually this is just one function; + // however, when we use FE_Q_DG0 on arbitrary nodes a + // parent support point might not be a child support + // point, and then we will get more than one nonzero + // value per row. Still, the values should sum up to 1 + const double val = val_extra_dim * evaluations1d[jj][0]; + const unsigned int child_dof = index_map_inverse[j+jj]; + if (std::fabs (val-1.) < eps) + this->restriction[ref-1][child] + (mother_dof,child_dof)=1.; + else if (std::fabs(val) > eps) + this->restriction[ref-1][child] + (mother_dof,child_dof)=val; + sum_check += val; + } FE_Q_DG0_Helper::increment_indices (j_indices, dofs1d); } - Assert (std::fabs(sum_check-1) < eps, - ExcInternalError()); - } + Assert (std::fabs(sum_check-1) < eps, + ExcInternalError()); } - } + } + } // for the discontinuous node we take the mean of all the child cells' // contributions. so there is also a diagonal entry equal to the inverse of // n_children. for (unsigned int ref=RefinementCase::cut_x; ref<=RefinementCase::isotropic_refinement; ++ref) - for (unsigned int child=0; - child::n_children(RefinementCase(ref)); - ++child) - { - for (unsigned int i=0; idofs_per_cell-1; ++i) - { - this->restriction[ref-1][child](i,this->dofs_per_cell-1)=0.; - this->restriction[ref-1][child](this->dofs_per_cell-1,i)=0.; - } - this->restriction[ref-1][child] - (this->dofs_per_cell-1,this->dofs_per_cell-1) = - 1./(double)GeometryInfo::n_children - (RefinementCase(ref)); - } - + for (unsigned int child=0; + child::n_children(RefinementCase(ref)); + ++child) + { + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + this->restriction[ref-1][child](i,this->dofs_per_cell-1)=0.; + this->restriction[ref-1][child](this->dofs_per_cell-1,i)=0.; + } + this->restriction[ref-1][child] + (this->dofs_per_cell-1,this->dofs_per_cell-1) = + 1./(double)GeometryInfo::n_children + (RefinementCase(ref)); + } + } - + //--------------------------------------------------------------------------- diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc index 02949471a2..1d91325670 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -371,8 +371,8 @@ GridGenerator::torus (Triangulation<2,3> &tria, template<> void GridGenerator::parallelogram (Triangulation<2> &tria, - const Point<2> (&corners)[2], - const bool colorize) + const Point<2> (&corners)[2], + const bool colorize) { std::vector > vertices (GeometryInfo<2>::vertices_per_cell); @@ -396,8 +396,8 @@ GridGenerator::parallelogram (Triangulation<2> &tria, template<> void GridGenerator::parallelogram (Triangulation<2> &tria, - const Tensor<2,2> &corners, - const bool colorize) + const Tensor<2,2> &corners, + const bool colorize) { // simply pass everything to the other function of same name const Point<2> x[2] = { corners[0], corners[1] }; @@ -417,8 +417,8 @@ GridGenerator::parallelogram (Triangulation<2> &tria, template void GridGenerator::parallelepiped (Triangulation &tria, - const Point (&corners) [dim], - const bool colorize) + const Point (&corners) [dim], + const bool colorize) { // Check that none of the user defined vertices overlap for (unsigned int i=0; i &tria, template void GridGenerator::subdivided_parallelepiped (Triangulation &tria, - const unsigned int n_subdivisions, - const Point (&corners) [dim], - const bool colorize) + const unsigned int n_subdivisions, + const Point (&corners) [dim], + const bool colorize) { // Equalise number of subdivisions in each dim-direction, heir // validity will be checked later @@ -493,17 +493,17 @@ GridGenerator::subdivided_parallelepiped (Triangulation &tria, n_subdivisions_[i] = n_subdivisions; // and call the function below - GridGenerator::subdivided_parallelepiped (tria, n_subdivisions_, - corners, - colorize); + GridGenerator::subdivided_parallelepiped (tria, n_subdivisions_, + corners, + colorize); } template void GridGenerator::subdivided_parallelepiped (Triangulation &tria, - const unsigned int ( n_subdivisions) [dim], - const Point (&corners) [dim], - const bool colorize) + const unsigned int ( n_subdivisions) [dim], + const Point (&corners) [dim], + const bool colorize) { // Zero n_subdivisions is the origin only, which makes no sense for (unsigned int i=0; i &tria, { case 1: for (unsigned int x=0; x<=n_subdivisions[0]; ++x) - points.push_back (Point (x*delta[0])); + points.push_back (Point (x*delta[0])); break; - + case 2: for (unsigned int y=0; y<=n_subdivisions[1]; ++y) - for (unsigned int x=0; x<=n_subdivisions[0]; ++x) - points.push_back (Point (x*delta[0] + y*delta[1])); + for (unsigned int x=0; x<=n_subdivisions[0]; ++x) + points.push_back (Point (x*delta[0] + y*delta[1])); break; case 3: for (unsigned int z=0; z<=n_subdivisions[2]; ++z) - for (unsigned int y=0; y<=n_subdivisions[1]; ++y) - for (unsigned int x=0; x<=n_subdivisions[0]; ++x) - points.push_back (Point (x*delta[0] + y*delta[1] + z*delta[2])); + for (unsigned int y=0; y<=n_subdivisions[1]; ++y) + for (unsigned int x=0; x<=n_subdivisions[0]; ++x) + points.push_back (Point (x*delta[0] + y*delta[1] + z*delta[2])); break; default: @@ -554,19 +554,19 @@ GridGenerator::subdivided_parallelepiped (Triangulation &tria, for (unsigned int i=0; i > cells (n_cells); - // Create fixed ordering of + // Create fixed ordering of switch (dim) { case 1: @@ -575,60 +575,60 @@ GridGenerator::subdivided_parallelepiped (Triangulation &tria, cells[x].vertices[0] = x; cells[x].vertices[1] = x+1; - // wipe material id + // wipe material id cells[x].material_id = 0; } break; case 2: - { - // Shorthand - const unsigned int n_dy = n_subdivisions[1]; - const unsigned int n_dx = n_subdivisions[0]; - - for (unsigned int y=0; y &tria, dof_handler.distribute_dofs(q1); CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (), - dof_handler.n_dofs ()); + dof_handler.n_dofs ()); DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern); c_sparsity_pattern.compress (); diff --git a/deal.II/source/lac/petsc_matrix_base.cc b/deal.II/source/lac/petsc_matrix_base.cc index a49bc7cb43..c8a3d456c1 100644 --- a/deal.II/source/lac/petsc_matrix_base.cc +++ b/deal.II/source/lac/petsc_matrix_base.cc @@ -602,7 +602,7 @@ namespace PETScWrappers // Set options PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD, - format); + format); // Write to screen MatView (matrix, PETSC_VIEWER_STDOUT_WORLD); diff --git a/deal.II/source/lac/petsc_parallel_vector.cc b/deal.II/source/lac/petsc_parallel_vector.cc index 612be27193..888f451a5a 100644 --- a/deal.II/source/lac/petsc_parallel_vector.cc +++ b/deal.II/source/lac/petsc_parallel_vector.cc @@ -322,28 +322,28 @@ namespace PETScWrappers out.setf (std::ios::fixed, std::ios::floatfield); for ( unsigned int i = 0; - i < Utilities::MPI::n_mpi_processes(communicator); - i++) - { - // This is slow, but most likely only used to debug. - MPI_Barrier(communicator); - if (i == Utilities::MPI::this_mpi_process(communicator)) - { - if (across) - { - out << "[Proc" << i << " " << istart << "-" << iend-1 << "]" << ' '; - for (PetscInt i=0; ieps); - AssertThrow (ierr == 0, ExcSLEPcError(ierr)); - - // set eigenspectrum problem type (general/standard) - AssertThrow (opA, ExcSLEPcWrappersUsageError()); - if (opB) - ierr = EPSSetOperators (solver_data->eps, *opA, *opB); - else - ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL); - AssertThrow (ierr == 0, ExcSLEPcError(ierr)); - - // set runtime options - set_solver_type (solver_data->eps); + // reset solver dtaa + solver_data.reset (new SolverData()); + + // create eigensolver context and set operators + ierr = EPSCreate (mpi_communicator, &solver_data->eps); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // set eigenspectrum problem type (general/standard) + AssertThrow (opA, ExcSLEPcWrappersUsageError()); + if (opB) + ierr = EPSSetOperators (solver_data->eps, *opA, *opB); + else + ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // set runtime options + set_solver_type (solver_data->eps); } // set the initial vector(s) if any @@ -153,14 +153,14 @@ namespace SLEPcWrappers AssertThrow (ierr == 0, ExcSLEPcError(ierr)); // Set convergence test to be absolute - ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS); - AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); // Set the convergence test function - // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test, - // reinterpret_cast(&solver_control)); - // AssertThrow (ierr == 0, ExcSLEPcError(ierr)); - + // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test, + // reinterpret_cast(&solver_control)); + // AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // solve the eigensystem ierr = EPSSolve (solver_data->eps); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); @@ -170,9 +170,9 @@ namespace SLEPcWrappers #ifdef PETSC_USE_64BIT_INDICES reinterpret_cast(n_converged) #else - reinterpret_cast(n_converged) + reinterpret_cast(n_converged) #endif - ); + ); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); int n_iterations = 0; @@ -182,24 +182,24 @@ namespace SLEPcWrappers // complete eigenspectrum { // get the number of solver iterations - ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations); - AssertThrow (ierr == 0, ExcSLEPcError(ierr)); - + ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // get the residual norm of the most extreme eigenvalue ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm); - AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); // check the solver state const SolverControl::State state - = solver_control.check (n_iterations, residual_norm); + = solver_control.check (n_iterations, residual_norm); // get the solver state according to SLEPc get_solver_state (state); - + // and in case of failure: throw exception - if (solver_control.last_check () != SolverControl::success) - throw SolverControl::NoConvergence (solver_control.last_step (), - solver_control.last_value ()); + if (solver_control.last_check () != SolverControl::success) + throw SolverControl::NoConvergence (solver_control.last_step (), + solver_control.last_value ()); } } @@ -234,7 +234,7 @@ namespace SLEPcWrappers { if (solver_data.get () == 0) return NULL; - + return &solver_data->eps; } @@ -244,22 +244,22 @@ namespace SLEPcWrappers switch (state) { case ::dealii::SolverControl::iterate: - solver_data->reason = EPS_CONVERGED_ITERATING; - break; - + solver_data->reason = EPS_CONVERGED_ITERATING; + break; + case ::dealii::SolverControl::success: - solver_data->reason = static_cast(1); - break; - + solver_data->reason = static_cast(1); + break; + case ::dealii::SolverControl::failure: - if (solver_control.last_step() > solver_control.max_steps()) - solver_data->reason = EPS_DIVERGED_ITS; - else - solver_data->reason = EPS_DIVERGED_BREAKDOWN; - break; - + if (solver_control.last_step() > solver_control.max_steps()) + solver_data->reason = EPS_DIVERGED_ITS; + else + solver_data->reason = EPS_DIVERGED_BREAKDOWN; + break; + default: - Assert (false, ExcNotImplemented()); + Assert (false, ExcNotImplemented()); } } @@ -272,10 +272,10 @@ namespace SLEPcWrappers int SolverBase::convergence_test (EPS /*eps */, - PetscScalar /*kr */, - PetscScalar /*ki */, - PetscReal /*residual_norm */, - PetscReal */*estimated_error */, + PetscScalar /*kr */, + PetscScalar /*ki */, + PetscReal /*residual_norm */, + PetscReal */*estimated_error */, void */*solver_control_x*/) { // This function is undefined (future reference only). @@ -392,8 +392,8 @@ namespace SLEPcWrappers /* ---------------- Generalized Davidson ----------------- */ SolverGeneralizedDavidson::SolverGeneralizedDavidson (SolverControl &cn, - const MPI_Comm &mpi_communicator, - const AdditionalData &data) + const MPI_Comm &mpi_communicator, + const AdditionalData &data) : SolverBase (cn, mpi_communicator), additional_data (data) @@ -425,8 +425,8 @@ namespace SLEPcWrappers /* ------------------ Jacobi Davidson -------------------- */ SolverJacobiDavidson::SolverJacobiDavidson (SolverControl &cn, - const MPI_Comm &mpi_communicator, - const AdditionalData &data) + const MPI_Comm &mpi_communicator, + const AdditionalData &data) : SolverBase (cn, mpi_communicator), additional_data (data) diff --git a/deal.II/source/lac/trilinos_sparse_matrix.cc b/deal.II/source/lac/trilinos_sparse_matrix.cc index 821d2fbac4..5442c86c27 100644 --- a/deal.II/source/lac/trilinos_sparse_matrix.cc +++ b/deal.II/source/lac/trilinos_sparse_matrix.cc @@ -588,13 +588,13 @@ namespace TrilinosWrappers dealii_sparse_matrix.get_sparsity_pattern(); if (matrix.get() == 0 || - m() != n_rows || + m() != n_rows || n_nonzero_elements() != sparsity_pattern.n_nonzero_elements()) - { - SparsityPattern trilinos_sparsity; - trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern); - reinit (trilinos_sparsity); - } + { + SparsityPattern trilinos_sparsity; + trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern); + reinit (trilinos_sparsity); + } // fill the values. the same as above: go through all rows of the matrix, // and then all columns. since the sparsity patterns of the input matrix @@ -619,10 +619,10 @@ namespace TrilinosWrappers // optimized diagonal AssertDimension(it->column(), row); if (std::fabs(it->value()) > drop_tolerance) - { - values[col] = it->value(); - row_indices[col++] = it->column(); - } + { + values[col] = it->value(); + row_indices[col++] = it->column(); + } ++select_index; ++it; } diff --git a/deal.II/source/multigrid/mg_tools.cc b/deal.II/source/multigrid/mg_tools.cc index bd45b98853..331add324f 100644 --- a/deal.II/source/multigrid/mg_tools.cc +++ b/deal.II/source/multigrid/mg_tools.cc @@ -1636,7 +1636,7 @@ namespace MGTools std::map boundary_values; for (std::set::const_iterator p=boundary_dofs.begin(); - p != boundary_dofs.end(); ++p) + p != boundary_dofs.end(); ++p) boundary_values[*p] = 0; Vector dummy(matrix.m()); @@ -1674,7 +1674,7 @@ namespace MGTools std::map boundary_values; for (std::set::const_iterator p=boundary_dofs.begin(); - p != boundary_dofs.end(); ++p) + p != boundary_dofs.end(); ++p) boundary_values[*p] = 0; BlockVector dummy(matrix.n_block_rows()); diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc index 5ab18c9bcf..d6516f5435 100644 --- a/deal.II/source/numerics/matrix_tools.cc +++ b/deal.II/source/numerics/matrix_tools.cc @@ -1988,11 +1988,11 @@ namespace MatrixTools // set entries of this line to zero except for the diagonal // entry - for (typename SparseMatrix::iterator - p = matrix.begin(dof_number); - p != matrix.end(dof_number); ++p) - if (p->column() != dof_number) - p->value() = 0.; + for (typename SparseMatrix::iterator + p = matrix.begin(dof_number); + p != matrix.end(dof_number); ++p) + if (p->column() != dof_number) + p->value() = 0.; // set right hand side to // wanted value: if main diagonal @@ -2045,8 +2045,8 @@ namespace MatrixTools // since that is the diagonal element and thus the present // row for (typename SparseMatrix::iterator - q = matrix.begin(dof_number)+1; - q != matrix.end(dof_number); ++q) + q = matrix.begin(dof_number)+1; + q != matrix.end(dof_number); ++q) { const unsigned int row = q->column(); @@ -2055,7 +2055,7 @@ namespace MatrixTools // (row,dof_number) bool (*comp)(const typename SparseMatrix::iterator::value_type p, const unsigned int column) - = &column_less_than::iterator>; + = &column_less_than::iterator>; const typename SparseMatrix::iterator p = Utilities::lower_bound(matrix.begin(row)+1, matrix.end(row), @@ -2289,10 +2289,10 @@ namespace MatrixTools // the diagonal block for (typename SparseMatrix::iterator q = (block_index.first == block_row ? - transpose_matrix.begin(block_index.second)+1 : - transpose_matrix.begin(block_index.second)); - q != transpose_matrix.end(block_index.second); - ++q) + transpose_matrix.begin(block_index.second)+1 : + transpose_matrix.begin(block_index.second)); + q != transpose_matrix.end(block_index.second); + ++q) { // get the number // of the column in @@ -2320,7 +2320,7 @@ namespace MatrixTools // sub-matrices bool (*comp)(const typename SparseMatrix::iterator::value_type p, const unsigned int column) - = &column_less_than::iterator>; + = &column_less_than::iterator>; typename SparseMatrix::iterator p = this_matrix.end(); -- 2.39.5