/* This file is part of the OpenLB library
*
* Copyright (C) 2013, 2014 Mathias J. Krause
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
*
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
/** \file
* Representation of a statistic for a parallel 3D geometry -- generic implementation.
*/
#include
#include
#include
#include
//#define PARALLEL_MODE_MPI
#include "geometry/superGeometry3D.h"
#include "geometry/superGeometryStatistics3D.h"
#include "utilities/vectorHelpers.h"
#include "core/vector.h"
using namespace olb::util;
namespace olb {
template
SuperGeometryStatistics3D::SuperGeometryStatistics3D(SuperGeometry3D* superGeometry)
: _superGeometry(superGeometry), _statisticsUpdateNeeded(true),
_overlap(superGeometry->getOverlap()), _nMaterials(0), _minOverMaterial(3,T(0)), _maxOverMaterial(3,T(0)),
clout(std::cout,"SuperGeometryStatistics3D")
{
}
template
SuperGeometryStatistics3D::SuperGeometryStatistics3D(SuperGeometryStatistics3D const& rhs)
: _superGeometry(rhs._superGeometry), _statisticsUpdateNeeded(true),
_overlap(rhs._superGeometry->getOverlap() ), _nMaterials(rhs._nMaterials),
_minOverMaterial(rhs._minOverMaterial), _maxOverMaterial(rhs._maxOverMaterial),
clout(std::cout,"SuperGeometryStatistics3D")
{
}
template
SuperGeometryStatistics3D& SuperGeometryStatistics3D::operator=(SuperGeometryStatistics3D const& rhs)
{
_superGeometry = rhs._superGeometry;
_statisticsUpdateNeeded = true;
_overlap = rhs._overlap;
_nMaterials = rhs._nMaterials;
_minOverMaterial = rhs._minOverMaterial;
_maxOverMaterial = rhs._maxOverMaterial;
return *this;
}
template
bool& SuperGeometryStatistics3D::getStatisticsStatus()
{
return _statisticsUpdateNeeded;
}
template
bool const & SuperGeometryStatistics3D::getStatisticsStatus() const
{
return _statisticsUpdateNeeded;
}
template
void SuperGeometryStatistics3D::update(bool verbose)
{
#ifdef PARALLEL_MODE_MPI
// This needs to be done with integers due to undesired behaviour with bool and LOR implicated by the MPI standard
int updateReallyNeededGlobal = 0;
if (_statisticsUpdateNeeded) {
updateReallyNeededGlobal = 1;
}
singleton::mpi().reduceAndBcast(updateReallyNeededGlobal, MPI_SUM);
if (updateReallyNeededGlobal>0) {
_statisticsUpdateNeeded = true;
}
//singleton::mpi().reduceAndBcast(_statisticsUpdateNeeded, MPI_LOR);
#endif
// check if update is really needed
if (_statisticsUpdateNeeded ) {
int updateReallyNeeded = 0;
for (int iCloc=0; iCloc<_superGeometry->getLoadBalancer().size(); iCloc++) {
if (_superGeometry->getBlockGeometry(iCloc).getStatistics().getStatisticsStatus() ) {
_superGeometry->getBlockGeometry(iCloc).getStatistics().update(false);
updateReallyNeeded++;
}
}
#ifdef PARALLEL_MODE_MPI
singleton::mpi().reduceAndBcast(updateReallyNeeded, MPI_SUM);
#endif
if (updateReallyNeeded==0) {
_statisticsUpdateNeeded = false;
// clout << "almost updated" << std::endl;
return;
}
// get total number of different materials right
_material2n = std::map();
_nMaterials = int();
for (int iCloc=0; iCloc<_superGeometry->getLoadBalancer().size(); iCloc++) {
if (_superGeometry->getBlockGeometry(iCloc).getStatistics(false).getNmaterials() > 0) {
_nMaterials = _superGeometry->getBlockGeometry(iCloc).getStatistics(false).getNmaterials();
}
}
#ifdef PARALLEL_MODE_MPI
singleton::mpi().reduceAndBcast(_nMaterials, MPI_SUM);
#endif
// store the number and min., max. possition for each rank
for (int iCloc=0; iCloc<_superGeometry->getLoadBalancer().size(); iCloc++) {
std::map material2n = _superGeometry->getBlockGeometry(iCloc).getStatistics(false).getMaterial2n();
std::map::iterator iter;
for (iter = material2n.begin(); iter != material2n.end(); iter++) {
if (iter->second!=0) {
std::vector minPhysR = _superGeometry->getBlockGeometry(iCloc).getStatistics(false).getMinPhysR(iter->first);
std::vector maxPhysR = _superGeometry->getBlockGeometry(iCloc).getStatistics(false).getMaxPhysR(iter->first);
if (_material2n.count(iter->first) == 0) {
_material2n[iter->first] = iter->second;
_material2min[iter->first] = minPhysR;
_material2max[iter->first] = maxPhysR;
//std::cout << iter->first<<":"<<_material2n[iter->first]<first] += iter->second;
for (int iDim=0; iDim<3; iDim++) {
if (_material2min[iter->first][iDim] > minPhysR[iDim]) {
_material2min[iter->first][iDim] = minPhysR[iDim];
}
if (_material2max[iter->first][iDim] < maxPhysR[iDim]) {
_material2max[iter->first][iDim] = maxPhysR[iDim];
}
}
//std::cout << iter->first<<":"<<_material2n[iter->first]<::iterator iMaterial;
for (iMaterial = _material2n.begin(); iMaterial != _material2n.end(); iMaterial++) {
materials[counter] = iMaterial->first;
materialCount[counter] = iMaterial->second;
for (int iDim=0; iDim<3; iDim++) {
materialMinR[3*counter + iDim] = _material2min[iMaterial->first][iDim];
materialMaxR[3*counter + iDim] = _material2max[iMaterial->first][iDim];
}
counter++;
}
for (int iRank=1; iRank minPhysR(3,T());
std::vector maxPhysR(3,T());
for (int iDim=0; iDim<3; iDim++) {
minPhysR[iDim] = materialMinRinBuf[3*iM + iDim];
maxPhysR[iDim] = materialMaxRinBuf[3*iM + iDim];
}
if (_material2n.count(materialsInBuf[iM]) == 0) {
_material2n[materialsInBuf[iM]] = materialCountInBuf[iM];
_material2min[materialsInBuf[iM]] = minPhysR;
_material2max[materialsInBuf[iM]] = maxPhysR;
} else {
_material2n[materialsInBuf[iM]] += materialCountInBuf[iM];
for (int iDim=0; iDim<3; iDim++) {
if (_material2min[materialsInBuf[iM]][iDim] > minPhysR[iDim]) {
_material2min[materialsInBuf[iM]][iDim] = minPhysR[iDim];
}
if (_material2max[materialsInBuf[iM]][iDim] < maxPhysR[iDim]) {
_material2max[materialsInBuf[iM]][iDim] = maxPhysR[iDim];
}
}
}
}
}
}
#endif
// update _minOverMaterial and _maxOverMaterial
typename std::map< int, std::vector >::const_iterator iter;
// find componentwise minimal extension over all material numbers
for ( int iDim = 0; iDim < 3; iDim++ ) {
// minimum
for ( iter = _material2min.begin(); iter != _material2min.end(); iter++ ) {
if ( iter->first != 0 ) { // only relevant material number are considered
if ( _minOverMaterial[iDim] > iter->second[iDim] ) {
_minOverMaterial[iDim] = iter->second[iDim];
}
}
}
// maximum
for ( iter = _material2max.begin(); iter != _material2max.end(); iter++ ) {
if ( iter->first != 0 ) { // only relevant material number are considered
if ( _maxOverMaterial[iDim] < iter->second[iDim] ) {
_maxOverMaterial[iDim] = iter->second[iDim];
}
}
}
}
//clout.setMultiOutput(true);
// print();
//clout.setMultiOutput(false);
if (verbose) {
clout << "updated" << std::endl;
}
_statisticsUpdateNeeded = false;
}
}
template
int SuperGeometryStatistics3D::getNmaterials()
{
update();
return _nMaterials;
}
template
int SuperGeometryStatistics3D::getNvoxel(int material)
{
update(true);
return _material2n[material];
}
template
int SuperGeometryStatistics3D::getNvoxel()
{
update();
int total = 0;
std::map::iterator iter;
for (iter = _material2n.begin(); iter != _material2n.end(); iter++) {
if (iter->first!=0) {
total+=iter->second;
}
}
return total;
}
template
std::vector SuperGeometryStatistics3D::getMinPhysR(int material)
{
update();
return _material2min[material];
}
template
std::vector SuperGeometryStatistics3D::getMinPhysR()
{
update();
return _minOverMaterial;
}
template
std::vector SuperGeometryStatistics3D::getMaxPhysR(int material)
{
update();
return _material2max[material];
}
template
std::vector SuperGeometryStatistics3D::getMaxPhysR()
{
update();
return _maxOverMaterial;
}
template
std::vector SuperGeometryStatistics3D::getPhysExtend(int material)
{
update();
std::vector extend;
for (int iDim = 0; iDim < 3; iDim++) {
extend.push_back(_material2max[material][iDim] - _material2min[material][iDim]);
}
return extend;
}
template
std::vector SuperGeometryStatistics3D::getPhysRadius(int material)
{
update();
std::vector radius;
for (int iDim=0; iDim<3; iDim++) {
radius.push_back((getMaxPhysR(material)[iDim] - getMinPhysR(material)[iDim])/2.);
}
return radius;
}
template
std::vector SuperGeometryStatistics3D::getCenterPhysR(int material)
{
update();
std::vector center;
for (int iDim=0; iDim<3; iDim++) {
center.push_back(getMinPhysR(material)[iDim] + getPhysRadius(material)[iDim]);
}
return center;
}
template
std::vector SuperGeometryStatistics3D::getType(int iC, int iX, int iY, int iZ)
{
update();
int iCloc=_superGeometry->getLoadBalancer().loc(iC);
std::vector discreteNormal = _superGeometry->getExtendedBlockGeometry(iCloc).getStatistics(false).getType(iX+_overlap, iY+_overlap, iZ+_overlap);
return discreteNormal;
}
template
std::vector SuperGeometryStatistics3D::computeNormal(int material)
{
update();
std::vector normal (3,int());
for (int iCloc=0; iCloc<_superGeometry->getLoadBalancer().size(); iCloc++) {
for (int iDim=0; iDim<3; iDim++) {
if (_superGeometry->getBlockGeometry(iCloc).getStatistics(false).getNvoxel(material)!=0) {
normal[iDim] += _superGeometry->getBlockGeometry(iCloc).getStatistics(false).computeNormal(material)[iDim]*_superGeometry->getBlockGeometry(iCloc).getStatistics(false).getNvoxel(material);
}
}
}
#ifdef PARALLEL_MODE_MPI
for (int iDim=0; iDim<3; iDim++) {
singleton::mpi().reduceAndBcast((normal[iDim]), MPI_SUM);
}
#endif
if (getNvoxel(material) == 0) {
std::cerr << "Unkown material number" << std::endl;
std::exit(-1);
}
for (int iDim=0; iDim<3; iDim++) {
normal[iDim] /= getNvoxel(material);
}
T norm = sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
if (norm>0.) {
normal[0]/=norm;
normal[1]/=norm;
normal[2]/=norm;
}
return normal;
}
template
std::vector SuperGeometryStatistics3D::computeDiscreteNormal(int material, T maxNorm)
{
update();
std::vector normal = computeNormal(material);
std::vector discreteNormal(3,int(0));
T smallestAngle = T(0);
for (int iX = -1; iX<=1; iX++) {
for (int iY = -1; iY<=1; iY++) {
for (int iZ = -1; iZ<=1; iZ++) {
T norm = sqrt(iX*iX+iY*iY+iZ*iZ);
if (norm>0.&& norm=smallestAngle) {
smallestAngle=angle;
discreteNormal[0] = iX;
discreteNormal[1] = iY;
discreteNormal[2] = iZ;
}
}
}
}
}
return discreteNormal;
}
template
T SuperGeometryStatistics3D::computeMaxPhysDistance( int material )
{
Vector vec(getMaxPhysR(material)[0] -getMinPhysR(material)[0], getMaxPhysR(material)[1] -getMinPhysR(material)[1], getMaxPhysR(material)[2] -getMinPhysR(material)[2]);
return vec.norm();
}
template
T SuperGeometryStatistics3D::computeMaxPhysDistance()
{
Vector vec(getMaxPhysR()[0] -getMinPhysR()[0], getMaxPhysR()[1] -getMinPhysR()[1], getMaxPhysR()[2] -getMinPhysR()[2]);
return vec.norm();
}
template
void SuperGeometryStatistics3D::print()
{
update();
std::map::iterator iter;
for (iter = _material2n.begin(); iter != _material2n.end(); iter++) {
clout << "materialNumber=" << iter->first
<< "; count=" << iter->second
<< "; minPhysR=(" << _material2min[iter->first][0] <<","<< _material2min[iter->first][1] <<","<< _material2min[iter->first][2] <<")"
<< "; maxPhysR=(" << _material2max[iter->first][0] <<","<< _material2max[iter->first][1] <<","<< _material2max[iter->first][2] <<")"
<< std::endl;
}
}
} // namespace olb