Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/filter/compute_extract.cc
1d1b2315bf45master

Last Change
Raw File

compute_extract.cc
View Options

/**
* @file compute_extract.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This compute extract data from point-like references
*
* @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 "compute_extract.hh"
#include "container.hh"
#include "factory_multiscale.hh"
#include "field.hh"
#include "filter_geometry.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"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
ComputeExtract::ComputeExtract(const std::string &name)
: LMObject(name), point_tags("tags:" + name) {
component = -1;
// register for computes
this->createInput("input");
this->createOutput("field");
this->createOutput("point_tags") = this->point_tags;
this->default_output = "field";
}
/* -------------------------------------------------------------------------- */
ComputeExtract::~ComputeExtract() {}
/* -------------------------------------------------------------------------- */
/* LMDESC EXTRACT
This computes generate an array of real containing the requested field
*/
/* LMEXAMPLE
COMPUTE disp EXTRACT INPUT md FIELD displacement COMPONENT 0
*/
inline void ComputeExtract::declareParams() {
FilterInterface::declareParams();
/* LMKEYWORD FIELD
Specify the field to extract
*/
this->parseKeyword("FIELD", field);
/* LMKEYWORD COMPONENT
Specify the coordinate to extract. Default is -1 and means all 3
coordinates are extracted
*/
this->parseKeyword("COMPONENT", component, -1);
}
/* -------------------------------------------------------------------------- */
template <typename Cont> enable_if_md<Cont> ComputeExtract::build(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
this->clear();
std::stringstream sstr;
sstr << this->getID() << ":" << field;
this->getOutput("field") = ContainerArray<Real>(sstr.str());
switch (field) {
case _displacement:
buildContainerField<Dim, _displacement>(cont);
break;
case _position0:
buildContainerField<Dim, _position0>(cont);
break;
case _position:
buildContainerField<Dim, _position>(cont);
break;
case _velocity:
buildContainerField<Dim, _velocity>(cont);
break;
case _force:
buildContainerField<Dim, _force>(cont);
break;
case _grouprank:
this->getOutputAsArray().push_back(cont.getCommGroup().getMyRank());
break;
case _mass:
buildContainerField<1, _mass>(cont);
break;
// case _epot:
// buildContainerField<1, _epot>(cont);
// break;
// case _tag:
// buildContainerField<1, _tag>(cont);
// break;
default:
LM_FATAL("unknown field type " << field);
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, FieldType ftype, typename Cont>
void ComputeExtract::buildContainerField(Cont &cont) {
this->point_tags.clear();
auto &array = this->getOutputAsArray("field");
array.clear();
if (component >= 0)
array.resize(0, 1);
else
array.resize(0, Dim);
for (auto &&at : cont) {
this->point_tags.push_back(at.tag());
decltype(auto) f = at.template field<ftype>();
DUMP("passing atom " << at.tag(), DBG_DETAIL);
if (component >= 0){
array.push_back(f[component]);
}
else {
#if __GNUC__ > 8
array.push_back(Vector<Dim>(f));
#else
array.push_back(f.operator Vector<Dim>());
#endif
}
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
enable_if_mesh<Cont> ComputeExtract::build(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
std::stringstream sstr;
sstr << this->getID() << ":" << field;
this->getOutput("field") = ContainerArray<Real>(sstr.str());
Cont::Ref::fields::foreach ([&](auto &&a) {
constexpr auto field_type = std::decay_t<decltype(a)>::value;
if (field_type == field) {
buildContainerNodalField<Dim, field_type>(cont);
}
});
Cont::RefElem::fields::foreach ([&](auto &&a) {
constexpr auto field_type = std::decay_t<decltype(a)>::value;
if (field_type == field) {
auto _field = make_field<field_type>(cont.getContainerElems());
buildContainerElementalField(_field);
}
});
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, FieldType ftype, typename Cont>
void ComputeExtract::buildContainerNodalField(Cont &cont) {
auto &array = this->getOutputAsArray();
array.clear();
if (component >= 0)
array.resize(0, 1);
else
array.resize(0, Dim);
for (auto &&at : cont) {
decltype(auto) f = at.template field<ftype>();
if (component >= 0) {
array.push_back(f[component]);
} else {
Vector<Dim, typename decltype(f)::value_type> _f = f;
Vector<Dim> __f = _f.template cast<Real>();
array.push_back(__f);
}
}
}
/* -------------------------------------------------------------------------- */
template <typename Field>
void ComputeExtract::buildContainerElementalField(Field &f) {
auto &array = this->getOutputAsArray();
array.clear();
UInt Dim = f.getDim();
UInt rows = array.rows();
if (component >= 0)
array.resize(rows, 1);
else
array.resize(rows, Dim);
for (auto &&el : f) {
if (component >= 0) {
array.push_back(el[component]);
} else {
array.push_back(ArrayView(el));
}
}
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_MAKE_CALL(ComputeExtract)
__END_LIBMULTISCALE__
c4science · Help