Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/filter/compute_ekin_mesh.cc
b05d5f2aba78features/coupler_compon…

Last Change
Raw File

compute_ekin_mesh.cc
View Options

/**
* @file compute_ekin_mesh.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Tue Dec 03 22:43:33 2013
*
* @brief This computes the kinetic energy of a set of DOFs
*
* @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 "lib_md.hh"
#include "lib_dd.hh"
#include "lib_continuum.hh"
#include "compute_ekin.hh"
#include "communicator.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename Cont>
void ComputeEKin<Cont>::build(Cont & cont){
static const UInt Dim = Cont::Dim;
if(based_on_nodes_flag){
buildBasedOnNodes<Dim>(cont);
return;
}
Real computed_ekin = 0.0;
typename Cont::ContainerElems::iterator it = cont.getContainerElems().getIterator();
typedef typename Cont::RefElem RefPoint;
for (RefPoint el = it.getFirst();!it.end();el = it.getNext()) {
computed_ekin += el.getEKin();
}
Communicator & comm = Communicator::getCommunicator();
CommGroup group = cont.getCommGroup();
comm.allReduce(&computed_ekin,1,group," kinetic energy reduce", OP_SUM);
this->empty();
this->setDim(1);
this->add(computed_ekin);
this->name_computed.clear();
this->name_computed.push_back("EKin");
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
template <UInt Dim>
void ComputeEKin<Cont>::buildBasedOnNodes(Cont & cont){
Real computed_ekin[3] = {0,0,0};
UInt cpt = cont.nbElem();
typename Cont::iterator it = cont.getIterator();
typedef typename Cont::Ref RefPoint;
for (RefPoint at = it.getFirst();!it.end();at = it.getNext()) {
Real vel[Dim];
at.getVelocities(vel);
for (UInt i = 0; i < Dim; ++i) computed_ekin[i] += 0.5*at.mass()*vel[i]*vel[i]*code_unit_system.mvv2e;
}
Communicator & comm = Communicator::getCommunicator();
CommGroup group = cont.getCommGroup();
comm.allReduce(computed_ekin,Dim,group,"kinetic energy reduce", OP_SUM);
comm.allReduce(&cpt,1,group,"kinetic energy counter reduce", OP_SUM);
Real ekin = 0.0;
for (UInt i = 0; i < Dim; ++i){
DUMP("ekin[" << i << "] = " << computed_ekin[i],DBG_INFO);
ekin += computed_ekin[i];
}
DUMP("nb_points = " << cpt,DBG_INFO);
this->empty();
Real nb_points = static_cast<Real>(cpt);
this->add(nb_points);
this->add(ekin);
for (UInt i = 0; i < Dim; ++i) this->add(computed_ekin[i]);
this->name_computed.clear();
this->name_computed.push_back("nbpoints");
this->name_computed.push_back("EKin");
this->name_computed.push_back("EKin_x");
if (Dim > 1) this->name_computed.push_back("EKin_y");
if (Dim == 3) this->name_computed.push_back("EKin_z");
this->setDim(Dim + 2);
DUMP("Kinetic energy based on nodes could be higher than the one computed from elements.", DBG_WARNING);
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_REF(ComputeEKin,LIST_CONTINUUM_MODEL);
__END_LIBMULTISCALE__
c4science · Help