allow different base_elements for compositional fields

include/aspect/introspection.h

@@ -246,18 +246,17 @@ namespace aspect
        * A structure that enumerates the base elements of the finite element
        * that correspond to each of the variables in this problem.
        *
-       * If there are compositional fields, they are all discretized with the
-       * same base element and, consequently, we only need a single index. If
-       * a variable does not exist in the problem (e.g., we do not have
-       * compositional fields), then the corresponding index is set to an
-       * invalid number.
+       * The indices here can be used to access the dealii::FiniteElement
+       * that describes the given variable. We support different finite
+       * elements for compositional fields, but we try to reuse the same
+       * element if possible.
        */
       struct BaseElements
       {
-        unsigned int       velocities;
-        unsigned int       pressure;
-        unsigned int       temperature;
-        unsigned int       compositional_fields;
+        unsigned int              velocities;
+        unsigned int              pressure;
+        unsigned int              temperature;
+        std::vector<unsigned int> compositional_fields;
       };
 
       /**
@@ -441,6 +440,7 @@ namespace aspect
          */
         IndexSet locally_owned_fluid_pressure_dofs;
       };
+
       /**
        * A variable that contains index sets describing which of the globally
        * enumerated degrees of freedom are owned by the current processor in a
@@ -465,6 +465,21 @@ namespace aspect
        * @}
        */
 
+      /**
+       * Return a vector with all used indices of base elements of the deal.II FiniteElement
+       * that are used for compositional fields. If several fields are the same type, they
+       * share base elements.
+      */
+      std::vector<unsigned int>
+      get_compositional_field_base_element_indices() const;
+
+      /**
+       * Return a vector with all compositional field indices that belong to a given
+       * base element index.
+      */
+      std::vector<unsigned int>
+      get_compositional_field_indices_with_base_element(const unsigned int base_element_index) const;
+
       /**
        * A function that gets the name of a compositional field as an input
        * parameter and returns its index. If the name is not found, an
@@ -593,6 +608,16 @@ namespace aspect
       bool
       is_stokes_component (const unsigned int component_index) const;
 
+      /**
+       * A function that gets a component index as an input
+       * parameter and returns if the component is one of the
+       * compositional fields.
+       *
+       * @param component_index The component index to check.
+       */
+      bool
+      is_composition_component (const unsigned int component_index) const;
+
     private:
       /**
        * A vector that stores the names of the compositional fields that will

source/mesh_refinement/compaction_length.cc

@@ -40,7 +40,12 @@ namespace aspect
       if (this->get_dof_handler().n_locally_owned_dofs() == 0)
         return;
 
-      const Quadrature<dim> quadrature(this->get_fe().base_element(this->introspection().base_elements.compositional_fields).get_unit_support_points());
+      // Use a quadrature in the support points of the porosity to compute the
+      // compaction length at:
+      const typename Simulator<dim>::AdvectionField porosity = Simulator<dim>::AdvectionField::composition(
+                                                                 this->introspection().find_composition_type(CompositionalFieldDescription::porosity));
+      const Quadrature<dim> quadrature(this->get_fe().base_element(porosity.base_element(this->introspection())).get_unit_support_points());
+
       FEValues<dim> fe_values (this->get_mapping(),
                                this->get_fe(),
                                quadrature,
@@ -65,7 +70,7 @@ namespace aspect
                         ExcMessage("Need MeltOutputs from the material model for computing the melt properties."));
 
             // for each composition dof, check if the compaction length exceeds the cell size
-            for (unsigned int i=0; i<this->get_fe().base_element(this->introspection().base_elements.compositional_fields).dofs_per_cell; ++i)
+            for (unsigned int i=0; i<in.n_evaluation_points(); ++i)
               {
                 const double compaction_length = std::sqrt((out.viscosities[i] + 4./3. * melt_out->compaction_viscosities[i])
                                                            * melt_out->permeabilities[i] / melt_out->fluid_viscosities[i]);

source/mesh_refinement/composition_gradient.cc

@@ -19,6 +19,7 @@
 */
 
 
+#include <aspect/simulator.h>
 #include <aspect/mesh_refinement/composition_gradient.h>
 
 #include <deal.II/base/quadrature_lib.h>
@@ -40,41 +41,46 @@ namespace aspect
       Vector<float> this_indicator (indicators.size());
       const double power = 1.0 + dim/2.0;
 
-      const Quadrature<dim> quadrature(this->get_fe().base_element(this->introspection().base_elements.compositional_fields).get_unit_support_points());
-      FEValues<dim> fe_values (this->get_mapping(),
-                               this->get_fe(),
-                               quadrature,
-                               update_quadrature_points | update_gradients);
+      for (const unsigned int k : this->introspection().get_compositional_field_base_element_indices())
+        {
+          const Quadrature<dim> quadrature (composition.base_element(k).get_unit_support_points());
+          const unsigned int dofs_per_cell = quadrature.size();
+          FEValues<dim> fe_values (this->get_mapping(),
+                                   this->get_fe(),
+                                   quadrature,
+                                   update_quadrature_points | update_gradients);
 
-      // the values of the compositional fields are stored as block vectors for each field
-      // we have to extract them in this structure
-      std::vector<Tensor<1,dim>> composition_gradients (quadrature.size());
+          // the values of the compositional fields are stored as block vectors for each field
+          // we have to extract them in this structure
+          std::vector<Tensor<1,dim>> composition_gradients (quadrature.size());
 
-      for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
-        {
           for (const auto &cell : this->get_dof_handler().active_cell_iterators())
             if (cell->is_locally_owned())
               {
                 const unsigned int idx = cell->active_cell_index();
                 fe_values.reinit(cell);
-                fe_values[this->introspection().extractors.compositional_fields[c]].get_function_gradients (this->get_solution(),
-                    composition_gradients);
-
-                // For each composition dof, write into the output vector the
-                // composition gradient. Note that quadrature points and dofs
-                // are enumerated in the same order.
-                for (unsigned int j=0; j<this->get_fe().base_element(this->introspection().base_elements.compositional_fields).dofs_per_cell; ++j)
-                  this_indicator[idx] += composition_gradients[j].norm();
-
-                // Scale gradient in each cell with the correct power of h. Otherwise,
-                // error indicators do not reduce when refined if there is a density
-                // jump. We need at least order 1 for the error not to grow when
-                // refining, so anything >1 should work. (note that the gradient
-                // itself scales like 1/h, so multiplying it with any factor h^s, s>1
-                // will yield convergence of the error indicators to zero as h->0)
-                this_indicator[idx] *= std::pow(cell->diameter(), power);
+
+                for (const unsigned int c : this->introspection().get_compositional_field_indices_with_base_element(k))
+                  {
+                    fe_values[composition.scalar_extractor(c)].get_function_gradients (this->get_solution(),
+                                                                                       composition_gradients);
+
+                    // For each composition dof, write into the output vector the
+                    // composition gradient. Note that quadrature points and dofs
+                    // are enumerated in the same order.
+                    for (unsigned int j=0; j<dofs_per_cell; ++j)
+                      this_indicator[idx] += composition_gradients[j].norm();
+
+                    // Scale gradient in each cell with the correct power of h. Otherwise,
+                    // error indicators do not reduce when refined if there is a density
+                    // jump. We need at least order 1 for the error not to grow when
+                    // refining, so anything >1 should work. (note that the gradient
+                    // itself scales like 1/h, so multiplying it with any factor h^s, s>1
+                    // will yield convergence of the error indicators to zero as h->0)
+                    this_indicator[idx] *= std::pow(cell->diameter(), power);
+                  }
+                indicators.add(composition_scaling_factors[c], this_indicator);
               }
-          indicators.add(composition_scaling_factors[c], this_indicator);
         }
     }
 

source/postprocess/composition_statistics.cc

@@ -37,10 +37,6 @@ namespace aspect
         return {"", ""};
 
       // create a quadrature formula based on the compositional element alone.
-      // be defensive about determining that a compositional field actually exists
-      AssertThrow (this->introspection().base_elements.compositional_fields
-                   != numbers::invalid_unsigned_int,
-                   ExcMessage("This postprocessor cannot be used without compositional fields."));
       const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_fields;
       const unsigned int n_q_points = quadrature_formula.size();
 

source/postprocess/max_depth_field.cc

@@ -36,7 +36,7 @@ namespace aspect
     {
       // create a quadrature formula based on the compositional element alone.
       // be defensive about determining that a compositional field actually exists
-      AssertThrow(this->introspection().base_elements.compositional_fields != numbers::invalid_unsigned_int,
+      AssertThrow(this->n_compositional_fields() > 0,
                   ExcMessage("This postprocessor cannot be used without compositional fields."));
       const Quadrature<dim> &quadrature_formula = this->introspection().quadratures.compositional_fields;
 

source/simulator/helper_functions.cc

@@ -174,7 +174,7 @@ namespace aspect
     if (this->is_temperature())
       return introspection.base_elements.temperature;
     else
-      return introspection.base_elements.compositional_fields;
+      return introspection.base_elements.compositional_fields[compositional_variable];
   }
 
   template <int dim>

source/simulator/introspection.cc

@@ -56,14 +56,13 @@ namespace aspect
      */
     template <int dim>
     typename Introspection<dim>::BlockIndices
-    setup_blocks (FEVariableCollection<dim> &fevs)
+    setup_blocks (FEVariableCollection<dim> &fevs, const unsigned int n_compositional_fields)
     {
       typename Introspection<dim>::BlockIndices b;
       b.velocities = fevs.variable("velocity").block_index;
       b.pressure = fevs.variable("pressure").block_index;
       b.temperature = fevs.variable("temperature").block_index;
 
-      unsigned int n_compositional_fields = fevs.variable("compositions").n_components();
       const unsigned int first_composition_block_index = fevs.variable("compositions").block_index;
       for (unsigned int i=0; i<n_compositional_fields; ++i)
         {
@@ -78,14 +77,15 @@ namespace aspect
 
     template <int dim>
     typename Introspection<dim>::BaseElements
-    setup_base_elements (FEVariableCollection<dim> &fevs)
+    setup_base_elements (FEVariableCollection<dim> &fevs, const unsigned int n_compositional_fields)
     {
       typename Introspection<dim>::BaseElements base_elements;
 
       base_elements.velocities = fevs.variable("velocity").base_index;
       base_elements.pressure = fevs.variable("pressure").base_index;
       base_elements.temperature = fevs.variable("temperature").base_index;
-      base_elements.compositional_fields = fevs.variable("compositions").base_index;
+      base_elements.compositional_fields = Utilities::possibly_extend_from_1_to_N(std::vector<unsigned int>({fevs.variable("compositions").base_index}),
+                                                                                  n_compositional_fields, "");
 
       return base_elements;
     }
@@ -247,9 +247,9 @@ namespace aspect
     use_discontinuous_composition_discretization (parameters.use_discontinuous_composition_discretization),
     component_indices (internal::setup_component_indices<dim>(*this)),
     n_blocks(FEVariableCollection<dim>::n_blocks()),
-    block_indices (internal::setup_blocks<dim>(*this)),
+    block_indices (internal::setup_blocks<dim>(*this, parameters.n_compositional_fields)),
     extractors (component_indices),
-    base_elements (internal::setup_base_elements<dim>(*this)),
+    base_elements (internal::setup_base_elements<dim>(*this, parameters.n_compositional_fields)),
     polynomial_degree (internal::setup_polynomial_degree<dim>(parameters)),
     quadratures (internal::setup_quadratures<dim>(parameters, ReferenceCells::get_hypercube<dim>())),
     face_quadratures (internal::setup_face_quadratures<dim>(parameters, ReferenceCells::get_hypercube<dim>())),
@@ -332,6 +332,37 @@ namespace aspect
 
 
 
+  template <int dim>
+  std::vector<unsigned int>
+  Introspection<dim>::get_compositional_field_base_element_indices() const
+  {
+    std::vector<unsigned int> result;
+    const unsigned int last = *this->base_elements.compositional_fields.end();
+    for (unsigned int i = this->base_elements.compositional_fields[0]; i<=last; ++i)
+      {
+        result.emplace_back(i);
+      }
+    return result;
+  }
+
+
+
+  template <int dim>
+  std::vector<unsigned int>
+  Introspection<dim>::get_compositional_field_indices_with_base_element(const unsigned int base_element_index) const
+  {
+    std::vector<unsigned int> result;
+
+    for (const auto idx : this->base_elements.compositional_fields)
+      {
+        if (idx == base_element_index)
+          result.emplace_back(idx);
+      }
+    return result;
+  }
+
+
+
   template <int dim>
   unsigned int
   Introspection<dim>::compositional_index_for_name (const std::string &name) const
@@ -487,6 +518,22 @@ namespace aspect
   }
 
 
+
+  template <int dim>
+  bool
+  Introspection<dim>::is_composition_component (const unsigned int component_index) const
+  {
+    // All compositions live at the end. Just to be sure, there are no other components
+    // in our system after compositional fields, right?
+    Assert(component_indices.compositional_fields[0] > component_indices.temperature
+           && component_indices.compositional_fields.last() == n_components-1, ExcInternalError());
+
+    if (component_index >= component_indices.compositional_fields[0])
+      return true;
+    else
+      return false;
+  }
+
 }
 
 

unit_tests/introspection.cc

@@ -67,4 +67,11 @@ TEST_CASE("Introspection::1")
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[0] == 3);
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[1] == 3);
   CHECK(introspection.block_indices.compositional_field_sparsity_pattern[2] == 3);
+
+  CHECK(introspection.base_elements.velocities == 0);
+  CHECK(introspection.base_elements.pressure == 1);
+  CHECK(introspection.base_elements.temperature == 2);
+  // all compositional fields are the same FE:
+  CHECK(introspection.base_elements.compositional_fields == std::vector<unsigned int>({3,3,3}));
+  CHECK(introspection.get_compositional_field_base_element_indices() == std::vector<unsigned int>({4}));
 }