From a041e7cb2ddadb30faabccfc1b60132384e3487b Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 27 Jul 2010 17:32:46 +0000 Subject: [PATCH] Added support for MappingQ of higher order also in codim one case. git-svn-id: https://svn.dealii.org/trunk@21577 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/mapping_q.h | 16 +-- deal.II/deal.II/source/fe/mapping_q.cc | 116 ++++++++++++++---- .../deal.II/source/fe/mapping_q_eulerian.cc | 12 ++ .../deal.II/source/grid/tria_boundary_lib.cc | 36 +++++- deal.II/deal.II/source/numerics/data_out.cc | 4 +- deal.II/doc/news/changes.h | 6 + 6 files changed, 155 insertions(+), 35 deletions(-) diff --git a/deal.II/deal.II/include/fe/mapping_q.h b/deal.II/deal.II/include/fe/mapping_q.h index 4dfe95605e..737add3c16 100644 --- a/deal.II/deal.II/include/fe/mapping_q.h +++ b/deal.II/deal.II/include/fe/mapping_q.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors +// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -226,7 +226,7 @@ class MappingQ : public MappingQ1 typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors) const ; /** * Implementation of the interface in @@ -241,7 +241,7 @@ class MappingQ : public MappingQ1 typename std::vector > &quadrature_points, std::vector &JxW_values, typename std::vector > &exterior_form, - typename std::vector > &normal_vectors) const ; + typename std::vector > &normal_vectors) const ; /** * For dim=2,3. Append the @@ -268,7 +268,7 @@ class MappingQ : public MappingQ1 */ virtual void add_line_support_points (const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + std::vector > &a) const; /** * For dim=3. Append the @@ -296,7 +296,7 @@ class MappingQ : public MappingQ1 */ virtual void add_quad_support_points(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + std::vector > &a) const; private: @@ -396,7 +396,7 @@ class MappingQ : public MappingQ1 * `mapping' report. */ void apply_laplace_vector(const Table<2,double> &lvs, - std::vector > &a) const; + std::vector > &a) const; /** * Computes the support points of @@ -404,7 +404,7 @@ class MappingQ : public MappingQ1 */ virtual void compute_mapping_support_points( const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + std::vector > &a) const; /** * Computes all support points of @@ -422,7 +422,7 @@ class MappingQ : public MappingQ1 */ void compute_support_points_laplace( const typename Triangulation::cell_iterator &cell, - std::vector > &a) const; + std::vector > &a) const; /** * Needed by the diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index ce211ca58b..82bc61a809 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -59,7 +59,7 @@ MappingQ::InternalData::memory_consumption () const #if deal_II_dimension == 1 // in 1d, it is irrelevant which polynomial degree to use, since all -// cells are scaled linearly +// cells are scaled linearly. Unless codimension is equal to two template<> MappingQ<1>::MappingQ (const unsigned int, const bool /*use_mapping_q_on_all_cells*/) @@ -88,7 +88,6 @@ MappingQ<1>::MappingQ (const MappingQ<1> &m): feq(degree) {} - template<> MappingQ<1>::~MappingQ () {} @@ -119,7 +118,8 @@ MappingQ::MappingQ (const unsigned int p, : degree(p), n_inner(Utilities::fixed_power(degree-1)), - n_outer((dim==2) ? 4+4*(degree-1) + n_outer((dim==1) ? 2 : + (dim==2) ? 4+4*(degree-1) :8+12*(degree-1)+6*(degree-1)*(degree-1)), tensor_pols(0), n_shape_functions(Utilities::fixed_power(degree+1)), @@ -127,7 +127,8 @@ MappingQ::MappingQ (const unsigned int p, lexicographic_to_hierarchic_numbering ( FiniteElementData (get_dpo_vector(degree), 1, degree))), - use_mapping_q_on_all_cells (use_mapping_q_on_all_cells), + use_mapping_q_on_all_cells (use_mapping_q_on_all_cells + || (dim != spacedim)), feq(degree) { // Construct the tensor product @@ -373,7 +374,7 @@ MappingQ::fill_fe_face_values ( std::vector > &quadrature_points, std::vector &JxW_values, std::vector > &exterior_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors) const { // convert data object to internal // data for this class. fails with @@ -435,7 +436,7 @@ MappingQ::fill_fe_subface_values (const typename Triangulation > &quadrature_points, std::vector &JxW_values, std::vector > &exterior_forms, - std::vector > &normal_vectors) const + std::vector > &normal_vectors) const { // convert data object to internal // data for this class. fails with @@ -842,7 +843,7 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const template void MappingQ::apply_laplace_vector(const Table<2,double> &lvs, - std::vector > &a) const + std::vector > &a) const { // check whether the data we need // is really available. if you fail @@ -873,7 +874,7 @@ MappingQ::apply_laplace_vector(const Table<2,double> &lvs, for (unsigned int unit_point=0; unit_point p; + Point p; for (unsigned int k=0; k void MappingQ::compute_mapping_support_points( const typename Triangulation::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { // if this is a cell for which we // want to compute the full @@ -912,7 +913,7 @@ MappingQ::compute_mapping_support_points( template void MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { // in any case, we need the // vertices first @@ -923,6 +924,10 @@ MappingQ::compute_support_points_laplace(const typename Triangulat if (degree>1) switch (dim) { + case 1: + Assert(spacedim == 2, ExcInternalError()); + add_line_support_points(cell, a); + break; case 2: // in 2d, add the // points on the four @@ -930,7 +935,10 @@ MappingQ::compute_support_points_laplace(const typename Triangulat // the exterior (outer) // points add_line_support_points (cell, a); - apply_laplace_vector (laplace_on_quad_vector,a); + if(dim != spacedim) + add_quad_support_points(cell, a); + else + apply_laplace_vector (laplace_on_quad_vector,a); break; case 3: @@ -966,15 +974,47 @@ MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, Assert (dim > 1, ExcImpossibleInDim(dim)); } + +template <> +void +MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator &cell, + std::vector > &a) const +{ + const unsigned int dim=1, spacedim=2; + // Ask for the mid point, if that's + // the only thing we need. + if (degree==2) + { + const Boundary * const boundary + = &(cell->get_triangulation().get_boundary(cell->material_id())); + a.push_back(boundary->get_new_point_on_line(cell)); + } + else + // otherwise call the more + // complicated functions and ask + // for inner points from the + // boundary description + { + std::vector > line_points (degree-1); + + const Boundary * const boundary + = &(cell->get_triangulation().get_boundary(cell->material_id())); + + boundary->get_intermediate_points_on_line (cell, line_points); + // Append all points + a.insert (a.end(), line_points.begin(), line_points.end()); + } +} + #endif template void MappingQ::add_line_support_points (const typename Triangulation::cell_iterator &cell, - std::vector > &a) const + std::vector > &a) const { - static const StraightBoundary straight_boundary; + static const StraightBoundary straight_boundary; // if we only need the midpoint, // then ask for it. if (degree==2) @@ -982,9 +1022,11 @@ MappingQ::add_line_support_points (const typename Triangulation::lines_per_cell; ++line_no) { const typename Triangulation::line_iterator line = cell->line(line_no); - const Boundary * const boundary - = (line->at_boundary() ? - &line->get_triangulation().get_boundary(line->boundary_indicator()) : + const Boundary * const boundary + = (line->at_boundary()? + &line->get_triangulation().get_boundary(line->boundary_indicator()): + (dim != spacedim) ? + &line->get_triangulation().get_boundary(cell->material_id()): &straight_boundary); a.push_back(boundary->get_new_point_on_line(line)); @@ -996,7 +1038,7 @@ MappingQ::add_line_support_points (const typename Triangulation > line_points (degree-1); + std::vector > line_points (degree-1); // loop over each of the lines, // and if it is at the @@ -1008,9 +1050,11 @@ MappingQ::add_line_support_points (const typename Triangulation::line_iterator line = cell->line(line_no); - const Boundary * const boundary - = (line->at_boundary() ? + const Boundary * const boundary + = (line->at_boundary()? &line->get_triangulation().get_boundary(line->boundary_indicator()) : + (dim != spacedim) ? + &line->get_triangulation().get_boundary(cell->material_id()) : &straight_boundary); boundary->get_intermediate_points_on_line (line, line_points); @@ -1208,12 +1252,30 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, #endif +#if deal_II_dimension == 2 + +template<> +void +MappingQ<2,3>:: +add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, + std::vector > &a) const +{ + std::vector > quad_points ((degree-1)*(degree-1)); + + cell->get_triangulation().get_boundary(cell->material_id()) + .get_intermediate_points_on_quad (cell, quad_points); + for (unsigned int i=0; i void MappingQ:: add_quad_support_points(const typename Triangulation::cell_iterator &, - std::vector > &) const + std::vector > &) const { Assert (dim > 2, ExcImpossibleInDim(dim)); } @@ -1312,8 +1374,8 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it mdata (dynamic_cast ( get_data(update_transformation_values, point_quadrature))); - mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells - || cell->has_boundary_lines()); + mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || + cell->has_boundary_lines()); typename MappingQ1::InternalData *p_data = (mdata->use_mapping_q1_on_current_cell ? @@ -1340,7 +1402,9 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // then a Newton iteration based on // the full MappingQ if we need // this - if (cell->has_boundary_lines() || use_mapping_q_on_all_cells) + if (cell->has_boundary_lines() || + use_mapping_q_on_all_cells || + (dim!=spacedim) ) { const Quadrature point_quadrature(p_unit); std::auto_ptr @@ -1351,7 +1415,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it mdata->use_mapping_q1_on_current_cell = false; - std::vector > &points = mdata->mapping_support_points; + std::vector > &points = mdata->mapping_support_points; compute_mapping_support_points (cell, points); this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit); @@ -1382,4 +1446,8 @@ MappingQ::clone () const // explicit instantiation template class MappingQ; +#if deal_II_dimension < 3 +template class MappingQ; +#endif + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/fe/mapping_q_eulerian.cc b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc index 4c85e5eb4a..2e35b3b463 100644 --- a/deal.II/deal.II/source/fe/mapping_q_eulerian.cc +++ b/deal.II/deal.II/source/fe/mapping_q_eulerian.cc @@ -214,4 +214,16 @@ template class MappingQEulerian >; template class MappingQEulerian; #endif +#if deal_II_dimension != 3 +template class MappingQEulerian, + deal_II_dimension+1>; + +# ifdef DEAL_II_USE_PETSC +template class MappingQEulerian; +# endif + +#endif + + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc index 2f13093807..ce03ff0286 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -699,6 +699,40 @@ HyperBallBoundary<3>::get_intermediate_points_on_quad ( #endif +#if deal_II_dimension == 2 + +template <> +void +HyperBallBoundary<2,3>::get_intermediate_points_on_quad ( + const Triangulation<2,3>::quad_iterator &quad, + std::vector > &points) const +{ + if (points.size()==1) + points[0]=get_new_point_on_quad(quad); + else + { + unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); + Assert(points.size()==m*m, ExcInternalError()); + + std::vector > lp0(m); + std::vector > lp1(m); + + get_intermediate_points_on_line(quad->line(0), lp0); + get_intermediate_points_on_line(quad->line(1), lp1); + + std::vector > lps(m); + for (unsigned int i=0; i void HyperBallBoundary::get_intermediate_points_on_quad ( diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index d56923130d..845d23e2e3 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -717,7 +717,7 @@ build_one_patch (const std::pair *cell_and_index, if (curved_cell_region==curved_inner_cells || (curved_cell_region==curved_boundary && cell_and_index->first->at_boundary())) { - Assert(patch.space_dim==dim, ExcInternalError()); + Assert(patch.space_dim==DH::space_dimension, ExcInternalError()); const std::vector > & q_points=fe_patch_values.get_quadrature_points(); // resize the patch.data member // in order to have enough memory @@ -731,7 +731,7 @@ build_one_patch (const std::pair *cell_and_index, // copy points to patch.data for (unsigned int i=0; i +
  • New: MappingQ and MappingQEulerian now support order > 1 also in + codimension one. Step-34 has been modified to show how this works. +
    + (Luca Heltai 2010/07/23-27) +

    +
  • New: QGaussOneOverR now has a new constructor for arbitrary quadrature points and not only the vertices of the reference cell.
    -- 2.39.5