// ---------------------------------------------------------------------
#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/std_cxx1x/array.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/table.h>
#include <deal.II/base/template_constraints.h>
Assert(face_1->at_boundary() && face_2->at_boundary(),
ExcMessage ("Faces for periodicity constraints must be on the boundary"));
- Assert(first_vector_components.size() == 0 || matrix.size(0) == matrix.size(1) == dim,
+ Assert(first_vector_components.empty()
+ || (matrix.size(0) == (int)dim && matrix.size(1) == (int)dim),
ExcMessage ("first_vector_components is nonempty, so matrix must be a rotation "
- "matrix exactly of of size dim"));
+ "matrix exactly of size dim"));
#ifdef DEBUG
{
// In the case that both faces have children, we loop over all
// children and apply make_periodicty_constrains recursively:
if (face_1->has_children() && face_2->has_children())
- {
- Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
- face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
- ExcNotImplemented());
+ {
+ Assert(face_1->n_children() == GeometryInfo<dim>::max_children_per_face &&
+ face_2->n_children() == GeometryInfo<dim>::max_children_per_face,
+ ExcNotImplemented());
- for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- {
- // Lookup the index for the second face
- unsigned int j;
- switch (dim)
- {
- case 2:
- j = lookup_table_2d[face_flip][i];
- break;
- case 3:
- j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
- break;
- default:
- AssertThrow(false, ExcNotImplemented());
- }
+ for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
+ {
+ // Lookup the index for the second face
+ unsigned int j;
+ switch (dim)
+ {
+ case 2:
+ j = lookup_table_2d[face_flip][i];
+ break;
+ case 3:
+ j = lookup_table_3d[face_orientation][face_flip][face_rotation][i];
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
- make_periodicity_constraints (face_1->child(i),
- face_2->child(j),
- constraint_matrix,
- component_mask,
- face_orientation,
- face_flip,
- face_rotation,
- matrix,
- first_vector_components);
- }
+ make_periodicity_constraints (face_1->child(i),
+ face_2->child(j),
+ constraint_matrix,
+ component_mask,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ matrix,
+ first_vector_components);
}
+ }
else
// otherwise at least one of the two faces is active and
// we need to enter the constraints
- {
- // Build up the transformation matrix:
+ {
+ // Build up the transformation matrix:
- if (first_vector_components.size() != 0 && matrix.size(0) == 0)
- {
- // We need to create a local transformation matrix...
- AssertThrow(false, ExcNotImplemented());
+// if (first_vector_components.size() != 0 && matrix.size(0) == 0)
+// {
+// // We need to create a local transformation matrix...
+// AssertThrow(false, ExcNotImplemented());
+//
+// // FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face));
+// }
- // FullMatrix<double>(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face));
- }
- else if (face_2->has_children() == false)
- {
- // TODO: This is horrible...
- FullMatrix<double> inverse;
- inverse.invert(matrix);
- set_periodicity_constraints(face_2, face_1,
- inverse,
- constraint_matrix,
- component_mask,
- face_orientation, face_flip, face_rotation);
- }
- else
+ FullMatrix<double> transformation;
+ if (!first_vector_components.empty() && matrix.size(0) == (int)dim)
+ {
+ // The matrix describes a rotation and we have to build a
+ // transformation matrix, we assume that for a 0° rotation
+ // we would have to build the identity matrix
+
+ const FiniteElement<dim> &fe1
+ = face_1->get_fe(face_1->nth_active_fe_index(0));
+
+ Quadrature<dim-1> quadrature (fe1.get_unit_face_support_points());
+
+ // now create the object with which we will generate the normal vectors
+ FEFaceValues<dim> fe_face_values
+ (fe1, quadrature, update_q_points | update_normal_vectors);
+
+ // have a vector that stores the location of each vector-dof tuple
+ // we want to rotate.
+ typedef std_cxx1x::array<unsigned int, dim> DoFTuple;
+
+ FullMatrix<double> rot_transformation
+ =IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face);
+
+ for (unsigned int i=0; i<fe1.dofs_per_face; ++i)
+ {
+ std::vector<unsigned int>::const_iterator comp_it
+ = std::find (first_vector_components.begin(),
+ first_vector_components.end(),
+ fe1.face_system_to_component_index(i).first);
+ if (comp_it != first_vector_components.end())
{
- set_periodicity_constraints(face_1, face_2,
- matrix,
- constraint_matrix,
- component_mask,
- face_orientation, face_flip, face_rotation);
+ const unsigned int first_vector_component = *comp_it;
+
+ // find corresponding other components of vector
+ DoFTuple vector_dofs;
+ vector_dofs[0] = i;
+
+ Assert(*comp_it+dim<=fe1.n_components(),
+ ExcMessage("Error: the finite element does not have enough components "
+ "to define rotated periodic boundaries."));
+
+ for (unsigned int k=0; k<fe1.dofs_per_face; ++k)
+ if ((k != i)
+ &&
+ (quadrature.point(k) == quadrature.point(i))
+ &&
+ (fe1.face_system_to_component_index(k).first >=
+ first_vector_component)
+ &&
+ (fe1.face_system_to_component_index(k).first <
+ first_vector_component + dim))
+ vector_dofs[fe1.face_system_to_component_index(k).first -
+ first_vector_component]
+ = k;
+
+ for(int i=0; i<dim; ++i)
+ {
+ rot_transformation[vector_dofs[i]][vector_dofs[i]]=0.;
+ for(int j=0; j<dim; ++j)
+ rot_transformation[vector_dofs[i]][vector_dofs[j]]=matrix[i][j];
+ }
}
+ }
+
+ transformation = rot_transformation;
+
}
+ else
+ transformation = matrix;
+
+ if (face_2->has_children() == false)
+ {
+ // TODO: This is horrible...
+ FullMatrix<double> inverse(transformation.size(0));
+ inverse.invert(transformation);
+ set_periodicity_constraints(face_2, face_1,
+ inverse,
+ constraint_matrix,
+ component_mask,
+ face_orientation, face_flip, face_rotation);
+ }
+ else
+ {
+ set_periodicity_constraints(face_1, face_2,
+ transformation,
+ constraint_matrix,
+ component_mask,
+ face_orientation, face_flip, face_rotation);
+ }
+ }
}
* An orthogonal equality test for points:
*
* point1 and point2 are considered equal, if
- * (point1 + offset) - point2
+ * matrix.(point1 + offset) - point2
* is parallel to the unit vector in <direction>
*/
template<int spacedim>
inline bool orthogonal_equality (const dealii::Point<spacedim> &point1,
const dealii::Point<spacedim> &point2,
const int direction,
- const dealii::Tensor<1,spacedim> &offset)
+ const dealii::Tensor<1,spacedim> &offset,
+ const FullMatrix<double> &matrix)
{
Assert (0<=direction && direction<spacedim,
ExcIndexRange (direction, 0, spacedim));
for (int i = 0; i < spacedim; ++i)
- {
- // Only compare coordinate-components != direction:
- if (i == direction)
- continue;
+ {
+ // Only compare coordinate-components != direction:
+ if (i == direction)
+ continue;
- if (fabs(point1(i) + offset[i] - point2(i)) > 1.e-10)
- return false;
- }
+ double transformed_p1_comp=0.;
+
+ for (int j = 0; j < spacedim; ++j)
+ transformed_p1_comp += matrix(i,j)*point1(j) + offset[j];
+
+ if (fabs(transformed_p1_comp - point2(i)) > 1.e-10)
+ return false;
+ }
+
return true;
}
const FaceIterator &face1,
const FaceIterator &face2,
const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+ const FullMatrix<double> &matrix)
{
static const int dim = FaceIterator::AccessorType::dimension;
it++)
{
if (orthogonal_equality(face1->vertex(i),face2->vertex(*it),
- direction, offset))
+ direction, offset, matrix))
{
matching[i] = *it;
face2_vertices.erase(it);
orthogonal_equality (const FaceIterator &face1,
const FaceIterator &face2,
const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
+ const FullMatrix<double> &matrix)
{
// Call the function above with a dummy orientation array
std::bitset<3> dummy;
- return orthogonal_equality (dummy, face1, face2, direction, offset);
+ return orthogonal_equality (dummy, face1, face2, direction, offset, matrix);
}
if (GridTools::orthogonal_equality(orientation,
cell1->face(face_idx1),
cell2->face(face_idx2),
- direction, offset))
+ direction, offset,
+ matrix))
{
// We have a match, so insert the matching pairs and
// remove the matched cell in pairs2 to speed up the
const X::active_face_iterator&,
const X::active_face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::face_iterator> (std::bitset<3> &,
const X::face_iterator&,
const X::face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::active_face_iterator> (const X::active_face_iterator&,
const X::active_face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
template
bool orthogonal_equality<X::face_iterator> (const X::face_iterator&,
const X::face_iterator&,
const int,
- const Tensor<1,deal_II_space_dimension> &);
+ const Tensor<1,deal_II_space_dimension> &,
+ const FullMatrix<double> &);
#if deal_II_dimension >= 2
projects / dealii-svn.git / commitdiff