Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/lib/atc/TransferLibrary.cpp
8d47d9dc3965efficient_neuronet

Last Change
Raw File

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);
}
}
}
}
}
}
c4science · Help