From: daniel.arndt Date: Fri, 18 Jul 2014 09:52:45 +0000 (+0000) Subject: First working version for rotated periodic boundaries, tested for non-hp X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dd3128cb604e0ff99a82b96766b1f15510df7b79;p=dealii-svn.git First working version for rotated periodic boundaries, tested for non-hp git-svn-id: https://svn.dealii.org/branches/branch_periodic_bc@33195 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 5844243ead..7ce0590d18 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -1175,7 +1175,7 @@ namespace GridTools const int direction, std::vector > &matched_pairs, const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(), - const FullMatrix &matrix = FullMatrix(), + const FullMatrix &matrix = FullMatrix(IdentityMatrix(CONTAINER::space_dimension)), const std::vector &first_vector_components = std::vector()); diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index d045093d15..a2136cbd96 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -15,6 +15,7 @@ // --------------------------------------------------------------------- #include +#include #include #include #include @@ -1912,9 +1913,10 @@ namespace DoFTools Assert(face_1->at_boundary() && face_2->at_boundary(), ExcMessage ("Faces for periodicity constraints must be on the boundary")); - Assert(first_vector_components.size() == 0 || matrix.size(0) == matrix.size(1) == dim, + Assert(first_vector_components.empty() + || (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim), ExcMessage ("first_vector_components is nonempty, so matrix must be a rotation " - "matrix exactly of of size dim")); + "matrix exactly of size dim")); #ifdef DEBUG { @@ -1960,71 +1962,142 @@ namespace DoFTools // In the case that both faces have children, we loop over all // children and apply make_periodicty_constrains recursively: if (face_1->has_children() && face_2->has_children()) - { - Assert(face_1->n_children() == GeometryInfo::max_children_per_face && - face_2->n_children() == GeometryInfo::max_children_per_face, - ExcNotImplemented()); + { + Assert(face_1->n_children() == GeometryInfo::max_children_per_face && + face_2->n_children() == GeometryInfo::max_children_per_face, + ExcNotImplemented()); - for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - { - // Lookup the index for the second face - unsigned int j; - switch (dim) - { - case 2: - j = lookup_table_2d[face_flip][i]; - break; - case 3: - j = lookup_table_3d[face_orientation][face_flip][face_rotation][i]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + { + // Lookup the index for the second face + unsigned int j; + switch (dim) + { + case 2: + j = lookup_table_2d[face_flip][i]; + break; + case 3: + j = lookup_table_3d[face_orientation][face_flip][face_rotation][i]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } - make_periodicity_constraints (face_1->child(i), - face_2->child(j), - constraint_matrix, - component_mask, - face_orientation, - face_flip, - face_rotation, - matrix, - first_vector_components); - } + make_periodicity_constraints (face_1->child(i), + face_2->child(j), + constraint_matrix, + component_mask, + face_orientation, + face_flip, + face_rotation, + matrix, + first_vector_components); } + } else // otherwise at least one of the two faces is active and // we need to enter the constraints - { - // Build up the transformation matrix: + { + // Build up the transformation matrix: - if (first_vector_components.size() != 0 && matrix.size(0) == 0) - { - // We need to create a local transformation matrix... - AssertThrow(false, ExcNotImplemented()); +// if (first_vector_components.size() != 0 && matrix.size(0) == 0) +// { +// // We need to create a local transformation matrix... +// AssertThrow(false, ExcNotImplemented()); +// +// // FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)); +// } - // FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)); - } - else if (face_2->has_children() == false) - { - // TODO: This is horrible... - FullMatrix inverse; - inverse.invert(matrix); - set_periodicity_constraints(face_2, face_1, - inverse, - constraint_matrix, - component_mask, - face_orientation, face_flip, face_rotation); - } - else + FullMatrix transformation; + if (!first_vector_components.empty() && matrix.size(0) == (int)dim) + { + // The matrix describes a rotation and we have to build a + // transformation matrix, we assume that for a 0° rotation + // we would have to build the identity matrix + + const FiniteElement &fe1 + = face_1->get_fe(face_1->nth_active_fe_index(0)); + + Quadrature quadrature (fe1.get_unit_face_support_points()); + + // now create the object with which we will generate the normal vectors + FEFaceValues fe_face_values + (fe1, quadrature, update_q_points | update_normal_vectors); + + // have a vector that stores the location of each vector-dof tuple + // we want to rotate. + typedef std_cxx1x::array DoFTuple; + + FullMatrix rot_transformation + =IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face); + + for (unsigned int i=0; i::const_iterator comp_it + = std::find (first_vector_components.begin(), + first_vector_components.end(), + fe1.face_system_to_component_index(i).first); + if (comp_it != first_vector_components.end()) { - set_periodicity_constraints(face_1, face_2, - matrix, - constraint_matrix, - component_mask, - face_orientation, face_flip, face_rotation); + const unsigned int first_vector_component = *comp_it; + + // find corresponding other components of vector + DoFTuple vector_dofs; + vector_dofs[0] = i; + + Assert(*comp_it+dim<=fe1.n_components(), + ExcMessage("Error: the finite element does not have enough components " + "to define rotated periodic boundaries.")); + + for (unsigned int k=0; k= + first_vector_component) + && + (fe1.face_system_to_component_index(k).first < + first_vector_component + dim)) + vector_dofs[fe1.face_system_to_component_index(k).first - + first_vector_component] + = k; + + for(int i=0; ihas_children() == false) + { + // TODO: This is horrible... + FullMatrix inverse(transformation.size(0)); + inverse.invert(transformation); + set_periodicity_constraints(face_2, face_1, + inverse, + constraint_matrix, + component_mask, + face_orientation, face_flip, face_rotation); + } + else + { + set_periodicity_constraints(face_1, face_2, + transformation, + constraint_matrix, + component_mask, + face_orientation, face_flip, face_rotation); + } + } } diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index a67a64cee0..6255822b2a 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -2211,26 +2211,33 @@ next_cell: * An orthogonal equality test for points: * * point1 and point2 are considered equal, if - * (point1 + offset) - point2 + * matrix.(point1 + offset) - point2 * is parallel to the unit vector in */ template inline bool orthogonal_equality (const dealii::Point &point1, const dealii::Point &point2, const int direction, - const dealii::Tensor<1,spacedim> &offset) + const dealii::Tensor<1,spacedim> &offset, + const FullMatrix &matrix) { Assert (0<=direction && direction 1.e-10) - return false; - } + double transformed_p1_comp=0.; + + for (int j = 0; j < spacedim; ++j) + transformed_p1_comp += matrix(i,j)*point1(j) + offset[j]; + + if (fabs(transformed_p1_comp - point2(i)) > 1.e-10) + return false; + } + return true; } @@ -2319,7 +2326,8 @@ next_cell: const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix) { static const int dim = FaceIterator::AccessorType::dimension; @@ -2339,7 +2347,7 @@ next_cell: it++) { if (orthogonal_equality(face1->vertex(i),face2->vertex(*it), - direction, offset)) + direction, offset, matrix)) { matching[i] = *it; face2_vertices.erase(it); @@ -2362,11 +2370,12 @@ next_cell: orthogonal_equality (const FaceIterator &face1, const FaceIterator &face2, const int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset, + const FullMatrix &matrix) { // Call the function above with a dummy orientation array std::bitset<3> dummy; - return orthogonal_equality (dummy, face1, face2, direction, offset); + return orthogonal_equality (dummy, face1, face2, direction, offset, matrix); } @@ -2409,7 +2418,8 @@ next_cell: if (GridTools::orthogonal_equality(orientation, cell1->face(face_idx1), cell2->face(face_idx2), - direction, offset)) + direction, offset, + matrix)) { // We have a match, so insert the matching pairs and // remove the matched cell in pairs2 to speed up the diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index d1c02d06ad..b28705ca04 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -291,26 +291,30 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const X::active_face_iterator&, const X::active_face_iterator&, const int, - const Tensor<1,deal_II_space_dimension> &); + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); template bool orthogonal_equality (std::bitset<3> &, const X::face_iterator&, const X::face_iterator&, const int, - const Tensor<1,deal_II_space_dimension> &); + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); template bool orthogonal_equality (const X::active_face_iterator&, const X::active_face_iterator&, const int, - const Tensor<1,deal_II_space_dimension> &); + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); template bool orthogonal_equality (const X::face_iterator&, const X::face_iterator&, const int, - const Tensor<1,deal_II_space_dimension> &); + const Tensor<1,deal_II_space_dimension> &, + const FullMatrix &); #if deal_II_dimension >= 2