peterrum / dealii

The development repository for the deal.II finite element library.
https://www.dealii.org/
Other
1 stars 0 forks source link

MappingPGeneric #16

Open peterrum opened 4 years ago

peterrum commented 4 years ago

References https://github.com/dealii/dealii/pull/10646.

The class FEValues::reinit() needs for a face (see https://github.com/dealii/dealii/blob/08727cc441b9cfd1081cb677fb75e3b62c8ba060/source/fe/fe_values.cc#L4809-L4835):

drwells commented 4 years ago

Here is the 'original' version of DataSetDescriptor::subface:

template <int dim>
typename QProjector<dim>::DataSetDescriptor
QProjector<dim>::DataSetDescriptor::
subface (const unsigned int face_no,
         const unsigned int subface_no,
         const bool         face_orientation,
         const bool         face_flip,
         const bool         face_rotation,
         const unsigned int n_quadrature_points)
{
  Assert (dim != 1, ExcInternalError());
  Assert (face_no < GeometryInfo<dim>::faces_per_cell,
          ExcInternalError());
                                   // the trick with +1 prevents that we get a
                                   // warning in 1d
  Assert (subface_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
          ExcInternalError());

                   // in 3d, we have to account for faces that
                   // have non-standard face orientation, flip
                   // and rotation. thus, we have to store
                   // _eight_ data sets per face or subface

                   // set up a table with the according offsets
                   // for non-standard orientation, first index:
                   // face_orientation (standard true=1), second
                   // index: face_flip (standard false=0), third
                   // index: face_rotation (standard false=0)
                   //
                   // note, that normally we should use the
                   // obvious offsets 0,1,2,3,4,5,6,7. However,
                   // prior to the changes enabling flipped and
                   // rotated faces, in many places of the
                   // library the convention was used, that the
                   // first dataset with offset 0 corresponds to
                   // a face in standard orientation. therefore
                   // we use the offsets 4,5,6,7,0,1,2,3 here to
                   // stick to that (implicit) convention
  static const unsigned int offset[2][2][2]=
    {{
                       // face_orientation=false; face_flip=false; face_rotation=false and true
      {4*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face,
       5*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face},
                       // face_orientation=false; face_flip=true;  face_rotation=false and true
      {6*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face,
       7*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face}},   
     {
                        // face_orientation=true;  face_flip=false; face_rotation=false and true
       {0*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face,
        1*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face},
                        // face_orientation=true;  face_flip=true;  face_rotation=false and true
       {2*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face,
        3*GeometryInfo<dim>::faces_per_cell*GeometryInfo<dim>::subfaces_per_face}}};

  switch (dim)
    {
      case 1:
      case 2:
            return ((face_no * GeometryInfo<dim>::subfaces_per_face +
                     subface_no)
                    * n_quadrature_points);
      case 3:
            return (((face_no * GeometryInfo<dim>::subfaces_per_face +
                      subface_no)
                     + offset[face_orientation][face_flip][face_rotation]
                     )
                    * n_quadrature_points);
      default:
            Assert (false, ExcInternalError());
    }
  return numbers::invalid_unsigned_int;              
}
drwells commented 4 years ago

See commit b3a7efd

drwells commented 4 years ago
template <int dim>
typename QProjector<dim>::DataSetDescriptor
QProjector<dim>::DataSetDescriptor::
face (const unsigned int face_no,
      const bool         face_orientation,
      const bool         face_flip,
      const bool         face_rotation,
      const unsigned int n_quadrature_points)
{
  Assert (dim != 1, ExcInternalError());
  Assert (face_no < GeometryInfo<dim>::faces_per_cell,
          ExcInternalError());

                   // in 3d, we have to account for faces that
                   // have non-standard face orientation, flip
                   // and rotation. thus, we have to store
                   // _eight_ data sets per face or subface

                   // set up a table with the according offsets
                   // for non-standard orientation, first index:
                   // face_orientation (standard true=1), second
                   // index: face_flip (standard false=0), third
                   // index: face_rotation (standard false=0)
                   //
                   // note, that normally we should use the
                   // obvious offsets 0,1,2,3,4,5,6,7. However,
                   // prior to the changes enabling flipped and
                   // rotated faces, in many places of the
                   // library the convention was used, that the
                   // first dataset with offset 0 corresponds to
                   // a face in standard orientation. therefore
                   // we use the offsets 4,5,6,7,0,1,2,3 here to
                   // stick to that (implicit) convention
  static const unsigned int offset[2][2][2]=
    {{{4*GeometryInfo<dim>::faces_per_cell, 5*GeometryInfo<dim>::faces_per_cell},    // face_orientation=false; face_flip=false; face_rotation=false and true
      {6*GeometryInfo<dim>::faces_per_cell, 7*GeometryInfo<dim>::faces_per_cell}},   // face_orientation=false; face_flip=true;  face_rotation=false and true
     {{0*GeometryInfo<dim>::faces_per_cell, 1*GeometryInfo<dim>::faces_per_cell},    // face_orientation=true;  face_flip=false; face_rotation=false and true
      {2*GeometryInfo<dim>::faces_per_cell, 3*GeometryInfo<dim>::faces_per_cell}}};  // face_orientation=true;  face_flip=true;  face_rotation=false and true

  switch (dim)
    {
      case 1:
      case 2:
            return face_no * n_quadrature_points;

      case 3:
            return ((face_no
             + offset[face_orientation][face_flip][face_rotation])
            * n_quadrature_points);

      default:
            Assert (false, ExcInternalError());
    }
  return numbers::invalid_unsigned_int;
}