From: heister Date: Mon, 3 Dec 2012 19:24:30 +0000 (+0000) Subject: merge from upstream X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c4d173a3554cec5488f223c91954b8401a024e7d;p=dealii-svn.git merge from upstream git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27733 0785d39b-7218-0410-832d-ea1e28bc413d --- c4d173a3554cec5488f223c91954b8401a024e7d diff --cc deal.II/examples/step-13/step-13.cc index 63e62cc64b,b85cb3d50a..b180dc8325 --- a/deal.II/examples/step-13/step-13.cc +++ b/deal.II/examples/step-13/step-13.cc @@@ -292,32 -203,19 +203,19 @@@ namespace Step1 void PointValueEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { - // First allocate a variable that - // will hold the point - // value. Initialize it with a - // value that is clearly bogus, - // so that if we fail to set it - // to a reasonable value, we will - // note at once. This may not be - // necessary in a function as - // small as this one, since we - // can easily see all possible - // paths of execution here, but - // it proved to be helpful for - // more complex cases, and so we - // employ this strategy here as - // well. + // First allocate a variable that will hold the point value. Initialize + // it with a value that is clearly bogus, so that if we fail to set it + // to a reasonable value, we will note at once. This may not be + // necessary in a function as small as this one, since we can easily see + // all possible paths of execution here, but it proved to be helpful for + // more complex cases, and so we employ this strategy here as well. double point_value = 1e20; - // Then loop over all cells and - // all their vertices, and check - // whether a vertex matches the - // evaluation point. If this is - // the case, then extract the - // point value, set a flag that - // we have found the point of + // Then loop over all cells and all their vertices, and check whether a + // vertex matches the evaluation point. If this is the case, then + // extract the point value, set a flag that we have found the point of // interest, and exit the loop. typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), diff --cc deal.II/examples/step-14/step-14.cc index c1d4613872,6242b228e3..3dc49df692 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@@ -216,19 -198,15 +198,15 @@@ namespace Step1 void PointXDerivativeEvaluation:: operator () (const DoFHandler &dof_handler, - const Vector &solution) const + const Vector &solution) const { - // This time initialize the - // return value with something - // useful, since we will have to - // add up a number of - // contributions and take the - // mean value afterwards... + // This time initialize the return value with something useful, since we + // will have to add up a number of contributions and take the mean value + // afterwards... double point_derivative = 0; - // ...then have some objects of - // which the meaning wil become - // clear below... + // ...then have some objects of which the meaning wil become clear + // below... QTrapez vertex_quadrature; FEValues fe_values (dof_handler.get_fe(), vertex_quadrature, diff --cc deal.II/examples/step-20/step-20.cc index fd816b8be1,49bb88bd09..c7c9fca010 --- a/deal.II/examples/step-20/step-20.cc +++ b/deal.II/examples/step-20/step-20.cc @@@ -191,15 -150,11 +150,11 @@@ namespace Step2 }; - // And then we also have to define - // these respective functions, of - // course. Given our discussion in - // the introduction of how the - // solution should look like, the - // following computations should be - // straightforward: + // And then we also have to define these respective functions, of + // course. Given our discussion in the introduction of how the solution + // should look like, the following computations should be straightforward: template - double RightHandSide::value (const Point & /*p*/, + double RightHandSide::value (const Point &/*p*/, const unsigned int /*component*/) const { return 0; diff --cc deal.II/examples/step-28/step-28.cc index 0363c2dad6,dc4a31b387..049d63568e --- a/deal.II/examples/step-28/step-28.cc +++ b/deal.II/examples/step-28/step-28.cc @@@ -1638,21 -1219,13 +1219,13 @@@ namespace Step2 // @sect5{Private member variables} - // Next, we have a few member - // variables. In particular, - // these are (i) a reference to - // the parameter object (owned by - // the main function of this - // program, and passed to the - // constructor of this class), - // (ii) an object describing the - // material parameters for the - // number of energy groups - // requested in the input file, - // and (iii) the finite element - // to be used by all energy - // groups: + // Next, we have a few member variables. In particular, these are (i) a + // reference to the parameter object (owned by the main function of this + // program, and passed to the constructor of this class), (ii) an object + // describing the material parameters for the number of energy groups + // requested in the input file, and (iii) the finite element to be used by + // all energy groups: - const Parameters ¶meters; + const Parameters ¶meters; const MaterialData material_data; FE_Q fe; diff --cc deal.II/examples/step-29/step-29.cc index 4e4bff6dd9,253a0c92f6..b97a0558b7 --- a/deal.II/examples/step-29/step-29.cc +++ b/deal.II/examples/step-29/step-29.cc @@@ -572,16 -400,12 +400,12 @@@ namespace Step2 - // The constructor takes the - // ParameterHandler object and stores - // it in a reference. It also - // initializes the DoF-Handler and - // the finite element system, which - // consists of two copies of the - // scalar Q1 field, one for $v$ and - // one for $w$: + // The constructor takes the ParameterHandler object and stores it in a + // reference. It also initializes the DoF-Handler and the finite element + // system, which consists of two copies of the scalar Q1 field, one for $v$ + // and one for $w$: template - UltrasoundProblem::UltrasoundProblem (ParameterHandler ¶m) + UltrasoundProblem::UltrasoundProblem (ParameterHandler ¶m) : prm(param), dof_handler(triangulation), diff --cc deal.II/examples/step-33/step-33.cc index 1721b8b098,77199f884b..6b4b5cd418 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@@ -765,26 -587,19 +587,19 @@@ namespace Step3 EulerEquations::Postprocessor:: compute_derived_quantities_vector (const std::vector > &uh, const std::vector > > &duh, - const std::vector > > & /*dduh*/, - const std::vector > & /*normals*/, - const std::vector > & /*evaluation_points*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &/*evaluation_points*/, std::vector > &computed_quantities) const { - // At the beginning of the function, let us - // make sure that all variables have the - // correct sizes, so that we can access - // individual vector elements without - // having to wonder whether we might read - // or write invalid elements; we also check - // that the duh vector only - // contains data if we really need it (the - // system knows about this because we say - // so in the - // get_needed_update_flags() - // function below). For the inner vectors, - // we check that at least the first element - // of the outer vector has the correct - // inner size: + // At the beginning of the function, let us make sure that all variables + // have the correct sizes, so that we can access individual vector + // elements without having to wonder whether we might read or write + // invalid elements; we also check that the duh vector only + // contains data if we really need it (the system knows about this because + // we say so in the get_needed_update_flags() function + // below). For the inner vectors, we check that at least the first element + // of the outer vector has the correct inner size: const unsigned int n_quadrature_points = uh.size(); if (do_schlieren_plot == true) diff --cc deal.II/examples/step-4/step-4.cc index 81fa0ad715,7aa256fa80..a39f1010f9 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@@ -437,21 -320,15 +320,15 @@@ void Step4::assemble_system ( FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - // Next, we again have to loop over all - // cells and assemble local contributions. - // Note, that a cell is a quadrilateral in - // two space dimensions, but a hexahedron - // in 3D. In fact, the - // active_cell_iterator data - // type is something different, depending - // on the dimension we are in, but to the - // outside world they look alike and you - // will probably never see a difference - // although the classes that this typedef - // stands for are in fact completely - // unrelated: + // Next, we again have to loop over all cells and assemble local + // contributions. Note, that a cell is a quadrilateral in two space + // dimensions, but a hexahedron in 3D. In fact, the + // active_cell_iterator data type is something different, + // depending on the dimension we are in, but to the outside world they look + // alike and you will probably never see a difference although the classes + // that this typedef stands for are in fact completely unrelated: typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); diff --cc deal.II/examples/step-5/step-5.cc index 0b059c108d,0430781c5c..7b36606830 --- a/deal.II/examples/step-5/step-5.cc +++ b/deal.II/examples/step-5/step-5.cc @@@ -442,61 -304,33 +304,33 @@@ void Step5::assemble_system ( FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); - // Here is one difference: for this - // program, we will again use a - // constant right hand side - // function and zero boundary - // values, but a variable - // coefficient. We have already - // declared the class that - // represents this coefficient - // above, so we only have to - // declare a corresponding object + // Here is one difference: for this program, we will again use a constant + // right hand side function and zero boundary values, but a variable + // coefficient. We have already declared the class that represents this + // coefficient above, so we only have to declare a corresponding object // here. // - // Then, below, we will ask the - // coefficient function object - // to compute the values of the - // coefficient at all quadrature - // points on one cell at once. The - // reason for this is that, if you - // look back at how we did this in - // step-4, you will realize that we - // called the function computing - // the right hand side value inside - // nested loops over all degrees of - // freedom and over all quadrature - // points, - // i.e. dofs_per_cell*n_q_points - // times. For the coefficient that - // is used inside the matrix, this - // would actually be - // dofs_per_cell*dofs_per_cell*n_q_points. On - // the other hand, the function - // will of course return the same - // value every time it is called - // with the same quadrature point, - // independently of what shape - // function we presently treat; - // secondly, these are virtual - // function calls, so are rather - // expensive. Obviously, there are - // only n_q_point different values, - // and we shouldn't call the - // function more often than - // that. Or, even better than this, - // compute all of these values at - // once, and get away with a single + // Then, below, we will ask the coefficient function object to + // compute the values of the coefficient at all quadrature points on one + // cell at once. The reason for this is that, if you look back at how we did + // this in step-4, you will realize that we called the function computing + // the right hand side value inside nested loops over all degrees of freedom + // and over all quadrature points, i.e. dofs_per_cell*n_q_points times. For + // the coefficient that is used inside the matrix, this would actually be + // dofs_per_cell*dofs_per_cell*n_q_points. On the other hand, the function + // will of course return the same value every time it is called with the + // same quadrature point, independently of what shape function we presently + // treat; secondly, these are virtual function calls, so are rather + // expensive. Obviously, there are only n_q_point different values, and we + // shouldn't call the function more often than that. Or, even better than + // this, compute all of these values at once, and get away with a single // function call per cell. // - // This is exactly what we are - // going to do. For this, we need - // some space to store the values - // in. We therefore also have to - // declare an array to hold these - // values: + // This is exactly what we are going to do. For this, we need some space to + // store the values in. We therefore also have to declare an array to hold + // these values: const Coefficient coefficient; std::vector coefficient_values (n_q_points); diff --cc deal.II/examples/step-9/step-9.cc index 3ae105a770,34c1acbb41..25b77c3c2b --- a/deal.II/examples/step-9/step-9.cc +++ b/deal.II/examples/step-9/step-9.cc @@@ -1459,23 -976,15 +976,15 @@@ namespace Step template void GradientEstimation::estimate (const DoFHandler &dof_handler, - const Vector &solution, + const Vector &solution, Vector &error_per_cell) { - // Before starting with the work, - // we check that the vector into - // which the results are written, - // has the right size. It is a - // common error that such - // parameters have the wrong size, - // but the resulting damage by not - // catching these errors are very - // subtle as they are usually - // corruption of data somewhere in - // memory. Often, the problems - // emerging from this are not - // reproducible, and we found that - // it is well worth the effort to + // Before starting with the work, we check that the vector into which the + // results are written, has the right size. It is a common error that such + // parameters have the wrong size, but the resulting damage by not + // catching these errors are very subtle as they are usually corruption of + // data somewhere in memory. Often, the problems emerging from this are + // not reproducible, and we found that it is well worth the effort to // check for such things. Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(), ExcInvalidVectorLength (error_per_cell.size(), diff --cc deal.II/include/deal.II/base/exceptions.h index ed9b1615e1,822a0f7fe7..d74bbb592b --- a/deal.II/include/deal.II/base/exceptions.h +++ b/deal.II/include/deal.II/base/exceptions.h @@@ -155,12 -134,11 +134,11 @@@ protected /** * Name of the exception and call sequence. */ - const char *exc; + const char *exc; /** - * A backtrace to the position - * where the problem happened, if - * the system supports this. + * A backtrace to the position where the problem happened, if the + * system supports this. */ char **stacktrace; diff --cc deal.II/include/deal.II/lac/chunk_sparse_matrix.h index b0bd665389,7f47346c8e..bfeeed0804 --- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h @@@ -374,10 -374,35 +374,35 @@@ public * is allowed to store zero * values in non-existent fields. */ - void add (const unsigned int i, - const unsigned int j, + void add (const types::global_dof_index i, + const types::global_dof_index j, const number value); + /** + * Add an array of values given by + * values in the given + * global matrix row at columns + * specified by col_indices in the + * sparse matrix. + * + * The optional parameter + * elide_zero_values can be + * used to specify whether zero + * values should be added anyway or + * these should be filtered away and + * only non-zero data is added. The + * default value is true, + * i.e., zero values won't be added + * into the matrix. + */ + template + void add (const unsigned int row, + const unsigned int n_cols, + const unsigned int *col_indices, + const number2 *values, + const bool elide_zero_values = true, + const bool col_indices_are_sorted = false); + /** * Multiply the entire matrix by a * fixed factor. diff --cc deal.II/source/dofs/dof_handler.cc index 721d0713c0,a91f82f6c2..923de03c7a --- a/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/source/dofs/dof_handler.cc @@@ -341,381 -336,19 +336,385 @@@ namespace interna } dof_handler.faces = new internal::DoFHandler::DoFFaces<3>; - dof_handler.faces->lines.dofs - .resize (dof_handler.tria->n_raw_lines() * - dof_handler.selected_fe->dofs_per_line, - DoFHandler<3,spacedim>::invalid_dof_index); - dof_handler.faces->quads.dofs - .resize (dof_handler.tria->n_raw_quads() * - dof_handler.selected_fe->dofs_per_quad, - DoFHandler<3,spacedim>::invalid_dof_index); + // avoid access to n_raw_lines when there are no cells + if (dof_handler.tria->n_cells() > 0) + { + dof_handler.faces->lines.dofs + .resize (dof_handler.tria->n_raw_lines() * + dof_handler.selected_fe->dofs_per_line, + DoFHandler<3,spacedim>::invalid_dof_index); + dof_handler.faces->quads.dofs + .resize (dof_handler.tria->n_raw_quads() * + dof_handler.selected_fe->dofs_per_quad, + DoFHandler<3,spacedim>::invalid_dof_index); + } } + + template + static + void reserve_space_mg (DoFHandler<1, spacedim> &dof_handler) + { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); + + const dealii::Triangulation<1, spacedim> &tria = dof_handler.get_tria (); + const unsigned int &dofs_per_line = dof_handler.get_fe ().dofs_per_line; + const unsigned int &n_levels = tria.n_levels (); + + for (unsigned int i = 0; i < n_levels; ++i) + { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<1>); + dof_handler.mg_levels.back ()->dof_object.dofs = std::vector (tria.n_raw_lines (i) * dofs_per_line, DoFHandler<1>::invalid_dof_index); + } + + const unsigned int &n_vertices = tria.n_vertices (); + + dof_handler.mg_vertex_dofs.resize (n_vertices); + + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); + + for (typename dealii::Triangulation<1, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) + { + const unsigned int level = cell->level (); + + for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) + { + const unsigned int vertex_index = cell->vertex_index (vertex); + + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; + + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } + + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) + { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], dof_handler.get_fe ().dofs_per_vertex); + } + + else + { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } + + template + static + void reserve_space_mg (DoFHandler<2, spacedim> &dof_handler) + { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); + + const dealii::FiniteElement<2, spacedim> &fe = dof_handler.get_fe (); + const dealii::Triangulation<2, spacedim> &tria = dof_handler.get_tria (); + const unsigned int &n_levels = tria.n_levels (); + + for (unsigned int i = 0; i < n_levels; ++i) + { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<2>); + dof_handler.mg_levels.back ()->dof_object.dofs = std::vector (tria.n_raw_quads (i) * fe.dofs_per_quad, DoFHandler<2>::invalid_dof_index); + } + + dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<2>; + dof_handler.mg_faces->lines.dofs = std::vector (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<2>::invalid_dof_index); + + const unsigned int &n_vertices = tria.n_vertices (); + + dof_handler.mg_vertex_dofs.resize (n_vertices); + + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); + + for (typename dealii::Triangulation<2, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) + { + const unsigned int level = cell->level (); + + for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) + { + const unsigned int vertex_index = cell->vertex_index (vertex); + + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; + + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } + + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) + { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex); + } + + else + { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } + + template + static + void reserve_space_mg (DoFHandler<3, spacedim> &dof_handler) + { + Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation")); + dof_handler.clear_mg_space (); + + const dealii::FiniteElement<3, spacedim> &fe = dof_handler.get_fe (); + const dealii::Triangulation<3, spacedim> &tria = dof_handler.get_tria (); + const unsigned int &n_levels = tria.n_levels (); + + for (unsigned int i = 0; i < n_levels; ++i) + { + dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<3>); + dof_handler.mg_levels.back ()->dof_object.dofs = std::vector (tria.n_raw_hexs (i) * fe.dofs_per_hex, DoFHandler<3>::invalid_dof_index); + } + + dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<3>; + dof_handler.mg_faces->lines.dofs = std::vector (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index); + dof_handler.mg_faces->quads.dofs = std::vector (tria.n_raw_quads () * fe.dofs_per_quad, DoFHandler<3>::invalid_dof_index); + + const unsigned int &n_vertices = tria.n_vertices (); + + dof_handler.mg_vertex_dofs.resize (n_vertices); + + std::vector max_level (n_vertices, 0); + std::vector min_level (n_vertices, n_levels); + + for (typename dealii::Triangulation<3, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) + { + const unsigned int level = cell->level (); + + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) + { + const unsigned int vertex_index = cell->vertex_index (vertex); + + if (min_level[vertex_index] > level) + min_level[vertex_index] = level; + + if (max_level[vertex_index] < level) + max_level[vertex_index] = level; + } + } + + for (unsigned int vertex = 0; vertex < n_vertices; ++vertex) + if (tria.vertex_used (vertex)) + { + Assert (min_level[vertex] < n_levels, ExcInternalError ()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex); + } + + else + { + Assert (min_level[vertex] == n_levels, ExcInternalError ()); + Assert (max_level[vertex] == 0, ExcInternalError ()); + dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0); + } + } + + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<1, spacedim>::cell_iterator &cell, types::global_dof_index next_free_dof) + { + const FiniteElement<1, spacedim> &fe = cell->get_fe (); + + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) + { + typename DoFHandler<1, spacedim>::cell_iterator neighbor = cell->neighbor (vertex); + + if (neighbor.state () == IteratorState::valid) + if (neighbor->user_flag_set () && (neighbor->level () == cell->level ())) + { + if (vertex == 0) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), 0, dof, neighbor->mg_vertex_dof_index (cell->level (), 1, dof)); + + else + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), 1, dof, neighbor->mg_vertex_dof_index (cell->level (), 0, dof)); + + continue; + } + + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); + } + + if (fe.dofs_per_line > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); + + cell->set_user_flag (); + return next_free_dof; + } + + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<2, spacedim>::cell_iterator &cell, types::global_dof_index next_free_dof) + { + const FiniteElement<2, spacedim> &fe = cell->get_fe (); + + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) + if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<2>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); + + if (fe.dofs_per_line > 0) + for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face) + { + typename DoFHandler<2, spacedim>::line_iterator line = cell->line (face); + + if (line->mg_dof_index (cell->level (), 0) == DoFHandler<2>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + line->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } + + if (fe.dofs_per_quad > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); + + cell->set_user_flag (); + return next_free_dof; + } + + template + static + types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<3, spacedim>::cell_iterator &cell, types::global_dof_index next_free_dof) + { + const FiniteElement<3, spacedim> &fe = cell->get_fe (); + + if (fe.dofs_per_vertex > 0) + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) + if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof) + cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++); + + if (fe.dofs_per_line > 0) + for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line) + { + typename DoFHandler<3, spacedim>::line_iterator line_it = cell->line (line); + + if (line_it->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof) + line_it->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } + + if (fe.dofs_per_quad > 0) + for (unsigned int face = 0; face < GeometryInfo<3>::quads_per_cell; ++face) + { + typename DoFHandler<3, spacedim>::quad_iterator quad = cell->quad (face); + + if (quad->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index) + for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof) + quad->set_mg_dof_index (cell->level (), dof, next_free_dof++); + } + + if (fe.dofs_per_hex > 0) + for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof) + cell->set_mg_dof_index (cell->level (), dof, next_free_dof++); + + cell->set_user_flag (); + return next_free_dof; + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<1, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<1> &mg_level, internal::DoFHandler::DoFFaces<1> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) + { + return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &, internal::DoFHandler::DoFFaces<2> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) + { + return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &mg_level, internal::DoFHandler::DoFFaces<2> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) + { + return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) + { + return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) + { + return mg_faces.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + types::global_dof_index get_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &mg_level, internal::DoFHandler::DoFFaces<3> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<3>) + { + return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index); + } + + template + static + void set_dof_index (const DoFHandler<1, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<1> &mg_level, internal::DoFHandler::DoFFaces<1> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) + { + mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } + + template + static + void set_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &, internal::DoFHandler::DoFFaces<2> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) + { + mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } + + template + static + void set_dof_index (const DoFHandler<2, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<2> &mg_level, internal::DoFHandler::DoFFaces<2> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) + { + mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } + + template + static + void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) + { + mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } + + template + static + void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &, internal::DoFHandler::DoFFaces<3> &mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) + { + mg_faces.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } + + template + static + void set_dof_index (const DoFHandler<3, spacedim> &dof_handler, internal::DoFHandler::DoFLevel<3> &mg_level, internal::DoFHandler::DoFFaces<3> &, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<3>) + { + mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index); + } }; } } diff --cc deal.II/source/dofs/dof_handler_policy.cc index ee8bdf5b73,8d6aa68f8f..0c04c9cc9e --- a/deal.II/source/dofs/dof_handler_policy.cc +++ b/deal.II/source/dofs/dof_handler_policy.cc @@@ -212,14 -212,10 +212,10 @@@ namespace interna /** - * Distribute degrees of - * freedom on all cells, or - * on cells with the - * correct subdomain_id if - * the corresponding - * argument is not equal to - * types::invalid_subdomain_id. Return - * the next free dof index. + * Distribute degrees of freedom on all cells, or on cells with the + * correct subdomain_id if the corresponding argument is not equal to + * types::invalid_subdomain_id. Return the total number of dofs - * returned. ++ * distributed. */ template static