Revision e2963065808465e78c566a62c31a67b267867bb6 authored by Dmitry Yu. Naumov on 19 October 2021, 09:43:49 UTC, committed by Dmitry Yu. Naumov on 19 October 2021, 09:43:49 UTC
[TH2M] Fixed a bug in using pressure, strain, stress dependent permeability models See merge request ogs/ogs!3821
HeuristicSearchLength.cpp
/**
* \file
* \date 2014-09-19
* \brief Implementation of heuristic search length strategy.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#include "HeuristicSearchLength.h"
#include "BaseLib/Logging.h"
#include "MeshLib/Elements/Element.h"
namespace MeshGeoToolsLib
{
HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh,
LengthType length_type)
: _mesh(mesh)
{
double sum(0.0);
double sum_of_sqr(0.0); // total length of edges
std::size_t n_sampling(0); // total length of edges squared
std::vector<MeshLib::Element*> const& elements(_mesh.getElements());
if (length_type == LengthType::Edge)
{
for (auto element : elements)
{
std::size_t const n_edges(element->getNumberOfEdges());
for (std::size_t k(0); k < n_edges; k++)
{
auto edge =
element->getEdge(k); // allocation inside getEdge().
double const len = edge->getContent();
delete edge;
sum += len;
sum_of_sqr += len * len;
}
n_sampling += n_edges;
}
}
else
{
for (const MeshLib::Element* e : elements)
{
auto const [min, max] =
MeshLib::computeSqrNodeDistanceRange(*e, true);
sum += std::sqrt(min);
sum_of_sqr += min;
}
n_sampling = _mesh.getNumberOfElements();
}
const double mean(sum / n_sampling);
const double variance((sum_of_sqr - (sum * sum) / n_sampling) /
(n_sampling - 1));
// Set the search length for the case of non-positive variance (which can
// happen due to numerics).
_search_length = mean / 2;
if (variance > 0)
{
if (variance < mean * mean / 4)
{
_search_length -= std::sqrt(variance);
}
else
{
_search_length = std::numeric_limits<double>::epsilon();
}
}
DBUG(
"[MeshNodeSearcher::MeshNodeSearcher] Calculated search length for "
"mesh '{:s}' is {:f}.",
_mesh.getName(), _search_length);
}
} // end namespace MeshGeoToolsLib
Computing file changes ...