From 38413ee720cb33e614421e8b88a147072bdecd72 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Tue, 8 May 2018 14:08:22 +0200
Subject: [PATCH] [PL] Implement getFlux() in HT process local asm.
---
ProcessLib/HT/HTFEM.h | 68 +++++++++++++++++++++++
ProcessLib/HT/HTLocalAssemblerInterface.h | 4 ++
ProcessLib/HT/MonolithicHTFEM.h | 53 ------------------
3 files changed, 72 insertions(+), 53 deletions(-)
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 27e67395b3e..3cc9cfca86f 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -88,6 +88,74 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
+ /// Computes the flux in the point \c pnt_local_coords that is given in
+ /// local coordinates using the values from \c local_x.
+ // TODO add time dependency
+ std::vector<double> getFlux(
+ MathLib::Point3d const& pnt_local_coords,
+ std::vector<double> const& local_x) const override
+ {
+ // eval dNdx and invJ at given point
+ using FemType =
+ NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+ FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+ &this->_element));
+
+ typename ShapeMatricesType::ShapeMatrices shape_matrices(
+ ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+ // Note: Axial symmetry is set to false here, because we only need dNdx
+ // here, which is not affected by axial symmetry.
+ fe.computeShapeFunctions(pnt_local_coords.getCoords(), shape_matrices,
+ GlobalDim, false);
+
+ // fetch permeability, viscosity, density
+ SpatialPosition pos;
+ pos.setElementID(this->_element.getID());
+
+ MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+ // local_x contains the local temperature and pressure values
+ NumLib::shapeFunctionInterpolate(
+ local_x, shape_matrices.N,
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)],
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::p)]);
+
+ // TODO remove following line if time dependency is implemented
+ double const t = 0.0;
+ auto const K =
+ this->_material_properties.porous_media_properties
+ .getIntrinsicPermeability(t, pos)
+ .getValue(t, pos, 0.0,
+ vars[static_cast<int>(
+ MaterialLib::Fluid::PropertyVariableType::T)]);
+
+ auto const mu = this->_material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+ GlobalDimMatrixType const K_over_mu = K / mu;
+
+ auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+ &local_x[local_x.size()/2], ShapeFunction::NPOINTS);
+ GlobalDimVectorType q =
+ -K_over_mu * shape_matrices.dNdx * p_nodal_values;
+
+ if (this->_material_properties.has_gravity)
+ {
+ auto const rho_w =
+ this->_material_properties.fluid_properties->getValue(
+ MaterialLib::Fluid::FluidPropertyType::Density, vars);
+ auto const b = this->_material_properties.specific_body_force;
+ q += K_over_mu * rho_w * b;
+ }
+ std::vector<double> flux;
+ flux.resize(3);
+ Eigen::Map<GlobalDimVectorType>(flux.data(), flux.size()) = q;
+
+ return flux;
+ }
+
protected:
MeshLib::Element const& _element;
HTMaterialProperties const& _material_properties;
diff --git a/ProcessLib/HT/HTLocalAssemblerInterface.h b/ProcessLib/HT/HTLocalAssemblerInterface.h
index 9aa542996e2..53117d69412 100644
--- a/ProcessLib/HT/HTLocalAssemblerInterface.h
+++ b/ProcessLib/HT/HTLocalAssemblerInterface.h
@@ -59,6 +59,10 @@ public:
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
+ virtual std::vector<double> getFlux(
+ MathLib::Point3d const& pnt_local_coords,
+ std::vector<double> const& local_x) const = 0;
+
protected:
// TODO: remove _coupled_solutions or move integration point data from
// local assembler class to a new class to make local assembler unique
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index ddbd94a5bca..33d7f6be096 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -184,59 +184,6 @@ public:
}
}
- /// Computes the flux in the point \c pnt_local_coords that is given in
- /// local coordinates using the values from \c local_x.
- // TODO add time dependency
- std::vector<double> getFlux(
- MathLib::Point3d const& pnt_local_coords,
- std::vector<double> const& local_x) const override
- {
- // eval dNdx and invJ at given point
- using FemType =
- NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
- FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
- &_element));
-
- typename ShapeMatricesType::ShapeMatrices shape_matrices(
- ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
-
- // Note: Axial symmetry is set to false here, because we only need dNdx
- // here, which is not affected by axial symmetry.
- fe.computeShapeFunctions(pnt_local_coords.getCoords(), shape_matrices,
- GlobalDim, false);
- std::vector<double> flux;
- flux.resize(3);
-
- // fetch permeability, viscosity, density
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- MaterialLib::Fluid::FluidProperty::ArrayType vars;
- double T_int_pt = 0.0;
- double p_int_pt = 0.0;
-
- // local_p, local_T?
- NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
- NumLib::shapeFunctionInterpolate(local_T, N, T_int_pt);
-
- // set the values into array vars
- vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =
- T_int_pt;
- vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
- p_int_pt;
-
- // TODO remove following line if time dependency is implemented
- double const t = 0.0;
- auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
-
- Eigen::Map<Eigen::RowVectorXd>(flux.data(), flux.size()) =
- -k * shape_matrices.dNdx *
- Eigen::Map<const Eigen::VectorXd>(local_x.data(), local_x.size());
-
- return flux;
- }
-
std::vector<double> const& getIntPtDarcyVelocity(
const double t,
GlobalVector const& current_solution,
--
GitLab