namespace internal
{
-
/**
* A class that provides all possible situations a face (in the current
* space dimension @p dim) might be subdivided into subfaces. For
/**
* Return the amount of memory occupied by an object of this type.
*/
- static std::size_t memory_consumption ();
+ static constexpr std::size_t memory_consumption ();
/**
* Exception.
RefinementCase<1>
RefinementCase<1>::cut_axis (const unsigned int i)
{
- const unsigned int dim = 1;
- Assert (i < dim, ExcIndexRange(i, 0, dim));
+ Assert (i < 1, ExcIndexRange(i, 0, 1));
- static const RefinementCase options[dim] = { RefinementPossibilities<1>::cut_x };
+ const RefinementCase options[1] = { RefinementPossibilities<1>::cut_x };
return options[i];
}
RefinementCase<2>
RefinementCase<2>::cut_axis (const unsigned int i)
{
- const unsigned int dim = 2;
- Assert (i < dim, ExcIndexRange(i, 0, dim));
+ Assert (i < 2, ExcIndexRange(i, 0, 2));
- static const RefinementCase options[dim] = { RefinementPossibilities<2>::cut_x,
- RefinementPossibilities<2>::cut_y
- };
+ const RefinementCase options[2] = { RefinementPossibilities<2>::cut_x,
+ RefinementPossibilities<2>::cut_y
+ };
return options[i];
}
RefinementCase<3>
RefinementCase<3>::cut_axis (const unsigned int i)
{
- const unsigned int dim = 3;
- Assert (i < dim, ExcIndexRange(i, 0, dim));
+ Assert (i < 3, ExcIndexRange(i, 0, 3));
- static const RefinementCase options[dim] = { RefinementPossibilities<3>::cut_x,
- RefinementPossibilities<3>::cut_y,
- RefinementPossibilities<3>::cut_z
- };
+ const RefinementCase options[3] = { RefinementPossibilities<3>::cut_x,
+ RefinementPossibilities<3>::cut_y,
+ RefinementPossibilities<3>::cut_z
+ };
return options[i];
}
template <int dim>
template <class Archive>
+inline
void RefinementCase<dim>::serialize (Archive &ar,
const unsigned int)
{
const RefinementCase<dim> /*refine_case*/)
{
- AssertThrow (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
return Point<dim>();
}
const unsigned int /*child_index*/,
const RefinementCase<dim> /*refine_case*/)
{
- AssertThrow (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
return Point<dim>();
}
return numbers::invalid_unsigned_int;
}
+template <>
+inline
+unsigned int
+GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert(vertex<GeometryInfo<3>::vertices_per_face,
+ ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
+
+ // set up a table to make sure that
+ // we handle non-standard faces correctly
+ //
+ // so set up a table that for each vertex (of
+ // a quad in standard position) describes
+ // which vertex to take
+ //
+ // first index: four vertices 0...3
+ //
+ // second index: face_orientation; 0:
+ // opposite normal, 1: standard
+ //
+ // third index: face_flip; 0: standard, 1:
+ // face rotated by 180 degrees
+ //
+ // forth index: face_rotation: 0: standard,
+ // 1: face rotated by 90 degrees
+
+ constexpr unsigned int vertex_translation[4][2][2][2] =
+ {
+ { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
+ { 3, 1 }
+ }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
+ { { 0, 2 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
+ { 3, 1 }
+ }
+ },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
+
+ { { { 2, 3 }, // vertex 1 ...
+ { 1, 0 }
+ },
+ { { 1, 0 },
+ { 2, 3 }
+ }
+ },
+
+ { { { 1, 0 }, // vertex 2 ...
+ { 2, 3 }
+ },
+ { { 2, 3 },
+ { 1, 0 }
+ }
+ },
+
+ { { { 3, 1 }, // vertex 3 ...
+ { 0, 2 }
+ },
+ { { 3, 1 },
+ { 0, 2 }
+ }
+ }
+ };
+
+ return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
+ const bool,
+ const bool,
+ const bool)
+{
+ Assert(dim>1, ExcImpossibleInDim(dim));
+ Assert(vertex<GeometryInfo<dim>::vertices_per_face,
+ ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
+ return vertex;
+}
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
+{
+ constexpr unsigned int n_children[RefinementCase<3>::cut_xyz+1]=
+ {0, 2, 2, 4, 2, 4, 4, 8};
+
+ return n_children[ref_case];
+}
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
+{
+ Assert(false, ExcImpossibleInDim(1));
+ return 0;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
+{
+ return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
+{
+ const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic+1]=
+ {0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
+ return nsubs[subface_case];
+}
+
+
+template <>
+inline
+double
+GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
+ const unsigned int)
+{
+ return 1;
+}
+
+
+template <>
+inline
+double
+GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
+ const unsigned int)
+{
+ double ratio=1;
+ switch (subface_case)
+ {
+ case internal::SubfaceCase<2>::case_none:
+ // Here, an
+ // Assert(false,ExcInternalError())
+ // would be the right
+ // choice, but
+ // unfortunately the
+ // current function is
+ // also called for faces
+ // without children (see
+ // tests/fe/mapping.cc).
+// Assert(false, ExcMessage("Face has no subfaces."));
+ // Furthermore, assign
+ // following value as
+ // otherwise the
+ // bits/volume_x tests
+ // break
+ ratio=1./GeometryInfo<2>::max_children_per_face;
+ break;
+ case internal::SubfaceCase<2>::case_x:
+ ratio=0.5;
+ break;
+ default:
+ // there should be no
+ // cases left
+ Assert(false, ExcInternalError());
+ break;
+ }
+
+ return ratio;
+}
+
+
+template <>
+inline
+double
+GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
+ const unsigned int subface_no)
+{
+ double ratio=1;
+ switch (subface_case)
+ {
+ case internal::SubfaceCase<3>::case_none:
+ // Here, an
+ // Assert(false,ExcInternalError())
+ // would be the right
+ // choice, but
+ // unfortunately the
+ // current function is
+ // also called for faces
+ // without children (see
+ // tests/bits/mesh_3d_16.cc). Add
+ // following switch to
+ // avoid diffs in
+ // tests/bits/mesh_3d_16
+ ratio=1./GeometryInfo<3>::max_children_per_face;
+ break;
+ case internal::SubfaceCase<3>::case_x:
+ case internal::SubfaceCase<3>::case_y:
+ ratio=0.5;
+ break;
+ case internal::SubfaceCase<3>::case_xy:
+ case internal::SubfaceCase<3>::case_x1y2y:
+ case internal::SubfaceCase<3>::case_y1x2x:
+ ratio=0.25;
+ break;
+ case internal::SubfaceCase<3>::case_x1y:
+ case internal::SubfaceCase<3>::case_y1x:
+ if (subface_no<2)
+ ratio=0.25;
+ else
+ ratio=0.5;
+ break;
+ case internal::SubfaceCase<3>::case_x2y:
+ case internal::SubfaceCase<3>::case_y2x:
+ if (subface_no==0)
+ ratio=0.5;
+ else
+ ratio=0.25;
+ break;
+ default:
+ // there should be no
+ // cases left
+ Assert(false, ExcInternalError());
+ break;
+ }
+
+ return ratio;
+}
+
+
+
+template <>
+RefinementCase<0>
+inline
+GeometryInfo<1>::face_refinement_case(const RefinementCase<1> &,
+ const unsigned int,
+ const bool,
+ const bool,
+ const bool)
+{
+ Assert(false, ExcImpossibleInDim(1));
+
+ return RefinementCase<0>::no_refinement;
+}
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<2>::face_refinement_case(const RefinementCase<2> &cell_refinement_case,
+ const unsigned int face_no,
+ const bool,
+ const bool,
+ const bool)
+{
+ const unsigned int dim=2;
+ Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
+ Assert(face_no<GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
+
+ const RefinementCase<dim-1>
+ ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
+ {
+ {
+ RefinementCase<dim-1>::no_refinement, // no_refinement
+ RefinementCase<dim-1>::no_refinement
+ },
+
+ {
+ RefinementCase<dim-1>::no_refinement,
+ RefinementCase<dim-1>::cut_x
+ },
+
+ {
+ RefinementCase<dim-1>::cut_x,
+ RefinementCase<dim-1>::no_refinement
+ },
+
+ {
+ RefinementCase<dim-1>::cut_x, // cut_xy
+ RefinementCase<dim-1>::cut_x
+ }
+ };
+
+ return ref_cases[cell_refinement_case][face_no/2];
+}
+
+
+template <>
+inline
+RefinementCase<2>
+GeometryInfo<3>::face_refinement_case(const RefinementCase<3> &cell_refinement_case,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool /*face_flip*/,
+ const bool face_rotation)
+{
+ const unsigned int dim=3;
+ Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
+ Assert(face_no<GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
+
+ const RefinementCase<dim-1>
+ ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
+ {
+ {
+ RefinementCase<dim-1>::no_refinement, // no_refinement
+ RefinementCase<dim-1>::no_refinement,
+ RefinementCase<dim-1>::no_refinement
+ },
+
+ {
+ RefinementCase<dim-1>::no_refinement, // cut_x
+ RefinementCase<dim-1>::cut_y,
+ RefinementCase<dim-1>::cut_x
+ },
+
+ {
+ RefinementCase<dim-1>::cut_x, // cut_y
+ RefinementCase<dim-1>::no_refinement,
+ RefinementCase<dim-1>::cut_y
+ },
+
+ {
+ RefinementCase<dim-1>::cut_x, // cut_xy
+ RefinementCase<dim-1>::cut_y,
+ RefinementCase<dim-1>::cut_xy
+ },
+
+ {
+ RefinementCase<dim-1>::cut_y, // cut_z
+ RefinementCase<dim-1>::cut_x,
+ RefinementCase<dim-1>::no_refinement
+ },
+
+ {
+ RefinementCase<dim-1>::cut_y, // cut_xz
+ RefinementCase<dim-1>::cut_xy,
+ RefinementCase<dim-1>::cut_x
+ },
+
+ {
+ RefinementCase<dim-1>::cut_xy, // cut_yz
+ RefinementCase<dim-1>::cut_x,
+ RefinementCase<dim-1>::cut_y
+ },
+
+ {
+ RefinementCase<dim-1>::cut_xy, // cut_xyz
+ RefinementCase<dim-1>::cut_xy,
+ RefinementCase<dim-1>::cut_xy
+ },
+ };
+
+ const RefinementCase<dim-1> ref_case=ref_cases[cell_refinement_case][face_no/2];
+
+ const RefinementCase<dim-1> flip[4]=
+ {
+ RefinementCase<dim-1>::no_refinement,
+ RefinementCase<dim-1>::cut_y,
+ RefinementCase<dim-1>::cut_x,
+ RefinementCase<dim-1>::cut_xy
+ };
+
+ // correct the ref_case for face_orientation
+ // and face_rotation. for face_orientation,
+ // 'true' is the default value whereas for
+ // face_rotation, 'false' is standard. If
+ // <tt>face_rotation==face_orientation</tt>,
+ // then one of them is non-standard and we
+ // have to swap cut_x and cut_y, otherwise no
+ // change is necessary. face_flip has no
+ // influence. however, in order to keep the
+ // interface consistent with other functions,
+ // we still include it as an argument to this
+ // function
+ return (face_orientation==face_rotation) ? flip[ref_case] : ref_case;
+}
+
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<1>::line_refinement_case(const RefinementCase<1> &cell_refinement_case,
+ const unsigned int line_no)
+{
+ (void)line_no;
+ const unsigned int dim = 1;
+ (void)dim;
+ Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
+ Assert(line_no<GeometryInfo<dim>::lines_per_cell,
+ ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
+
+ return cell_refinement_case;
+}
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<2>::line_refinement_case(const RefinementCase<2> &cell_refinement_case,
+ const unsigned int line_no)
+{
+ // Assertions are in face_refinement_case()
+ return face_refinement_case(cell_refinement_case, line_no);
+}
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<3>::line_refinement_case(const RefinementCase<3> &cell_refinement_case,
+ const unsigned int line_no)
+{
+ const unsigned int dim=3;
+ Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
+ ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
+ Assert(line_no<GeometryInfo<dim>::lines_per_cell,
+ ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
+
+ // array indicating, which simple refine
+ // case cuts a line in direction x, y or
+ // z. For example, cut_y and everything
+ // containing cut_y (cut_xy, cut_yz,
+ // cut_xyz) cuts lines, which are in y
+ // direction.
+ const RefinementCase<dim>
+ cut_one[dim] =
+ {
+ RefinementCase<dim>::cut_x,
+ RefinementCase<dim>::cut_y,
+ RefinementCase<dim>::cut_z
+ };
+
+ // order the direction of lines
+ // 0->x, 1->y, 2->z
+ const unsigned int direction[lines_per_cell]=
+ {1,1,0,0,1,1,0,0,2,2,2,2};
+
+ return ((cell_refinement_case & cut_one[direction[line_no]]) ?
+ RefinementCase<1>::cut_x : RefinementCase<1>::no_refinement);
+}
+
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(const RefinementCase<0> &,
+ const unsigned int,
+ const bool,
+ const bool,
+ const bool)
+{
+ const unsigned int dim = 1;
+ Assert(false, ExcImpossibleInDim(dim));
+
+ return RefinementCase<dim>::no_refinement;
+}
+
+
+template <>
+inline
+RefinementCase<2>
+GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(const RefinementCase<1> &face_refinement_case,
+ const unsigned int face_no,
+ const bool,
+ const bool,
+ const bool)
+{
+ const unsigned int dim = 2;
+ Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
+ ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
+ Assert(face_no<GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
+
+ if (face_refinement_case==RefinementCase<dim>::cut_x)
+ return (face_no/2) ? RefinementCase<dim>::cut_x : RefinementCase<dim>::cut_y;
+ else
+ return RefinementCase<dim>::no_refinement;
+}
+
+
+template <>
+inline
+RefinementCase<3>
+GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(const RefinementCase<2> &face_refinement_case,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool /*face_flip*/,
+ const bool face_rotation)
+{
+ const unsigned int dim=3;
+ Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
+ ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
+ Assert(face_no<GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
+
+ const RefinementCase<2> flip[4]=
+ {
+ RefinementCase<2>::no_refinement,
+ RefinementCase<2>::cut_y,
+ RefinementCase<2>::cut_x,
+ RefinementCase<2>::cut_xy
+ };
+
+ // correct the face_refinement_case for
+ // face_orientation and face_rotation. for
+ // face_orientation, 'true' is the default
+ // value whereas for face_rotation, 'false'
+ // is standard. If
+ // <tt>face_rotation==face_orientation</tt>,
+ // then one of them is non-standard and we
+ // have to swap cut_x and cut_y, otherwise no
+ // change is necessary. face_flip has no
+ // influence. however, in order to keep the
+ // interface consistent with other functions,
+ // we still include it as an argument to this
+ // function
+ const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_refinement_case] : face_refinement_case;
+
+ const RefinementCase<dim> face_to_cell[3][4]=
+ {
+ {
+ RefinementCase<dim>::no_refinement, // faces 0 and 1
+ RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
+ RefinementCase<dim>::cut_z,
+ RefinementCase<dim>::cut_yz
+ },
+
+ {
+ RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are "exchanged on faces 2 and 3")
+ RefinementCase<dim>::cut_z,
+ RefinementCase<dim>::cut_x,
+ RefinementCase<dim>::cut_xz
+ },
+
+ {
+ RefinementCase<dim>::no_refinement, // faces 4 and 5
+ RefinementCase<dim>::cut_x,
+ RefinementCase<dim>::cut_y,
+ RefinementCase<dim>::cut_xy
+ }
+ };
+
+ return face_to_cell[face_no/2][std_face_ref];
+}
+
+
+
+template <>
+inline
+RefinementCase<1>
+GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
+{
+ (void)line_no;
+ Assert(line_no==0, ExcIndexRange(line_no,0,1));
+
+ return RefinementCase<1>::cut_x;
+}
+
+
+template <>
+inline
+RefinementCase<2>
+GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
+{
+ const unsigned int dim = 2;
+ (void)dim;
+ Assert(line_no<GeometryInfo<dim>::lines_per_cell,
+ ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
+
+ return (line_no/2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
+}
+
+
+template <>
+inline
+RefinementCase<3>
+GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
+{
+ const unsigned int dim=3;
+ Assert(line_no<GeometryInfo<dim>::lines_per_cell,
+ ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
+
+ const RefinementCase<dim> ref_cases[6]=
+ {
+ RefinementCase<dim>::cut_y, // lines 0 and 1
+ RefinementCase<dim>::cut_x, // lines 2 and 3
+ RefinementCase<dim>::cut_y, // lines 4 and 5
+ RefinementCase<dim>::cut_x, // lines 6 and 7
+ RefinementCase<dim>::cut_z, // lines 8 and 9
+ RefinementCase<dim>::cut_z
+ }; // lines 10 and 11
+
+ return ref_cases[line_no/2];
+}
+
+
+
+
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert(vertex<GeometryInfo<3>::vertices_per_face,
+ ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
+
+ // set up a table to make sure that
+ // we handle non-standard faces correctly
+ //
+ // so set up a table that for each vertex (of
+ // a quad in standard position) describes
+ // which vertex to take
+ //
+ // first index: four vertices 0...3
+ //
+ // second index: face_orientation; 0:
+ // opposite normal, 1: standard
+ //
+ // third index: face_flip; 0: standard, 1:
+ // face rotated by 180 degrees
+ //
+ // forth index: face_rotation: 0: standard,
+ // 1: face rotated by 90 degrees
+
+ const unsigned int vertex_translation[4][2][2][2] =
+ {
+ { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
+ { 3, 1 }
+ }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
+ { { 0, 1 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
+ { 3, 2 }
+ }
+ },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
+
+ { { { 2, 3 }, // vertex 1 ...
+ { 1, 0 }
+ },
+ { { 1, 3 },
+ { 2, 0 }
+ }
+ },
+
+ { { { 1, 0 }, // vertex 2 ...
+ { 2, 3 }
+ },
+ { { 2, 0 },
+ { 1, 3 }
+ }
+ },
+
+ { { { 3, 1 }, // vertex 3 ...
+ { 0, 2 }
+ },
+ { { 3, 2 },
+ { 0, 1 }
+ }
+ }
+ };
+
+ return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
+ const bool,
+ const bool,
+ const bool)
+{
+ Assert(dim>1, ExcImpossibleInDim(dim));
+ Assert(vertex<GeometryInfo<dim>::vertices_per_face,
+ ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
+ return vertex;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert(line<GeometryInfo<3>::lines_per_face,
+ ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
+
+
+ // make sure we handle
+ // non-standard faces correctly
+ //
+ // so set up a table that for each line (of a
+ // quad) describes which line to take
+ //
+ // first index: four lines 0...3
+ //
+ // second index: face_orientation; 0:
+ // opposite normal, 1: standard
+ //
+ // third index: face_flip; 0: standard, 1:
+ // face rotated by 180 degrees
+ //
+ // forth index: face_rotation: 0: standard,
+ // 1: face rotated by 90 degrees
+
+ const unsigned int line_translation[4][2][2][2] =
+ {
+ { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
+ { 3, 1 }
+ }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
+ { { 0, 3 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
+ { 1, 2 }
+ }
+ },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
+
+ { { { 3, 1 }, // line 1 ...
+ { 2, 0 }
+ },
+ { { 1, 2 },
+ { 0, 3 }
+ }
+ },
+
+ { { { 0, 3 }, // line 2 ...
+ { 1, 2 }
+ },
+ { { 2, 0 },
+ { 3, 1 }
+ }
+ },
+
+ { { { 1, 2 }, // line 3 ...
+ { 0, 3 }
+ },
+ { { 3, 1 },
+ { 2, 0 }
+ }
+ }
+ };
+
+ return line_translation[line][face_orientation][face_flip][face_rotation];
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
+ const bool,
+ const bool,
+ const bool)
+{
+ Assert(false, ExcNotImplemented());
+ return line;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert(line<GeometryInfo<3>::lines_per_face,
+ ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
+
+
+ // make sure we handle
+ // non-standard faces correctly
+ //
+ // so set up a table that for each line (of a
+ // quad) describes which line to take
+ //
+ // first index: four lines 0...3
+ //
+ // second index: face_orientation; 0:
+ // opposite normal, 1: standard
+ //
+ // third index: face_flip; 0: standard, 1:
+ // face rotated by 180 degrees
+ //
+ // forth index: face_rotation: 0: standard,
+ // 1: face rotated by 90 degrees
+
+ const unsigned int line_translation[4][2][2][2] =
+ {
+ { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
+ { 3, 1 }
+ }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
+ { { 0, 2 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
+ { 1, 3 }
+ }
+ },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
+
+ { { { 3, 1 }, // line 1 ...
+ { 2, 0 }
+ },
+ { { 1, 3 },
+ { 0, 2 }
+ }
+ },
+
+ { { { 0, 3 }, // line 2 ...
+ { 1, 2 }
+ },
+ { { 2, 1 },
+ { 3, 0 }
+ }
+ },
+
+ { { { 1, 2 }, // line 3 ...
+ { 0, 3 }
+ },
+ { { 3, 0 },
+ { 2, 1 }
+ }
+ }
+ };
+
+ return line_translation[line][face_orientation][face_flip][face_rotation];
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
+ const bool,
+ const bool,
+ const bool)
+{
+ Assert(false, ExcNotImplemented());
+ return line;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<1>::child_cell_on_face (const RefinementCase<1> &,
+ const unsigned int face,
+ const unsigned int subface,
+ const bool, const bool, const bool,
+ const RefinementCase<0> &)
+{
+ (void)subface;
+ Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (subface<max_children_per_face,
+ ExcIndexRange(subface, 0, max_children_per_face));
+
+ return face;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<2>::child_cell_on_face (const RefinementCase<2> &ref_case,
+ const unsigned int face,
+ const unsigned int subface,
+ const bool /*face_orientation*/,
+ const bool face_flip,
+ const bool /*face_rotation*/,
+ const RefinementCase<1> &)
+{
+ Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (subface<max_children_per_face,
+ ExcIndexRange(subface, 0, max_children_per_face));
+
+ // always return the child adjacent to the specified
+ // subface. if the face of a cell is not refined, don't
+ // throw an assertion but deliver the child adjacent to
+ // the face nevertheless, i.e. deliver the child of
+ // this cell adjacent to the subface of a possibly
+ // refined neighbor. this simplifies setting neighbor
+ // information in execute_refinement.
+ const unsigned int
+ subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell][max_children_per_face] =
+ {
+ {
+ // Normal orientation (face_flip = false)
+ {{0,0},{1,1},{0,1},{0,1}}, // cut_x
+ {{0,1},{0,1},{0,0},{1,1}}, // cut_y
+ {{0,2},{1,3},{0,1},{2,3}} // cut_xy, i.e., isotropic
+ },
+ {
+ // Flipped orientation (face_flip = true)
+ {{0,0},{1,1},{1,0},{1,0}}, // cut_x
+ {{1,0},{1,0},{0,0},{1,1}}, // cut_y
+ {{2,0},{3,1},{1,0},{3,2}} // cut_xy, i.e., isotropic
+ }
+ };
+
+ return subcells[face_flip][ref_case-1][face][subface];
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::child_cell_on_face (const RefinementCase<3> &ref_case,
+ const unsigned int face,
+ const unsigned int subface,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const RefinementCase<2> &face_ref_case)
+{
+ const unsigned int dim = 3;
+
+ Assert (ref_case>RefinementCase<dim-1>::no_refinement, ExcMessage("Cell has no children."));
+ Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (subface<GeometryInfo<dim-1>::n_children(face_ref_case) ||
+ (subface==0 && face_ref_case==RefinementCase<dim-1>::no_refinement),
+ ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
+
+ // invalid number used for invalid cases,
+ // e.g. when the children are more refined at
+ // a given face than the face itself
+ const unsigned int e = numbers::invalid_unsigned_int;
+
+ // the whole process of finding a child cell
+ // at a given subface considering the
+ // possibly anisotropic refinement cases of
+ // the cell and the face as well as
+ // orientation, flip and rotation of the face
+ // is quite complicated. thus, we break it
+ // down into several steps.
+
+ // first step: convert the given face refine
+ // case to a face refine case concerning the
+ // face in standard orientation (, flip and
+ // rotation). This only affects cut_x and
+ // cut_y
+ const RefinementCase<dim-1> flip[4]=
+ {
+ RefinementCase<dim-1>::no_refinement,
+ RefinementCase<dim-1>::cut_y,
+ RefinementCase<dim-1>::cut_x,
+ RefinementCase<dim-1>::cut_xy
+ };
+ // for face_orientation, 'true' is the
+ // default value whereas for face_rotation,
+ // 'false' is standard. If
+ // <tt>face_rotation==face_orientation</tt>,
+ // then one of them is non-standard and we
+ // have to swap cut_x and cut_y, otherwise no
+ // change is necessary.
+ const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_ref_case] : face_ref_case;
+
+ // second step: convert the given subface
+ // index to the one for a standard face
+ // respecting face_orientation, face_flip and
+ // face_rotation
+
+ // first index: face_ref_case
+ // second index: face_orientation
+ // third index: face_flip
+ // forth index: face_rotation
+ // fifth index: subface index
+ const unsigned int subface_exchange[4][2][2][2][4]=
+ {
+ // no_refinement (subface 0 stays 0,
+ // all others are invalid)
+ { { { {0,e,e,e},
+ {0,e,e,e}
+ },
+ { {0,e,e,e},
+ {0,e,e,e}
+ }
+ },
+ { { {0,e,e,e},
+ {0,e,e,e}
+ },
+ { {0,e,e,e},
+ {0,e,e,e}
+ }
+ }
+ },
+ // cut_x (here, if the face is only
+ // rotated OR only falsely oriented,
+ // then subface 0 of the non-standard
+ // face does NOT correspond to one of
+ // the subfaces of a standard
+ // face. Thus we indicate the subface
+ // which is located at the lower left
+ // corner (the origin of the face's
+ // local coordinate system) with
+ // '0'. The rest of this issue is
+ // taken care of using the above
+ // conversion to a 'standard face
+ // refine case')
+ { { { {0,1,e,e},
+ {0,1,e,e}
+ },
+ { {1,0,e,e},
+ {1,0,e,e}
+ }
+ },
+ { { {0,1,e,e},
+ {0,1,e,e}
+ },
+ { {1,0,e,e},
+ {1,0,e,e}
+ }
+ }
+ },
+ // cut_y (the same applies as for
+ // cut_x)
+ { { { {0,1,e,e},
+ {1,0,e,e}
+ },
+ { {1,0,e,e},
+ {0,1,e,e}
+ }
+ },
+ { { {0,1,e,e},
+ {1,0,e,e}
+ },
+ { {1,0,e,e},
+ {0,1,e,e}
+ }
+ }
+ },
+ // cut_xyz: this information is
+ // identical to the information
+ // returned by
+ // GeometryInfo<3>::real_to_standard_face_vertex()
+ { { { {0,2,1,3}, // face_orientation=false, face_flip=false, face_rotation=false, subfaces 0,1,2,3
+ {2,3,0,1}
+ }, // face_orientation=false, face_flip=false, face_rotation=true, subfaces 0,1,2,3
+ { {3,1,2,0}, // face_orientation=false, face_flip=true, face_rotation=false, subfaces 0,1,2,3
+ {1,0,3,2}
+ }
+ }, // face_orientation=false, face_flip=true, face_rotation=true, subfaces 0,1,2,3
+ { { {0,1,2,3}, // face_orientation=true, face_flip=false, face_rotation=false, subfaces 0,1,2,3
+ {1,3,0,2}
+ }, // face_orientation=true, face_flip=false, face_rotation=true, subfaces 0,1,2,3
+ { {3,2,1,0}, // face_orientation=true, face_flip=true, face_rotation=false, subfaces 0,1,2,3
+ {2,0,3,1}
+ }
+ }
+ }
+ };// face_orientation=true, face_flip=true, face_rotation=true, subfaces 0,1,2,3
+
+ const unsigned int std_subface=subface_exchange
+ [face_ref_case]
+ [face_orientation]
+ [face_flip]
+ [face_rotation]
+ [subface];
+ Assert (std_subface!=e, ExcInternalError());
+
+ // third step: these are the children, which
+ // can be found at the given subfaces of an
+ // isotropically refined (standard) face
+ //
+ // first index: (refinement_case-1)
+ // second index: face_index
+ // third index: subface_index (isotropic refinement)
+ const unsigned int
+ iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell][max_children_per_face] =
+ {
+ // cut_x
+ { {0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
+ {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
+ {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
+ {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
+ {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
+ {0, 1, 0, 1}
+ }, // face 5, subfaces 0,1,2,3
+ // cut_y
+ { {0, 1, 0, 1},
+ {0, 1, 0, 1},
+ {0, 0, 0, 0},
+ {1, 1, 1, 1},
+ {0, 0, 1, 1},
+ {0, 0, 1, 1}
+ },
+ // cut_xy
+ { {0, 2, 0, 2},
+ {1, 3, 1, 3},
+ {0, 0, 1, 1},
+ {2, 2, 3, 3},
+ {0, 1, 2, 3},
+ {0, 1, 2, 3}
+ },
+ // cut_z
+ { {0, 0, 1, 1},
+ {0, 0, 1, 1},
+ {0, 1, 0, 1},
+ {0, 1, 0, 1},
+ {0, 0, 0, 0},
+ {1, 1, 1, 1}
+ },
+ // cut_xz
+ { {0, 0, 1, 1},
+ {2, 2, 3, 3},
+ {0, 1, 2, 3},
+ {0, 1, 2, 3},
+ {0, 2, 0, 2},
+ {1, 3, 1, 3}
+ },
+ // cut_yz
+ { {0, 1, 2, 3},
+ {0, 1, 2, 3},
+ {0, 2, 0, 2},
+ {1, 3, 1, 3},
+ {0, 0, 1, 1},
+ {2, 2, 3, 3}
+ },
+ // cut_xyz
+ { {0, 2, 4, 6},
+ {1, 3, 5, 7},
+ {0, 4, 1, 5},
+ {2, 6, 3, 7},
+ {0, 1, 2, 3},
+ {4, 5, 6, 7}
+ }
+ };
+
+ // forth step: check, whether the given face
+ // refine case is valid for the given cell
+ // refine case. this is the case, if the
+ // given face refine case is at least as
+ // refined as the face is for the given cell
+ // refine case
+
+ // note, that we are considering standard
+ // face refinement cases here and thus must
+ // not pass the given orientation, flip and
+ // rotation flags
+ if ((std_face_ref & face_refinement_case(ref_case, face))
+ == face_refinement_case(ref_case, face))
+ {
+ // all is fine. for anisotropic face
+ // refine cases, select one of the
+ // isotropic subfaces which neighbors the
+ // same child
+
+ // first index: (standard) face refine case
+ // second index: subface index
+ const unsigned int equivalent_iso_subface[4][4]=
+ {
+ {0,e,e,e}, // no_refinement
+ {0,3,e,e}, // cut_x
+ {0,3,e,e}, // cut_y
+ {0,1,2,3}
+ }; // cut_xy
+
+ const unsigned int equ_std_subface
+ =equivalent_iso_subface[std_face_ref][std_subface];
+ Assert (equ_std_subface!=e, ExcInternalError());
+
+ return iso_children[ref_case-1][face][equ_std_subface];
+ }
+ else
+ {
+ // the face_ref_case was too coarse,
+ // throw an error
+ Assert(false,
+ ExcMessage("The face RefineCase is too coarse "
+ "for the given cell RefineCase."));
+ }
+ // we only get here in case of an error
+ return e;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<4>::child_cell_on_face (const RefinementCase<4> &,
+ const unsigned int,
+ const unsigned int,
+ const bool, const bool, const bool,
+ const RefinementCase<3> &)
+{
+ Assert(false, ExcNotImplemented());
+ return numbers::invalid_unsigned_int;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<1>::face_to_cell_lines (const unsigned int face,
+ const unsigned int line,
+ const bool, const bool, const bool)
+{
+ (void)face;
+ (void)line;
+ Assert (face+1<faces_per_cell+1, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (line+1<lines_per_face+1, ExcIndexRange(line, 0, lines_per_face));
+
+ // There is only a single line, so
+ // it must be this.
+ return 0;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<2>::face_to_cell_lines (const unsigned int face,
+ const unsigned int line,
+ const bool, const bool, const bool)
+{
+ (void)line;
+ Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
+
+ // The face is a line itself.
+ return face;
+}
+
+
+
+template <>
+inline
+unsigned int
+GeometryInfo<3>::face_to_cell_lines (const unsigned int face,
+ const unsigned int line,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
+ Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
+
+ const unsigned
+ lines[faces_per_cell][lines_per_face] = {{8,10, 0, 4}, // left face
+ {9,11, 1, 5}, // right face
+ {2, 6, 8, 9}, // front face
+ {3, 7,10,11}, // back face
+ {0, 1, 2, 3}, // bottom face
+ {4, 5, 6, 7}
+ };// top face
+ return lines[face][real_to_standard_face_line(line,
+ face_orientation,
+ face_flip,
+ face_rotation)];
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::face_to_cell_lines (const unsigned int,
+ const unsigned int,
+ const bool, const bool, const bool)
+{
+ Assert(false, ExcNotImplemented());
+ return numbers::invalid_unsigned_int;
+}
+
+
+
+template <int dim>
+inline
+unsigned int
+GeometryInfo<dim>::face_to_cell_vertices (const unsigned int face,
+ const unsigned int vertex,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation)
+{
+ return child_cell_on_face(RefinementCase<dim>::isotropic_refinement, face, vertex,
+ face_orientation, face_flip, face_rotation);
+}
+
+
+
+template <int dim>
+inline
+Point<dim>
+GeometryInfo<dim>::project_to_unit_cell (const Point<dim> &q)
+{
+ Point<dim> p = q;
+ for (unsigned int i=0; i<dim; i++)
+ if (p[i] < 0.) p[i] = 0.;
+ else if (p[i] > 1.) p[i] = 1.;
+
+ return p;
+}
+
+
+
+template <int dim>
+inline
+double
+GeometryInfo<dim>::distance_to_unit_cell (const Point<dim> &p)
+{
+ double result = 0.0;
+
+ for (unsigned int i=0; i<dim; i++)
+ if ((-p[i]) > result)
+ result = -p[i];
+ else if ((p[i]-1.) > result)
+ result = (p[i] - 1.);
+
+ return result;
+}
+
+
+
+template <int dim>
+inline
+double
+GeometryInfo<dim>::
+d_linear_shape_function (const Point<dim> &xi,
+ const unsigned int i)
+{
+ Assert (i < GeometryInfo<dim>::vertices_per_cell,
+ ExcIndexRange (i, 0, GeometryInfo<dim>::vertices_per_cell));
+
+ switch (dim)
+ {
+ case 1:
+ {
+ const double x = xi[0];
+ switch (i)
+ {
+ case 0:
+ return 1-x;
+ case 1:
+ return x;
+ }
+ break;
+ }
+
+ case 2:
+ {
+ const double x = xi[0];
+ const double y = xi[1];
+ switch (i)
+ {
+ case 0:
+ return (1-x)*(1-y);
+ case 1:
+ return x*(1-y);
+ case 2:
+ return (1-x)*y;
+ case 3:
+ return x*y;
+ }
+ break;
+ }
+
+ case 3:
+ {
+ const double x = xi[0];
+ const double y = xi[1];
+ const double z = xi[2];
+ switch (i)
+ {
+ case 0:
+ return (1-x)*(1-y)*(1-z);
+ case 1:
+ return x*(1-y)*(1-z);
+ case 2:
+ return (1-x)*y*(1-z);
+ case 3:
+ return x*y*(1-z);
+ case 4:
+ return (1-x)*(1-y)*z;
+ case 5:
+ return x*(1-y)*z;
+ case 6:
+ return (1-x)*y*z;
+ case 7:
+ return x*y*z;
+ }
+ break;
+ }
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ return -1e9;
+}
+
+
+
+template <>
+Tensor<1,1>
+inline
+GeometryInfo<1>::
+d_linear_shape_function_gradient (const Point<1> &,
+ const unsigned int i)
+{
+ Assert (i < GeometryInfo<1>::vertices_per_cell,
+ ExcIndexRange (i, 0, GeometryInfo<1>::vertices_per_cell));
+
+ switch (i)
+ {
+ case 0:
+ return Point<1>(-1.);
+ case 1:
+ return Point<1>(1.);
+ }
+
+ return Point<1>(-1e9);
+}
+
+
+
+template <>
+Tensor<1,2>
+inline
+GeometryInfo<2>::
+d_linear_shape_function_gradient (const Point<2> &xi,
+ const unsigned int i)
+{
+ Assert (i < GeometryInfo<2>::vertices_per_cell,
+ ExcIndexRange (i, 0, GeometryInfo<2>::vertices_per_cell));
+
+ const double x = xi[0];
+ const double y = xi[1];
+ switch (i)
+ {
+ case 0:
+ return Point<2>(-(1-y),-(1-x));
+ case 1:
+ return Point<2>(1-y,-x);
+ case 2:
+ return Point<2>(-y, 1-x);
+ case 3:
+ return Point<2>(y,x);
+ }
+ return Point<2> (-1e9, -1e9);
+}
+
+
+
+template <>
+Tensor<1,3>
+inline
+GeometryInfo<3>::
+d_linear_shape_function_gradient (const Point<3> &xi,
+ const unsigned int i)
+{
+ Assert (i < GeometryInfo<3>::vertices_per_cell,
+ ExcIndexRange (i, 0, GeometryInfo<3>::vertices_per_cell));
+
+ const double x = xi[0];
+ const double y = xi[1];
+ const double z = xi[2];
+ switch (i)
+ {
+ case 0:
+ return Point<3>(-(1-y)*(1-z),
+ -(1-x)*(1-z),
+ -(1-x)*(1-y));
+ case 1:
+ return Point<3>((1-y)*(1-z),
+ -x*(1-z),
+ -x*(1-y));
+ case 2:
+ return Point<3>(-y*(1-z),
+ (1-x)*(1-z),
+ -(1-x)*y);
+ case 3:
+ return Point<3>(y*(1-z),
+ x*(1-z),
+ -x*y);
+ case 4:
+ return Point<3>(-(1-y)*z,
+ -(1-x)*z,
+ (1-x)*(1-y));
+ case 5:
+ return Point<3>((1-y)*z,
+ -x*z,
+ x*(1-y));
+ case 6:
+ return Point<3>(-y*z,
+ (1-x)*z,
+ (1-x)*y);
+ case 7:
+ return Point<3>(y*z, x*z, x*y);
+ }
+
+ return Point<3> (-1e9, -1e9, -1e9);
+}
+
+
+
+template <int dim>
+inline
+Tensor<1,dim>
+GeometryInfo<dim>::
+d_linear_shape_function_gradient (const Point<dim> &,
+ const unsigned int)
+{
+ Assert (false, ExcNotImplemented());
+ return Tensor<1,dim>();
+}
+
+
+unsigned int
+inline
+GeometryInfo<0>::n_children(const RefinementCase<0> &)
+{
+ return 0;
+}
+
+
+namespace internal
+{
+ namespace GeometryInfoHelper
+ {
+ // wedge product of a single
+ // vector in 2d: we just have to
+ // rotate it by 90 degrees to the
+ // right
+ inline
+ Tensor<1,2>
+ wedge_product (const Tensor<1,2> (&derivative)[1])
+ {
+ Tensor<1,2> result;
+ result[0] = derivative[0][1];
+ result[1] = -derivative[0][0];
+
+ return result;
+ }
+
+
+ // wedge product of 2 vectors in
+ // 3d is the cross product
+ inline
+ Tensor<1,3>
+ wedge_product (const Tensor<1,3> (&derivative)[2])
+ {
+ return cross_product_3d (derivative[0], derivative[1]);
+ }
+
+
+ // wedge product of dim vectors
+ // in dim-d: that's the
+ // determinant of the matrix
+ template <int dim>
+ inline
+ Tensor<0,dim>
+ wedge_product (const Tensor<1,dim> (&derivative)[dim])
+ {
+ Tensor<2,dim> jacobian;
+ for (unsigned int i=0; i<dim; ++i)
+ jacobian[i] = derivative[i];
+
+ return determinant (jacobian);
+ }
+ }
+}
+
+
+template <int dim>
+template <int spacedim>
+inline
+void
+GeometryInfo<dim>::
+alternating_form_at_vertices
+#ifndef DEAL_II_CONSTEXPR_BUG
+(const Point<spacedim> (&vertices)[vertices_per_cell],
+ Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
+#else
+(const Point<spacedim> *vertices,
+ Tensor<spacedim-dim,spacedim> *forms)
+#endif
+{
+ // for each of the vertices,
+ // compute the alternating form
+ // of the mapped unit
+ // vectors. consider for
+ // example the case of a quad
+ // in spacedim==3: to do so, we
+ // need to see how the
+ // infinitesimal vectors
+ // (d\xi_1,0) and (0,d\xi_2)
+ // are transformed into
+ // spacedim-dimensional space
+ // and then form their cross
+ // product (i.e. the wedge product
+ // of two vectors). to this end, note
+ // that
+ // \vec x = sum_i \vec v_i phi_i(\vec xi)
+ // so the transformed vectors are
+ // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
+ // and
+ // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
+ // which boils down to the columns
+ // of the 3x2 matrix \grad_\xi x(\xi)
+ //
+ // a similar reasoning would
+ // hold for all dim,spacedim
+ // pairs -- we only have to
+ // compute the wedge product of
+ // the columns of the
+ // derivatives
+ for (unsigned int i=0; i<vertices_per_cell; ++i)
+ {
+ Tensor<1,spacedim> derivatives[dim];
+
+ for (unsigned int j=0; j<vertices_per_cell; ++j)
+ {
+ const Tensor<1,dim> grad_phi_j
+ = d_linear_shape_function_gradient (unit_cell_vertex(i),
+ j);
+ for (unsigned int l=0; l<dim; ++l)
+ derivatives[l] += vertices[j] * grad_phi_j[l];
+ }
+
+ forms[i] = internal::GeometryInfoHelper::wedge_product (derivatives);
+ }
+}
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
DEAL_II_NAMESPACE_OPEN
-unsigned int
-GeometryInfo<0>::n_children(const RefinementCase<0> &)
-{
- return 0;
-}
-
-
-
template <int dim> constexpr unsigned int GeometryInfo<dim>::max_children_per_cell;
template <int dim> constexpr unsigned int GeometryInfo<dim>::faces_per_cell;
template <int dim> constexpr unsigned int GeometryInfo<dim>::max_children_per_face;
GeometryInfo<0>::dx_to_deal;
-using namespace numbers;
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::n_children(const RefinementCase<dim> &ref_case)
-{
- static const unsigned int n_children[RefinementCase<3>::cut_xyz+1]=
- {0, 2, 2, 4, 2, 4, 4, 8};
-
- return n_children[ref_case];
-}
-
-
-template <>
-unsigned int
-GeometryInfo<1>::n_subfaces(const internal::SubfaceCase<1> &)
-{
- Assert(false, ExcImpossibleInDim(1));
- return 0;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<2>::n_subfaces(const internal::SubfaceCase<2> &subface_case)
-{
- return (subface_case == internal::SubfaceCase<2>::case_x) ? 2 : 0;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::n_subfaces(const internal::SubfaceCase<3> &subface_case)
-{
- static const unsigned int nsubs[internal::SubfaceCase<3>::case_isotropic+1]=
- {0, 2, 3, 3, 4, 2, 3, 3, 4, 4};
- return nsubs[subface_case];
-}
-
-
-template <>
-double
-GeometryInfo<1>::subface_ratio(const internal::SubfaceCase<1> &,
- const unsigned int)
-{
- return 1;
-}
-
-
-template <>
-double
-GeometryInfo<2>::subface_ratio(const internal::SubfaceCase<2> &subface_case,
- const unsigned int)
-{
- const unsigned int dim=2;
-
- double ratio=1;
- switch (subface_case)
- {
- case internal::SubfaceCase<dim>::case_none:
- // Here, an
- // Assert(false,ExcInternalError())
- // would be the right
- // choice, but
- // unfortunately the
- // current function is
- // also called for faces
- // without children (see
- // tests/fe/mapping.cc).
-// Assert(false, ExcMessage("Face has no subfaces."));
- // Furthermore, assign
- // following value as
- // otherwise the
- // bits/volume_x tests
- // break
- ratio=1./GeometryInfo<dim>::max_children_per_face;
- break;
- case internal::SubfaceCase<dim>::case_x:
- ratio=0.5;
- break;
- default:
- // there should be no
- // cases left
- Assert(false, ExcInternalError());
- break;
- }
-
- return ratio;
-}
-
-
-template <>
-double
-GeometryInfo<3>::subface_ratio(const internal::SubfaceCase<3> &subface_case,
- const unsigned int subface_no)
-{
- const unsigned int dim=3;
-
- double ratio=1;
- switch (subface_case)
- {
- case internal::SubfaceCase<dim>::case_none:
- // Here, an
- // Assert(false,ExcInternalError())
- // would be the right
- // choice, but
- // unfortunately the
- // current function is
- // also called for faces
- // without children (see
- // tests/bits/mesh_3d_16.cc). Add
- // following switch to
- // avoid diffs in
- // tests/bits/mesh_3d_16
- ratio=1./GeometryInfo<dim>::max_children_per_face;
- break;
- case internal::SubfaceCase<dim>::case_x:
- case internal::SubfaceCase<dim>::case_y:
- ratio=0.5;
- break;
- case internal::SubfaceCase<dim>::case_xy:
- case internal::SubfaceCase<dim>::case_x1y2y:
- case internal::SubfaceCase<dim>::case_y1x2x:
- ratio=0.25;
- break;
- case internal::SubfaceCase<dim>::case_x1y:
- case internal::SubfaceCase<dim>::case_y1x:
- if (subface_no<2)
- ratio=0.25;
- else
- ratio=0.5;
- break;
- case internal::SubfaceCase<dim>::case_x2y:
- case internal::SubfaceCase<dim>::case_y2x:
- if (subface_no==0)
- ratio=0.5;
- else
- ratio=0.25;
- break;
- default:
- // there should be no
- // cases left
- Assert(false, ExcInternalError());
- break;
- }
-
- return ratio;
-}
-
-
-
-template <>
-RefinementCase<0>
-GeometryInfo<1>::face_refinement_case(const RefinementCase<1> &,
- const unsigned int,
- const bool,
- const bool,
- const bool)
-{
- Assert(false, ExcImpossibleInDim(1));
-
- return RefinementCase<0>::no_refinement;
-}
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<2>::face_refinement_case(const RefinementCase<2> &cell_refinement_case,
- const unsigned int face_no,
- const bool,
- const bool,
- const bool)
-{
- const unsigned int dim=2;
- Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
- ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
- Assert(face_no<GeometryInfo<dim>::faces_per_cell,
- ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
- static const RefinementCase<dim-1>
- ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
- {
- {
- RefinementCase<dim-1>::no_refinement, // no_refinement
- RefinementCase<dim-1>::no_refinement
- },
-
- {
- RefinementCase<dim-1>::no_refinement,
- RefinementCase<dim-1>::cut_x
- },
-
- {
- RefinementCase<dim-1>::cut_x,
- RefinementCase<dim-1>::no_refinement
- },
-
- {
- RefinementCase<dim-1>::cut_x, // cut_xy
- RefinementCase<dim-1>::cut_x
- }
- };
-
- return ref_cases[cell_refinement_case][face_no/2];
-}
-
-
-template <>
-RefinementCase<2>
-GeometryInfo<3>::face_refinement_case(const RefinementCase<3> &cell_refinement_case,
- const unsigned int face_no,
- const bool face_orientation,
- const bool /*face_flip*/,
- const bool face_rotation)
-{
- const unsigned int dim=3;
- Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
- ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
- Assert(face_no<GeometryInfo<dim>::faces_per_cell,
- ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
- static const RefinementCase<dim-1>
- ref_cases[RefinementCase<dim>::isotropic_refinement+1][GeometryInfo<dim>::faces_per_cell/2]=
- {
- {
- RefinementCase<dim-1>::no_refinement, // no_refinement
- RefinementCase<dim-1>::no_refinement,
- RefinementCase<dim-1>::no_refinement
- },
-
- {
- RefinementCase<dim-1>::no_refinement, // cut_x
- RefinementCase<dim-1>::cut_y,
- RefinementCase<dim-1>::cut_x
- },
-
- {
- RefinementCase<dim-1>::cut_x, // cut_y
- RefinementCase<dim-1>::no_refinement,
- RefinementCase<dim-1>::cut_y
- },
-
- {
- RefinementCase<dim-1>::cut_x, // cut_xy
- RefinementCase<dim-1>::cut_y,
- RefinementCase<dim-1>::cut_xy
- },
-
- {
- RefinementCase<dim-1>::cut_y, // cut_z
- RefinementCase<dim-1>::cut_x,
- RefinementCase<dim-1>::no_refinement
- },
-
- {
- RefinementCase<dim-1>::cut_y, // cut_xz
- RefinementCase<dim-1>::cut_xy,
- RefinementCase<dim-1>::cut_x
- },
-
- {
- RefinementCase<dim-1>::cut_xy, // cut_yz
- RefinementCase<dim-1>::cut_x,
- RefinementCase<dim-1>::cut_y
- },
-
- {
- RefinementCase<dim-1>::cut_xy, // cut_xyz
- RefinementCase<dim-1>::cut_xy,
- RefinementCase<dim-1>::cut_xy
- },
- };
-
- const RefinementCase<dim-1> ref_case=ref_cases[cell_refinement_case][face_no/2];
-
- static const RefinementCase<dim-1> flip[4]=
- {
- RefinementCase<dim-1>::no_refinement,
- RefinementCase<dim-1>::cut_y,
- RefinementCase<dim-1>::cut_x,
- RefinementCase<dim-1>::cut_xy
- };
-
- // correct the ref_case for face_orientation
- // and face_rotation. for face_orientation,
- // 'true' is the default value whereas for
- // face_rotation, 'false' is standard. If
- // <tt>face_rotation==face_orientation</tt>,
- // then one of them is non-standard and we
- // have to swap cut_x and cut_y, otherwise no
- // change is necessary. face_flip has no
- // influence. however, in order to keep the
- // interface consistent with other functions,
- // we still include it as an argument to this
- // function
- return (face_orientation==face_rotation) ? flip[ref_case] : ref_case;
-}
-
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<1>::line_refinement_case(const RefinementCase<1> &cell_refinement_case,
- const unsigned int line_no)
-{
- (void)line_no;
- const unsigned int dim = 1;
- (void)dim;
- Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
- ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
- Assert(line_no<GeometryInfo<dim>::lines_per_cell,
- ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
-
- return cell_refinement_case;
-}
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<2>::line_refinement_case(const RefinementCase<2> &cell_refinement_case,
- const unsigned int line_no)
-{
- // Assertions are in face_refinement_case()
- return face_refinement_case(cell_refinement_case, line_no);
-}
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<3>::line_refinement_case(const RefinementCase<3> &cell_refinement_case,
- const unsigned int line_no)
-{
- const unsigned int dim=3;
- Assert(cell_refinement_case<RefinementCase<dim>::isotropic_refinement+1,
- ExcIndexRange(cell_refinement_case, 0, RefinementCase<dim>::isotropic_refinement+1));
- Assert(line_no<GeometryInfo<dim>::lines_per_cell,
- ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
-
- // array indicating, which simple refine
- // case cuts a line in direction x, y or
- // z. For example, cut_y and everything
- // containing cut_y (cut_xy, cut_yz,
- // cut_xyz) cuts lines, which are in y
- // direction.
- static const RefinementCase<dim>
- cut_one[dim] =
- {
- RefinementCase<dim>::cut_x,
- RefinementCase<dim>::cut_y,
- RefinementCase<dim>::cut_z
- };
-
- // order the direction of lines
- // 0->x, 1->y, 2->z
- static const unsigned int direction[lines_per_cell]=
- {1,1,0,0,1,1,0,0,2,2,2,2};
-
- return ((cell_refinement_case & cut_one[direction[line_no]]) ?
- RefinementCase<1>::cut_x : RefinementCase<1>::no_refinement);
-}
-
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<1>::min_cell_refinement_case_for_face_refinement(const RefinementCase<0> &,
- const unsigned int,
- const bool,
- const bool,
- const bool)
-{
- const unsigned int dim = 1;
- Assert(false, ExcImpossibleInDim(dim));
-
- return RefinementCase<dim>::no_refinement;
-}
-
-
-template <>
-RefinementCase<2>
-GeometryInfo<2>::min_cell_refinement_case_for_face_refinement(const RefinementCase<1> &face_refinement_case,
- const unsigned int face_no,
- const bool,
- const bool,
- const bool)
-{
- const unsigned int dim = 2;
- Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
- ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
- Assert(face_no<GeometryInfo<dim>::faces_per_cell,
- ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
- if (face_refinement_case==RefinementCase<dim>::cut_x)
- return (face_no/2) ? RefinementCase<dim>::cut_x : RefinementCase<dim>::cut_y;
- else
- return RefinementCase<dim>::no_refinement;
-}
-
-
-template <>
-RefinementCase<3>
-GeometryInfo<3>::min_cell_refinement_case_for_face_refinement(const RefinementCase<2> &face_refinement_case,
- const unsigned int face_no,
- const bool face_orientation,
- const bool /*face_flip*/,
- const bool face_rotation)
-{
- const unsigned int dim=3;
- Assert(face_refinement_case<RefinementCase<dim-1>::isotropic_refinement+1,
- ExcIndexRange(face_refinement_case, 0, RefinementCase<dim-1>::isotropic_refinement+1));
- Assert(face_no<GeometryInfo<dim>::faces_per_cell,
- ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
- static const RefinementCase<2> flip[4]=
- {
- RefinementCase<2>::no_refinement,
- RefinementCase<2>::cut_y,
- RefinementCase<2>::cut_x,
- RefinementCase<2>::cut_xy
- };
-
- // correct the face_refinement_case for
- // face_orientation and face_rotation. for
- // face_orientation, 'true' is the default
- // value whereas for face_rotation, 'false'
- // is standard. If
- // <tt>face_rotation==face_orientation</tt>,
- // then one of them is non-standard and we
- // have to swap cut_x and cut_y, otherwise no
- // change is necessary. face_flip has no
- // influence. however, in order to keep the
- // interface consistent with other functions,
- // we still include it as an argument to this
- // function
- const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_refinement_case] : face_refinement_case;
-
- static const RefinementCase<dim> face_to_cell[3][4]=
- {
- {
- RefinementCase<dim>::no_refinement, // faces 0 and 1
- RefinementCase<dim>::cut_y, // cut_x in face 0 means cut_y for the cell
- RefinementCase<dim>::cut_z,
- RefinementCase<dim>::cut_yz
- },
-
- {
- RefinementCase<dim>::no_refinement, // faces 2 and 3 (note that x and y are "exchanged on faces 2 and 3")
- RefinementCase<dim>::cut_z,
- RefinementCase<dim>::cut_x,
- RefinementCase<dim>::cut_xz
- },
-
- {
- RefinementCase<dim>::no_refinement, // faces 4 and 5
- RefinementCase<dim>::cut_x,
- RefinementCase<dim>::cut_y,
- RefinementCase<dim>::cut_xy
- }
- };
-
- return face_to_cell[face_no/2][std_face_ref];
-}
-
-
-
-template <>
-RefinementCase<1>
-GeometryInfo<1>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
-{
- (void)line_no;
- Assert(line_no==0, ExcIndexRange(line_no,0,1));
-
- return RefinementCase<1>::cut_x;
-}
-
-
-template <>
-RefinementCase<2>
-GeometryInfo<2>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
-{
- const unsigned int dim = 2;
- (void)dim;
- Assert(line_no<GeometryInfo<dim>::lines_per_cell,
- ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
-
- return (line_no/2) ? RefinementCase<2>::cut_x : RefinementCase<2>::cut_y;
-}
-
-
-template <>
-RefinementCase<3>
-GeometryInfo<3>::min_cell_refinement_case_for_line_refinement(const unsigned int line_no)
-{
- const unsigned int dim=3;
- Assert(line_no<GeometryInfo<dim>::lines_per_cell,
- ExcIndexRange(line_no, 0, GeometryInfo<dim>::lines_per_cell));
-
- static const RefinementCase<dim> ref_cases[6]=
- {
- RefinementCase<dim>::cut_y, // lines 0 and 1
- RefinementCase<dim>::cut_x, // lines 2 and 3
- RefinementCase<dim>::cut_y, // lines 4 and 5
- RefinementCase<dim>::cut_x, // lines 6 and 7
- RefinementCase<dim>::cut_z, // lines 8 and 9
- RefinementCase<dim>::cut_z
- }; // lines 10 and 11
-
- return ref_cases[line_no/2];
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::standard_to_real_face_vertex(const unsigned int vertex,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- Assert(vertex<GeometryInfo<3>::vertices_per_face,
- ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
-
- // set up a table to make sure that
- // we handle non-standard faces correctly
- //
- // so set up a table that for each vertex (of
- // a quad in standard position) describes
- // which vertex to take
- //
- // first index: four vertices 0...3
- //
- // second index: face_orientation; 0:
- // opposite normal, 1: standard
- //
- // third index: face_flip; 0: standard, 1:
- // face rotated by 180 degrees
- //
- // forth index: face_rotation: 0: standard,
- // 1: face rotated by 90 degrees
-
- static const unsigned int vertex_translation[4][2][2][2] =
- {
- { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
- { 3, 1 }
- }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
- { { 0, 2 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
- { 3, 1 }
- }
- },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
-
- { { { 2, 3 }, // vertex 1 ...
- { 1, 0 }
- },
- { { 1, 0 },
- { 2, 3 }
- }
- },
-
- { { { 1, 0 }, // vertex 2 ...
- { 2, 3 }
- },
- { { 2, 3 },
- { 1, 0 }
- }
- },
-
- { { { 3, 1 }, // vertex 3 ...
- { 0, 2 }
- },
- { { 3, 1 },
- { 0, 2 }
- }
- }
- };
-
- return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::standard_to_real_face_vertex(const unsigned int vertex,
- const bool,
- const bool,
- const bool)
-{
- Assert(dim>1, ExcImpossibleInDim(dim));
- Assert(vertex<GeometryInfo<dim>::vertices_per_face,
- ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
- return vertex;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::real_to_standard_face_vertex(const unsigned int vertex,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- Assert(vertex<GeometryInfo<3>::vertices_per_face,
- ExcIndexRange(vertex,0,GeometryInfo<3>::vertices_per_face));
-
- // set up a table to make sure that
- // we handle non-standard faces correctly
- //
- // so set up a table that for each vertex (of
- // a quad in standard position) describes
- // which vertex to take
- //
- // first index: four vertices 0...3
- //
- // second index: face_orientation; 0:
- // opposite normal, 1: standard
- //
- // third index: face_flip; 0: standard, 1:
- // face rotated by 180 degrees
- //
- // forth index: face_rotation: 0: standard,
- // 1: face rotated by 90 degrees
-
- static const unsigned int vertex_translation[4][2][2][2] =
- {
- { { { 0, 2 }, // vertex 0, face_orientation=false, face_flip=false, face_rotation=false and true
- { 3, 1 }
- }, // vertex 0, face_orientation=false, face_flip=true, face_rotation=false and true
- { { 0, 1 }, // vertex 0, face_orientation=true, face_flip=false, face_rotation=false and true
- { 3, 2 }
- }
- },// vertex 0, face_orientation=true, face_flip=true, face_rotation=false and true
-
- { { { 2, 3 }, // vertex 1 ...
- { 1, 0 }
- },
- { { 1, 3 },
- { 2, 0 }
- }
- },
-
- { { { 1, 0 }, // vertex 2 ...
- { 2, 3 }
- },
- { { 2, 0 },
- { 1, 3 }
- }
- },
-
- { { { 3, 1 }, // vertex 3 ...
- { 0, 2 }
- },
- { { 3, 2 },
- { 0, 1 }
- }
- }
- };
-
- return vertex_translation[vertex][face_orientation][face_flip][face_rotation];
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::real_to_standard_face_vertex(const unsigned int vertex,
- const bool,
- const bool,
- const bool)
-{
- Assert(dim>1, ExcImpossibleInDim(dim));
- Assert(vertex<GeometryInfo<dim>::vertices_per_face,
- ExcIndexRange(vertex,0,GeometryInfo<dim>::vertices_per_face));
- return vertex;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::standard_to_real_face_line(const unsigned int line,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- Assert(line<GeometryInfo<3>::lines_per_face,
- ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
-
-
- // make sure we handle
- // non-standard faces correctly
- //
- // so set up a table that for each line (of a
- // quad) describes which line to take
- //
- // first index: four lines 0...3
- //
- // second index: face_orientation; 0:
- // opposite normal, 1: standard
- //
- // third index: face_flip; 0: standard, 1:
- // face rotated by 180 degrees
- //
- // forth index: face_rotation: 0: standard,
- // 1: face rotated by 90 degrees
-
- static const unsigned int line_translation[4][2][2][2] =
- {
- { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
- { 3, 1 }
- }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
- { { 0, 3 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
- { 1, 2 }
- }
- },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
-
- { { { 3, 1 }, // line 1 ...
- { 2, 0 }
- },
- { { 1, 2 },
- { 0, 3 }
- }
- },
-
- { { { 0, 3 }, // line 2 ...
- { 1, 2 }
- },
- { { 2, 0 },
- { 3, 1 }
- }
- },
-
- { { { 1, 2 }, // line 3 ...
- { 0, 3 }
- },
- { { 3, 1 },
- { 2, 0 }
- }
- }
- };
-
- return line_translation[line][face_orientation][face_flip][face_rotation];
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::standard_to_real_face_line(const unsigned int line,
- const bool,
- const bool,
- const bool)
-{
- Assert(false, ExcNotImplemented());
- return line;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::real_to_standard_face_line(const unsigned int line,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- Assert(line<GeometryInfo<3>::lines_per_face,
- ExcIndexRange(line,0,GeometryInfo<3>::lines_per_face));
-
-
- // make sure we handle
- // non-standard faces correctly
- //
- // so set up a table that for each line (of a
- // quad) describes which line to take
- //
- // first index: four lines 0...3
- //
- // second index: face_orientation; 0:
- // opposite normal, 1: standard
- //
- // third index: face_flip; 0: standard, 1:
- // face rotated by 180 degrees
- //
- // forth index: face_rotation: 0: standard,
- // 1: face rotated by 90 degrees
-
- static const unsigned int line_translation[4][2][2][2] =
- {
- { { { 2, 0 }, // line 0, face_orientation=false, face_flip=false, face_rotation=false and true
- { 3, 1 }
- }, // line 0, face_orientation=false, face_flip=true, face_rotation=false and true
- { { 0, 2 }, // line 0, face_orientation=true, face_flip=false, face_rotation=false and true
- { 1, 3 }
- }
- },// line 0, face_orientation=true, face_flip=true, face_rotation=false and true
-
- { { { 3, 1 }, // line 1 ...
- { 2, 0 }
- },
- { { 1, 3 },
- { 0, 2 }
- }
- },
-
- { { { 0, 3 }, // line 2 ...
- { 1, 2 }
- },
- { { 2, 1 },
- { 3, 0 }
- }
- },
-
- { { { 1, 2 }, // line 3 ...
- { 0, 3 }
- },
- { { 3, 0 },
- { 2, 1 }
- }
- }
- };
-
- return line_translation[line][face_orientation][face_flip][face_rotation];
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::real_to_standard_face_line(const unsigned int line,
- const bool,
- const bool,
- const bool)
-{
- Assert(false, ExcNotImplemented());
- return line;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<1>::child_cell_on_face (const RefinementCase<1> &,
- const unsigned int face,
- const unsigned int subface,
- const bool, const bool, const bool,
- const RefinementCase<0> &)
-{
- (void)subface;
- Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
- Assert (subface<max_children_per_face,
- ExcIndexRange(subface, 0, max_children_per_face));
-
- return face;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<2>::child_cell_on_face (const RefinementCase<2> &ref_case,
- const unsigned int face,
- const unsigned int subface,
- const bool /*face_orientation*/,
- const bool face_flip,
- const bool /*face_rotation*/,
- const RefinementCase<1> &)
-{
- Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
- Assert (subface<max_children_per_face,
- ExcIndexRange(subface, 0, max_children_per_face));
-
- // always return the child adjacent to the specified
- // subface. if the face of a cell is not refined, don't
- // throw an assertion but deliver the child adjacent to
- // the face nevertheless, i.e. deliver the child of
- // this cell adjacent to the subface of a possibly
- // refined neighbor. this simplifies setting neighbor
- // information in execute_refinement.
- static const unsigned int
- subcells[2][RefinementCase<2>::isotropic_refinement][faces_per_cell][max_children_per_face] =
- {
- {
- // Normal orientation (face_flip = false)
- {{0,0},{1,1},{0,1},{0,1}}, // cut_x
- {{0,1},{0,1},{0,0},{1,1}}, // cut_y
- {{0,2},{1,3},{0,1},{2,3}} // cut_xy, i.e., isotropic
- },
- {
- // Flipped orientation (face_flip = true)
- {{0,0},{1,1},{1,0},{1,0}}, // cut_x
- {{1,0},{1,0},{0,0},{1,1}}, // cut_y
- {{2,0},{3,1},{1,0},{3,2}} // cut_xy, i.e., isotropic
- }
- };
-
- return subcells[face_flip][ref_case-1][face][subface];
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::child_cell_on_face (const RefinementCase<3> &ref_case,
- const unsigned int face,
- const unsigned int subface,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation,
- const RefinementCase<2> &face_ref_case)
-{
- const unsigned int dim = 3;
-
- Assert (ref_case>RefinementCase<dim-1>::no_refinement, ExcMessage("Cell has no children."));
- Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
- Assert (subface<GeometryInfo<dim-1>::n_children(face_ref_case) ||
- (subface==0 && face_ref_case==RefinementCase<dim-1>::no_refinement),
- ExcIndexRange(subface, 0, GeometryInfo<2>::n_children(face_ref_case)));
-
- // invalid number used for invalid cases,
- // e.g. when the children are more refined at
- // a given face than the face itself
- static const unsigned int e=invalid_unsigned_int;
-
- // the whole process of finding a child cell
- // at a given subface considering the
- // possibly anisotropic refinement cases of
- // the cell and the face as well as
- // orientation, flip and rotation of the face
- // is quite complicated. thus, we break it
- // down into several steps.
-
- // first step: convert the given face refine
- // case to a face refine case concerning the
- // face in standard orientation (, flip and
- // rotation). This only affects cut_x and
- // cut_y
- static const RefinementCase<dim-1> flip[4]=
- {
- RefinementCase<dim-1>::no_refinement,
- RefinementCase<dim-1>::cut_y,
- RefinementCase<dim-1>::cut_x,
- RefinementCase<dim-1>::cut_xy
- };
- // for face_orientation, 'true' is the
- // default value whereas for face_rotation,
- // 'false' is standard. If
- // <tt>face_rotation==face_orientation</tt>,
- // then one of them is non-standard and we
- // have to swap cut_x and cut_y, otherwise no
- // change is necessary.
- const RefinementCase<dim-1> std_face_ref = (face_orientation==face_rotation) ? flip[face_ref_case] : face_ref_case;
-
- // second step: convert the given subface
- // index to the one for a standard face
- // respecting face_orientation, face_flip and
- // face_rotation
-
- // first index: face_ref_case
- // second index: face_orientation
- // third index: face_flip
- // forth index: face_rotation
- // fifth index: subface index
- static const unsigned int subface_exchange[4][2][2][2][4]=
- {
- // no_refinement (subface 0 stays 0,
- // all others are invalid)
- { { { {0,e,e,e},
- {0,e,e,e}
- },
- { {0,e,e,e},
- {0,e,e,e}
- }
- },
- { { {0,e,e,e},
- {0,e,e,e}
- },
- { {0,e,e,e},
- {0,e,e,e}
- }
- }
- },
- // cut_x (here, if the face is only
- // rotated OR only falsely oriented,
- // then subface 0 of the non-standard
- // face does NOT correspond to one of
- // the subfaces of a standard
- // face. Thus we indicate the subface
- // which is located at the lower left
- // corner (the origin of the face's
- // local coordinate system) with
- // '0'. The rest of this issue is
- // taken care of using the above
- // conversion to a 'standard face
- // refine case')
- { { { {0,1,e,e},
- {0,1,e,e}
- },
- { {1,0,e,e},
- {1,0,e,e}
- }
- },
- { { {0,1,e,e},
- {0,1,e,e}
- },
- { {1,0,e,e},
- {1,0,e,e}
- }
- }
- },
- // cut_y (the same applies as for
- // cut_x)
- { { { {0,1,e,e},
- {1,0,e,e}
- },
- { {1,0,e,e},
- {0,1,e,e}
- }
- },
- { { {0,1,e,e},
- {1,0,e,e}
- },
- { {1,0,e,e},
- {0,1,e,e}
- }
- }
- },
- // cut_xyz: this information is
- // identical to the information
- // returned by
- // GeometryInfo<3>::real_to_standard_face_vertex()
- { { { {0,2,1,3}, // face_orientation=false, face_flip=false, face_rotation=false, subfaces 0,1,2,3
- {2,3,0,1}
- }, // face_orientation=false, face_flip=false, face_rotation=true, subfaces 0,1,2,3
- { {3,1,2,0}, // face_orientation=false, face_flip=true, face_rotation=false, subfaces 0,1,2,3
- {1,0,3,2}
- }
- }, // face_orientation=false, face_flip=true, face_rotation=true, subfaces 0,1,2,3
- { { {0,1,2,3}, // face_orientation=true, face_flip=false, face_rotation=false, subfaces 0,1,2,3
- {1,3,0,2}
- }, // face_orientation=true, face_flip=false, face_rotation=true, subfaces 0,1,2,3
- { {3,2,1,0}, // face_orientation=true, face_flip=true, face_rotation=false, subfaces 0,1,2,3
- {2,0,3,1}
- }
- }
- }
- };// face_orientation=true, face_flip=true, face_rotation=true, subfaces 0,1,2,3
-
- const unsigned int std_subface=subface_exchange
- [face_ref_case]
- [face_orientation]
- [face_flip]
- [face_rotation]
- [subface];
- Assert (std_subface!=e, ExcInternalError());
-
- // third step: these are the children, which
- // can be found at the given subfaces of an
- // isotropically refined (standard) face
- //
- // first index: (refinement_case-1)
- // second index: face_index
- // third index: subface_index (isotropic refinement)
- static const unsigned int
- iso_children[RefinementCase<dim>::cut_xyz][faces_per_cell][max_children_per_face] =
- {
- // cut_x
- { {0, 0, 0, 0}, // face 0, subfaces 0,1,2,3
- {1, 1, 1, 1}, // face 1, subfaces 0,1,2,3
- {0, 0, 1, 1}, // face 2, subfaces 0,1,2,3
- {0, 0, 1, 1}, // face 3, subfaces 0,1,2,3
- {0, 1, 0, 1}, // face 4, subfaces 0,1,2,3
- {0, 1, 0, 1}
- }, // face 5, subfaces 0,1,2,3
- // cut_y
- { {0, 1, 0, 1},
- {0, 1, 0, 1},
- {0, 0, 0, 0},
- {1, 1, 1, 1},
- {0, 0, 1, 1},
- {0, 0, 1, 1}
- },
- // cut_xy
- { {0, 2, 0, 2},
- {1, 3, 1, 3},
- {0, 0, 1, 1},
- {2, 2, 3, 3},
- {0, 1, 2, 3},
- {0, 1, 2, 3}
- },
- // cut_z
- { {0, 0, 1, 1},
- {0, 0, 1, 1},
- {0, 1, 0, 1},
- {0, 1, 0, 1},
- {0, 0, 0, 0},
- {1, 1, 1, 1}
- },
- // cut_xz
- { {0, 0, 1, 1},
- {2, 2, 3, 3},
- {0, 1, 2, 3},
- {0, 1, 2, 3},
- {0, 2, 0, 2},
- {1, 3, 1, 3}
- },
- // cut_yz
- { {0, 1, 2, 3},
- {0, 1, 2, 3},
- {0, 2, 0, 2},
- {1, 3, 1, 3},
- {0, 0, 1, 1},
- {2, 2, 3, 3}
- },
- // cut_xyz
- { {0, 2, 4, 6},
- {1, 3, 5, 7},
- {0, 4, 1, 5},
- {2, 6, 3, 7},
- {0, 1, 2, 3},
- {4, 5, 6, 7}
- }
- };
-
- // forth step: check, whether the given face
- // refine case is valid for the given cell
- // refine case. this is the case, if the
- // given face refine case is at least as
- // refined as the face is for the given cell
- // refine case
-
- // note, that we are considering standard
- // face refinement cases here and thus must
- // not pass the given orientation, flip and
- // rotation flags
- if ((std_face_ref & face_refinement_case(ref_case, face))
- == face_refinement_case(ref_case, face))
- {
- // all is fine. for anisotropic face
- // refine cases, select one of the
- // isotropic subfaces which neighbors the
- // same child
-
- // first index: (standard) face refine case
- // second index: subface index
- static const unsigned int equivalent_iso_subface[4][4]=
- {
- {0,e,e,e}, // no_refinement
- {0,3,e,e}, // cut_x
- {0,3,e,e}, // cut_y
- {0,1,2,3}
- }; // cut_xy
-
- const unsigned int equ_std_subface
- =equivalent_iso_subface[std_face_ref][std_subface];
- Assert (equ_std_subface!=e, ExcInternalError());
-
- return iso_children[ref_case-1][face][equ_std_subface];
- }
- else
- {
- // the face_ref_case was too coarse,
- // throw an error
- Assert(false,
- ExcMessage("The face RefineCase is too coarse "
- "for the given cell RefineCase."));
- }
- // we only get here in case of an error
- return e;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<4>::child_cell_on_face (const RefinementCase<4> &,
- const unsigned int,
- const unsigned int,
- const bool, const bool, const bool,
- const RefinementCase<3> &)
-{
- Assert(false, ExcNotImplemented());
- return invalid_unsigned_int;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<1>::face_to_cell_lines (const unsigned int face,
- const unsigned int line,
- const bool, const bool, const bool)
-{
- (void)face;
- (void)line;
- Assert (face+1<faces_per_cell+1, ExcIndexRange(face, 0, faces_per_cell));
- Assert (line+1<lines_per_face+1, ExcIndexRange(line, 0, lines_per_face));
-
- // There is only a single line, so
- // it must be this.
- return 0;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<2>::face_to_cell_lines (const unsigned int face,
- const unsigned int line,
- const bool, const bool, const bool)
-{
- (void)line;
- Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
- Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
-
- // The face is a line itself.
- return face;
-}
-
-
-
-template <>
-unsigned int
-GeometryInfo<3>::face_to_cell_lines (const unsigned int face,
- const unsigned int line,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- Assert (face<faces_per_cell, ExcIndexRange(face, 0, faces_per_cell));
- Assert (line<lines_per_face, ExcIndexRange(line, 0, lines_per_face));
-
- static const unsigned
- lines[faces_per_cell][lines_per_face] = {{8,10, 0, 4}, // left face
- {9,11, 1, 5}, // right face
- {2, 6, 8, 9}, // front face
- {3, 7,10,11}, // back face
- {0, 1, 2, 3}, // bottom face
- {4, 5, 6, 7}
- };// top face
- return lines[face][real_to_standard_face_line(line,
- face_orientation,
- face_flip,
- face_rotation)];
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::face_to_cell_lines (const unsigned int,
- const unsigned int,
- const bool, const bool, const bool)
-{
- Assert(false, ExcNotImplemented());
- return invalid_unsigned_int;
-}
-
-
-
-template <int dim>
-unsigned int
-GeometryInfo<dim>::face_to_cell_vertices (const unsigned int face,
- const unsigned int vertex,
- const bool face_orientation,
- const bool face_flip,
- const bool face_rotation)
-{
- return child_cell_on_face(RefinementCase<dim>::isotropic_refinement, face, vertex,
- face_orientation, face_flip, face_rotation);
-}
-
-
-
-template <int dim>
-Point<dim>
-GeometryInfo<dim>::project_to_unit_cell (const Point<dim> &q)
-{
- Point<dim> p = q;
- for (unsigned int i=0; i<dim; i++)
- if (p[i] < 0.) p[i] = 0.;
- else if (p[i] > 1.) p[i] = 1.;
-
- return p;
-}
-
-
-
-template <int dim>
-double
-GeometryInfo<dim>::distance_to_unit_cell (const Point<dim> &p)
-{
- double result = 0.0;
-
- for (unsigned int i=0; i<dim; i++)
- if ((-p[i]) > result)
- result = -p[i];
- else if ((p[i]-1.) > result)
- result = (p[i] - 1.);
-
- return result;
-}
-
-
-
-template <int dim>
-double
-GeometryInfo<dim>::
-d_linear_shape_function (const Point<dim> &xi,
- const unsigned int i)
-{
- Assert (i < GeometryInfo<dim>::vertices_per_cell,
- ExcIndexRange (i, 0, GeometryInfo<dim>::vertices_per_cell));
-
- switch (dim)
- {
- case 1:
- {
- const double x = xi[0];
- switch (i)
- {
- case 0:
- return 1-x;
- case 1:
- return x;
- }
- break;
- }
-
- case 2:
- {
- const double x = xi[0];
- const double y = xi[1];
- switch (i)
- {
- case 0:
- return (1-x)*(1-y);
- case 1:
- return x*(1-y);
- case 2:
- return (1-x)*y;
- case 3:
- return x*y;
- }
- break;
- }
-
- case 3:
- {
- const double x = xi[0];
- const double y = xi[1];
- const double z = xi[2];
- switch (i)
- {
- case 0:
- return (1-x)*(1-y)*(1-z);
- case 1:
- return x*(1-y)*(1-z);
- case 2:
- return (1-x)*y*(1-z);
- case 3:
- return x*y*(1-z);
- case 4:
- return (1-x)*(1-y)*z;
- case 5:
- return x*(1-y)*z;
- case 6:
- return (1-x)*y*z;
- case 7:
- return x*y*z;
- }
- break;
- }
-
- default:
- Assert (false, ExcNotImplemented());
- }
- return -1e9;
-}
-
-
-
-template <>
-Tensor<1,1>
-GeometryInfo<1>::
-d_linear_shape_function_gradient (const Point<1> &,
- const unsigned int i)
-{
- Assert (i < GeometryInfo<1>::vertices_per_cell,
- ExcIndexRange (i, 0, GeometryInfo<1>::vertices_per_cell));
-
- switch (i)
- {
- case 0:
- return Point<1>(-1.);
- case 1:
- return Point<1>(1.);
- }
-
- return Point<1>(-1e9);
-}
-
-
-
-template <>
-Tensor<1,2>
-GeometryInfo<2>::
-d_linear_shape_function_gradient (const Point<2> &xi,
- const unsigned int i)
-{
- Assert (i < GeometryInfo<2>::vertices_per_cell,
- ExcIndexRange (i, 0, GeometryInfo<2>::vertices_per_cell));
-
- const double x = xi[0];
- const double y = xi[1];
- switch (i)
- {
- case 0:
- return Point<2>(-(1-y),-(1-x));
- case 1:
- return Point<2>(1-y,-x);
- case 2:
- return Point<2>(-y, 1-x);
- case 3:
- return Point<2>(y,x);
- }
- return Point<2> (-1e9, -1e9);
-}
-
-
-
-template <>
-Tensor<1,3>
-GeometryInfo<3>::
-d_linear_shape_function_gradient (const Point<3> &xi,
- const unsigned int i)
-{
- Assert (i < GeometryInfo<3>::vertices_per_cell,
- ExcIndexRange (i, 0, GeometryInfo<3>::vertices_per_cell));
-
- const double x = xi[0];
- const double y = xi[1];
- const double z = xi[2];
- switch (i)
- {
- case 0:
- return Point<3>(-(1-y)*(1-z),
- -(1-x)*(1-z),
- -(1-x)*(1-y));
- case 1:
- return Point<3>((1-y)*(1-z),
- -x*(1-z),
- -x*(1-y));
- case 2:
- return Point<3>(-y*(1-z),
- (1-x)*(1-z),
- -(1-x)*y);
- case 3:
- return Point<3>(y*(1-z),
- x*(1-z),
- -x*y);
- case 4:
- return Point<3>(-(1-y)*z,
- -(1-x)*z,
- (1-x)*(1-y));
- case 5:
- return Point<3>((1-y)*z,
- -x*z,
- x*(1-y));
- case 6:
- return Point<3>(-y*z,
- (1-x)*z,
- (1-x)*y);
- case 7:
- return Point<3>(y*z, x*z, x*y);
- }
-
- return Point<3> (-1e9, -1e9, -1e9);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-GeometryInfo<dim>::
-d_linear_shape_function_gradient (const Point<dim> &,
- const unsigned int)
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,dim>();
-}
-
-
-
-
-
-namespace internal
-{
- namespace GeometryInfo
- {
- // wedge product of a single
- // vector in 2d: we just have to
- // rotate it by 90 degrees to the
- // right
- inline
- Tensor<1,2>
- wedge_product (const Tensor<1,2> (&derivative)[1])
- {
- Tensor<1,2> result;
- result[0] = derivative[0][1];
- result[1] = -derivative[0][0];
-
- return result;
- }
-
-
- // wedge product of 2 vectors in
- // 3d is the cross product
- inline
- Tensor<1,3>
- wedge_product (const Tensor<1,3> (&derivative)[2])
- {
- return cross_product_3d (derivative[0], derivative[1]);
- }
-
-
- // wedge product of dim vectors
- // in dim-d: that's the
- // determinant of the matrix
- template <int dim>
- inline
- Tensor<0,dim>
- wedge_product (const Tensor<1,dim> (&derivative)[dim])
- {
- Tensor<2,dim> jacobian;
- for (unsigned int i=0; i<dim; ++i)
- jacobian[i] = derivative[i];
-
- return determinant (jacobian);
- }
- }
-}
-
-
-template <int dim>
-template <int spacedim>
-void
-GeometryInfo<dim>::
-alternating_form_at_vertices
-#ifndef DEAL_II_CONSTEXPR_BUG
-(const Point<spacedim> (&vertices)[vertices_per_cell],
- Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
-#else
-(const Point<spacedim> *vertices,
- Tensor<spacedim-dim,spacedim> *forms)
-#endif
-{
- // for each of the vertices,
- // compute the alternating form
- // of the mapped unit
- // vectors. consider for
- // example the case of a quad
- // in spacedim==3: to do so, we
- // need to see how the
- // infinitesimal vectors
- // (d\xi_1,0) and (0,d\xi_2)
- // are transformed into
- // spacedim-dimensional space
- // and then form their cross
- // product (i.e. the wedge product
- // of two vectors). to this end, note
- // that
- // \vec x = sum_i \vec v_i phi_i(\vec xi)
- // so the transformed vectors are
- // [x(\xi+(d\xi_1,0))-x(\xi)]/d\xi_1
- // and
- // [x(\xi+(0,d\xi_2))-x(\xi)]/d\xi_2
- // which boils down to the columns
- // of the 3x2 matrix \grad_\xi x(\xi)
- //
- // a similar reasoning would
- // hold for all dim,spacedim
- // pairs -- we only have to
- // compute the wedge product of
- // the columns of the
- // derivatives
- for (unsigned int i=0; i<vertices_per_cell; ++i)
- {
- Tensor<1,spacedim> derivatives[dim];
-
- for (unsigned int j=0; j<vertices_per_cell; ++j)
- {
- const Tensor<1,dim> grad_phi_j
- = d_linear_shape_function_gradient (unit_cell_vertex(i),
- j);
- for (unsigned int l=0; l<dim; ++l)
- derivatives[l] += vertices[j] * grad_phi_j[l];
- }
-
- forms[i] = internal::GeometryInfo::wedge_product (derivatives);
- }
-}
-
-
template struct GeometryInfo<1>;
template struct GeometryInfo<2>;
template struct GeometryInfo<3>;
#endif
;
-
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff