Homec4science
Diffusion LibMultiScale (features/cadd) Browse

libmultiscale/src/coupling/point_association_inline_impl.hh
0666acaf17c9features/cadd

Last Change
Raw File

point_association_inline_impl.hh
View Options

/**
* @file point_association_inline_impl.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
*
* @date Mon Jun 16 17:13:26 2014
*
* @brief computes the association between point degrees of freedom of
* different domains
*
* @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 "point_association.hh"
#include <string>
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA, typename ContB>
inline void PointAssociation::mixField(ContA &pointsAlist, ContB &pointsBlist) {
if (this->in_group_A && !this->in_group_B)
LM_FATAL("not a parallel routine yet");
if (pointsAlist->size() != pointsBlist->size())
LM_FATAL("inconsistency here");
constexpr UInt Dim = ContA::Dim;
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
// make the average in the field_buffer thing
UInt index = 0;
UInt index_weight = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] = weightsA[index_weight] * _field[dir];
}
++index_weight;
}
index = 0;
index_weight = 0;
for (auto &&point : pointsBlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] += weightsB[index_weight] * _field[dir];
}
++index_weight;
}
// apply the result to the fields
index = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
_field[dir] = this->field_buffer[index];
}
}
index = 0;
for (auto &&point : pointsBlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
_field[dir] = this->field_buffer[index];
}
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA>
inline void PointAssociation::prepareSynchronisation(ContA &pointsAlist) {
if (ftype == _displacement) {
this->prepareSynchronisationDisplacement(pointsAlist);
return;
}
constexpr UInt Dim = ContA::Dim;
if (this->in_group_A) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
UInt index = 0;
for (auto &&point : pointsAlist) {
auto &&_field = point.template field<ftype>();
for (UInt dir = 0; dir < Dim; ++dir, ++index) {
this->field_buffer[index] = _field[dir];
}
}
} else if (this->in_group_B) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
}
}
/* -------------------------------------------------------------------------- */
// specialize this for ComputeTrueDisplacement
template <typename ContA>
inline void
PointAssociation::prepareSynchronisationDisplacement(ContA &pointsAlist) {
constexpr UInt Dim = ContA::Dim;
if (this->in_group_A) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
DUMP(this->getID() << ": bounding box of " << pointsAlist.getID() << " is "
<< pointsAlist->getBoundingBox(),
DBG_INFO);
this->true_disp.buildManual(this->pointsAlist);
UInt index = 0;
for (auto &&disp : this->true_disp.getArray()) {
this->field_buffer[index] = disp;
DUMP(this->getID() << ": prepareSyncDisp.template field<" << _displacement
<< ">(" << index % Dim << ") = " << disp,
DBG_DETAIL);
++index;
}
} else if (this->in_group_B) {
this->field_buffer.resize(Dim * pointsAlist->size(), 0);
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype> inline void PointAssociation::synchronise() {
std::stringstream sstr;
sstr << "Sync " << ftype << " for object " << this->getID();
this->field_buffer.resize(this->local_points, spatial_dimension);
this->distributeVectorA2B(sstr.str(), this->field_buffer);
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype, typename ContA, typename ContB>
inline void PointAssociation::applyBuffer(ContA &pointsAlist,
ContB &pointsBlist) {
if (ftype == _position) {
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0)
this->template applyBufferPosition<false, false, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 0)
this->template applyBufferPosition<true, false, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 0)
this->template applyBufferPosition<true, true, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 1)
this->template applyBufferPosition<true, true, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 1)
this->template applyBufferPosition<true, false, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 0)
this->template applyBufferPosition<false, true, false>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 1)
this->template applyBufferPosition<false, true, true>(pointsAlist,
pointsBlist);
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 1)
this->template applyBufferPosition<false, false, true>(pointsAlist,
pointsBlist);
return;
}
// constexpr UInt Dim = ContainerA::Dim;
for (UInt point_idx = 0; point_idx < pointsBlist->size(); ++point_idx) {
auto point = pointsBlist->get(point_idx);
DUMP(this->getID() << ": apply buffer for point " << point, DBG_DETAIL);
auto &&field_vals = this->field_buffer(point_idx);
auto &&_field = point.template field<ftype>();
_field = field_vals;
// for (UInt dir = 0; dir < Dim; ++dir) {
// _field[dir] = field_vals[dir];
// DUMP(this->getID() << ": point.template field<" << ftype << ">(" << dir
// << ") = " << field_vals[dir],
// DBG_DETAIL);
// }
}
}
/* -------------------------------------------------------------------------- */
template <typename Vec1, typename Vec2, typename Vec3, UInt Dim, bool pbc_x,
bool pbc_y, bool pbc_z>
inline void fix_pbc(Vec1 &pos, Vec2 &xmin, Vec3 &d_size) {
for (UInt i = 0; i < Dim; ++i) {
if ((pbc_x && i == 0) || (pbc_y && i == 1) || (pbc_z && i == 2)) {
int offset = floor((pos[i] - xmin[i]) / d_size[i]);
pos[i] -= offset * d_size[i];
}
}
}
/* -------------------------------------------------------------------------- */
template <bool pbc_x, bool pbc_y, bool pbc_z, typename ContA, typename ContB>
inline void PointAssociation::applyBufferPosition(ContA &pointsAlist,
ContB &pointsBlist) {
DUMP(this->getID() << ": bounding box of " << this->pointsBlist.getID()
<< " is " << pointsBlist->getBoundingBox(),
DBG_INFO);
constexpr UInt Dim = ContA::Dim;
Cube &bbox = pointsBlist->getBoundingBox();
auto &&d_size = bbox.getSize();
auto &&xmin = bbox.getXmin<Dim>();
for (UInt point_idx = 0; point_idx < pointsBlist->size(); ++point_idx) {
auto point = pointsBlist->get(point_idx);
Vector<Dim> field_vals = this->field_buffer.row(point_idx);
DUMP(this->getID() << ": apply buffer position for point " << point,
DBG_DETAIL);
fix_pbc<Vector<Dim>, Vector<Dim>, Vector<3, Real>, Dim, pbc_x, pbc_y,
pbc_z>(field_vals, xmin, d_size);
auto &&_field = point.template field<_position>();
for (UInt dir = 0; dir < Dim; ++dir) {
_field[dir] = field_vals[dir];
DUMP(this->getID() << ": point.template field<" << _position << ">("
<< dir << ") = " << field_vals[dir],
DBG_DETAIL);
}
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype> inline void PointAssociation::enforceField() {
if (this->in_group_A) {
STARTTIMER("Point Assoc. Synch.");
this->prepareSynchronisation<ftype>();
STOPTIMER("Point Assoc. Synch.");
}
this->synchronise<ftype>();
if (this->in_group_B) {
STARTTIMER("Apply Assoc. Synch.");
this->applyBuffer<ftype>();
STOPTIMER("Apply Assoc. Synch.");
}
}
/* -------------------------------------------------------------------------- */
template <FieldType ftype1, FieldType ftype2, typename ContA, typename ContB>
inline void PointAssociation::zeroFields(ContA &pointsAlist,
ContB &pointsBlist) {
constexpr UInt Dim = ContA::Dim;
for (UInt point_idx = 0; point_idx < pointsBlist->size(); ++point_idx) {
auto point = pointsBlist->get(point_idx);
auto &&f1 = point.template field<ftype1>();
auto &&f2 = point.template field<ftype1>();
for (UInt dir = 0; dir < Dim; ++dir) {
f1[dir] = 0.;
f2[dir] = 0.;
}
}
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
c4science · Help