[PL] Interpolate process solutions to int_pts.

ProcessLib/ComponentTransport/ComponentTransportFEM.h

@@ -916,6 +916,21 @@ public:
         return flux;
     }
 
+    std::vector<double> interpolateNodalValuesToIntegrationPoints(
+        std::vector<double> const& local_x) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> interpolated_values(n_integration_points);
+        for (unsigned ip(0); ip < n_integration_points; ++ip)
+        {
+            NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N,
+                                             interpolated_values[ip]);
+        }
+        return interpolated_values;
+    }
+
     std::vector<double> const& getInterpolatedLocalSolution(
         const double /*t*/,
         std::vector<GlobalVector*> const& int_pt_x,

ProcessLib/ComponentTransport/ComponentTransportProcess.cpp

@@ -174,6 +174,77 @@ void ComponentTransportProcess::
         _local_assemblers, pv.getActiveElementIDs(), _coupled_solutions);
 }
 
+std::vector<GlobalVector>
+ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
+    std::vector<GlobalVector*> const& nodal_values_vectors) const
+{
+    // Result is for each process a vector of integration point values for each
+    // element stored consecutively.
+    auto interpolateNodalValuesToIntegrationPoints =
+        [](std::size_t mesh_item_id,
+           LocalAssemblerInterface& local_assembler,
+           std::vector<
+               std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+               dof_tables,
+           std::vector<GlobalVector*> const& x,
+           std::vector<std::vector<double>>& int_pt_x) {
+            for (unsigned process_id = 0; process_id < x.size(); ++process_id)
+            {
+                auto const& dof_table = dof_tables[process_id].get();
+                auto const indices =
+                    NumLib::getIndices(mesh_item_id, dof_table);
+                auto const local_x = x[process_id]->get(indices);
+
+                std::vector<double> const interpolated_values =
+                    local_assembler.interpolateNodalValuesToIntegrationPoints(
+                        local_x);
+
+                // For each element (mesh_item_id) concatenate the integration
+                // point values.
+                int_pt_x[process_id].insert(int_pt_x[process_id].end(),
+                                            interpolated_values.begin(),
+                                            interpolated_values.end());
+            }
+        };
+
+    // Same dof table for each primary variable.
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    dof_tables.reserve(nodal_values_vectors.size());
+    std::generate_n(std::back_inserter(dof_tables), nodal_values_vectors.size(),
+                    [&]() { return std::ref(*_local_to_global_index_map); });
+
+    std::vector<std::vector<double>> integration_point_values_vecs(
+        nodal_values_vectors.size());
+
+    GlobalExecutor::executeDereferenced(
+        interpolateNodalValuesToIntegrationPoints, _local_assemblers,
+        dof_tables, nodal_values_vectors, integration_point_values_vecs);
+
+    // Convert from std::vector<std::vector<double>> to
+    // std::vector<GlobalVector>
+    std::vector<GlobalVector> integration_point_values_vectors;
+    integration_point_values_vectors.reserve(nodal_values_vectors.size());
+    for (auto const& integration_point_values_vec :
+         integration_point_values_vecs)
+    {
+        GlobalIndexType const size = integration_point_values_vec.size();
+        // New vector of size (elements x integration_points_per_element)
+        GlobalVector integration_point_values_vector(size);
+
+        // Copy one by one.
+        for (GlobalIndexType i = 0; i < size; ++i)
+        {
+            integration_point_values_vector.set(
+                i, integration_point_values_vec[i]);
+        }
+        integration_point_values_vectors.push_back(
+            std::move(integration_point_values_vector));
+    }
+
+    return integration_point_values_vectors;
+}
+
 void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
     const double t,
     std::vector<GlobalVector*> const& integration_point_values_vectors,

ProcessLib/ComponentTransport/ComponentTransportProcess.h

@@ -117,6 +117,9 @@ public:
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
 
+    std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+        std::vector<GlobalVector*> const& nodal_values_vectors) const override;
+
     void extrapolateIntegrationPointValuesToNodes(
         const double t,
         std::vector<GlobalVector*> const& integration_point_values_vectors,

ProcessLib/LocalAssemblerInterface.h

@@ -98,6 +98,12 @@ public:
                              GlobalVector const& x, double const t,
                              double const dt, bool const use_monolithic_scheme);
 
+    virtual std::vector<double> interpolateNodalValuesToIntegrationPoints(
+        std::vector<double> const& /*local_x*/)
+    {
+        return {};
+    }
+
     /// Computes the flux in the point \c p_local_coords that is given in local
     /// coordinates using the values from \c local_x.
     /// Fits to monolithic scheme.

ProcessLib/Process.h

@@ -101,6 +101,12 @@ public:
     {
     }
 
+    virtual std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+        std::vector<GlobalVector*> const& /*nodal_values_vectors*/) const
+    {
+        return {};
+    }
+
     virtual void extrapolateIntegrationPointValuesToNodes(
         const double /*t*/,
         std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,

ProcessLib/TimeLoop.cpp

@@ -453,6 +453,11 @@ void TimeLoop::initialize()
     {
         BaseLib::RunTime time_phreeqc;
         time_phreeqc.start();
+
+        auto& pcs = _per_process_data[0]->process;
+        auto const interpolated_process_solutions =
+            pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
         _chemical_solver_interface->executeInitialCalculation(
             _process_solutions);
         INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
@@ -815,6 +820,11 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
         // space and localized chemical equilibrium between solutes.
         BaseLib::RunTime time_phreeqc;
         time_phreeqc.start();
+
+        auto& pcs = _per_process_data[0]->process;
+        auto const interpolated_process_solutions =
+            pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
         _chemical_solver_interface->doWaterChemistryCalculation(
             _process_solutions, dt);
         INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());