From: leicht Date: Fri, 9 Feb 2007 07:05:32 +0000 (+0000) Subject: Adjust QProjector to create all eight (instead of two) face-quadrature rules originat... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0b0c12670faa3a8385d9f28103d885c8abb7e507;p=dealii-svn.git Adjust QProjector to create all eight (instead of two) face-quadrature rules originating from all combinations of face_orientation, face_flip and face_rotation. This is a partial merge from branch_general_meshes. git-svn-id: https://svn.dealii.org/trunk@14447 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/qprojector.h b/deal.II/base/include/base/qprojector.h index 700d692adb..df3fdb3c72 100644 --- a/deal.II/base/include/base/qprojector.h +++ b/deal.II/base/include/base/qprojector.h @@ -169,7 +169,7 @@ class QProjector * later. * * @note In 3D, this function - * produces two sets of + * produces eight sets of * quadrature points for each * face, in order to cope * possibly different @@ -279,16 +279,15 @@ class QProjector static DataSetDescriptor cell (); /** - * Static function to - * generate an offset object - * for a given face of a cell - * with the given face - * orientation. This function - * of course is only allowed - * if dim>=2, and the - * face orientation is - * ignored if the space - * dimension equals 2. + * Static function to generate an + * offset object for a given face of a + * cell with the given face + * orientation, flip and rotation. This + * function of course is only allowed + * if dim>=2, and the face + * orientation, flip and rotation are + * ignored if the space dimension + * equals 2. * * The last argument denotes * the number of quadrature @@ -303,19 +302,20 @@ class QProjector DataSetDescriptor face (const unsigned int face_no, const bool face_orientation, + const bool face_flip, + const bool face_rotation, const unsigned int n_quadrature_points); /** - * Static function to - * generate an offset object - * for a given subface of a - * cell with the given face - * orientation. This function - * of course is only allowed - * if dim>=2, and the - * face orientation is - * ignored if the space - * dimension equals 2. + * Static function to generate an + * offset object for a given subface of + * a cell with the given face + * orientation, flip and rotation. This + * function of course is only allowed + * if dim>=2, and the face + * orientation, flip and rotation are + * ignored if the space dimension + * equals 2. * * The last argument denotes * the number of quadrature @@ -331,6 +331,8 @@ class QProjector subface (const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, + const bool face_flip, + const bool face_rotation, const unsigned int n_quadrature_points); /** @@ -382,6 +384,25 @@ class QProjector * orientations. */ static Quadrature<2> reflect (const Quadrature<2> &q); + + /** + * Given a quadrature object in + * 2d, rotate all quadrature + * points by @p n_times * 90 degrees + * counterclockwise + * and return them with their + * original weights. + * + * This function is necessary for + * projecting a 2d quadrature + * rule onto the faces of a 3d + * cube, since there we need all + * rotations to account for + * face_flip and face_rotation + * of non-standard faces. + */ + static Quadrature<2> rotate (const Quadrature<2> &q, + const unsigned int n_times); }; /*@}*/ diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 282ddd9bce..92c9f8939b 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -402,6 +402,46 @@ QProjector::reflect (const Quadrature<2> &q) } +template +Quadrature<2> +QProjector::rotate (const Quadrature<2> &q, + const unsigned int n_times) +{ + std::vector > q_points (q.n_quadrature_points); + std::vector weights (q.n_quadrature_points); + for (unsigned int i=0; i (q_points, weights); +} + + template <> void QProjector<1>::project_to_face (const Quadrature<0> &, @@ -805,49 +845,56 @@ QProjector<3>::project_to_all_faces (const SubQuadrature &quadrature) { const unsigned int dim = 3; - const SubQuadrature q_reflected = reflect (quadrature); + SubQuadrature q_reflected=reflect (quadrature); + SubQuadrature q[8]= + {quadrature, + rotate (quadrature,1), + rotate (quadrature,2), + rotate (quadrature,3), + q_reflected, + rotate (q_reflected,3), + rotate (q_reflected,2), + rotate (q_reflected,1)}; + + const unsigned int n_points = quadrature.n_quadrature_points, n_faces = GeometryInfo::faces_per_cell; // first fix quadrature points std::vector > q_points; - q_points.reserve(n_points * n_faces * 2); + q_points.reserve(n_points * n_faces * 8); std::vector > help(n_points); - // project to each face and append - // results - for (unsigned int face=0; face weights; - weights.reserve (n_points * n_faces * 2); - for (unsigned int face=0; face(q_points, weights); @@ -915,7 +962,16 @@ QProjector<3>::project_to_all_subfaces (const SubQuadrature &quadrature) { const unsigned int dim = 3; - const SubQuadrature q_reflected = reflect (quadrature); + SubQuadrature q_reflected=reflect (quadrature); + SubQuadrature q[8]= + {quadrature, + rotate (quadrature,1), + rotate (quadrature,2), + rotate (quadrature,3), + q_reflected, + rotate (q_reflected,3), + rotate (q_reflected,2), + rotate (q_reflected,1)}; const unsigned int n_points = quadrature.n_quadrature_points, n_faces = GeometryInfo::faces_per_cell, @@ -923,46 +979,39 @@ QProjector<3>::project_to_all_subfaces (const SubQuadrature &quadrature) // first fix quadrature points std::vector > q_points; - q_points.reserve (n_points * n_faces * subfaces_per_face * 2); + q_points.reserve (n_points * n_faces * subfaces_per_face * 8); std::vector > help(n_points); - // project to each face and copy - // results - for (unsigned int face=0; face weights; - weights.reserve (n_points * n_faces * subfaces_per_face * 2); - for (unsigned int face=0; face(q_points, weights); @@ -1033,30 +1082,51 @@ typename QProjector::DataSetDescriptor QProjector::DataSetDescriptor:: face (const unsigned int face_no, const bool face_orientation, + const bool face_flip, + const bool face_rotation, const unsigned int n_quadrature_points) { Assert (dim != 1, ExcInternalError()); Assert (face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + // in 3d, we have to account for faces that + // have non-standard face orientation, flip + // and rotation. thus, we have to store + // _eight_ data sets per face or subface + + // set up a table with the according offsets + // for non-standard orientation, first index: + // face_orientation (standard true=1), second + // index: face_flip (standard false=0), third + // index: face_rotation (standard false=0) + // + // note, that normally we should use the + // obvious offsets 0,1,2,3,4,5,6,7. However, + // prior to the changes enabling flipped and + // rotated faces, in many places of the + // library the convention was used, that the + // first dataset with offset 0 corresponds to + // a face in standard orientation. therefore + // we use the offsets 4,5,6,7,0,1,2,3 here to + // stick to that (implicit) convention + static const unsigned int offset[2][2][2]= + {{{4*GeometryInfo::faces_per_cell, 5*GeometryInfo::faces_per_cell}, // face_orientation=false; face_flip=false; face_rotation=false and true + {6*GeometryInfo::faces_per_cell, 7*GeometryInfo::faces_per_cell}}, // face_orientation=false; face_flip=true; face_rotation=false and true + {{0*GeometryInfo::faces_per_cell, 1*GeometryInfo::faces_per_cell}, // face_orientation=true; face_flip=false; face_rotation=false and true + {2*GeometryInfo::faces_per_cell, 3*GeometryInfo::faces_per_cell}}}; // face_orientation=true; face_flip=true; face_rotation=false and true + switch (dim) { case 1: case 2: return face_no * n_quadrature_points; - // in 3d, we have to - // account for faces - // that have reverse - // orientation. thus, - // we have to store - // _two_ data sets per - // face or subface + case 3: - return ((face_no + - (face_orientation == true ? - 0 : GeometryInfo::faces_per_cell)) - * n_quadrature_points); + return ((face_no + + offset[face_orientation][face_flip][face_rotation]) + * n_quadrature_points); default: Assert (false, ExcInternalError()); @@ -1072,16 +1142,54 @@ QProjector::DataSetDescriptor:: subface (const unsigned int face_no, const unsigned int subface_no, const bool face_orientation, + const bool face_flip, + const bool face_rotation, const unsigned int n_quadrature_points) { Assert (dim != 1, ExcInternalError()); Assert (face_no < GeometryInfo::faces_per_cell, ExcInternalError()); - // the trick with +1 prevents - // that we get a warning in 1d + // the trick with +1 prevents that we get a + // warning in 1d Assert (subface_no+1 < GeometryInfo::subfaces_per_face+1, ExcInternalError()); + // in 3d, we have to account for faces that + // have non-standard face orientation, flip + // and rotation. thus, we have to store + // _eight_ data sets per face or subface + + // set up a table with the according offsets + // for non-standard orientation, first index: + // face_orientation (standard true=1), second + // index: face_flip (standard false=0), third + // index: face_rotation (standard false=0) + // + // note, that normally we should use the + // obvious offsets 0,1,2,3,4,5,6,7. However, + // prior to the changes enabling flipped and + // rotated faces, in many places of the + // library the convention was used, that the + // first dataset with offset 0 corresponds to + // a face in standard orientation. therefore + // we use the offsets 4,5,6,7,0,1,2,3 here to + // stick to that (implicit) convention + static const unsigned int offset[2][2][2]= + {{ + // face_orientation=false; face_flip=false; face_rotation=false and true + {4*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face, + 5*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face}, + // face_orientation=false; face_flip=true; face_rotation=false and true + {6*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face, + 7*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face}}, + { + // face_orientation=true; face_flip=false; face_rotation=false and true + {0*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face, + 1*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face}, + // face_orientation=true; face_flip=true; face_rotation=false and true + {2*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face, + 3*GeometryInfo::faces_per_cell*GeometryInfo::subfaces_per_face}}}; + switch (dim) { case 1: @@ -1089,15 +1197,10 @@ subface (const unsigned int face_no, return ((face_no * GeometryInfo::subfaces_per_face + subface_no) * n_quadrature_points); - - // for 3d, same as above: case 3: return (((face_no * GeometryInfo::subfaces_per_face + subface_no) - + (face_orientation == true ? - 0 : - GeometryInfo::faces_per_cell * - GeometryInfo::subfaces_per_face) + + offset[face_orientation][face_flip][face_rotation] ) * n_quadrature_points); default: