* 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 <typename FaceIterator, typename number>
void
const bool face_rotation = false,
const FullMatrix<double> & matrix = FullMatrix<double>(),
const std::vector<unsigned int> &first_vector_components =
- std::vector<unsigned int>());
+ std::vector<unsigned int>(),
+ const number periodicity_factor = 1.);
AffineConstraints<number> & constraints,
const ComponentMask & component_mask = ComponentMask(),
const std::vector<unsigned int> &first_vector_components =
- std::vector<unsigned int>());
+ std::vector<unsigned int>(),
+ const number periodicity_factor = 1.);
const types::boundary_id b_id2,
const unsigned int direction,
AffineConstraints<number> &constraints,
- const ComponentMask & component_mask = ComponentMask());
+ const ComponentMask & component_mask = ComponentMask(),
+ const number periodicity_factor = 1.);
const types::boundary_id b_id,
const unsigned int direction,
AffineConstraints<number> &constraints,
- const ComponentMask & component_mask = ComponentMask());
+ const ComponentMask & component_mask = ComponentMask(),
+ const number periodicity_factor = 1.);
/**
* @}
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;
component_mask,
face_orientation,
face_flip,
- face_rotation);
+ face_rotation,
+ periodicity_factor);
}
return;
}
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)
}
is_identity_constrained = true;
target = jj;
- constraint_factor = entry;
+ constraint_factor = entry * periodicity_factor;
}
}
// 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
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 */
}
const bool face_flip,
const bool face_rotation,
const FullMatrix<double> & matrix,
- const std::vector<unsigned int> & first_vector_components)
+ const std::vector<unsigned int> & first_vector_components,
+ const number periodicity_factor)
{
static const int dim = FaceIterator::AccessorType::dimension;
static const int spacedim = FaceIterator::AccessorType::space_dimension;
face_flip,
face_rotation,
matrix,
- first_vector_components);
+ first_vector_components,
+ periodicity_factor);
}
}
else
component_mask,
face_orientation,
face_flip,
- face_rotation);
+ face_rotation,
+ periodicity_factor);
}
else
{
component_mask,
face_orientation,
face_flip,
- face_rotation);
+ face_rotation,
+ periodicity_factor);
}
}
else
face_orientation ?
face_rotation ^ face_flip :
face_flip,
- face_rotation);
+ face_rotation,
+ periodicity_factor);
}
}
}
& periodic_faces,
AffineConstraints<number> & constraints,
const ComponentMask & component_mask,
- const std::vector<unsigned int> &first_vector_components)
+ const std::vector<unsigned int> &first_vector_components,
+ const number periodicity_factor)
{
// Loop over all periodic faces...
for (auto &pair : periodic_faces)
pair.orientation[1],
pair.orientation[2],
pair.matrix,
- first_vector_components);
+ first_vector_components,
+ periodicity_factor);
}
}
const types::boundary_id b_id2,
const unsigned int direction,
dealii::AffineConstraints<number> &constraints,
- const ComponentMask &component_mask)
+ const ComponentMask &component_mask,
+ const number periodicity_factor)
{
static const int space_dim = DoFHandlerType::space_dimension;
(void)space_dim;
make_periodicity_constraints<DoFHandlerType>(matched_faces,
constraints,
- component_mask);
+ component_mask,
+ std::vector<unsigned int>(),
+ periodicity_factor);
}
const types::boundary_id b_id,
const unsigned int direction,
AffineConstraints<number> &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;
make_periodicity_constraints<DoFHandlerType>(matched_faces,
constraints,
- component_mask);
+ component_mask,
+ std::vector<unsigned int>(),
+ periodicity_factor);
}
bool,
bool,
const FullMatrix<double> &,
- const std::vector<unsigned int> &);
+ const std::vector<unsigned int> &,
+ const double);
#ifdef DEAL_II_WITH_COMPLEX_VALUES
template void DoFTools::make_periodicity_constraints(
bool,
bool,
const FullMatrix<double> &,
- const std::vector<unsigned int> &);
+ const std::vector<unsigned int> &,
+ std::complex<double>);
#endif
template void
GridTools::PeriodicFacePair<DH<deal_II_dimension>::cell_iterator>> &,
AffineConstraints<double> &,
const ComponentMask &,
- const std::vector<unsigned int> &);
+ const std::vector<unsigned int> &,
+ const double);
#ifdef DEAL_II_WITH_COMPLEX_VALUES
template void DoFTools::make_periodicity_constraints<DH<deal_II_dimension>,
GridTools::PeriodicFacePair<DH<deal_II_dimension>::cell_iterator>> &,
AffineConstraints<std::complex<double>> &,
const ComponentMask &,
- const std::vector<unsigned int> &);
+ const std::vector<unsigned int> &,
+ const std::complex<double>);
#endif
template void DoFTools::make_periodicity_constraints(
const types::boundary_id,
const unsigned int,
AffineConstraints<double> &,
- const ComponentMask &);
+ const ComponentMask &,
+ const double);
#ifdef DEAL_II_WITH_COMPLEX_VALUES
template void DoFTools::make_periodicity_constraints(
const types::boundary_id,
const unsigned int,
AffineConstraints<std::complex<double>> &,
- const ComponentMask &);
+ const ComponentMask &,
+ const std::complex<double>);
#endif
template void DoFTools::make_periodicity_constraints(
const types::boundary_id,
const unsigned int,
AffineConstraints<double> &,
- const ComponentMask &);
+ const ComponentMask &,
+ const double);
#ifdef DEAL_II_WITH_COMPLEX_VALUES
template void DoFTools::make_periodicity_constraints(
const types::boundary_id,
const unsigned int,
AffineConstraints<std::complex<double>> &,
- const ComponentMask &);
+ const ComponentMask &,
+ const std::complex<double>);
#endif
}
projects / dealii.git / commitdiff