From: leicht Date: Fri, 9 Feb 2007 07:13:14 +0000 (+0000) Subject: Use new interfaces to GeometryInfo and QProjector where necessary. Furthermore, let... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=149eca76d35488eb84954fb67570c9b2493a96a8;p=dealii-svn.git Use new interfaces to GeometryInfo and QProjector where necessary. Furthermore, let the FiniteElements tell, how dofs on faces and lines in 3D have to be permuted if the face has non-standard orientation, flip or rotation and/or the lines have non-standard orientation. Let DoFAccessor use this information if asked for the dofs on a cell. Currently this is only implemented for FE_Q and FESystem, for the other FEs there is a TODO. This is the last partial merge from branch_general_meshes. git-svn-id: https://svn.dealii.org/trunk@14448 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h index 612dd8e571..4660850f7a 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -481,9 +481,22 @@ DoFObjectAccessor<2,DH>::get_dof_indices (std::vector &dof_indices for (unsigned int vertex=0; vertex<4; ++vertex) for (unsigned int d=0; dvertex_dof_index(vertex,d,fe_index); + // now copy dof numbers from the line. for + // lines with the wrong orientation (which + // might occur in 3d), we have already made + // sure that we're ok by picking the correct + // vertices (this happens automatically in + // the vertex() function). however, if the + // line is in wrong orientation, we look at + // it in flipped orientation and we will have + // to adjust the shape function indices that + // we see to correspond to the correct + // (face-local) ordering. for (unsigned int line=0; line<4; ++line) for (unsigned int d=0; dline(line)->dof_index(d,fe_index); + *next++ = this->line(line)->dof_index(this->dof_handler->get_fe()[fe_index]. + adjust_line_dof_index_for_line_orientation(d, + this->line_orientation(line)),fe_index); for (unsigned int d=0; ddof_index(d,fe_index); } @@ -585,29 +598,42 @@ DoFObjectAccessor<3,DH>::get_dof_indices (std::vector &dof_indices for (unsigned int vertex=0; vertex<8; ++vertex) for (unsigned int d=0; dvertex_dof_index(vertex,d,fe_index); + // now copy dof numbers from the line. for + // lines with the wrong orientation, we have + // already made sure that we're ok by picking + // the correct vertices (this happens + // automatically in the vertex() + // function). however, if the line is in + // wrong orientation, we look at it in + // flipped orientation and we will have to + // adjust the shape function indices that we + // see to correspond to the correct + // (cell-local) ordering. for (unsigned int line=0; line<12; ++line) for (unsigned int d=0; dline(line)->dof_index(d,fe_index); - + *next++ = this->line(line)->dof_index(this->dof_handler->get_fe()[fe_index]. + adjust_line_dof_index_for_line_orientation(d, + this->line_orientation(line)),fe_index); // now copy dof numbers from the face. for // faces with the wrong orientation, we // have already made sure that we're ok by // picking the correct lines and vertices // (this happens automatically in the - // line() and vertex() functions. however, + // line() and vertex() functions). however, // if the face is in wrong orientation, we // look at it in flipped orientation and we // will have to adjust the shape function // indices that we see to correspond to the - // correct (cell-local) ordering. + // correct (cell-local) ordering. The same + // applies, if the face_rotation or + // face_orientation is non-standard for (unsigned int quad=0; quad<6; ++quad) - if (this->face_orientation(quad)) - for (unsigned int d=0; dquad(quad)->dof_index(d,fe_index); - else - for (unsigned int d=0; dquad(quad)->dof_index(this->dof_handler->get_fe()[fe_index]. - adjust_quad_dof_index_for_face_orientation(d),fe_index); + for (unsigned int d=0; dquad(quad)->dof_index(this->dof_handler->get_fe()[fe_index]. + adjust_quad_dof_index_for_face_orientation(d, + this->face_orientation(quad), + this->face_flip(quad), + this->face_rotation(quad)),fe_index); for (unsigned int d=0; ddof_index(d,fe_index); } diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index e3cf286baa..2c61f473a1 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1201,14 +1201,33 @@ class FiniteElement : public Subscriptor, * faces (quads) have to be permuted in * order to be combined with the correct * shape functions. Given a local dof @p - * index on a quad, return the local - * index, if the face has non-standard - * face_orientation. In 2D and 1D there - * is no need for permutation and - * consequently, and exception is thrown. - */ - unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index) const; - + * index on a quad, return the local index, + * if the face has non-standard + * face_orientation, face_flip or + * face_rotation. In 2D and 1D there is no + * need for permutation and consequently + * an exception is thrown. + */ + unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) const; + + /** + * For lines with non-standard + * line_orientation in 3D, the dofs on + * lines have to be permuted in order to be + * combined with the correct shape + * functions. Given a local dof @p index on + * a line, return the local index, if the + * line has non-standard + * line_orientation. In 2D and 1D there is + * no need for permutation, so the given + * index is simply returned. + */ + unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index, + const bool line_orientation) const; + /** * Return in which of the vector * components of this finite @@ -2079,8 +2098,35 @@ class FiniteElement : public Subscriptor, * i.e. old_index + shift = * new_index. In 2D and 1D there is * no need for permutation so the vector is + * empty. In 3D it has the size of + * dofs_per_quad * 8 , where 8 is + * the number of orientations, a face can + * be in (all comibinations of the three + * bool flags face_orientation, face_flip + * and face_rotation). + * + * The standard implementation fills this + * with zeros, i.e. no permuatation at + * all. Derived finite element classes have + * to fill this Table with the correct + * values. + */ + Table<2,int> adjust_quad_dof_index_for_face_orientation_table; + + /** + * For lines with non-standard + * line_orientation in 3D, the dofs on + * lines have to be permuted in + * order to be combined with the correct + * shape functions. Given a local dof @p + * index on a line, return the shift in the + * local index, if the line has + * non-standard line_orientation, + * i.e. old_index + shift = + * new_index. In 2D and 1D there is + * no need for permutation so the vector is * empty. In 3D it has the size of @p - * dofs_per_quad. + * dofs_per_line. * * The standard implementation fills this * with zeros, i.e. no permuatation at @@ -2088,7 +2134,7 @@ class FiniteElement : public Subscriptor, * to fill this vector with the correct * values. */ - std::vector adjust_quad_dof_index_for_face_orientation_table; + std::vector adjust_line_dof_index_for_line_orientation_table; private: /** diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 6eebc85d31..de5c14aae9 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -557,7 +557,9 @@ class FiniteElementData */ unsigned int face_to_cell_index (const unsigned int face_index, const unsigned int face, - const bool orientation = true) const; + const bool face_orientation = true, + const bool face_flip = false, + const bool face_rotation = false) const; /** * @deprecated This function is diff --git a/deal.II/deal.II/include/fe/fe_poly.templates.h b/deal.II/deal.II/include/fe/fe_poly.templates.h index 450e1715a1..1f9c0d586a 100644 --- a/deal.II/deal.II/include/fe/fe_poly.templates.h +++ b/deal.II/deal.II/include/fe/fe_poly.templates.h @@ -294,7 +294,10 @@ FE_Poly::fill_fe_face_values (const Mapping &ma const typename QProjector::DataSetDescriptor dsd; const typename QProjector::DataSetDescriptor offset - = dsd.face (face, cell->face_orientation(face), + = dsd.face (face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), quadrature.n_quadrature_points); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -345,8 +348,11 @@ FE_Poly::fill_fe_subface_values (const Mapping const typename QProjector::DataSetDescriptor dsd; const typename QProjector::DataSetDescriptor offset - = dsd.subface (face, subface, cell->face_orientation(face), - quadrature.n_quadrature_points); + = dsd.subface (face, subface, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), + quadrature.n_quadrature_points); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index 78ad32b969..de75f523ca 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -249,29 +249,42 @@ DoFCellAccessor >::update_cell_dof_indices_cache () const for (unsigned int vertex=0; vertex<8; ++vertex) for (unsigned int d=0; dline(line)->dof_index(d); - + *next++ = this->line(line)->dof_index(this->dof_handler->get_fe(). + adjust_line_dof_index_for_line_orientation(d, + this->line_orientation(line))); // now copy dof numbers from the face. for // faces with the wrong orientation, we // have already made sure that we're ok by // picking the correct lines and vertices // (this happens automatically in the - // line() and vertex() functions. however, + // line() and vertex() functions). however, // if the face is in wrong orientation, we // look at it in flipped orientation and we // will have to adjust the shape function // indices that we see to correspond to the - // correct (cell-local) ordering. + // correct (cell-local) ordering. The same + // applies, if the face_rotation or + // face_orientation is non-standard for (unsigned int quad=0; quad<6; ++quad) - if (this->face_orientation(quad)) - for (unsigned int d=0; dquad(quad)->dof_index(d); - else - for (unsigned int d=0; dquad(quad)->dof_index(this->dof_handler->get_fe(). - adjust_quad_dof_index_for_face_orientation(d)); + for (unsigned int d=0; dquad(quad)->dof_index(this->dof_handler->get_fe(). + adjust_quad_dof_index_for_face_orientation(d, + this->face_orientation(quad), + this->face_flip(quad), + this->face_rotation(quad))); for (unsigned int d=0; d::FiniteElement ( FiniteElementData (fe_data), cached_primitivity(false), adjust_quad_dof_index_for_face_orientation_table (dim == 3 ? - this->dofs_per_quad : 0, - 0), + this->dofs_per_quad : 0 , + dim==3 ? 8 : 0), + adjust_line_dof_index_for_line_orientation_table (dim == 3 ? + this->dofs_per_line : 0), system_to_base_table(this->dofs_per_cell), face_system_to_base_table(this->dofs_per_face), component_to_base_table (this->components, @@ -206,6 +208,7 @@ FiniteElement::FiniteElement ( // changed by constructor of // derived class. first_block_of_base_table.resize(1,0); + adjust_quad_dof_index_for_face_orientation_table.fill(0); } @@ -338,21 +341,27 @@ FiniteElement::component_to_block_index (const unsigned int index) const template unsigned int -FiniteElement::adjust_quad_dof_index_for_face_orientation (const unsigned int) const +FiniteElement::adjust_quad_dof_index_for_face_orientation (const unsigned int index, + const bool, + const bool, + const bool) const { // general template for 1D and 2D: not // implemented. in fact, the function // shouldn't even be called unless we are // in 3d, so throw an internal error Assert (false, ExcInternalError()); - return deal_II_numbers::invalid_unsigned_int; + return index; } #if deal_II_dimension == 3 template <> unsigned int -FiniteElement<3>::adjust_quad_dof_index_for_face_orientation (const unsigned int index) const +FiniteElement<3>::adjust_quad_dof_index_for_face_orientation (const unsigned int index, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) const { // adjust dofs on 3d faces if the face is // flipped. note that we query a table that @@ -364,9 +373,41 @@ FiniteElement<3>::adjust_quad_dof_index_for_face_orientation (const unsigned int // the function should also not have been // called Assert (indexdofs_per_quad, ExcIndexRange(index,0,this->dofs_per_quad)); - Assert (adjust_quad_dof_index_for_face_orientation_table.size()==this->dofs_per_quad, + Assert (adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, ExcInternalError()); - return index+adjust_quad_dof_index_for_face_orientation_table[index]; + return index+adjust_quad_dof_index_for_face_orientation_table(index,4*face_orientation+2*face_flip+face_rotation); +} + +#endif + + + +template +unsigned int +FiniteElement::adjust_line_dof_index_for_line_orientation (const unsigned int index, + const bool) const +{ + // general template for 1D and 2D: do + // nothing. Do not throw an Assertion, + // however, in order to allow to call this + // function in 2D as well + return index; +} + +#if deal_II_dimension == 3 + +template <> +unsigned int +FiniteElement<3>::adjust_line_dof_index_for_line_orientation (const unsigned int index, + const bool line_orientation) const +{ + Assert (indexdofs_per_line, ExcIndexRange(index,0,this->dofs_per_line)); + Assert (adjust_line_dof_index_for_line_orientation_table.size()==this->dofs_per_line, + ExcInternalError()); + if (line_orientation) + return index; + else + return index+adjust_line_dof_index_for_line_orientation_table[index]; } #endif diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc index 34c68f8b0c..ce9c2b20c9 100644 --- a/deal.II/deal.II/source/fe/fe_abf.cc +++ b/deal.II/deal.II/source/fe/fe_abf.cc @@ -32,8 +32,9 @@ #include -//TODO: implement the adjust_quad_dof_index_for_face_orientation_table field, -//and write a test similar to bits/face_orientation_and_fe_q_01 +//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and +//adjust_line_dof_index_for_line_orientation_table fields, and write tests +//similar to bits/face_orientation_and_fe_q_* DEAL_II_NAMESPACE_OPEN @@ -686,8 +687,8 @@ FE_ABF::interpolate( double n_orient = (double) GeometryInfo::unit_normal_orientation[face]; for (unsigned int fp=0; fp < n_face_points; ++fp) { - // TODO: Check what the face_orientation has to be in 3D - unsigned int k = QProjector::DataSetDescriptor::face (face, false, n_face_points); + // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D + unsigned int k = QProjector::DataSetDescriptor::face (face, false, false, false, n_face_points); for (unsigned int i=0; i::unit_normal_direction[face]+offset); @@ -746,8 +747,8 @@ FE_ABF::interpolate( double n_orient = (double) GeometryInfo::unit_normal_orientation[face]; for (unsigned int fp=0; fp < n_face_points; ++fp) { - // TODO: Check what the face_orientation has to be in 3D - unsigned int k = QProjector::DataSetDescriptor::face (face, false, n_face_points); + // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D + unsigned int k = QProjector::DataSetDescriptor::face (face, false, false, false, n_face_points); for (unsigned int i=0; i::unit_normal_direction[face]][k + fp]; diff --git a/deal.II/deal.II/source/fe/fe_data.cc b/deal.II/deal.II/source/fe/fe_data.cc index 7f6c9e04da..2457fdaaa1 100644 --- a/deal.II/deal.II/source/fe/fe_data.cc +++ b/deal.II/deal.II/source/fe/fe_data.cc @@ -103,7 +103,9 @@ unsigned int FiniteElementData:: face_to_cell_index (const unsigned int face_index, const unsigned int face, - const bool orientation) const + const bool face_orientation, + const bool face_flip, + const bool face_rotation) const { Assert (face_index < this->dofs_per_face, ExcIndexRange(face_index, 0, this->dofs_per_face)); @@ -116,7 +118,10 @@ face_to_cell_index (const unsigned int face_index, // Vertex number on the face const unsigned int face_vertex = face_index / this->dofs_per_vertex; return face_index % this->dofs_per_vertex - + GeometryInfo::face_to_cell_vertices(face, face_vertex, orientation) + + GeometryInfo::face_to_cell_vertices(face, face_vertex, + face_orientation, + face_flip, + face_rotation) * this->dofs_per_vertex; } // Else, DoF on a line? @@ -127,7 +132,10 @@ face_to_cell_index (const unsigned int face_index, // Line number on the face const unsigned int face_line = index / this->dofs_per_line; return this->first_line_index + index % this->dofs_per_line - + GeometryInfo::face_to_cell_lines(face, face_line, orientation) + + GeometryInfo::face_to_cell_lines(face, face_line, + face_orientation, + face_flip, + face_rotation) * this->dofs_per_line; } // Else, DoF is on a quad diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index b2d0aa8cfa..ee6842a65e 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -24,6 +24,8 @@ #include +//TODO: implement the adjust_line_dof_index_for_line_orientation_table field, +//and write a test similar to bits/face_orientation_and_fe_q_02 DEAL_II_NAMESPACE_OPEN @@ -1130,7 +1132,10 @@ FE_Nedelec::fill_fe_face_values (const Mapping &mapp // faces are stored contiguously) const typename QProjector::DataSetDescriptor offset = (QProjector::DataSetDescriptor:: - face (face, cell->face_orientation(face), + face (face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), quadrature.n_quadrature_points)); // get the flags indicating the @@ -1157,7 +1162,7 @@ FE_Nedelec::fill_fe_face_values (const Mapping &mapp // ways Assert (fe_data.shape_values[0].size() == GeometryInfo::faces_per_cell * n_q_points * - (dim == 3 ? 2 : 1), + (dim == 3 ? 8 : 1), ExcInternalError()); std::vector > shape_values (n_q_points); @@ -1190,7 +1195,7 @@ FE_Nedelec::fill_fe_face_values (const Mapping &mapp // ways Assert (fe_data.shape_gradients[0].size() == GeometryInfo::faces_per_cell * n_q_points * - (dim == 3 ? 2 : 1), + (dim == 3 ? 8 : 1), ExcInternalError()); std::vector > shape_grads1 (n_q_points); @@ -1276,8 +1281,11 @@ FE_Nedelec::fill_fe_subface_values (const Mapping &m // faces are stored contiguously) const typename QProjector::DataSetDescriptor offset = (QProjector::DataSetDescriptor:: - subface (face, subface, cell->face_orientation(face), - quadrature.n_quadrature_points)); + subface (face, subface, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), + quadrature.n_quadrature_points)); // get the flags indicating the // fields that have to be filled diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index ba4fc66431..d0c82ab304 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -536,7 +536,10 @@ FE_PolyTensor::fill_fe_face_values ( const typename QProjector::DataSetDescriptor dsd; const typename QProjector::DataSetDescriptor offset - = dsd.face (face, cell->face_orientation(face), + = dsd.face (face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), n_q_points); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); @@ -717,8 +720,11 @@ FE_PolyTensor::fill_fe_subface_values ( const typename QProjector::DataSetDescriptor dsd; const typename QProjector::DataSetDescriptor offset - = dsd.subface (face, subface, cell->face_orientation(face), - n_q_points); + = dsd.subface (face, subface, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face), + n_q_points); const UpdateFlags flags(fe_data.update_once | fe_data.update_each); diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 9b4908e85e..7de581d316 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -810,18 +810,66 @@ template <> void FE_Q<3>::initialize_quad_dof_index_permutation () { - Assert (adjust_quad_dof_index_for_face_orientation_table.size()==this->dofs_per_quad, + Assert (adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, ExcInternalError()); - const unsigned int points = this->degree-1; - Assert(points*points==this->dofs_per_quad, ExcInternalError()); - + const unsigned int n=this->degree-1; + Assert(n*n==this->dofs_per_quad, ExcInternalError()); + + // alias for the table to fill + Table<2,int> &data=this->adjust_quad_dof_index_for_face_orientation_table; + + // the dofs on a face are connected to a n x + // n matrix. for example, for degree==4 we + // have the following dofs on a quad + + // ___________ + // | | + // | 6 7 8 | + // | | + // | 3 4 5 | + // | | + // | 0 1 2 | + // |___________| + // + // we have dof_no=i+n*j with index i in + // x-direction and index j in y-direction + // running from 0 to n-1. to extract i and j + // we can use i=dof_no%n and j=dof_no/n. i + // and j can then be used to construct the + // rotated and mirrored numbers. + + for (unsigned int local=0; localdofs_per_quad; ++local) // face support points are in lexicographic // ordering with x running fastest. invert // that (y running fastest) - this->adjust_quad_dof_index_for_face_orientation_table[local] - = (local%points)*points + local/points - local; + { + unsigned int i=local%n, + j=local/n; + + // face_orientation=false, face_flip=false, face_rotation=false + data(local,0)=j + i *n - local; + // face_orientation=false, face_flip=false, face_rotation=true + data(local,1)=(n-1-i) + j *n - local; + // face_orientation=false, face_flip=true, face_rotation=false + data(local,2)=(n-1-j) + (n-1-i)*n - local; + // face_orientation=false, face_flip=true, face_rotation=true + data(local,3)=i + (n-1-j)*n - local; + // face_orientation=true, face_flip=false, face_rotation=false + data(local,4)=0; + // face_orientation=true, face_flip=false, face_rotation=true + data(local,5)=j + (n-1-i)*n - local; + // face_orientation=true, face_flip=true, face_rotation=false + data(local,6)=(n-1-i) + (n-1-j)*n - local; + // face_orientation=true, face_flip=true, face_rotation=true + data(local,7)=(n-1-j) + i *n - local; + } + + // aditionally initialize reordering of line + // dofs + for (unsigned int i=0; idofs_per_line; ++i) + this->adjust_line_dof_index_for_line_orientation_table[i]=this->dofs_per_line-1-i - i; } #endif diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index e2faf68cf1..0ebed3355e 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -16,8 +16,9 @@ #include #include -//TODO: implement the adjust_quad_dof_index_for_face_orientation_table field, -//and write a test similar to bits/face_orientation_and_fe_q_01 +//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and +//adjust_line_dof_index_for_line_orientation_table fields, and write tests +//similar to bits/face_orientation_and_fe_q_* DEAL_II_NAMESPACE_OPEN diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index 859d86879e..7b89c2ac5d 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -27,8 +27,9 @@ #include #include -//TODO: implement the adjust_quad_dof_index_for_face_orientation_table field, -//and write a test similar to bits/face_orientation_and_fe_q_01 +//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and +//adjust_line_dof_index_for_line_orientation_table fields, and write tests +//similar to bits/face_orientation_and_fe_q_* DEAL_II_NAMESPACE_OPEN @@ -213,7 +214,7 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) { // Enter the support point // into the vector - this->generalized_support_points[current] = faces.point(current); + this->generalized_support_points[current] = faces.point(current+QProjector::DataSetDescriptor::face(0,true,false,false,n_face_points)); } } diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc index 957902c8f3..084ae9b768 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -25,8 +25,9 @@ #include -//TODO: implement the adjust_quad_dof_index_for_face_orientation_table field, -//and write a test similar to bits/face_orientation_and_fe_q_01 +//TODO: implement the adjust_quad_dof_index_for_face_orientation_table and +//adjust_line_dof_index_for_line_orientation_table fields, and write tests +//similar to bits/face_orientation_and_fe_q_* DEAL_II_NAMESPACE_OPEN @@ -328,7 +329,12 @@ FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) for (unsigned int k=0; kdofs_per_face*GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = faces.point(k); + this->generalized_support_points[k] = faces.point(k+QProjector + ::DataSetDescriptor::face(0, + true, + false, + false, + this->dofs_per_face)); current = this->dofs_per_face*GeometryInfo::faces_per_cell; } diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index ef70cb174b..2e68834440 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -1866,24 +1866,44 @@ FESystem<3>::initialize_quad_dof_index_permutation () { // the array into which we want to write // should have the correct size already. - Assert (adjust_quad_dof_index_for_face_orientation_table.size()==this->dofs_per_quad, + Assert (adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, ExcInternalError()); // to obtain the shifts for this composed - // element, concatenate the shift vectors of + // element, copy the shift information of // the base elements unsigned int index = 0; for (unsigned int b=0; b &temp + const Table<2,int> &temp = this->base_element(b).adjust_quad_dof_index_for_face_orientation_table; for (unsigned int c=0; cdofs_per_quad, + ExcInternalError()); + + // aditionally compose the permutation + // information for lines + Assert (adjust_line_dof_index_for_line_orientation_table.size()==this->dofs_per_line, + ExcInternalError()); + index = 0; + for (unsigned int b=0; b &temp2 + = this->base_element(b).adjust_line_dof_index_for_line_orientation_table; + for (unsigned int c=0; cdofs_per_line, ExcInternalError()); } diff --git a/deal.II/deal.II/source/fe/mapping_cartesian.cc b/deal.II/deal.II/source/fe/mapping_cartesian.cc index 0e8094359b..5a114c3bc5 100644 --- a/deal.II/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/deal.II/source/fe/mapping_cartesian.cc @@ -205,11 +205,15 @@ MappingCartesian::compute_fill (const typename Triangulation::cell_ite // called from FEFaceValues QProjector::DataSetDescriptor::face (face_no, cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), quadrature_points.size()) : // called from FESubfaceValues QProjector::DataSetDescriptor::subface (face_no, sub_no, cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), quadrature_points.size()) )); diff --git a/deal.II/deal.II/source/fe/mapping_q.cc b/deal.II/deal.II/source/fe/mapping_q.cc index 42af4fedd7..3b0f604b6d 100644 --- a/deal.II/deal.II/source/fe/mapping_q.cc +++ b/deal.II/deal.II/source/fe/mapping_q.cc @@ -383,7 +383,10 @@ MappingQ::fill_fe_face_values (const typename Triangulation::cell_iter this->compute_fill_face (cell, face_no, false, n_q_points, QProjector::DataSetDescriptor:: - face (face_no, cell->face_orientation(face_no), + face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), n_q_points), q.get_weights(), *p_data, @@ -446,6 +449,8 @@ MappingQ::fill_fe_subface_values (const typename Triangulation::cell_i QProjector::DataSetDescriptor:: subface (face_no, sub_no, cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), n_q_points), q.get_weights(), *p_data, @@ -1023,12 +1028,17 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // select the correct mappings // for the present face - const bool face_orientation = cell->face_orientation(face_no); + const bool face_orientation = cell->face_orientation(face_no), + face_flip = cell->face_flip (face_no), + face_rotation = cell->face_rotation (face_no); // some sanity checks up front for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( - GeometryInfo<3>::face_to_cell_vertices(face_no, i, face_orientation)), + GeometryInfo<3>::face_to_cell_vertices(face_no, i, + face_orientation, + face_flip, + face_rotation)), ExcInternalError()); // indices of the lines that @@ -1037,7 +1047,8 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // face_to_cell_lines for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( - face_no, i, face_orientation)), ExcInternalError()); + face_no, i, face_orientation, face_flip, face_rotation)), + ExcInternalError()); // if face at boundary, then // ask boundary object to @@ -1088,13 +1099,16 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // sort the points into b for (unsigned int i=0; i::face_to_cell_vertices(face_no, i, face_orientation)]; + b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i, + face_orientation, + face_flip, + face_rotation)]; for (unsigned int i=0; i::face_to_cell_lines( - face_no, i, face_orientation)*(degree-1)+j]; + face_no, i, face_orientation, face_flip, face_rotation)*(degree-1)+j]; // Now b includes the // right order of diff --git a/deal.II/deal.II/source/fe/mapping_q1.cc b/deal.II/deal.II/source/fe/mapping_q1.cc index e715a28f3a..991f94258a 100644 --- a/deal.II/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/deal.II/source/fe/mapping_q1.cc @@ -724,6 +724,8 @@ MappingQ1::fill_fe_face_values (const typename Triangulation::cell_ite n_q_points, DataSetDescriptor::face (face_no, cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), n_q_points), q.get_weights(), data, @@ -757,8 +759,10 @@ MappingQ1::fill_fe_subface_values (const typename Triangulation::cell_ compute_fill_face (cell, face_no, true, n_q_points, DataSetDescriptor::subface (face_no, sub_no, - cell->face_orientation(face_no), - n_q_points), + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + n_q_points), q.get_weights(), data, quadrature_points,