Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/coupling/bridging_atomic_continuum.cc
0b9417d37655features/coupler_compon…

Last Change
Raw File

bridging_atomic_continuum.cc
View Options

/**
* @file bridging_atomic_continuum.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief Bridging object between atomistic and finite elements
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
//#define TIMER
/* -------------------------------------------------------------------------- */
#include "bridging_atomic_continuum.hh"
#include "compute_extract.hh"
#include "lib_bridging.hh"
#include "lm_common.hh"
#include "reference_manager.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename DomainA, typename DomainC>
BridgingAtomicContinuum<DomainA, DomainC>::BridgingAtomicContinuum(
const std::string &name, DomainA &A, DomainC &C)
: LMObject(name), Bridging<typename DomainA::ContainerAtoms,
typename DomainC::ContainerPoints>(
name, A.getContainer(), C.getContainer()),
domA(A), domC(C), comm_group_continuum(domC.getCommGroup()),
comm_group_atomic(domA.getCommGroup()) {
if (domA.getContainer().hasRefManager())
domA.getContainer().getRefManager()->addSubSet(this->pointList);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
BridgingAtomicContinuum<DomainA, DomainC>::~BridgingAtomicContinuum() {}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::clearAll() {
Bridging<typename DomainA::ContainerPoints,
typename DomainC::ContainerPoints>::clearAll();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::updateForMigration() {
Bridging<typename DomainA::ContainerPoints,
typename DomainC::ContainerPoints>::updateForMigration();
if (this->check_coherency) {
checkPositionCoherency();
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::init() {
Bridging<typename DomainA::ContainerPoints,
typename DomainC::ContainerPoints>::init();
DUMP(this->getID() << " : attach duo vector", DBG_MESSAGE);
if (this->in_group_A) {
try {
this->contA.getRefManager()->attachObject(this->getDuoVector(),
this->pointList);
DUMP(this->getID() << " : attaching Parent::duo_vector", DBG_INFO);
} catch (LibMultiScaleException &e) {
}
}
MPI_Barrier(MPI_COMM_WORLD);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::averageDOFsOnElements(
ContainerArray<Real> &data_mesh, FieldType field) {
LM_FATAL("to adapt");
if (comm_group_atomic == comm_group_continuum) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
extracted_field.build(*this->pointList);
this->averageOnElements(extracted_field.getArray(), data_mesh);
return;
}
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
if (this->is_in_atomic)
extracted_field.build(*this->pointList);
LM_TOIMPLEMENT;
// extracted_field.resize(Parent::local_points * Dim);
/* do communication */
this->communicateBufferAtom2Continuum(extracted_field.getArray());
if (this->is_in_continuum)
this->averageOnElements(extracted_field.getArray(), data_mesh);
}
/* -------------------------------------------------------------------------- */
// template <typename DomainA, typename DomainC,UInt Dim>
// template <FieldType field>
// void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
// buildAtomicDOFsVector(std::vector<Real> & v_out){
// IteratorAtomsSubset it = this->pointList->begin();
// UInt nmax = this->pointList->size();
// UInt i =0;
// Real * v_out_ptr = &v_out[0];
// for (RefAtom at = it.getFirst(); !it.end() ;
// at = it.getNext() , ++i, v_out_ptr+= Dim) {
// LM_ASSERT(i < nmax,"overflow detected");
// at.template getField<field>(v_out_ptr);
// }
// }
// /* --------------------------------------------------------------------------
// */
// template <typename DomainA, typename DomainC,UInt Dim>
// void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
// buildAtomicDOFsVector(std::vector<Real> & v_out,
// FieldType field){
// switch(field){
// case _position0: buildAtomicDOFsVector<_position0 >(v_out); break;
// case _position: buildAtomicDOFsVector<_position >(v_out); break;
// case _displacement: buildAtomicDOFsVector<_displacement>(v_out); break;
// case _velocity: buildAtomicDOFsVector<_velocity >(v_out); break;
// case _force: buildAtomicDOFsVector<_force >(v_out);
// break;
// default: LM_FATAL("not handled field");
// }
// }
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::leastSquareAtomicData(
ContainerArray<Real> &data_mesh, ContainerArray<Real> &data_atom) {
if (comm_group_atomic == comm_group_continuum) {
this->solveLeastSquare(data_mesh, data_atom);
return;
}
/* do communication */
this->communicateBufferAtom2Continuum(data_atom);
if (this->is_in_continuum) {
this->solveLeastSquare(data_mesh, data_atom);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::leastSquareAtomicDOFs(
ContainerArray<Real> &data_mesh, FieldType field) {
LM_FATAL("to adapt");
if (comm_group_atomic == comm_group_continuum) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
extracted_field.build<typename DomainA::ContainerSubset>(this->pointList);
this->solveLeastSquare(data_mesh, extracted_field.getArray());
return;
}
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
if (this->is_in_atomic)
extracted_field.build<typename DomainA::ContainerSubset>(this->pointList);
LM_TOIMPLEMENT;
// extracted_field.resize(Parent::local_points * Dim);
this->leastSquareAtomicData(data_mesh, extracted_field.getArray());
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::buildAtomicDOFsNoiseVector(
std::vector<Real> &, FieldType field) {
ComputeExtract extracted_field("ComputeExtract:" + this->getID());
extracted_field.setParam("FIELD", field);
if (this->is_in_atomic)
extracted_field.build(*this->pointList);
if (this->is_in_continuum) {
this->interpolateAtomicDOFs(field, extracted_field.getArray());
for (UInt i = 0; i < this->local_points * DomainA::ContainerPoints::Dim;
++i)
LM_TOIMPLEMENT;
// extracted_field[i] *= -1;
}
/* do communication */
this->synchSumBuffer(extracted_field.getArray());
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
void BridgingAtomicContinuum<DomainA, DomainC>::checkPositionCoherency() {
constexpr UInt Dim = DomainA::ContainerPoints::Dim;
if (this->in_group_B) {
// interpolate atomic sites to check everything is ok from this
auto &&field = domC.getField(_position0);
UInt *assoc_tmp = &this->pointToElement[0];
for (UInt i = 0; i < this->local_points; ++i) {
Vector<Dim> X;
Vector<Dim> pos = VectorView<Dim>(&this->positions[Dim * i]);
X = this->getShapeMatrix().interpolate(i, field, assoc_tmp[i]);
for (UInt k = 0; k < Dim; ++k) {
if (X[k] != pos[k])
if (fabs(pos[k] - X[k]) > 1e-5) {
searchAtomInLocalPool(
static_cast<typename DomainA::ContainerSubset &>(
this->pointList),
pos);
DUMP("mismatch in positions (X), atom "
<< i << " dist is " << pos[k] - X[k]
<< " : probably a redistribution trouble check "
"migrations\n"
<< "computed position is " << X[0] << " " << X[1] << " "
<< X[2] << " position was " << pos[0] << " " << pos[1]
<< " " << pos[2],
DBG_MESSAGE);
}
}
}
}
this->distributeVectorB2A("temp_positions", this->positions);
if (this->in_group_A) {
UInt i = 0;
UInt global_flag = 0;
for (auto &&at : this->pointList) {
UInt local_flag = 0;
Vector<Dim> pos = VectorView<Dim>(&this->positions[Dim * i]);
Vector<Dim> at_pos = at.position0();
for (UInt k = 0; k < Dim; ++k) {
if (pos[k] != at_pos[k])
if (!local_flag && fabs(pos[k] - at_pos[k]) > 1e-5) {
searchAtomInLocalPool(this->pointList, pos);
DUMP("mismatch in positions (X), atom "
<< i << " dist is "
<< this->positions[Dim * i] - at.position0()[0]
<< " : probably a redistribution trouble check migrations"
<< std::endl
<< "ref is " << at << "\nand send position is " << pos[0]
<< " " << pos[1] << " " << pos[2],
DBG_MESSAGE);
local_flag = 1;
global_flag = 1;
}
}
++i;
}
if (global_flag) {
this->getDuoVector().print("buggy-redistrib");
LM_FATAL("abort due to migration problem");
}
}
DUMP("position coherency OK", DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC>
template <typename Cont, typename Vec>
void BridgingAtomicContinuum<DomainA, DomainC>::searchAtomInLocalPool(
Cont &container, Vec &x) {
UInt i = 0;
for (auto &&at : container) {
#ifdef TRACE_ATOM
auto pos = at.position0();
#endif
IF_COORD_EQUALS(pos, x) {
DUMP("actually searched atom is at ref " << at << " bridge index " << i,
DBG_MESSAGE);
break;
}
++i;
}
if (i == container.size())
DUMP("atom " << x[0] << " " << x[1] << " " << x[2] << " "
<< " not found locally : probably has migrated"
<< " or is simply lost !"
<< " NOT GOOD",
DBG_MESSAGE);
i = 0;
for (auto &&at : domA.getContainer()) {
#ifdef TRACE_ATOM
auto pos = at.position0();
#endif
IF_COORD_EQUALS(pos, x) {
DUMP("actually searched atom is at ref " << at << " bridge index " << i,
DBG_MESSAGE);
break;
}
++i;
}
if (i == domA.getContainer().size())
DUMP("atom not found locally : probably has migrated"
<< "or is simply lost !"
<< " NOT GOOD",
DBG_MESSAGE);
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(BridgingAtomicContinuum)
__END_LIBMULTISCALE__
c4science · Help