Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/template_cadd.cc
4577fdf5306cmaster

Last Change
Raw File

template_cadd.cc
View Options

/**
* @file template_cadd.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
*
* @date Thu Sep 18 14:51:07 2014
*
* @brief Simple CADD-like coupler
*
* @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 "template_cadd.hh"
#include "lib_bridging.hh"
#include <fstream>
#include <cmath>
#include <assert.h>
#include <unistd.h>
#include <algorithm>
#include <vector>
#include <utility>
#include "geometry_manager.hh"
#include "reference_manager_interface.hh"
#include "filter_manager.hh"
/* -------------------------------------------------------------------------- */
//added by guillaume
#include "spatial_grid_libmultiscale.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
typedef std::pair<double, double> point;
template<unsigned int dimension>
inline Real norm2(const Real * x1, const Real * x2){
Real distance = 0.0;
for (unsigned int i = 0; i < dimension; ++i) {
distance += square(x1[i]-x2[i]);
}
return distance;
}
inline Real sign(Real a)
{
if (a<0) return -1;
else if (a>0) return 1;
return 0;
}
const unsigned int dim = 2;
template <typename DomainDD, typename DomainMD,UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::init(){
// attach atom containers to refmanager
if (this->is_in_atomic){
ReferenceManagerInterface &ref_manager = this->domC.getContainer().getRefManager();
ref_manager.addSubSet("pad atoms",this->pad_atoms);
ref_manager.addSubSet("detection atoms",this->detection_atoms);
}
//extract nodes from domC.getContainerNodes();
IteratorAtoms it = this->domC.getContainer().getIterator();
IteratorDDNodes nodes = this->domA.getContainerNodes().getIterator();
Geometry* pad = GeometryManager::getManager().getGeometry(pad_geom);
Geometry* md = GeometryManager::getManager().getGeometry(md_geom);
Cube & md_cube = dynamic_cast<Cube&>(*md);
Geometry* detect = GeometryManager::getManager().getGeometry(detect_geom);
if (!dynamic_cast<Cube*>(md)) LM_FATAL("md geometry must be a cube");
/* -------------------------------------------------------------------------- */
//added by jcho
// read all the provided template files
for (unsigned int i = 0; i < this->num_templates; ++i) {
TwoDTemplate<RefAtom, IteratorAtomsSubset, false> * single_template
= new TwoDTemplate<RefAtom, IteratorAtomsSubset, false>();
this->core_templates.push_back(single_template);
this->core_templates[i]->init(this->templatefilename_list[i]);
}
/* -------------------------------------------------------------------------- */
UInt index = 0;
for (RefDDNode node = nodes.getFirst(); !nodes.end(); node = nodes.getNext()){
if (node.getConstraint() == 7) {
this->fictitious_nodes.add(node);
}
}
// added by jcho
// add pad nodes by computing intersection points
double x,y,z,nbr_x,nbr_z,nbr_y;
//RefDDNode nbr;
UInt nc = 0;
md_xrange[0] = md_cube.getXmin()[0];
md_xrange[1] = md_cube.getXmax()[0];
md_yrange[0] = md_cube.getXmin()[1];
md_yrange[1] = md_cube.getXmax()[1];
md_zrange[0] = md_cube.getXmin()[2];
md_zrange[1] = md_cube.getXmax()[2];
double pts0[2];
pts0[0] = md_cube.getXmax()[0];
pts0[1] = md_cube.getXmin()[2];
double pts1[2];
pts1[0] = md_cube.getXmax()[0];
pts1[1] = md_cube.getXmax()[2];
double pts2[2];
pts2[0] = md_cube.getXmin()[0];
pts2[1] = md_cube.getXmax()[2];
double pts3[2];
pts3[0] = md_cube.getXmin()[0];
pts3[1] = md_cube.getXmin()[2];
double ptsx,ptsz,nextidx;
std::vector<int> copied_indices;
double current_DD[3];
double next_DD[3];
double dx,dy,dz;
double bx,by,bz,br;
double dangle,dangle_deg,distance,cangle;
double lx,ly,lz,lr;
double maxseg;
for (RefDDNode node = nodes.getFirst(); !nodes.end(); node = nodes.getNext()){
x = node.position(0);
y = node.position(1);
z = node.position(2);
index = node.getIndex();
maxseg = node.getmaxSeg();
// hard code: y is normal to slip direction!
// we do not care y-coordinate whethere intersect MD boundaries or not
current_DD[0] = x;
current_DD[1] = z;
current_DD[2] = y;
nc = node.numNbrs();
if(nc>=3){
LM_FATAL("node (x:" << x << ", y:" << y << ", z:" << z << ")"
<<" has more than 2 neighbors ("<< nc <<" neighbors)."
<<" Currently, we only support upto 2 neighbors.");
}
for (unsigned int n = 0; n < nc; n++){
nbr_x = node.getNeighborNode(n,0);
nbr_z = node.getNeighborNode(n,2);
nbr_y = node.getNeighborNode(n,1);
nextidx = node.getNeighborNodeIndex(n);
// hard code: y is normal to slip direction!
// we do not care y-coordinate whethere intersect MD boundaries or not
next_DD[0] = nbr_x;
next_DD[1] = nbr_z;
next_DD[2] = nbr_y;
dx = next_DD[0] - current_DD[0];
dy = next_DD[1] - current_DD[1];
dz = next_DD[2] - current_DD[2];
lr = sqrt(dx*dx + dy*dy + dz*dz);
if(lr > maxseg){
// the connectivity over the periodic boundary condition of DD box
}
else if( (std::find(copied_indices.begin(),copied_indices.end(),index) ==copied_indices.end()) &&
(std::find(copied_indices.begin(),copied_indices.end(),nextidx)==copied_indices.end())) {
// intersect the MD boundary: pts0-pts1 || pts1-pts2 || pts2-pts3 || pts3-pts0
if ( (intersectTwoD(current_DD,next_DD,pts0,pts1, &ptsx, &ptsz)) ||
(intersectTwoD(current_DD,next_DD,pts1,pts2, &ptsx, &ptsz)) ||
(intersectTwoD(current_DD,next_DD,pts2,pts3, &ptsx, &ptsz)) ||
(intersectTwoD(current_DD,next_DD,pts3,pts0, &ptsx, &ptsz)) ){
copied_indices.push_back(index);
lx = dx/lr; ly = dy/lr; lz = dz/lr;
bx = node.burgersVectorX(n);
by = node.burgersVectorY(n);
bz = node.burgersVectorZ(n);
br = sqrt(bx*bx+by*by+bz*bz);
bx/=br; by/=br; bz/=br; /* unit vector along Burgers vector */
dangle = fabs(bx*lx+by*ly+bz*lz);
dangle_deg = acos(dangle)*180.0/3.14159265359;
distance = fabs(90.0 - dangle_deg);
cangle = distance + 90.0;
std::vector<double> temp;
temp.push_back(ptsx);
temp.push_back(0.0);
temp.push_back(ptsz);
temp.push_back(cangle);
this->pad_dd.push_back(temp);
}
}
}
}
for (RefDDNode node = nodes.getFirst(); !nodes.end(); node = nodes.getNext()){
index = node.getIndex();
if(std::find(copied_indices.begin(),copied_indices.end(),index)!=copied_indices.end()){
this->pad_nodes.add(node);
}
}
if(this->pad_nodes.size() != this->pad_dd.size()){
LM_FATAL("something wrong during detecting pad_dd nodes" );
}
for (RefAtom at = it.getFirst(); !it.end(); at = it.getNext()){
DUMP("Testing atom at position " << at.position0(0)
<< " " << at.position0(1)
<< " " << at.position0(2),DBG_ALL);
bool inboth = false;
if (pad->contains(at.position0(0),at.position0(1),at.position0(2))){
DUMP("Accepting atom at position " << at.position0(0),DBG_DETAIL);
this->pad_atoms.add(at);
inboth = true;
} else {
inboth = false;
}
if(detect->contains(at.position0(0),at.position0(1),at.position0(2))) {
this->detection_atoms.add(at);
for (unsigned int i = 0; i < 3; ++i) {
assert(!isnan(at.position(i)));
}
inboth &= true;
} else {
inboth = false;
}
if (inboth){
LM_FATAL("An atom is both in the dislocation detection zone AND the pad"
<<" , which is conceptually inadmissible. Make sure, "
<< "your pad_geom = " << pad_geom << " and detection_geom = "
<< detect_geom << " are NONoverlapping geometries");
}
}
UInt local_points = this->pad_atoms.nbElem();
if (!local_points){
DUMP("We found no atom in the pad",DBG_INFO);
}
DUMP("We found " << local_points << " atoms in the bridging pad",DBG_INFO_STARTUP);
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD,UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::coupling(CouplingStage stage)
{
if (stage != COUPLING_STEP3) return;// WILL MAYBE CHANGE!!!!!
//! Main coupling calls
STARTTIMER("dislo_detection");
this->applyDetection();
STOPTIMER("dislo_detection");
STARTTIMER("template_application");
this->applyTemplateToPad();
STOPTIMER("template_application");
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::applyDetection()
{
if (current_step % this->detect_freq == 0){
// For the moment, detect MD dislocation core using centro-symmetry,
// but it need to be replaced to DXA method (Max's work).
this->detect_cores_by_compute();
if (this->is_in_dd){
// if(this->num_cores != this->fictitious_nodes.size()){
// // do something ... insertion algorithm(?)
// LM_FATAL("not match number of cores MD:" << this->num_cores << " DD:" <<this->fictitious_nodes.size());
// }
double mindist,dist;
double ndx = 0.0; double ndy = 0.0; double ndz = 0.0;
IteratorDDNodesSubset it_fict = this->fictitious_nodes.getIterator();
for (RefDDNode nd = it_fict.getFirst(); !it_fict.end(); nd = it_fict.getNext()) {
mindist = 100.0;
for (int i=0; i<this->num_cores; ++i){
dist = sqrt( (nd.position(0)-this->pos_core[i][0])*(nd.position(0)-this->pos_core[i][0]) +
(nd.position(1)-this->pos_core[i][1])*(nd.position(1)-this->pos_core[i][1]) +
(nd.position(2)-this->pos_core[i][2])*(nd.position(2)-this->pos_core[i][2]) );
if(dist < mindist){
mindist = dist;
ndx = this->pos_core[i][0];
ndy = this->pos_core[i][1];
ndz = this->pos_core[i][2];
}
}
nd.position(0) = ndx;
nd.position(1) = ndy;
nd.position(2) = ndz;
}
}
}
}
/* -------------------------------------------------------------------------- */
//added by jcho
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::detect_cores_by_compute()
{
UInt num_global = 0;
if (this->is_in_atomic){
double csd_val;
int atom_index;
std::vector<ComputeInterface*> computes_ptr;
computes_ptr.resize(compute_list.size());
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
FilterID id = compute_list[comp];
ComputeInterface * my_compute = dynamic_cast<ComputeInterface*>(FilterManager::getManager().getObject(id));
if (!my_compute)
LM_FATAL("invalid compute passed to dumper text " << compute_list[0]);
computes_ptr[comp] = my_compute;
my_compute->build();
ComputeCompatibility<> * out_compute = new ComputeCompatibility<>("",computes_ptr[comp]->getDim());
out_computes.push_back(out_compute);
}
std::vector<std::vector<double> > dislocated_atoms;
IteratorAtomsSubset it_detection_atoms = this->detection_atoms.getIterator();
for (RefAtom at = it_detection_atoms.getFirst();
!it_detection_atoms.end();
at = it_detection_atoms.getNext()){
atom_index = at.getIndex();
csd_val = computes_ptr[0]->get(atom_index);
if ((this->epsilon[0]<=csd_val) && (csd_val<this->epsilon[1])){
std::vector<double> coordinate;
coordinate.push_back(at.position(0));
coordinate.push_back(at.position(1));
coordinate.push_back(at.position(2));
dislocated_atoms.push_back(coordinate);
}
}
std::vector<std::vector<double> > dislocations;
cores_by_distance(dislocated_atoms, &dislocations);
this->num_cores = dislocations.size();
for (int i=0; i<dislocations.size(); ++i){
Quantity <Length,3> temp;
temp[0] = dislocations[i][0];
temp[1] = dislocations[i][1];
temp[2] = dislocations[i][2];
this->pos_core.push_back(temp);
}
num_global = dislocations.size();
DUMP("total " << num_global << " dislocations are detected in MD processor.",DBG_MESSAGE);
if (Communicator::getCommunicator().groupRank(lm_my_proc_id,
this->comm_group_atomic) == 0)
Communicator::getCommunicator().send(&num_global,1,0,this->comm_group_dd,"number of dislocation cores");
}
else{
Communicator::getCommunicator().receive(&num_global,1,0,this->comm_group_atomic,"number of dislocation cores");
DUMP("total " << num_global << " dislocations are received in DD processor.",DBG_MESSAGE);
}
}
/* -------------------------------------------------------------------------- */
//added by jcho
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::cores_by_distance(std::vector<std::vector<double> >& atoms,
std::vector<std::vector<double> > * dislocations)
{
bool finish;
std::vector<std::vector<double> > clusters_x;
std::vector<std::vector<double> > clusters_y;
std::vector<std::vector<double> > clusters_z;
for (int i = 0; i < atoms.size(); ++i){
std::vector<double> coordinate_x;
std::vector<double> coordinate_y;
std::vector<double> coordinate_z;
coordinate_x.push_back(atoms[i][0]);
coordinate_y.push_back(atoms[i][1]);
coordinate_z.push_back(atoms[i][2]);
clusters_x.push_back(coordinate_x);
clusters_y.push_back(coordinate_y);
clusters_z.push_back(coordinate_z);
}
finish = false;
while (finish == false){
finish = all_sorted(clusters_x, clusters_y, clusters_z);
clustering_between_groups(&clusters_x,&clusters_y,&clusters_z);
}
for (int i = 0; i < clusters_x.size(); ++i){
if(clusters_x[i].size()>1){
std::vector<double> coordinate;
double center_x;
double center_y;
double center_z;
for (int j=0; j < clusters_x[i].size(); ++j){
center_x += clusters_x[i][j];
center_y += clusters_y[i][j];
center_z += clusters_z[i][j];
}
center_x /= clusters_x[i].size();
center_y /= clusters_y[i].size();
center_z /= clusters_z[i].size();
coordinate.push_back(center_x);
coordinate.push_back(center_y);
coordinate.push_back(center_z);
dislocations->push_back(coordinate);
}
}
}
/* -------------------------------------------------------------------------- */
//added by jcho
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::clustering_between_groups(std::vector<std::vector<double> > * clusters_x,
std::vector<std::vector<double> > * clusters_y,
std::vector<std::vector<double> > * clusters_z)
{
std::vector<int> visited_index;
std::vector<std::vector<double> > modified_clusters_x;
std::vector<std::vector<double> > modified_clusters_y;
std::vector<std::vector<double> > modified_clusters_z;
std::vector<double> test_master_group_x;
std::vector<double> test_master_group_y;
std::vector<double> test_master_group_z;
std::vector<double> test_slave_group_x;
std::vector<double> test_slave_group_y;
std::vector<double> test_slave_group_z;
std::vector<double> a_cluster_x;
std::vector<double> a_cluster_y;
std::vector<double> a_cluster_z;
std::vector<double> distances;
double rcut,temp;
rcut = this->cluster_radius;
for (int i=0; i < clusters_x[0].size(); ++i){
if (std::find(visited_index.begin(), visited_index.end(), i) == visited_index.end()){
visited_index.push_back(i);
test_master_group_x = clusters_x[0][i];
test_master_group_y = clusters_y[0][i];
test_master_group_z = clusters_z[0][i];
a_cluster_x = test_master_group_x;
a_cluster_y = test_master_group_y;
a_cluster_z = test_master_group_z;
for (int j=i+1; j < clusters_x[0].size(); ++j){
if (std::find(visited_index.begin(), visited_index.end(), j) == visited_index.end()){
test_slave_group_x = clusters_x[0][j];
test_slave_group_y = clusters_y[0][j];
test_slave_group_z = clusters_z[0][j];
distances = distances_two_groups(test_master_group_x,test_master_group_y,test_master_group_z,
test_slave_group_x, test_slave_group_y, test_slave_group_z);
temp = *min_element(distances.begin(), distances.end());
if(temp < rcut){
visited_index.push_back(j);
for (int z=0; z < test_slave_group_x.size(); ++z){
a_cluster_x.push_back(test_slave_group_x[z]);
a_cluster_y.push_back(test_slave_group_y[z]);
a_cluster_z.push_back(test_slave_group_z[z]);
}
}
}
}
modified_clusters_x.push_back(a_cluster_x);
modified_clusters_y.push_back(a_cluster_y);
modified_clusters_z.push_back(a_cluster_z);
}
}
clusters_x[0] = modified_clusters_x;
clusters_y[0] = modified_clusters_y;
clusters_z[0] = modified_clusters_z;
}
/* -------------------------------------------------------------------------- */
//added by jcho
template <typename DomainDD, typename DomainMD, UInt Dim>
bool CADD_template<DomainDD,DomainMD,Dim>::all_sorted(std::vector<std::vector<double> >& groups_x,
std::vector<std::vector<double> >& groups_y,
std::vector<std::vector<double> >& groups_z)
{
double mindist,rcut,temp;
std::vector<double> current_group_x,next_group_x;
std::vector<double> current_group_y,next_group_y;
std::vector<double> current_group_z,next_group_z;
std::vector<double> distances;
mindist = 100.0;
rcut = this->cluster_radius;
for (int i = 0; i < groups_x.size(); ++i){
current_group_x = groups_x[i];
current_group_y = groups_y[i];
current_group_z = groups_z[i];
for (int j = i+1; j < groups_x.size(); ++j){
next_group_x = groups_x[j];
next_group_y = groups_y[j];
next_group_z = groups_z[j];
distances = distances_two_groups(current_group_x,current_group_y,current_group_z,
next_group_x, next_group_y, next_group_z);
temp = *min_element(distances.begin(), distances.end());
if(temp < mindist){
mindist = temp;
}
}
}
if(mindist <= rcut){
return false;
}
else{
return true;
}
}
/* -------------------------------------------------------------------------- */
//added by jcho
template <typename DomainDD, typename DomainMD, UInt Dim>
std::vector<double>
CADD_template<DomainDD,DomainMD,Dim>::distances_two_groups(std::vector<double> & group1_x,
std::vector<double> & group1_y,
std::vector<double> & group1_z,
std::vector<double> & group2_x,
std::vector<double> & group2_y,
std::vector<double> & group2_z)
{
std::vector<double> dists;
double dist;
for (int i=0; i < group1_x.size(); ++i){
for (int j=0; j < group2_x.size(); ++j){
dist = sqrt((group1_x[i]-group2_x[j])*(group1_x[i]-group2_x[j])+
(group1_y[i]-group2_y[j])*(group1_y[i]-group2_y[j])+
(group1_z[i]-group2_z[j])*(group1_z[i]-group2_z[j]));
dists.push_back(dist);
}
}
return dists;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::applyTemplateToPad()
{
int position; //0: pure core template, 1; core template + barnett, 2: barnett
double pos_bar[3];
double barnett_disp[] = {0.0, 0.0, 0.0};
double factor,angle;
double two_angles[2];
int two_indices[2];
UInt counter = 0;
// if (this->is_in_atomic) {
std::vector<ComputeInterface*> computes_ptr;
computes_ptr.resize(1);
IteratorAtomsSubset it_pad_atoms = this->pad_atoms.getIterator();
std::vector<Real> allDisplacements(3*this->pad_atoms.size());
if (current_step % detect_freq == 0){
if(barnett_compute.size()){
FilterID id = barnett_compute[0];
ComputeInterface * my_compute = dynamic_cast<ComputeInterface*>(FilterManager::getManager().getObject(id));
if (!my_compute)
LM_FATAL("invalid compute for barnett displacements " << barnett_compute[0]);
computes_ptr[0] = my_compute;
DUMP("start to compute Barnett displacement fields.",DBG_MESSAGE);
my_compute->build();
DUMP("finish to compute Barnett displacement fields.",DBG_MESSAGE);
ComputeCompatibility<> * out_compute = new ComputeCompatibility<>("",computes_ptr[0]->getDim());
out_computes.push_back(out_compute);
}
else{
//barnett_disp_for_all_padatoms(it_pad_atoms, allDisplacements);
}
counter = 0;
typename ComputeInterface::iterator comp_vals = computes_ptr[0]->getIterator();
for(double val = comp_vals.getFirst(); !comp_vals.end(); val = comp_vals.getNext()) {
allDisplacements[counter] = val;
counter++;
}
}
if (this->is_in_atomic) {
counter = 0;
DUMP("start to constrain pad atoms.",DBG_MESSAGE);
for (RefAtom at = it_pad_atoms.getFirst(); !it_pad_atoms.end(); at = it_pad_atoms.getNext()){
for (unsigned int i = 0; i < 3; ++i) {
at.velocity(i) = 0.0; //constraint the movements of pad atoms
barnett_disp[i] = allDisplacements[counter+i];//extract the computed barnett displacement fields
}
counter += 3;
if (current_step % detect_freq == 0){
position_of_atom(at.position0(0), at.position0(1), at.position0(2), &position, pos_bar, &factor, &angle);
compute_angle_ratio(angle, two_indices, two_angles);
applyDisplacement(at, angle, pos_bar,position, factor, this->detect_freq, two_indices, two_angles, barnett_disp);
}
}
DUMP("finish to constrain pad atoms.",DBG_MESSAGE);
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::position_of_atom(double x,
double y,
double z,
int *position,
double *pos_bar,
double *factor,
double *angle)
{
int atom_region,core_region;
double pad_ddx,pad_ddz,pad_ddy;
double temp_pos_bar[3];
double temp_factor;
double pts0[2],pts1[2],pts2[2],pts3[2],center[2];
pts0[0] = this->md_xrange[1]*2.0; pts0[1] = this->md_zrange[0]*2.0;
pts1[0] = this->md_xrange[1]*2.0; pts1[1] = this->md_zrange[1]*2.0;
pts2[0] = this->md_xrange[0]*2.0; pts2[1] = this->md_zrange[1]*2.0;
pts3[0] = this->md_xrange[0]*2.0; pts3[1] = this->md_zrange[0]*2.0;
center[0] = (this->md_xrange[1]+this->md_xrange[1]+this->md_xrange[0]+this->md_xrange[0])/4.0;
center[1] = (this->md_zrange[0]+this->md_zrange[1]+this->md_zrange[1]+this->md_zrange[0])/4.0;
atom_region = evaluate_edge(x,y,z,pts0,pts1,pts2,pts3,center);
*factor = 1.0;
pos_bar[0] = x;
pos_bar[1] = y;
pos_bar[2] = z;
for (int i=0; i < this->pad_dd.size(); ++i){
pad_ddx = pad_dd[i][0];
pad_ddy = pad_dd[i][1];
pad_ddz = pad_dd[i][2];
*angle = pad_dd[i][3];
core_region = evaluate_edge(pad_ddx,pad_ddy,pad_ddz,pts0,pts1,pts2,pts3,center);
if(atom_region != core_region){
*position = 2;
}
else{
if( (core_region==1) || (core_region==3)) {
temp_pos_bar[0] = z-pad_dd[i][2];
temp_pos_bar[1] = y-pad_dd[i][1];
temp_pos_bar[2] = x-pad_dd[i][0];
}
else {
temp_pos_bar[0] = x-pad_dd[i][0];
temp_pos_bar[1] = y-pad_dd[i][1];
temp_pos_bar[2] = z-pad_dd[i][2];
}
if (in_side_ellipse( temp_pos_bar, this->rmax, this->rmin, &temp_factor)){
if(temp_factor <= this->factor){
*position = 0;
pos_bar[0] = temp_pos_bar[0];
pos_bar[1] = temp_pos_bar[1];
pos_bar[2] = temp_pos_bar[2];
*factor = temp_factor;
return;
}
else{
*position = 1;
pos_bar[0] = temp_pos_bar[0];
pos_bar[1] = temp_pos_bar[1];
pos_bar[2] = temp_pos_bar[2];
*factor = temp_factor;
}
}
else{
if(*position != 1){
*position = 2;
}
else{
}
}
}
}
return;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::compute_angle_ratio(double angle, int *two_indices,double *two_angles)
{
int j;
if(angle < 90.0){
angle = 180.0-angle;
}
else if ((180.0<angle) && (angle<=270.0)){
angle = 360.0-angle;
}
else if ((270.0<angle) && (angle<=360.0)){
angle = angle-180.0;
}
for (j=0; j < this->num_templates-1; ++j){
two_angles[0] = this->core_templates[j]->get_theta();
two_angles[1] = this->core_templates[j+1]->get_theta();
two_indices[0] = j;
two_indices[1] = j+1;
if((angle-two_angles[0])*(angle-two_angles[1])<=0.0){
break;}
}
return;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
bool CADD_template<DomainDD,DomainMD,Dim>::in_side_ellipse(double pos[], double rmax, double rmin, double *factor)
{
double eq,val;
eq = ( (pos[0]/rmax)*(pos[0]/rmax)+(pos[1]/rmin)*(pos[1]/rmin) );
val = 1.0 - eq;
if(val >= 0.0){
*factor = sqrt(eq);
return true;
}
else{
*factor = sqrt(eq);
return false;
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
int CADD_template<DomainDD,DomainMD,Dim>::evaluate_edge(double x,
double y,
double z,
double pts1[],
double pts2[],
double pts3[],
double pts4[],
double cent[])
{
double point[2];
point[0] = x;
point[1] = z;
if(is_in_triangle(cent,pts1,pts2,point)){
return 1;
}
else if(is_in_triangle(cent,pts2,pts3,point)){
return 2;
}
else if(is_in_triangle(cent,pts3,pts4,point)){
return 3;
}
else if(is_in_triangle(cent,pts4,pts1,point)){
return 4;
}
int * p = NULL;
return *p;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::shear_strain(double pos[],
double disp[],
double mu,
double sigma_yx,
double sigma_yz)
{
double eps_yx,eps_yz;
eps_yx = sigma_yx/mu;
eps_yz = sigma_yz/mu;
disp[0] = eps_yx/pos[1];
disp[0] = 0.0;
disp[0] = eps_yz/pos[1];
}
/* -------------------------------------------------------------------------- */
template<typename base>
void alloc(base**& table, base*& container, unsigned int nrows, unsigned int stride)
{
table = new base*[nrows];
container = new base[stride*nrows];
for (unsigned int i = 0; i < nrows; ++i) {
table[i] = container+i*stride;
}
}
/* -------------------------------------------------------------------------- */
////////////////////////////////////////////////////////////////////////
template <typename RefAtom, typename IteratorAtoms,bool use_template>
void TwoDTemplate<RefAtom, IteratorAtoms, use_template>::init(const std::string & filename)
{
// check search_dim
if ((this->search_dim > 2) || (this->search_dim < 1)) {
LM_THROW("SEARCH_DIM can only be 1 or 2");
}
this->atoms_around_core = new ContainerAtomsSubset;
alloc(this->vertices, this->vertices_cont, this->search_dim+1, 2);
std::ifstream infile;
// infile.exceptions ( ifstream::failbit);
infile.open(filename.c_str());
char tmp;
enum read_state {NOT_READING, READ_MATERIAL, READ_ATOMS, READ_CONNECTIVITIES};
unsigned int buflen = 1024;
read_state state=NOT_READING;
if (!infile.good()){
char buf[256];
char * path = buf;
UInt size = 0;
path = getcwd(buf,size);
LM_FATAL("The template file \""<< filename.c_str()
<< "\" could NOT be opened, rdstate = " << std::hex
<< infile.rdstate() << std::endl
<< " badbit = " << infile.bad()<<std::endl
<< " failbit = " << infile.fail()<<std::endl
<< " eofbit = " << infile.eof()<<std::endl
<< " goodbit = " << infile.good()<<std::endl
<< "The current working directory is \""
<< path << "\"" << std::endl)
}
while (infile.good()){
tmp = infile.peek();
switch (tmp){
case '#':
break;
default:
if (state==NOT_READING){
state = READ_MATERIAL;
infile >> this->poisson >> this->burgers >> this->theta;
}
else if (state == READ_MATERIAL){
state = READ_ATOMS;
infile >> this->num_nodes;
alloc(this->coords, this->coords_cont, this->num_nodes, 2);
for (unsigned int i = 0; i < this->num_nodes; ++i) {
infile.ignore(buflen, '\n');
infile>>this->coords[i][0]>>this->coords[i][1];
}
}
else if (state == READ_ATOMS){
state = READ_CONNECTIVITIES;
infile>>this->num_elems;
alloc(this->connectivity, this->connectivity_cont, this->num_elems, 3);
alloc(this->shape, this->shape_cont, this->num_elems, 9);
alloc(this->x_l_transition, this->x_l_transition_cont, this->num_elems, 9);
for (unsigned int i = 0; i < this->num_elems; ++i) {
infile.ignore(buflen, '\n');
infile>>this->connectivity[i][0]
>>this->connectivity[i][1]
>>this->connectivity[i][2];
infile.ignore(buflen, '\n');
infile>>this->shape[i][0]
>>this->shape[i][1]
>>this->shape[i][2]
>>this->shape[i][3]
>>this->shape[i][4]
>>this->shape[i][5]
>>this->shape[i][6]
>>this->shape[i][7]
>>this->shape[i][8];
infile.ignore(buflen, '\n');
infile>>this->x_l_transition[i][0]
>>this->x_l_transition[i][1]
>>this->x_l_transition[i][2]
>>this->x_l_transition[i][3]
>>this->x_l_transition[i][4]
>>this->x_l_transition[i][5]
>>this->x_l_transition[i][6]
>>this->x_l_transition[i][7]
>>this->x_l_transition[i][8];
}
}
break;
}
infile.ignore(buflen, '\n');
}
infile.close();
}
inline void template_disp(Real a[], Real b[], Real c[], Real Su[], Real SxInv[], Real atom_coords[], Real disp[])
{
UInt index[3][3] = {{0,1,2},{3,4,5},{6,7,8}};
for (UInt i = 0; i < 3; ++i) {
disp[i]=0.0;
for (UInt j = 0; j < 3; ++j) {
for (UInt k = 0; k < 3; ++k) {
disp[i] += Su[index[i][j]] * SxInv[index[j][k]] * atom_coords[k];
}
}
}
}
template <typename RefAtom, typename IteratorAtoms,bool use_template>
void TwoDTemplate<RefAtom, IteratorAtoms, use_template>::coreDispAtom(double pos_core[], double disp[])
{
int elem_num = findElement(pos_core);
if(elem_num != -1){
double atom_coord[] = {1.0, pos_core[0], pos_core[1]};
template_disp(this->coords[this->connectivity[elem_num][0]],
this->coords[this->connectivity[elem_num][1]],
this->coords[this->connectivity[elem_num][2]],
this->shape[elem_num], this->x_l_transition[elem_num],
atom_coord, disp);
}
else{
LM_FATAL("cannot find element!");
}
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::barnett_disp_for_all_padatoms(IteratorAtomsSubset it_pad_atoms, std::vector<Real> & allDisplacements){
UInt take_paradis_ptr = 0;
DomainDDInterface *_ptr;
if (this->is_in_dd){
DomainInterface *domDD = DomainMultiScale::getManager().getObject(domdd);
*_ptr = dynamic_cast<DomainDDInterface*>(domDD);
Communicator::getCommunicator().send(&_ptr,1,0,this->comm_group_dd,"paradis pointer");
}
else{
Communicator::getCommunicator().receive(&_ptr,1,0,this->comm_group_atomic,"paradis pointer");
}
UInt num_pad_atoms;
if (this->is_in_atomic) {
std::vector<Real> allPositions(3*this->pad_atoms.size());
int count = 0;
UInt total_count = 0;
for (RefAtom at = it_pad_atoms.getFirst(); !it_pad_atoms.end(); at = it_pad_atoms.getNext()){
allPositions[count] = at.position0(0);
allPositions[count+1] = at.position0(1);
allPositions[count+2] = at.position0(2);
count += 3;
total_count += 1;
}
_ptr->computeDisplacements(allPositions,allDisplacements);
num_pad_atoms = total_count;
Communicator::getCommunicator().send(&num_pad_atoms,1,0,this->comm_group_dd,"total number of local pad atoms");
}
else{
Communicator::getCommunicator().receive(&num_pad_atoms,1,0,this->comm_group_atomic,"total number of local pad atoms");
}
return;
}
/* -------------------------------------------------------------------------- */
template <typename DomainDD, typename DomainMD, UInt Dim>
void CADD_template<DomainDD,DomainMD,Dim>::applyDisplacement(RefAtom at,
double angle,
double* pos_bar,
int position,
double f,
int detect_freq,
int *two_indices,
double *two_angles,
double disp_barnett[])
{
STARTTIMER("application");
Real disp_templates[] = {0.0, 0.0, 0.0};
Real disp_strain[] = {0.0, 0.0, 0.0};
Real atom_coords[3];
// Real dx[] = {0.0, 0.0};
Real dz[] = {0.0, 0.0};
// Real total_disp[] = {0.0, 0.0, 0.0};
double factor_a = -1.0/this->factor*(f-1.0);
double factor_b = 1.0/this->factor*(f-(1.0-this->factor));
for (unsigned int i = 0; i < 3; ++i) {
atom_coords[i] = at.position0(i);
}
shear_strain(atom_coords,disp_strain,this->mu,this->shear_stress_yx,this->shear_stress_yz);
if ((position == 0) || (position == 1)){
for (int j=0; j < 2; ++j) {
Real disp_template[] = {0.0, 0.0, 0.0};
this->core_templates[two_indices[j]]->coreDispAtom(pos_bar, disp_template);
// dx[j] = disp_template[2]*sin(two_angles[j]*M_PI/180.0) + disp_template[0]*cos(two_angles[j]*M_PI/180.0);
dz[j] = disp_template[0]*sin(two_angles[j]*M_PI/180.0) + disp_template[2]*cos(two_angles[j]*M_PI/180.0);
}
disp_templates[0] = disp_barnett[0];//(dx[1]-dx[0])/(two_angles[1]-two_angles[0])*(angle-two_angles[0])+dx[0];
disp_templates[1] = disp_barnett[1];
disp_templates[2] = (dz[1]-dz[0])/(two_angles[1]-two_angles[0])*(angle-two_angles[0])+dz[0];
}
if(position == 0){
at.position(0) = at.position0(0) + disp_templates[0] + disp_strain[0];
at.position(1) = at.position0(1) + disp_templates[1] + disp_strain[1];
at.position(2) = at.position0(2) + disp_templates[2] + disp_strain[2];
}
else if(position ==1){
at.position(0) = at.position0(0) + disp_templates[0]*factor_a + disp_barnett[0]*factor_b + disp_strain[0];
at.position(1) = at.position0(1) + disp_templates[1]*factor_a + disp_barnett[1]*factor_b + disp_strain[1];
at.position(2) = at.position0(2) + disp_templates[2]*factor_a + disp_barnett[2]*factor_b + disp_strain[2];
}
else if(position ==2){
at.position(0) = at.position0(0) + disp_barnett[0] + disp_strain[0];
at.position(1) = at.position0(1) + disp_barnett[1] + disp_strain[1];
at.position(2) = at.position0(2) + disp_barnett[2] + disp_strain[2];
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC CADDTEMPLATE
CADDTEMPLATE is a simple CADD-like coupler which can deal with straing egde
dislocations normal to the interface (i.e., 1;3400;0cno elastic problem)
*/
/* LMEXAMPLE
COUPLING_CODE
*/
template <typename DomainDD, typename DomainMD,unsigned int Dim>
void CADD_template<DomainDD,DomainMD,Dim>::declareParams(){
/* LMKEYWORD PADGEOM
Identifier of the pad geometry
*/
this->parseKeyword("PADGEOM",pad_geom);
/* LMKEYWORD SLIP_PLANE
Specifies the normal vector to the slip plane
*/
this->parseVectorKeyword("SLIP_PLANE", 3, this->normal_quant);
/* LMKEYWORD MDGEOM
Identifier of the MD geometry
*/
this->parseKeyword("MDGEOM", this->md_geom);
/* LMKEYWORD DETECTGEOM
Identifier of the geometry in which dislos are detected
*/
this->parseKeyword("DETECTGEOM", this->detect_geom);
/* LMKEYWORD GRID_SIZEX
Specifies the size of grid cells in X directions
*/
this->parseKeyword("GRID_SIZEX", this->grid_division[0]);
/* LMKEYWORD GRID_SIZEY
Specifies the size of grid cells in Y directions
*/
this->parseKeyword("GRID_SIZEY", this->grid_division[1]);
/* LMKEYWORD GRID_SIZEZ
Specifies the number of grid cells in the Z direction
*/
this->parseKeyword("GRID_SIZEZ", this->grid_division[2]);
/* LMKEYWORD DETECTION_FREQUENCY
Specifies the frequency at which detection computations are performed
*/
this->parseKeyword("DETECTION_FREQUENCY", this->detect_freq);
/* LMKEYWORD DETECT_COMPUTE
Add computes to input to filter
*/
this->parseKeyword("DETECTION_COMPUTE",compute_list);
/* LMKEYWORD BARNETT_COMPUTE
Add computes having Barnett displacement fields
*/
this->parseKeyword("BARNETT_COMPUTE",barnett_compute);
/* LMKEYWORD EPSILON
Specifies the lower and upper bounds for the compute result to
identify dislocations.
*/
this->parseVectorKeyword("EPSILON",2,epsilon);
/* LMKEYWORD CLUSTER_RADIUS
Specifies the radius to identify dislocation core.
*/
this->parseKeyword("CLUSTER_RADIUS",cluster_radius);
/* LMKEYWORD NUM_TEMPLATES
Specifies the number of templates provided.
*/
this->parseKeyword("NUM_TEMPLATES", this->num_templates);
/* LMKEYWORD NUM_TEMPLATES
Specifies the template file.
*/
this->parseKeyword("TEMPLATE_FILE", this->templatefilename_list);
/* LMKEYWORD RMAX
Specifies the major radius of the region where core template applied
*/
this->parseKeyword("RMAX", this->rmax);
/* LMKEYWORD RMIN
Specifies the minor radius of the region where core template applied
*/
this->parseKeyword("RMIN", this->rmin);
/* LMKEYWORD RATIO_OF_PURE_CORE_TEMPLATE
Specifies the ratio of region where only pure core template is applied
with respect to ellipse defined by RMAX and RMIN
where
*/
this->parseKeyword("RATIO_OF_PURE_CORE_TEMPLATE", this->factor);
/* LMKEYWORD SHEAR_STRESS_YX
Specifies the shear sterss in X direction (\sigma_{yx}).
*/
this->parseKeyword("SHEAR_STRESS_YX", this->shear_stress_yx);
/* LMKEYWORD SHEAR_STRESS_YZ
Specifies the shear sterss in YZ direction (\sigma_{yz}).
*/
this->parseKeyword("SHEAR_STRESS_YZ", this->shear_stress_yz);
/* LMKEYWORD SHEAR_MODULUS
Specifies the shear modulus in Burgers vector direction.
mega-Pascal unit
*/
this->parseKeyword("SHEAR_MODULUS", this->mu);
/* LMKEYWORD DOMAIN_DD
Specifies the dd domain id
*/
this->parseKeyword("DOMAIN_DD", this->domdd);
this->addSubParsableObject(this->core_template);
}
/* -------------------------------------------------------------------------- */
template <typename RefAtom, typename IteratorAtoms, bool use_template>
void TwoDTemplate<RefAtom, IteratorAtoms, use_template>::declareParams() {
}
#ifdef LIBMULTISCALE_USE_LAMMPS
template class TwoDTemplate<RefLammps<2>, ContainerArray<RefLammps<2> >::iterator, true >;
template class TwoDTemplate<RefLammps<3>, ContainerArray<RefLammps<3> >::iterator, true >;
#endif
DECLARE_DD_ATOMIC_TEMPLATE(CADD_template);
__END_LIBMULTISCALE__
c4science · Help