{"action":"created","comment":{"url":"https://api.github.com/repos/dealii/dealii/pulls/comments/80804936","pull_request_review_id":1834319,"id":80804936,"diff_hunk":"@@ -0,0 +1,927 @@\n+// ---------------------------------------------------------------------\n+//\n+// Copyright (C) 2016 by the deal.II authors\n+//\n+// This file is part of the deal.II library.\n+//\n+// The deal.II library is free software; you can use it, redistribute\n+// it, and/or modify it under the terms of the GNU Lesser General\n+// Public License as published by the Free Software Foundation; either\n+// version 2.1 of the License, or (at your option) any later version.\n+// The full text of the license can be found in the file LICENSE at\n+// the top level of the deal.II distribution.\n+//\n+// ---------------------------------------------------------------------\n+#include <deal.II/fe/fe_enriched.h>\n+#ifdef DEAL_II_WITH_CXX11\n+\n+#include <deal.II/fe/fe_tools.h>\n+\n+DEAL_II_NAMESPACE_OPEN\n+\n+namespace\n+{\n+  /**\n+   * Auxiliary function to create multiplicity vector from input enrichment functions.\n+   */\n+  template <typename T>\n+  std::vector<unsigned int>\n+  build_multiplicities(const std::vector<std::vector<T > > &functions )\n+  {\n+    std::vector<unsigned int> multiplicities;\n+    multiplicities.push_back(1); // the first one is non-enriched FE\n+    for (unsigned int i = 0; i < functions.size(); i++)\n+      multiplicities.push_back(functions[i].size());\n+\n+    return multiplicities;\n+  }\n+\n+\n+  /**\n+   * Auxiliary function to build FiniteElement's vector\n+   */\n+  template <int dim, int spacedim>\n+  std::vector< const FiniteElement< dim, spacedim > * >\n+  build_fes(const FiniteElement<dim,spacedim> *fe_base,\n+            const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched)\n+  {\n+    std::vector< const FiniteElement< dim, spacedim > * > fes;\n+    fes.push_back(fe_base);\n+    for (unsigned int i = 0; i < fe_enriched.size(); i++)\n+      fes.push_back(fe_enriched[i]);\n+\n+    return fes;\n+  }\n+\n+\n+  /**\n+   * Auxiliary function which check consistency of the input parameters.\n+   * Returns true if everything is ok.\n+   */\n+  template <int dim, int spacedim>\n+  bool\n+  consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes,\n+                     const std::vector< unsigned int > &multiplicities,\n+                     const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)\n+  {\n+    AssertThrow(fes.size() > 0,\n+                ExcMessage(\"FEs size should be >=1\"));\n+    AssertThrow(fes.size() == multiplicities.size(),\n+                ExcMessage(\"FEs and multiplicities should have the same size\"));\n+\n+    AssertThrow (functions.size() == fes.size() - 1,\n+                 ExcDimensionMismatch(functions.size(), fes.size()-1));\n+\n+    AssertThrow(multiplicities[0] == 1,\n+                ExcMessage(\"First multiplicity should be 1\"));\n+\n+    const unsigned int n_comp_base = fes[0]->n_components();\n+\n+    // start from fe=1 as 0th is always non-enriched FE.\n+    for (unsigned int fe=1; fe < fes.size(); fe++)\n+      {\n+        const FE_Nothing<dim> *fe_nothing = dynamic_cast<const FE_Nothing<dim>*>(fes[fe]);\n+        if (fe_nothing)\n+          AssertThrow (fe_nothing->is_dominating(),\n+                       ExcMessage(\"Only dominating FE_Nothing can be used in FE_Enriched\"));\n+\n+        AssertThrow (fes[fe]->n_components() == n_comp_base,\n+                     ExcMessage(\"All elements must have the same number of components\"));\n+      }\n+    return true;\n+  }\n+\n+\n+  /**\n+   * Auxiliary function which determines whether the FiniteElement will be enriched.\n+   */\n+  template <int dim, int spacedim>\n+  bool\n+  check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes)\n+  {\n+    // start from fe=1 as 0th is always non-enriched FE.\n+    for (unsigned int fe=1; fe < fes.size(); fe++)\n+      if (dynamic_cast<const FE_Nothing<dim>*>(fes[fe]) == 0)\n+        // this is not FE_Nothing => there will be enrichment\n+        return true;\n+\n+    return false;\n+  }\n+}\n+\n+\n+template <int dim, int spacedim>\n+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> &fe_base)\n+  :\n+  FE_Enriched<dim,spacedim>(fe_base,\n+                            FE_Nothing<dim,spacedim>(fe_base.n_components(),true),\n+                            NULL)\n+{\n+}\n+\n+\n+template <int dim, int spacedim>\n+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> &fe_base,\n+                                        const FiniteElement<dim,spacedim> &fe_enriched,\n+                                        const Function<spacedim>      *enrichment_function)\n+  :\n+  FE_Enriched<dim,spacedim>\n+  (&fe_base, { &fe_enriched },\n+{\n+  {\n+    [=] (const typename Triangulation<dim, spacedim>::cell_iterator &) -> const Function<spacedim> *\n+    {\n+      return enrichment_function;\n+    }\n+  }\n+})\n+{\n+}\n+\n+\n+template <int dim, int spacedim>\n+FE_Enriched<dim,spacedim>::FE_Enriched (const FiniteElement<dim,spacedim> *fe_base,\n+                                        const std::vector<const FiniteElement<dim,spacedim> * > &fe_enriched,\n+                                        const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)\n+  :\n+  FE_Enriched<dim,spacedim> (build_fes(fe_base,fe_enriched),\n+                             build_multiplicities(functions),\n+                             functions)\n+{\n+}\n+\n+\n+template <int dim, int spacedim>\n+FE_Enriched<dim,spacedim>::FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes,\n+                                        const std::vector< unsigned int > &multiplicities,\n+                                        const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > > &functions)\n+  :\n+  FiniteElement<dim,spacedim> (FETools::Compositing::multiply_dof_numbers(fes,multiplicities,false),\n+                               FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities),\n+                               FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)),\n+  enrichments(functions),\n+  is_enriched(check_if_enriched(fes)),\n+  fe_system(fes,multiplicities)\n+{\n+  // descriptive error are thrown within the function.\n+  Assert(consistency_check(fes,multiplicities,functions),\n+         ExcInternalError());\n+\n+  initialize(fes, multiplicities);\n+\n+  // resize to be consistent with all FEs used to construct the FE_Enriched,\n+  // even though we will never use the 0th element.\n+  base_no_mult_local_enriched_dofs.resize(fes.size());\n+  for (unsigned int fe=1; fe < fes.size(); fe++)\n+    base_no_mult_local_enriched_dofs[fe].resize(multiplicities[fe]);\n+\n+  Assert (base_no_mult_local_enriched_dofs.size() == this->n_base_elements(),\n+          ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),\n+                               this->n_base_elements()));\n+\n+  // build the map: (base_no, base_m) -> vector of local element DoFs\n+  for (unsigned int system_index=0; system_index<this->dofs_per_cell;\n+       ++system_index)\n+    {\n+      const unsigned int base_no = this->system_to_base_table[system_index].first.first;\n+      if (base_no == 0) // 0th is always non-enriched FE\n+        continue;\n+\n+      const unsigned int base_m  = this->system_to_base_table[system_index].first.second;\n+      const unsigned int\n+      base_index = this->system_to_base_table[system_index].second;\n+\n+      Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(),\n+              ExcMessage(\"Size mismatch for base_no_mult_local_enriched_dofs: \"\n+                         \"base_index = \" + std::to_string(base_index) +\n+                         \"; base_no = \" + std::to_string(base_no) +\n+                         \"; base_m = \" + std::to_string(base_m) +\n+                         \"; system_index = \" + std::to_string(system_index)));\n+\n+      Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(),\n+              ExcDimensionMismatch(base_m,\n+                                   base_no_mult_local_enriched_dofs[base_no].size()));\n+\n+      base_no_mult_local_enriched_dofs[base_no][base_m].push_back(system_index);\n+    }\n+\n+  // make sure that local_enriched_dofs.size() is correct, that is equals to DoFs\n+  // per cell of the corresponding FE.\n+  for (unsigned int base_no = 1; base_no < base_no_mult_local_enriched_dofs.size(); base_no++)\n+    {\n+      for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)\n+        Assert ( base_no_mult_local_enriched_dofs[base_no][m].size() == fes[base_no]->dofs_per_cell,\n+                 ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no][m].size(),\n+                                      fes[base_no]->dofs_per_cell));\n+    }\n+}\n+\n+\n+template <int dim, int spacedim>\n+const std::vector<std::vector<std::function<const Function<spacedim> *(const typename Triangulation<dim, spacedim>::cell_iterator &) > > >\n+FE_Enriched<dim,spacedim>::get_enrichments() const\n+{\n+  return enrichments;\n+}\n+\n+\n+template <int dim, int spacedim>\n+double\n+FE_Enriched<dim,spacedim>::shape_value(const unsigned int   i,\n+                                       const Point< dim > &p) const\n+{\n+  Assert(!is_enriched, ExcMessage(\"For enriched finite elements shape_value() can not be defined on the reference element.\"));\n+  return fe_system.shape_value(i,p);\n+}\n+\n+\n+template <int dim, int spacedim>\n+FiniteElement<dim,spacedim> *\n+FE_Enriched<dim,spacedim>::clone() const\n+{\n+  std::vector< const FiniteElement< dim, spacedim > * > fes;\n+  std::vector< unsigned int > multiplicities;\n+\n+  for (unsigned int i=0; i<this->n_base_elements(); i++)\n+    {\n+      fes.push_back( & base_element(i) );\n+      multiplicities.push_back(this->element_multiplicity(i) );\n+    }\n+\n+  return new FE_Enriched<dim,spacedim>(fes, multiplicities, get_enrichments());\n+}\n+\n+\n+template <int dim, int spacedim>\n+UpdateFlags\n+FE_Enriched<dim,spacedim>::requires_update_flags (const UpdateFlags flags) const\n+{\n+  UpdateFlags out = fe_system.requires_update_flags(flags);\n+\n+  if (is_enriched)\n+    {\n+      // if we ask for values or gradients, then we would need quadrature points\n+      if (flags & (update_values | update_gradients))\n+        out |= update_quadrature_points;\n+\n+      // if need gradients, add update_values due to product rule\n+      if (out & update_gradients)\n+        out |= update_values;\n+    }\n+\n+  Assert (!(flags & update_3rd_derivatives),\n+          ExcNotImplemented());\n+\n+  return out;\n+}\n+\n+\n+template <int dim, int spacedim>\n+template <int dim_1>\n+typename FiniteElement<dim,spacedim>::InternalDataBase *\n+FE_Enriched<dim,spacedim>::setup_data (std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase> fes_data,\n+                                       const UpdateFlags      flags,\n+                                       const Quadrature<dim_1> &quadrature) const\n+{\n+  Assert ((dynamic_cast<typename FESystem<dim,spacedim>::InternalData *> (fes_data.get()) != NULL),\n+          ExcInternalError());\n+  typename FESystem<dim,spacedim>::InternalData *data_fesystem =\n+    static_cast<typename FESystem<dim,spacedim>::InternalData *> (fes_data.get());\n+\n+  // FESystem::InternalData will be aggregated (owned) by\n+  // our InternalData.\n+  fes_data.release();\n+  InternalData *data = new InternalData(std::unique_ptr<typename FESystem<dim,spacedim>::InternalData>(data_fesystem));\n+\n+  // copy update_each from FESystem data:\n+  data->update_each = data_fesystem->update_each;\n+\n+  // resize cache array according to requested flags\n+  data->enrichment.resize(this->n_base_elements());\n+\n+  const unsigned int n_q_points = quadrature.size();\n+\n+  for (unsigned int base=0; base < this->n_base_elements(); ++base)\n+    {\n+      data->enrichment[base].resize(this->element_multiplicity(base));\n+      for (unsigned int m = 0; m < this->element_multiplicity(base); ++m)\n+        {\n+          if (flags & update_values)\n+            data->enrichment[base][m].values.resize(n_q_points);\n+\n+          if (flags & update_gradients)\n+            data->enrichment[base][m].gradients.resize(n_q_points);\n+\n+          if (flags & update_hessians)\n+            data->enrichment[base][m].hessians.resize(n_q_points);\n+        }\n+    }\n+\n+  return data;\n+}\n+\n+\n+template <int dim, int spacedim>\n+typename FiniteElement<dim,spacedim>::InternalDataBase *\n+FE_Enriched<dim,spacedim>::get_face_data (const UpdateFlags      update_flags,\n+                                          const Mapping<dim,spacedim>    &mapping,\n+                                          const Quadrature<dim-1> &quadrature,\n+                                          internal::FEValues::FiniteElementRelatedData< dim, spacedim >        &output_data) const\n+{\n+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_face_data(update_flags,mapping,quadrature,output_data)),\n+                    update_flags,\n+                    quadrature);\n+}\n+\n+\n+template <int dim, int spacedim>\n+typename FiniteElement<dim,spacedim>::InternalDataBase *\n+FE_Enriched<dim,spacedim>::get_subface_data (const UpdateFlags      update_flags,\n+                                             const Mapping<dim,spacedim>    &mapping,\n+                                             const Quadrature<dim-1> &quadrature,\n+                                             dealii::internal::FEValues::FiniteElementRelatedData<dim, spacedim> &output_data) const\n+{\n+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_subface_data(update_flags,mapping,quadrature,output_data)),\n+                    update_flags,\n+                    quadrature);\n+}\n+\n+\n+template <int dim, int spacedim>\n+typename FiniteElement<dim,spacedim>::InternalDataBase *\n+FE_Enriched<dim,spacedim>::get_data (const UpdateFlags      flags,\n+                                     const Mapping<dim,spacedim>    &mapping,\n+                                     const Quadrature<dim> &quadrature,\n+                                     internal::FEValues::FiniteElementRelatedData< dim, spacedim >   &output_data) const\n+{\n+  return setup_data(std::unique_ptr<typename FiniteElement<dim,spacedim>::InternalDataBase>(fe_system.get_data(flags,mapping,quadrature,output_data)),\n+                    flags,\n+                    quadrature);\n+}\n+\n+\n+template <int dim, int spacedim>\n+void FE_Enriched<dim,spacedim>::initialize (const std::vector<const FiniteElement<dim,spacedim>*> &fes,\n+                                            const std::vector<unsigned int> &multiplicities)\n+{\n+  Assert (fes.size() == multiplicities.size(),\n+          ExcDimensionMismatch (fes.size(), multiplicities.size()) );\n+\n+  // Note that we need to skip every fe with multiplicity 0 in the following block of code\n+  this->base_to_block_indices.reinit(0, 0);\n+\n+  for (unsigned int i=0; i<fes.size(); i++)\n+    if (multiplicities[i]>0)\n+      this->base_to_block_indices.push_back( multiplicities[i] );\n+\n+  {\n+    // If the system is not primitive, these have not been initialized by\n+    // FiniteElement\n+    this->system_to_component_table.resize(this->dofs_per_cell);\n+    this->face_system_to_component_table.resize(this->dofs_per_face);\n+\n+    FETools::Compositing::build_cell_tables(this->system_to_base_table,\n+                                            this->system_to_component_table,\n+                                            this->component_to_base_table,\n+                                            *this,\n+                                            false);\n+\n+    FETools::Compositing::build_face_tables(this->face_system_to_base_table,\n+                                            this->face_system_to_component_table,\n+                                            *this,\n+                                            false);\n+\n+  }\n+\n+  // restriction and prolongation matrices are built on demand\n+\n+  // now set up the interface constraints for h-refinement.\n+  // take them from fe_system:\n+  this->interface_constraints = fe_system.interface_constraints;\n+\n+  // if we just wrap another FE (i.e. use FE_Nothing as a second FE)\n+  // then it makes sense to have support points.\n+  // However, functions like interpolate_boundary_values() need all FEs inside\n+  // FECollection to be able to provide support points irrespectively whether\n+  // this FE sits on the boundary or not. Thus for moment just copy support\n+  // points from fe system:\n+  {\n+    this->unit_support_points      = fe_system.unit_support_points;\n+    this->unit_face_support_points = fe_system.unit_face_support_points;\n+  }\n+\n+  // take adjust_quad_dof_index_for_face_orientation_table from FESystem:\n+  {\n+    this->adjust_line_dof_index_for_line_orientation_table = fe_system.adjust_line_dof_index_for_line_orientation_table;\n+  }\n+}\n+\n+\n+template <int dim, int spacedim>\n+std::string\n+FE_Enriched<dim,spacedim>::get_name () const\n+{\n+  std::ostringstream namebuf;\n+\n+  namebuf << \"FE_Enriched<\"\n+          << Utilities::dim_string(dim,spacedim)\n+          << \">[\";\n+  for (unsigned int i=0; i< this->n_base_elements(); ++i)\n+    {\n+      namebuf << base_element(i).get_name();\n+      if (this->element_multiplicity(i) != 1)\n+        namebuf << '^' << this->element_multiplicity(i);\n+      if (i != this->n_base_elements()-1)\n+        namebuf << '-';\n+    }\n+  namebuf << ']';\n+\n+  return namebuf.str();\n+}\n+\n+\n+template <int dim, int spacedim>\n+const FiniteElement<dim,spacedim> &\n+FE_Enriched<dim,spacedim>::base_element (const unsigned int index) const\n+{\n+  return fe_system.base_element(index);\n+}\n+\n+\n+template <int dim, int spacedim>\n+void\n+FE_Enriched<dim,spacedim>::fill_fe_values (const typename Triangulation< dim, spacedim >::cell_iterator &cell,\n+                                           const CellSimilarity::Similarity cell_similarity,\n+                                           const Quadrature< dim > &quadrature,\n+                                           const Mapping< dim, spacedim > &mapping,\n+                                           const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,\n+                                           const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,\n+                                           const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,\n+                                           internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data\n+                                          ) const\n+{\n+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,\n+          ExcInternalError());\n+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);\n+\n+  // call FESystem's method to fill everything without enrichment function\n+  fe_system.fill_fe_values(cell,\n+                           cell_similarity,\n+                           quadrature,\n+                           mapping,\n+                           mapping_internal,\n+                           mapping_data,\n+                           *fe_data.fesystem_data,\n+                           output_data);\n+\n+  if (is_enriched)\n+    multiply_by_enrichment(quadrature,\n+                           fe_data,\n+                           mapping_data,\n+                           cell,\n+                           output_data);\n+}\n+\n+\n+template <int dim, int spacedim>\n+void\n+FE_Enriched<dim,spacedim>::fill_fe_face_values\n+(const typename Triangulation< dim, spacedim >::cell_iterator &cell,\n+ const unsigned int face_no,\n+ const Quadrature< dim-1 > &quadrature,\n+ const Mapping< dim, spacedim > &mapping,\n+ const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,\n+ const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,\n+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,\n+ internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data\n+) const\n+{\n+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,\n+          ExcInternalError());\n+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);\n+\n+  // call FESystem's method to fill everything without enrichment function\n+  fe_system.fill_fe_face_values(cell,\n+                                face_no,\n+                                quadrature,\n+                                mapping,\n+                                mapping_internal,\n+                                mapping_data,\n+                                *fe_data.fesystem_data,\n+                                output_data);\n+\n+  if (is_enriched)\n+    multiply_by_enrichment(quadrature,\n+                           fe_data,\n+                           mapping_data,\n+                           cell,\n+                           output_data);\n+}\n+\n+\n+template <int dim, int spacedim>\n+void\n+FE_Enriched<dim,spacedim>::fill_fe_subface_values\n+(const typename Triangulation< dim, spacedim >::cell_iterator &cell,\n+ const unsigned int face_no,\n+ const unsigned int sub_no,\n+ const Quadrature< dim-1 > &quadrature,\n+ const Mapping< dim, spacedim > &mapping,\n+ const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal,\n+ const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data,\n+ const typename FiniteElement<dim,spacedim>::InternalDataBase &fe_internal,\n+ internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data\n+) const\n+{\n+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != NULL,\n+          ExcInternalError());\n+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);\n+\n+  // call FESystem's method to fill everything without enrichment function\n+  fe_system.fill_fe_subface_values(cell,\n+                                   face_no,\n+                                   sub_no,\n+                                   quadrature,\n+                                   mapping,\n+                                   mapping_internal,\n+                                   mapping_data,\n+                                   *fe_data.fesystem_data,\n+                                   output_data);\n+\n+  if (is_enriched)\n+    multiply_by_enrichment(quadrature,\n+                           fe_data,\n+                           mapping_data,\n+                           cell,\n+                           output_data);\n+\n+}\n+\n+\n+template <int dim, int spacedim>\n+template <int dim_1>\n+void\n+FE_Enriched<dim,spacedim>::multiply_by_enrichment\n+(const Quadrature<dim_1> &quadrature,\n+ const InternalData &fe_data,\n+ const internal::FEValues::MappingRelatedData<dim,spacedim> &mapping_data,\n+ const typename Triangulation< dim, spacedim >::cell_iterator &cell,\n+ internal::FEValues::FiniteElementRelatedData<dim,spacedim> &output_data) const\n+{\n+  // mapping_data will contain quadrature points on the real element.\n+  // fe_internal is needed to get update flags\n+  // finally, output_data should store all the results we need.\n+\n+  // Either dim_1==dim\n+  // (fill_fe_values) or dim_1==dim-1\n+  // (fill_fe_(sub)face_values)\n+  Assert(dim_1==dim || dim_1==dim-1, ExcInternalError());\n+  const UpdateFlags flags = fe_data.update_each;\n+//  const UpdateFlags flags(dim_1==dim ?\n+//                          fe_internal.current_update_flags() :\n+//                          fe_internal.update_flags);\n+\n+  const unsigned int n_q_points  = quadrature.size();\n+\n+  // First, populate output_data object (that shall hold everything requested such\n+  // as shape value, gradients, hessians, etc) from each base element.\n+  // That is almost identical to FESystem::compute_fill_one_base(),\n+  // the difference being that we do it irrespectively of cell_similarity\n+  // and use base_fe_data.update_flags\n+\n+  // TODO: do we need it only for dim_1 == dim (i.e. fill_fe_values)?\n+  if (dim_1 == dim)\n+    for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)\n+      {\n+        const FiniteElement<dim,spacedim> &\n+        base_fe      = base_element(base_no);\n+        typename FiniteElement<dim,spacedim>::InternalDataBase &\n+        base_fe_data = fe_data.get_fe_data(base_no);\n+        internal::FEValues::FiniteElementRelatedData<dim,spacedim> &\n+        base_data    = fe_data.get_fe_output_object(base_no);\n+\n+        const UpdateFlags base_flags = base_fe_data.update_each;\n+\n+        for (unsigned int system_index=0; system_index<this->dofs_per_cell;\n+             ++system_index)\n+          if (this->system_to_base_table[system_index].first.first == base_no)\n+            {\n+              const unsigned int\n+              base_index = this->system_to_base_table[system_index].second;\n+              Assert (base_index<base_fe.dofs_per_cell, ExcInternalError());\n+\n+              // now copy. if the shape function is primitive, then there\n+              // is only one value to be copied, but for non-primitive\n+              // elements, there might be more values to be copied\n+              //\n+              // so, find out from which index to take this one value, and\n+              // to which index to put\n+              unsigned int out_index = 0;\n+              for (unsigned int i=0; i<system_index; ++i)\n+                out_index += this->n_nonzero_components(i);\n+              unsigned int in_index = 0;\n+              for (unsigned int i=0; i<base_index; ++i)\n+                in_index += base_fe.n_nonzero_components(i);\n+\n+              // then loop over the number of components to be copied\n+              Assert (this->n_nonzero_components(system_index) ==\n+                      base_fe.n_nonzero_components(base_index),\n+                      ExcInternalError());\n+              for (unsigned int s=0; s<this->n_nonzero_components(system_index); ++s)\n+                {\n+                  if (base_flags & update_values)\n+                    for (unsigned int q=0; q<n_q_points; ++q)\n+                      output_data.shape_values[out_index+s][q] =\n+                        base_data.shape_values(in_index+s,q);\n+\n+                  if (base_flags & update_gradients)\n+                    for (unsigned int q=0; q<n_q_points; ++q)\n+                      output_data.shape_gradients[out_index+s][q] =\n+                        base_data.shape_gradients[in_index+s][q];\n+\n+                  if (base_flags & update_hessians)\n+                    for (unsigned int q=0; q<n_q_points; ++q)\n+                      output_data.shape_hessians[out_index+s][q] =\n+                        base_data.shape_hessians[in_index+s][q];\n+\n+                }\n+            }\n+      }\n+\n+  Assert (base_no_mult_local_enriched_dofs.size() == fe_data.enrichment.size(),\n+          ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(),\n+                               fe_data.enrichment.size()));\n+  // calculate hessians, gradients and values for each function\n+  for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)\n+    {\n+      Assert (base_no_mult_local_enriched_dofs[base_no].size() == fe_data.enrichment[base_no].size(),\n+              ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no].size(),\n+                                   fe_data.enrichment[base_no].size()));\n+      for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)\n+        {\n+          if (flags & update_hessians)\n+            {\n+              Assert (fe_data.enrichment[base_no][m].hessians.size() == n_q_points,\n+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].hessians.size(),\n+                                           n_q_points));\n+              for (unsigned int q=0; q<n_q_points; q++)\n+                fe_data.enrichment[base_no][m].hessians[q]  = enrichments[base_no-1][m](cell)->hessian (mapping_data.quadrature_points[q]);\n+            }\n+\n+          if (flags & update_gradients)\n+            {\n+              Assert (fe_data.enrichment[base_no][m].gradients.size() == n_q_points,\n+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].gradients.size(),\n+                                           n_q_points));\n+              for (unsigned int q=0; q<n_q_points; q++)\n+                fe_data.enrichment[base_no][m].gradients[q] = enrichments[base_no-1][m](cell)->gradient(mapping_data.quadrature_points[q]);\n+            }\n+\n+          if (flags & update_values)\n+            {\n+              Assert (fe_data.enrichment[base_no][m].values.size() == n_q_points,\n+                      ExcDimensionMismatch(fe_data.enrichment[base_no][m].values.size(),\n+                                           n_q_points));\n+              for (unsigned int q=0; q<n_q_points; q++)\n+                fe_data.enrichment[base_no][m].values[q]    = enrichments[base_no-1][m](cell)->value   (mapping_data.quadrature_points[q]);\n+            }\n+        }\n+    }\n+\n+  // Finally, update the standard data stored in output_data\n+  // by expanding the product rule for enrichment function.\n+  // note that the order if important, namely\n+  // output_data.shape_XYZ contains values of standard FEM and we overwrite\n+  // it with the updated one in the following order: hessians -> gradients -> values\n+  if (flags & update_hessians)\n+    {\n+      Tensor<2, spacedim> grad_grad, grad_grad_r;\n+      Tensor<2,spacedim,double> sym_grad_grad, sym_grad_grad_r;\n+      Tensor<1,spacedim> grad_r;\n+      for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++)\n+        {\n+          for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++)\n+            for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++)\n+              {\n+                const unsigned int &enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i];\n+                for (unsigned int q=0; q<n_q_points; ++q)\n+                  {\n+                    grad_grad =  outer_product(output_data.shape_gradients[enriched_dof][q],  fe_data.enrichment[base_no][m].gradients[q]);\n+                    sym_grad_grad = symmetrize (grad_grad) * 2.0; // symmetrize does [s+s^T]/2","path":"source/fe/fe_enriched.cc","position":null,"original_position":710,"commit_id":"62a2da37079dad53799a791b34f0f2cc03c98b8d","original_commit_id":"0ad35105c8fac12217dd3717b454062684a6180d","user":{"login":"bangerth","id":7504421,"avatar_url":"https://avatars.githubusercontent.com/u/7504421?v=3","gravatar_id":"","url":"https://api.github.com/users/bangerth","html_url":"https://github.com/bangerth","followers_url":"https://api.github.com/users/bangerth/followers","following_url":"https://api.github.com/users/bangerth/following{/other_user}","gists_url":"https://api.github.com/users/bangerth/gists{/gist_id}","starred_url":"https://api.github.com/users/bangerth/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bangerth/subscriptions","organizations_url":"https://api.github.com/users/bangerth/orgs","repos_url":"https://api.github.com/users/bangerth/repos","events_url":"https://api.github.com/users/bangerth/events{/privacy}","received_events_url":"https://api.github.com/users/bangerth/received_events","type":"User","site_admin":false},"body":"I see. Might be worth expanding on in the documentation here or in the class documentation. In fact, writing out all of those (first and second) derivatives may be useful in the class doc because it would be useful to state there anyway that the `Function` object you pass in needs to be able to provide first and second derivatives if you want first or second derivatives of the enriched finite element.","created_at":"2016-09-27T21:57:02Z","updated_at":"2016-09-27T21:57:02Z","html_url":"https://github.com/dealii/dealii/pull/3115#discussion_r80804936","pull_request_url":"https://api.github.com/repos/dealii/dealii/pulls/3115","_links":{"self":{"href":"https://api.github.com/repos/dealii/dealii/pulls/comments/80804936"},"html":{"href":"https://github.com/dealii/dealii/pull/3115#discussion_r80804936"},"pull_request":{"href":"https://api.github.com/repos/dealii/dealii/pulls/3115"}}},"pull_request":{"url":"https://api.github.com/repos/dealii/dealii/pulls/3115","id":85237254,"html_url":"https://github.com/dealii/dealii/pull/3115","diff_url":"https://github.com/dealii/dealii/pull/3115.diff","patch_url":"https://github.com/dealii/dealii/pull/3115.patch","issue_url":"https://api.github.com/repos/dealii/dealii/issues/3115","number":3115,"state":"open","locked":false,"title":"FE_Enriched","user":{"login":"davydden","id":8023934,"avatar_url":"https://avatars.githubusercontent.com/u/8023934?v=3","gravatar_id":"","url":"https://api.github.com/users/davydden","html_url":"https://github.com/davydden","followers_url":"https://api.github.com/users/davydden/followers","following_url":"https://api.github.com/users/davydden/following{/other_user}","gists_url":"https://api.github.com/users/davydden/gists{/gist_id}","starred_url":"https://api.github.com/users/davydden/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/davydden/subscriptions","organizations_url":"https://api.github.com/users/davydden/orgs","repos_url":"https://api.github.com/users/davydden/repos","events_url":"https://api.github.com/users/davydden/events{/privacy}","received_events_url":"https://api.github.com/users/davydden/received_events","type":"User","site_admin":false},"body":"a new FE to do partition of unity method, where enriched and non-enriched elements can be combined.\r\n\r\n- [ ] squash\r\n- [ ] add an entry in `changes.h`\r\n- [ ] change `we store the casted result in` -> `we store the cast result in`\r\n- [x] remove `Threads::Mutex`\r\n- [x] remove destructor.\r\n- [x] instantiate for `dim != spacedim`\r\n- [ ] use references to FEs? (or at least for the base FE) https://github.com/dealii/dealii/pull/3115/files#r80339138\r\n- [x] require c++11\r\n- [x] add a note on disjoint domains and different functions for different cells\r\n- [x] redraw `doc/doxygen/images/fe_enriched_1d.png`","created_at":"2016-09-14T14:23:06Z","updated_at":"2016-09-27T21:57:02Z","closed_at":null,"merged_at":null,"merge_commit_sha":"c817547040edd5835e7c4ae0a32de902cd3426d5","assignee":null,"assignees":[],"milestone":null,"commits_url":"https://api.github.com/repos/dealii/dealii/pulls/3115/commits","review_comments_url":"https://api.github.com/repos/dealii/dealii/pulls/3115/comments","review_comment_url":"https://api.github.com/repos/dealii/dealii/pulls/comments{/number}","comments_url":"https://api.github.com/repos/dealii/dealii/issues/3115/comments","statuses_url":"https://api.github.com/repos/dealii/dealii/statuses/62a2da37079dad53799a791b34f0f2cc03c98b8d","head":{"label":"davydden:feature/fe_enriched","ref":"feature/fe_enriched","sha":"62a2da37079dad53799a791b34f0f2cc03c98b8d","user":{"login":"davydden","id":8023934,"avatar_url":"https://avatars.githubusercontent.com/u/8023934?v=3","gravatar_id":"","url":"https://api.github.com/users/davydden","html_url":"https://github.com/davydden","followers_url":"https://api.github.com/users/davydden/followers","following_url":"https://api.github.com/users/davydden/following{/other_user}","gists_url":"https://api.github.com/users/davydden/gists{/gist_id}","starred_url":"https://api.github.com/users/davydden/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/davydden/subscriptions","organizations_url":"https://api.github.com/users/davydden/orgs","repos_url":"https://api.github.com/users/davydden/repos","events_url":"https://api.github.com/users/davydden/events{/privacy}","received_events_url":"https://api.github.com/users/davydden/received_events","type":"User","site_admin":false},"repo":{"id":22635610,"name":"dealii","full_name":"davydden/dealii","owner":{"login":"davydden","id":8023934,"avatar_url":"https://avatars.githubusercontent.com/u/8023934?v=3","gravatar_id":"","url":"https://api.github.com/users/davydden","html_url":"https://github.com/davydden","followers_url":"https://api.github.com/users/davydden/followers","following_url":"https://api.github.com/users/davydden/following{/other_user}","gists_url":"https://api.github.com/users/davydden/gists{/gist_id}","starred_url":"https://api.github.com/users/davydden/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/davydden/subscriptions","organizations_url":"https://api.github.com/users/davydden/orgs","repos_url":"https://api.github.com/users/davydden/repos","events_url":"https://api.github.com/users/davydden/events{/privacy}","received_events_url":"https://api.github.com/users/davydden/received_events","type":"User","site_admin":false},"private":false,"html_url":"https://github.com/davydden/dealii","description":"The development repository for the deal.II finite element library.","fork":true,"url":"https://api.github.com/repos/davydden/dealii","forks_url":"https://api.github.com/repos/davydden/dealii/forks","keys_url":"https://api.github.com/repos/davydden/dealii/keys{/key_id}","collaborators_url":"https://api.github.com/repos/davydden/dealii/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/davydden/dealii/teams","hooks_url":"https://api.github.com/repos/davydden/dealii/hooks","issue_events_url":"https://api.github.com/repos/davydden/dealii/issues/events{/number}","events_url":"https://api.github.com/repos/davydden/dealii/events","assignees_url":"https://api.github.com/repos/davydden/dealii/assignees{/user}","branches_url":"https://api.github.com/repos/davydden/dealii/branches{/branch}","tags_url":"https://api.github.com/repos/davydden/dealii/tags","blobs_url":"https://api.github.com/repos/davydden/dealii/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/davydden/dealii/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/davydden/dealii/git/refs{/sha}","trees_url":"https://api.github.com/repos/davydden/dealii/git/trees{/sha}","statuses_url":"https://api.github.com/repos/davydden/dealii/statuses/{sha}","languages_url":"https://api.github.com/repos/davydden/dealii/languages","stargazers_url":"https://api.github.com/repos/davydden/dealii/stargazers","contributors_url":"https://api.github.com/repos/davydden/dealii/contributors","subscribers_url":"https://api.github.com/repos/davydden/dealii/subscribers","subscription_url":"https://api.github.com/repos/davydden/dealii/subscription","commits_url":"https://api.github.com/repos/davydden/dealii/commits{/sha}","git_commits_url":"https://api.github.com/repos/davydden/dealii/git/commits{/sha}","comments_url":"https://api.github.com/repos/davydden/dealii/comments{/number}","issue_comment_url":"https://api.github.com/repos/davydden/dealii/issues/comments{/number}","contents_url":"https://api.github.com/repos/davydden/dealii/contents/{+path}","compare_url":"https://api.github.com/repos/davydden/dealii/compare/{base}...{head}","merges_url":"https://api.github.com/repos/davydden/dealii/merges","archive_url":"https://api.github.com/repos/davydden/dealii/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/davydden/dealii/downloads","issues_url":"https://api.github.com/repos/davydden/dealii/issues{/number}","pulls_url":"https://api.github.com/repos/davydden/dealii/pulls{/number}","milestones_url":"https://api.github.com/repos/davydden/dealii/milestones{/number}","notifications_url":"https://api.github.com/repos/davydden/dealii/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/davydden/dealii/labels{/name}","releases_url":"https://api.github.com/repos/davydden/dealii/releases{/id}","deployments_url":"https://api.github.com/repos/davydden/dealii/deployments","created_at":"2014-08-05T08:06:18Z","updated_at":"2015-01-10T19:48:48Z","pushed_at":"2016-09-27T21:08:29Z","git_url":"git://github.com/davydden/dealii.git","ssh_url":"git@github.com:davydden/dealii.git","clone_url":"https://github.com/davydden/dealii.git","svn_url":"https://github.com/davydden/dealii","homepage":"http://www.dealii.org/","size":147791,"stargazers_count":0,"watchers_count":0,"language":"C++","has_issues":false,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":0,"mirror_url":null,"open_issues_count":0,"forks":0,"open_issues":0,"watchers":0,"default_branch":"master"}},"base":{"label":"dealii:master","ref":"master","sha":"2863c8df4cc3c3f94739bf929292e704a55fdf15","user":{"login":"dealii","id":7506070,"avatar_url":"https://avatars.githubusercontent.com/u/7506070?v=3","gravatar_id":"","url":"https://api.github.com/users/dealii","html_url":"https://github.com/dealii","followers_url":"https://api.github.com/users/dealii/followers","following_url":"https://api.github.com/users/dealii/following{/other_user}","gists_url":"https://api.github.com/users/dealii/gists{/gist_id}","starred_url":"https://api.github.com/users/dealii/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/dealii/subscriptions","organizations_url":"https://api.github.com/users/dealii/orgs","repos_url":"https://api.github.com/users/dealii/repos","events_url":"https://api.github.com/users/dealii/events{/privacy}","received_events_url":"https://api.github.com/users/dealii/received_events","type":"Organization","site_admin":false},"repo":{"id":21943980,"name":"dealii","full_name":"dealii/dealii","owner":{"login":"dealii","id":7506070,"avatar_url":"https://avatars.githubusercontent.com/u/7506070?v=3","gravatar_id":"","url":"https://api.github.com/users/dealii","html_url":"https://github.com/dealii","followers_url":"https://api.github.com/users/dealii/followers","following_url":"https://api.github.com/users/dealii/following{/other_user}","gists_url":"https://api.github.com/users/dealii/gists{/gist_id}","starred_url":"https://api.github.com/users/dealii/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/dealii/subscriptions","organizations_url":"https://api.github.com/users/dealii/orgs","repos_url":"https://api.github.com/users/dealii/repos","events_url":"https://api.github.com/users/dealii/events{/privacy}","received_events_url":"https://api.github.com/users/dealii/received_events","type":"Organization","site_admin":false},"private":false,"html_url":"https://github.com/dealii/dealii","description":"The development repository for the deal.II finite element library.","fork":false,"url":"https://api.github.com/repos/dealii/dealii","forks_url":"https://api.github.com/repos/dealii/dealii/forks","keys_url":"https://api.github.com/repos/dealii/dealii/keys{/key_id}","collaborators_url":"https://api.github.com/repos/dealii/dealii/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/dealii/dealii/teams","hooks_url":"https://api.github.com/repos/dealii/dealii/hooks","issue_events_url":"https://api.github.com/repos/dealii/dealii/issues/events{/number}","events_url":"https://api.github.com/repos/dealii/dealii/events","assignees_url":"https://api.github.com/repos/dealii/dealii/assignees{/user}","branches_url":"https://api.github.com/repos/dealii/dealii/branches{/branch}","tags_url":"https://api.github.com/repos/dealii/dealii/tags","blobs_url":"https://api.github.com/repos/dealii/dealii/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/dealii/dealii/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/dealii/dealii/git/refs{/sha}","trees_url":"https://api.github.com/repos/dealii/dealii/git/trees{/sha}","statuses_url":"https://api.github.com/repos/dealii/dealii/statuses/{sha}","languages_url":"https://api.github.com/repos/dealii/dealii/languages","stargazers_url":"https://api.github.com/repos/dealii/dealii/stargazers","contributors_url":"https://api.github.com/repos/dealii/dealii/contributors","subscribers_url":"https://api.github.com/repos/dealii/dealii/subscribers","subscription_url":"https://api.github.com/repos/dealii/dealii/subscription","commits_url":"https://api.github.com/repos/dealii/dealii/commits{/sha}","git_commits_url":"https://api.github.com/repos/dealii/dealii/git/commits{/sha}","comments_url":"https://api.github.com/repos/dealii/dealii/comments{/number}","issue_comment_url":"https://api.github.com/repos/dealii/dealii/issues/comments{/number}","contents_url":"https://api.github.com/repos/dealii/dealii/contents/{+path}","compare_url":"https://api.github.com/repos/dealii/dealii/compare/{base}...{head}","merges_url":"https://api.github.com/repos/dealii/dealii/merges","archive_url":"https://api.github.com/repos/dealii/dealii/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/dealii/dealii/downloads","issues_url":"https://api.github.com/repos/dealii/dealii/issues{/number}","pulls_url":"https://api.github.com/repos/dealii/dealii/pulls{/number}","milestones_url":"https://api.github.com/repos/dealii/dealii/milestones{/number}","notifications_url":"https://api.github.com/repos/dealii/dealii/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/dealii/dealii/labels{/name}","releases_url":"https://api.github.com/repos/dealii/dealii/releases{/id}","deployments_url":"https://api.github.com/repos/dealii/dealii/deployments","created_at":"2014-07-17T14:00:26Z","updated_at":"2016-09-24T21:10:45Z","pushed_at":"2016-09-27T21:08:30Z","git_url":"git://github.com/dealii/dealii.git","ssh_url":"git@github.com:dealii/dealii.git","clone_url":"https://github.com/dealii/dealii.git","svn_url":"https://github.com/dealii/dealii","homepage":"http://www.dealii.org/","size":168264,"stargazers_count":172,"watchers_count":172,"language":"C++","has_issues":true,"has_downloads":true,"has_wiki":true,"has_pages":false,"forks_count":185,"mirror_url":null,"open_issues_count":185,"forks":185,"open_issues":185,"watchers":172,"default_branch":"master"}},"_links":{"self":{"href":"https://api.github.com/repos/dealii/dealii/pulls/3115"},"html":{"href":"https://github.com/dealii/dealii/pull/3115"},"issue":{"href":"https://api.github.com/repos/dealii/dealii/issues/3115"},"comments":{"href":"https://api.github.com/repos/dealii/dealii/issues/3115/comments"},"review_comments":{"href":"https://api.github.com/repos/dealii/dealii/pulls/3115/comments"},"review_comment":{"href":"https://api.github.com/repos/dealii/dealii/pulls/comments{/number}"},"commits":{"href":"https://api.github.com/repos/dealii/dealii/pulls/3115/commits"},"statuses":{"href":"https://api.github.com/repos/dealii/dealii/statuses/62a2da37079dad53799a791b34f0f2cc03c98b8d"}}}}
{"actor":{"display_login":"bangerth"}}