From: bangerth Date: Mon, 20 May 2013 14:08:26 +0000 (+0000) Subject: Make make_periodicity_constraints work with locally refined meshes as well. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=60c93f650d2bf83cfe868ef5c8d789e2aba261bb;p=dealii-svn.git Make make_periodicity_constraints work with locally refined meshes as well. git-svn-id: https://svn.dealii.org/trunk@29525 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index bcc47bed17..3159feed5a 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -128,6 +128,13 @@ this function.
    +
  1. 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. +
    +(Wolfgang Bangerth, 2013/05/20) +
  2. +
  3. 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 diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 765bf03568..dc30e7863a 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -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 diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index 4d55e79e4f..ef49025893 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -2827,31 +2827,39 @@ namespace DoFTools Assert (face_2->n_children() == GeometryInfo::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 child_transformation (dofs_per_face, dofs_per_face); + FullMatrix subface_interpolation (dofs_per_face, dofs_per_face); for (unsigned int c=0; cn_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 &fe = face_1->get_fe( - face_1_index); + const FiniteElement &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; ihas_children()) - set_periodicity_constraints(face_2, face_1, - FullMatrix(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(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(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(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); + } } diff --git a/tests/bits/periodicity_01.cc b/tests/bits/periodicity_01.cc index 8b5e34114f..1bba307672 100644 --- a/tests/bits/periodicity_01.cc +++ b/tests/bits/periodicity_01.cc @@ -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 index 0000000000..56dc51eec4 --- /dev/null +++ b/tests/bits/periodicity_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +template +void test () +{ + deallog << dim << "D" << std::endl; + + // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell + Triangulation triangulation; + std::vector repetitions (dim, 1); + repetitions[0] = 2; + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + Point(), + (dim == 2 ? + Point(2,1) : + Point(2,1,1))); + triangulation.begin_active()->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + + FE_Q fe(1); + DoFHandler 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 index 0000000000..a80f23950c --- /dev/null +++ b/tests/bits/periodicity_02/cmp/generic @@ -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 index 0000000000..e7a4a68138 --- /dev/null +++ b/tests/bits/periodicity_03.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +template +void test () +{ + deallog << dim << "D" << std::endl; + + // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell twice + Triangulation triangulation; + std::vector repetitions (dim, 1); + repetitions[0] = 2; + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + Point(), + (dim == 2 ? + Point(2,1) : + Point(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 fe(1); + DoFHandler 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 index 0000000000..80becf73fa --- /dev/null +++ b/tests/bits/periodicity_03/cmp/generic @@ -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 index 0000000000..fc6fc614a9 --- /dev/null +++ b/tests/bits/periodicity_04.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +template +void test () +{ + deallog << dim << "D" << std::endl; + + // create a 2x1 (or 2x1x1) mesh and refine the leftmost cell twice + Triangulation triangulation; + std::vector repetitions (dim, 1); + repetitions[0] = 2; + GridGenerator::subdivided_hyper_rectangle (triangulation, + repetitions, + Point(), + (dim == 2 ? + Point(2,1) : + Point(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 fe(FE_Q(1),2); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs (fe); + DoFRenumbering::component_wise (dof_handler); + + std::vector 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 index 0000000000..80becf73fa --- /dev/null +++ b/tests/bits/periodicity_04/cmp/generic @@ -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