]> https://gitweb.dealii.org/ - dealii.git/commitdiff
DoFTools: restructure set_periodicity_constraints
authorMatthias Maier <tamiko@43-1.org>
Fri, 18 Jan 2019 17:14:09 +0000 (11:14 -0600)
committerMatthias Maier <tamiko@43-1.org>
Fri, 18 Jan 2019 19:37:58 +0000 (13:37 -0600)
This commit restructures the function set_periodicity_constraints. Over
the last years (and after fixing about a dozen periodicity related bugs)
this function has gotten out of control. The aim for this restructuring
is to put back some order into the function. This is achieved by:

 - Merging loops and if clauses that were redundant.

 - Minor code cleanup.

 - Switching to a different idiom:

   if(... special case ...)
     {
       // do special case
       continue;
     }

   // next special case

   This saves about 5 (!) levels of indentation and makes the control
   flow significantly more obvious to (a) comment on and (b) understand
   again later.

This commit does not introduce a functional change.

source/dofs/dof_tools_constraints.cc

index 6550f7adc32022bbe9a1972da07ee146a11736a1..5ad8b268970c7b782a3155cd75ad8f55e537490a 100644 (file)
@@ -1825,28 +1825,29 @@ namespace DoFTools
 
       Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
 
-      // if face_2 does have children, then we need to iterate over them
+      // If face_2 does have children, then we need to iterate over these
+      // children and set periodic constraints in the inverse direction:
+
       if (face_2->has_children())
         {
           Assert(face_2->n_children() ==
                    GeometryInfo<dim>::max_children_per_face,
                  ExcNotImplemented());
+
           const unsigned int dofs_per_face =
             face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face;
           FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
-          FullMatrix<double> subface_interpolation(dofs_per_face,
-                                                   dofs_per_face);
+          FullMatrix<double> subface_interp(dofs_per_face, dofs_per_face);
+
           for (unsigned int c = 0; c < face_2->n_children(); ++c)
             {
               // get the interpolation matrix recursively from the one that
               // interpolated from face_1 to face_2 by multiplying from the left
               // with the one that interpolates from face_2 to its child
-              face_1->get_fe(face_1->nth_active_fe_index(0))
-                .get_subface_interpolation_matrix(
-                  face_1->get_fe(face_1->nth_active_fe_index(0)),
-                  c,
-                  subface_interpolation);
-              subface_interpolation.mmult(child_transformation, transformation);
+              const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
+              fe.get_subface_interpolation_matrix(fe, c, subface_interp);
+              subface_interp.mmult(child_transformation, transformation);
+
               set_periodicity_constraints(face_1,
                                           face_2->child(c),
                                           child_transformation,
@@ -1856,362 +1857,362 @@ namespace DoFTools
                                           face_flip,
                                           face_rotation);
             }
+          return;
         }
-      else
-        // both faces are active. we need to match the corresponding DoFs of
-        // both faces
+
+      //
+      // If we reached this point then both faces are active. Now all
+      // that is left is to 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_2_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<types::global_dof_index> dofs_1(dofs_per_face);
+      std::vector<types::global_dof_index> 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);
+
+      // If either of the two faces has an invalid dof index, stop. This is
+      // so that there is no attempt to match artificial cells of parallel
+      // distributed triangulations.
+      //
+      // While it seems like we ought to be able to avoid even calling
+      // set_periodicity_constraints for artificial faces, this situation
+      // can arise when a face that is being made periodic is only
+      // partially touched by the local subdomain.
+      // make_periodicity_constraints will be called recursively even for
+      // the section of the face that is not touched by the local
+      // subdomain.
+      //
+      // Until there is a better way to determine if the cells that
+      // neighbor a face are artificial, we simply test to see if the face
+      // does not have a valid dof initialization.
+
+      for (unsigned int i = 0; i < dofs_per_face; i++)
+        if (dofs_1[i] == numbers::invalid_dof_index ||
+            dofs_2[i] == numbers::invalid_dof_index)
+          {
+            return;
+          }
+
+      // 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 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_2_index),
-            ExcMessage(
-              "Matching periodic cells need to use the same finite element"));
+          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;
+        }
 
-          const FiniteElement<dim, spacedim> &fe = face_1->get_fe(face_1_index);
+      //
+      // Loop over all dofs on face 2 and constrain them against all
+      // matching dofs on face 1:
+      //
 
-          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."));
+      for (unsigned int i = 0; i < dofs_per_face; ++i)
+        {
+          // Obey the component mask
+          if ((component_mask.n_selected_components(fe.n_components()) !=
+               fe.n_components()) &&
+              !component_mask[fe.face_system_to_component_index(i).first])
+            continue;
 
-          const unsigned int dofs_per_face = fe.dofs_per_face;
+          // As mentioned in the comment above, we need to be careful about
+          // treating identity constraints differently. consequently, find
+          // out whether this dof 'i' will be identity constrained
+          //
+          // To check whether this is the case, first see whether there
+          // are any weights other than 0 and 1, then in a first stage
+          // make sure that if so there is only one weight equal to 1,
+          //
+          // afterwards do the same for constraints of type dof1=-dof2.
 
-          std::vector<types::global_dof_index> dofs_1(dofs_per_face);
-          std::vector<types::global_dof_index> dofs_2(dofs_per_face);
+          // FIXME: Refactor into one go.
 
-          face_1->get_dof_indices(dofs_1, face_1_index);
-          face_2->get_dof_indices(dofs_2, face_2_index);
+          bool         is_identity_constrained = true;
+          const double eps                     = 1.e-13;
+          for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
+            if (std::abs(transformation(i, jj)) > eps &&
+                std::abs(transformation(i, jj) - 1.) > eps)
+              {
+                is_identity_constrained = false;
+                break;
+              }
 
-          for (unsigned int i = 0; i < dofs_per_face; i++)
+          unsigned int identity_constraint_target =
+            numbers::invalid_unsigned_int;
+          if (is_identity_constrained == true)
             {
-              if (dofs_1[i] == numbers::invalid_dof_index ||
-                  dofs_2[i] == numbers::invalid_dof_index)
+              bool one_identity_found = false;
+              for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
+                if (std::abs(transformation(i, jj) - 1.) < eps)
+                  {
+                    if (one_identity_found == false)
+                      {
+                        one_identity_found         = true;
+                        identity_constraint_target = jj;
+                      }
+                    else
+                      {
+                        is_identity_constrained = false;
+                        identity_constraint_target =
+                          numbers::invalid_unsigned_int;
+                        break;
+                      }
+                  }
+            }
+
+          bool         is_inverse_constrained = !is_identity_constrained;
+          unsigned int inverse_constraint_target =
+            numbers::invalid_unsigned_int;
+          if (is_inverse_constrained)
+            for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
+              if (std::abs(transformation(i, jj)) > eps &&
+                  std::abs(transformation(i, jj) + 1.) > eps)
                 {
-                  /* If either of these faces have no indices, stop.  This is so
-                   * that there is no attempt to match artificial cells of
-                   * parallel distributed triangulations.
-                   *
-                   * While it seems like we ought to be able to avoid even
-                   * calling set_periodicity_constraints for artificial faces,
-                   * this situation can arise when a face that is being made
-                   * periodic is only partially touched by the local subdomain.
-                   * make_periodicity_constraints will be called recursively
-                   * even for the section of the face that is not touched by the
-                   * local subdomain.
-                   *
-                   * Until there is a better way to determine if the cells that
-                   * neighbor a face are artificial, we simply test to see if
-                   * the face does not have a valid dof initialization.
-                   */
-                  return;
+                  is_inverse_constrained = false;
+                  break;
                 }
+          if (is_inverse_constrained)
+            {
+              bool one_identity_found = false;
+              for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
+                if (std::abs(transformation(i, jj) + 1) < eps)
+                  {
+                    if (one_identity_found == false)
+                      {
+                        one_identity_found        = true;
+                        inverse_constraint_target = jj;
+                      }
+                    else
+                      {
+                        is_inverse_constrained = false;
+                        inverse_constraint_target =
+                          numbers::invalid_unsigned_int;
+                        break;
+                      }
+                  }
             }
 
-          // Well, this is a hack:
+          const unsigned int target = is_identity_constrained ?
+                                        identity_constraint_target :
+                                        inverse_constraint_target;
+
+          // Fix up a rare corner case:
           //
-          // 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)
+          // Find out whether the current dof (living on face 2) also
+          // exists on face 1. If this is the case then both faces share
+          // the same dof and we are in one of two situations:
+          //  - We are about to enter an identity constraint of the dof to
+          //    itself. In this case simply do nothing.
+          //  - Otherwise, we force the dof to zero.
+          {
+            bool continue_with_next_dof = false;
+            for (unsigned int j = 0; j < dofs_per_face; ++j)
+              if (dofs_2[i] == dofs_1[j])
+                {
+                  if (!(is_identity_constrained && target == i))
+                    affine_constraints.add_line(dofs_2[i]);
+                  continue_with_next_dof = true;
+                }
+
+            if (continue_with_next_dof)
+              continue;
+          }
+
+          // Next, we work on all constraints that are not identity
+          // constraints, i.e., constraints that involve an interpolation
+          // step that constrains the current dof (on face 2) to more than
+          // one dof on face 1.
+
+          if (!is_identity_constrained && !is_inverse_constrained)
             {
-              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;
+              // The current dof is already constrained. There is nothing
+              // left to do.
+              if (affine_constraints.is_constrained(dofs_2[i]))
+                continue;
+
+              // Enter the constraint piece by piece:
+              affine_constraints.add_line(dofs_2[i]);
+
+              for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
+                {
+                  // Query the correct face_index on face_1 respecting
+                  // the given orientation:
+                  const unsigned int j =
+                    cell_to_rotated_face_index[fe.face_to_cell_index(
+                      jj, 0, face_orientation, face_flip, face_rotation)];
+
+                  if (std::abs(transformation(i, jj)) > eps)
+                    affine_constraints.add_entry(dofs_2[i],
+                                                 dofs_1[j],
+                                                 transformation(i, jj));
+                }
+
+              // Continue with next dof.
+              continue;
             }
 
-          // loop over all dofs on face 2 and constrain them against the ones on
-          // face 1
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            if ((component_mask.n_selected_components(fe.n_components()) ==
-                 fe.n_components()) ||
-                component_mask[fe.face_system_to_component_index(i).first])
-              {
-                // as mentioned in the comment above this function, we need to
-                // be careful about treating identity constraints differently.
-                // consequently, find out whether this dof 'i' will be identity
-                // constrained
-                //
-                // to check whether this is the case, first see whether there
-                // are any weights other than 0 and 1, then in a first stage
-                // make sure that if so there is only one weight equal to 1
-                //
-                // afterwards do the same for constraints of type dof1=-dof2
-                bool         is_identity_constrained = true;
-                const double eps                     = 1.e-13;
-                for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
-                  if (std::abs(transformation(i, jj)) > eps &&
-                      std::abs(transformation(i, jj) - 1.) > eps)
-                    {
-                      is_identity_constrained = false;
-                      break;
-                    }
-                unsigned int identity_constraint_target =
-                  numbers::invalid_unsigned_int;
-                if (is_identity_constrained == true)
-                  {
-                    bool one_identity_found = false;
-                    for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
-                      if (std::abs(transformation(i, jj) - 1.) < eps)
-                        {
-                          if (one_identity_found == false)
-                            {
-                              one_identity_found         = true;
-                              identity_constraint_target = jj;
-                            }
-                          else
-                            {
-                              is_identity_constrained = false;
-                              identity_constraint_target =
-                                numbers::invalid_unsigned_int;
-                              break;
-                            }
-                        }
-                  }
+          //
+          // We are left with an equality constraint.
+          //
+
+          // Query the correct face_index on face_1 respecting the
+          // given orientation:
+          const unsigned int j =
+            cell_to_rotated_face_index[fe.face_to_cell_index(
+              target, 0, face_orientation, face_flip, face_rotation)];
 
-                bool         is_inverse_constrained = !is_identity_constrained;
-                unsigned int inverse_constraint_target =
-                  numbers::invalid_unsigned_int;
-                if (is_inverse_constrained)
-                  for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
-                    if (std::abs(transformation(i, jj)) > eps &&
-                        std::abs(transformation(i, jj) + 1.) > eps)
+          if (affine_constraints.is_constrained(dofs_2[i]))
+            {
+              // If the two aren't already identity constrained (whichever
+              // way around) or already identical (in case of rotated
+              // periodicity constraints), then enter the constraint.
+              // Otherwise, there is nothing for us to do.
+              bool enter_constraint = false;
+              // see if this would add an identity constraint
+              // cycle
+              if (!affine_constraints.is_constrained(dofs_1[j]))
+                {
+                  types::global_dof_index new_dof = dofs_2[i];
+                  while (new_dof != dofs_1[j])
+                    if (affine_constraints.is_constrained(new_dof))
                       {
-                        is_inverse_constrained = false;
+                        const std::vector<
+                          std::pair<types::global_dof_index, double>>
+                          *constraint_entries =
+                            affine_constraints.get_constraint_entries(new_dof);
+                        if (constraint_entries->size() == 1)
+                          new_dof = (*constraint_entries)[0].first;
+                        else
+                          {
+                            enter_constraint = true;
+                            break;
+                          }
+                      }
+                    else
+                      {
+                        enter_constraint = true;
                         break;
                       }
-                if (is_inverse_constrained)
-                  {
-                    bool one_identity_found = false;
-                    for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
-                      if (std::abs(transformation(i, jj) + 1) < eps)
-                        {
-                          if (one_identity_found == false)
-                            {
-                              one_identity_found        = true;
-                              inverse_constraint_target = jj;
-                            }
-                          else
-                            {
-                              is_inverse_constrained = false;
-                              inverse_constraint_target =
-                                numbers::invalid_unsigned_int;
-                              break;
-                            }
-                        }
-                  }
-
-                const unsigned int target = is_identity_constrained ?
-                                              identity_constraint_target :
-                                              inverse_constraint_target;
+                }
 
-                // find out whether this dof also exists on face 1 if this is
-                // true and the constraint is no identity constraint to itself,
-                // set it to zero
-                bool constraint_set = false;
-                for (unsigned int j = 0; j < dofs_per_face; ++j)
-                  {
-                    if (dofs_2[i] == dofs_1[j])
-                      if (!(is_identity_constrained && target == i))
+              if (enter_constraint)
+                {
+                  affine_constraints.add_line(dofs_1[j]);
+                  affine_constraints.add_entry(
+                    dofs_1[j], dofs_2[i], is_identity_constrained ? 1.0 : -1.0);
+                }
+            }
+          else
+            {
+              // if the two aren't already identity constrained
+              // (whichever way around) or already identical (in
+              // case of rotated periodicity constraints), then
+              // enter the constraint. Otherwise there is nothing
+              // for us to do
+              bool enter_constraint = false;
+              if (!affine_constraints.is_constrained(dofs_1[j]))
+                {
+                  if (dofs_2[i] != dofs_1[j])
+                    enter_constraint = true;
+                }
+              else // dofs_1[j] is constrained, is it identity or
+                   // inverse constrained?
+                {
+                  const std::vector<std::pair<types::global_dof_index, double>>
+                    *constraint_entries =
+                      affine_constraints.get_constraint_entries(dofs_1[j]);
+                  if (constraint_entries->size() == 1 &&
+                      (*constraint_entries)[0].first == dofs_2[i])
+                    {
+                      if ((is_identity_constrained &&
+                           std::abs((*constraint_entries)[0].second - 1) >
+                             eps) ||
+                          (is_inverse_constrained &&
+                           std::abs((*constraint_entries)[0].second + 1) > eps))
                         {
+                          // this pair of constraints means that
+                          // both dofs have to be constrained to
+                          // 0.
                           affine_constraints.add_line(dofs_2[i]);
-                          constraint_set = true;
                         }
-                  }
-
-                if (!constraint_set)
-                  {
-                    // now treat constraints, either as an equality constraint
-                    // or as a sequence of constraints
-                    if (is_identity_constrained || is_inverse_constrained)
-                      {
-                        // Query the correct face_index on face_1 respecting the
-                        // given orientation:
-                        const unsigned int j =
-                          cell_to_rotated_face_index[fe.face_to_cell_index(
-                            target,
-                            0, /* It doesn't really matter, just assume
-                                * we're on the first face...
-                                */
-                            face_orientation,
-                            face_flip,
-                            face_rotation)];
-
-                        if (affine_constraints.is_constrained(dofs_2[i]))
+                    }
+                  else
+                    {
+                      // see if this would add an identity
+                      // constraint cycle
+                      types::global_dof_index new_dof = dofs_1[j];
+                      while (new_dof != dofs_2[i])
+                        if (affine_constraints.is_constrained(new_dof))
                           {
-                            // if the two aren't already identity constrained
-                            // (whichever way around) or already identical (in
-                            // case of rotated periodicity constraints), then
-                            // enter the constraint. otherwise there is nothing
-                            // for us still to do
-                            bool enter_constraint = false;
-                            // see if this would add an identity constraint
-                            // cycle
-                            if (!affine_constraints.is_constrained(dofs_1[j]))
+                            const std::vector<
+                              std::pair<types::global_dof_index, double>>
+                              *constraint_entries =
+                                affine_constraints.get_constraint_entries(
+                                  new_dof);
+                            if (constraint_entries->size() == 1)
+                              new_dof = (*constraint_entries)[0].first;
+                            else
                               {
-                                types::global_dof_index new_dof = dofs_2[i];
-                                while (new_dof != dofs_1[j])
-                                  if (affine_constraints.is_constrained(
-                                        new_dof))
-                                    {
-                                      const std::vector<
-                                        std::pair<types::global_dof_index,
-                                                  double>> *constraint_entries =
-                                        affine_constraints
-                                          .get_constraint_entries(new_dof);
-                                      if (constraint_entries->size() == 1)
-                                        new_dof =
-                                          (*constraint_entries)[0].first;
-                                      else
-                                        {
-                                          enter_constraint = true;
-                                          break;
-                                        }
-                                    }
-                                  else
-                                    {
-                                      enter_constraint = true;
-                                      break;
-                                    }
-                              }
-
-                            if (enter_constraint)
-                              {
-                                affine_constraints.add_line(dofs_1[j]);
-                                affine_constraints.add_entry(
-                                  dofs_1[j],
-                                  dofs_2[i],
-                                  is_identity_constrained ? 1.0 : -1.0);
+                                enter_constraint = true;
+                                break;
                               }
                           }
                         else
                           {
-                            // if the two aren't already identity constrained
-                            // (whichever way around) or already identical (in
-                            // case of rotated periodicity constraints), then
-                            // enter the constraint. otherwise there is nothing
-                            // for us still to do
-                            bool enter_constraint = false;
-                            if (!affine_constraints.is_constrained(dofs_1[j]))
-                              {
-                                if (dofs_2[i] != dofs_1[j])
-                                  enter_constraint = true;
-                              }
-                            else // dofs_1[j] is constrained, is it identity or
-                                 // inverse constrained?
-                              {
-                                const std::vector<
-                                  std::pair<types::global_dof_index, double>>
-                                  *constraint_entries =
-                                    affine_constraints.get_constraint_entries(
-                                      dofs_1[j]);
-                                if (constraint_entries->size() == 1 &&
-                                    (*constraint_entries)[0].first == dofs_2[i])
-                                  {
-                                    if ((is_identity_constrained &&
-                                         std::abs(
-                                           (*constraint_entries)[0].second -
-                                           1) > eps) ||
-                                        (is_inverse_constrained &&
-                                         std::abs(
-                                           (*constraint_entries)[0].second +
-                                           1) > eps))
-                                      {
-                                        // this pair of constraints means that
-                                        // both dofs have to be constrained to
-                                        // 0.
-                                        affine_constraints.add_line(dofs_2[i]);
-                                      }
-                                  }
-                                else
-                                  {
-                                    // see if this would add an identity
-                                    // constraint cycle
-                                    types::global_dof_index new_dof = dofs_1[j];
-                                    while (new_dof != dofs_2[i])
-                                      if (affine_constraints.is_constrained(
-                                            new_dof))
-                                        {
-                                          const std::vector<std::pair<
-                                            types::global_dof_index,
-                                            double>> *constraint_entries =
-                                            affine_constraints
-                                              .get_constraint_entries(new_dof);
-                                          if (constraint_entries->size() == 1)
-                                            new_dof =
-                                              (*constraint_entries)[0].first;
-                                          else
-                                            {
-                                              enter_constraint = true;
-                                              break;
-                                            }
-                                        }
-                                      else
-                                        {
-                                          enter_constraint = true;
-                                          break;
-                                        }
-                                  }
-                              }
-
-                            if (enter_constraint)
-                              {
-                                affine_constraints.add_line(dofs_2[i]);
-                                affine_constraints.add_entry(
-                                  dofs_2[i],
-                                  dofs_1[j],
-                                  is_identity_constrained ? 1.0 : -1.0);
-                              }
-                          }
-                      }
-                    else if (!affine_constraints.is_constrained(dofs_2[i]))
-                      {
-                        // this is just a regular constraint. enter it piece by
-                        // piece
-                        affine_constraints.add_line(dofs_2[i]);
-                        for (unsigned int jj = 0; jj < dofs_per_face; ++jj)
-                          {
-                            // Query the correct face_index on face_1 respecting
-                            // the given orientation:
-                            const unsigned int j =
-                              cell_to_rotated_face_index[fe.face_to_cell_index(
-                                jj,
-                                0,
-                                face_orientation,
-                                face_flip,
-                                face_rotation)];
-
-                            // And finally constrain the two DoFs respecting
-                            // component_mask:
-                            if (transformation(i, jj) != 0)
-                              affine_constraints.add_entry(
-                                dofs_2[i], dofs_1[j], transformation(i, jj));
+                            enter_constraint = true;
+                            break;
                           }
-                      }
-                  }
-              }
-        }
+                    }
+                }
+
+              if (enter_constraint)
+                {
+                  affine_constraints.add_line(dofs_2[i]);
+                  affine_constraints.add_entry(
+                    dofs_2[i], dofs_1[j], is_identity_constrained ? 1.0 : -1.0);
+                }
+            }
+        } /* for dofs_per_face */
     }
 
 

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.