From b94df41b63caf404d541f8b64ae7e5210a6a07a5 Mon Sep 17 00:00:00 2001 From: buerg Date: Mon, 10 Jun 2013 19:53:19 +0000 Subject: [PATCH] Add support for faces with non-standard orientation for lowest-order FE_Nedelec<3>. git-svn-id: https://svn.dealii.org/trunk@29803 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_poly_tensor.cc | 118 +++++++++++++++++++++++----- 1 file changed, 97 insertions(+), 21 deletions(-) diff --git a/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/source/fe/fe_poly_tensor.cc index 8d5562849c..7f29e909db 100644 --- a/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/source/fe/fe_poly_tensor.cc @@ -37,9 +37,9 @@ namespace * given cell. */ void - get_face_sign_change (const Triangulation<1>::cell_iterator &, - const unsigned int , - std::vector &face_sign) + get_face_sign_change_rt (const Triangulation<1>::cell_iterator &, + const unsigned int , + std::vector &face_sign) { // nothing to do in 1d std::fill (face_sign.begin (), face_sign.end (), 1.0); @@ -48,9 +48,9 @@ namespace void - get_face_sign_change (const Triangulation<2>::cell_iterator &cell, - const unsigned int dofs_per_face, - std::vector &face_sign) + get_face_sign_change_rt (const Triangulation<2>::cell_iterator &cell, + const unsigned int dofs_per_face, + std::vector &face_sign) { const unsigned int dim = 2; const unsigned int spacedim = 2; @@ -83,13 +83,72 @@ namespace void - get_face_sign_change (const Triangulation<3>::cell_iterator &/*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) + get_face_sign_change_rt (const Triangulation<3>::cell_iterator & cell, + const unsigned int dofs_per_face, + std::vector &face_sign) { std::fill (face_sign.begin (), face_sign.end (), 1.0); //TODO: think about what it would take here } + + void + get_face_sign_change_nedelec (const Triangulation<1>::cell_iterator &, + const unsigned int , + std::vector &face_sign) + { + // nothing to do in 1d + std::fill (face_sign.begin (), face_sign.end (), 1.0); + } + + + + void + get_face_sign_change_nedelec (const Triangulation<2>::cell_iterator &cell, + const unsigned int dofs_per_face, + std::vector &face_sign) + { + std::fill (face_sign.begin (), face_sign.end (), 1.0); +//TODO: think about what it would take here + } + + + void + get_face_sign_change_nedelec (const Triangulation<3>::cell_iterator & cell, + const unsigned int dofs_per_face, + std::vector &face_sign) + { + const unsigned int dim = 3; + std::fill (face_sign.begin (), face_sign.end (), 1.0); + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + { + if (!(cell->face (f)->at_boundary ())) + { + const bool face_orientation = cell->face_orientation (f); + const bool face_flip = cell->face_flip (f); + const bool face_rotation = cell->face_rotation (f); +//TODO: This is probably only going to work for those elements for which all dofs are face dofs + if (face_flip) + { + face_sign[GeometryInfo::face_to_cell_lines (f, 2, face_orientation, face_flip, face_rotation)] = -1.0; + face_sign[GeometryInfo::face_to_cell_lines (f, 3, face_orientation, face_flip, face_rotation)] = -1.0; + + if (!face_rotation) + { + face_sign[GeometryInfo::face_to_cell_lines (f, 0, face_orientation, face_flip, face_rotation)] = -1.0; + face_sign[GeometryInfo::face_to_cell_lines (f, 1, face_orientation, face_flip, face_rotation)] = -1.0; + } + } + + else + if (face_rotation) + { + face_sign[GeometryInfo::face_to_cell_lines (f, 0, face_orientation, face_flip, face_rotation)] = -1.0; + face_sign[GeometryInfo::face_to_cell_lines (f, 1, face_orientation, face_flip, face_rotation)] = -1.0; + } + } + } + } } @@ -383,13 +442,19 @@ FE_PolyTensor::fill_fe_values ( // Compute eventual sign changes depending on the neighborhood // between two faces. std::vector sign_change (this->dofs_per_cell, 1.0); - get_face_sign_change (cell, this->dofs_per_face, sign_change); - + + if (mapping_type == mapping_raviart_thomas) + get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + + else + if (mapping_type == mapping_nedelec) + get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); // for Piola mapping, the similarity // concept cannot be used because of // possible sign changes from one cell to // the next. - if ( (mapping_type == mapping_piola) || (mapping_type == mapping_raviart_thomas) ) + if ( (mapping_type == mapping_piola) || (mapping_type == mapping_raviart_thomas) + || (mapping_type == mapping_nedelec)) if (cell_similarity == CellSimilarity::translation) cell_similarity = CellSimilarity::none; @@ -444,7 +509,7 @@ FE_PolyTensor::fill_fe_values ( for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_values(first+d,k) = shape_values[k][d]; + data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; break; } @@ -530,8 +595,7 @@ FE_PolyTensor::fill_fe_values ( for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_gradients[first + d][k] = shape_grads1[k][d]; - // then copy over to target: + data.shape_gradients[first + d][k] = sign_change[i] * shape_grads1[k][d]; break; } @@ -591,7 +655,13 @@ FE_PolyTensor::fill_fe_face_values ( // Compute eventual sign changes depending // on the neighborhood between two faces. std::vector sign_change (this->dofs_per_cell, 1.0); - get_face_sign_change (cell, this->dofs_per_face, sign_change); + + if (mapping_type == mapping_raviart_thomas) + get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + + else + if (mapping_type == mapping_nedelec) + get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { @@ -646,7 +716,7 @@ FE_PolyTensor::fill_fe_face_values ( for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_values(first+d,k) = shape_values[k][d]; + data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; break; } @@ -732,8 +802,8 @@ FE_PolyTensor::fill_fe_face_values ( for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_gradients[first + d][k] = shape_grads1[k][d]; - // then copy over to target: + data.shape_gradients[first + d][k] = sign_change[i] * shape_grads1[k][d]; + break; } @@ -795,7 +865,13 @@ FE_PolyTensor::fill_fe_subface_values ( // Compute eventual sign changes depending // on the neighborhood between two faces. std::vector sign_change (this->dofs_per_cell, 1.0); - get_face_sign_change (cell, this->dofs_per_face, sign_change); + + if (mapping_type == mapping_raviart_thomas) + get_face_sign_change_rt (cell, this->dofs_per_face, sign_change); + + else + if (mapping_type == mapping_nedelec) + get_face_sign_change_nedelec (cell, this->dofs_per_face, sign_change); for (unsigned int i=0; idofs_per_cell; ++i) { @@ -849,7 +925,7 @@ FE_PolyTensor::fill_fe_subface_values ( for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_values(first+d,k) = shape_values[k][d]; + data.shape_values(first+d,k) = sign_change[i] * shape_values[k][d]; break; } -- 2.39.5