]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
merge from upstream
authorheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Dec 2012 19:24:30 +0000 (19:24 +0000)
committerheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Dec 2012 19:24:30 +0000 (19:24 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27733 0785d39b-7218-0410-832d-ea1e28bc413d

39 files changed:
1  2 
deal.II/examples/step-13/step-13.cc
deal.II/examples/step-14/step-14.cc
deal.II/examples/step-18/step-18.cc
deal.II/examples/step-20/step-20.cc
deal.II/examples/step-21/step-21.cc
deal.II/examples/step-22/step-22.cc
deal.II/examples/step-23/step-23.cc
deal.II/examples/step-28/step-28.cc
deal.II/examples/step-29/step-29.cc
deal.II/examples/step-31/step-31.cc
deal.II/examples/step-32/step-32.cc
deal.II/examples/step-33/step-33.cc
deal.II/examples/step-35/step-35.cc
deal.II/examples/step-37/step-37.cc
deal.II/examples/step-4/step-4.cc
deal.II/examples/step-42/step-42.cc
deal.II/examples/step-43/step-43.cc
deal.II/examples/step-46/step-46.cc
deal.II/examples/step-47/step-47.cc
deal.II/examples/step-5/step-5.cc
deal.II/examples/step-6/step-6.cc
deal.II/examples/step-9/step-9.cc
deal.II/include/deal.II/base/exceptions.h
deal.II/include/deal.II/base/types.h
deal.II/include/deal.II/dofs/dof_accessor.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/include/deal.II/grid/tria_objects.h
deal.II/include/deal.II/lac/chunk_sparse_matrix.h
deal.II/include/deal.II/lac/petsc_solver.h
deal.II/include/deal.II/lac/sparse_matrix.h
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/dofs/dof_handler.cc
deal.II/source/dofs/dof_handler_policy.cc
deal.II/source/dofs/dof_tools.cc
deal.II/source/dofs/dof_tools.inst.in
deal.II/source/fe/fe_q.cc
deal.II/source/lac/constraint_matrix.cc
deal.II/source/lac/petsc_solver.cc

index 63e62cc64bae6f2e2f4cf446b644e2708646b422,b85cb3d50a0a129499e1ee0b7ecda1e40aa29ea1..b180dc832547002dd249773488b8cad759b185d4
@@@ -292,32 -203,19 +203,19 @@@ namespace Step1
      void
      PointValueEvaluation<dim>::
      operator () (const DoFHandler<dim> &dof_handler,
 -                 const Vector<double> &solution) const
 +                 const Vector<double>  &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<dim>::active_cell_iterator
        cell = dof_handler.begin_active(),
index c1d4613872548835754f131d47447347976ff445,6242b228e3d37e5ac894894ec8f24c0c8a53be29..3dc49df6928786b9d082f1d50fac157ffff7fe76
@@@ -216,19 -198,15 +198,15 @@@ namespace Step1
      void
      PointXDerivativeEvaluation<dim>::
      operator () (const DoFHandler<dim> &dof_handler,
 -                 const Vector<double> &solution) const
 +                 const Vector<double>  &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<dim>  vertex_quadrature;
        FEValues<dim> fe_values (dof_handler.get_fe(),
                                 vertex_quadrature,
Simple merge
index fd816b8be1a573ee5691e16e898869f6d68dcd48,49bb88bd09c133a432b1072ce7ba4e9f63bc010f..c7c9fca01042a6cd323455b85c437aef351d2cd7
@@@ -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 <int dim>
 -  double RightHandSide<dim>::value (const Point<dim>  & /*p*/,
 +  double RightHandSide<dim>::value (const Point<dim>  &/*p*/,
                                      const unsigned int /*component*/) const
    {
      return 0;
Simple merge
Simple merge
Simple merge
index 0363c2dad69a11b94291487252ab9d42528b81b7,dc4a31b38764ad9ca20814b667f87b37ce17fdbb..049d63568e5c2789781e05f4a73425110413febe
@@@ -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 &parameters;
 +    const Parameters  &parameters;
      const MaterialData material_data;
      FE_Q<dim>          fe;
  
index 4e4bff6dd9bec4f435bf263a287a7aeb097e1062,253a0c92f6e8ad87dc65ab43a59d5c9d076a6995..b97a0558b7bcbf4125278640f1394f0294421846
@@@ -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 <int dim>
 -  UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler &param)
 +  UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler  &param)
      :
      prm(param),
      dof_handler(triangulation),
Simple merge
Simple merge
index 1721b8b098f8f5ecd5ed72bf6386dd0166e84096,77199f884b12b95241ce7375272948d24bea0546..6b4b5cd4182c1111ff4aa84c2069ce95e16d8592
@@@ -765,26 -587,19 +587,19 @@@ namespace Step3
    EulerEquations<dim>::Postprocessor::
    compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
                                       const std::vector<std::vector<Tensor<1,dim> > > &duh,
 -                                     const std::vector<std::vector<Tensor<2,dim> > > & /*dduh*/,
 -                                     const std::vector<Point<dim> >                  & /*normals*/,
 -                                     const std::vector<Point<dim> >                  & /*evaluation_points*/,
 +                                     const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
 +                                     const std::vector<Point<dim> >                  &/*normals*/,
 +                                     const std::vector<Point<dim> >                  &/*evaluation_points*/,
                                       std::vector<Vector<double> >                    &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 <code>duh</code> vector only
-     // contains data if we really need it (the
-     // system knows about this because we say
-     // so in the
-     // <code>get_needed_update_flags()</code>
-     // 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 <code>duh</code> vector only
+     // contains data if we really need it (the system knows about this because
+     // we say so in the <code>get_needed_update_flags()</code> 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)
Simple merge
Simple merge
index 81fa0ad7157ea33c9cafbb7385efea7b15725a6e,7aa256fa808245b86258febcebd1e2ba827992bd..a39f1010f9d5b1a12caee1d7ea7635eaf821704a
@@@ -437,21 -320,15 +320,15 @@@ void Step4<dim>::assemble_system (
    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double>       cell_rhs (dofs_per_cell);
  
 -  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 +  std::vector<types::global_dof_index> 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
-   // <code>active_cell_iterator</code> 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
+   // <code>active_cell_iterator</code> 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<dim>::active_cell_iterator
    cell = dof_handler.begin_active(),
    endc = dof_handler.end();
Simple merge
Simple merge
Simple merge
Simple merge
index 0b059c108d5fed0b0c6a5d92e81a8179653a0142,0430781c5cd93433fb42930b9d6552775661a73a..7b366068304412a0c8af59e1911f604599abe0c1
@@@ -442,61 -304,33 +304,33 @@@ void Step5<dim>::assemble_system (
    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double>       cell_rhs (dofs_per_cell);
  
 -  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 +  std::vector<types::global_dof_index> 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
-   // <code>coefficient</code> 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 <code>coefficient</code> 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<dim> coefficient;
    std::vector<double>    coefficient_values (n_q_points);
  
Simple merge
index 3ae105a770b68b72ee3efb9b99f3125c382840d8,34c1acbb41aeef387753bb3180df14ba4011ceda..25b77c3c2b1b7d05c15205df46e4f830dcd65080
@@@ -1459,23 -976,15 +976,15 @@@ namespace Step
    template <int dim>
    void
    GradientEstimation::estimate (const DoFHandler<dim> &dof_handler,
 -                                const Vector<double> &solution,
 +                                const Vector<double>  &solution,
                                  Vector<float>         &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(),
index ed9b1615e1f708b6f28ff769ec19fb4c4843bc24,822a0f7fe757eb5134fd9d516acec85f53675905..d74bbb592bfd9e99628227acfefbf441d0ed33bf
@@@ -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;
  
Simple merge
index b0bd665389951d2de86302fe2e5facd8d53f5853,7f47346c8e68f0936b38e15a935a6ad0849be87c..bfeeed0804e039e41b6273bcf5908910941cca0a
@@@ -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
+    * <tt>values</tt> in the given
+    * global matrix row at columns
+    * specified by col_indices in the
+    * sparse matrix.
+    *
+    * The optional parameter
+    * <tt>elide_zero_values</tt> 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 <tt>true</tt>,
+    * i.e., zero values won't be added
+    * into the matrix.
+    */
+   template <typename number2>
+   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.
index 721d0713c036d962406c727010c8f44020b85269,a91f82f6c296e84dfd6e06cbec4e452e494d553e..923de03c7a320c288c463ecf88cd52973c87ff60
@@@ -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<int spacedim>
 +      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<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
 +        std::vector<unsigned int> 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<int spacedim>
 +      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<unsigned int> (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<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
 +        std::vector<unsigned int> 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<int spacedim>
 +      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<unsigned int> (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<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index);
 +        dof_handler.mg_faces->quads.dofs = std::vector<unsigned int> (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<unsigned int> max_level (n_vertices, 0);
 +        std::vector<unsigned int> 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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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<int spacedim>
 +      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);
 +      }
      };
    }
  }
index ee8bdf5b73b1e27aa5b0dec8681fde27d33c04e7,8d6aa68f8ff12011d9716eb60426ceb0798d74ae..0c04c9cc9e0ee53fdd6e687f231e5c191abf4f88
@@@ -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 <int dim, int spacedim>
          static
Simple merge
Simple merge
Simple merge
Simple merge

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.