]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Refactor a bit, in preparation to dealing with periodic faces not on
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 May 2013 22:10:32 +0000 (22:10 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 May 2013 22:10:32 +0000 (22:10 +0000)
the same refinement level.

git-svn-id: https://svn.dealii.org/trunk@29521 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/dofs/dof_tools.cc

index 45a2e0406f9cb2ee51a3d28835af00fd222c59dc..675ba22fd0bd2197815c059fc053286cdeb5d15f 100644 (file)
@@ -2786,7 +2786,120 @@ namespace DoFTools
 
 
 
-  template<typename FaceIterator>
+  namespace
+  {
+    // enter constraints for periodicity into the given ConstraintMatrix object.
+    // this function is called when at least one of the two face iterators corresponds
+    // to an active object without further children
+    //
+    // @param transformation A matrix that maps degrees of freedom from one face
+    // to another. If the DoFs on the two faces are supposed to match exactly, then
+    // the matrix so provided will be the identity matrix. if face 2 is once refined
+    // from face 1, then the matrix needs to be the interpolation matrix from a face
+    // to this particular child
+    template <typename FaceIterator>
+    void
+    set_periodicity_constraints (const FaceIterator                          &face_1,
+                                 const typename identity<FaceIterator>::type &face_2,
+                                 const FullMatrix<double>                    &transformation,
+                                 dealii::ConstraintMatrix                    &constraint_matrix,
+                                 const ComponentMask                         &component_mask,
+                                 const bool                                   face_orientation,
+                                 const bool                                   face_flip,
+                                 const bool                                   face_rotation)
+    {
+      static const int dim      = FaceIterator::AccessorType::dimension;
+      static const int spacedim = FaceIterator::AccessorType::space_dimension;
+
+      // we should be in the case were both faces are active
+      // and have no children:
+      Assert (!face_1->has_children() && !face_2->has_children(),
+              ExcNotImplemented());
+
+      Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1,
+              ExcInternalError());
+
+      // ... then we match the corresponding DoFs of both faces ...
+      const unsigned int face_1_index = face_1->nth_active_fe_index(0);
+      const unsigned int face_2_index = face_2->nth_active_fe_index(0);
+      Assert (face_1->get_fe(face_1_index) == face_2->get_fe(face_1_index),
+              ExcMessage ("Matching periodic cells need to use the same finite element"));
+
+      const FiniteElement<dim,spacedim> &fe = face_1->get_fe(face_1_index);
+
+      Assert (component_mask.represents_n_components(fe.n_components()),
+              ExcMessage ("The number of components in the mask has to be either "
+                          "zero or equal to the number of components in the finite "
+                          "element."));
+
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+
+      std::vector<unsigned int> dofs_1(dofs_per_face);
+      std::vector<unsigned int> dofs_2(dofs_per_face);
+
+      face_1->get_dof_indices(dofs_1, face_1_index);
+      face_2->get_dof_indices(dofs_2, face_2_index);
+
+
+      // Well, this is a hack:
+      //
+      // There is no
+      //   face_to_face_index(face_index,
+      //                      face_orientation,
+      //                      face_flip,
+      //                      face_rotation)
+      // function in FiniteElementData, so we have to use
+      //   face_to_cell_index(face_index, face
+      //                      face_orientation,
+      //                      face_flip,
+      //                      face_rotation)
+      // But this will give us an index on a cell - something we cannot work
+      // with directly. But luckily we can match them back :-]
+
+      std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
+
+      // Build up a cell to face index for face_2:
+      for (unsigned int i = 0; i < dofs_per_face; ++i)
+        {
+          const unsigned int cell_index =
+            fe.face_to_cell_index(i,
+                                  0, /* It doesn't really matter, just assume
+                                      * we're on the first face...
+                                      */
+                                  true, false, false // default orientation
+                                 );
+          cell_to_rotated_face_index[cell_index] = i;
+        }
+
+      for (unsigned int i = 0; i < dofs_per_face; ++i)
+        {
+          // Query the correct face_index on face_2 respecting the given
+          // orientation:
+          const unsigned int j =
+            cell_to_rotated_face_index[fe.face_to_cell_index(i,
+                                                             0, /* It doesn't really matter, just assume
+                                                                 * we're on the first face...
+                                                                 */
+                                                             face_orientation,
+                                                             face_flip,
+                                                             face_rotation)];
+
+          // And finally constrain the two DoFs respecting component_mask:
+          if (component_mask.n_selected_components(fe.n_components()) == fe.n_components()
+              || component_mask[fe.face_system_to_component_index(i).first])
+            {
+              if (!constraint_matrix.is_constrained(dofs_1[i]))
+                {
+                  constraint_matrix.add_line(dofs_1[i]);
+                  constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0);
+                }
+            }
+        }
+    }
+  }
+
+
+  template <typename FaceIterator>
   void
   make_periodicity_constraints (const FaceIterator                          &face_1,
                                 const typename identity<FaceIterator>::type &face_2,
@@ -2882,93 +2995,15 @@ namespace DoFTools
                                           face_flip,
                                           face_rotation);
           }
-        return;
-      }
-
-    // .. otherwise we should be in the case were both faces are active
-    // and have no children ..
-    Assert (!face_1->has_children() && !face_2->has_children(),
-            ExcNotImplemented());
-
-    Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1,
-            ExcInternalError());
-
-    // .. then we match the corresponding DoFs of both faces ..
-    const unsigned int face_1_index = face_1->nth_active_fe_index(0);
-    const unsigned int face_2_index = face_2->nth_active_fe_index(0);
-    Assert (face_1->get_fe(face_1_index) == face_2->get_fe(face_1_index),
-            ExcMessage ("Matching periodic cells need to use the same finite element"));
-
-    const dealii::FiniteElement<dim> &fe = face_1->get_fe(face_1_index);
-
-    Assert (component_mask.represents_n_components(fe.n_components()),
-            ExcMessage ("The number of components in the mask has to be either "
-                        "zero or equal to the number of components in the finite "
-                        "element."));
-
-    const unsigned int dofs_per_face = fe.dofs_per_face;
-
-    std::vector<unsigned int> dofs_1(dofs_per_face);
-    std::vector<unsigned int> dofs_2(dofs_per_face);
-
-    face_1->get_dof_indices(dofs_1, face_1_index);
-    face_2->get_dof_indices(dofs_2, face_2_index);
-
-
-    // Well, this is a hack:
-    //
-    // There is no
-    //   face_to_face_index(face_index,
-    //                      face_orientation,
-    //                      face_flip,
-    //                      face_rotation)
-    // function in FiniteElementData, so we have to use
-    //   face_to_cell_index(face_index, face
-    //                      face_orientation,
-    //                      face_flip,
-    //                      face_rotation)
-    // But this will give us an index on a cell - something we cannot work
-    // with directly. But luckily we can match them back :-]
-
-    std::map<unsigned int, unsigned int> cell_to_rotated_face_index;
-
-    // Build up a cell to face index for face_2:
-    for (unsigned int i = 0; i < dofs_per_face; ++i)
-      {
-        const unsigned int cell_index =
-          fe.face_to_cell_index(i,
-                                0, /* It doesn't really matter, just assume
-                                    * we're on the first face...
-                                    */
-                                true, false, false // default orientation
-                               );
-        cell_to_rotated_face_index[cell_index] = i;
-      }
-
-    for (unsigned int i = 0; i < dofs_per_face; ++i)
-      {
-        // Query the correct face_index on face_2 respecting the given
-        // orientation:
-        const unsigned int j =
-          cell_to_rotated_face_index[fe.face_to_cell_index(i,
-                                                           0, /* It doesn't really matter, just assume
-                                                               * we're on the first face...
-                                                               */
-                                                           face_orientation,
-                                                           face_flip,
-                                                           face_rotation)];
-
-        // And finally constrain the two DoFs respecting component_mask:
-        if (component_mask.n_selected_components(fe.n_components()) == fe.n_components()
-            || component_mask[fe.face_system_to_component_index(i).first])
-          {
-            if (!constraint_matrix.is_constrained(dofs_1[i]))
-              {
-                constraint_matrix.add_line(dofs_1[i]);
-                constraint_matrix.add_entry(dofs_1[i], dofs_2[j], 1.0);
-              }
-          }
       }
+    else
+      // otherwise at least one of the two faces is active and
+      // we need to enter the constraints
+      set_periodicity_constraints(face_1, face_2,
+                                 FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)),
+                                 constraint_matrix,
+                                 component_mask,
+                                  face_orientation, face_flip, face_rotation);
   }
 
 

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.