]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Make make_periodicity_constraints work with locally refined meshes as well.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 May 2013 14:08:26 +0000 (14:08 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 May 2013 14:08:26 +0000 (14:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@29525 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/source/dofs/dof_tools.cc
tests/bits/periodicity_01.cc
tests/bits/periodicity_02.cc [new file with mode: 0644]
tests/bits/periodicity_02/cmp/generic [new file with mode: 0644]
tests/bits/periodicity_03.cc [new file with mode: 0644]
tests/bits/periodicity_03/cmp/generic [new file with mode: 0644]
tests/bits/periodicity_04.cc [new file with mode: 0644]
tests/bits/periodicity_04/cmp/generic [new file with mode: 0644]

index bcc47bed1708111a5e27a36139a382bbf1e4f0e8..3159feed5af22747b2b49de05ebd3874f35cd474 100644 (file)
@@ -128,6 +128,13 @@ this function.
 
 <ol>
 
+<li> Improved: DoFTools::make_periodicity_constraints now also works
+for meshes where the refinement level of the two sides of the domain
+is not the same, i.e., one side is more refined than the other.
+<br>
+(Wolfgang Bangerth, 2013/05/20)
+</li>
+
 <li> Improved: Through the fields DataOutBase::VtkFlags::time and
 DataOutBase::VtkFlags::cycle, it is now possible to encode the time and/or
 cycle within a nonlinear or other iteration in VTK and VTU files written
index 765bf035684780a5b2ddc03d7e186f387991c0fc..dc30e7863a47378cbec2fa22eb3f2b828c037e97 100644 (file)
@@ -954,12 +954,14 @@ namespace DoFTools
    *
    * Otherwise, if @p face_1 and @p face_2 are not active faces, this
    * function loops recursively over the children of @p face_1 and @p face_2.
-   *
-   * For this to work @p face_1 and @p face_2 must have the same
-   * refinement history, i.e. either @p face_1 and @p face_2 must be active
-   * faces or must be isotropically refined and have the same number of
-   * child faces that recursively obey this rule. (The anisotropic case
-   * is not yet implemented.)
+   * If only one of the two faces is active, then we recursively iterate
+   * over the children of the non-active ones and make sure that the
+   * solution function on the refined side equals that on the non-refined
+   * face in much the same way as we enforce hanging node constraints
+   * at places where differently refined cells come together. (However,
+   * unlike hanging nodes, we do not enforce the requirement that there
+   * be only a difference of one refinement level between the two sides
+   * of the domain you would like to be periodic).
    *
    * This routine only constrains DoFs that are not already constrained.
    * If this routine encounters a DoF that already is constrained (for
@@ -967,9 +969,6 @@ namespace DoFTools
    * constraint (dofs the entry is constrained to, inhomogeneities) is
    * kept and nothing happens.
    *
-   * Furthermore, no DoFs belonging to (or belonging to any descendant
-   * of) @p face_2 get constrained or get marked as being constrained.
-   *
    * The flags in the @p component_mask (see @ref GlossComponentMask)
    * denote which components of the finite element space shall be
    * constrained with periodic boundary conditions. If it is left as
@@ -1107,7 +1106,7 @@ namespace DoFTools
    * boundary with faces belonging to the second boundary with the help
    * of @p orthogonal_equality.
    *
-   * If this matching is successfull it constrains all DoFs associated
+   * If this matching is successful it constrains all DoFs associated
    * with the 'first' boundary to the respective DoFs of the 'second'
    * boundary respecting the relative orientation of the two faces.
    *
@@ -1117,9 +1116,6 @@ namespace DoFTools
    * constraint (dofs the entry is constrained to, inhomogeneities) is
    * kept and nothing happens.
    *
-   * Furthermore, no DoFs belonging to the 'second' boundary get
-   * constrained or get marked as being constrained.
-   *
    * The flags in the last parameter, @p component_mask (see @ref
    * GlossComponentMask) denote which components of the finite element space
    * shall be constrained with periodic boundary conditions. If it is left
index 4d55e79e4ffd350be7677c6450e3582f31888046..ef49025893a81a15147b820dd07aa38d6bc7010b 100644 (file)
@@ -2827,31 +2827,39 @@ namespace DoFTools
           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;
+            = 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);
           for (unsigned int c=0; c<face_2->n_children(); ++c)
             {
-//TODO: Need to get interpolation matrix right
+              // 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);
               set_periodicity_constraints(face_1, face_2->child(c),
-                  transformation,
-                  constraint_matrix, component_mask,
-                  face_orientation, face_flip, face_rotation);
+                                          child_transformation,
+                                          constraint_matrix, component_mask,
+                                          face_orientation, face_flip, face_rotation);
             }
         }
       else
-      // both faces are active. we need to match the corresponding DoFs of both faces
+        // both faces are active. we need 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_1_index),
-              ExcMessage ("Matching periodic cells need to use the same finite element"));
+          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);
+          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."));
+                 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;
 
@@ -2884,34 +2892,33 @@ namespace DoFTools
               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
-              );
+                                                                    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])
+          // loop over all dofs on face 2 and constrain them again the ones on face 1
+          for (unsigned int i=0; i<dofs_per_face; ++i)
+            if (!constraint_matrix.is_constrained(dofs_2[i]))
+              if (component_mask[fe.face_system_to_component_index(i).first])
                 {
-                  if (!constraint_matrix.is_constrained(dofs_1[i]))
+                  constraint_matrix.add_line(dofs_2[i]);
+                  for (unsigned int jj=0; jj<dofs_per_face; ++jj)
                     {
-                      constraint_matrix.add_line(dofs_1[i]);
-                      constraint_matrix.add_entry(dofs_1[i], dofs_2[j],
-                          1.0);
+                      // 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(jj, 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 (transformation(i,j) != 0)
+                        constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
+                                                    transformation(i,j));
                     }
                 }
-            }
         }
     }
   }
@@ -3017,18 +3024,20 @@ namespace DoFTools
     else
       // otherwise at least one of the two faces is active and
       // we need to enter the constraints
-      if (face_1->has_children())
-        set_periodicity_constraints(face_2, face_1,
-                                    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);
-      else
-        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);
+      {
+        if (face_2->has_children() == false)
+          set_periodicity_constraints(face_2, face_1,
+                                      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);
+        else
+          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);
+      }
   }
 
 
index 8b5e34114fca37b431c0dab91c26d7e5e8515fd3..1bba307672509e6212ee26389b5cf78f6278ec33 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$ 
 //
-//    Copyright (C) 2002, 2003, 2004, 2005, 2010, 2013 by the deal.II authors and Anna Schneebeli
+//    Copyright (C) 2002, 2003, 2004, 2005, 2010, 2013 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
diff --git a/tests/bits/periodicity_02.cc b/tests/bits/periodicity_02.cc
new file mode 100644 (file)
index 0000000..56dc51e
--- /dev/null
@@ -0,0 +1,78 @@
+//----------------------------  periodicity_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2002, 2003, 2004, 2005, 2010, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  periodicity_02.cc  ---------------------------
+
+
+// check periodic boundary conditions for a simple enough case where we know
+// the exact set of constraints
+//
+// this test simply uses two hypercubes, refines one of them and matches the
+// faces at the far ends
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <iomanip>
+#include <fstream>
+
+
+
+template <int dim>
+void test () 
+{
+  deallog << dim << "D" << std::endl;
+  
+  // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell
+  Triangulation<dim> triangulation;
+  std::vector<unsigned int> repetitions (dim, 1);
+  repetitions[0] = 2;
+  GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                            repetitions,
+                                            Point<dim>(),
+                                            (dim == 2 ?
+                                             Point<dim>(2,1) :
+                                             Point<dim>(2,1,1)));
+  triangulation.begin_active()->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+  
+  FE_Q<dim>          fe(1);
+  DoFHandler<dim>    dof_handler (triangulation);
+  dof_handler.distribute_dofs (fe);
+
+  ConstraintMatrix cm;
+  DoFTools::make_periodicity_constraints (dof_handler.begin(0)->face(0),
+                                         (++dof_handler.begin(0))->face(1),
+                                         cm);
+  cm.print (deallog.get_file_stream());
+}
+
+    
+
+int main () 
+{
+  std::ofstream logfile("periodicity_02/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+  return 0;
+}
diff --git a/tests/bits/periodicity_02/cmp/generic b/tests/bits/periodicity_02/cmp/generic
new file mode 100644 (file)
index 0000000..a80f239
--- /dev/null
@@ -0,0 +1,23 @@
+
+DEAL::2D
+    4 1:  1.00000
+    6 1:  0.500000
+    6 3:  0.500000
+    9 3:  1.00000
+DEAL::3D
+    8 1:  1.00000
+    10 1:  0.500000
+    10 3:  0.500000
+    12 1:  0.500000
+    12 5:  0.500000
+    14 1:  0.250000
+    14 3:  0.250000
+    14 5:  0.250000
+    14 7:  0.250000
+    19 3:  1.00000
+    21 3:  0.500000
+    21 7:  0.500000
+    24 5:  1.00000
+    26 5:  0.500000
+    26 7:  0.500000
+    29 7:  1.00000
diff --git a/tests/bits/periodicity_03.cc b/tests/bits/periodicity_03.cc
new file mode 100644 (file)
index 0000000..e7a4a68
--- /dev/null
@@ -0,0 +1,80 @@
+//----------------------------  periodicity_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2002, 2003, 2004, 2005, 2010, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  periodicity_03.cc  ---------------------------
+
+
+// check periodic boundary conditions for a simple enough case where we know
+// the exact set of constraints
+//
+// this test simply uses two hypercubes, refines one of them twice and matches
+// the faces at the far ends. this requires recursing into children
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <iomanip>
+#include <fstream>
+
+
+
+template <int dim>
+void test () 
+{
+  deallog << dim << "D" << std::endl;
+  
+  // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell twice
+  Triangulation<dim> triangulation;
+  std::vector<unsigned int> repetitions (dim, 1);
+  repetitions[0] = 2;
+  GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                            repetitions,
+                                            Point<dim>(),
+                                            (dim == 2 ?
+                                             Point<dim>(2,1) :
+                                             Point<dim>(2,1,1)));
+  triangulation.begin_active()->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+  triangulation.begin_active(1)->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+  
+  FE_Q<dim>          fe(1);
+  DoFHandler<dim>    dof_handler (triangulation);
+  dof_handler.distribute_dofs (fe);
+
+  ConstraintMatrix cm;
+  DoFTools::make_periodicity_constraints (dof_handler.begin(0)->face(0),
+                                         (++dof_handler.begin(0))->face(1),
+                                         cm);
+  cm.print (deallog.get_file_stream());
+}
+
+    
+
+int main () 
+{
+  std::ofstream logfile("periodicity_03/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+  return 0;
+}
diff --git a/tests/bits/periodicity_03/cmp/generic b/tests/bits/periodicity_03/cmp/generic
new file mode 100644 (file)
index 0000000..80becf7
--- /dev/null
@@ -0,0 +1,41 @@
+
+DEAL::2D
+    10 1:  1.00000
+    12 1:  0.750000
+    12 3:  0.250000
+    7 1:  0.500000
+    7 3:  0.500000
+    8 3:  1.00000
+DEAL::3D
+    30 1:  1.00000
+    32 1:  0.750000
+    32 3:  0.250000
+    34 1:  0.750000
+    34 5:  0.250000
+    36 1:  0.562500
+    36 3:  0.187500
+    36 5:  0.187500
+    36 7:  0.0625000
+    15 1:  0.500000
+    15 3:  0.500000
+    42 1:  0.375000
+    42 3:  0.375000
+    42 5:  0.125000
+    42 7:  0.125000
+    22 1:  0.500000
+    22 5:  0.500000
+    46 1:  0.375000
+    46 3:  0.125000
+    46 5:  0.375000
+    46 7:  0.125000
+    18 1:  0.250000
+    18 3:  0.250000
+    18 5:  0.250000
+    18 7:  0.250000
+    16 3:  1.00000
+    19 3:  0.500000
+    19 7:  0.500000
+    23 5:  1.00000
+    25 5:  0.500000
+    25 7:  0.500000
+    28 7:  1.00000
diff --git a/tests/bits/periodicity_04.cc b/tests/bits/periodicity_04.cc
new file mode 100644 (file)
index 0000000..fc6fc61
--- /dev/null
@@ -0,0 +1,88 @@
+//----------------------------  periodicity_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2002, 2003, 2004, 2005, 2010, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  periodicity_04.cc  ---------------------------
+
+
+// check periodic boundary conditions for a simple enough case where we know
+// the exact set of constraints
+//
+// this test simply uses two hypercubes, refines one of them twice and matches
+// the faces at the far ends. this requires recursing into children
+//
+// compared to the _03 test, we also set a component mask
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <iomanip>
+#include <fstream>
+
+
+
+template <int dim>
+void test () 
+{
+  deallog << dim << "D" << std::endl;
+  
+  // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell twice
+  Triangulation<dim> triangulation;
+  std::vector<unsigned int> repetitions (dim, 1);
+  repetitions[0] = 2;
+  GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                            repetitions,
+                                            Point<dim>(),
+                                            (dim == 2 ?
+                                             Point<dim>(2,1) :
+                                             Point<dim>(2,1,1)));
+  triangulation.begin_active()->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+  triangulation.begin_active(1)->set_refine_flag ();
+  triangulation.execute_coarsening_and_refinement ();
+
+  FESystem<dim> fe(FE_Q<dim>(1),2);
+  DoFHandler<dim>    dof_handler (triangulation);
+  dof_handler.distribute_dofs (fe);
+  DoFRenumbering::component_wise (dof_handler);
+
+  std::vector<bool> mask(2, true);
+  mask[1] = false;
+  ConstraintMatrix cm;
+  DoFTools::make_periodicity_constraints (dof_handler.begin(0)->face(0),
+                                         (++dof_handler.begin(0))->face(1),
+                                         cm,
+                                         mask);
+  cm.print (deallog.get_file_stream());
+}
+
+    
+
+int main () 
+{
+  std::ofstream logfile("periodicity_04/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<2>();
+  test<3>();
+  return 0;
+}
diff --git a/tests/bits/periodicity_04/cmp/generic b/tests/bits/periodicity_04/cmp/generic
new file mode 100644 (file)
index 0000000..80becf7
--- /dev/null
@@ -0,0 +1,41 @@
+
+DEAL::2D
+    10 1:  1.00000
+    12 1:  0.750000
+    12 3:  0.250000
+    7 1:  0.500000
+    7 3:  0.500000
+    8 3:  1.00000
+DEAL::3D
+    30 1:  1.00000
+    32 1:  0.750000
+    32 3:  0.250000
+    34 1:  0.750000
+    34 5:  0.250000
+    36 1:  0.562500
+    36 3:  0.187500
+    36 5:  0.187500
+    36 7:  0.0625000
+    15 1:  0.500000
+    15 3:  0.500000
+    42 1:  0.375000
+    42 3:  0.375000
+    42 5:  0.125000
+    42 7:  0.125000
+    22 1:  0.500000
+    22 5:  0.500000
+    46 1:  0.375000
+    46 3:  0.125000
+    46 5:  0.375000
+    46 7:  0.125000
+    18 1:  0.250000
+    18 3:  0.250000
+    18 5:  0.250000
+    18 7:  0.250000
+    16 3:  1.00000
+    19 3:  0.500000
+    19 7:  0.500000
+    23 5:  1.00000
+    25 5:  0.500000
+    25 7:  0.500000
+    28 7:  1.00000

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.