TransferLibrary.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, Jul 13, 23:27

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 ~ phiM^{-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);
	}
	}
	}
	}
	}
	}

TransferLibrary.cppNo OneTemporaryActions

File Metadata

TransferLibrary.cppView Options

Event Timeline

TransferLibrary.cpp
No OneTemporary
Actions

TransferLibrary.cpp
View Options