From 7831a7977fba4e1cd29a5cfb6f27cc29d2dbbd70 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Wed, 5 Jun 2024 17:49:48 -0400 Subject: [PATCH] allow different base_elements for compositional fields --- benchmarks/solitary_wave/solitary_wave.cc | 9 +- include/aspect/introspection.h | 47 +++- .../mesh_refinement/composition_gradient.cc | 63 ++--- source/postprocess/composition_statistics.cc | 4 - source/postprocess/max_depth_field.cc | 2 +- source/simulator/helper_functions.cc | 227 +++++++++--------- source/simulator/introspection.cc | 72 +++++- unit_tests/introspection.cc | 8 + 8 files changed, 266 insertions(+), 166 deletions(-) diff --git a/benchmarks/solitary_wave/solitary_wave.cc b/benchmarks/solitary_wave/solitary_wave.cc index 202350a76dc..cbfba28ca6c 100644 --- a/benchmarks/solitary_wave/solitary_wave.cc +++ b/benchmarks/solitary_wave/solitary_wave.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -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::AdvectionField porosity = Simulator::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 quadrature_formula (this->get_fe().base_element(this->introspection().base_elements.compositional_fields).degree+1); + + const QGauss quadrature_formula (this->get_fe().base_element(porosity.base_element(this->introspection())).degree+1); const unsigned int n_q_points = quadrature_formula.size(); FEValues fe_values (this->get_mapping(), diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h index 163e1d01538..dd573427745 100644 --- a/include/aspect/introspection.h +++ b/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 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 + 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 + 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 diff --git a/source/mesh_refinement/composition_gradient.cc b/source/mesh_refinement/composition_gradient.cc index ec313514af0..6d7ad428952 100644 --- a/source/mesh_refinement/composition_gradient.cc +++ b/source/mesh_refinement/composition_gradient.cc @@ -19,6 +19,7 @@ */ +#include #include #include @@ -37,44 +38,50 @@ namespace aspect "compositional fields are active!")); indicators = 0; - Vector this_indicator (indicators.size()); const double power = 1.0 + dim/2.0; - const Quadrature quadrature(this->get_fe().base_element(this->introspection().base_elements.compositional_fields).get_unit_support_points()); - FEValues 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 quadrature (this->get_fe().base_element(base_element_index).get_unit_support_points()); + const unsigned int dofs_per_cell = quadrature.size(); + FEValues 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> 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> composition_gradients (quadrature.size()); - for (unsigned int c=0; cn_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; jget_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::AdvectionField composition = Simulator::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; j1 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); } } diff --git a/source/postprocess/composition_statistics.cc b/source/postprocess/composition_statistics.cc index 748d22e1c17..d0073fcf11c 100644 --- a/source/postprocess/composition_statistics.cc +++ b/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 &quadrature_formula = this->introspection().quadratures.compositional_fields; const unsigned int n_q_points = quadrature_formula.size(); diff --git a/source/postprocess/max_depth_field.cc b/source/postprocess/max_depth_field.cc index fdd9f43d6bf..d0cba1d2a56 100644 --- a/source/postprocess/max_depth_field.cc +++ b/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 &quadrature_formula = this->introspection().quadratures.compositional_fields; diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index f100bede553..a5192eda657 100644 --- a/source/simulator/helper_functions.cc +++ b/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 @@ -184,7 +184,11 @@ namespace aspect if (this->is_temperature()) return introspection.extractors.temperature; else - return introspection.extractors.compositional_fields[compositional_variable]; + { + Assert(compositional_variable < introspection.n_compositional_fields, + ExcMessage("Invalid AdvectionField.")); + return introspection.extractors.compositional_fields[compositional_variable]; + } } template @@ -1686,55 +1690,80 @@ namespace aspect pcout << " Solving composition reactions... " << std::flush; - // make one fevalues for the composition, and one for the temperature (they might use different finite elements) - const Quadrature quadrature_C(dof_handler.get_fe().base_element(introspection.base_elements.compositional_fields).get_unit_support_points()); - FEValues fe_values_C (*mapping, - dof_handler.get_fe(), - quadrature_C, - update_quadrature_points | update_values | update_gradients); - - std::vector local_dof_indices (dof_handler.get_fe().dofs_per_cell); - MaterialModel::MaterialModelInputs in_C(quadrature_C.size(), introspection.n_compositional_fields); - MaterialModel::MaterialModelOutputs out_C(quadrature_C.size(), introspection.n_compositional_fields); - HeatingModel::HeatingModelOutputs heating_model_outputs_C(quadrature_C.size(), introspection.n_compositional_fields); + // We want to compute reactions in each support point for all fields (compositional fields and temperature). The reaction + // rate for an individual field depends on the values of all other fields, so we have to step them forward in time together. + // The rates comes from the material and heating model, otherwise we have a simple ODE in each point on each cell to solve. + // + // So far so good. Except that fields can have different Finite Element discretizations (degree, continuous/discontinuous) + // and will have different support points. We solve this by computing the union of all support points and evaluating all fields + // in these points for every cell (if necessary by interpolation). Then we solve the ODEs on each cell together and + // write back all values that correspond to support points of that particular field. Field values we computed that are + // not actually support points, we just throw away. + + // First compute the union of all support points (temperature and compositions): + std::vector> unique_support_points; + // For each field (T or composition) store where each support point is found in the list of indices of support points + // unique_support_points. This means index=support_point_index_by_field[0][3] is the index of the third support point + // of the temperature field and unique_support_points[index] is its location. + std::vector> support_point_index_by_field(1+introspection.n_compositional_fields); + + // Fill in the support_point_index_by_field data structure. This is a bit complicated... + { + // first fill in temperature: + unique_support_points = dof_handler.get_fe().base_element(introspection.base_elements.temperature).get_unit_support_points(); + for (unsigned int i=0; i::AdvectionField::composition(c); + const std::vector> &points = dof_handler.get_fe().base_element(composition.base_element(introspection)).get_unit_support_points(); - // temperature element - const Quadrature quadrature_T(dof_handler.get_fe().base_element(introspection.base_elements.temperature).get_unit_support_points()); + for (unsigned int i=0; i fe_values_T (*mapping, - dof_handler.get_fe(), - quadrature_T, - update_quadrature_points | update_values | update_gradients); + const Quadrature combined_support_points(unique_support_points); + FEValues fe_values (*mapping, + dof_handler.get_fe(), + combined_support_points, + update_quadrature_points | update_values | update_gradients); - MaterialModel::MaterialModelInputs in_T(quadrature_T.size(), introspection.n_compositional_fields); - MaterialModel::MaterialModelOutputs out_T(quadrature_T.size(), introspection.n_compositional_fields); - HeatingModel::HeatingModelOutputs heating_model_outputs_T(quadrature_T.size(), introspection.n_compositional_fields); + std::vector local_dof_indices (dof_handler.get_fe().dofs_per_cell); + MaterialModel::MaterialModelInputs in(combined_support_points.size(), introspection.n_compositional_fields); + MaterialModel::MaterialModelOutputs out(combined_support_points.size(), introspection.n_compositional_fields); + HeatingModel::HeatingModelOutputs heating_model_outputs(combined_support_points.size(), introspection.n_compositional_fields); // add reaction rate outputs - material_model->create_additional_named_outputs(out_C); - material_model->create_additional_named_outputs(out_T); - - MaterialModel::ReactionRateOutputs *reaction_rate_outputs_C - = out_C.template get_additional_output>(); + material_model->create_additional_named_outputs(out); - MaterialModel::ReactionRateOutputs *reaction_rate_outputs_T - = out_T.template get_additional_output>(); + MaterialModel::ReactionRateOutputs *reaction_rate_outputs + = out.template get_additional_output>(); - AssertThrow(reaction_rate_outputs_C != nullptr && reaction_rate_outputs_T != nullptr, + AssertThrow(reaction_rate_outputs != nullptr, ExcMessage("You are trying to use the operator splitting solver scheme, " "but the material model you use does not support operator splitting " "(it does not create ReactionRateOutputs, which are required for this " "solver scheme).")); // some heating models require the additional outputs - heating_model_manager.create_additional_material_model_inputs_and_outputs(in_C, out_C); - heating_model_manager.create_additional_material_model_inputs_and_outputs(in_T, out_T); + heating_model_manager.create_additional_material_model_inputs_and_outputs(in, out); // Make a loop first over all cells, than over all reaction time steps, and then over // all degrees of freedom in each element to compute the reactions. This is possible @@ -1746,114 +1775,88 @@ namespace aspect // interface between cells), as we loop over all cells, and then over all degrees of freedom // on each cell. Although this means we do some additional work, the results are still // correct, as we never read from distributed_vector inside the loop over all cells. - // We initialize the material model inputs objects in_T and in_C using the solution vector + // We initialize the material model inputs object in using the solution vector // on every cell, compute the update, and then on every cell put the result into the // distributed_vector vector. Only after the loop over all cells do we copy distributed_vector // back onto the solution vector. // So even though we touch some DoF more than once, we always start from the same value, compute the // same value, and then overwrite the same value in distributed_vector. - // TODO: make this more efficient. + // TODO: make this even more efficient. for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned()) { - fe_values_C.reinit (cell); - in_C.reinit(fe_values_C, cell, introspection, solution); - - if (temperature_and_composition_use_same_fe == false) - { - fe_values_T.reinit (cell); - in_T.reinit(fe_values_T, cell, introspection, solution); - } + fe_values.reinit (cell); + in.reinit(fe_values, cell, introspection, solution); - std::vector> accumulated_reactions_C (quadrature_C.size(),std::vector (introspection.n_compositional_fields)); - std::vector accumulated_reactions_T (quadrature_T.size()); + std::vector> accumulated_reactions_C (combined_support_points.size(), std::vector (introspection.n_compositional_fields)); + std::vector accumulated_reactions_T (combined_support_points.size()); // Make the reaction time steps: We have to update the values of compositional fields and the temperature. - // Because temperature and composition might use different finite elements, we loop through their elements - // separately, and update the temperature and the compositions for both. // We can reuse the same material model inputs and outputs structure for each reaction time step. // We store the computed updates to temperature and composition in a separate (accumulated_reactions) vector, // so that we can later copy it over to the solution vector. for (unsigned int i=0; ifill_additional_material_model_inputs(in_C, solution, fe_values_C, introspection); + material_model->fill_additional_material_model_inputs(in, solution, fe_values, introspection); - material_model->evaluate(in_C, out_C); - heating_model_manager.evaluate(in_C, out_C, heating_model_outputs_C); + material_model->evaluate(in, out); + heating_model_manager.evaluate(in, out, heating_model_outputs); - for (unsigned int j=0; jreaction_rates[j][c]; - accumulated_reactions_C[j][c] += reaction_time_step_size * reaction_rate_outputs_C->reaction_rates[j][c]; + in.composition[j][c] = in.composition[j][c] + + reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c]; + accumulated_reactions_C[j][c] += reaction_time_step_size * reaction_rate_outputs->reaction_rates[j][c]; } - in_C.temperature[j] = in_C.temperature[j] - + reaction_time_step_size * heating_model_outputs_C.rates_of_temperature_change[j]; + in.temperature[j] = in.temperature[j] + + reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j]; - if (temperature_and_composition_use_same_fe) - accumulated_reactions_T[j] += reaction_time_step_size * heating_model_outputs_C.rates_of_temperature_change[j]; - } - - if (!temperature_and_composition_use_same_fe) - { - // loop over temperature element - material_model->fill_additional_material_model_inputs(in_T, solution, fe_values_T, introspection); - - material_model->evaluate(in_T, out_T); - heating_model_manager.evaluate(in_T, out_T, heating_model_outputs_T); - - for (unsigned int j=0; jreaction_rates[j][c]; - } + accumulated_reactions_T[j] += reaction_time_step_size * heating_model_outputs.rates_of_temperature_change[j]; } } cell->get_dof_indices (local_dof_indices); + const unsigned int component_idx_T = introspection.component_indices.temperature; - // copy reaction rates and new values for the compositional fields - for (unsigned int j=0; j=component_idx_T) // ignore velocity, pressure, etc. { - if (temperature_and_composition_use_same_fe) - distributed_vector(local_dof_indices[temperature_idx]) = in_C.temperature[j]; - else - distributed_vector(local_dof_indices[temperature_idx]) = in_T.temperature[j]; - - distributed_reaction_vector(local_dof_indices[temperature_idx]) = accumulated_reactions_T[j]; + // We found a DoF that belongs to component component_idx, which is a temperature or compositional + // field. That means we want to find where this DoF in the computed reactions above to copy it + // back into the global solution vector. + + // These two variables tell us the how-manyth shape function of which field (and therefore + // field) this DoF is: + const unsigned int index_within = comp_pair.second; + const unsigned int field_index = component_idx-component_idx_T; + // Now we can look up in the support_point_index_by_field data structure where this support + // point is in the list of unique_support_points (and in the Quadrature): + const unsigned int idx = support_point_index_by_field[field_index][index_within]; + + // The final step is grabbing the value from the reaction computation and write it into + // the global vector (if we own it, of course): + if (dof_handler.locally_owned_dofs().is_element(local_dof_indices[dof_idx])) + { + // temperatures and compositions are stored differently: + if (component_idx == component_idx_T) + { + distributed_vector(local_dof_indices[dof_idx]) = in.temperature[idx]; + distributed_reaction_vector(local_dof_indices[dof_idx]) = accumulated_reactions_T[idx]; + } + else + { + const unsigned int composition = field_index-1; // 0 is temperature... + distributed_vector(local_dof_indices[dof_idx]) = in.composition[idx][composition]; + distributed_reaction_vector(local_dof_indices[dof_idx]) = accumulated_reactions_C[idx][composition]; + } + } } } } diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc index c1378342368..15ddb5c4ae6 100644 --- a/source/simulator/introspection.cc +++ b/source/simulator/introspection.cc @@ -21,7 +21,7 @@ #include #include - +#include #include #include @@ -56,14 +56,13 @@ namespace aspect */ template typename Introspection::BlockIndices - setup_blocks (FEVariableCollection &fevs) + setup_blocks (FEVariableCollection &fevs, const unsigned int n_compositional_fields) { typename Introspection::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 typename Introspection::BaseElements - setup_base_elements (FEVariableCollection &fevs) + setup_base_elements (FEVariableCollection &fevs, const unsigned int n_compositional_fields) { typename Introspection::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({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(*this)), n_blocks(FEVariableCollection::n_blocks()), - block_indices (internal::setup_blocks(*this)), + block_indices (internal::setup_blocks(*this, parameters.n_compositional_fields)), extractors (component_indices), - base_elements (internal::setup_base_elements(*this)), + base_elements (internal::setup_base_elements(*this, parameters.n_compositional_fields)), polynomial_degree (internal::setup_polynomial_degree(parameters)), quadratures (internal::setup_quadratures(parameters, ReferenceCells::get_hypercube())), face_quadratures (internal::setup_face_quadratures(parameters, ReferenceCells::get_hypercube())), @@ -333,6 +333,46 @@ namespace aspect + template + std::vector + Introspection::get_compositional_field_base_element_indices() const + { + // We are assigning base elements in order, so the first compositional field + // gives us the first base element index. Then we find the largest index + // in the vector. This is necessary, because the fields could have type A,B,A. + std::vector result; + const unsigned int first = this->base_elements.compositional_fields[0]; + const unsigned int last = *std::max_element(this->base_elements.compositional_fields.begin(), this->base_elements.compositional_fields.end()); + for (unsigned int i = first; i<=last; ++i) + result.emplace_back(i); + return result; + } + + + + template + std::vector + Introspection::get_compositional_field_indices_with_base_element(const unsigned int base_element_index) const + { + std::vector result; + Assert(base_element_index <= get_compositional_field_base_element_indices().back(), + ExcMessage("Invalid base_element_index specified.")); + + unsigned int idx = 0; + for (const auto base_idx : this->base_elements.compositional_fields) + { + if (base_idx == base_element_index) + result.emplace_back(idx); + + ++idx; + } + + Assert(result.size() > 0, ExcInternalError("There should be at least one compositional field for a valid base element.")); + return result; + } + + + template unsigned int Introspection::compositional_index_for_name (const std::string &name) const @@ -353,7 +393,7 @@ namespace aspect Introspection::name_for_compositional_index (const unsigned int index) const { // make sure that what we get here is really an index of one of the compositional fields - AssertIndexRange(index,composition_names.size()); + AssertIndexRange(index, composition_names.size()); return composition_names[index]; } @@ -488,6 +528,22 @@ namespace aspect } + + template + bool + Introspection::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.back() == n_components-1, ExcInternalError()); + + if (component_index >= component_indices.compositional_fields[0]) + return true; + else + return false; + } + } diff --git a/unit_tests/introspection.cc b/unit_tests/introspection.cc index 29c44c37490..7d7f0fd6131 100644 --- a/unit_tests/introspection.cc +++ b/unit_tests/introspection.cc @@ -67,4 +67,12 @@ 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 have the same FE: + CHECK(introspection.base_elements.compositional_fields == std::vector({3,3,3})); + CHECK(introspection.get_compositional_field_base_element_indices() == std::vector({3})); + CHECK(introspection.get_compositional_field_indices_with_base_element(3) == std::vector({0,1,2})); }