]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
First working version for rotated periodic boundaries, tested for non-hp branch_periodic_bc
authordaniel.arndt <daniel.arndt@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jul 2014 09:52:45 +0000 (09:52 +0000)
committerdaniel.arndt <daniel.arndt@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jul 2014 09:52:45 +0000 (09:52 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_periodic_bc@33195 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/grid_tools.h
deal.II/source/dofs/dof_tools_constraints.cc
deal.II/source/grid/grid_tools.cc
deal.II/source/grid/grid_tools.inst.in

index 5844243ead70a5bb00417bbb30ea46604d3dd228..7ce0590d1818c9f9934b47b3a72d6fa1bfaf4ac4 100644 (file)
@@ -1175,7 +1175,7 @@ namespace GridTools
    const int                direction,
    std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
    const Tensor<1,CONTAINER::space_dimension> &offset = dealii::Tensor<1,CONTAINER::space_dimension>(),
-   const FullMatrix<double> &matrix = FullMatrix<double>(),
+   const FullMatrix<double> &matrix = FullMatrix<double>(IdentityMatrix(CONTAINER::space_dimension)),
    const std::vector<unsigned int> &first_vector_components = std::vector<unsigned int>());
 
 
index d045093d1527fc21442812fce1926e87ef9441fb..a2136cbd9668c9fcfafeed946aa6050a2dd20fad 100644 (file)
@@ -15,6 +15,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/multithread_info.h>
+#include <deal.II/base/std_cxx1x/array.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/base/table.h>
 #include <deal.II/base/template_constraints.h>
@@ -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<dim>::max_children_per_face &&
-               face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
-               ExcNotImplemented());
+    {
+      Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
+              face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
+              ExcNotImplemented());
 
-        for (unsigned int i = 0; i < GeometryInfo<dim>::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<dim>::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<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face));
+//           }
 
-            // FullMatrix<double>(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<double> inverse;
-            inverse.invert(matrix);
-            set_periodicity_constraints(face_2, face_1,
-                                        inverse,
-                                        constraint_matrix,
-                                        component_mask,
-                                        face_orientation, face_flip, face_rotation);
-          }
-        else
+      FullMatrix<double> 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<dim> &fe1
+          = face_1->get_fe(face_1->nth_active_fe_index(0));
+
+        Quadrature<dim-1> quadrature (fe1.get_unit_face_support_points());
+
+        // now create the object with which we will generate the normal vectors
+        FEFaceValues<dim> 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<unsigned int, dim> DoFTuple;
+        
+        FullMatrix<double> rot_transformation
+          =IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face);
+            
+        for (unsigned int i=0; i<fe1.dofs_per_face; ++i)
+        {
+          std::vector<unsigned int>::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<fe1.dofs_per_face; ++k)
+              if ((k != i)
+                  &&
+                  (quadrature.point(k) == quadrature.point(i))
+                  &&
+                  (fe1.face_system_to_component_index(k).first >=
+                    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; i<dim; ++i)
+            {
+              rot_transformation[vector_dofs[i]][vector_dofs[i]]=0.;
+              for(int j=0; j<dim; ++j)
+                rot_transformation[vector_dofs[i]][vector_dofs[j]]=matrix[i][j];
+            }
           }
+        }
+
+        transformation = rot_transformation;
+
       }
+      else
+        transformation = matrix;
+        
+      if (face_2->has_children() == false)
+      {
+        // TODO: This is horrible...
+        FullMatrix<double> 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);
+      }
+    }
   }
 
 
index a67a64cee050331abaff71e28861cb4c1008122a..6255822b2a9f80594f70b79f1ca2e96cd99b06a5 100644 (file)
@@ -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 <direction>
    */
   template<int spacedim>
   inline bool orthogonal_equality (const dealii::Point<spacedim> &point1,
                                    const dealii::Point<spacedim> &point2,
                                    const int                     direction,
-                                   const dealii::Tensor<1,spacedim> &offset)
+                                   const dealii::Tensor<1,spacedim> &offset,
+                                   const FullMatrix<double>      &matrix)
   {
     Assert (0<=direction && direction<spacedim,
             ExcIndexRange (direction, 0, spacedim));
     for (int i = 0; i < spacedim; ++i)
-      {
-        // Only compare coordinate-components != direction:
-        if (i == direction)
-          continue;
+    {
+      // Only compare coordinate-components != direction:
+      if (i == direction)
+        continue;
 
-        if (fabs(point1(i) + offset[i] - point2(i)) > 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<double> &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<double> &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
index d1c02d06adb12305bc434842f8f7a3df7085ffb8..b28705ca04f8d67d989b1d7beb11855f36f5972d 100644 (file)
@@ -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<double> &);
 
     template
     bool orthogonal_equality<X::face_iterator> (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<double> &);
 
     template
     bool orthogonal_equality<X::active_face_iterator> (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<double> &);
 
     template
     bool orthogonal_equality<X::face_iterator> (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<double> &);
 
     #if deal_II_dimension >= 2
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.