Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/filter_point.cc
1537c5d55d4dmaster

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 "lm_common.hh"
#include "filter_point.hh"
#include <list>
#include "lib_md.hh"
#include "lib_dd.hh"
#include "lib_continuum.hh"
#include "ref_point_data.hh"
#include <fstream>
#include <sstream>
#include "communicator.hh"
#include "compute_compatibility.hh"
#include "compute_interface.hh"
#include "dumper_text.hh"
#include "dumper.hh"
/* -------------------------------------------------------------------------- */
__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 {
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){
Quantity<Length,3> xmax = cube.getXmax();
Quantity<Length,3> xmin = cube.getXmin();
for (UInt k = 0; k < Dim; ++k) {
xmax[k] += 2.*tolerance;
xmin[k] -= 2.*tolerance;
}
cube.setDimensions(xmin,xmax);
UInt _n = positions.size()/Dim;
for (UInt i = 0; i < _n; ++i) {
const Real * 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]);
}
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
UInt FilterPoint<Cont>::getNbSharedPoints() {
UInt counter = 0;
typename Cont::ContainerSubset::iterator points = this->getIterator(this->dt);
typename Cont::ContainerSubset::Ref point;
for (point = points.getFirst(); !points.end(); point = points.getNext()) {
if ((point.getDOFType() & dt_local_master) == false) {
counter++;
}
}
return counter;
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void FilterPoint<Cont>::build(Cont & cont) {
std::vector<Real> positions_copy;
static const UInt Dim = Cont::Dim;
Cube cube(Dim);
typename Cont::iterator points = cont.getIterator(this->dt);
typename Cont::Ref point;
for (point = points.getFirst(); !points.end(); point = points.getNext()) {
Real pos[Dim];
point.getPositions0(pos);
cube.extendBoundingBox(pos);
}
setPositionsCopy<Dim>(this->positions, cube,this->tolerance,positions_copy);
DUMP(this->getID() << ": positions_copy "
<< positions_copy.size() << " "
<< positions.size(),DBG_INFO);
this->empty();
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);
for (point = points.getFirst(); !points.end(); point = points.getNext()) {
if (findPointInSet(point, indices, positions_copy, this->tolerance)) {
this->add(point);
}
}
//remove dummy point
indices.pop_back();
for (UInt i = 0; i < Dim ; ++i) {
positions_copy.pop_back();
}
nb_points = this->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();
comm.allReduce(&nb_points,1,group,"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(),Dim);
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 < this->getArray().size(); ++i) {
found_positions.push_back(this->getArray()[i].position0(0));
found_positions.push_back(this->getArray()[i].position0(1));
found_positions.push_back(this->getArray()[i].position0(2));
}
found_positions.setCommGroup(group);
DumperText<ComputeInterface> 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 " << this->size() << " points");
}
dimension = Dim;
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
FilterPoint<Cont>::FilterPoint(const FilterID name, FilterInterface & f)
:Filter<Cont>(name,f),path(""), tolerance(1e-5), dimension(0){
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
FilterPoint<Cont>::FilterPoint(const FilterID name, DomainInterface & d)
:Filter<Cont>(name,d),path(""), dimension(0){
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
FilterPoint<Cont>::FilterPoint(const FilterID name, UInt dim):
Filter<Cont>(name,dim), dimension(0){
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
FilterPoint<Cont>::~FilterPoint(){
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void FilterPoint<Cont>::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);
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
UInt FilterPoint<Cont>::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 */
template <typename Cont>
void FilterPoint<Cont>::declareParams(){
Filter<Cont>::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_FILTER(FilterPoint,LIST_ATOM_MODEL);
DECLARE_FILTER(FilterPoint,LIST_CONTINUUM_MODEL);
DECLARE_FILTER(FilterPoint,LIST_DD_MODEL);
DECLARE_FILTER_REFPOINT(FilterPoint);
DECLARE_FILTER_GENERIC_MESH(FilterPoint);
__END_LIBMULTISCALE__
c4science · Help