From df31dc6e41f2075691148eccc42f538fc1f705cd Mon Sep 17 00:00:00 2001 From: Daniel Garcia-Sanchez Date: Thu, 31 Oct 2019 21:39:00 +0100 Subject: [PATCH] Add Bloch periodic conditions --- include/deal.II/dofs/dof_tools.h | 28 +++++++-- source/dofs/dof_tools_constraints.cc | 71 ++++++++++++++++++----- source/dofs/dof_tools_constraints.inst.in | 24 +++++--- 3 files changed, 94 insertions(+), 29 deletions(-) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index f04bb8d077..79cf3de253 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -1150,10 +1150,24 @@ namespace DoFTools * This function makes sure that identity constraints don't create cycles * in @p constraints. * + * @p periodicity_factor can be used to to implement Bloch periodic conditions + * (a.k.a. phase shift periodic conditions) of the form + * $\psi(\mathbf{r})=e^{-i\mathbf{k}\cdot\mathbf{r}}u(\mathbf{r})$ + * where $u$ is periodic with the same periodicity as the crystal lattice and + * and $\mathbf{k}$ is the wavevector, see + * [https://en.wikipedia.org/wiki/Bloch_wave](https://en.wikipedia.org/wiki/Bloch_wave). + * The solution at @p face_2 is equal to the solution at @p face_1 times + * @p periodicity_factor. For example, if the solution at @p face_1 is + * $\psi(0)$ and $\mathbf{d} is the corresponding point on @p face_2, then + * the solution at @p face_2 should be + * $\psi(d) = \psi(0)e^{-i \mathbf{k}\cdot \mathbf{d}}$. This condition can be + * implemented using + * $\mathrm{periodicity\_factor}=e^{-i \mathbf{k}\cdot \mathbf{d}}$. + * * Detailed information can be found in the see * @ref GlossPeriodicConstraints "Glossary entry on periodic boundary conditions". * - * @author Matthias Maier, 2012 - 2015 + * @author Matthias Maier, Daniel Garcia-Sanchez 2012 - 2019 */ template void @@ -1167,7 +1181,8 @@ namespace DoFTools const bool face_rotation = false, const FullMatrix & matrix = FullMatrix(), const std::vector &first_vector_components = - std::vector()); + std::vector(), + const number periodicity_factor = 1.); @@ -1201,7 +1216,8 @@ namespace DoFTools AffineConstraints & constraints, const ComponentMask & component_mask = ComponentMask(), const std::vector &first_vector_components = - std::vector()); + std::vector(), + const number periodicity_factor = 1.); @@ -1244,7 +1260,8 @@ namespace DoFTools const types::boundary_id b_id2, const unsigned int direction, AffineConstraints &constraints, - const ComponentMask & component_mask = ComponentMask()); + const ComponentMask & component_mask = ComponentMask(), + const number periodicity_factor = 1.); @@ -1280,7 +1297,8 @@ namespace DoFTools const types::boundary_id b_id, const unsigned int direction, AffineConstraints &constraints, - const ComponentMask & component_mask = ComponentMask()); + const ComponentMask & component_mask = ComponentMask(), + const number periodicity_factor = 1.); /** * @} diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 344b3b618c..a1a619d2d3 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1823,7 +1823,8 @@ namespace DoFTools const ComponentMask & component_mask, const bool face_orientation, const bool face_flip, - const bool face_rotation) + const bool face_rotation, + const number periodicity_factor) { static const int dim = FaceIterator::AccessorType::dimension; static const int spacedim = FaceIterator::AccessorType::space_dimension; @@ -1864,7 +1865,8 @@ namespace DoFTools component_mask, face_orientation, face_flip, - face_rotation); + face_rotation, + periodicity_factor); } return; } @@ -1978,7 +1980,7 @@ namespace DoFTools bool is_identity_constrained = false; unsigned int target = numbers::invalid_unsigned_int; - double constraint_factor = 1.; + number constraint_factor = periodicity_factor; constexpr double eps = 1.e-13; for (unsigned int jj = 0; jj < dofs_per_face; ++jj) @@ -1996,7 +1998,7 @@ namespace DoFTools } is_identity_constrained = true; target = jj; - constraint_factor = entry; + constraint_factor = entry * periodicity_factor; } } @@ -2095,7 +2097,7 @@ namespace DoFTools // In case of a dependency cycle we, either // - do nothing if cycle_constraint_factor == 1. In this case all - // degrees + // degrees // of freedom are already periodically constrained, // - otherwise, force all dofs to zero (by setting dof_left to // zero). The reasoning behind this is the fact that @@ -2114,6 +2116,30 @@ namespace DoFTools affine_constraints.add_entry(dof_left, dof_right, constraint_factor); + // The number 1e10 in the assert below is arbitrary. If the + // absolute value of constraint_factor is too large, then probably + // the absolute value of periodicity_factor is too large or too + // small. This would be equivalent to an evanescent wave that has + // a very small wavelength. A quick calculation shows that if + // |periodicity_factor| > 1e10 -> |np.exp(ikd)|> 1e10, therefore k + // is imaginary (evanescent wave) and the evanescent wavelength is + // 0.27 times smaller than the dimension of the structure, + // lambda=((2*pi)/log(1e10))*d. Imaginary wavenumbers can be + // interesting in some cases + // (https://doi.org/10.1103/PhysRevA.94.033813).In order to + // implement the case of in which the wavevector can be imaginary + // it would be necessary to rewrite this function and the dof + // ordering method should be modified. + // Let's take the following constraint a*x1 + b*x2 = 0. You could + // just always pick x1 = b/a*x2, but in practice this is not so + // stable if a could be a small number -- intended to be zero, but + // just very small due to roundoff. Of course, constraining x2 in + // terms of x1 has the same problem. So one chooses x1 = b/a*x2 if + // |b|<|a|, and x2 = a/b*x1 if |a|<|b|. + Assert( + std::abs(constraint_factor) < 1e10, + ExcMessage( + "The periodicity constraint is too large. The parameter periodicity_factor might be too large or too small.")); } } /* for dofs_per_face */ } @@ -2230,7 +2256,8 @@ namespace DoFTools const bool face_flip, const bool face_rotation, const FullMatrix & matrix, - const std::vector & first_vector_components) + const std::vector & first_vector_components, + const number periodicity_factor) { static const int dim = FaceIterator::AccessorType::dimension; static const int spacedim = FaceIterator::AccessorType::space_dimension; @@ -2368,7 +2395,8 @@ namespace DoFTools face_flip, face_rotation, matrix, - first_vector_components); + first_vector_components, + periodicity_factor); } } else @@ -2405,7 +2433,8 @@ namespace DoFTools component_mask, face_orientation, face_flip, - face_rotation); + face_rotation, + periodicity_factor); } else { @@ -2419,7 +2448,8 @@ namespace DoFTools component_mask, face_orientation, face_flip, - face_rotation); + face_rotation, + periodicity_factor); } } else @@ -2441,7 +2471,8 @@ namespace DoFTools face_orientation ? face_rotation ^ face_flip : face_flip, - face_rotation); + face_rotation, + periodicity_factor); } } } @@ -2456,7 +2487,8 @@ namespace DoFTools & periodic_faces, AffineConstraints & constraints, const ComponentMask & component_mask, - const std::vector &first_vector_components) + const std::vector &first_vector_components, + const number periodicity_factor) { // Loop over all periodic faces... for (auto &pair : periodic_faces) @@ -2480,7 +2512,8 @@ namespace DoFTools pair.orientation[1], pair.orientation[2], pair.matrix, - first_vector_components); + first_vector_components, + periodicity_factor); } } @@ -2495,7 +2528,8 @@ namespace DoFTools const types::boundary_id b_id2, const unsigned int direction, dealii::AffineConstraints &constraints, - const ComponentMask &component_mask) + const ComponentMask &component_mask, + const number periodicity_factor) { static const int space_dim = DoFHandlerType::space_dimension; (void)space_dim; @@ -2516,7 +2550,9 @@ namespace DoFTools make_periodicity_constraints(matched_faces, constraints, - component_mask); + component_mask, + std::vector(), + periodicity_factor); } @@ -2527,7 +2563,8 @@ namespace DoFTools const types::boundary_id b_id, const unsigned int direction, AffineConstraints &constraints, - const ComponentMask & component_mask) + const ComponentMask & component_mask, + const number periodicity_factor) { static const int dim = DoFHandlerType::dimension; static const int space_dim = DoFHandlerType::space_dimension; @@ -2551,7 +2588,9 @@ namespace DoFTools make_periodicity_constraints(matched_faces, constraints, - component_mask); + component_mask, + std::vector(), + periodicity_factor); } diff --git a/source/dofs/dof_tools_constraints.inst.in b/source/dofs/dof_tools_constraints.inst.in index 8b321f8948..6f6966eb1e 100644 --- a/source/dofs/dof_tools_constraints.inst.in +++ b/source/dofs/dof_tools_constraints.inst.in @@ -58,7 +58,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) bool, bool, const FullMatrix &, - const std::vector &); + const std::vector &, + const double); #ifdef DEAL_II_WITH_COMPLEX_VALUES template void DoFTools::make_periodicity_constraints( @@ -70,7 +71,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) bool, bool, const FullMatrix &, - const std::vector &); + const std::vector &, + std::complex); #endif template void @@ -79,7 +81,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) GridTools::PeriodicFacePair::cell_iterator>> &, AffineConstraints &, const ComponentMask &, - const std::vector &); + const std::vector &, + const double); #ifdef DEAL_II_WITH_COMPLEX_VALUES template void DoFTools::make_periodicity_constraints, @@ -88,7 +91,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) GridTools::PeriodicFacePair::cell_iterator>> &, AffineConstraints> &, const ComponentMask &, - const std::vector &); + const std::vector &, + const std::complex); #endif template void DoFTools::make_periodicity_constraints( @@ -97,7 +101,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const types::boundary_id, const unsigned int, AffineConstraints &, - const ComponentMask &); + const ComponentMask &, + const double); #ifdef DEAL_II_WITH_COMPLEX_VALUES template void DoFTools::make_periodicity_constraints( @@ -106,7 +111,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const types::boundary_id, const unsigned int, AffineConstraints> &, - const ComponentMask &); + const ComponentMask &, + const std::complex); #endif template void DoFTools::make_periodicity_constraints( @@ -114,7 +120,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const types::boundary_id, const unsigned int, AffineConstraints &, - const ComponentMask &); + const ComponentMask &, + const double); #ifdef DEAL_II_WITH_COMPLEX_VALUES template void DoFTools::make_periodicity_constraints( @@ -122,7 +129,8 @@ for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS) const types::boundary_id, const unsigned int, AffineConstraints> &, - const ComponentMask &); + const ComponentMask &, + const std::complex); #endif } -- 2.39.5