allow different base_elements for compositional fields

benchmarks/solitary_wave/solitary_wave.cc

@@ -22,6 +22,7 @@
 #include <aspect/postprocess/interface.h>
 #include <aspect/gravity_model/interface.h>
 #include <aspect/geometry_model/interface.h>
+#include <aspect/simulator.h>
 #include <aspect/simulator_access.h>
 #include <aspect/global.h>
 
@@ -828,13 +829,13 @@ namespace aspect
       AssertThrow(this->introspection().compositional_name_exists("porosity"),
                   ExcMessage("Postprocessor Solitary Wave only works if there is a compositional field called porosity."));
       const unsigned int porosity_index = this->introspection().compositional_index_for_name("porosity");
+      const typename Simulator<dim>::AdvectionField porosity = Simulator<dim>::AdvectionField::composition(porosity_index);
 
       // 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->introspection().n_compositional_fields > 0,
                    ExcMessage("This postprocessor cannot be used without compositional fields."));
-      const QGauss<dim> quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.compositional_fields).degree+1);
+
+      const QGauss<dim> quadrature_formula (this->get_fe().base_element(porosity.base_element(this->introspection())).degree+1);
       const unsigned int n_q_points = quadrature_formula.size();
 
       FEValues<dim> fe_values (this->get_mapping(),

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,25 @@ 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. If you have no compositional fields, the vector has length 0.
+       * If all compositional fields have the same finite element space, the length is 1.
+       */
+      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 as returned by get_compositional_field_base_element_indices().
+       * The indices returned are therefore between 0 and n_compositional_fields-1.
+       * If you have a single compositional field, this function returns {0} when passing
+       * in base_element_index=0.
+       */
+      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 +612,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/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>
@@ -37,44 +38,50 @@ namespace aspect
                                "compositional fields are active!"));
 
       indicators = 0;
-      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 base_element_index : this->introspection().get_compositional_field_base_element_indices())
+        {
+          const Quadrature<dim> quadrature (this->get_fe().base_element(base_element_index).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(base_element_index))
+                  {
+                    const typename Simulator<dim>::AdvectionField composition = Simulator<dim>::AdvectionField::composition(c);
+                    fe_values[composition.scalar_extractor(this->introspection())].get_function_gradients (this->get_solution(),
+                        composition_gradients);
+
+                    // Some up the indicators for this composition on this cell. Note that quadrature points and dofs
+                    // are enumerated in the same order.
+                    double this_indicator = 0.0;
+                    for (unsigned int j=0; j<dofs_per_cell; ++j)
+                      this_indicator += 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 *= std::pow(cell->diameter(), power);
+
+                    indicators[idx] += 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;