--- /dev/null
+New: The tutorial examples now all use an interpolation formula with a degree depending on the degree of the finite element, which can be changed dynamically.
+<br>
+(Roland Richter, 2019/04/25)
this->triangulation->n_active_cells());
KellyErrorEstimator<dim>::estimate(
this->dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(this->fe->degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
this->solution,
estimated_error_per_cell);
this->triangulation->n_active_cells());
KellyErrorEstimator<dim>::estimate(
this->dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(this->fe->degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
this->solution,
estimated_error_per_cell);
// Initialize a <code>FEValues</code> object with a quadrature formula,
// have abbreviations for the number of quadrature points and shape
// functions...
- QGauss<dim> quadrature(4);
+ QGauss<dim> quadrature(dof_handler.get_fe().degree + 1);
FEValues<dim> fe_values(dof_handler.get_fe(),
quadrature,
update_gradients | update_quadrature_points |
template <int dim>
void MinimalSurfaceProblem<dim>::assemble_system()
{
- const QGauss<dim> quadrature_formula(3);
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
system_matrix = 0;
system_rhs = 0;
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
present_solution,
estimated_error_per_cell);
evaluation_point = present_solution;
evaluation_point.add(alpha, newton_update);
- const QGauss<dim> quadrature_formula(3);
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_gradients | update_quadrature_points |
DEAL:cg::Starting value 0.510691
DEAL:cg::Convergence step 6 value 4.59193e-14
DEAL::Cycle 1
-DEAL:: Number of active cells: 41
-DEAL:: Number of degrees of freedom: 52 (by level: 8, 25, 41)
-DEAL:cg::Starting value 0.455356
-DEAL:cg::Convergence step 8 value 3.09682e-13
+DEAL:: Number of active cells: 44
+DEAL:: Number of degrees of freedom: 55 (by level: 8, 25, 45)
+DEAL:cg::Starting value 0.440678
+DEAL:cg::Convergence step 8 value 1.99419e-13
DEAL::Cycle 2
-DEAL:: Number of active cells: 80
-DEAL:: Number of degrees of freedom: 100 (by level: 8, 25, 61, 52)
-DEAL:cg::Starting value 0.394469
-DEAL:cg::Convergence step 9 value 1.96993e-13
+DEAL:: Number of active cells: 86
+DEAL:: Number of degrees of freedom: 105 (by level: 8, 25, 69, 49)
+DEAL:cg::Starting value 0.371855
+DEAL:cg::Convergence step 9 value 1.13984e-13
DEAL::Cycle 3
-DEAL:: Number of active cells: 161
-DEAL:: Number of degrees of freedom: 190 (by level: 8, 25, 77, 160)
-DEAL:cg::Starting value 0.322156
-DEAL:cg::Convergence step 9 value 2.94418e-13
+DEAL:: Number of active cells: 170
+DEAL:: Number of degrees of freedom: 200 (by level: 8, 25, 77, 174)
+DEAL:cg::Starting value 0.318967
+DEAL:cg::Convergence step 9 value 2.62112e-13
DEAL::Cycle 4
-DEAL:: Number of active cells: 311
-DEAL:: Number of degrees of freedom: 364 (by level: 8, 25, 86, 227, 174)
-DEAL:cg::Starting value 0.279667
-DEAL:cg::Convergence step 10 value 3.45746e-13
+DEAL:: Number of active cells: 332
+DEAL:: Number of degrees of freedom: 388 (by level: 8, 25, 86, 231, 204)
+DEAL:cg::Starting value 0.276534
+DEAL:cg::Convergence step 10 value 1.69562e-13
DEAL::Cycle 5
-DEAL:: Number of active cells: 593
-DEAL:: Number of degrees of freedom: 667 (by level: 8, 25, 89, 231, 490, 96)
-DEAL:cg::Starting value 0.215917
-DEAL:cg::Convergence step 10 value 1.03758e-13
+DEAL:: Number of active cells: 632
+DEAL:: Number of degrees of freedom: 714 (by level: 8, 25, 89, 231, 514, 141)
+DEAL:cg::Starting value 0.215300
+DEAL:cg::Convergence step 10 value 6.47463e-13
DEAL::Cycle 6
-DEAL:: Number of active cells: 1127
-DEAL:: Number of degrees of freedom: 1251 (by level: 8, 25, 89, 274, 760, 417, 178)
-DEAL:cg::Starting value 0.185906
-DEAL:cg::Convergence step 10 value 3.40351e-13
+DEAL:: Number of active cells: 1202
+DEAL:: Number of degrees of freedom: 1332 (by level: 8, 25, 89, 282, 771, 435, 257)
+DEAL:cg::Starting value 0.175848
+DEAL:cg::Convergence step 10 value 1.80664e-13
DEAL::Cycle 7
-DEAL:: Number of active cells: 2144
-DEAL:: Number of degrees of freedom: 2359 (by level: 8, 25, 89, 308, 779, 1262, 817)
-DEAL:cg::Starting value 0.141519
-DEAL:cg::Convergence step 10 value 5.74965e-13
+DEAL:: Number of active cells: 2288
+DEAL:: Number of degrees of freedom: 2511 (by level: 8, 25, 89, 318, 779, 1420, 829, 30)
+DEAL:cg::Starting value 0.136724
+DEAL:cg::Convergence step 11 value 9.73331e-14
</pre>
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
Vector<float> local_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(2),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
localized_solution,
local_error_per_cell,
: triangulation(MPI_COMM_WORLD)
, fe(FE_Q<dim>(1), dim)
, dof_handler(triangulation)
- , quadrature_formula(2)
+ , quadrature_formula(fe.degree + 1)
, present_time(0.0)
, present_timestep(1.0)
, end_time(10.0)
Vector<float> error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(2),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
incremental_displacement,
error_per_cell,
matrix_u.reinit(sparsity_pattern);
matrix_v.reinit(sparsity_pattern);
- MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+ MatrixCreator::create_mass_matrix(dof_handler,
+ QGauss<dim>(fe.degree + 1),
+ mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
laplace_matrix);
// The rest of the function is spent on setting vector sizes to the
VectorTools::project(dof_handler,
constraints,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
InitialValuesU<dim>(),
old_solution_u);
VectorTools::project(dof_handler,
constraints,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
InitialValuesV<dim>(),
old_solution_v);
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
- MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+ MatrixCreator::create_mass_matrix(dof_handler,
+ QGauss<dim>(fe.degree + 1),
+ mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
laplace_matrix);
// The second difference, as mentioned, to step-23 is that we need to
// do something only if that particular face is at the boundary of the
// domain. Like this:
{
- const QGauss<dim - 1> quadrature_formula(3);
+ const QGauss<dim - 1> quadrature_formula(fe.degree + 1);
FEFaceValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_JxW_values);
VectorTools::project(dof_handler,
constraints,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
InitialValuesP<dim>(),
old_solution_p);
old_solution_v = 0;
mass_matrix.reinit(sparsity_pattern);
laplace_matrix.reinit(sparsity_pattern);
- MatrixCreator::create_mass_matrix(dof_handler, QGauss<dim>(3), mass_matrix);
+ MatrixCreator::create_mass_matrix(dof_handler,
+ QGauss<dim>(fe.degree + 1),
+ mass_matrix);
MatrixCreator::create_laplace_matrix(dof_handler,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
laplace_matrix);
solution.reinit(dof_handler.n_dofs());
Vector<double> &nl_term) const
{
nl_term = 0;
- const QGauss<dim> quadrature_formula(3);
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_JxW_values |
const Vector<double> &new_data,
SparseMatrix<double> &nl_matrix) const
{
- QGauss<dim> quadrature_formula(3);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_JxW_values |
constraints.close();
VectorTools::project(dof_handler,
constraints,
- QGauss<dim>(3),
+ QGauss<dim>(fe.degree + 1),
InitialValues<dim>(1, time),
solution);
}
// used. Since our bilinear form involves boundary integrals on
// $\Gamma_2$, we also need a quadrature rule for surface integration on
// the faces, which are $dim-1$ dimensional:
- QGauss<dim> quadrature_formula(2);
- QGauss<dim - 1> face_quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
+ QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
const unsigned int n_q_points = quadrature_formula.size(),
n_face_q_points = face_quadrature_formula.size(),
@code
std::cout << "Mean value: "
<< VectorTools::compute_mean_value (dof_handler,
- QGauss<2>(3),
+ QGauss<2>(fe.degree + 1),
solution,
0)
<< std::endl;
// 2D. This quadrature formula integrates polynomials of degrees up to three
// exactly (in 1D). It is easy to check that this is sufficient for the
// present problem:
- QGauss<2> quadrature_formula(2);
+ QGauss<2> quadrature_formula(fe.degree + 1);
// And we initialize the object which we have briefly talked about above. It
// needs to be told which finite element we want to use, and the quadrature
// points and their weights (jointly described by a Quadrature object). As
: mapping()
, fe(FE_Q<dim>(1), EulerEquations<dim>::n_components)
, dof_handler(triangulation)
- , quadrature(2)
- , face_quadrature(2)
+ , quadrature(fe.degree + 1)
+ , face_quadrature(fe.degree + 1)
, verbose_cout(std::cout, false)
{
ParameterHandler prm;
template <int dim>
void EigenvalueProblem<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
template <int dim>
void Step4<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
// We wanted to have a non-constant right hand side, so we use an object of
// the class declared above to generate the necessary data. Since this right
{
TimerOutput::Scope t(computing_timer, "assembly");
- const QGauss<dim> quadrature_formula(3);
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
locally_relevant_solution,
estimated_error_per_cell);
template <int dim>
void Step5<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
std::map<types::boundary_id, const Function<dim> *>
neumann_boundary;
KellyErrorEstimator<dim>::estimate(dof_handler_local,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
neumann_boundary,
solution_local,
estimated_error_per_cell,
FEValuesExtractors::Vector velocities(0);
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
locally_relevant_solution,
estimated_error_per_cell,
template <int dim>
void Step6<dim>::assemble_system()
{
- const QGauss<dim> quadrature_formula(3);
+ const QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
Number of degrees of freedom: 15912
cycle cells dofs L2 H1 Linfty
- 0 64 81 1.576e-01 1.418e+00 2.707e-01
- 1 124 157 4.285e-02 1.285e+00 1.469e-01
- 2 280 341 1.593e-02 7.909e-01 8.034e-02
- 3 577 690 9.359e-03 5.096e-01 2.784e-02
- 4 1099 1264 2.865e-03 3.038e-01 9.822e-03
- 5 2191 2452 1.480e-03 2.106e-01 5.679e-03
- 6 4165 4510 6.907e-04 1.462e-01 2.338e-03
- 7 7915 8440 4.743e-04 1.055e-01 1.442e-03
- 8 15196 15912 1.920e-04 7.468e-02 7.259e-04
+ 0 64 81 1.840e+00 2.858e+00 1.835e+00
+ 1 124 157 5.190e-02 1.200e+00 1.344e-01
+ 2 280 341 1.439e-02 7.892e-01 7.554e-02
+ 3 577 690 8.627e-03 5.061e-01 2.805e-02
+ 4 1099 1264 3.217e-03 3.030e-01 1.073e-02
+ 5 2191 2452 1.445e-03 2.097e-01 5.073e-03
+ 6 4165 4510 8.387e-04 1.460e-01 2.013e-03
+ 7 7915 8440 7.051e-04 1.053e-01 1.804e-03
+ 8 15196 15912 2.774e-04 7.463e-02 6.911e-04
Solving with Q1 elements, global refinement
===========================================
Number of degrees of freedom: 16641
cycle cells dofs L2 H1 Linfty
- 0 64 81 1.576e-01 1.418e+00 2.707e-01
- 1 256 289 4.280e-02 1.285e+00 1.444e-01
- 2 1024 1089 1.352e-02 7.556e-01 7.772e-02
- 3 4096 4225 3.423e-03 3.822e-01 2.332e-02
- 4 16384 16641 8.586e-04 1.917e-01 6.097e-03
-
-n cells H1 L2
-0 64 1.418e+00 - - 1.576e-01 - -
-1 256 1.285e+00 1.10 0.14 4.280e-02 3.68 1.88
-2 1024 7.556e-01 1.70 0.77 1.352e-02 3.17 1.66
-3 4096 3.822e-01 1.98 0.98 3.423e-03 3.95 1.98
-4 16384 1.917e-01 1.99 1.00 8.586e-04 3.99 2.00
+ 0 64 81 1.840e+00 2.858e+00 1.835e+00
+ 1 256 289 3.570e-02 1.199e+00 1.307e-01
+ 2 1024 1089 1.192e-02 7.565e-01 7.168e-02
+ 3 4096 4225 3.047e-03 3.823e-01 2.128e-02
+ 4 16384 16641 7.660e-04 1.917e-01 5.554e-03
+
+n cells H1 L2
+0 64 2.858e+00 - - 1.840e+00 - -
+1 256 1.199e+00 2.38 1.25 3.570e-02 51.54 5.69
+2 1024 7.565e-01 1.58 0.66 1.192e-02 2.99 1.58
+3 4096 3.823e-01 1.98 0.98 3.047e-03 3.91 1.97
+4 16384 1.917e-01 1.99 1.00 7.660e-04 3.98 1.99
Solving with Q2 elements, global refinement
===========================================
5 2059 9223 7.738e-05 1.974e-02 7.270e-04
6 3913 17887 2.925e-05 8.772e-03 1.463e-04
7 7441 33807 1.024e-05 4.121e-03 8.567e-05
- 8 14212 64731 3.761e-06 2.108e-03 2.167e-05
+ 8 14212 64731 3.761e-06 2.108e-03 2.167e-05
@endcode
Go ahead and run the program with higher order elements (Q3, Q4, ...). You
will notice that assertions in several parts of the code will trigger (for
-example in the generation of the filename for the data output). After fixing
-these you will not see the correct convergence orders that the theory
-predicts. This is because the orders for the quadrature formulas are
-hard-coded in this program and this order is not enough for higher order
-discretizations. What is a good way to pick the orders dynamically?
+example in the generation of the filename for the data output). You might have to address these,
+but it should not be very hard to get the program to work!
<h4> Convergence Comparison </h4>
template <int dim>
void HelmholtzProblem<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(3);
- QGauss<dim - 1> face_quadrature_formula(3);
+ QGauss<dim> quadrature_formula(fe->degree + 1);
+ QGauss<dim - 1> face_quadrature_formula(fe->degree + 1);
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_face_q_points = face_quadrature_formula.size();
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(fe->degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
solution,
Solution<dim>(),
difference_per_cell,
- QGauss<dim>(3),
+ QGauss<dim>(fe->degree + 1),
VectorTools::L2_norm);
const double L2_error =
VectorTools::compute_global_error(triangulation,
solution,
Solution<dim>(),
difference_per_cell,
- QGauss<dim>(3),
+ QGauss<dim>(fe->degree + 1),
VectorTools::H1_seminorm);
const double H1_error =
VectorTools::compute_global_error(triangulation,
// points. Since this depends quite sensitively on the quadrature rule
// being used, and since we would like to avoid false results due to
// super-convergence effects at some points, we use a special quadrature
- // rule that is obtained by iterating the trapezoidal rule five times in
- // each space direction. Note that the constructor of the QIterated class
+ // rule that is obtained by iterating the trapezoidal rule by the degree of
+ // of the finite element times two plus one in each space direction.
+ // Note that the constructor of the QIterated class
// takes a one-dimensional quadrature rule and a number that tells it how
// often it shall use this rule in each space direction.
//
// from the L infinity errors on each cell with a call to
// VectorTools::compute_global_error.
const QTrapez<1> q_trapez;
- const QIterated<dim> q_iterated(q_trapez, 5);
+ const QIterated<dim> q_iterated(q_trapez, fe->degree * 2 + 1);
VectorTools::integrate_difference(dof_handler,
solution,
Solution<dim>(),
template <int dim>
void ElasticProblem<dim>::assemble_system()
{
- QGauss<dim> quadrature_formula(2);
+ QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(2),
+ QGauss<dim - 1>(fe.degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
projects / dealii.git / commitdiff