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(),
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,
};
- // 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;
// @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<dim> fe;
- // 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 ¶m)
+ UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler ¶m)
:
prm(param),
dof_handler(triangulation),
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)
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();
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);
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(),
/**
* 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;
* 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.
}
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);
+ }
};
}
}
/**
- * 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
projects / dealii-svn.git / commitdiff