Commit 3ae1fb81 authored by Lars Bilke's avatar Lars Bilke
Browse files

Merge branch 'SmallFixes' into 'master'

Small fixes

See merge request ogs/ogs!3425
parents c096b07b 88447b6c
......@@ -239,7 +239,7 @@ MeshLib::Node* createNode(std::stringstream& sstr)
bool parseNodes(std::ifstream& in,
std::vector<MeshLib::Node*>& nodes,
std::map<std::size_t, std::size_t>& node_id_map,
MeshLib::Properties& mesh_prop)
MeshLib::Properties const& mesh_prop)
{
NodeType t = NodeType::UNSPECIFIED;
std::streampos pos = in.tellg();
......
......@@ -16,21 +16,21 @@
namespace MathLib
{
void applyKnownSolution(EigenMatrix &A_, EigenVector &b_, EigenVector &/*x*/,
void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &/*x*/,
const std::vector<EigenMatrix::IndexType> &vec_knownX_id,
const std::vector<double> &vec_knownX_x, double /*penalty_scaling*/)
{
using SpMat = EigenMatrix::RawMatrixType;
static_assert(SpMat::IsRowMajor, "matrix is assumed to be row major!");
auto &A = A_.getRawMatrix();
auto &b = b_.getRawVector();
auto &A_eigen = A.getRawMatrix();
auto &b_eigen = b.getRawVector();
// A(k, j) = 0.
// A_eigen(k, j) = 0.
// set row to zero
for (auto row_id : vec_knownX_id)
{
for (SpMat::InnerIterator it(A,row_id); it; ++it) {
for (SpMat::InnerIterator it(A_eigen,row_id); it; ++it) {
if (it.col() != decltype(it.col())(row_id))
{
it.valueRef() = 0.0;
......@@ -38,14 +38,14 @@ void applyKnownSolution(EigenMatrix &A_, EigenVector &b_, EigenVector &/*x*/,
}
}
SpMat AT = A.transpose();
SpMat AT = A_eigen.transpose();
for (std::size_t ix=0; ix<vec_knownX_id.size(); ix++)
{
SpMat::Index const row_id = vec_knownX_id[ix];
auto const x = vec_knownX_x[ix];
// b_i -= A(i,k)*val, i!=k
// b_i -= A_eigen(i,k)*val, i!=k
// set column to zero, subtract from rhs
for (SpMat::InnerIterator it(AT, row_id); it; ++it)
{
......@@ -54,20 +54,20 @@ void applyKnownSolution(EigenMatrix &A_, EigenVector &b_, EigenVector &/*x*/,
continue;
}
b[it.col()] -= it.value()*x;
b_eigen[it.col()] -= it.value()*x;
it.valueRef() = 0.0;
}
auto& c = AT.coeffRef(row_id, row_id);
if (c != 0.0) {
b[row_id] = x * c;
b_eigen[row_id] = x * c;
} else {
b[row_id] = x;
b_eigen[row_id] = x;
c = 1.0;
}
}
A = AT.transpose();
A_eigen = AT.transpose();
}
} // namespace MathLib
......@@ -24,19 +24,20 @@ class EigenVector;
*
* @param A Coefficient matrix
* @param b RHS vector
* @param _vec_knownX_id a vector of known solution entry IDs
* @param _vec_knownX_x a vector of known solutions
* @param vec_knownX_id a vector of known solution entry IDs
* @param vec_knownX_x a vector of known solutions
* @param penalty_scaling value for scaling some matrix and right hand side
* entries to enforce some conditions
* entries to enforce some conditions, value ignored in the current
* implementation
*/
void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &/*x*/,
const std::vector<EigenMatrix::IndexType> &_vec_knownX_id,
const std::vector<double> &_vec_knownX_x, double penalty_scaling = 1e+10);
const std::vector<EigenMatrix::IndexType> &vec_knownX_id,
const std::vector<double> &vec_knownX_x, double penalty_scaling = 1e+10);
inline
void applyKnownSolution(Eigen::MatrixXd const &A, Eigen::VectorXd &b, Eigen::VectorXd &/*x*/,
const std::vector<Eigen::MatrixXd::Index> &_vec_knownX_id,
const std::vector<double> &_vec_knownX_x, double penalty_scaling = 1e+10)
inline void applyKnownSolution(
Eigen::MatrixXd const& A, Eigen::VectorXd const& b, Eigen::VectorXd& /*x*/,
const std::vector<Eigen::MatrixXd::Index>& _vec_knownX_id,
const std::vector<double>& _vec_knownX_x, double penalty_scaling = 1e+10)
{
(void) A; (void) b; (void) _vec_knownX_id; (void) _vec_knownX_x;
(void) penalty_scaling;
......
......@@ -86,7 +86,7 @@ public:
/// Checks if a property vector with given type \c T, \c name, \c
/// mesh_item_type, and \c number_of_components exists.
template <typename T>
bool existsPropertyVector(std::string const& property_name,
bool existsPropertyVector(std::string const& name,
MeshItemType const mesh_item_type,
int const number_of_components) const;
......
......@@ -41,7 +41,8 @@ public:
_bhe(bhe),
_py_bc_object(py_bc_object)
{
const auto g_idx_T_out = static_cast<int>(in_out_global_indices.second);
const auto g_idx_T_out =
static_cast<int>(_in_out_global_indices.second);
// store the bc node ids to BHE network dataframe
std::get<3>(_py_bc_object.dataframe_network).emplace_back(g_idx_T_out);
......
......@@ -351,7 +351,6 @@ void Output::doOutputAlways(Process const& process,
DBUG("Found {:d} nodes for output at mesh '{:s}'.", nodes.size(),
mesh.getName());
MeshLib::MeshSubset mesh_subset(mesh, nodes);
std::vector<std::unique_ptr<NumLib::LocalToGlobalIndexMap>>
mesh_dof_tables;
mesh_dof_tables.reserve(x.size());
......@@ -359,7 +358,7 @@ void Output::doOutputAlways(Process const& process,
{
mesh_dof_tables.push_back(
process.getDOFTable(i).deriveBoundaryConstrainedMap(
std::move(mesh_subset)));
MeshLib::MeshSubset{mesh, nodes}));
}
std::vector<NumLib::LocalToGlobalIndexMap const*>
mesh_dof_table_pointers;
......
......@@ -12,10 +12,7 @@
#include <memory>
#include "MaterialLib/MPL/VariableType.h"
#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
#include "MathLib/KelvinVector.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
namespace ProcessLib
{
......
......@@ -15,7 +15,6 @@
#include "MeshLib/Elements/Utils.h"
#include "NumLib/DOF/ComputeSparsityPattern.h"
#include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/RichardsMechanics/CreateLocalAssemblers.h"
#include "RichardsMechanicsFEM.h"
#include "RichardsMechanicsProcessData.h"
......
......@@ -31,22 +31,22 @@ public:
PointGenerator point_generator1 = PointGenerator(
MathLib::Point3d(std::array<double, 3>{{0.0, 0.0, 0.0}}), 1.0);
SymmSegmentGenerator segment_generator1 = SymmSegmentGenerator{point_generator1,
std::bind(ac::reflect, point_generator1.center, std::placeholders::_1)};
[&](auto p){ return ac::reflect(point_generator1.center, p); }};
PointGenerator point_generator2 = PointGenerator(
MathLib::Point3d(std::array<double, 3>{{2.0, 0.0, 0.0}}), 1.0);
SymmSegmentGenerator segment_generator2 = SymmSegmentGenerator{point_generator2,
std::bind(ac::reflect, point_generator2.center, std::placeholders::_1)};
[&](auto p){ return ac::reflect(point_generator2.center, p); }};
Eigen::Vector3d const translation_vector1 = {2, 2, 0};
PairSegmentGenerator pair_segment_generator1 = PairSegmentGenerator{
segment_generator1,
std::bind(ac::translate, translation_vector1, std::placeholders::_1)};
[&](auto p) { return ac::translate(translation_vector1, p); }};
Eigen::Vector3d const translation_vector2 = {0, 0, 0};
PairSegmentGenerator pair_segment_generator2 = PairSegmentGenerator{
segment_generator1,
std::bind(ac::translate, translation_vector2, std::placeholders::_1)};
[&](auto p) { return ac::translate(translation_vector2, p); }};
ac::gtest_reporter gtest_reporter;
};
......
......@@ -24,8 +24,9 @@ public:
PointGenerator point_generator = PointGenerator(
MathLib::Point3d(std::array<double, 3>{{0.0, 0.0, 0.0}}), 1.0);
SymmSegmentGenerator segment_generator = SymmSegmentGenerator{point_generator,
std::bind(ac::reflect, point_generator.center, std::placeholders::_1)};
SymmSegmentGenerator segment_generator = SymmSegmentGenerator{
point_generator,
[&](auto p) { return ac::reflect(point_generator.center, p); }};
ac::gtest_reporter gtest_reporter;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment