https://gitlab.opengeosys.org/ogs/ogs.git
Tip revision: a8a2144fa41d836d81c2bcb0570897aec75788a6 authored by Dominik Kern on 25 March 2021, 10:52:55 UTC
use unit mD precisely and correct nondim time
use unit mD precisely and correct nondim time
Tip revision: a8a2144
MeshElementGrid.cpp
/*
* \brief Definition of the class MeshElementGrid.
*
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
*/
#include "MeshElementGrid.h"
#include <algorithm>
#include <bitset>
#include <cmath>
#include <memory>
#include "../Mesh.h"
#include "../Node.h"
#include "../Elements/Element.h"
#include "GeoLib/GEOObjects.h"
namespace MeshLib
{
MeshElementGrid::MeshElementGrid(MeshLib::Mesh const& mesh)
: _aabb{mesh.getNodes().cbegin(), mesh.getNodes().cend()},
_n_steps({{1, 1, 1}})
{
auto getDimensions =
[](MathLib::Point3d const& min, MathLib::Point3d const& max)
{
std::bitset<3> dim; // all bits are set to zero.
for (std::size_t k(0); k < 3; ++k) {
double const tolerance(
std::nexttoward(max[k],std::numeric_limits<double>::max())-max[k]);
if (std::abs(max[k] - min[k]) > tolerance)
{
dim[k] = true;
}
}
return dim;
};
MathLib::Point3d const& min_pnt(_aabb.getMinPoint());
MathLib::Point3d const& max_pnt(_aabb.getMaxPoint());
auto const dim = getDimensions(min_pnt, max_pnt);
std::array<double, 3> delta{{ max_pnt[0] - min_pnt[0],
max_pnt[1] - min_pnt[1], max_pnt[2] - min_pnt[2] }};
const std::size_t n_eles(mesh.getNumberOfElements());
const std::size_t n_eles_per_cell(100);
// *** condition: n_eles / n_cells < n_eles_per_cell
// where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
// *** with _n_steps[0] = ceil(pow(n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2]), 1/3.)));
// _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
// _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
auto sc_ceil = [](double v){
return static_cast<std::size_t>(std::ceil(v));
};
switch (dim.count()) {
case 3: // 3d case
_n_steps[0] = sc_ceil(std::cbrt(
n_eles*delta[0]*delta[0]/(n_eles_per_cell*delta[1]*delta[2])));
_n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
_n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
break;
case 2: // 2d cases
if (dim[0] && dim[2]) { // 2d case: xz plane, y = const
_n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[2])));
_n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]);
}
else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const
_n_steps[0] = sc_ceil(std::sqrt(n_eles*delta[0]/(n_eles_per_cell*delta[1])));
_n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]);
}
else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const
_n_steps[1] = sc_ceil(std::sqrt(n_eles*delta[1]/(n_eles_per_cell*delta[2])));
_n_steps[2] = sc_ceil(n_eles * delta[2] / (n_eles_per_cell*delta[1]));
}
break;
case 1: // 1d cases
for (std::size_t k(0); k<3; ++k) {
if (dim[k]) {
_n_steps[k] = sc_ceil(static_cast<double>(n_eles)/n_eles_per_cell);
}
}
}
// some frequently used expressions to fill the vector of elements per grid
// cell
for (std::size_t k(0); k<3; k++) {
_step_sizes[k] = delta[k] / _n_steps[k];
_inverse_step_sizes[k] = 1.0 / _step_sizes[k];
}
_elements_in_grid_box.resize(_n_steps[0]*_n_steps[1]*_n_steps[2]);
sortElementsInGridCells(mesh);
}
MathLib::Point3d const& MeshElementGrid::getMinPoint() const
{
return _aabb.getMinPoint();
}
MathLib::Point3d const& MeshElementGrid::getMaxPoint() const
{
return _aabb.getMaxPoint();
}
void MeshElementGrid::sortElementsInGridCells(MeshLib::Mesh const& mesh)
{
for (auto const element : mesh.getElements())
{
if (! sortElementInGridCells(*element)) {
OGS_FATAL("Sorting element (id={:d}) into mesh element grid.",
element->getID());
}
}
}
bool MeshElementGrid::sortElementInGridCells(MeshLib::Element const& element)
{
std::array<std::size_t,3> min{};
std::array<std::size_t,3> max{};
std::pair<bool, std::array<std::size_t, 3>> c(
getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(0)))));
if (c.first) {
min = c.second;
max = min;
} else {
return false;
}
for (std::size_t k(1); k<element.getNumberOfNodes(); ++k) {
// compute coordinates of the grid for each node of the element
c = getGridCellCoordinates(*(static_cast<MathLib::Point3d const*>(element.getNode(k))));
if (!c.first)
{
return false;
}
for (std::size_t j(0); j < 3; ++j)
{
if (min[j] > c.second[j])
{
min[j] = c.second[j];
}
if (max[j] < c.second[j])
{
max[j] = c.second[j];
}
}
}
const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
// If a node of an element is almost equal to the upper right point of the
// AABB the grid cell coordinates computed by getGridCellCoordintes() could
// be to large (due to numerical errors). The following lines ensure that
// the grid cell coordinates are in the valid range.
for (std::size_t k(0); k < 3; ++k)
{
max[k] = std::min(_n_steps[k] - 1, max[k]);
}
// insert the element into the grid cells
for (std::size_t i(min[0]); i<=max[0]; i++) {
for (std::size_t j(min[1]); j<=max[1]; j++) {
for (std::size_t k(min[2]); k<=max[2]; k++) {
_elements_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(&element);
}
}
}
return true;
}
std::pair<bool, std::array<std::size_t, 3>>
MeshElementGrid::getGridCellCoordinates(MathLib::Point3d const& p) const
{
bool valid(true);
std::array<std::size_t, 3> coords{};
for (std::size_t k(0); k<3; ++k) {
const double d(p[k]-_aabb.getMinPoint()[k]);
if (d < 0.0) {
valid = false;
coords[k] = 0;
} else if (_aabb.getMaxPoint()[k] <= p[k]) {
valid = false;
coords[k] = _n_steps[k]-1;
} else {
coords[k] = static_cast<std::size_t>(d * _inverse_step_sizes[k]);
}
}
return std::make_pair(valid, coords);
}
} // namespace MeshLib