const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
... // assemble local contributions and copy them into global object
- }
+ }
</pre>
</code>
<h3>Complications with linear systems for hp discretizations</h3>
+<h4>Creating the sparsity pattern</h4>
+
One of the problems with $hp$ methods is that the high polynomial degree of
shape functions together with the large number of constrained degrees of
freedom leads to matrices with large numbers of nonzero entries in some
order and $hp$ elements in 3d, the CompressedSetSparsityPattern has a definite
edge. We will therefore use it later when we build the sparsity pattern in
this tutorial program.
+
+
+<h4>Eliminating constrained degrees of freedom</h4>
+
+A second problem particular to $hp$ methods arises because we have so many
+constrained degrees of freedom: typically up to about one third of all degrees
+of freedom (in 3d) are constrained because they either belong to cells with
+hanging nodes or because they are on cells adjacent to cells with a higher or
+lower polynomial degree. This is much more than the fraction of constrained
+degrees of freedom in non-$hp$ mode.
+
+In previous programs, for example in step-6, we have dealt with eliminating
+these constrained degrees of freedom by first building an object of type
+ConstraintMatrix that stores the constraints, and then "condensing" away these
+degrees of freedom first from the sparsity pattern. The same scheme is used
+for the matrix and right hand side, which are also both first built then
+condensed.
+
+Dealing with sparsity patterns and matrices in this way turns out to be
+inefficient because the effort necessary is at least ${\cal O}(N \log N)$ in
+the number of unknowns, whereas an ideal finite element program would of
+course only have algorithms that are linear in the number of unknowns. The
+solution to this problem is to use the ConstraintMatrix object already at the
+time of creating the sparsity pattern, or while copying local contributions
+into the global matrix object (or global vector).
+
+So, instead of the code snippet (taken from @ref step_6 "step-6")
+<code>
+<pre>
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+ hanging_node_constraints.condense (sparsity_pattern);
+</pre>
+</code>
+we now use the following:
+<code>
+<pre>
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
+ hanging_node_constraints);
+</pre>
+</code>
+
+In a similar vein, we have to slightly modify the way we copy local
+contributions into global matrices and vectors. In previous tutorial programs,
+we have always used a process like this:
+<code>
+<pre>
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ ... // assemble local contributions them into global object
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+</pre>
+</code>
+We now replace copying and later condensing into one step, as already shown in
+@ref step_17 "step-17" and @ref step_18 "step-18":
+<code>
+<pre>
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ ... // assemble local contributions them into global object
+
+ cell->get_dof_indices (local_dof_indices);
+
+ hanging_node_constraints
+ .distribute_local_to_global (cell_matrix,
+ local_dof_indices,
+ system_matrix);
+
+ hanging_node_constraints
+ .distribute_local_to_global (cell_rhs,
+ local_dof_indices,
+ system_rhs);
+ }
+</pre>
+</code>
+Essentially, what the
+<code>hanging_node_constraints.distribute_local_to_global</code> calls do is
+to implement the same loop as before, but whenever we hit a constrained degree
+of freedom, the function does the right thing and already condenses it away.
+
+Using this technique of eliminating constrained nodes already when
+transferring local contributions into the global objects, we avoid the problem
+of having to go back later and change these objects. Timing these operations
+shows that this makes the overall algorithms faster.
// directly build the sparsity pattern, but
// first create an intermediate object that
// we later copy into the right data
- // structure. This is as explained in the
- // introduction of this program.
+ // structure. In another slight deviation, we
+ // do not first build the sparsity pattern
+ // then condense away constrained degrees of
+ // freedom, but pass the constraint matrix
+ // object directly to the function that
+ // builds the sparsity pattern. Both of these
+ // steps are explained in the introduction of
+ // this program.
//
// The second change, maybe hidden in plain
// sight, is that the dof_handler variable
// @sect4{LaplaceProblem::assemble_system}
+
+ // This is the function that assembles the
+ // global matrix and right hand side vector
+ // from the local contributions of each
+ // cell. Its main working is as has been
+ // described in many of the tutorial programs
+ // before. The significant deviations are the
+ // ones necessary for $hp$ finite element
+ // methods. In particular, that we need to
+ // use a collection of FEValues object
+ // (implemented through the hp::FEValues
+ // class), and that we have to eliminate
+ // constrained degrees of freedom already
+ // when copying local contributions into
+ // global objects. Both of these are
+ // explained in detail in the introduction of
+ // this program.
template <int dim>
void LaplaceProblem<dim>::assemble_system ()
{
system_rhs);
}
+ // After the steps already discussed above,
+ // all we still have to do is to treat
+ // Dirichlet boundary values
+ // correctly. This works in exactly the
+ // same way as for non-$hp$ objects:
std::map<unsigned int,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
0,
system_rhs);
}
+
+
+ // @sect4{LaplaceProblem::solve}
+
+ // The function solving the linear system is
+ // entirely unchanged from previous
+ // examples. We simply try to reduce the
+ // initial residual (which equals the $l_2$
+ // norm of the right hand side) by a certain
+ // factor:
template <int dim>
void LaplaceProblem<dim>::solve ()
{
-template <int dim>
-void
-LaplaceProblem<dim>::
-estimate_smoothness (Vector<float> &smoothness_indicators) const
-{
- const unsigned int N = max_degree;
-
- // form all the Fourier vectors
- // that we want to
- // consider. exclude k=0 to avoid
- // problems with |k|^{-mu} and also
- // logarithms of |k|
- std::vector<Tensor<1,dim> > k_vectors;
- std::vector<unsigned int> k_vectors_magnitude;
- switch (dim)
- {
- case 2:
- {
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int j=0; j<N; ++j)
- if (!((i==0) && (j==0))
- &&
- (i*i + j*j < N*N))
- {
- k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
- deal_II_numbers::PI * j));
- k_vectors_magnitude.push_back (i*i+j*j);
- }
-
- break;
- }
-
- case 3:
- {
- for (unsigned int i=0; i<N; ++i)
- for (unsigned int j=0; j<N; ++j)
- for (unsigned int k=0; k<N; ++k)
- if (!((i==0) && (j==0) && (k==0))
- &&
- (i*i + j*j + k*k < N*N))
- {
- k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
- deal_II_numbers::PI * j,
- deal_II_numbers::PI * k));
- k_vectors_magnitude.push_back (i*i+j*j+k*k);
- }
-
- break;
- }
-
- default:
- Assert (false, ExcNotImplemented());
- }
-
- const unsigned n_fourier_modes = k_vectors.size();
- std::vector<double> ln_k (n_fourier_modes);
- for (unsigned int i=0; i<n_fourier_modes; ++i)
- ln_k[i] = std::log (k_vectors[i].norm());
-
-
- // assemble the matrices that do
- // the Fourier transforms for each
- // of the finite elements we deal
- // with. note that these matrices
- // are complex-valued, so we can't
- // use FullMatrix
- QGauss<1> base_quadrature (2);
- QIterated<dim> quadrature (base_quadrature, N);
-
- std::vector<Table<2,std::complex<double> > >
- fourier_transform_matrices (fe_collection.size());
- for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
- {
- fourier_transform_matrices[fe].reinit (n_fourier_modes,
- fe_collection[fe].dofs_per_cell);
-
- for (unsigned int k=0; k<n_fourier_modes; ++k)
- for (unsigned int i=0; i<fe_collection[fe].dofs_per_cell; ++i)
- {
- std::complex<double> sum = 0;
- for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
- {
- const Point<dim> x_q = quadrature.point(q);
- sum += std::exp(std::complex<double>(0,1) *
- (k_vectors[k] * x_q)) *
- fe_collection[fe].shape_value(i,x_q) *
- quadrature.weight(q);
- }
- fourier_transform_matrices[fe](k,i)
- = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
- }
- }
-
- // the next thing is to loop over
- // all cells and do our work there,
- // i.e. to locally do the Fourier
- // transform and estimate the decay
- // coefficient
- std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
- Vector<double> local_dof_values;
-
- typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (unsigned int index=0; cell!=endc; ++cell, ++index)
- {
- local_dof_values.reinit (cell->get_fe().dofs_per_cell);
- cell->get_dof_values (solution, local_dof_values);
-
- // first compute the Fourier
- // transform of the local
- // solution
- std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
- for (unsigned int f=0; f<n_fourier_modes; ++f)
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- fourier_coefficients[f] +=
- fourier_transform_matrices[cell->active_fe_index()](f,i)
- *
- local_dof_values(i);
-
- // enter the Fourier
- // coefficients into a map with
- // the k-magnitudes, to make
- // sure that we get only the
- // largest magnitude for each
- // value of |k|
- std::map<unsigned int, double> k_to_max_U_map;
- for (unsigned int f=0; f<n_fourier_modes; ++f)
- if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
- k_to_max_U_map.end())
- ||
- (k_to_max_U_map[k_vectors_magnitude[f]] <
- std::abs (fourier_coefficients[f])))
- k_to_max_U_map[k_vectors_magnitude[f]]
- = std::abs (fourier_coefficients[f]);
-
- // now we have to calculate the
- // various contributions to the
- // formula for mu. we'll only
- // take those fourier
- // coefficients with the
- // largest value for a given
- // |k|
- double sum_1 = 0,
- sum_ln_k = 0,
- sum_ln_k_square = 0,
- sum_ln_U = 0,
- sum_ln_U_ln_k = 0;
- for (unsigned int f=0; f<n_fourier_modes; ++f)
- if (k_to_max_U_map[k_vectors_magnitude[f]] ==
- std::abs (fourier_coefficients[f]))
- {
- sum_1 += 1;
- sum_ln_k += ln_k[f];
- sum_ln_k_square += ln_k[f]*ln_k[f];
- sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
- sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) * ln_k[f];
- }
-
- const double mu
- = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
- *
- (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
-
- smoothness_indicators(index) = mu - 1.*dim/2;
- }
-}
-
-
+ // @sect4{LaplaceProblem::postprocess}
template <int dim>
void LaplaceProblem<dim>::postprocess (const unsigned int cycle)
}
+ // @sect4{LaplaceProblem::create_coarse_grid}
template <>
void LaplaceProblem<2>::create_coarse_grid ()
{
+ // @sect4{LaplaceProblem::run}
template <int dim>
void LaplaceProblem<dim>::run ()
{
}
}
+
+
+template <int dim>
+void
+LaplaceProblem<dim>::
+estimate_smoothness (Vector<float> &smoothness_indicators) const
+{
+ const unsigned int N = max_degree;
+
+ // form all the Fourier vectors
+ // that we want to
+ // consider. exclude k=0 to avoid
+ // problems with |k|^{-mu} and also
+ // logarithms of |k|
+ std::vector<Tensor<1,dim> > k_vectors;
+ std::vector<unsigned int> k_vectors_magnitude;
+ switch (dim)
+ {
+ case 2:
+ {
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ if (!((i==0) && (j==0))
+ &&
+ (i*i + j*j < N*N))
+ {
+ k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
+ deal_II_numbers::PI * j));
+ k_vectors_magnitude.push_back (i*i+j*j);
+ }
+
+ break;
+ }
+
+ case 3:
+ {
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ for (unsigned int k=0; k<N; ++k)
+ if (!((i==0) && (j==0) && (k==0))
+ &&
+ (i*i + j*j + k*k < N*N))
+ {
+ k_vectors.push_back (Point<dim>(deal_II_numbers::PI * i,
+ deal_II_numbers::PI * j,
+ deal_II_numbers::PI * k));
+ k_vectors_magnitude.push_back (i*i+j*j+k*k);
+ }
+
+ break;
+ }
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ const unsigned n_fourier_modes = k_vectors.size();
+ std::vector<double> ln_k (n_fourier_modes);
+ for (unsigned int i=0; i<n_fourier_modes; ++i)
+ ln_k[i] = std::log (k_vectors[i].norm());
+
+
+ // assemble the matrices that do
+ // the Fourier transforms for each
+ // of the finite elements we deal
+ // with. note that these matrices
+ // are complex-valued, so we can't
+ // use FullMatrix
+ QGauss<1> base_quadrature (2);
+ QIterated<dim> quadrature (base_quadrature, N);
+
+ std::vector<Table<2,std::complex<double> > >
+ fourier_transform_matrices (fe_collection.size());
+ for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
+ {
+ fourier_transform_matrices[fe].reinit (n_fourier_modes,
+ fe_collection[fe].dofs_per_cell);
+
+ for (unsigned int k=0; k<n_fourier_modes; ++k)
+ for (unsigned int i=0; i<fe_collection[fe].dofs_per_cell; ++i)
+ {
+ std::complex<double> sum = 0;
+ for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
+ {
+ const Point<dim> x_q = quadrature.point(q);
+ sum += std::exp(std::complex<double>(0,1) *
+ (k_vectors[k] * x_q)) *
+ fe_collection[fe].shape_value(i,x_q) *
+ quadrature.weight(q);
+ }
+ fourier_transform_matrices[fe](k,i)
+ = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
+ }
+ }
+
+ // the next thing is to loop over
+ // all cells and do our work there,
+ // i.e. to locally do the Fourier
+ // transform and estimate the decay
+ // coefficient
+ std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
+ Vector<double> local_dof_values;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (unsigned int index=0; cell!=endc; ++cell, ++index)
+ {
+ local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+ cell->get_dof_values (solution, local_dof_values);
+
+ // first compute the Fourier
+ // transform of the local
+ // solution
+ std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ fourier_coefficients[f] +=
+ fourier_transform_matrices[cell->active_fe_index()](f,i)
+ *
+ local_dof_values(i);
+
+ // enter the Fourier
+ // coefficients into a map with
+ // the k-magnitudes, to make
+ // sure that we get only the
+ // largest magnitude for each
+ // value of |k|
+ std::map<unsigned int, double> k_to_max_U_map;
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
+ k_to_max_U_map.end())
+ ||
+ (k_to_max_U_map[k_vectors_magnitude[f]] <
+ std::abs (fourier_coefficients[f])))
+ k_to_max_U_map[k_vectors_magnitude[f]]
+ = std::abs (fourier_coefficients[f]);
+
+ // now we have to calculate the
+ // various contributions to the
+ // formula for mu. we'll only
+ // take those fourier
+ // coefficients with the
+ // largest value for a given
+ // |k|
+ double sum_1 = 0,
+ sum_ln_k = 0,
+ sum_ln_k_square = 0,
+ sum_ln_U = 0,
+ sum_ln_U_ln_k = 0;
+ for (unsigned int f=0; f<n_fourier_modes; ++f)
+ if (k_to_max_U_map[k_vectors_magnitude[f]] ==
+ std::abs (fourier_coefficients[f]))
+ {
+ sum_1 += 1;
+ sum_ln_k += ln_k[f];
+ sum_ln_k_square += ln_k[f]*ln_k[f];
+ sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
+ sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) * ln_k[f];
+ }
+
+ const double mu
+ = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
+ *
+ (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
+
+ smoothness_indicators(index) = mu - 1.*dim/2;
+ }
+}
+
+
+ // @sect3{The main function}
+
+ // The main function is again verbatim what
+ // we had before: wrap creating and running
+ // an object of the main class into a
+ // <code>try</code> block and catch whatever
+ // exceptions are thrown, thereby producing
+ // meaningful output if anything should go
+ // wrong:
int main ()
{
try
projects / dealii-svn.git / commitdiff