Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/filter/filter_point.cc
d8e9d1cc530clammps_2020

Last Change
Raw File

filter_point.cc
View Options

/**
* @file filter_point.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Tue Jul 29 17:03:26 2014
*
* @brief This allows to filter a set of DOFs by matching a list of coordinates
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "filter_point.hh"
#include "compute_interface.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "lm_communicator.hh"
#include "ref_point_data.hh"
#include <fstream>
#include <list>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename Ref>
inline bool findPointInSet(Ref &point, std::list<UInt> &indices,
const std::vector<Real> &positions,
const Real &tolerance) {
// static const UInt Dim = Ref::Dim;
std::list<UInt>::iterator index = indices.begin();
std::list<UInt>::iterator end_index = indices.end();
Real distance = 0;
do {
LM_TOIMPLEMENT;
// distance = point.computeInitialDistance(&positions.at(Dim * (*index)));
++index;
} while ((distance > tolerance) && (index != end_index));
if (index == end_index) {
return false;
} else {
indices.erase(--index);
return true;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void setPositionsCopy(const std::vector<Real> &positions, Cube cube,
Real tolerance, std::vector<Real> &positions_copy) {
auto &&xmax = cube.getXmax();
auto &&xmin = cube.getXmin();
for (UInt k = 0; k < Dim; ++k) {
xmax[k] += 2. * tolerance;
xmin[k] -= 2. * tolerance;
}
cube.setDimensions(xmin.cast<PhysicalScalar<Length>>(),
xmax.cast<PhysicalScalar<Length>>());
UInt _n = positions.size() / Dim;
for (UInt i = 0; i < _n; ++i) {
Vector<Dim> tmp(&positions[i * Dim]);
bool test = static_cast<Geometry &>(cube).template contains<Dim>(tmp);
if (test) {
for (UInt k = 0; k < Dim; ++k)
positions_copy.push_back(positions[i * Dim + k]);
}
}
}
/* -------------------------------------------------------------------------- */
UInt FilterPoint::getNbSharedPoints() {
LM_TOIMPLEMENT;
// UInt counter = 0;
// auto it = this->begin(this->dt);
// auto end = this->end(this->dt);
// for (; it != end; ++it) {
// auto &&point = *it;
// if ((point.getDOFType() & dt_local_master) == false) {
// counter++;
// }
// }
// return counter;
}
/* -------------------------------------------------------------------------- */
template <typename Cont> void FilterPoint::build([[gnu::unused]] Cont &cont) {
LM_TOIMPLEMENT;
// auto &output =
// allocAndGetCastedOutput<typename Cont::ContainerSubset>(this->getID(),
// this->getID());
// std::vector<Real> positions_copy;
// static const UInt Dim = Cont::Dim;
// Cube cube(Dim);
// auto it = cont.begin(this->dt);
// auto end = cont.end(this->dt);
// for (; it != end; ++it) {
// auto &&point = *it;
// auto &&pos = point.position0();
// cube.extendBoundingBox(pos);
// }
// setPositionsCopy<Dim>(this->positions, cube, this->tolerance,
// positions_copy); DUMP(this->getID() << ": positions_copy " <<
// positions_copy.size() << " "
// << positions.size(),
// DBG_INFO);
// output.clear();
// std::list<UInt> indices;
// for (UInt i = 0; i < Dim; ++i) {
// // This is a dummy point I add for the findPointInSet algo, cleaned later
// positions_copy.push_back(lm_real_max);
// }
// UInt nb_points = positions_copy.size() / Dim;
// for (UInt i = 0; i < nb_points; ++i) {
// indices.push_back(i);
// }
// DUMP("number of indices = " << indices.size(), DBG_INFO);
// it = cont.begin(this->dt);
// for (; it != end; ++it) {
// auto &&point = *it;
// if (findPointInSet(point, indices, positions_copy, this->tolerance)) {
// output.push_back(point);
// }
// }
// // remove dummy point
// indices.pop_back();
// for (UInt i = 0; i < Dim; ++i) {
// positions_copy.pop_back();
// }
// nb_points = output.size();
// if (this->dt != dt_local_master) {
// UInt shared = this->getNbSharedPoints();
// DUMP("locally found " << nb_points << " total of which " << shared
// << " are shared.",
// DBG_DETAIL);
// nb_points -= shared;
// }
// // Communicator &comm = Communicator::getCommunicator();
// CommGroup group = cont.getCommGroup();
// group.allReduce(&nb_points, 1, "DeviationOperator reduce res", OP_SUM);
// if (nb_points == 0) {
// LM_FATAL(this->getID() << ": Could not match any point out of "
// << positions.size() / Dim);
// }
// if (nb_points != positions.size() / Dim) {
// std::stringstream sstr;
// sstr << this->getID() << "-found_pos-" << lm_my_proc_id;
// ComputeCompatibility found_positions(sstr.str());
// // found_positions.name_computed.resize(Dim);
// // for (UInt k = 0; k < Dim; ++k) {
// // std::stringstream sstr;
// // sstr << "pos" << k;
// // found_positions.name_computed[k] = sstr.str();
// // }
// for (UInt i = 0; i < output.size(); ++i) {
// auto &&point = output[i];
// auto &&pos = point.position0();
// found_positions.push_back(pos[0]);
// found_positions.push_back(pos[1]);
// found_positions.push_back(pos[2]);
// }
// found_positions.setCommGroup(group);
// LM_TOIMPLEMENT;
// // DumperText error_positions(sstr.str(), found_positions);
// // error_positions.setParam("SEPARATOR", std::string(","));
// // error_positions.setParam("RANK", comm.groupRank(lm_my_proc_id,
// group));
// // error_positions.manualDump();
// LM_FATAL(this->getID() << ": Could match " << nb_points - 1 << " out of "
// << positions.size() / Dim << " positions"
// << " locally I found " << output.size()
// << " points");
// }
// dimension = Dim;
}
/* -------------------------------------------------------------------------- */
FilterPoint::FilterPoint(const LMID &name)
: LMObject(name), path(""), tolerance(1e-5), dimension(0) {}
/* -------------------------------------------------------------------------- */
FilterPoint::~FilterPoint() {}
/* -------------------------------------------------------------------------- */
void FilterPoint::init() {
// yeah, could be done more elegantly
std::ifstream infile(this->path.c_str());
if (!infile.good()) {
LM_FATAL("The atom file " << path.c_str()
<< " could NOT be opened, rdstate = " << std::hex
<< infile.rdstate() << std::endl
<< " badbit = " << infile.bad() << std::endl
<< " failbit = " << infile.fail() << std::endl
<< " eofbit = " << infile.eof() << std::endl
<< " goodbit = " << infile.good() << std::endl);
}
char tmp;
while (infile.good()) {
tmp = infile.peek();
switch (tmp) {
case '#':
infile.ignore(1024, '\n');
break;
default:
while (infile.good()) {
this->positions.push_back(0);
infile >> this->positions.back();
}
break;
}
}
if (this->positions.size() == 0)
LM_FATAL("no valid positions found in file " << this->path.c_str());
this->positions.pop_back();
DUMP("found " << this->positions.size() / 3 << " points", DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
UInt FilterPoint::getNbPointsInFile() {
return this->positions.size() / dimension;
}
/* -------------------------------------------------------------------------- */
/* LMDESC POINT
This filter allows to extract a set of points and/or
elements contained in a provided point
*/
/* LMHERITANCE filter */
/* LMEXAMPLE FILTER fpoint POINT INPUT md PATH ./my-point-list.txt TOL 1e-5 */
void FilterPoint::declareParams() {
FilterInterface::declareParams();
/* LMKEYWORD PATH
Specify the file containing the points to filter
*/
this->parseKeyword("PATH", path);
/* LMKEYWORD TOL
Specify the file containing the points to filter
*/
this->parseKeyword("TOL", tolerance, 1e-5);
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(FilterPoint)
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help