Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/filter/filter_geometry.cc
d8e9d1cc530clammps_2020

Last Change
Raw File

filter_geometry.cc
View Options

/**
* @file filter_geometry.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Thu Jul 24 14:22:22 2014
*
* @brief This filter allows to select DOFs based on geometrical criteria
*
* @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_geometry.hh"
#include "compute_interface.hh"
#include "factory_multiscale.hh"
#include "geometry_manager.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "reference_manager_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
/* LMDESC GEOM
This filter allows to extract a set of points and/or
elements contained in a provided geometry
*/
/* LMHERITANCE filter */
/* LMEXAMPLE
FILTER topAtoms GEOM INPUT md GEOMETRY top
*/
/* -------------------------------------------------------------------------- */
void FilterGeometry::init() {
std::vector<ComputeInterface *> computes_ptr;
computes_ptr.resize(compute_list.size());
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
LMID id = compute_list[comp];
auto &my_compute =
FilterManager::getManager().getCastedObject<ComputeInterface>(id);
computes_ptr[comp] = &my_compute;
my_compute.compute();
LM_TOIMPLEMENT;
// ComputeCompatibility *out_compute =
// new ComputeCompatibility("", computes_ptr[comp]->getDim());
// out_computes.push_back(out_compute);
}
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
LMID id = compute_list[comp];
std::stringstream out_compute_name;
out_compute_name << id << ":" << this->getID();
// this->out_computes[comp]->name_computed.push_back(out_compute_name.str());
this->out_computes[comp]->setID(out_compute_name.str());
}
}
/* -------------------------------------------------------------------------- */
void FilterGeometry::declareParams() {
FilterInterface::declareParams();
/* LMKEYWORD GEOMETRY
Specify the geometry on which the filtering must proceed
*/
this->parseKeyword("GEOMETRY", geom);
/* LMKEYWORD CURRENT_POSITION
Specify that the selection should be based on the
current position as opposed to the initial position
*/
this->parseTag("CURRENT_POSITION", current_position_flag, false);
/* LMKEYWORD DONOTRESETBOUNDINGBOX
Flag to use the bounding box of complete domain
*/
this->parseTag("DONOTRESETBOUNDINGBOX", do_not_reset_boundingbox_flag, false);
/* LMKEYWORD CENTROID
Specify that the selection should be based on the
centriod of the element
*/
this->parseTag("CENTROID", centroid_flag, false);
/* LMKEYWORD BASED_ON_NODES
Specify that the selection should be based on the
nodes of the element
*/
this->parseTag("BASED_ON_NODES", based_on_nodes_flag, false);
/* LMKEYWORD SORT
Flag to sort the container elements
*/
this->parseTag("SORT", sort_flag, false);
/* LMKEYWORD DO_NOT_REGISTER
used in point association
*/
this->parseTag("DO_NOT_REGISTER", do_not_register, false);
/* LMKEYWORD ADD_COMPUTE
Add computes to input to filter
*/
this->parseKeyword("ADD_COMPUTE", compute_list, std::vector<LMID>{});
}
/* -------------------------------------------------------------------------- */
template <typename Ref> class ComparatorRef {
public:
bool operator()(const Ref &r1, const Ref &r2) {
auto pos1 = r1.position0();
auto pos2 = r2.position0();
for (UInt i = 0; i < pos1.size(); ++i) {
if (pos1[i] != pos2[i])
return pos1[i] < pos2[i];
}
return r1.position0()[0] < r2.position0()[0];
if (r1.position0()[1] != r2.position0()[1])
return r1.position0()[1] < r2.position0()[1];
return r1.position0()[2] < r2.position0()[2];
}
bool operator()(Ref &r1, Ref &r2) {
auto pos1 = r1.position0();
auto pos2 = r2.position0();
for (UInt i = 0; i < pos1.size(); ++i) {
if (pos1[i] != pos2[i])
return pos1[i] < pos2[i];
}
return r1.position0()[0] < r2.position0()[0];
if (r1.position0()[1] != r2.position0()[1])
return r1.position0()[1] < r2.position0()[1];
return r1.position0()[2] < r2.position0()[2];
}
};
/* -------------------------------------------------------------------------- */
FilterGeometry::~FilterGeometry() {}
/* -------------------------------------------------------------------------- */
FilterGeometry::FilterGeometry(const std::string &name)
: LMObject(name), geom(invalidGeom), current_position_flag(false),
centroid_flag(false), based_on_nodes_flag(false),
do_not_reset_boundingbox_flag(false), do_not_register(false),
sort_flag(false) {
this->createInput("input");
this->createOutput("references");
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
enable_if_point<Cont> FilterGeometry::build(Cont &cont) {
auto &output = allocOutput<typename Cont::ContainerSubset>("references");
constexpr UInt Dim = Cont::Dim;
if (output.size() && !current_position_flag)
return;
if (!do_not_reset_boundingbox_flag) {
output.getBoundingBox().resetDimensions();
}
output.clear();
auto g = GeometryManager::getManager().getGeometry(geom);
std::vector<ComputeInterface *> computes_ptr;
computes_ptr.resize(compute_list.size());
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
this->out_computes[comp]->clear();
LMID id = compute_list[comp];
auto &my_compute =
FilterManager::getManager().getCastedObject<ComputeInterface>(id);
computes_ptr[comp] = &my_compute;
my_compute.compute();
LM_TOIMPLEMENT;
// this->out_computes[comp]->copyReleaseInfo(cont);
// this->out_computes[comp]->copyContainerInfo(cont);
}
int counter = 0;
for (auto &&at : cont) {
Vector<Dim> positions;
positions.setZero();
if (current_position_flag) {
positions = at.position();
} else {
positions = at.position0();
}
DUMP("Testing point at position " << at.position()[0], DBG_ALL);
if (g->contains<Dim>(positions)) {
if (!do_not_reset_boundingbox_flag)
output.getBoundingBox().extendBoundingBox(positions);
DUMP("Accepting atom " << at, DBG_INFO);
output.push_back(at);
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
LM_TOIMPLEMENT;
// this->out_computes[comp]->push_back(computes_ptr[comp]->get(counter));
}
}
counter++;
}
#ifndef LM_OPTIMIZED
UInt nb = output.size();
#endif // LM_OPTIMIZED
DUMP("The filter contains " << nb << " atoms in geometry " << *g, DBG_INFO);
if (!this->do_not_register) {
if (!current_position_flag) {
output.addSubSet(output);
}
}
if (sort_flag) {
// std::vector<typename Cont::Ref> &arr = output.getArray();
// do sort
std::sort(output.begin(), output.end(),
ComparatorRef<typename Cont::Ref>());
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
enable_if_mesh<Cont> FilterGeometry::build(Cont &cont) {
auto &output =
allocOutput<typename Cont::ContainerSubset>("references", this->getID());
constexpr UInt Dim = Cont::Dim;
if (output.size() && !current_position_flag)
return;
if (!do_not_reset_boundingbox_flag)
output.getBoundingBox().resetDimensions();
if (based_on_nodes_flag) {
buildBasedOnNodes<Dim>(cont);
return;
}
auto g = GeometryManager::getGeometry(geom);
output.clear();
auto &contElems = cont.getContainerElems();
auto &contNodes = cont.getContainerNodes();
auto &contElemsFiltered = output.getContainerElems();
for (auto &&el : contElems) {
auto &&connectivity = el.globalIndexes();
UInt nb_nodes = connectivity.size();
bool test = false;
if (centroid_flag) {
Vector<Dim> centroid;
centroid.setZero();
for (UInt i = 0; i < nb_nodes; ++i) {
auto nd = contNodes.get(connectivity[i]);
auto X = nd.position0();
centroid += X;
centroid /= nb_nodes;
test = test | g->contains(centroid);
}
} else {
for (UInt i = 0; i < nb_nodes; ++i) {
auto nd = contNodes.get(connectivity[i]);
test = test | g->contains(nd.position0());
}
}
if (test) {
contElemsFiltered.push_back(el);
}
}
// recompute connectivity in altered connectivity
output.computeAlteredConnectivity(contNodes);
// reshape the geometry
if (!do_not_reset_boundingbox_flag) {
for (auto &&nd : output.getContainerNodes()) {
auto X = nd.position0();
output.getBoundingBox().extendBoundingBox(X);
}
}
}
/* --------------------------------------------------------------------------
*/
template <UInt Dim, typename Cont>
void FilterGeometry::buildBasedOnNodes(Cont &cont) {
auto &output =
allocOutput<typename Cont::ContainerSubset>("references", this->getID());
// using type = typename Cont::ContainerSubset;
// using ptr_type = std::shared_ptr<type>;
// allocAndGetCastedOutput<ptr_type>
// try {
// getCastedOutput<ptr_type>();
// } catch (ArgumentContainer::UnallocatedPointer &e) {
// this->getOutput() =
// make_argument(std::make_shared<type>(this->getID()));
// }
// auto &output = *getCastedOutput<ptr_type>();
output.clear();
auto g = GeometryManager::getManager().getGeometry(geom);
for (auto &&at : cont) {
auto positions0 = at.position0();
DUMP("Testing point at position " << at.position0(), DBG_ALL);
if (g->contains<Dim>(positions0)) {
if (!do_not_reset_boundingbox_flag)
output.getBoundingBox().extendBoundingBox(positions0);
DUMP("Accepting point " << at, DBG_ALL);
output.push_back(at);
}
}
#ifndef LM_OPTIMIZED
UInt nb = output.size();
#endif // LM_OPTIMIZED
DUMP("le filter a trouve " << nb << " atoms in geometry " << *g, DBG_INFO);
}
/* --------------------------------------------------------------------------
*/
DECLARE_COMPUTE_MAKE_CALL(FilterGeometry)
/* --------------------------------------------------------------------------
*/
__END_LIBMULTISCALE__
c4science · Help