From ebabe2ecf8bbf592a2fd12c5e78136758b9bd592 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Thu, 15 Mar 2001 19:05:58 +0000 Subject: [PATCH] Remove alternative_normals_computation, hence could remove overload of compute_face_data, compute_fill_face and update_each functions. Make use_mapping_q_on_all_cells static const. Change straight_boundary from static member variable to static function variables. git-svn-id: https://svn.dealii.org/trunk@4224 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/mapping_q.cc | 275 +++++++------------------ 1 file changed, 75 insertions(+), 200 deletions(-) diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 960fba2980..d5f5fcff52 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -25,11 +25,13 @@ #include +template +const bool MappingQ::use_mapping_q_on_all_cells; + template MappingQ::InternalData::InternalData (const unsigned int n_shape_functions): MappingQ1::InternalData(n_shape_functions), -//TODO: in 1d, use_mapping_q1_on_current_cell is always false, but should be true. use_mapping_q1_on_current_cell(false), mapping_q1_data(1 << dim) { @@ -50,11 +52,7 @@ MappingQ<1>::MappingQ (const unsigned int): n_outer(0), tensor_pols(0), n_shape_functions(2), - renumber(0), -//TODO: why have two ways to compute? if they both work, choose one and remove the other - alternative_normals_computation(false), -//TODO: remove use_mapping_q_on_all_cells as it is set to false in the constructor and never set again - use_mapping_q_on_all_cells(false) + renumber(0) {} @@ -86,11 +84,7 @@ MappingQ::MappingQ (const unsigned int p): :8+12*(degree-1)+6*(degree-1)*(degree-1)), tensor_pols(0), n_shape_functions(power(degree+1,dim)), - renumber(0), -//TODO: why have two ways to compute? if they both work, choose one and remove the other - alternative_normals_computation(false), -//TODO: remove use_mapping_q_on_all_cells as it is set to false in the constructor and never set again - use_mapping_q_on_all_cells(false) + renumber(0) { // Construct the tensor product // polynomials used as shape @@ -191,21 +185,6 @@ MappingQ::compute_shapes_virtual (const std::vector > &unit_poin -template -UpdateFlags -MappingQ::update_each (const UpdateFlags in) const -{ - UpdateFlags out=MappingQ1::update_each(in); - - if (in & update_normal_vectors) - if (alternative_normals_computation) - out |= update_covariant_transformation; - - return out; -} - - - template Mapping::InternalDataBase* MappingQ::get_data (const UpdateFlags update_flags, @@ -232,10 +211,9 @@ MappingQ::get_face_data (const UpdateFlags update_flags, compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) - MappingQ1::compute_face_data (update_flags, q, - quadrature.n_quadrature_points, - data->mapping_q1_data); - + compute_face_data (update_flags, q, + quadrature.n_quadrature_points, + data->mapping_q1_data); return data; } @@ -251,71 +229,14 @@ MappingQ::get_subface_data (const UpdateFlags update_flags, compute_face_data (update_flags, q, quadrature.n_quadrature_points, *data); if (!use_mapping_q_on_all_cells) - MappingQ1::compute_face_data (update_flags, q, - quadrature.n_quadrature_points, - data->mapping_q1_data); - + compute_face_data (update_flags, q, + quadrature.n_quadrature_points, + data->mapping_q1_data); return data; } -template -void -MappingQ::compute_face_data (UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - MappingQ1::InternalData &mapping_q1_data) const -{ - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - InternalData &data = dynamic_cast (mapping_q1_data); - - MappingQ1::compute_face_data(update_flags, q, - n_original_q_points, data); - -//TODO: externalize this to a proper template function. or: scrap the whole function for 1d since it doesn't make much sense anyway? -#if (deal_II_dimension>1) - if ((data.update_flags & update_normal_vectors) - && alternative_normals_computation) - { - const unsigned int nfaces = GeometryInfo::faces_per_cell; - data.unit_normals.resize(nfaces); - std::vector > n(nfaces); - if (dim==2) - { - n[0][1]=-1; - n[1][0]=1; - n[2][1]=1; - n[3][0]=-1; - } - else if (dim==3) - { - n[0][1]=-1; - n[1][1]=1; - n[2][2]=-1; - n[3][0]=1; - n[4][2]=1; - n[5][0]=-1; - } - else - Assert(false, ExcNotImplemented()); - - for (unsigned int i=0; i void MappingQ::fill_fe_values (const DoFHandler::cell_iterator &cell, @@ -341,12 +262,14 @@ MappingQ::fill_fe_values (const DoFHandler::cell_iterator &cell, // depending on this result, use // this or the other data object // for the mapping + MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) - MappingQ1::fill_fe_values(cell, q, data.mapping_q1_data, - quadrature_points, JxW_values); + p_data=&data.mapping_q1_data; else - MappingQ1::fill_fe_values(cell, q, data, - quadrature_points, JxW_values); + p_data=&data; + + MappingQ1::fill_fe_values(cell, q, *p_data, + quadrature_points, JxW_values); } @@ -373,7 +296,7 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato // treated by a reduced Q1 mapping, // e.g. if the cell is in the // interior of the domain -//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell? +//TODO: (comment needed) shouldn't we ask whether the face is at the boundary, rather than the cell? data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -383,18 +306,17 @@ MappingQ::fill_fe_face_values (const typename DoFHandler::cell_iterato // depending on this result, use // this or the other data object // for the mapping + MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) - MappingQ1::compute_fill_face (cell, face_no, false, - npts, offset, q.get_weights(), - data.mapping_q1_data, - quadrature_points, JxW_values, - exterior_forms, normal_vectors); + p_data=&data.mapping_q1_data; else - compute_fill_face (cell, face_no, false, - npts, offset, q.get_weights(), - data, - quadrature_points, JxW_values, - exterior_forms, normal_vectors); + p_data=&data; + + compute_fill_face (cell, face_no, false, + npts, offset, q.get_weights(), + *p_data, + quadrature_points, JxW_values, + exterior_forms, normal_vectors); } @@ -421,7 +343,7 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter // treated by a reduced Q1 mapping, // e.g. if the cell is in the // interior of the domain -//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell? +//TODO: (comment needed) shouldn't we ask whether the face is at the boundary, rather than the cell? data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -432,18 +354,17 @@ MappingQ::fill_fe_subface_values (const typename DoFHandler::cell_iter // depending on this result, use // this or the other data object // for the mapping + MappingQ1::InternalData *p_data=0; if (data.use_mapping_q1_on_current_cell) - MappingQ1::compute_fill_face (cell, face_no, true, - npts, offset, q.get_weights(), - data.mapping_q1_data, - quadrature_points, JxW_values, - exterior_forms, normal_vectors); + p_data=&data.mapping_q1_data; else - compute_fill_face (cell, face_no, true, - npts, offset, q.get_weights(), - data, - quadrature_points, JxW_values, - exterior_forms, normal_vectors); + p_data=&data; + + compute_fill_face (cell, face_no, true, + npts, offset, q.get_weights(), + *p_data, + quadrature_points, JxW_values, + exterior_forms, normal_vectors); } @@ -508,7 +429,8 @@ MappingQ::set_laplace_on_quad_vector(std::vector > &loq // not precomputed, then do so now if (dim==2) compute_laplace_vector(loqvs); -//TODO: what if dim==3? + else + Assert(false, ExcNotImplemented()); } // the sum of weights of the points @@ -756,8 +678,11 @@ MappingQ::compute_mapping_support_points( if (use_mapping_q_on_all_cells || cell->has_boundary_lines()) compute_support_points_laplace(cell, a); // compute_support_points_simple(cell, a); -//TODO: can we delete the previous line? - +//TODO: (later) can we delete the previous line? + // keep last line for easy + // switching between the two + // possible cases of computing the + // support points. else // otherwise: use a Q1 mapping // for which the mapping shape @@ -802,7 +727,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulation: // points located on // the boundary faces add_line_support_points (cell, a); - add_face_support_points (cell, a); + add_quad_support_points (cell, a); apply_laplace_vector (laplace_on_hex_vector, a); break; @@ -814,14 +739,14 @@ MappingQ::compute_support_points_laplace(const typename Triangulation: -//TODO: remove the following function altogether +//TODO: (later) remove the following function altogether template void MappingQ::compute_support_points_simple(const typename Triangulation::cell_iterator &cell, std::vector > &a) const { + Assert(a.size()==0, ExcInternalError()); // the vertices first -// TODO: check for size of 'a' first? for (unsigned int i=0; i::vertices_per_cell; ++i) a.push_back(cell->vertex(i)); @@ -856,10 +781,6 @@ MappingQ::compute_support_points_simple(const typename Triangulation:: // the midpoint between // the middle point and // the 4 vertices - - // TODO: better position of - // points: transform them by - // a Q2 transformation. for (unsigned int i=0; i::vertices_per_cell; ++i) a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.); break; @@ -919,6 +840,7 @@ void MappingQ::add_line_support_points (const Triangulation::cell_iterator &cell, std::vector > &a) const { + static const StraightBoundary straight_boundary; // if we only need the midpoint, // then ask for it. if (degree==2) @@ -968,10 +890,10 @@ MappingQ::add_line_support_points (const Triangulation::cell_iterator #if deal_II_dimension==3 -//TODO: rename function to add_quad_support_points, to unify notation + template<> void -MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell, +MappingQ<3>::add_quad_support_points(const Triangulation<3>::cell_iterator &cell, std::vector > &a) const { const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, @@ -994,6 +916,14 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell {9,5,10,1}, {2,10,6,11}, {8,7,11,3}}; + + static StraightBoundary<3> straight_boundary; + // used if face quad at boundary or + // entirely in the interior of the + // domain + std::vector > quad_points ((degree-1)*(degree-1)); + // used if only one line of face quad is at boundary + std::vector > b(4*degree); // loop over all faces and collect points on them @@ -1017,9 +947,6 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell // on it if (face->at_boundary()) { -//TODO: move this variable out of the inner loop - std::vector > quad_points ((degree-1)*(degree-1)); - face->get_triangulation().get_boundary(face->boundary_indicator()) .get_intermediate_points_on_quad (face, quad_points); a.insert (a.end(), quad_points.begin(), quad_points.end()); @@ -1043,14 +970,14 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell // separately if (lines_at_boundary>0) { + // b is of size + // 4*degree, make sure + // that this is the + // right size + Assert(b.size()==vertices_per_face+lines_per_face*(degree-1), + ExcDimensionMismatch(b.size(), vertices_per_face+lines_per_face*(degree-1))); + // sort the points into b -//TODO: this is not thread-safe!!! b might be used for objects with -//TODO: different degrees at the same time! - static std::vector > b; - b.resize(4*degree); - Assert(4*degree==vertices_per_face+lines_per_face*(degree-1), - ExcDimensionMismatch(4*degree, - vertices_per_face+lines_per_face*(degree-1))); for (unsigned int i=0; i::add_face_support_points(const Triangulation<3>::cell_iterator &cell for (unsigned int j=0; j::add_face_support_points(const Triangulation<3>::cell_iterator &cell // intermediate points // from a straight // boundary object -//TODO: move this variable out of the loop - std::vector > quad_points ((degree-1)*(degree-1)); - straight_boundary.get_intermediate_points_on_quad (face, quad_points); a.insert (a.end(), quad_points.begin(), quad_points.end()); } @@ -1085,17 +1014,16 @@ MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell #endif + template void -MappingQ::add_face_support_points(const typename Triangulation::cell_iterator &, +MappingQ::add_quad_support_points(const typename Triangulation::cell_iterator &, std::vector > &) const { Assert(false, ExcInternalError()); } - - #if deal_II_dimension==3 template <> @@ -1103,6 +1031,8 @@ void MappingQ<3>::fill_quad_support_points_simple (const Triangulation<3>::cell_iterator &cell, std::vector > &a) const { + static StraightBoundary<3> straight_boundary; + const Boundary<3> *boundary = 0; std::vector > quad_points; @@ -1131,61 +1061,6 @@ MappingQ::fill_quad_support_points_simple (const Triangulation::cell_i {} - -//TODO: remove call of cross_product for dim==2 -#if deal_II_dimension==2 - -void cross_product (Tensor<1,2> &, const Tensor<1,2> &, const Tensor<1,2> &) -{ - Assert(false, ExcInternalError()); -} - -#endif - -template -void -MappingQ::compute_fill_face (const typename DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const bool is_subface, - const unsigned int npts, - const unsigned int offset, - const std::vector &weights, - MappingQ1::InternalData &mapping_q1_data, - std::vector > &quadrature_points, - std::vector &JxW_values, - std::vector > &boundary_forms, - std::vector > &normal_vectors) const -{ - MappingQ1::compute_fill_face (cell, face_no, is_subface, - npts, - offset, - weights, - mapping_q1_data, - quadrature_points, - JxW_values, - boundary_forms, - normal_vectors); - - const UpdateFlags update_flags(mapping_q1_data.current_update_flags()); - - if ((update_flags & update_normal_vectors) - && alternative_normals_computation) - { - InternalData *data_ptr = dynamic_cast (&mapping_q1_data); - Assert(data_ptr!=0, ExcInternalError()); - InternalData &data=*data_ptr; - - transform_covariant(normal_vectors, - data.unit_normals[face_no], - data, 0); - - for (unsigned int i=0; i void MappingQ::transform_covariant (std::vector > &dst, -- 2.39.5