From 11d7b7b52312616e20489e5c5e575ca6f1f10c3c Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 30 Mar 2017 10:53:16 +0200
Subject: [PATCH] [PL/HC] Include computation of Darcy vel.
---
ProcessLib/HC/HCFEM.h | 52 +++++++++++++++++++++++++++++++++++++
ProcessLib/HC/HCProcess.cpp | 10 +++++++
ProcessLib/HC/HCProcess.h | 4 +++
3 files changed, 66 insertions(+)
diff --git a/ProcessLib/HC/HCFEM.h b/ProcessLib/HC/HCFEM.h
index 11c2c65aa7c..7e8a10f813d 100644
--- a/ProcessLib/HC/HCFEM.h
+++ b/ProcessLib/HC/HCFEM.h
@@ -257,6 +257,58 @@ public:
}
}
+ 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.viscosity_model->getValue(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_density->getValue(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
{
diff --git a/ProcessLib/HC/HCProcess.cpp b/ProcessLib/HC/HCProcess.cpp
index cb84dc12efb..ded20583b71 100644
--- a/ProcessLib/HC/HCProcess.cpp
+++ b/ProcessLib/HC/HCProcess.cpp
@@ -98,6 +98,16 @@ void HCProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, coupling_term);
}
+void HCProcess::computeSecondaryVariableConcrete(
+ double const t, GlobalVector const& x,
+ StaggeredCouplingTerm const& coupling_term)
+{
+ DBUG("Compute the Darcy velocity for HCProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &HCLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
+}
+
} // namespace HC
} // namespace ProcessLib
diff --git a/ProcessLib/HC/HCProcess.h b/ProcessLib/HC/HCProcess.h
index 4e72ad050a5..1addab3e2d1 100644
--- a/ProcessLib/HC/HCProcess.h
+++ b/ProcessLib/HC/HCProcess.h
@@ -57,6 +57,10 @@ public:
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,
--
GitLab