Page MenuHomec4science
Files F87712642

bridging_atomic_continuum.cc
No OneTemporary
Actions

File Metadata

Created
Mon, Oct 14, 10:56

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 "lm_common.hh"
#include "bridging_atomic_continuum.hh"
#include "lib_bridging.hh"
#include "compute_extract.hh"
#include "reference_manager.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename DomainA, typename DomainC,UInt Dim>
BridgingAtomicContinuum<DomainA,DomainC,Dim>::
BridgingAtomicContinuum(const std::string & name,
DomainA & A,DomainC & C):
Bridging<typename DomainA::ContainerAtoms,
typename DomainC::ContainerPoints,Dim>(name,
A.getContainer(),
C.getContainer()),
domA(A),domC(C) {
if (domA.getContainer().hasRefManager())
domA.getContainer().getRefManager().addSubSet("atomes_rec_subset",Parent::pointList);
field_sizes[0] = Dim;
field_sizes[1] = Dim;
field_sizes[2] = Dim;
field_sizes[3] = 1;
comm_group_continuum = domC.getGroupID();
comm_group_atomic = domA.getGroupID();
is_in_continuum = Communicator::getCommunicator().amIinGroup(comm_group_continuum);
is_in_atomic = Communicator::getCommunicator().amIinGroup(comm_group_atomic);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
BridgingAtomicContinuum<DomainA,DomainC,Dim>::~BridgingAtomicContinuum(){
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::clearAll(){
Parent::clearAll();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::updateForMigration(){
Parent::updateForMigration();
#ifndef LM_OPTIMIZED
checkPositionCoherency();
#endif //LM_OPTIMIZED
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::init(){
Parent::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,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::projectAtomicDOFsOnMesh(FieldType field){
if(is_in_continuum && !this->total_points)
DUMP("Warning : atomic container is empty: "
<< "cannot proceed atomic projection (check geometries ?!) "
<< this->getID(),DBG_WARNING);
if (!Parent::local_points) return;
Parent::allocateBuffer(Parent::local_points,Dim);
if(is_in_continuum) this->interpolateAtomicDOFs(field);
this->distributeVectorB2A("correction_fictifs",Parent::buffer,Dim);
if(is_in_atomic) setAtomicDOFs(field);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::projectAtomicDOFsOnMesh(FieldType field, std::vector<Real> & buffer){
if(is_in_continuum && !this->total_points)
DUMP("Warning : atomic container is empty: "
<< "cannot proceed atomic projection (check geometries ?!)",DBG_WARNING);
if(is_in_continuum) this->interpolateAtomicDOFs(field, buffer);
this->distributeVectorB2A("correction_fictifs",buffer,Dim);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::interpolateAtomicDOFs(FieldType field_code){
interpolateAtomicDOFs(field_code,Parent::buffer);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
interpolateAtomicDOFs(FieldType field_code, std::vector<Real> & buffer){
if (Parent::local_points == 0){
DUMP("Warning : atomic container is empty: "
<< "cannot proceed atomic projection (check geometries ?!)",DBG_WARNING);
return;
}
_Vec_ field;
if (field_code == _displacement) field = domC.getU();
else if (field_code == _velocity) field = domC.getV();
else LM_FATAL("field not handled yet " << field_code);
for (UInt i = 0 ; i < Parent::local_points ; ++i)
{
LM_ASSERT(i*Dim < buffer.size(),"overflow detected");
buffer[i*Dim] = this->getShapeMatrix().interpolate(i,field,Parent::assoc[i],0);
if (Dim > 1) buffer[i*Dim+1] = this->getShapeMatrix().interpolate(i,field,Parent::assoc[i],1);
if (Dim == 3)buffer[i*Dim+2] = this->getShapeMatrix().interpolate(i,field,Parent::assoc[i],2);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::setAtomicDOFs(FieldType field_code){
if (Parent::local_points == 0){
DUMP("Warning : atomic container is empty: "
<< "cannot proceed atomic projection (check geometries ?!)",DBG_WARNING);
return;
}
if (Parent::buffer.empty()){
LM_FATAL("function InterpolateAtomicDOFs should be called first! exiting");
}
typename BridgingAtomicContinuum<DomainA,DomainC,Dim>::IteratorAtomsSubset
it = Parent::pointList.getIterator();
typename BridgingAtomicContinuum<DomainA,DomainC,Dim>::RefAtom at = it.getFirst();
UInt at_index = 0;
// Real xmin,xmax,ymin,ymax,zmin,zmax;
// domA.getDomainDimensions(xmin,xmax,
// ymin,ymax,
// zmin,zmax);
// Real xsize = xmax - xmin;
// Real ysize = ymax - ymin;
for (at = it.getFirst(); !it.end() ; ++at_index, at = it.getNext())
{
at.force(0) = 0;
if (field_code == _velocity) at.velocity(0) = Parent::buffer[at_index*Dim];
if (field_code == _displacement){
at.position(0) = at.position0(0) + Parent::buffer[at_index*Dim];
//if (at.position(0) >= xmax) at.position(0) -= xsize;
//else if (at.position(0) < xmin) at.position(0) += xsize;
//****************************************************************
// Guillaume needs to include PBC treatment
//****************************************************************
at.velocity(0) = 0;
}
if (Dim > 1){
at.force(1) = 0;
if (field_code == _velocity) at.velocity(1) = Parent::buffer[at_index*Dim+1];
if (field_code == _displacement){
at.position(1) = at.position0(1) + Parent::buffer[at_index*Dim+1];
at.velocity(1) = 0;
//if (at.position(1) >= ymax) at.position(1) -= ysize;
//else if (at.position(1) < ymin) at.position(1) += ysize;
//****************************************************************
// Guillaume needs to include PBC treatment
//****************************************************************
}
}
if (Dim == 3){
at.force(2) = 0;
if (field_code == _velocity) at.velocity(2) = Parent::buffer[at_index*Dim+2];
if (field_code == _displacement) {
at.position(2) = at.position0(2) + Parent::buffer[at_index*Dim+2];
at.velocity(2) = 0;
}
}
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
averageDOFsOnElements(std::vector<Real> & data_mesh,
FieldType field){
LM_FATAL("to adapt");
if (comm_group_atomic == comm_group_continuum){
ComputeExtract<ContainerAtomsSubset> extracted_field("ComputeExtract:"+this->getID(),Dim);
extracted_field.setParam("FIELD",field);
extracted_field.build(Parent::pointList);
Parent::averageOnElements(extracted_field,data_mesh,Dim);
return;
}
ComputeExtract<ContainerAtomsSubset> extracted_field("ComputeExtract:"+this->getID(),Dim);
extracted_field.setParam("FIELD",field);
if(is_in_atomic)
extracted_field.build(Parent::pointList);
extracted_field.resize(Parent::local_points*Dim);
/* do communication */
this->communicateBufferAtom2Continuum(Dim,extracted_field);
if(is_in_continuum)
Parent::averageOnElements(extracted_field,data_mesh,Dim);
}
/* -------------------------------------------------------------------------- */
// template <typename DomainA, typename DomainC,UInt Dim>
// template <FieldType field>
// void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
// buildAtomicDOFsVector(std::vector<Real> & v_out){
// IteratorAtomsSubset it = Parent::pointList->getIterator();
// UInt nmax = Parent::pointList->nbElem();
// 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,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
leastSquareAtomicData(std::vector<Real> & data_mesh,
std::vector<Real> & data_atom,UInt stride){
if (comm_group_atomic == comm_group_continuum){
Parent::solveLeastSquare(data_mesh,data_atom,stride);
return;
}
/* do communication */
this->communicateBufferAtom2Continuum(stride,data_atom);
if(is_in_continuum){
Parent::solveLeastSquare(data_mesh,data_atom,stride);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
leastSquareAtomicDOFs(std::vector<Real> & data_mesh,
FieldType field){
LM_FATAL("to adapt");
if (comm_group_atomic == comm_group_continuum){
ComputeExtract<ContainerAtomsSubset> extracted_field("ComputeExtract:"+this->getID(),Dim);
extracted_field.setParam("FIELD",field);
extracted_field.build(Parent::pointList);
Parent::solveLeastSquare(data_mesh,extracted_field,Dim);
return;
}
ComputeExtract<ContainerAtomsSubset> extracted_field("ComputeExtract:"+this->getID(),Dim);
extracted_field.setParam("FIELD",field);
if(is_in_atomic)
extracted_field.build(Parent::pointList);
extracted_field.resize(Parent::local_points*Dim);
this->leastSquareAtomicData(data_mesh,extracted_field,Dim);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
buildAtomicDOFsNoiseVector(std::vector<Real> & v_out,
FieldType field){
ComputeExtract<ContainerAtomsSubset> extracted_field("ComputeExtract:"+this->getID(),Dim);
extracted_field.setParam("FIELD",field);
if(is_in_atomic)
extracted_field.build(Parent::pointList);
if(is_in_continuum){
this->interpolateAtomicDOFs(field,extracted_field);
for (UInt i = 0 ; i < Parent::local_points*Dim ; ++i)
extracted_field[i] *= -1;
}
/* do communication */
this->synchSumBuffer(Dim,extracted_field);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
attachVector(std::vector<Real> & tab,UInt stride){
if (domA.getContainer().hasRefManager())
domA.getContainer().getRefManager().attachVector(tab,Parent::pointList,stride);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::checkPositionCoherency(){
if(this->in_group_B){
//interpolate atomic sites to check everything is ok from this
typename DomainC::_Vec_ & field = domC.getP0();
UInt * assoc_tmp = &Parent::assoc[0];
for (UInt i = 0 ; i < Parent::local_points ; ++i){
{
Real x[3] = {0.0,0.0,0.0};
Real pos[3] = {0.0,0.0,0.0};
for (UInt k = 0; k < Dim; ++k) pos[k] = Parent::positions[Dim*i+k];
for (UInt k = 0; k < Dim; ++k){
x[k] = this->getShapeMatrix().interpolate(i,field,assoc_tmp[i],k);
if (x[k] != pos[k])
if (fabs(pos[k] - x[k]) > 1e-5) {
searchAtomInLocalPool(Parent::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",Parent::positions,Dim);
if(this->in_group_A){
IteratorAtomsSubset it = Parent::pointList.getIterator();
RefAtom at = it.getFirst();
UInt i = 0;
UInt global_flag = 0;
for (at = it.getFirst(); !it.end();at = it.getNext(),++i){
UInt local_flag = 0;
Real pos[3] = {0.0,0.0,0.0};
Real at_pos[3] = {0.0,0.0,0.0};
for (UInt k = 0; k < Dim; ++k) pos[k] = Parent::positions[Dim*i+k];
at.getPositions0(at_pos);
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(Parent::pointList,pos);
DUMP("mismatch in positions (X), atom " << i << " dist is "
<< Parent::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;
}
}
}
if (global_flag) {
Parent::getDuoVector().print("buggy-redistrib");
LM_FATAL("abort due to migration problem");
}
}
DUMP("position coherency OK",DBG_INFO);
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC,UInt Dim>
void BridgingAtomicContinuum<DomainA,DomainC,Dim>::
searchAtomInLocalPool(typename Parent::ContainerPointsSubset & container,Real * x){
UInt i = 0;
{
IteratorAtomsSubset it = container.getIterator();
RefAtom at = it.getFirst();
for (at = it.getFirst(), i = 0; !it.end() ; at = it.getNext(), ++i){
Real pos[3] = {0.0,0.0,0.0};
at.getPositions0(pos);
IF_COORD_EQUALS(pos,x) {
DUMP("actually searched atom is at ref " << at << " bridge index "
<< i,DBG_MESSAGE);
break;
}
}
if (it.end()) DUMP("atom " << x[0] << " " << x[1] << " " << x[2] << " "
<< " not found locally : probably has migrated"
<< " or is simply lost !"
<< " NOT GOOD",DBG_MESSAGE);
}
{
ContainerAtoms & c = domA.getContainer();
IteratorAtoms it = c.getIterator();
RefAtom at = it.getFirst();
for (at = it.getFirst(),i = 0; !it.end() ; at = it.getNext(), ++i){
Real pos[3] = {0.0,0.0,0.0};
at.getPositions0(pos);
IF_COORD_EQUALS(pos,x) {
DUMP("actually searched atom is at ref " << at << " bridge index "
<< i,DBG_MESSAGE);
break;
}
}
if (it.end()) DUMP("atom not found locally : probably has migrated"
<< "or is simply lost !"
<< " NOT GOOD",DBG_MESSAGE);
}
}
/* -------------------------------------------------------------------------- */
DECLARE_ATOMIC_CONTINUUM_TEMPLATE(BridgingAtomicContinuum)
__END_LIBMULTISCALE__

Event Timeline

c4science · Help