<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
*
* 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
* 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
* 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.
*
* 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
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;
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));
}
}
- }
}
}
}
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);
+ }
}
// $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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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
projects / dealii-svn.git / commitdiff