From 419409bd0b13317173c9ae130867fd4aed41016e Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 30 Mar 2017 09:26:35 +0200
Subject: [PATCH] [PL/HC] Change handling of body force vector.
---
ProcessLib/HC/CreateHCProcess.cpp | 10 +++++++++-
ProcessLib/HC/HCFEM.h | 12 +++++++++---
ProcessLib/HC/HCProcessData.h | 4 ++--
3 files changed, 20 insertions(+), 6 deletions(-)
diff --git a/ProcessLib/HC/CreateHCProcess.cpp b/ProcessLib/HC/CreateHCProcess.cpp
index 078fcfee26c..ac617993cdc 100644
--- a/ProcessLib/HC/CreateHCProcess.cpp
+++ b/ProcessLib/HC/CreateHCProcess.cpp
@@ -116,14 +116,22 @@ std::unique_ptr<Process> createHCProcess(
"decay_rate", parameters, 1);
// Specific body force parameter.
- Eigen::Vector3d specific_body_force;
+ Eigen::VectorXd specific_body_force;
std::vector<double> const b =
//! \ogs_file_param{prj__processes__process__HC__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());
+ }
HCProcessData process_data{
std::move(porous_media_properties),
diff --git a/ProcessLib/HC/HCFEM.h b/ProcessLib/HC/HCFEM.h
index 6464581e7b9..11c2c65aa7c 100644
--- a/ProcessLib/HC/HCFEM.h
+++ b/ProcessLib/HC/HCFEM.h
@@ -144,7 +144,7 @@ public:
auto p_nodal_values =
Eigen::Map<const NodalVectorType>(&local_x[num_nodes], num_nodes);
- auto const & b = _process_data.specific_body_force.head(GlobalDim);
+ auto const& b = _process_data.specific_body_force;
MaterialLib::Fluid::FluidProperty::ArrayType vars;
@@ -212,7 +212,12 @@ public:
GlobalDimMatrixType const K_over_mu = K / mu;
GlobalDimVectorType const velocity =
- -perm_over_visc * (dNdx * p_nodal_values - density_water * b);
+ _process_data.has_gravity
+ ? GlobalDimVectorType(
+ -perm_over_visc *
+ (dNdx * p_nodal_values - density_water * b))
+ : GlobalDimVectorType(-perm_over_visc * dNdx *
+ p_nodal_values);
double const velocity_magnitude = velocity.norm();
GlobalDimMatrixType const& I(
@@ -245,7 +250,8 @@ public:
w * N.transpose() * porosity * retardation_factor * N;
Kpp.noalias() += w * dNdx.transpose() * perm_over_visc * dNdx;
Mpp.noalias() += w * N.transpose() * specific_storage * N;
- Bp += w * density_water * dNdx.transpose() * perm_over_visc * b;
+ if (_process_data.has_gravity)
+ Bp += w * density_water * dNdx.transpose() * perm_over_visc * b;
/* with Oberbeck-Boussing assumption density difference only exists
* in buoyancy effects */
}
diff --git a/ProcessLib/HC/HCProcessData.h b/ProcessLib/HC/HCProcessData.h
index d92119210ff..f83c9c55ccf 100644
--- a/ProcessLib/HC/HCProcessData.h
+++ b/ProcessLib/HC/HCProcessData.h
@@ -41,7 +41,7 @@ struct HCProcessData
ProcessLib::Parameter<double> const& solute_dispersivity_transverse_,
ProcessLib::Parameter<double> const& retardation_factor_,
ProcessLib::Parameter<double> const& decay_rate_,
- Eigen::Vector3d const& specific_body_force_,
+ Eigen::VectorXd const& specific_body_force_,
bool const has_gravity_)
: porous_media_properties(std::move(porous_media_properties_)),
viscosity_model(std::move(viscosity_model_)),
@@ -92,7 +92,7 @@ struct HCProcessData
Parameter<double> const& solute_dispersivity_transverse;
Parameter<double> const& retardation_factor;
Parameter<double> const& decay_rate;
- Eigen::Vector3d const specific_body_force;
+ Eigen::VectorXd const specific_body_force;
bool const has_gravity;
};
--
GitLab