Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81631454
TransferLibrary.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 7, 03:25
Size
71 KB
Mime Type
text/x-c
Expires
Mon, Sep 9, 03:25 (2 d)
Engine
blob
Format
Raw Data
Handle
20601050
Attached To
rLAMMPS lammps
TransferLibrary.cpp
View Options
// ATC headers
#include "TransferLibrary.h"
#include "ATC_Coupling.h"
#include "PrescribedDataManager.h"
#include "LinearSolver.h"
#include "PerAtomQuantityLibrary.h"
#include "KernelFunction.h"
#include "MoleculeSet.h"
//#include <typeinfo>
#include <set>
#include <sstream>
#include <utility>
#include <vector>
using std::set;
using std::map;
using std::string;
using std::stringstream;
using std::pair;
using std::vector;
namespace ATC {
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodalAtomVolume
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
NodalAtomVolume::NodalAtomVolume(ATC_Method * atc,
SPAR_MAN * shapeFunction) :
atc_(atc),
shapeFunction_(shapeFunction),
lammpsInterface_(atc->lammps_interface()),
feEngine_(atc->fe_engine()),
tol_(1.e-10)
{
shapeFunction_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void NodalAtomVolume::reset_quantity() const
{
// solve equation \sum_a N_Ia \sum_J N_Ja dV_J = \int_Omega N_I dV
// form left-hand side
int nNodes = shapeFunction_->nCols();
SPAR_MAT lhs(nNodes,nNodes);
atc_->compute_consistent_md_mass_matrix(shapeFunction_->quantity(),lhs);
// form right-hand side
_rhsMatrix_.resize(nNodes,nNodes);
feEngine_->compute_lumped_mass_matrix(_rhsMatrix_);
_rhs_.resize(nNodes);
_rhs_.copy(_rhsMatrix_.ptr(),_rhsMatrix_.size(),1);
// change entries for all entries if no atoms in shape function support
double totalVolume = _rhs_.sum();
double averageVolume = averaging_operation(totalVolume);
_scale_.resize(nNodes);
for (int i = 0; i < nNodes; i++) {
if ((abs(lhs(i,i)) > 0.))
_scale_(i) = 1.;
else
_scale_(i) = 0.;
}
lhs.row_scale(_scale_);
for (int i = 0; i < nNodes; i++) {
if (_scale_(i) < 0.5) {
lhs.set(i,i,1.);
_rhs_(i) = averageVolume;
}
}
lhs.compress();
// solve equation
LinearSolver solver(lhs, ATC::LinearSolver::ITERATIVE_SOLVE_SYMMETRIC, true);
solver.set_max_iterations(lhs.nRows());
solver.set_tolerance(tol_);
quantity_.reset(nNodes,1);
CLON_VEC tempQuantity(quantity_,CLONE_COL,0);
solver.solve(tempQuantity,_rhs_);
}
//--------------------------------------------------------
// averaging_operation
//--------------------------------------------------------
double NodalAtomVolume::averaging_operation(const double totalVolume) const
{
int nLocal[1] = {shapeFunction_->nRows()};
int nGlobal[1] = {0};
lammpsInterface_->int_allsum(nLocal,nGlobal,1);
return totalVolume/(double(nGlobal[0]));
}
//--------------------------------------------------------
// overloading operation to get number of rows
//--------------------------------------------------------
int NodalAtomVolume::nRows() const
{
return atc_->num_nodes();
}
//--------------------------------------------------------
// overloading operation to get number of columns
//--------------------------------------------------------
int NodalAtomVolume::nCols() const
{
return 1;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodalVolume
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// averaging_operation
//--------------------------------------------------------
double NodalVolume::averaging_operation(const double totalVolume) const
{
int nNodes = shapeFunction_->nCols();
return totalVolume/nNodes;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodalAtomVolumeElement
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
NodalAtomVolumeElement::NodalAtomVolumeElement(ATC_Method * atc,
SPAR_MAN * shapeFunction,
PerAtomQuantity<int> * atomElement) :
atc_(atc),
shapeFunction_(shapeFunction),
atomElement_(atomElement),
feEngine_(atc->fe_engine()),
tol_(1.e-10)
{
shapeFunction_->register_dependence(this);
if (!atomElement_) {
InterscaleManager & interscaleManager = atc_->interscale_manager();
atomElement_ = interscaleManager.per_atom_int_quantity("AtomElement");
}
atomElement_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void NodalAtomVolumeElement::reset_quantity() const
{
// Using analyses by G. Wagner and J. Templeton, weights ~ phi*M^{-1}*V
// where phi are the dimensionless shape/weighting functions,
// M is the "mass" matrix M_IJ,
// V is the vector of nodal and element volumes
//
//
// form atom-element shape functions and elemental volumes
const FE_Mesh * feMesh = feEngine_->fe_mesh();
int nElts = feMesh->num_elements();
int nNodes = shapeFunction_->nCols();
int nLocal = shapeFunction_->nRows();
// form "mass" matrix M_IJ (I,J = 0, 1, ..., nNodes+nElts-1)
int neSize = nNodes+nElts;
const INT_ARRAY & atomElement(atomElement_->quantity());
SPAR_MAT nodEltShpFcnMatrix(nLocal,neSize);
const SPAR_MAT & shapeFunction(shapeFunction_->quantity());
for(int a = 0; a < nLocal; a++) {
for(int I = 0; I < nNodes; I++) {
nodEltShpFcnMatrix.set(a,I,shapeFunction(a,I));
}
int thisCol = nNodes+atomElement(a,0);
nodEltShpFcnMatrix.set(a,thisCol,1);
}
SPAR_MAT neMassMatrix(neSize,neSize);
atc_->compute_consistent_md_mass_matrix(nodEltShpFcnMatrix,neMassMatrix);
// form vector of nodal and elemental volumes
_nodeVolumesMatrix_.resize(nNodes,nNodes);
feEngine_->compute_lumped_mass_matrix(_nodeVolumesMatrix_);
_nodeVolumes_.resize(nNodes);
_nodeVolumes_.copy(_nodeVolumesMatrix_.ptr(),_nodeVolumesMatrix_.size(),1);
DENS_VEC _nodeElementVolumes_(neSize);
for(int I = 0; I < nNodes; I++) {
_nodeElementVolumes_(I) = _nodeVolumes_(I);
}
double averageEltVolume = 0.0;
for(int E = 0; E < nElts; E++) {
double minx, maxx, miny, maxy, minz, maxz;
feMesh->element_coordinates(E,_nodalCoords_);
minx = _nodalCoords_(0,0); maxx = _nodalCoords_(0,0);
miny = _nodalCoords_(1,0); maxy = _nodalCoords_(1,0);
minz = _nodalCoords_(2,0); maxz = _nodalCoords_(2,0);
for (int j = 1; j < _nodalCoords_.nCols(); ++j) {
if (_nodalCoords_(0,j)<minx) minx = _nodalCoords_(0,j);
if (_nodalCoords_(0,j)>maxx) maxx = _nodalCoords_(0,j);
if (_nodalCoords_(1,j)<miny) miny = _nodalCoords_(1,j);
if (_nodalCoords_(1,j)>maxy) maxy = _nodalCoords_(1,j);
if (_nodalCoords_(2,j)<minz) minz = _nodalCoords_(2,j);
if (_nodalCoords_(2,j)>maxz) maxz = _nodalCoords_(2,j);
}
_nodeElementVolumes_(nNodes+E) = (maxx-minx)*(maxy-miny)*(maxz-minz);
averageEltVolume += (maxx-minx)*(maxy-miny)*(maxz-minz);
}
averageEltVolume /= nElts;
// correct entries of mass matrix if no atoms in shape function support
double totalNodalVolume = _nodeVolumes_.sum();
double averageNodalVolume = totalNodalVolume/nNodes;
_scale_.resize(neSize);
for (int i = 0; i < neSize; i++) {
if ((abs(neMassMatrix(i,i)) > 0.)) {
_scale_(i) = 1.;
} else {
printf("No atoms are in support of node/element %i\n",i);
_scale_(i) = 0.;
}
}
neMassMatrix.row_scale(_scale_);
for (int i = 0; i < neSize; i++) {
if (_scale_(i) < 0.5) {
neMassMatrix.set(i,i,1.);
if (i < nNodes) {
_nodeElementVolumes_(i) = averageNodalVolume;
} else {
_nodeElementVolumes_(i) = averageEltVolume;
}
}
}
neMassMatrix.compress();
// solve equation
LinearSolver solver(neMassMatrix, ATC::LinearSolver::ITERATIVE_SOLVE_SYMMETRIC, true);
solver.set_max_iterations(neMassMatrix.nRows());
double myTol = 1.e-10;
solver.set_tolerance(myTol);
quantity_.resize(neSize,0);
CLON_VEC tempQuantity(quantity_,CLONE_COL,0);
solver.solve(tempQuantity,_nodeElementVolumes_);
}
//--------------------------------------------------------
// overloading operation to get number of rows
//--------------------------------------------------------
int NodalAtomVolumeElement::nRows() const
{
return atc_->num_nodes();
}
//--------------------------------------------------------
// overloading operation to get number of columns
//--------------------------------------------------------
int NodalAtomVolumeElement::nCols() const
{
return 1;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class AtomTypeElement
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
AtomTypeElement::AtomTypeElement(ATC_Coupling * atc,
PerAtomQuantity<int> * atomElement) :
atomElement_(atomElement),
nElts_((atc->fe_engine())->num_elements())
{
if (!atomElement_) {
atomElement_ = (atc->interscale_manager()).per_atom_int_quantity("AtomElement");
}
atomElement_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void AtomTypeElement::reset_quantity() const
{
// determine which elements contain internal atoms
quantity_.resize(nElts_,1);
_quantityLocal_.resize(nElts_,1);
_quantityLocal_ = 0;
const INT_ARRAY & atomElement(atomElement_->quantity());
for (int i = 0; i < atomElement_->nRows(); ++i) {
_quantityLocal_(atomElement(i,0),0) = 1;
}
// swap contributions
lammpsInterface_->logical_or(_quantityLocal_.ptr(),
quantity_.ptr(),nElts_);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ElementMask
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ElementMask::ElementMask(ATC_Coupling * atc,
MatrixDependencyManager<DenseMatrix, int> * hasInternal,
MatrixDependencyManager<DenseMatrix, int> * hasGhost) :
hasInternal_(hasInternal),
hasGhost_(hasGhost),
feEngine_(atc->fe_engine())
{
if (!hasInternal_) {
hasInternal_ = (atc->interscale_manager()).dense_matrix_int("ElementHasInternal");
}
if (!hasGhost_) {
hasGhost_ = (atc->interscale_manager()).dense_matrix_int("ElementHasGhost");
}
hasInternal_->register_dependence(this);
if (hasGhost_) hasGhost_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void ElementMask::reset_quantity() const
{
const INT_ARRAY & hasInternal(hasInternal_->quantity());
int nElts = hasInternal.size();
quantity_.resize(nElts,1);
if (hasGhost_) {
const INT_ARRAY & hasGhost(hasGhost_->quantity());
for (int i = 0; i < nElts; ++i) {
quantity_(i,0) = !hasInternal(i,0) || hasGhost(i,0);
}
}
else {
for (int i = 0; i < nElts; ++i) {
quantity_(i,0) = !hasInternal(i,0);
}
}
const set<int> & nullElements = feEngine_->null_elements();
set<int>::const_iterator iset;
for (iset = nullElements.begin(); iset != nullElements.end(); iset++) {
int ielem = *iset;
quantity_(ielem,0) = false;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class AtomElementMask
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
AtomElementMask::AtomElementMask(ATC_Coupling * atc,
MatrixDependencyManager<DenseMatrix, int> * hasAtoms) :
hasAtoms_(hasAtoms),
feEngine_(atc->fe_engine())
{
if (!hasAtoms_) {
hasAtoms_ = (atc->interscale_manager()).dense_matrix_int("ElementHasInternal");
}
if (!hasAtoms_) {
throw ATC_Error("AtomElementMask::AtomElementMask - no element has atoms transfer provided");
}
hasAtoms_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void AtomElementMask::reset_quantity() const
{
const INT_ARRAY & hasAtoms(hasAtoms_->quantity());
int nElts = hasAtoms.size();
quantity_.resize(nElts,1);
for (int i = 0; i < nElts; ++i) {
quantity_(i,0) = hasAtoms(i,0);
}
// this seems to cause problems because many materials end up being null
const set<int> & nullElements = feEngine_->null_elements();
set<int>::const_iterator iset;
for (iset = nullElements.begin(); iset != nullElements.end(); iset++) {
int ielem = *iset;
quantity_(ielem,0) = false;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ElementMaskNodeSet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ElementMaskNodeSet::ElementMaskNodeSet(ATC_Coupling * atc,
SetDependencyManager<int> * nodeSet) :
nodeSet_(nodeSet),
feMesh_((atc->fe_engine())->fe_mesh())
{
nodeSet_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void ElementMaskNodeSet::reset_quantity() const
{
quantity_.resize(feMesh_->num_elements(),1);
quantity_ = false;
// get the maximal element set corresponding to those nodes
set<int> elementSet;
feMesh_->nodeset_to_maximal_elementset(nodeSet_->quantity(),elementSet);
set<int>::const_iterator iset;
for (iset = elementSet.begin(); iset != elementSet.end(); iset++) {
quantity_(*iset,0) = true;
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodalGeometryType
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
NodalGeometryType::NodalGeometryType(ATC_Coupling * atc,
MatrixDependencyManager<DenseMatrix, int> * hasInternal,
MatrixDependencyManager<DenseMatrix, int> * hasGhost) :
hasInternal_(hasInternal),
hasGhost_(hasGhost),
feEngine_(atc->fe_engine()),
nNodes_(atc->num_nodes()),
nElts_((atc->fe_engine())->num_elements())
{
if (!hasInternal_) {
hasInternal_ = (atc->interscale_manager()).dense_matrix_int("ElementHasInternal");
}
if (!hasGhost_) {
hasGhost_ = (atc->interscale_manager()).dense_matrix_int("ElementHasGhost");
}
hasInternal_->register_dependence(this);
if (hasGhost_) hasGhost_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void NodalGeometryType::reset_quantity() const
{
const INT_ARRAY & hasInternal(hasInternal_->quantity());
_nodesInternal_.resize(nNodes_);
_nodesInternal_ = 0;
_nodesGhost_.reset(nNodes_);
_nodesGhost_ = 0;
Array<int> nodes;
vector<int> myElems = feEngine_->fe_mesh()->owned_elts();
if (hasGhost_) {
const INT_ARRAY & hasGhost(hasGhost_->quantity()) ;
// iterate through all elements owned by this processor
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
if (hasInternal(ielem,0) || hasGhost(ielem,0)) {
feEngine_->element_connectivity(ielem,nodes);
for (int j = 0; j < nodes.size(); j++) {
if (hasInternal(ielem,0)) {
_nodesInternal_(nodes(j)) = 1;
}
if (hasGhost(ielem,0)) {
_nodesGhost_(nodes(j)) = 1;
}
}
}
}
// sum up partial result arrays
lammpsInterface_->logical_or(MPI_IN_PLACE, _nodesInternal_.ptr(), _nodesInternal_.size());
lammpsInterface_->logical_or(MPI_IN_PLACE, _nodesGhost_.ptr(), _nodesGhost_.size());
}
else {
// iterate through all elements owned by this processor
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
if (hasInternal(ielem,0)) {
feEngine_->element_connectivity(ielem,nodes);
for (int j = 0; j < nodes.size(); j++) {
_nodesInternal_(nodes(j)) = 1;
}
}
}
// sum up partial result arrays
lammpsInterface_->logical_or(MPI_IN_PLACE, _nodesInternal_.ptr(), _nodesInternal_.size());
}
quantity_.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (_nodesInternal_(i) && _nodesGhost_(i)) {
quantity_(i,0) = BOUNDARY;
}
else if (_nodesInternal_(i)) {
quantity_(i,0) = MD_ONLY;
}
else {
quantity_(i,0) = FE_ONLY;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodalGeometryTypeElementSet
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
NodalGeometryTypeElementSet::NodalGeometryTypeElementSet(ATC_Coupling * atc,
MatrixDependencyManager<DenseMatrix, int> * hasInternal) :
hasInternal_(hasInternal),
feEngine_(atc->fe_engine()),
nNodes_(atc->num_nodes()),
nElts_((atc->fe_engine())->num_elements())
{
if (!hasInternal_) {
hasInternal_ = (atc->interscale_manager()).dense_matrix_int("ElementHasInternal");
}
if (!hasInternal_) {
throw ATC_Error("NodalGeometryTypeElementSet: No ElementHasInternal object provided or exists");
}
hasInternal_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void NodalGeometryTypeElementSet::reset_quantity() const
{
const INT_ARRAY & hasInternal(hasInternal_->quantity());
_nodesInternal_.resize(nNodes_);
_nodesInternal_ = 0;
_nodesGhost_.reset(nNodes_);
_nodesGhost_ = 0;
Array<int> nodes;
vector<int> myElems = feEngine_->fe_mesh()->owned_elts();
// iterate through all elements owned by this processor
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
if (hasInternal(ielem,0)) {
feEngine_->element_connectivity(ielem,nodes);
for (int j = 0; j < nodes.size(); j++) {
_nodesInternal_(nodes(j)) = 1;
}
}
else {
feEngine_->element_connectivity(ielem,nodes);
for (int j = 0; j < nodes.size(); j++) {
_nodesGhost_(nodes(j)) = 1;
}
}
}
// sum up partial result arrays
lammpsInterface_->logical_or(MPI_IN_PLACE, _nodesInternal_.ptr(), _nodesInternal_.size());
lammpsInterface_->logical_or(MPI_IN_PLACE, _nodesGhost_.ptr(), _nodesGhost_.size());
quantity_.resize(nNodes_,1);
for (int i = 0; i < nNodes_; i++) {
if (_nodesInternal_(i) && _nodesGhost_(i)) {
quantity_(i,0) = BOUNDARY;
}
else if (_nodesInternal_(i)) {
quantity_(i,0) = MD_ONLY;
}
else {
quantity_(i,0) = FE_ONLY;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NodeToSubset
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
NodeToSubset::NodeToSubset(ATC_Method * atc,
SetDependencyManager<int> * subsetNodes) :
LargeToSmallMap(),
atc_(atc),
subsetNodes_(subsetNodes)
{
subsetNodes_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void NodeToSubset::reset_quantity() const
{
int nNodes = atc_->num_nodes();
const set<int> & subsetNodes(subsetNodes_->quantity());
quantity_.resize(nNodes,1);
size_ = 0;
for (int i = 0; i < nNodes; i++) {
if (subsetNodes.find(i) != subsetNodes.end()) {
quantity_(i,0) = size_;
size_++;
}
else {
quantity_(i,0) = -1;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class SubsetToNode
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
SubsetToNode::SubsetToNode(NodeToSubset * nodeToSubset) :
nodeToSubset_(nodeToSubset)
{
nodeToSubset_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void SubsetToNode::reset_quantity() const
{
const INT_ARRAY & nodeToSubset(nodeToSubset_->quantity());
int nNodes = nodeToSubset.nRows();
int count = 0;
quantity_.resize(nodeToSubset_->size(),1);
for (int i = 0; i < nNodes; i++) {
if (nodeToSubset(i,0) > -1) {
quantity_(count,0) = i;
count++;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class ReducedSparseMatrix
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
ReducedSparseMatrix::ReducedSparseMatrix(ATC_Method * atc,
SPAR_MAN * source,
LargeToSmallAtomMap * map) :
SparseMatrixTransfer<double>(),
source_(source),
map_(map)
{
source_->register_dependence(this);
map_->register_dependence(this);
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
ReducedSparseMatrix::~ReducedSparseMatrix()
{
source_->remove_dependence(this);
map_->remove_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void ReducedSparseMatrix::reset_quantity() const
{
const SPAR_MAT & source(source_->quantity());
const INT_ARRAY & map(map_->quantity());
quantity_.reset(source.nRows(),source.nCols());
for (int i = 0; i < source.nRows(); i++) {
int idx = map(i,0);
if (idx > -1) {
source.row(i,_row_,_index_);
for (int j = 0; j < _row_.size(); j++) {
quantity_.set(i,_index_(j),_row_(j));
}
}
}
quantity_.compress();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class RowMappedSparseMatrix
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void RowMappedSparseMatrix::reset_quantity() const
{
const SPAR_MAT & source(source_->quantity());
const INT_ARRAY & map(map_->quantity());
quantity_.reset(map_->size(),source.nCols());
for (int i = 0; i < source.nRows(); i++) {
int idx = map(i,0);
if (idx > -1) {
source.row(i,_row_,_index_);
for (int j = 0; j < _row_.size(); j++) {
quantity_.set(idx,_index_(j),_row_(j));
}
}
}
quantity_.compress();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class RowMappedSparseMatrixVector
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
RowMappedSparseMatrixVector::RowMappedSparseMatrixVector(ATC_Method * atc,
VectorDependencyManager<SPAR_MAT * > * source,
LargeToSmallAtomMap * map) :
VectorTransfer<SPAR_MAT * >(),
source_(source),
map_(map)
{
source_->register_dependence(this);
map_->register_dependence(this);
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
RowMappedSparseMatrixVector::~RowMappedSparseMatrixVector()
{
source_->remove_dependence(this);
map_->remove_dependence(this);
for (unsigned i = 0; i < quantity_.size(); ++i) {
if (quantity_[i]) delete quantity_[i];
}
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void RowMappedSparseMatrixVector::reset_quantity() const
{
const vector<SPAR_MAT * > & source(source_->quantity());
const INT_ARRAY & map(map_->quantity());
for (unsigned i = 0; i < quantity_.size(); ++i) {
if (quantity_[i]) delete quantity_[i];
}
quantity_.resize(source.size(),NULL);
for (unsigned i = 0; i < source.size(); i++) {
quantity_[i] = new SPAR_MAT(map_->size(),source[i]->nCols());
}
for (unsigned i = 0; i < source.size(); i++) {
for (int j = 0; j < source[i]->nRows(); j++) {
int idx = map(j,0);
if (idx > -1) {
source[i]->row(j,_row_,_index_);
for (int k = 0; k < _row_.size(); k++) {
quantity_[i]->set(idx,_index_(k),_row_(k));
}
}
}
}
for (unsigned i = 0; i < source.size(); i++) {
quantity_[i]->compress();
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class MappedDiagonalMatrix
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
MappedDiagonalMatrix::MappedDiagonalMatrix(ATC_Method * atc,
DIAG_MAN * source,
LargeToSmallAtomMap * map) :
DiagonalMatrixTransfer<double>(),
source_(source),
map_(map)
{
source_->register_dependence(this);
map_->register_dependence(this);
}
//--------------------------------------------------------
// Destructor
//--------------------------------------------------------
MappedDiagonalMatrix::~MappedDiagonalMatrix()
{
source_->remove_dependence(this);
map_->remove_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void MappedDiagonalMatrix::reset_quantity() const
{
const DIAG_MAT & source(source_->quantity());
const INT_ARRAY & map(map_->quantity());
quantity_.resize(map_->size(),map_->size());
for (int i = 0; i < source.nRows(); i++) {
int idx = map(i,0);
if (idx > -1) {
quantity_(idx,idx) = source(i,i);
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class MappedQuantity
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
MappedQuantity::MappedQuantity(ATC_Method * atc,
DENS_MAN * source,
LargeToSmallMap * map) :
DenseMatrixTransfer<double>(),
source_(source),
map_(map)
{
source_->register_dependence(this);
map_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void MappedQuantity::reset_quantity() const
{
const DENS_MAT & source(source_->quantity());
const INT_ARRAY & map(map_->quantity());
quantity_.resize(map_->size(),source.nCols());
for (int i = 0; i < source.nRows(); i++) {
int idx = map(i,0);
if (idx > -1) {
for (int j = 0; j < source.nCols(); j++) {
quantity_(idx,j) = source(i,j);
}
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class RegulatedNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
RegulatedNodes::RegulatedNodes(ATC_Coupling * atc,
FieldName fieldName,
MatrixDependencyManager<DenseMatrix, int> * nodalGeometryType) :
SetTransfer<int>(),
nodalGeometryType_(nodalGeometryType),
atc_(atc),
feEngine_(atc->fe_engine()),
prescribedDataManager_(atc->prescribed_data_manager())
{
if (!nodalGeometryType_) {
nodalGeometryType_ = (atc_->interscale_manager()).dense_matrix_int("NodalGeometryType");
}
if (nodalGeometryType_) {
nodalGeometryType_->register_dependence(this);
}
else {
throw ATC_Error("TransferLibrary::RegulatedNodes - No Nodal Geometry Type provided");
}
const map<FieldName,int> & atcFieldSizes(atc_->field_sizes());
if (fieldName == NUM_TOTAL_FIELDS) {
fieldSizes_ = atcFieldSizes;
}
else {
map<FieldName,int>::const_iterator fs_iter = atcFieldSizes.find(fieldName);
fieldSizes_.insert(pair<FieldName,int>(fieldName,fs_iter->second));
}
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void RegulatedNodes::reset_quantity() const
{
const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
quantity_.clear();
for (int i = 0; i < nodeType.size(); ++i) {
if (nodeType(i,0) != FE_ONLY) {
quantity_.insert(i);
}
}
}
//--------------------------------------------------------
// insert_boundary_nodes
//--------------------------------------------------------
void RegulatedNodes::insert_boundary_nodes() const
{
const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
for (int i = 0; i < nodeType.size(); ++i) {
if (nodeType(i,0) == BOUNDARY) {
quantity_.insert(i);
}
}
}
//--------------------------------------------------------
// insert_boundary_faces
//--------------------------------------------------------
void RegulatedNodes::insert_boundary_faces() const
{
const set<string> & boundaryFaceNames(atc_->boundary_face_names());
set<string>::const_iterator iter;
set<int>::const_iterator nodeIter;
for (iter = boundaryFaceNames.begin(); iter != boundaryFaceNames.end(); iter++) {
set<int> nodeSet;
feEngine_->fe_mesh()->faceset_to_nodeset(*iter,nodeSet);
for (nodeIter = nodeSet.begin(); nodeIter != nodeSet.end(); nodeIter++) {
quantity_.insert(*nodeIter);
}
}
}
//--------------------------------------------------------
// insert_fixed_nodes
//--------------------------------------------------------
void RegulatedNodes::insert_fixed_nodes() const
{
const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
map<FieldName,int>::const_iterator fs_iter;
for (int i = 0; i < nodeType.size(); ++i) {
bool isFixed = false;
for (fs_iter = fieldSizes_.begin(); fs_iter != fieldSizes_.end(); fs_iter++) {
for (int j = 0; j < fs_iter->second; j++) {
isFixed = prescribedDataManager_->is_fixed(i,fs_iter->first,j);
if (isFixed) break;
}
}
if (isFixed && ((nodeType(i,0)==MD_ONLY) || (nodeType(i,0)==BOUNDARY))) {
quantity_.insert(i);
}
}
}
//--------------------------------------------------------
// insert_face_fluxes
//--------------------------------------------------------
void RegulatedNodes::insert_face_fluxes() const
{
const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
set<int>::const_iterator inode;
map<FieldName,int>::const_iterator fs_iter;
for (fs_iter = fieldSizes_.begin(); fs_iter != fieldSizes_.end(); fs_iter++) {
for (int j = 0; j < fs_iter->second; j++) {
set<int> faceFluxNodes = prescribedDataManager_->flux_face_nodes(fs_iter->first,j);
for (inode = faceFluxNodes.begin(); inode != faceFluxNodes.end(); inode++) {
if (((nodeType(*inode,0)==MD_ONLY) || (nodeType(*inode,0)==BOUNDARY))) {
quantity_.insert(*inode);
}
}
}
}
}
//--------------------------------------------------------
// insert_element_fluxes
//--------------------------------------------------------
void RegulatedNodes::insert_element_fluxes() const
{
const INT_ARRAY & nodeType(nodalGeometryType_->quantity());
set<int>::const_iterator inode;
map<FieldName,int>::const_iterator fs_iter;
for (fs_iter = fieldSizes_.begin(); fs_iter != fieldSizes_.end(); fs_iter++) {
for (int j = 0; j < fs_iter->second; j++) {
set<int> elementFluxNodes = prescribedDataManager_->flux_element_nodes(fs_iter->first,j);
for (inode = elementFluxNodes.begin(); inode != elementFluxNodes.end(); inode++) {
if (((nodeType(*inode,0)==MD_ONLY) || (nodeType(*inode,0)==BOUNDARY))) {
quantity_.insert(*inode);
}
}
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class BoundaryNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void BoundaryNodes::reset_quantity() const
{
quantity_.clear();
// a) they are a boundary node
RegulatedNodes::insert_boundary_nodes();
// b) they are in a specified boundary faceset
RegulatedNodes::insert_boundary_faces();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class FluxNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void FluxNodes::reset_quantity() const
{
quantity_.clear();
// a) they have a fixed face flux
RegulatedNodes::insert_face_fluxes();
// b) they have a fixed element flux
RegulatedNodes::insert_element_fluxes();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class FluxBoundaryNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void FluxBoundaryNodes::reset_quantity() const
{
FluxNodes::reset_quantity();
// a) they are a boundary node
RegulatedNodes::insert_boundary_nodes();
// b) they are in a specified boundary faceset
RegulatedNodes::insert_boundary_faces();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class AllRegulatedNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void AllRegulatedNodes::reset_quantity() const
{
FluxBoundaryNodes::reset_quantity();
// a) they are a fixed node
RegulatedNodes::insert_fixed_nodes();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class FixedNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void FixedNodes::reset_quantity() const
{
quantity_.clear();
// a) they are a fixed node
RegulatedNodes::insert_fixed_nodes();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class FixedBoundaryNodes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void FixedBoundaryNodes::reset_quantity() const
{
FixedNodes::reset_quantity();
// a) they are a boundary node
RegulatedNodes::insert_boundary_nodes();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class DenseMatrixQuotient
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
DenseMatrixQuotient::DenseMatrixQuotient(DENS_MAN* matrixNumerator,
DENS_MAN* matrixDenominator):
matrixNumerator_(matrixNumerator),
matrixDenominator_(matrixDenominator)
{
matrixNumerator_->register_dependence(this);
matrixDenominator_->register_dependence(this);
}
//--------------------------------------------------------
// Reset_quantity
//--------------------------------------------------------
void DenseMatrixQuotient::reset_quantity() const
{
quantity_ = matrixNumerator_->quantity();
quantity_.divide_zero_safe(matrixDenominator_->quantity());
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class DenseMatrixSum
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
DenseMatrixSum::DenseMatrixSum(DENS_MAN* matrixOne,
DENS_MAN* matrixTwo):
matrixOne_(matrixOne), matrixTwo_(matrixTwo)
{
matrixOne_->register_dependence(this);
matrixTwo_->register_dependence(this);
}
//--------------------------------------------------------
// Reset_quantity
//--------------------------------------------------------
void DenseMatrixSum::reset_quantity() const
{
quantity_ = matrixOne_->quantity();
quantity_ += (matrixTwo_->quantity());
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class DenseMatrixDelta
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
DenseMatrixDelta::DenseMatrixDelta(DENS_MAN* matrix,
DENS_MAT* reference):
matrix_(matrix), reference_(reference)
{
matrix_->register_dependence(this);
}
//--------------------------------------------------------
// Reset_quantity
//--------------------------------------------------------
void DenseMatrixDelta::reset_quantity() const
{
quantity_ = matrix_->quantity();
quantity_ -= *reference_;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class PointToElementMap
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
PointToElementMap::PointToElementMap(ATC_Method * atc,
MatrixDependencyManager<DenseMatrix, double> * pointPositions) :
DenseMatrixTransfer<int>(),
pointPositions_(pointPositions),
feMesh_((atc->fe_engine())->fe_mesh())
{
pointPositions_->register_dependence(this);
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
PointToElementMap::~PointToElementMap()
{
pointPositions_->remove_dependence(this);
}
//--------------------------------------------------------
// reset
//--------------------------------------------------------
void PointToElementMap::reset_quantity() const
{
const DENS_MAT & position(pointPositions_->quantity());
int nsd = position.nCols();
int nRows = position.nRows();
quantity_.resize(nRows,nsd);
DENS_VEC coords(nsd);
for (int i = 0; i < nRows; i++) {
for (int j = 0; j < nsd; j++) {
coords(j) = position(i,j);
}
quantity_(i,0) = feMesh_->map_to_element(coords);
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class Interpolant
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// constructor
//--------------------------------------------------------
Interpolant::Interpolant(ATC_Method * atc,
MatrixDependencyManager<DenseMatrix, int>* pointToElementMap,
DENS_MAN* pointPositions) :
SparseMatrixTransfer<double>(),
pointToElementMap_(pointToElementMap),
pointPositions_(pointPositions),
feEngine_(atc->fe_engine())
{
pointToElementMap_->register_dependence(this);
pointPositions_->register_dependence(this);
}
//--------------------------------------------------------
// reset
//--------------------------------------------------------
void Interpolant::reset_quantity() const
{
const DENS_MAT & positions(pointPositions_->quantity());
const INT_ARRAY & elements(pointToElementMap_->quantity());
if (positions.nRows() > 0) {
feEngine_->evaluate_shape_functions(positions,
elements,
this->quantity_);
}
else {
quantity_.resize(0,feEngine_->num_nodes());
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class InterpolantGradient
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// constructor
//--------------------------------------------------------
InterpolantGradient::InterpolantGradient(ATC_Method * atc,
MatrixDependencyManager<DenseMatrix, int>* pointToElementMap,
DENS_MAN* pointPositions) :
VectorTransfer<SPAR_MAT * >(),
pointToElementMap_(pointToElementMap),
pointPositions_(pointPositions),
feEngine_(atc->fe_engine())
{
pointToElementMap_->register_dependence(this);
pointPositions_->register_dependence(this);
quantity_.resize(atc->nsd(),NULL);
for (int i = 0; i < atc->nsd(); ++i) {
quantity_[i] = new SPAR_MAT();
}
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
InterpolantGradient::~InterpolantGradient() {
pointToElementMap_->remove_dependence(this);
pointPositions_->remove_dependence(this);
for (unsigned i = 0; i < quantity_.size(); ++i) {
if (quantity_[i]) delete quantity_[i];
}
}
//--------------------------------------------------------
// reset_quantity()
//--------------------------------------------------------
void InterpolantGradient::reset_quantity() const
{
const DENS_MAT & positions(pointPositions_->quantity());
const INT_ARRAY & elements(pointToElementMap_->quantity());
if (positions.nRows() > 0) {
feEngine_->evaluate_shape_function_derivatives(positions,
elements,
this->quantity_);
}
else {
for (unsigned i = 0; i < quantity_.size(); ++i) {
(this->quantity_)[i]->resize(0,feEngine_->num_nodes());
}
}
}
//--------------------------------------------------------
// Class PerAtomShapeFunctionGradient
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// constructor
//--------------------------------------------------------
PerAtomShapeFunctionGradient::PerAtomShapeFunctionGradient(ATC_Method * atc,
MatrixDependencyManager<DenseMatrix, int>* atomToElementMap,
DENS_MAN* atomPositions,
const string & tag,
AtomType atomType) :
VectorTransfer<SPAR_MAT * >(),
atomToElementMap_(atomToElementMap),
atomPositions_(atomPositions),
feEngine_(atc->fe_engine())
{
InterscaleManager & interscaleManager(atc->interscale_manager());
if (!atomToElementMap_) {
atomToElementMap_ = interscaleManager.per_atom_int_quantity("AtomElement");
}
if (!atomPositions_) {
atomPositions_ = interscaleManager.per_atom_quantity("AtomicCoarseGrainingPositions");
}
atomToElementMap_->register_dependence(this);
atomPositions_->register_dependence(this);
// storage container
matrices_.resize(atc->nsd(),NULL);
for (int i = 0; i < atc->nsd(); ++i) {
matrices_[i] = new AtcAtomSparseMatrix<double>(atc,feEngine_->num_nodes(),
feEngine_->num_nodes_per_element(),
atomType);
stringstream myint;
myint << i;
interscaleManager.add_per_atom_sparse_matrix(matrices_[i],
tag+myint.str());
matrices_[i]->register_dependence(this);
}
quantity_.resize(atc->nsd(),NULL);
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
PerAtomShapeFunctionGradient::~PerAtomShapeFunctionGradient() {
atomToElementMap_->remove_dependence(this);
atomPositions_->remove_dependence(this);
for (unsigned i = 0; i < matrices_.size(); ++i) {
matrices_[i]->remove_dependence(this);
}
}
//--------------------------------------------------------
// reset_quantity()
//--------------------------------------------------------
void PerAtomShapeFunctionGradient::reset_quantity() const
{
const DENS_MAT & positions(atomPositions_->quantity());
const INT_ARRAY & elements(atomToElementMap_->quantity());
for (unsigned i = 0; i < quantity_.size(); ++i) {
(this->quantity_)[i] = & matrices_[i]->set_quantity();
}
if (positions.nRows() > 0) {
feEngine_->evaluate_shape_function_derivatives(positions,
elements,
this->quantity_);
}
else {
for (unsigned i = 0; i < quantity_.size(); ++i) {
(this->quantity_)[i]->resize(0,feEngine_->num_nodes());
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class InterpolantSmallMolecule
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// constructor
//--------------------------------------------------------
InterpolantSmallMolecule::InterpolantSmallMolecule(ATC_Method * atc,
MatrixDependencyManager<DenseMatrix, int>* moleculeToElementMap,
DENS_MAN* moleculePositions,
MoleculeSet * moleculeSet) :
Interpolant(atc,moleculeToElementMap,moleculePositions),
moleculeSet_(moleculeSet)
{
moleculeSet_->register_dependence(this);
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
InterpolantSmallMolecule::~InterpolantSmallMolecule()
{
moleculeSet_->remove_dependence(this);
}
//--------------------------------------------------------
// reset
//--------------------------------------------------------
void InterpolantSmallMolecule::reset_quantity() const
{
Interpolant::reset_quantity();
// scale rows by fraction of molecules on this proc
_fractions_.resize((this->quantity_).nRows());
for (int i = 0; i < moleculeSet_->local_molecule_count(); i++)
_fractions_(i) = moleculeSet_->local_fraction(i);
(this->quantity_).row_scale(_fractions_);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyKernelAccumulation
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyKernelAccumulation::OnTheFlyKernelAccumulation(ATC_Method * atc,
PerAtomQuantity<double> * source,
KernelFunction* kernelFunction,
DENS_MAN* atomCoarseGrainingPositions):
atc_(atc),
source_(source),
kernelFunction_(kernelFunction),
atomCoarseGrainingPositions_(atomCoarseGrainingPositions),
feMesh_((atc_->fe_engine())->fe_mesh())
{
source_->register_dependence(this);
atomCoarseGrainingPositions_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyKernelAccumulation::reset_quantity() const
{
const DENS_MAT & positions(atomCoarseGrainingPositions_->quantity());
const DENS_MAT & source(source_->quantity());
int nNodes = feMesh_->num_nodes_unique();
quantity_.resize(nNodes,source.nCols());
_quantityLocal_.reset(nNodes,source.nCols());
if (source.nRows()>0) {
DENS_VEC xI(positions.nCols()),xa(positions.nCols()),xaI(positions.nCols());
double val;
for (int i = 0; i < nNodes; i++) {
xI = feMesh_->nodal_coordinates(i);
for (int j = 0; j < positions.nRows(); j++) {
for (int k = 0; k < positions.nCols(); ++k) {
xa(k) = positions(j,k);
}
xaI = xa - xI;
atc_->lammps_interface()->periodicity_correction(xaI.ptr());
val = kernelFunction_->value(xaI);
if (val > 0) {
for (int k=0; k < source.nCols(); k++) {
_quantityLocal_(i,k) += val*source(j,k);
}
}
}
}
}
// accumulate across processors
int count = quantity_.nRows()*quantity_.nCols();
lammpsInterface_->allsum(_quantityLocal_.ptr(),quantity_.ptr(),count);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyKernelAccumulationNormalized
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyKernelAccumulationNormalized::OnTheFlyKernelAccumulationNormalized(ATC_Method * atc,
PerAtomQuantity<double> * source,
KernelFunction* kernelFunction,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN * weights):
OnTheFlyKernelAccumulation(atc,source,kernelFunction,atomCoarseGrainingPositions),
weights_(weights)
{
weights_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyKernelAccumulationNormalized::reset_quantity() const
{
OnTheFlyKernelAccumulation::reset_quantity();
quantity_ *= weights_->quantity();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyKernelAccumulationNormalizedReferenced
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyKernelAccumulationNormalizedReferenced::OnTheFlyKernelAccumulationNormalizedReferenced(ATC_Method * atc,
PerAtomQuantity<double> * source,
KernelFunction* kernelFunction,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN* weights,
DENS_MAN * reference):
OnTheFlyKernelAccumulationNormalized(atc,source,kernelFunction,atomCoarseGrainingPositions,weights),
reference_(reference)
{
reference_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyKernelAccumulationNormalizedReferenced::reset_quantity() const
{
OnTheFlyKernelAccumulationNormalized::reset_quantity();
quantity_ -= reference_->quantity();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyKernelAccumulationNormalizedScaled
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyKernelAccumulationNormalizedScaled::OnTheFlyKernelAccumulationNormalizedScaled(ATC_Method * atc,
PerAtomQuantity<double> * source,
KernelFunction* kernelFunction,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN* weights,
const double scale):
OnTheFlyKernelAccumulationNormalized(atc,source,kernelFunction,atomCoarseGrainingPositions,weights),
scale_(scale)
{
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyKernelAccumulationNormalizedScaled::reset_quantity() const
{
OnTheFlyKernelAccumulationNormalized::reset_quantity();
quantity_ *= scale_;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class AccumulantWeights
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
AccumulantWeights::AccumulantWeights(SPAR_MAN* accumulant):
DiagonalMatrixTransfer<double>(),
accumulant_(accumulant)
{
accumulant_->register_dependence(this);
}
//--------------------------------------------------------
// Reset_quantity
//--------------------------------------------------------
void AccumulantWeights::reset_quantity() const
{
const SPAR_MAT accumulant(accumulant_->quantity());
int nNodes = accumulant.nCols();
// get summation of atoms per node
_localWeights_.reset(nNodes); _weights_.resize(nNodes);
if (accumulant.nRows()>0) {
_localWeights_ = (accumulant_->quantity()).col_sum();
}
lammpsInterface_->allsum(_localWeights_.ptr(),_weights_.ptr(),nNodes);
// assign weights
quantity_.resize(nNodes,nNodes);
for (int i = 0; i < nNodes; i++) {
if (_weights_(i) > 0.) {
quantity_(i,i) = 1./_weights_(i);
}
else {
quantity_(i,i) = 0.;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyKernelWeights
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyKernelWeights::OnTheFlyKernelWeights(DENS_MAN* weights):
DiagonalMatrixTransfer<double>(),
weights_(weights)
{
weights_->register_dependence(this);
}
//--------------------------------------------------------
// Reset_quantity
//--------------------------------------------------------
void OnTheFlyKernelWeights::reset_quantity() const
{
const DENS_MAT & weights(weights_->quantity());
int nNodes = weights.nRows();
// assign weights
quantity_.resize(nNodes,nNodes);
for (int i = 0; i < nNodes; i++) {
if (weights(i,0) > 0.) {
quantity_(i,i) = 1./weights(i,0);
}
else {
quantity_(i,i) = 0.;
}
}
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class KernelInverseVolumes
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
KernelInverseVolumes::KernelInverseVolumes(ATC_Method * atc,
KernelFunction* kernelFunction):
kernelFunction_(kernelFunction),
feMesh_((atc->fe_engine())->fe_mesh())
{
// do nothing
}
//--------------------------------------------------------
// multiplication by transpose
//--------------------------------------------------------
void KernelInverseVolumes::reset_quantity() const
{
int nNodes = feMesh_->num_nodes_unique();
quantity_.resize(nNodes,nNodes);
quantity_ = kernelFunction_->inv_vol();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyMeshAccumulation
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyMeshAccumulation::OnTheFlyMeshAccumulation(ATC_Method * atc,
PerAtomQuantity<double> * source,
DENS_MAN* atomCoarseGrainingPositions):
atc_(atc),
source_(source),
atomCoarseGrainingPositions_(atomCoarseGrainingPositions),
feMesh_((atc_->fe_engine())->fe_mesh())
{
source_->register_dependence(this);
atomCoarseGrainingPositions_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyMeshAccumulation::reset_quantity() const
{
const DENS_MAT & positions(atomCoarseGrainingPositions_->quantity());
const DENS_MAT & source(source_->quantity());
int nNodes = feMesh_->num_nodes_unique();
int nodesPerElement = feMesh_->num_nodes_per_element();
Array<int> node_list(nodesPerElement);
DENS_VEC shp(nodesPerElement);
quantity_.resize(nNodes,source.nCols());
_quantityLocal_.reset(nNodes,source.nCols());
DENS_VEC xj(atc_->nsd());
if (source.nRows()>0) {
for (int j = 0; j < source.nRows(); j++) {
for (int k = 0; k < atc_->nsd(); k++) {
xj(k) = positions(j,k);
}
feMesh_->shape_functions(xj,shp,node_list);
for (int I = 0; I < nodesPerElement; I++) {
int inode = node_list(I);
for (int k = 0; k < source.nCols(); k++) {
//quantity_(inode,k) += shp(I)*source(j,k);
_quantityLocal_(inode,k) += shp(I)*source(j,k);
}
}
}
}
// accumulate across processors
int count = quantity_.nRows()*quantity_.nCols();
lammpsInterface_->allsum(_quantityLocal_.ptr(),quantity_.ptr(),count);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyMeshAccumulationNormalized
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyMeshAccumulationNormalized::OnTheFlyMeshAccumulationNormalized(ATC_Method * atc,
PerAtomQuantity<double> * source,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN * weights):
OnTheFlyMeshAccumulation(atc,source,atomCoarseGrainingPositions),
weights_(weights)
{
weights_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyMeshAccumulationNormalized::reset_quantity() const
{
OnTheFlyMeshAccumulation::reset_quantity();
quantity_ *= weights_->quantity();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyMeshAccumulationNormalizedReferenced
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyMeshAccumulationNormalizedReferenced::OnTheFlyMeshAccumulationNormalizedReferenced(ATC_Method * atc,
PerAtomQuantity<double> * source,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN* weights,
DENS_MAN * reference):
OnTheFlyMeshAccumulationNormalized(atc,source,atomCoarseGrainingPositions,weights),
reference_(reference)
{
reference_->register_dependence(this);
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyMeshAccumulationNormalizedReferenced::reset_quantity() const
{
OnTheFlyMeshAccumulationNormalized::reset_quantity();
quantity_ -= reference_->quantity();
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyMeshAccumulationNormalizedScaled
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyMeshAccumulationNormalizedScaled::OnTheFlyMeshAccumulationNormalizedScaled(ATC_Method * atc,
PerAtomQuantity<double> * source,
DENS_MAN* atomCoarseGrainingPositions,
DIAG_MAN* weights,
const double scale):
OnTheFlyMeshAccumulationNormalized(atc,source,atomCoarseGrainingPositions,weights),
scale_(scale)
{
}
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyMeshAccumulationNormalizedScaled::reset_quantity() const
{
OnTheFlyMeshAccumulationNormalized::reset_quantity();
quantity_ *= scale_;
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class NativeShapeFunctionGradient
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// constructor
//--------------------------------------------------------
NativeShapeFunctionGradient::NativeShapeFunctionGradient(ATC_Method * atc) :
VectorTransfer<SPAR_MAT * >(),
feEngine_(atc->fe_engine())
{
quantity_.resize(atc->nsd(),NULL);
for (int i = 0; i < atc->nsd(); ++i) {
quantity_[i] = new SPAR_MAT();
}
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
NativeShapeFunctionGradient::~NativeShapeFunctionGradient() {
for (unsigned i = 0; i < quantity_.size(); ++i) {
if (quantity_[i]) delete quantity_[i];
}
}
//--------------------------------------------------------
// reset_quantity()
//--------------------------------------------------------
void NativeShapeFunctionGradient::reset_quantity() const
{
feEngine_->compute_gradient_matrix(quantity_);
}
//--------------------------------------------------------
//--------------------------------------------------------
// Class OnTheFlyShapeFunctionProlongation
//--------------------------------------------------------
//--------------------------------------------------------
//--------------------------------------------------------
// Constructor
//--------------------------------------------------------
OnTheFlyShapeFunctionProlongation::OnTheFlyShapeFunctionProlongation(ATC_Method * atc,
DENS_MAN * source,
DENS_MAN * atomCoarseGrainingPositions):
FeToAtomTransfer(atc,source),
atomCoarseGrainingPositions_(atomCoarseGrainingPositions),
feMesh_((atc->fe_engine())->fe_mesh())
{
atomCoarseGrainingPositions_->register_dependence(this);
}
//--------------------------------------------------------
// destructor
//--------------------------------------------------------
OnTheFlyShapeFunctionProlongation::~OnTheFlyShapeFunctionProlongation()
{
atomCoarseGrainingPositions_->remove_dependence(this);
};
//--------------------------------------------------------
// reset_quantity
//--------------------------------------------------------
void OnTheFlyShapeFunctionProlongation::reset() const
{
if (need_reset()) {
PerAtomQuantity<double>::reset();
const DENS_MAT & positions(atomCoarseGrainingPositions_->quantity());
const DENS_MAT & source(source_->quantity());
int nodesPerElement = feMesh_->num_nodes_per_element();
Array<int> node_list(nodesPerElement);
DENS_VEC shp(nodesPerElement);
DENS_VEC xj(positions.nCols());
int nLocal = positions.nRows();
quantity_ = 0.;
for (int j = 0; j < nLocal; j++) {
for (int k = 0; k < source.nCols(); k++) {
xj(k) = positions(j,k);
}
feMesh_->shape_functions(xj,shp,node_list);
for (int I = 0; I < nodesPerElement; I++) {
int inode = node_list(I);
for (int k = 0; k < source.nCols(); k++) {
quantity_(j,k) += shp(I)*source(inode,k);
}
}
}
}
}
}
Event Timeline
Log In to Comment