/**
* \file
* \copyright
* Copyright (c) 2012-2023, 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 "MaterialLib/Adsorption/Adsorption.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "NumLib/Function/Interpolation.h"
#include "TESLocalAssembler.h"
#include "TESReactionAdaptor.h"
namespace
{
enum class MatOutType
{
OGS5,
PYTHON
};
const MatOutType MATRIX_OUTPUT_FORMAT = MatOutType::PYTHON;
// TODO move to some location in the OGS core.
template <typename Mat>
void ogs5OutMat(const Mat& mat)
{
for (unsigned r = 0; r < mat.rows(); ++r)
{
switch (MATRIX_OUTPUT_FORMAT)
{
case MatOutType::OGS5:
if (r != 0)
{
std::printf("\n");
}
std::printf("|");
break;
case MatOutType::PYTHON:
if (r != 0)
{
std::printf(",\n");
}
std::printf("[");
break;
}
for (unsigned c = 0; c < mat.cols(); ++c)
{
switch (MATRIX_OUTPUT_FORMAT)
{
case MatOutType::OGS5:
std::printf(" %.16e", mat(r, c));
break;
case MatOutType::PYTHON:
if (c != 0)
{
std::printf(",");
}
std::printf(" %23.16g", mat(r, c));
break;
}
}
switch (MATRIX_OUTPUT_FORMAT)
{
case MatOutType::OGS5:
std::printf(" | ");
break;
case MatOutType::PYTHON:
std::printf(" ]");
break;
}
}
std::printf("\n");
}
template <typename Vec>
void ogs5OutVec(const Vec& vec)
{
for (unsigned r = 0; r < vec.size(); ++r)
{
switch (MATRIX_OUTPUT_FORMAT)
{
case MatOutType::OGS5:
if (r != 0)
{
std::printf("\n");
}
std::printf("| %.16e | ", vec[r]);
break;
case MatOutType::PYTHON:
if (r != 0)
{
std::printf(",\n");
}
std::printf("[ %23.16g ]", vec[r]);
break;
}
}
std::printf("\n");
}
} // anonymous namespace
namespace ProcessLib
{
namespace TES
{
template <typename ShapeFunction_, int GlobalDim>
TESLocalAssembler<ShapeFunction_, GlobalDim>::TESLocalAssembler(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
NumLib::GenericIntegrationMethod const& integration_method,
bool const is_axially_symmetric,
AssemblyParams const& asm_params)
: _element(e),
_integration_method(integration_method),
_shape_matrices(NumLib::initShapeMatrices<ShapeFunction,
ShapeMatricesType, GlobalDim>(
e, is_axially_symmetric, _integration_method)),
_d(asm_params, _integration_method.getNumberOfPoints(), GlobalDim)
{
}
template <typename ShapeFunction_, int GlobalDim>
void TESLocalAssembler<ShapeFunction_, GlobalDim>::assemble(
double const /*t*/, double const /*dt*/, std::vector<double> const& local_x,
std::vector<double> const& /*local_x_prev*/,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data)
{
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 * NODAL_DOF);
auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
local_K_data, local_matrix_size, local_matrix_size);
auto local_b = MathLib::createZeroedVector<NodalVectorType>(local_b_data,
local_matrix_size);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
_d.preEachAssemble();
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
auto const& sm = _shape_matrices[ip];
auto const& wp = _integration_method.getWeightedPoint(ip);
auto const weight = wp.getWeight();
_d.assembleIntegrationPoint(ip, local_x, sm, weight, local_M, local_K,
local_b);
}
if (_d.getAssemblyParameters().output_element_matrices)
{
std::puts("### Element: ?");
std::puts("---Velocity of water");
for (auto const& vs : _d.getData().velocity)
{
std::printf("| ");
for (auto v : vs)
{
std::printf("%23.16e ", v);
}
std::printf("|\n");
}
std::printf("\n---Mass matrix: \n");
ogs5OutMat(local_M);
std::printf("\n");
std::printf("---Laplacian + Advective + Content matrix: \n");
ogs5OutMat(local_K);
std::printf("\n");
std::printf("---RHS: \n");
ogs5OutVec(local_b);
std::printf("\n");
}
}
template <typename ShapeFunction_, int GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, GlobalDim>::getIntPtSolidDensity(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& /*cache*/) const
{
return _d.getData().solid_density;
}
template <typename ShapeFunction_, int GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, GlobalDim>::getIntPtLoading(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const
{
auto const rho_SR = _d.getData().solid_density;
auto const rho_SR_dry = _d.getAssemblyParameters().rho_SR_dry;
cache.clear();
cache.reserve(rho_SR.size());
transform(cbegin(rho_SR), cend(rho_SR), back_inserter(cache),
[&](auto const rho) { return rho / rho_SR_dry - 1.0; });
return cache;
}
template <typename ShapeFunction_, int GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, GlobalDim>::getIntPtReactionDampingFactor(
const double /*t*/, std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const
{
auto const fac = _d.getData().reaction_adaptor->getReactionDampingFactor();
auto const num_integration_points = _d.getData().solid_density.size();
cache.clear();
cache.resize(num_integration_points, fac);
return cache;
}
template <typename ShapeFunction_, int GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, GlobalDim>::getIntPtReactionRate(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& /*cache*/) const
{
return _d.getData().reaction_rate;
}
template <typename ShapeFunction_, int GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, GlobalDim>::getIntPtDarcyVelocity(
const double /*t*/,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const
{
auto const n_integration_points = _integration_method.getNumberOfPoints();
constexpr int process_id = 0; // monolithic scheme
auto const indices =
NumLib::getIndices(_element.getID(), *dof_table[process_id]);
assert(!indices.empty());
auto const local_x = x[process_id]->get(indices);
// local_x is ordered by component, local_x_mat is row major
auto const local_x_mat = MathLib::toMatrix<
Eigen::Matrix<double, NODAL_DOF, Eigen::Dynamic, Eigen::RowMajor>>(
local_x, NODAL_DOF, ShapeFunction_::NPOINTS);
cache.clear();
auto cache_mat = MathLib::createZeroedMatrix<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
cache, GlobalDim, n_integration_points);
for (unsigned i = 0; i < n_integration_points; ++i)
{
double p, T, x;
NumLib::shapeFunctionInterpolate(local_x, _shape_matrices[i].N, p, T,
x);
const double eta_GR = fluid_viscosity(p, T, x);
auto const& k = _d.getAssemblyParameters().solid_perm_tensor;
cache_mat.col(i).noalias() =
k * (_shape_matrices[i].dNdx *
local_x_mat.row(COMPONENT_ID_PRESSURE).transpose()) /
-eta_GR;
}
return cache;
}
template <typename ShapeFunction_, int GlobalDim>
bool TESLocalAssembler<ShapeFunction_, GlobalDim>::checkBounds(
std::vector<double> const& local_x,
std::vector<double> const& local_x_prev_ts)
{
return _d.getReactionAdaptor().checkBounds(local_x, local_x_prev_ts);
}
} // namespace TES
} // namespace ProcessLib