diff --git a/dune/gdt/spaces/cg/interface.hh b/dune/gdt/spaces/cg/interface.hh index b7a830ae568940f7cf7c7a5c18a03741ab5fc4e9..02bc1f3bf0b11a2545302b02c7823485007b9f9a 100644 --- a/dune/gdt/spaces/cg/interface.hh +++ b/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