[spaces.cg.interface] add local Dirichlet DoFs for simplicial Lagrange elements

dune/gdt/spaces/cg/interface.hh

@@ -52,6 +52,7 @@ public:
   using typename BaseType::PatternType;
 
 private:
+  typedef DSC::FieldVector<DomainFieldType, dimDomain> StuffDomainType;
   static const constexpr RangeFieldType compare_tolerance_ = 1e-13;
 
 public:
@@ -166,6 +167,75 @@ public:
     return localDirichletDofs;
   } // ... local_dirichlet_DoFs_order_1(...)
 
+
+  std::set<size_t> local_dirichlet_DoFs_simplicial_lagrange_elements(const EntityType& entity,
+                                                                     const BoundaryInfoType& boundaryInfo) const
+  {
+    if (!entity.type().isSimplex())
+      DUNE_THROW(NotImplemented, "Only implemented for simplex elements!");
+    // check
+    assert(this->grid_view().indexSet().contains(entity));
+    // prepare
+    std::set<size_t> localDirichletDofs;
+    std::vector<DomainType> dirichlet_vertices;
+    // get all dirichlet vertices of this entity, therefore
+    // * loop over all intersections
+    const auto intersection_it_end = this->grid_view().iend(entity);
+    for (auto intersection_it = this->grid_view().ibegin(entity); intersection_it != intersection_it_end;
+         ++intersection_it) {
+      const auto& intersection = *intersection_it;
+      std::vector<StuffDomainType> dirichlet_vertices_intersection;
+      // only work on dirichlet ones, actual dirichlet intersections + process boundaries for parallel runs
+      if (boundaryInfo.dirichlet(intersection) || (!intersection.neighbor() && !intersection.boundary())) {
+        // and get the vertices of the intersection
+        const auto geometry = intersection.geometry();
+        for (auto cc : DSC::valueRange(geometry.corners())) {
+          dirichlet_vertices_intersection.emplace_back(entity.geometry().local(geometry.corner(cc)));
+          dirichlet_vertices.emplace_back(entity.geometry().local(geometry.corner(cc)));
+        }
+      } // only work on dirichlet ones
+      // for higher polynomial orders, add points associated to DoFs within the intersection
+      // calculate by using lagrange grid {x = \sum_{j=0}^d \lambda_j a_j | \sum_j lambda_j = 1}, where a_j are the
+      // vertices of the entity and \lambda_j \in {\frac{m}{polOrder} | m = 0, ... , polOrder}
+      std::vector<double> possible_coefficients(polOrder < 1 ? 0 : polOrder - 1);
+      for (int m = 0; m < polOrder; ++m)
+        possible_coefficients[m] = m / polOrder;
+      std::set<std::vector<double>> possible_coefficient_vectors;
+      possible_convex_combination_coefficients(
+          possible_coefficient_vectors, possible_coefficients, dirichlet_vertices_intersection.size());
+      for (const auto& coefficient_vector : possible_coefficient_vectors)
+        dirichlet_vertices.emplace_back(std::inner_product(dirichlet_vertices_intersection.begin(),
+                                                           dirichlet_vertices_intersection.end(),
+                                                           coefficient_vector.begin(),
+                                                           StuffDomainType(0)));
+    } // loop over all intersections
+
+    // find the corresponding basis functions
+    const auto basis = this->base_function_set(entity);
+    typedef typename BaseType::BaseFunctionSetType::RangeType RangeType;
+    std::vector<RangeType> tmp_basis_values(basis.size(), RangeType(0));
+    for (size_t cc = 0; cc < dirichlet_vertices.size(); ++cc) {
+      // find the basis function that evaluates to one here (has to be only one per range dimension!)
+      basis.evaluate(dirichlet_vertices[cc], tmp_basis_values);
+      size_t ones     = 0;
+      size_t zeros    = 0;
+      size_t failures = 0;
+      for (size_t jj = 0; jj < basis.size(); ++jj) {
+        for (size_t rr = 0; rr < dimRange; ++rr) {
+          if (std::abs(tmp_basis_values[jj][rr] - RangeFieldType(1)) < compare_tolerance_) {
+            localDirichletDofs.insert(jj);
+            ++ones;
+          } else if (std::abs(tmp_basis_values[jj][rr]) < compare_tolerance_)
+            ++zeros;
+          else
+            ++failures;
+        }
+      }
+      assert(ones == dimRange && zeros == ((basis.size() - 1) * dimRange) && failures == 0 && "This must not happen!");
+    }
+    return localDirichletDofs;
+  } // ... local_dirichlet_DoFs_simplicial_lagrange_elements(...)
+
   using BaseType::compute_pattern;
 
   template <class G, class S, size_t d, size_t r, size_t rC>
@@ -196,6 +266,34 @@ public:
     }
   } // ... local_constraints(..., Constraints::Dirichlet< ... > ...)
   /** @} */
+
+private:
+  void possible_convex_combination_coefficients(std::set<std::vector<double>>& vectors_in,
+                                                const std::vector<double>& possible_coefficients,
+                                                const size_t final_size) const
+  {
+    if (vectors_in.empty()) {
+      if (final_size != 0)
+        for (const auto& coeff : possible_coefficients)
+          vectors_in.insert(std::vector<double>(1, coeff));
+    } else if (vectors_in.begin()->size() != final_size) {
+      std::set<std::vector<double>> vectors_out;
+      for (auto& vec : vectors_in) {
+        for (const auto& coeff : possible_coefficients) {
+          if ((vec.size() != final_size - 1
+               && DSC::FloatCmp::le(std::accumulate(vec.begin(), vec.end(), 0.0) + coeff, 1.0))
+              || (vec.size() == final_size - 1
+                  && DSC::FloatCmp::eq(std::accumulate(vec.begin(), vec.end(), 0.0) + coeff, 1.0))) {
+            std::vector<double> vec_copy = vec;
+            vec_copy.push_back(coeff);
+            vectors_out.insert(vec_copy);
+          }
+        }
+      }
+      vectors_in = vectors_out;
+      possible_convex_combination_coefficients(vectors_in, possible_coefficients, final_size);
+    } // if (...)
+  } // ... possible_convex_combination_coefficients(...)
 }; // class CGInterface