diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 60954b084caadd2969f0f314b7e9bba1f74c22de..a038048709e3d97481938d59a05ca7d014e48216 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -34,8 +34,8 @@
#include "ProcessLib/UncoupledProcessesTimeLoop.h"
-#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
#include "ProcessLib/ComponentTransport/CreateComponentTransportProcess.h"
+#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
#include "ProcessLib/HT/CreateHTProcess.h"
#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
#include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
@@ -43,6 +43,7 @@
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
#include "ProcessLib/PhaseField/CreatePhaseFieldProcess.h"
+#include "ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h"
#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
@@ -420,6 +421,14 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
break;
}
}
+ else if (type == "RichardsComponentTransport")
+ {
+ process = ProcessLib::RichardsComponentTransport::
+ createRichardsComponentTransportProcess(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
+ }
else if (type == "SMALL_DEFORMATION")
{
switch (_mesh_vec[0]->getDimension())
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index aa8acef9f7287e25d49f6a699b533136a99b9d1b..5d0e0d3d03a805efb36f2c498edd2c2512f5db32 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -4,6 +4,7 @@ APPEND_SOURCE_FILES(SOURCES)
APPEND_SOURCE_FILES(SOURCES BoundaryCondition)
APPEND_SOURCE_FILES(SOURCES CalculateSurfaceFlux)
APPEND_SOURCE_FILES(SOURCES ComponentTransport)
+APPEND_SOURCE_FILES(SOURCES RichardsComponentTransport)
APPEND_SOURCE_FILES(SOURCES Deformation)
APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
APPEND_SOURCE_FILES(SOURCES HeatConduction)
diff --git a/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb65c752e570e541fde15b981ea2d1fa3c74490e
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.cpp
@@ -0,0 +1,89 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreatePorousMediaProperties.h"
+
+#include "BaseLib/reorderVector.h"
+
+#include "MaterialLib/PorousMedium/Permeability/createPermeabilityModel.h"
+#include "MaterialLib/PorousMedium/Porosity/createPorosityModel.h"
+#include "MaterialLib/PorousMedium/Storage/createStorageModel.h"
+
+#include "MeshLib/Mesh.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+PorousMediaProperties createPorousMediaProperties(
+ MeshLib::Mesh& mesh, BaseLib::ConfigTree const& porous_medium_configs)
+{
+ DBUG("Create PorousMediaProperties.");
+
+ std::vector<Eigen::MatrixXd> intrinsic_permeability_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+ porosity_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+ storage_models;
+
+ std::vector<int> mat_ids;
+ for (auto const& porous_medium_config :
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium}
+ porous_medium_configs.getConfigSubtreeList("porous_medium"))
+ {
+ //! \ogs_file_attr{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__id}
+ auto const id = porous_medium_config.getConfigAttribute<int>("id");
+ mat_ids.push_back(id);
+
+ auto const& porosity_config =
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__porosity}
+ porous_medium_config.getConfigSubtree("porosity");
+ porosity_models.emplace_back(
+ MaterialLib::PorousMedium::createPorosityModel(porosity_config));
+
+ // Configuration for the intrinsic permeability
+ auto const& permeability_config =
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__permeability}
+ porous_medium_config.getConfigSubtree("permeability");
+ intrinsic_permeability_models.emplace_back(
+ MaterialLib::PorousMedium::createPermeabilityModel(
+ permeability_config));
+
+ // Configuration for the specific storage.
+ auto const& storage_config =
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium__porous_medium__storage}
+ porous_medium_config.getConfigSubtree("storage");
+ storage_models.emplace_back(
+ MaterialLib::PorousMedium::createStorageModel(storage_config));
+ }
+
+ BaseLib::reorderVector(intrinsic_permeability_models, mat_ids);
+ BaseLib::reorderVector(porosity_models, mat_ids);
+ BaseLib::reorderVector(storage_models, mat_ids);
+
+ std::vector<int> material_ids(mesh.getNumberOfElements());
+ if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
+ {
+ auto const& mesh_material_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+ material_ids.reserve(mesh_material_ids->size());
+ std::copy(mesh_material_ids->cbegin(), mesh_material_ids->cend(),
+ material_ids.begin());
+ }
+
+ return PorousMediaProperties{std::move(porosity_models),
+ std::move(intrinsic_permeability_models),
+ std::move(storage_models),
+ std::move(material_ids)};
+}
+
+} // namespace ComponentTransport
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..73e0ece1062340d303ccfb5839b813c8ebfe9e07
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreatePorousMediaProperties.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "PorousMediaProperties.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+PorousMediaProperties createPorousMediaProperties(
+ MeshLib::Mesh& mesh, BaseLib::ConfigTree const& porous_media_config);
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..50a3377cd626089d55f2baebe075339a42309a3f
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
@@ -0,0 +1,153 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateRichardsComponentTransportProcess.h"
+
+#include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
+
+#include "ProcessLib/Parameter/ConstantParameter.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "CreatePorousMediaProperties.h"
+#include "RichardsComponentTransportProcess.h"
+#include "RichardsComponentTransportProcessData.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+std::unique_ptr<Process> createRichardsComponentTransportProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{prj__processes__process__type}
+ config.checkConfigParameter("type", "RichardsComponentTransport");
+
+ DBUG("Create RichardsComponentTransportProcess.");
+
+ // Process variable.
+
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ auto process_variables = findProcessVariables(
+ variables, pv_config,
+ {
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__process_variables__concentration}
+ "concentration",
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__process_variables__pressure}
+ "pressure"});
+
+ auto const& porous_medium_configs =
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__porous_medium}
+ config.getConfigSubtree("porous_medium");
+ PorousMediaProperties porous_media_properties{
+ createPorousMediaProperties(mesh, porous_medium_configs)};
+
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+
+ auto fluid_properties =
+ MaterialLib::Fluid::createFluidProperties(fluid_config);
+
+ // Parameter for the density of the fluid.
+ auto& fluid_reference_density= findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__fluid_reference_density}
+ "fluid_reference_density", parameters, 1);
+ DBUG("Use \'%s\' as fluid_reference_density parameter.",
+ fluid_reference_density.name.c_str());
+
+ // Parameter for the longitudinal solute dispersivity.
+ auto const& molecular_diffusion_coefficient = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__molecular_diffusion_coefficient}
+ "molecular_diffusion_coefficient", parameters, 1);
+ DBUG("Use \'%s\' as molecular diffusion coefficient parameter.",
+ molecular_diffusion_coefficient.name.c_str());
+
+ // Parameter for the longitudinal solute dispersivity.
+ auto const& solute_dispersivity_longitudinal = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__solute_dispersivity_longitudinal}
+ "solute_dispersivity_longitudinal", parameters, 1);
+ DBUG("Use \'%s\' as longitudinal solute dispersivity parameter.",
+ solute_dispersivity_longitudinal.name.c_str());
+
+ // Parameter for the transverse solute dispersivity.
+ auto const& solute_dispersivity_transverse = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__solute_dispersivity_transverse}
+ "solute_dispersivity_transverse", parameters, 1);
+ DBUG("Use \'%s\' as transverse solute dispersivity parameter.",
+ solute_dispersivity_transverse.name.c_str());
+
+ // Parameter for the retardation factor.
+ auto const& retardation_factor = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__retardation_factor}
+ "retardation_factor", parameters, 1);
+
+ // Parameter for the decay rate.
+ auto const& decay_rate = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__RichardsComponentTransport__decay_rate}
+ "decay_rate", parameters, 1);
+
+ // Specific body force parameter.
+ Eigen::VectorXd specific_body_force;
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__RichardsComponentTransport__specific_body_force}
+ config.getConfigParameter<std::vector<double>>("specific_body_force");
+ assert(b.size() > 0 && b.size() < 4);
+ if (b.size() < mesh.getDimension())
+ OGS_FATAL(
+ "specific body force (gravity vector) has %d components, mesh "
+ "dimension is %d",
+ b.size(), mesh.getDimension());
+ bool const has_gravity = MathLib::toVector(b).norm() > 0;
+ if (has_gravity)
+ {
+ specific_body_force.resize(b.size());
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+ }
+
+ RichardsComponentTransportProcessData process_data{
+ std::move(porous_media_properties),
+ fluid_reference_density,
+ std::move(fluid_properties),
+ molecular_diffusion_coefficient,
+ solute_dispersivity_longitudinal,
+ solute_dispersivity_transverse,
+ retardation_factor,
+ decay_rate,
+ specific_body_force,
+ has_gravity};
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"RichardsComponentTransport_concentration_pressure"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ return std::make_unique<RichardsComponentTransportProcess>(
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller));
+}
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..539b88eea044cd36b74b1d419e71e0d32b320bb7
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+std::unique_ptr<Process> createRichardsComponentTransportProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c89ce7c999ee640201bff1271a45bfe28a9c8679
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.cpp
@@ -0,0 +1,57 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "PorousMediaProperties.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+int PorousMediaProperties::getMaterialID(SpatialPosition const& pos) const
+{
+ int const element_id = pos.getElementID().get();
+ return _material_ids[element_id];
+}
+
+MaterialLib::PorousMedium::Porosity const& PorousMediaProperties::getPorosity(
+ double /*t*/, SpatialPosition const& pos) const
+{
+ return *_porosity_models[getMaterialID(pos)];
+}
+
+Eigen::MatrixXd const& PorousMediaProperties::getIntrinsicPermeability(
+ double /*t*/, SpatialPosition const& pos) const
+{
+ return _intrinsic_permeability_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::Storage const&
+PorousMediaProperties::getSpecificStorage(double /*t*/,
+ SpatialPosition const& pos) const
+{
+ return *_specific_storage_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::CapillaryPressureSaturation const&
+PorousMediaProperties::getCapillaryPressureSaturationModel(
+ double /*t*/, SpatialPosition const& pos) const
+{
+ return *_capillary_pressure_saturation_models[getMaterialID(pos)];
+}
+
+MaterialLib::PorousMedium::RelativePermeability const&
+PorousMediaProperties::getRelativePermeability(
+ double /*t*/, SpatialPosition const& pos) const
+{
+ return *_relative_permeability_models[getMaterialID(pos)];
+}
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..12dd178163fb0fa5dc5a11e048548c9e59778308
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/PorousMediaProperties.h
@@ -0,0 +1,77 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+#include <Eigen/Dense>
+#include "BaseLib/reorderVector.h"
+
+
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+
+class PorousMediaProperties
+{
+public:
+ PorousMediaProperties(
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&&
+ porosity_models,
+ std::vector<Eigen::MatrixXd>&& intrinsic_permeability_models,
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&&
+ specific_storage_models,
+ std::vector<int>&& material_ids)
+ : _porosity_models(std::move(porosity_models)),
+ _intrinsic_permeability_models(
+ std::move(intrinsic_permeability_models)),
+ _specific_storage_models(std::move(specific_storage_models)),
+ _material_ids(std::move(material_ids))
+ {
+ }
+
+ PorousMediaProperties(PorousMediaProperties&& other)
+ : _porosity_models(std::move(other._porosity_models)),
+ _intrinsic_permeability_models(other._intrinsic_permeability_models),
+ _specific_storage_models(std::move(other._specific_storage_models)),
+ _material_ids(other._material_ids)
+ {
+ }
+
+ MaterialLib::PorousMedium::Porosity const& getPorosity(
+ double t, SpatialPosition const& pos) const;
+
+ Eigen::MatrixXd const& getIntrinsicPermeability(
+ double t, SpatialPosition const& pos) const;
+
+ MaterialLib::PorousMedium::Storage const& getSpecificStorage(
+ double t, SpatialPosition const& pos) const;
+
+private:
+ int getMaterialID(SpatialPosition const& pos) const;
+private:
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+ _porosity_models;
+ std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+ _specific_storage_models;
+ std::vector<int> _material_ids;
+};
+
+}
+}
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..c23ac9fb2c6f3d9d82c87f7349ea27a92f6eda47
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
@@ -0,0 +1,379 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+
+#include "RichardsComponentTransportProcessData.h"
+#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
+struct IntegrationPointData final
+{
+ IntegrationPointData(NodalRowVectorType const& N_,
+ GlobalDimNodalMatrixType const& dNdx_,
+ double const& integration_weight_)
+ : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
+ {
+ }
+
+ NodalRowVectorType const N;
+ GlobalDimNodalMatrixType const dNdx;
+ double const integration_weight;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+const unsigned NUM_NODAL_DOF = 2;
+
+class RichardsComponentTransportLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class LocalAssemblerData
+ : public RichardsComponentTransportLocalAssemblerInterface
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalMatrixType = typename ShapeMatricesType::template MatrixType<
+ NUM_NODAL_DOF * ShapeFunction::NPOINTS,
+ NUM_NODAL_DOF * ShapeFunction::NPOINTS>;
+ using LocalVectorType =
+ typename ShapeMatricesType::template VectorType<NUM_NODAL_DOF *
+ ShapeFunction::NPOINTS>;
+
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using GlobalDimNodalMatrixType =
+ typename ShapeMatricesType::GlobalDimNodalMatrixType;
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+
+public:
+ LocalAssemblerData(
+ MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ RichardsComponentTransportProcessData const& process_data)
+ : _element(element),
+ _process_data(process_data),
+ _integration_method(integration_order),
+ _darcy_velocities(
+ GlobalDim,
+ std::vector<double>(_integration_method.getNumberOfPoints()))
+ {
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+ (void)local_matrix_size;
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ _ip_data.reserve(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data.emplace_back(
+ shape_matrices[ip].N, shape_matrices[ip].dNdx,
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ shape_matrices[ip].integralMeasure *
+ shape_matrices[ip].detJ);
+ }
+ }
+
+ void assemble(double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data) override
+ {
+ auto const local_matrix_size = local_x.size();
+ // This assertion is valid only if all nodal d.o.f. use the same shape
+ // matrices.
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<LocalVectorType>(
+ local_b_data, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ auto const num_nodes = ShapeFunction::NPOINTS;
+ auto p_nodal_values =
+ Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
+
+ auto const& b = _process_data.specific_body_force;
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ GlobalDimMatrixType const& I(
+ GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
+
+ auto KCC = local_K.template block<num_nodes, num_nodes>(0, 0);
+ auto MCC = local_M.template block<num_nodes, num_nodes>(0, 0);
+ auto Kpp =
+ local_K.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Mpp =
+ local_M.template block<num_nodes, num_nodes>(num_nodes, num_nodes);
+ auto Bp = local_b.template block<num_nodes, 1>(num_nodes, 0);
+
+ for (std::size_t ip(0); ip < n_integration_points; ++ip)
+ {
+ pos.setIntegrationPoint(ip);
+
+ // \todo the argument to getValue() has to be changed for non
+ // constant storage model
+ auto const specific_storage =
+ _process_data.porous_media_properties.getSpecificStorage(t, pos)
+ .getValue(0.0);
+
+ auto const& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+ auto const& w = ip_data.integration_weight;
+
+ double C_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ // Order matters: First C, then p!
+ NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
+
+ // \todo the first argument has to be changed for non constant
+ // porosity model
+ auto const porosity =
+ _process_data.porous_media_properties.getPorosity(t, pos)
+ .getValue(0.0, C_int_pt);
+
+ auto const retardation_factor =
+ _process_data.retardation_factor(t, pos)[0];
+
+ auto const& solute_dispersivity_transverse =
+ _process_data.solute_dispersivity_transverse(t, pos)[0];
+ auto const& solute_dispersivity_longitudinal =
+ _process_data.solute_dispersivity_longitudinal(t, pos)[0];
+
+ // Use the fluid density model to compute the density
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+ auto const density = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const& decay_rate = _process_data.decay_rate(t, pos)[0];
+ auto const& molecular_diffusion_coefficient =
+ _process_data.molecular_diffusion_coefficient(t, pos)[0];
+
+ auto const& K =
+ _process_data.porous_media_properties.getIntrinsicPermeability(
+ t, pos);
+ // Use the viscosity model to compute the viscosity
+ auto const mu = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ GlobalDimVectorType const velocity =
+ _process_data.has_gravity
+ ? GlobalDimVectorType(-K_over_mu *
+ (dNdx * p_nodal_values - density * b))
+ : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
+
+ double const velocity_magnitude = velocity.norm();
+ GlobalDimMatrixType const hydrodynamic_dispersion =
+ velocity_magnitude != 0.0
+ ? GlobalDimMatrixType(
+ (porosity * molecular_diffusion_coefficient +
+ solute_dispersivity_transverse *
+ velocity_magnitude) *
+ I +
+ (solute_dispersivity_longitudinal -
+ solute_dispersivity_transverse) /
+ velocity_magnitude * velocity *
+ velocity.transpose())
+ : GlobalDimMatrixType(
+ (porosity * molecular_diffusion_coefficient +
+ solute_dispersivity_transverse *
+ velocity_magnitude) *
+ I);
+
+ // matrix assembly
+ KCC.noalias() +=
+ (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
+ N.transpose() * velocity.transpose() * dNdx +
+ N.transpose() * decay_rate * porosity * retardation_factor *
+ N) *
+ w;
+ MCC.noalias() +=
+ w * N.transpose() * porosity * retardation_factor * N;
+ Kpp.noalias() += w * dNdx.transpose() * K_over_mu * dNdx;
+ Mpp.noalias() += w * N.transpose() * specific_storage * N;
+ if (_process_data.has_gravity)
+ Bp += w * density * dNdx.transpose() * K_over_mu * b;
+ /* with Oberbeck-Boussing assumption density difference only exists
+ * in buoyancy effects */
+ }
+ }
+
+ void computeSecondaryVariableConcrete(
+ double const t, std::vector<double> const& local_x) override
+ {
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ auto const& K =
+ _process_data.porous_media_properties.getIntrinsicPermeability(t,
+ pos);
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ auto const mu = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+ &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ auto const& ip_data = _ip_data[ip];
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+
+ GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+ if (_process_data.has_gravity)
+ {
+ double C_int_pt = 0.0;
+ double p_int_pt = 0.0;
+ NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt,
+ p_int_pt);
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+
+ auto const rho_w = _process_data.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const b = _process_data.specific_body_force;
+ // here it is assumed that the vector b is directed 'downwards'
+ velocity += K_over_mu * rho_w * b;
+ }
+
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ _darcy_velocities[d][ip] = velocity[d];
+ }
+ }
+ }
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _ip_data[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 0);
+ return _darcy_velocities[0];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 1);
+ return _darcy_velocities[1];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 2);
+ return _darcy_velocities[2];
+ }
+
+private:
+ MeshLib::Element const& _element;
+ RichardsComponentTransportProcessData const& _process_data;
+
+ IntegrationMethod const _integration_method;
+ std::vector<
+ IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>,
+ Eigen::aligned_allocator<
+ IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
+ _ip_data;
+ std::vector<std::vector<double>> _darcy_velocities;
+};
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1310718ea9145adf891b4a2838c22727b22f117d
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
@@ -0,0 +1,115 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "RichardsComponentTransportProcess.h"
+
+#include <cassert>
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+RichardsComponentTransportProcess::RichardsComponentTransportProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ RichardsComponentTransportProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), std::move(named_function_caller)),
+ _process_data(std::move(process_data))
+{
+}
+
+void RichardsComponentTransportProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
+ ProcessLib::createLocalAssemblers<LocalAssemblerData>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_x",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
+ &RichardsComponentTransportLocalAssemblerInterface::
+ getIntPtDarcyVelocityX));
+
+ if (mesh.getDimension() > 1)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_y",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &RichardsComponentTransportLocalAssemblerInterface::
+ getIntPtDarcyVelocityY));
+ }
+ if (mesh.getDimension() > 2)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_z",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &RichardsComponentTransportLocalAssemblerInterface::
+ getIntPtDarcyVelocityZ));
+ }
+}
+
+void RichardsComponentTransportProcess::assembleConcreteProcess(
+ const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble RichardsComponentTransportProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ *_local_to_global_index_map, t, x, M, K, b, _coupling_term);
+}
+
+void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac)
+{
+ DBUG("AssembleWithJacobian RichardsComponentTransportProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+ dx_dx, M, K, b, Jac, _coupling_term);
+}
+
+void RichardsComponentTransportProcess::computeSecondaryVariableConcrete(
+ double const t, GlobalVector const& x,
+ StaggeredCouplingTerm const& coupling_term)
+{
+ DBUG("Compute the Darcy velocity for RichardsComponentTransportProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &RichardsComponentTransportLocalAssemblerInterface::
+ computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
+}
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib
+
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..d1967fea7d7869596b251c893f82a178c609854b
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
@@ -0,0 +1,88 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "RichardsComponentTransportFEM.h"
+#include "RichardsComponentTransportProcessData.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace RichardsComponentTransport
+{
+/**
+ * # RichardsComponentTransport process
+ *
+ * The implementation uses a monolithic approach, i.e., both processes
+ * are assembled within one global system of equations.
+ *
+ * ## Process Coupling
+ *
+ * The advective term of the concentration equation is given by the Richards
+ * flow process, i.e., the concentration distribution depends on
+ * darcy velocity of the Richards flow process. On the other hand the
+ * concentration dependencies of the viscosity and density in the groundwater
+ * flow couples the unsaturated H process to the C process.
+ *
+ * \note At the moment there is not any coupling by source or sink terms, i.e.,
+ * the coupling is implemented only by density changes due to concentration
+ * changes in the buoyance term in the groundwater flow. This coupling schema is
+ * referred to as the Boussinesq approximation.
+ * */
+class RichardsComponentTransportProcess final : public Process
+{
+public:
+ RichardsComponentTransportProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
+ jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&&
+ process_variables,
+ RichardsComponentTransportProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller);
+
+ //! \name ODESystem interface
+ //! @{
+
+ bool isLinear() const override { return false; }
+ //! @}
+
+ void computeSecondaryVariableConcrete(
+ double const t, GlobalVector const& x,
+ StaggeredCouplingTerm const& coupling_term) override;
+
+private:
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override;
+
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+ RichardsComponentTransportProcessData _process_data;
+
+ std::vector<
+ std::unique_ptr<RichardsComponentTransportLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..97c4c58a803b9714ecd38469f84cdf201a4036c8
--- /dev/null
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcessData.h
@@ -0,0 +1,94 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+#include "PorousMediaProperties.h"
+
+namespace ProcessLib
+{
+template <typename ReturnType>
+struct Parameter;
+
+namespace RichardsComponentTransport
+{
+struct RichardsComponentTransportProcessData
+{
+ RichardsComponentTransportProcessData(
+ PorousMediaProperties&& porous_media_properties_,
+ ProcessLib::Parameter<double> const& fluid_reference_density_,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperties>&&
+ fluid_properties_,
+ ProcessLib::Parameter<double> const& molecular_diffusion_coefficient_,
+ ProcessLib::Parameter<double> const& solute_dispersivity_longitudinal_,
+ ProcessLib::Parameter<double> const& solute_dispersivity_transverse_,
+ ProcessLib::Parameter<double> const& retardation_factor_,
+ ProcessLib::Parameter<double> const& decay_rate_,
+ Eigen::VectorXd const& specific_body_force_,
+ bool const has_gravity_)
+ : porous_media_properties(std::move(porous_media_properties_)),
+ fluid_reference_density(fluid_reference_density_),
+ fluid_properties(std::move(fluid_properties_)),
+ molecular_diffusion_coefficient(molecular_diffusion_coefficient_),
+ solute_dispersivity_longitudinal(solute_dispersivity_longitudinal_),
+ solute_dispersivity_transverse(solute_dispersivity_transverse_),
+ retardation_factor(retardation_factor_),
+ decay_rate(decay_rate_),
+ specific_body_force(specific_body_force_),
+ has_gravity(has_gravity_)
+ {
+ }
+
+ RichardsComponentTransportProcessData(
+ RichardsComponentTransportProcessData&& other)
+ : porous_media_properties(std::move(other.porous_media_properties)),
+ fluid_reference_density(other.fluid_reference_density),
+ fluid_properties(other.fluid_properties.release()),
+ molecular_diffusion_coefficient(
+ other.molecular_diffusion_coefficient),
+ solute_dispersivity_longitudinal(
+ other.solute_dispersivity_longitudinal),
+ solute_dispersivity_transverse(other.solute_dispersivity_transverse),
+ retardation_factor(other.retardation_factor),
+ decay_rate(other.decay_rate),
+ specific_body_force(other.specific_body_force),
+ has_gravity(other.has_gravity)
+ {
+ }
+
+ //! Copies are forbidden.
+ RichardsComponentTransportProcessData(
+ RichardsComponentTransportProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(RichardsComponentTransportProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(RichardsComponentTransportProcessData&&) = delete;
+
+ PorousMediaProperties porous_media_properties;
+ Parameter<double> const& fluid_reference_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperties> fluid_properties;
+ Parameter<double> const& molecular_diffusion_coefficient;
+ Parameter<double> const& solute_dispersivity_longitudinal;
+ Parameter<double> const& solute_dispersivity_transverse;
+ Parameter<double> const& retardation_factor;
+ Parameter<double> const& decay_rate;
+ Eigen::VectorXd const specific_body_force;
+ bool const has_gravity;
+};
+
+} // namespace RichardsComponentTransport
+} // namespace ProcessLib