]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add functionality to support multiple mappings per element in FE_PolyTensor
authorgrahambenharper <grahambenharper@gmail.com>
Sun, 7 Jul 2019 03:30:22 +0000 (21:30 -0600)
committergrahambenharper <grahambenharper@gmail.com>
Sun, 7 Jul 2019 03:30:22 +0000 (21:30 -0600)
include/deal.II/fe/fe_dg_vector.templates.h
include/deal.II/fe/fe_poly_tensor.h
source/fe/fe_abf.cc
source/fe/fe_bdm.cc
source/fe/fe_dg_vector.cc
source/fe/fe_nedelec.cc
source/fe/fe_poly_tensor.cc
source/fe/fe_raviart_thomas.cc
source/fe/fe_raviart_thomas_nodal.cc
source/fe/fe_rt_bubbles.cc

index 8aa11a5a29f0cb049b8fe68ae52ecb1d997a2a28..a0893d2e47ef324e6d9eb334d4c58d4b7b51e66b 100644 (file)
@@ -43,7 +43,7 @@ FE_DGVector<PolynomialType, dim, spacedim>::FE_DGVector(const unsigned int deg,
       std::vector<ComponentMask>(PolynomialType::compute_n_pols(deg),
                                  ComponentMask(dim, true)))
 {
-  this->mapping_type                   = map;
+  this->mapping_type                   = {map};
   const unsigned int polynomial_degree = this->tensor_degree();
 
   QGauss<dim> quadrature(polynomial_degree + 1);
index ecff871dc7bbcddbf1e8da3dd15eac611addffac..c755a6f2a6dbf392f25558efda3e7a9e2f45235c 100644 (file)
@@ -207,9 +207,12 @@ public:
 protected:
   /**
    * The mapping type to be used to map shape functions from the reference
-   * cell to the mesh cell.
+   * cell to the mesh cell. If this vector is length one, the same mapping
+   * will be applied to all shape functions. If the vector size is equal to
+   * the finite element dofs per cell, then each shape function will be mapped
+   * according to the corresponding entry in the vector.
    */
-  MappingType mapping_type;
+  std::vector<MappingType> mapping_type;
 
 
   /* NOTE: The following function has its definition inlined into the class
@@ -247,11 +250,30 @@ protected:
     // initialize fields only if really
     // necessary. otherwise, don't
     // allocate memory
+
+    const bool update_transformed_shape_values =
+      (std::find_if(this->mapping_type.begin(),
+                    this->mapping_type.end(),
+                    [](const MappingType t) { return t != mapping_none; }) !=
+       this->mapping_type.end());
+
+    const bool update_transformed_shape_grads =
+      (std::find_if(this->mapping_type.begin(),
+                    this->mapping_type.end(),
+                    [](const MappingType t) {
+                      return (t == mapping_raviart_thomas ||
+                              t == mapping_piola || t == mapping_nedelec ||
+                              t == mapping_contravariant);
+                    }) != this->mapping_type.end());
+
+    const bool update_transformed_shape_hessian_tensors =
+      update_transformed_shape_values;
+
     if (update_flags & update_values)
       {
         values.resize(this->dofs_per_cell);
         data.shape_values.reinit(this->dofs_per_cell, n_q_points);
-        if (mapping_type != mapping_none)
+        if (update_transformed_shape_values)
           data.transformed_shape_values.resize(n_q_points);
       }
 
@@ -261,10 +283,7 @@ protected:
         data.shape_grads.reinit(this->dofs_per_cell, n_q_points);
         data.transformed_shape_grads.resize(n_q_points);
 
-        if ((mapping_type == mapping_raviart_thomas) ||
-            (mapping_type == mapping_piola) ||
-            (mapping_type == mapping_nedelec) ||
-            (mapping_type == mapping_contravariant))
+        if (update_transformed_shape_grads)
           data.untransformed_shape_grads.resize(n_q_points);
       }
 
@@ -273,7 +292,7 @@ protected:
         grad_grads.resize(this->dofs_per_cell);
         data.shape_grad_grads.reinit(this->dofs_per_cell, n_q_points);
         data.transformed_shape_hessians.resize(n_q_points);
-        if (mapping_type != mapping_none)
+        if (update_transformed_shape_hessian_tensors)
           data.untransformed_shape_hessian_tensors.resize(n_q_points);
       }
 
index 9efdc6e053777dffab38404a128641ae158e1b5c..229344046fd913430593398cf3f82f3e01f478bc 100644 (file)
@@ -60,7 +60,7 @@ FE_ABF<dim>::FE_ABF(const unsigned int deg)
   Assert(dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_raviart_thomas;
+  this->mapping_type = {mapping_raviart_thomas};
   // First, initialize the
   // generalized support points and
   // quadrature weights, since they
index 2b5335d940149949a23220d45ee41285008a5ceb..0dfd5374f6773ef2a043220879be1c8e82e73bbf 100644 (file)
@@ -56,7 +56,7 @@ FE_BDM<dim>::FE_BDM(const unsigned int deg)
 
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_bdm;
+  this->mapping_type = {mapping_bdm};
   // These must be done first, since
   // they change the evaluation of
   // basis functions
index 229c009ddb2a2e7060d7536d3f23199d1430fa2b..9409a368cbdcad3dbcde5d72ad890d0d98dae91b 100644 (file)
@@ -27,7 +27,7 @@ DEAL_II_NAMESPACE_OPEN
 
 template <int dim, int spacedim>
 FE_DGNedelec<dim, spacedim>::FE_DGNedelec(const unsigned int p)
-  : FE_DGVector<PolynomialsNedelec<dim>, dim, spacedim>(p, mapping_nedelec)
+  : FE_DGVector<PolynomialsNedelec<dim>, dim, spacedim>(p, {mapping_nedelec})
 {}
 
 
@@ -52,7 +52,7 @@ template <int dim, int spacedim>
 FE_DGRaviartThomas<dim, spacedim>::FE_DGRaviartThomas(const unsigned int p)
   : FE_DGVector<PolynomialsRaviartThomas<dim>, dim, spacedim>(
       p,
-      mapping_raviart_thomas)
+      {mapping_raviart_thomas})
 {}
 
 
@@ -77,7 +77,7 @@ FE_DGRaviartThomas<dim, spacedim>::get_name() const
 
 template <int dim, int spacedim>
 FE_DGBDM<dim, spacedim>::FE_DGBDM(const unsigned int p)
-  : FE_DGVector<PolynomialsBDM<dim>, dim, spacedim>(p, mapping_bdm)
+  : FE_DGVector<PolynomialsBDM<dim>, dim, spacedim>(p, {mapping_bdm})
 {}
 
 
index d12afe8a953b1f84bbe21d21ecaf21c0d9c36037..ad175ac7084c32cec514071f5f977c9c0c42f284 100644 (file)
@@ -87,7 +87,7 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
 
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_nedelec;
+  this->mapping_type = {mapping_nedelec};
   // First, initialize the
   // generalized support points and
   // quadrature weights, since they
index 6b72947c5591afd1010129b1a9254e7455a24674..bc1d0844a30af4c20a9a72c49da169d352f6ea71 100644 (file)
@@ -49,30 +49,29 @@ namespace internal
        * implements an algorithm that determines those DoFs that need this
        * sign change for a given cell.
        */
+      template <int spacedim>
       void
       get_face_sign_change_rt(const dealii::Triangulation<1>::cell_iterator &,
-                              const unsigned int,
-                              std::vector<double> &face_sign)
+                              const FiniteElement<1, spacedim> &,
+                              const std::vector<MappingType> &,
+                              std::vector<double> &)
       {
         // nothing to do in 1d
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
       }
 
 
 
+      //      template<int spacedim>
       void
       get_face_sign_change_rt(
         const dealii::Triangulation<2>::cell_iterator &cell,
-        const unsigned int                             dofs_per_face,
+        const FiniteElement<2, 2> &                    fe,
+        const std::vector<MappingType> &               mapping_type,
         std::vector<double> &                          face_sign)
       {
         const unsigned int dim      = 2;
         const unsigned int spacedim = 2;
 
-        // Default is no sign
-        // change. I.e. multiply by one.
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
-
         for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2;
              f < GeometryInfo<dim>::faces_per_cell;
              ++f)
@@ -84,14 +83,24 @@ namespace internal
                 const unsigned int nn = cell->neighbor_face_no(f);
 
                 if (nn < GeometryInfo<dim>::faces_per_cell / 2)
-                  for (unsigned int j = 0; j < dofs_per_face; ++j)
+                  for (unsigned int j = 0; j < fe.dofs_per_face; ++j)
                     {
-                      Assert(f * dofs_per_face + j < face_sign.size(),
+                      const unsigned int cell_j = fe.face_to_cell_index(j, f);
+                      // something in FiniteElement or FiniteElementData base
+                      // classes
+
+                      Assert(f * fe.dofs_per_face + j < face_sign.size(),
+                             ExcInternalError());
+                      Assert(mapping_type.size() == 1 ||
+                               cell_j < mapping_type.size(),
                              ExcInternalError());
 
                       // TODO: This is probably only going to work for those
                       // elements for which all dofs are face dofs
-                      face_sign[f * dofs_per_face + j] = -1.0;
+                      if ((mapping_type.size() > 1 ?
+                             mapping_type[cell_j] :
+                             mapping_type[0]) == mapping_raviart_thomas)
+                        face_sign[f * fe.dofs_per_face + j] = -1.0;
                     }
               }
           }
@@ -99,51 +108,57 @@ namespace internal
 
 
 
+      template <int spacedim>
       void
       get_face_sign_change_rt(
         const dealii::Triangulation<3>::cell_iterator & /*cell*/,
-        const unsigned int /*dofs_per_face*/,
-        std::vector<double> &face_sign)
+        const FiniteElement<3, spacedim> & /*fe*/,
+        const std::vector<MappingType> & /*mapping_type*/,
+        std::vector<double> & /*face_sign*/)
       {
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
         // TODO: think about what it would take here
       }
 
+
+
+      template <int spacedim>
       void
       get_face_sign_change_nedelec(
         const dealii::Triangulation<1>::cell_iterator & /*cell*/,
-        const unsigned int /*dofs_per_face*/,
-        std::vector<double> &face_sign)
+        const FiniteElement<1, spacedim> & /*fe*/,
+        const std::vector<MappingType> & /*mapping_type*/,
+        std::vector<double> & /*face_sign*/)
       {
         // nothing to do in 1d
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
       }
 
 
 
+      template <int spacedim>
       void
       get_face_sign_change_nedelec(
         const dealii::Triangulation<2>::cell_iterator & /*cell*/,
-        const unsigned int /*dofs_per_face*/,
-        std::vector<double> &face_sign)
+        const FiniteElement<2, spacedim> & /*fe*/,
+        const std::vector<MappingType> & /*mapping_type*/,
+        std::vector<double> & /*face_sign*/)
       {
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
         // TODO: think about what it would take here
       }
 
 
+      template <int spacedim>
       void
       get_face_sign_change_nedelec(
         const dealii::Triangulation<3>::cell_iterator &cell,
-        const unsigned int /*dofs_per_face*/,
-        std::vector<double> &face_sign)
+        const FiniteElement<3, spacedim> & /*fe*/,
+        const std::vector<MappingType> &mapping_type,
+        std::vector<double> &           face_sign)
       {
         const unsigned int dim = 3;
-        std::fill(face_sign.begin(), face_sign.end(), 1.0);
         // TODO: This is probably only going to work for those elements for
         // which all dofs are face dofs
         for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
-          if (!(cell->line_orientation(l)))
+          if (!(cell->line_orientation(l)) & mapping_type[0] == mapping_nedelec)
             face_sign[l] = -1.0;
       }
     } // namespace
@@ -161,7 +176,7 @@ FE_PolyTensor<PolynomialType, dim, spacedim>::FE_PolyTensor(
   : FiniteElement<dim, spacedim>(fe_data,
                                  restriction_is_additive_flags,
                                  nonzero_components)
-  , mapping_type(MappingType::mapping_none)
+  , mapping_type({MappingType::mapping_none})
   , poly_space(PolynomialType(degree))
 {
   cached_point(0) = -1;
@@ -364,18 +379,23 @@ FE_PolyTensor<PolynomialType, dim, spacedim>::fill_fe_values(
   // between two faces.
   std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
 
-  if (mapping_type == mapping_raviart_thomas)
-    internal::FE_PolyTensor::get_face_sign_change_rt(cell,
-                                                     this->dofs_per_face,
-                                                     fe_data.sign_change);
-  else if (mapping_type == mapping_nedelec)
-    internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
-                                                          this->dofs_per_face,
-                                                          fe_data.sign_change);
+  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
+                                                   *this,
+                                                   this->mapping_type,
+                                                   fe_data.sign_change);
+
+  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
+                                                        *this,
+                                                        this->mapping_type,
+                                                        fe_data.sign_change);
 
 
   for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
     {
+      MappingType mapping_type =
+        (this->mapping_type.size() > 1 ? this->mapping_type[i] :
+                                         this->mapping_type[0]);
+
       const unsigned int first =
         output_data.shape_function_to_row_table[i * this->n_components() +
                                                 this->get_nonzero_components(i)
@@ -940,18 +960,22 @@ FE_PolyTensor<PolynomialType, dim, spacedim>::fill_fe_face_values(
   // on the neighborhood between two faces.
   std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
 
-  if (mapping_type == mapping_raviart_thomas)
-    internal::FE_PolyTensor::get_face_sign_change_rt(cell,
-                                                     this->dofs_per_face,
-                                                     fe_data.sign_change);
+  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
+                                                   *this,
+                                                   this->mapping_type,
+                                                   fe_data.sign_change);
 
-  else if (mapping_type == mapping_nedelec)
-    internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
-                                                          this->dofs_per_face,
-                                                          fe_data.sign_change);
+  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
+                                                        *this,
+                                                        this->mapping_type,
+                                                        fe_data.sign_change);
 
   for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
     {
+      const MappingType mapping_type =
+        (this->mapping_type.size() > 1 ? this->mapping_type[i] :
+                                         this->mapping_type[0]);
+
       const unsigned int first =
         output_data.shape_function_to_row_table[i * this->n_components() +
                                                 this->get_nonzero_components(i)
@@ -1568,18 +1592,22 @@ FE_PolyTensor<PolynomialType, dim, spacedim>::fill_fe_subface_values(
   // on the neighborhood between two faces.
   std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0);
 
-  if (mapping_type == mapping_raviart_thomas)
-    internal::FE_PolyTensor::get_face_sign_change_rt(cell,
-                                                     this->dofs_per_face,
-                                                     fe_data.sign_change);
+  internal::FE_PolyTensor::get_face_sign_change_rt(cell,
+                                                   *this,
+                                                   this->mapping_type,
+                                                   fe_data.sign_change);
 
-  else if (mapping_type == mapping_nedelec)
-    internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
-                                                          this->dofs_per_face,
-                                                          fe_data.sign_change);
+  internal::FE_PolyTensor::get_face_sign_change_nedelec(cell,
+                                                        *this,
+                                                        this->mapping_type,
+                                                        fe_data.sign_change);
 
   for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
     {
+      MappingType mapping_type =
+        (this->mapping_type.size() > 1 ? this->mapping_type[i] :
+                                         this->mapping_type[0]);
+
       const unsigned int first =
         output_data.shape_function_to_row_table[i * this->n_components() +
                                                 this->get_nonzero_components(i)
@@ -2140,88 +2168,95 @@ FE_PolyTensor<PolynomialType, dim, spacedim>::requires_update_flags(
 {
   UpdateFlags out = update_default;
 
-  switch (mapping_type)
+  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
     {
-      case mapping_none:
-        {
-          if (flags & update_values)
-            out |= update_values;
-
-          if (flags & update_gradients)
-            out |= update_gradients | update_values |
-                   update_jacobian_pushed_forward_grads;
-
-          if (flags & update_hessians)
-            out |= update_hessians | update_values | update_gradients |
-                   update_jacobian_pushed_forward_grads |
-                   update_jacobian_pushed_forward_2nd_derivatives;
-          break;
-        }
+      MappingType mapping_type =
+        (this->mapping_type.size() > 1 ? this->mapping_type[i] :
+                                         this->mapping_type[0]);
 
-      case mapping_raviart_thomas:
-      case mapping_piola:
+      switch (mapping_type)
         {
-          if (flags & update_values)
-            out |= update_values | update_piola;
-
-          if (flags & update_gradients)
-            out |= update_gradients | update_values | update_piola |
-                   update_jacobian_pushed_forward_grads |
-                   update_covariant_transformation |
-                   update_contravariant_transformation;
-
-          if (flags & update_hessians)
-            out |= update_hessians | update_piola | update_values |
-                   update_gradients | update_jacobian_pushed_forward_grads |
-                   update_jacobian_pushed_forward_2nd_derivatives |
-                   update_covariant_transformation;
-
-          break;
-        }
+          case mapping_none:
+            {
+              if (flags & update_values)
+                out |= update_values;
+
+              if (flags & update_gradients)
+                out |= update_gradients | update_values |
+                       update_jacobian_pushed_forward_grads;
+
+              if (flags & update_hessians)
+                out |= update_hessians | update_values | update_gradients |
+                       update_jacobian_pushed_forward_grads |
+                       update_jacobian_pushed_forward_2nd_derivatives;
+              break;
+            }
+          case mapping_raviart_thomas:
+          case mapping_piola:
+            {
+              if (flags & update_values)
+                out |= update_values | update_piola;
+
+              if (flags & update_gradients)
+                out |= update_gradients | update_values | update_piola |
+                       update_jacobian_pushed_forward_grads |
+                       update_covariant_transformation |
+                       update_contravariant_transformation;
+
+              if (flags & update_hessians)
+                out |= update_hessians | update_piola | update_values |
+                       update_gradients | update_jacobian_pushed_forward_grads |
+                       update_jacobian_pushed_forward_2nd_derivatives |
+                       update_covariant_transformation;
+
+              break;
+            }
 
-      case mapping_contravariant:
-        {
-          if (flags & update_values)
-            out |= update_values | update_piola;
-
-          if (flags & update_gradients)
-            out |= update_gradients | update_values |
-                   update_jacobian_pushed_forward_grads |
-                   update_covariant_transformation |
-                   update_contravariant_transformation;
-
-          if (flags & update_hessians)
-            out |= update_hessians | update_piola | update_values |
-                   update_gradients | update_jacobian_pushed_forward_grads |
-                   update_jacobian_pushed_forward_2nd_derivatives |
-                   update_covariant_transformation;
-
-          break;
-        }
 
-      case mapping_nedelec:
-      case mapping_covariant:
-        {
-          if (flags & update_values)
-            out |= update_values | update_covariant_transformation;
+          case mapping_contravariant:
+            {
+              if (flags & update_values)
+                out |= update_values | update_piola;
+
+              if (flags & update_gradients)
+                out |= update_gradients | update_values |
+                       update_jacobian_pushed_forward_grads |
+                       update_covariant_transformation |
+                       update_contravariant_transformation;
+
+              if (flags & update_hessians)
+                out |= update_hessians | update_piola | update_values |
+                       update_gradients | update_jacobian_pushed_forward_grads |
+                       update_jacobian_pushed_forward_2nd_derivatives |
+                       update_covariant_transformation;
+
+              break;
+            }
 
-          if (flags & update_gradients)
-            out |= update_gradients | update_values |
-                   update_jacobian_pushed_forward_grads |
-                   update_covariant_transformation;
+          case mapping_nedelec:
+          case mapping_covariant:
+            {
+              if (flags & update_values)
+                out |= update_values | update_covariant_transformation;
 
-          if (flags & update_hessians)
-            out |= update_hessians | update_values | update_gradients |
-                   update_jacobian_pushed_forward_grads |
-                   update_jacobian_pushed_forward_2nd_derivatives |
-                   update_covariant_transformation;
+              if (flags & update_gradients)
+                out |= update_gradients | update_values |
+                       update_jacobian_pushed_forward_grads |
+                       update_covariant_transformation;
 
-          break;
-        }
+              if (flags & update_hessians)
+                out |= update_hessians | update_values | update_gradients |
+                       update_jacobian_pushed_forward_grads |
+                       update_jacobian_pushed_forward_2nd_derivatives |
+                       update_covariant_transformation;
 
-      default:
-        {
-          Assert(false, ExcNotImplemented());
+              break;
+            }
+
+          default:
+            {
+              Assert(false, ExcNotImplemented());
+            }
         }
     }
 
index 779623ee622af92d97fb0611f1ad989db3db20c5..0aead6d54e1ed104dca5d2e152ecb536dc96696a 100644 (file)
@@ -60,7 +60,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas(const unsigned int deg)
   Assert(dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_raviart_thomas;
+  this->mapping_type = {mapping_raviart_thomas};
   // First, initialize the
   // generalized support points and
   // quadrature weights, since they
index 736815808003ef610cf0db05753d7939ab844091..a7f40dfc9847a16b9c96bef0f36775d1da68c4a9 100644 (file)
@@ -52,7 +52,7 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int deg)
   Assert(dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_raviart_thomas;
+  this->mapping_type = {mapping_raviart_thomas};
   // First, initialize the
   // generalized support points and
   // quadrature weights, since they
index 240e9786988264752e71704eaccad667cfb61f50..a5ab697f275d70820c485686b6545d714c3f625a 100644 (file)
@@ -55,7 +55,7 @@ FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
       "Lowest order RT_Bubbles element is degree 1, but you requested for degree 0"));
   const unsigned int n_dofs = this->dofs_per_cell;
 
-  this->mapping_type = mapping_raviart_thomas;
+  this->mapping_type = {mapping_raviart_thomas};
   // Initialize support points and quadrature weights
   initialize_support_points(deg);
   // Compute the inverse node matrix to get

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.