tensor_pols(0),
n_shape_functions(2),
renumber(0),
- use_mapping_q_on_all_cells (false)
+ use_mapping_q_on_all_cells (false),
+ feq(degree)
{}
tensor_pols(0),
n_shape_functions(2),
renumber(0),
- use_mapping_q_on_all_cells (m.use_mapping_q_on_all_cells)
+ use_mapping_q_on_all_cells (m.use_mapping_q_on_all_cells),
+ feq(degree)
{}
tensor_pols(0),
n_shape_functions(Utilities::fixed_power<dim>(degree+1)),
renumber(0),
- use_mapping_q_on_all_cells (use_mapping_q_on_all_cells)
+ use_mapping_q_on_all_cells (use_mapping_q_on_all_cells),
+ feq(degree)
{
// Construct the tensor product
// polynomials used as shape
tensor_pols(0),
n_shape_functions(mapping.n_shape_functions),
renumber(mapping.renumber),
- use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells)
+ use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells),
+ feq(degree)
{
tensor_pols=new TensorProductPolynomials<dim> (*mapping.tensor_pols);
laplace_on_quad_vector=mapping.laplace_on_quad_vector;
&straight_boundary);
boundary->get_intermediate_points_on_line (line, line_points);
- a.insert (a.end(), line_points.begin(), line_points.end());
+ if (dim==3)
+ {
+ // in 3D, lines might be in wrong
+ // orientation. if so, reverse
+ // the vector
+ if (cell->line_orientation(line_no))
+ a.insert (a.end(), line_points.begin(), line_points.end());
+ else
+ a.insert (a.end(), line_points.rbegin(), line_points.rend());
+ }
+ else
+ // in 2D, lines always have the
+ // correct orientation. simply
+ // append all points
+ a.insert (a.end(), line_points.begin(), line_points.end());
+
}
}
}
face_flip = cell->face_flip (face_no),
face_rotation = cell->face_rotation (face_no);
+#ifdef DEBUG
// some sanity checks up front
for (unsigned int i=0; i<vertices_per_face; ++i)
Assert(face->vertex_index(i)==cell->vertex_index(
Assert(face->line(i)==cell->line(GeometryInfo<3>::face_to_cell_lines(
face_no, i, face_orientation, face_flip, face_rotation)),
ExcInternalError());
+#endif
// if face at boundary, then
// ask boundary object to
{
face->get_triangulation().get_boundary(face->boundary_indicator())
.get_intermediate_points_on_quad (face, quad_points);
- a.insert (a.end(), quad_points.begin(), quad_points.end());
+ // in 3D, the orientation, flip and
+ // rotation of the face might not
+ // match what we expect here, namely
+ // the standard orientation. thus
+ // reorder points accordingly. since
+ // a Mapping uses the same shape
+ // function as an FEQ, we can ask a
+ // FEQ to do the reordering for us.
+ for (unsigned int i=0; i<quad_points.size(); ++i)
+ a.push_back(quad_points[feq.adjust_quad_dof_index_for_face_orientation(i,
+ face_orientation,
+ face_flip,
+ face_rotation)]);
}
else
{
ExcDimensionMismatch(b.size(),
vertices_per_face+lines_per_face*(degree-1)));
- // sort the points into b
+ // sort the points into b. We
+ // used access from the cell (not
+ // from the face) to fill b, so
+ // we can assume a standard face
+ // orientation. Doing so, the
+ // calculated points will be in
+ // standard orientation as well.
for (unsigned int i=0; i<vertices_per_face; ++i)
- b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i,
- face_orientation,
- face_flip,
- face_rotation)];
+ b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)];
for (unsigned int i=0; i<lines_per_face; ++i)
for (unsigned int j=0; j<degree-1; ++j)
b[vertices_per_face+i*(degree-1)+j]=
a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
- face_no, i, face_orientation, face_flip, face_rotation)*(degree-1)+j];
+ face_no, i)*(degree-1)+j];
- // Now b includes the
- // right order of
- // support points on
- // the quad to apply
- // the laplace vector
+ // Now b includes the support
+ // points on the quad and we can
+ // apply the laplace vector
apply_laplace_vector(laplace_on_quad_vector, b);
Assert(b.size()==4*degree+(degree-1)*(degree-1),
ExcDimensionMismatch(b.size(), 4*degree+(degree-1)*(degree-1)));
// from a straight
// boundary object
straight_boundary.get_intermediate_points_on_quad (face, quad_points);
- a.insert (a.end(), quad_points.begin(), quad_points.end());
+ // in 3D, the orientation, flip
+ // and rotation of the face might
+ // not match what we expect here,
+ // namely the standard
+ // orientation. thus reorder
+ // points accordingly. since a
+ // Mapping uses the same shape
+ // function as an FEQ, we can ask
+ // a FEQ to do the reordering for
+ // us.
+ for (unsigned int i=0; i<quad_points.size(); ++i)
+ a.push_back(quad_points[feq.adjust_quad_dof_index_for_face_orientation(i,
+ face_orientation,
+ face_flip,
+ face_rotation)]);
}
}
}