From: Luca Heltai Date: Wed, 28 Jul 2010 13:18:59 +0000 (+0000) Subject: Added support for higher order mapping in step-34 X-Git-Tag: v8.0.0~5760 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=81317d71c1f5fb283665f68bd713d321836f4364;p=dealii.git Added support for higher order mapping in step-34 git-svn-id: https://svn.dealii.org/trunk@21581 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-34/doc/results.dox b/deal.II/examples/step-34/doc/results.dox index 01e7b38fe6..3c40a63fc7 100644 --- a/deal.II/examples/step-34/doc/results.dox +++ b/deal.II/examples/step-34/doc/results.dox @@ -154,12 +154,42 @@ equal size generates a regular polygon with N faces, whose angles are exactly $\pi-\frac {2\pi}{N}$, therefore the error we commit should be exactly $\frac 12 - (\frac 12 -\frac 1N) = \frac 1N$. In fact this is a very good indicator that we are performing the singular integrals in -an appropiate manner. +an appropriate manner. The error in the approximation of the potential $\phi$ is largely due to approximation of the domain. A much better approximation could be obtained by using higher order mappings. +If we modify the main() function, setting fe_degree and mapping_degree +to two, and raise the order of the quadrature formulas in +the parameter file, we obtain the following convergence table for the +two dimensional simulation + +@verbatim +cycle cells dofs L2(phi) Linfty(alpha) + 0 20 40 5.404e-05 - 2.306e-04 - + 1 40 80 3.578e-06 3.92 1.738e-05 3.73 + 2 80 160 2.479e-07 3.85 1.253e-05 0.47 + 3 160 320 1.856e-08 3.74 7.670e-06 0.71 +@endverbatim + +and + +@verbatim +cycle cells dofs L2(phi) Linfty(alpha) + 0 24 98 9.187e-03 - 8.956e-03 - + 1 96 386 3.991e-04 4.52 1.182e-03 2.92 + 2 384 1538 2.113e-05 4.24 1.499e-04 2.98 + 3 1536 6146 1.247e-06 4.08 1.896e-05 2.98 +@endverbatim + +for the three dimensional case. As we can see, convergence results are +much better with higher order mapping, mainly due to a better +resolution of the curved geometry. Notice that, given the same number +of degrees of freedom, for example in step 3 of the Q1 case and step 2 +of Q2 case in the three dimensional simulation, the error is roughly +three orders of magnitude lower. + The result of running these computations is a bunch of output files that we can pass to our visualization program of choice. The output files are of two kind: the potential on the boundary diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index b85e66d9b2..8d993b2b70 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -48,7 +48,7 @@ #include #include -#include +#include #include #include @@ -141,7 +141,8 @@ template class BEMProblem { public: - BEMProblem(); + BEMProblem(const unsigned int fe_degree = 1, + const unsigned int mapping_degree = 1); void run(); @@ -296,6 +297,19 @@ class BEMProblem void output_results(const unsigned int cycle); + // To allow for dimension + // independent programming, we + // specialize this single + // function to extract the + // singular quadrature formula + // needed to integrate the + // singular kernels in the + // interior of the cells. + const Quadrature & get_singular_quadrature( + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int index) const; + + // The usual deal.II classes can // be used for boundary element // methods by specifying the @@ -317,10 +331,21 @@ class BEMProblem // usual finite element classes // that we saw in all previous // examples. + // + // The class is constructed in a + // way to allow for arbitrary + // order of approximation of both + // the domain (through high order + // mapping) and the finite + // element space. The order of + // the finite element space and + // of the mapping can be selected + // in the constructor of the class. Triangulation tria; FE_Q fe; DoFHandler dh; + MappingQ mapping; // In BEM methods, the matrix // that is generated is @@ -348,6 +373,7 @@ class BEMProblem // from a point $\mathbf x$) at // the support points of our // shape functions. + Vector phi; Vector alpha; @@ -355,6 +381,7 @@ class BEMProblem // to output errors in the exact // solution and in the computed // alphas. + ConvergenceTable convergence_table; // The following variables are @@ -390,28 +417,17 @@ class BEMProblem // to define the order of the // singular quadrature rule. // - // Notice that the pointer given - // below for the quadrature rule - // is only used for non singular - // integrals. Whenever the - // integral is singular, then - // only the degree of the - // quadrature pointer is used, - // and the integration is a - // special one (see the - // assemble_system() function - // below for further details). - // // We also define a couple of // parameters which are used in // case we wanted to extend the // solution to the entire domain. + Functions::ParsedFunction wind; Functions::ParsedFunction exact_solution; - std_cxx1x::shared_ptr > quadrature; unsigned int singular_quadrature_order; - + std_cxx1x::shared_ptr > quadrature; + SolverControl solver_control; unsigned int n_cycles; @@ -451,11 +467,13 @@ class BEMProblem // is static, and has no knowledge of // the number of components. template -BEMProblem::BEMProblem() +BEMProblem::BEMProblem(const unsigned int fe_degree, + const unsigned int mapping_degree) : - fe(1), + fe(fe_degree), dh(tria), - wind(dim) + wind(dim), + mapping(mapping_degree, true) {} @@ -716,6 +734,7 @@ void BEMProblem::read_domain() GridIn gi; gi.attach_triangulation (tria); gi.read_ucd (in); + tria.set_boundary(1, boundary); } @@ -750,38 +769,11 @@ void BEMProblem::refine_and_resize() // of this program, assembling the // matrix that corresponds to the // boundary integral equation. - // - // At the beginning, we create the - // singular quadratures for the three - // dimensional problem (note that a - // 3d boundary integral problem - // requires a 2d quadrature - // formula!), since in this case they - // only depend on the reference - // element. This quadrature is a - // standard Gauss quadrature formula - // reparametrized in such a way that - // allows one to integrate - // singularities of the kind $1/R$ - // centered at one of the - // vertices. Here we define a vector - // of four such quadratures (one per - // vertex of the two dimensional - // cells for a surface in 3d) that - // will be used later on; note, - // however, that these objects will - // only be used in the three - // dimensional case. template void BEMProblem::assemble_system() { - std::vector > sing_quadratures_3d; - for(unsigned int i=0; i<4; ++i) - sing_quadratures_3d.push_back - (QGaussOneOverR<2>(singular_quadrature_order, i, true)); - - // Next, we initialize an FEValues + // First we initialize an FEValues // object with the quadrature // formula for the integration of // the kernel in non singular @@ -791,7 +783,7 @@ void BEMProblem::assemble_system() // precise, since the functions we // are integrating are not // polynomial functions. - FEValues fe_v(fe, *quadrature, + FEValues fe_v(mapping, fe, *quadrature, update_values | update_cell_normal_vectors | update_quadrature_points | @@ -826,28 +818,15 @@ void BEMProblem::assemble_system() // collocation points, which are // the support points of the $i$th // basis function, while $j$ runs - // on inner integration points. We - // perform the following check to - // ensure that we are not trying to - // use this code for high order - // elements. It will only work with - // Q1 elements, that is, for - // fe.dofs_per_cell == - // GeometryInfo::vertices_per_cell. - AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, - ExcMessage("The code in this function can only be used for " - "the usual Q1 elements.")); - - // Now that we have checked that - // the number of vertices is equal - // to the number of degrees of - // freedom, we construct a vector + // on inner integration points. + + // We construct a vector // of support points which will be // used in the local integrations: std::vector > support_points(dh.n_dofs()); - DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, - dh, support_points); + DoFTools::map_dofs_to_support_points( mapping, dh, support_points); + // After doing so, we can start the // integration loop over all cells, // where we first initialize the @@ -870,7 +849,6 @@ void BEMProblem::assemble_system() const std::vector > &normals = fe_v.get_cell_normal_vectors(); wind.vector_value_list(q_points, cell_wind); - // We then form the integral over // the current cell for all // degrees of freedom (note that @@ -960,213 +938,20 @@ void BEMProblem::assemble_system() // functions against a // singular weight on the // reference cell. - // Notice that singular - // integration requires a - // careful selection of - // the quadrature - // rules. In particular - // the deal.II library - // provides quadrature - // rules which are - // taylored for - // logarithmic - // singularities - // (QGaussLog, - // QGaussLogR), as well - // as for 1/R - // singularities - // (QGaussOneOverR). - // - // Singular integration - // is typically obtained - // by constructing - // weighted quadrature - // formulas with singular - // weights, so that it is - // possible to write - // - // \f[ - // \int_K f(x) s(x) dx = \sum_{i=1}^N w_i f(q_i) - // \f] - // - // where $s(x)$ is a given - // singularity, and the weights - // and quadrature points - // $w_i,q_i$ are carefully - // selected to make the formula - // above an equality for a - // certain class of functions - // $f(x)$. - // - // In all the finite - // element examples we - // have seen so far, the - // weight of the - // quadrature itself - // (namely, the function - // $s(x)$), was always - // constantly equal to 1. - // For singular - // integration, we have - // two choices: we can - // use the definition - // above, factoring out - // the singularity from - // the integrand (i.e., - // integrating $f(x)$ - // with the special - // quadrature rule), or - // we can ask the - // quadrature rule to - // "normalize" the - // weights $w_i$ with - // $s(q_i)$: - // - // \f[ - // \int_K f(x) s(x) dx = - // \int_K g(x) dx = \sum_{i=1}^N \frac{w_i}{s(q_i)} g(q_i) - // \f] - // - // We use this second - // option, through the @p - // factor_out_singularity - // parameter of both - // QGaussLogR and - // QGaussOneOverR. - // - // These integrals are - // somewhat delicate, - // especially in two - // dimensions, due to the - // transformation from - // the real to the - // reference cell, where - // the variable of - // integration is scaled - // with the determinant - // of the transformation. - // - // In two dimensions this - // process does not - // result only in a - // factor appearing as a - // constant factor on the - // entire integral, but - // also on an additional - // integral alltogether - // that needs to be - // evaluated: - // - // \f[ - // \int_0^1 f(x)\ln(x/\alpha) dx = - // \int_0^1 f(x)\ln(x) dx - \int_0^1 f(x) \ln(\alpha) dx. - // \f] - // - // This process is taken care of by - // the constructor of the QGaussLogR - // class, which adds additional - // quadrature points and weights to - // take into consideration also the - // second part of the integral. - // - // A similar reasoning - // should be done in the - // three dimensional - // case, since the - // singular quadrature is - // taylored on the - // inverse of the radius - // $r$ in the reference - // cell, while our - // singular function - // lives in real space, - // however in the three - // dimensional case - // everything is simpler - // because the - // singularity scales - // linearly with the - // determinant of the - // transformation. This - // allows us to build the - // singular two - // dimensional quadrature - // rules once and for all - // outside the loop over - // all cells, using only - // a pointer where needed. - // - // Notice that in one - // dimensional - // integration this is - // not possible, since we - // need to know the - // scaling parameter for - // the quadrature, which - // is not known a - // priori. Here, the - // quadrature rule itself - // depends also on the - // size of the current - // cell. For this reason, - // it is necessary to - // create a new - // quadrature for each - // singular - // integration. Since we - // create it using the - // new operator of C++, - // we also need to - // destroy it using the - // dual of new: - // delete. This is done - // at the end, and only - // if dim == 2. // - // Putting all this into a - // dimension independent - // framework requires a little - // trick. The problem is that, - // depending on dimension, we'd - // like to either assign a - // QGaussLogR<1> or a - // QGaussOneOverR<2> to a - // Quadrature. C++ - // doesn't allow this right - // away, and neither is a - // static_cast - // possible. However, we can - // attempt a dynamic_cast: the - // implementation will then - // look up at run time whether - // the conversion is possible - // (which we know it - // is) and if that isn't the - // case simply return a null - // pointer. To be sure we can - // then add a safety check at - // the end: + // The correct quadrature + // formula is selected by + // the + // get_singular_quadrature + // function, which is + // explained in detail below. Assert(singular_index != numbers::invalid_unsigned_int, ExcInternalError()); + + const Quadrature & singular_quadrature = + get_singular_quadrature(cell, singular_index); - const Quadrature * - singular_quadrature - = (dim == 2 - ? - dynamic_cast*>( - new QGaussLogR<1>(singular_quadrature_order, - Point<1>((double)singular_index), - 1./cell->measure(), true)) - : - (dim == 3 - ? - dynamic_cast*>( - &sing_quadratures_3d[singular_index]) - : - 0)); - Assert(singular_quadrature, ExcInternalError()); - - FEValues fe_v_singular (fe, *singular_quadrature, + FEValues fe_v_singular (mapping, fe, singular_quadrature, update_jacobians | update_values | update_cell_normal_vectors | @@ -1174,7 +959,7 @@ void BEMProblem::assemble_system() fe_v_singular.reinit(cell); - std::vector > singular_cell_wind( (*singular_quadrature).size(), + std::vector > singular_cell_wind( singular_quadrature.size(), Vector(dim) ); const std::vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); @@ -1182,7 +967,7 @@ void BEMProblem::assemble_system() wind.vector_value_list(singular_q_points, singular_cell_wind); - for(unsigned int q=0; qsize(); ++q) + for(unsigned int q=0; q R = singular_q_points[q] - support_points[i]; double normal_wind = 0; @@ -1201,8 +986,6 @@ void BEMProblem::assemble_system() fe_v_singular.JxW(q) ); } } - if(dim==2) - delete singular_quadrature; } // Finally, we need to add @@ -1274,10 +1057,10 @@ template void BEMProblem::compute_errors(const unsigned int cycle) { Vector difference_per_cell (tria.n_active_cells()); - VectorTools::integrate_difference (dh, phi, + VectorTools::integrate_difference (mapping, dh, phi, exact_solution, difference_per_cell, - QGauss<(dim-1)>(3), + QGauss<(dim-1)>(2*fe.degree+1), VectorTools::L2_norm); const double L2_error = difference_per_cell.l2_norm(); @@ -1316,6 +1099,161 @@ void BEMProblem::compute_errors(const unsigned int cycle) } + // Singular integration requires a + // careful selection of the + // quadrature rules. In particular + // the deal.II library provides + // quadrature rules which are + // taylored for logarithmic + // singularities (QGaussLog, + // QGaussLogR), as well as for 1/R + // singularities (QGaussOneOverR). + // + // Singular integration is typically + // obtained by constructing weighted + // quadrature formulas with singular + // weights, so that it is possible to + // write + // + // \f[ + // \int_K f(x) s(x) dx = \sum_{i=1}^N w_i f(q_i) + // \f] + // + // where $s(x)$ is a given + // singularity, and the weights and + // quadrature points $w_i,q_i$ are + // carefully selected to make the + // formula above an equality for a + // certain class of functions $f(x)$. + // + // In all the finite element examples + // we have seen so far, the weight of + // the quadrature itself (namely, the + // function $s(x)$), was always + // constantly equal to 1. For + // singular integration, we have two + // choices: we can use the definition + // above, factoring out the + // singularity from the integrand + // (i.e., integrating $f(x)$ with the + // special quadrature rule), or we + // can ask the quadrature rule to + // "normalize" the weights $w_i$ with + // $s(q_i)$: + // + // \f[ + // \int_K f(x) s(x) dx = + // \int_K g(x) dx = \sum_{i=1}^N \frac{w_i}{s(q_i)} g(q_i) + // \f] + // + // We use this second option, through + // the @p factor_out_singularity + // parameter of both QGaussLogR and + // QGaussOneOverR. + // + // These integrals are somewhat + // delicate, especially in two + // dimensions, due to the + // transformation from the real to + // the reference cell, where the + // variable of integration is scaled + // with the determinant of the + // transformation. + // + // In two dimensions this process + // does not result only in a factor + // appearing as a constant factor on + // the entire integral, but also on + // an additional integral alltogether + // that needs to be evaluated: + // + // \f[ + // \int_0^1 f(x)\ln(x/\alpha) dx = + // \int_0^1 f(x)\ln(x) dx - \int_0^1 f(x) \ln(\alpha) dx. + // \f] + // + // This process is taken care of by + // the constructor of the QGaussLogR + // class, which adds additional + // quadrature points and weights to + // take into consideration also the + // second part of the integral. + // + // A similar reasoning should be done + // in the three dimensional case, + // since the singular quadrature is + // taylored on the inverse of the + // radius $r$ in the reference cell, + // while our singular function lives + // in real space, however in the + // three dimensional case everything + // is simpler because the singularity + // scales linearly with the + // determinant of the + // transformation. This allows us to + // build the singular two dimensional + // quadrature rules only once and, + // reuse them over all cells. + // + // In the one dimensional singular + // integration this is not possible, + // since we need to know the scaling + // parameter for the quadrature, + // which is not known a priori. Here, + // the quadrature rule itself depends + // also on the size of the current + // cell. For this reason, it is + // necessary to create a new + // quadrature for each singular + // integration. + // + // The different quadrature rules are + // built inside the + // get_singular_quadrature, which is + // specialized for dim=2 and dim=3, + // and they are retrieved inside the + // assemble_system function. The + // index given as an argument is the + // index of the unit support point + // where the singularity is located. + +template<> +const Quadrature<2> & BEMProblem<3>::get_singular_quadrature( + const DoFHandler<2,3>::active_cell_iterator &, + const unsigned int index) const +{ + Assert(index < fe.dofs_per_cell, + ExcIndexRange(0, fe.dofs_per_cell, index)); + + static std::vector > quadratures; + if(quadratures.size() == 0) + for(unsigned int i=0; i(singular_quadrature_order, + fe.get_unit_support_points()[i], + true)); + return quadratures[index]; +} + + +template<> +const Quadrature<1> & BEMProblem<2>::get_singular_quadrature( + const DoFHandler<1,2>::active_cell_iterator &cell, + const unsigned int index) const +{ + Assert(index < fe.dofs_per_cell, + ExcIndexRange(0, fe.dofs_per_cell, index)); + + static Quadrature<1> * q_pointer = NULL; + if(q_pointer) delete q_pointer; + + q_pointer = new QGaussLogR<1>(singular_quadrature_order, + fe.get_unit_support_points()[index], + 1./cell->measure(), true); + return (*q_pointer); +} + + + // @sect4{BEMProblem::compute_exterior_solution} // We'd like to also know something @@ -1365,7 +1303,7 @@ void BEMProblem::compute_exterior_solution() endc = dh.end(); - FEValues fe_v(fe, *quadrature, + FEValues fe_v(mapping, fe, *quadrature, update_values | update_cell_normal_vectors | update_quadrature_points | @@ -1380,8 +1318,8 @@ void BEMProblem::compute_exterior_solution() std::vector > local_wind(n_q_points, Vector(dim) ); std::vector > external_support_points(external_dh.n_dofs()); - DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, - external_dh, external_support_points); + DoFTools::map_dofs_to_support_points(StaticMappingQ1::mapping, + external_dh, external_support_points); for(cell = dh.begin_active(); cell != endc; ++cell) { @@ -1446,7 +1384,9 @@ void BEMProblem::output_results(const unsigned int cycle) dataout.attach_dof_handler(dh); dataout.add_data_vector(phi, "phi"); dataout.add_data_vector(alpha, "alpha"); - dataout.build_patches(); + dataout.build_patches(mapping, + mapping.get_degree(), + DataOut >::curved_inner_cells); std::string filename = ( Utilities::int_to_string(dim) + "d_boundary_solution_" + @@ -1518,11 +1458,14 @@ int main () { try { + unsigned int degree = 1; + unsigned int mapping_degree = 1; + deallog.depth_console (3); - BEMProblem<2> laplace_problem_2d; + BEMProblem<2> laplace_problem_2d(degree, mapping_degree); laplace_problem_2d.run(); - BEMProblem<3> laplace_problem_3d; + BEMProblem<3> laplace_problem_3d(degree, mapping_degree); laplace_problem_3d.run(); } catch (std::exception &exc)