ATC_TransferHardy.cpp
f70511be8768lammm-devel

ATC_TransferHardy.cpp
View Options

			// ATC_Transfer headers ?
			#include "ATC_TransferHardy.h"
			#include "ATC_Error.h"
			#include "FE_Engine.h"
			#include "ATC_HardyKernel.h"
			#include "LammpsInterface.h"
			#include "Quadrature.h"

			// Other Headers
			#include <vector>
			#include <map>
			#include <set>
			#include <utility>
			#include <fstream>
			#include <exception>

			static int sgn(double x) { return (int) ((x>0) - (x<0)); }
			static int rnd(double x) { return (int) (x+sgn(x)*0.5); }

			using namespace std;

			namespace ATC {

			ATC_TransferHardy::ATC_TransferHardy(string groupName,
			string matParamFile)
			: ATC_Transfer(),
			setRefPE_(false),
			setRefPEvalue_(false),
			xPointer_(NULL),
			kernelFunction_(NULL),
			cauchyBornStress_(NULL),
			bondOnTheFly_(false),
			kernelOnTheFly_(false),
			useAtomicShapeFunctions_(true)
			{
			// assign default "internal" group
			int igroup = lammpsInterface_->find_group(groupName.c_str());
			groupbit_ \|= lammpsInterface_->group_bit(igroup);
			igroups_.insert(igroup);

			nTypes_ = lammpsInterface_->ntypes();


			// Defaults
			simTime_ = 0.;

			fieldFlags_.reset(NUM_HARDY_FIELDS);
			fieldFlags_ = true;
			gradFlags_.reset(NUM_HARDY_FIELDS);
			gradFlags_ = false;
			rateFlags_.reset(NUM_HARDY_FIELDS);
			rateFlags_ = false;

			fieldSizes_.reset(NUM_HARDY_FIELDS);
			fieldSizes_ = 1;
			fieldSizes_(HARDY_DISPLACEMENT) =3;
			fieldSizes_(HARDY_MOMENTUM) =3;
			fieldSizes_(HARDY_VELOCITY) =3;
			fieldSizes_(HARDY_PROJECTED_VELOCITY) =3;
			fieldSizes_(HARDY_HEAT_FLUX) =3;
			fieldSizes_(HARDY_STRESS) =9; // LAGRANGIAN
			fieldSizes_(HARDY_ESHELBY_STRESS) =9;
			fieldSizes_(HARDY_CAUCHY_BORN_STRESS) =6;
			fieldSizes_(HARDY_TRANSFORMED_STRESS) =9;
			fieldSizes_(HARDY_TYPE_CONCENTRATION) =nTypes_;

			atomicOutputMask_.insert("displacement");
			atomicOutputMask_.insert("velocity");
			atomicOutputMask_.insert("temperature");
			atomicOutputMask_.insert("potential_energy");
			atomicOutputMask_.insert("force");
			atomicOutputMask_.insert("virial");
			atomicOutputMask_.insert("centrosymmetry");

			set_line_quadrature(line_ngauss,line_xg,line_wg);

			int computeID;
			// create compute for pe/atom
			computeID = lammpsInterface_->atomPE_create();
			if (lammpsInterface_->comm_rank() == 0 ) {
			cout << "atomPE compute created with ID: " << computeID << "\n" << flush;
			}

			}

			//-------------------------------------------------------------------
			ATC_TransferHardy::~ATC_TransferHardy()
			{
			// Deallocate
			if (kernelFunction_) delete kernelFunction_;
			}

			//-------------------------------------------------------------------
			// called before the beginning of a "run"
			void ATC_TransferHardy::initialize()
			{
			// Base class initalizations
			if (!initialized_) {
			feEngine_->initialize();
			nNodes_ = feEngine_->get_nNodes();
			nsd_ = feEngine_->get_nsd();
			// set ref positions
			set_xref();

			// set nlocal & internal to atom map
			reset_nlocal();
			if (useAtomicShapeFunctions_) {
			try {
			reset_shape_functions();
			}
			catch(bad_alloc&) {
			if (lammpsInterface_->comm_rank() == 0) {
			cout << " ATC:: insufficient memory for reset_shape_functions() to execute\n";
			throw;
			}
			}
			set_atomic_weights();
			}

			NodeVolumes_.reset(nNodes_,nNodes_);
			invNodeVolumes_.reset(nNodes_,nNodes_);
			feEngine_->compute_lumped_mass_matrix(NodeVolumes_);
			invNodeVolumes_ = inv(NodeVolumes_);
			}

			// positions
			set_xPointer();

			// set up workspaces
			if (!initialized_) {

			nNodesGlobal_ = feEngine_->get_feMesh()->get_nNodes();

			// get periodicity and box bounds/lengths
			lammpsInterface_->get_box_periodicity(periodicity[0],
			periodicity[1],periodicity[2]);
			lammpsInterface_->get_box_bounds(box_bounds[0][0],box_bounds[1][0],
			box_bounds[0][1],box_bounds[1][1],
			box_bounds[0][2],box_bounds[1][2]);

			for (int k = 0; k < 3; k++) {
			box_length[k] = box_bounds[1][k] - box_bounds[0][k];
			}

			if (atomToElementMapType_ == EULERIAN) {
			fieldSizes_(HARDY_STRESS) =6;
			}

			// PE compute
			lammpsInterface_->atomPE_init();

			// ground state for PE
			nodalRefPotentialEnergy_.reset(nNodes_,1);

			// copy requested fieldFlags_ settings to outputFlags_
			outputFlags_ = fieldFlags_;

			// perform consistency check on fieldFlags
			check_fieldFlags_consistency();

			// check whether single_enable==0 for stress/heat flux calculation
			if (fieldFlags_(HARDY_STRESS) \|\| fieldFlags_(HARDY_HEAT_FLUX)) {
			if (lammpsInterface_->single_enable()==0) {
			throw ATC_Error(0,"Calculation of Hardy stress field not possible with selected pair type.");
			}
			}

			// size arrays for requested/required Hardy fields
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			string name;
			int size;
			if (fieldFlags_(index)) {
			hardy_field_to_string(index,name);
			size = fieldSizes_(index);
			hardyData_ [name].reset(nNodes_,size);
			hardyDataOld_ [name].reset(nNodes_,size);
			filteredHardyData_[name].reset(nNodes_,size);
			}
			}

			// compute table of Hardy bond functions for processors that own atoms
			if ((! bondOnTheFly_)
			&& ( ( fieldFlags_(HARDY_STRESS)
			\|\| fieldFlags_(HARDY_ESHELBY_STRESS)
			\|\| fieldFlags_(HARDY_HEAT_FLUX) ) ) ) {
			try {
			compute_bond_matrix();
			}
			catch(bad_alloc&) {
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "\n ATC:: stress/heat_flux will be computed on-the-fly\n";
			bondOnTheFly_ = true;
			}
			}
			}

			// compute shape function matrix for kernel functions
			if (kernelFunction_ && (! kernelOnTheFly_)) {
			try{
			compute_kernel_matrix();
			}
			catch(bad_alloc&) {
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "\n ATC:: kernel will be computed on-the-fly\n";
			kernelOnTheFly_ = true;
			}
			}
			}

			// compute table of gradient of Hardy functions
			if (fieldFlags_(HARDY_ESHELBY_STRESS)
			\|\| fieldFlags_(HARDY_CAUCHY_BORN_STRESS)
			\|\| fieldFlags_(HARDY_VACANCY_CONCENTRATION)) {
			gradFlags_(HARDY_DISPLACEMENT) = true;
			}
			gradientTable_ = vector<SPAR_MAT>(nsd_);
			if (gradFlags_.has_member(true)) {
			compute_gradient_matrix();
			}

			// construct & initialize filter
			init_filter();
			}
			initialized_ = true;

			// addstep needs to be done _every_ initialize not just first
			lammpsInterface_->atomPE_addstep(lammpsInterface_->ntimestep()+1);

			if (lammpsInterface_->comm_rank() == 0) {
			cout << " ATC:: conversion factor for energy/vol -> stress "
			<< lammpsInterface_->nktv2p() << "\n";
			cout << " ATC:: cutoff radius "
			<< lammpsInterface_->pair_cutoff() << "\n";
			}
			}

			//-------------------------------------------------------------------
			// sets initial values of filtered quantities
			void ATC_TransferHardy::init_filter()
			{
			timeFilters_.reset(NUM_HARDY_FIELDS);
			sampleCounter_ = 0;
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			string name;
			hardy_field_to_string(index,name);
			filteredHardyData_[name] = 0.0;
			timeFilters_(index) = timeFilterManager_.construct();
			}
			}


			//-------------------------------------------------------------------
			// called after the end of a "run"
			void ATC_TransferHardy::finish()
			{
			// base class
			ATC_Transfer::finish();
			}

			//-------------------------------------------------------------------
			// this is the parser
			bool ATC_TransferHardy::modify(int narg, char **arg)
			{
			bool match = false;
			double val, dval;

			// check to see if it is a transfer class command
			if (strcmp(arg[0],"transfer")==0) {
			/*! \page man_hardy_fields fix_modify AtC transfer fields
			\section syntax
			fix_modify AtC transfer fields <all \| none> \n
			fix_modify AtC transfer fields <add \| delete> <list_of_fields> \n
			- all \| none (keyword) = output all or no fields \n
			- add \| delete (keyword) = add or delete the listed output fields \n
			- fields (keyword) = \n
			density : mass per unit volume \n
			displacement : displacement vector \n
			momentum : momentum per unit volume \n
			velocity : defined by momentum divided by density \n
			projected_velocity : simple kernel estimation of atomic velocities \n
			temperature : temperature derived from the relative atomic kinetic energy (as done by Hardy) \n
			kinetic_temperature : temperature derived from the full kinetic energy \n
			energy : total energy (potential + kinetic) per unit volume \n
			number_density : simple kernel estimation of number of atoms per unit volume \n
			stress :
			Cauchy stress tensor for eulerian analysis (atom_element_map), or
			1st Piola-Kirchhoff stress tensor for lagrangian analysis \n
			transformed_stress :
			1st Piola-Kirchhoff stress tensor for eulerian analysis (atom_element_map), or
			Cauchy stress tensor for lagrangian analysis \n
			heat_flux : spatial heat flux vector for eulerian,
			or referential heat flux vector for lagrangian \n
			\section examples
			<TT> fix_modify AtC transfer fields add velocity temperature </TT>
			\section description
			Allows modification of the fields calculated and output by the Hardy
			transfer class. The commands are cumulative, e.g.\n
			<TT> fix_modify AtC transfer fields none </TT> \n
			followed by \n
			<TT> fix_modify AtC transfer fields add velocity temperature </TT> \n
			will only output the velocity and temperature fields.
			\section restrictions
			Must be used with the hardy AtC transfer, see \ref man_fix_atc.
			Currently, the stress and heat flux formulas are only correct for
			central force potentials, e.g. Lennard-Jones and EAM
			but not Stillinger-Weber.
			\section related
			See \ref man_hardy_gradients , \ref man_hardy_rates and \ref man_hardy_computes
			\section default
			All fields are output by default
			*/
			if (strcmp(arg[1],"fields")==0) {
			if (strcmp(arg[2],"all")==0) {
			fieldFlags_ = true;
			match = true;
			}
			else if (strcmp(arg[2],"none")==0) {
			fieldFlags_ = false;
			match = true;
			}
			else if (strcmp(arg[2],"add")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			fieldFlags_(field_name) = true; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			else if (strcmp(arg[2],"delete")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			fieldFlags_(field_name) = false; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			}

			/*! \page man_hardy_gradients fix_modify AtC transfer gradients
			\section syntax
			fix_modify AtC transfer gradients <add \| delete> <list_of_fields> \n
			- add \| delete (keyword) = add or delete the calculation of gradients for the listed output fields \n
			- fields (keyword) = \n
			gradients can be calculated for all fields listed in \ref man_hardy_fields

			\section examples
			<TT> fix_modify AtC transfer gradients add temperature velocity stress </TT> \n
			<TT> fix_modify AtC transfer gradients delete velocity </TT> \n
			\section description
			Requests calculation and ouput of gradients of the fields from the Hardy
			transfer class. These gradients will be with regard to spatial or material
			coordinate for eulerian or lagrangian analysis, respectively, as specified by
			atom_element_map (see \ref man_atom_element_map )
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			No gradients are calculated by default
			*/
			else if (strcmp(arg[1],"gradients")==0) {
			if (strcmp(arg[2],"add")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			gradFlags_(field_name) = true; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			else if (strcmp(arg[2],"delete")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			gradFlags_(field_name) = false; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			}

			/*! \page man_hardy_rates fix_modify AtC transfer rates
			\section syntax
			fix_modify AtC transfer rates <add \| delete> <list_of_fields> \n
			- add \| delete (keyword) = add or delete the calculation of rates (time derivatives) for the listed output fields \n
			- fields (keyword) = \n
			rates can be calculated for all fields listed in \ref man_hardy_fields

			\section examples
			<TT> fix_modify AtC transfer rates add temperature velocity stress </TT> \n
			<TT> fix_modify AtC transfer rates delete stress </TT> \n
			\section description
			Requests calculation and ouput of rates (time derivatives) of the fields from the Hardy
			transfer class. For eulerian analysis (see \ref man_atom_element_map ), these rates
			are the partial time derivatives of the nodal fields, not the full (material) time
			derivatives. \n
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			No rates are calculated by default
			*/
			else if (strcmp(arg[1],"rates")==0) {
			if (strcmp(arg[2],"add")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			rateFlags_(field_name) = true; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			else if (strcmp(arg[2],"delete")==0) {
			hardyFieldName field_name;
			for (int i = 3; i < narg; ++i) {
			if (string_to_hardy_field(arg[i],field_name)) {
			rateFlags_(field_name) = false; }
			else { throw ATC_Error(0,"unsupported Hardy field"); }
			}
			match = true;
			}
			}

			/*! \page man_hardy_computes fix_modify AtC transfer computes
			\section syntax
			fix_modify AtC transfer computes <add \| delete> [per-atom compute id] <volume \| number> \n
			- add \| delete (keyword) = add or delete the calculation of an equivalent continuum field
			for the specified per-atom compute as volume or number density quantity \n
			- per-atom compute id = name/id for per-atom compute,
			fields can be calculated for all per-atom computes available from LAMMPS \n
			- volume \| number (keyword) = field created is a per-unit-volume quantity
			or a per-atom quantity as weighted by kernel functions \n

			\section examples
			<TT> compute virial all stress/atom </TT> \n
			<TT> fix_modify AtC transfer computes add virial volume </TT> \n
			<TT> fix_modify AtC transfer computes delete virial </TT> \n
			\n
			<TT> compute centrosymmetry all centro/atom </TT> \n
			<TT> fix_modify AtC transfer computes add centrosymmetry number </TT> \n
			\section description
			Calculates continuum fields corresponding to specified per-atom computes created by LAMMPS \n
			\section restrictions
			Must be used with the hardy AtC transfer ( see \ref man_fix_atc ) \n
			Per-atom compute must be specified before corresponding continuum field can be requested \n
			\section related
			See manual page for compute
			\section default
			No defaults exist for this command
			*/
			else if (strcmp(arg[1],"computes")==0) {
			if (strcmp(arg[2],"add")==0) {
			string tag(arg[3]);
			int icompute = lammpsInterface_->find_compute(tag.c_str());
			if (icompute < 0)
			throw ATC_Error(0,"Could not find compute "+tag);
			int normalization = NO_NORMALIZATION;
			if (narg > 4) {
			if (strcmp(arg[4],"volume")==0) {
			normalization = VOLUME_NORMALIZATION;
			}
			else if (strcmp(arg[4],"number")==0) {
			normalization = NUMBER_NORMALIZATION;
			}
			}
			computes_[tag] = normalization;
			match = true;
			}
			else if (strcmp(arg[2],"delete")==0) {
			string tag(arg[3]);
			if (computes_.find(tag) != computes_.end()) {
			computes_.erase(tag);
			}
			else {
			throw ATC_Error(0,tag+" compute is not in list");
			}
			match = true;
			}
			}

			/*! \page man_hardy_kernel fix_modify AtC kernel
			\section syntax
			fix_modify AtC transfer kernel <type> <parameters>
			- type (keyword) = mesh, step, cell, cubic_cylinder, cubic_sphere, quartic_cylinder, quartic_sphere \n
			- parameters :\n
			mesh = none\n
			step = radius (double) \n
			cell = hx, hy, hz (double) or h (double) \n
			cubic_cylinder = radius (double) \n
			cubic_sphere = radius (double) \n
			quartic_cylinder = radius (double) \n
			quartic_sphere = radius (double) \n
			\section examples
			fix_modify AtC transfer kernel cell 1.0 1.0 1.0
			\section description

			\section restrictions
			Must be used with the hardy AtC transfer \n
			For cylinder kernel types, cylindrical axis is assumed to be in z-direction \n
			( see \ref man_fix_atc )
			\section related
			\section default
			Default to the mesh based kernel
			*/
			else if (strcmp(arg[1],"kernel")==0) {
			if (kernelFunction_) delete kernelFunction_;
			if (strcmp(arg[2],"mesh")==0) {
			match = true;
			useAtomicShapeFunctions_ = true;
			}
			else if (strcmp(arg[2],"step")==0) {
			double parameters[1] = {atof(arg[3])};
			kernelFunction_ = new ATC_HardyKernelStep(1,parameters);
			match = true;
			}
			else if (strcmp(arg[2],"cell")==0) {
			double parameters[3];
			parameters[0] = parameters[1] = parameters[2] = atof(arg[3]);
			if (narg > 5) {
			for (int i = 1; i < 3; i++) { parameters[i] = atof(arg[3+i]); }
			}
			kernelFunction_ = new ATC_HardyKernelCell(2,parameters);
			match = true;
			}
			else if (strcmp(arg[2],"cubic_cylinder")==0) {
			double parameters[1] = {atof(arg[3])}; // cutoff radius
			kernelFunction_ = new ATC_HardyKernelCubicCyl(1,parameters);
			match = true;
			}
			else if (strcmp(arg[2],"cubic_sphere")==0) {
			double parameters[1] = {atof(arg[3])}; // cutoff radius
			kernelFunction_ = new ATC_HardyKernelCubicSphere(1,parameters);
			match = true;
			}
			else if (strcmp(arg[2],"quartic_cylinder")==0) {
			double parameters[1] = {atof(arg[3])}; // cutoff radius
			kernelFunction_ = new ATC_HardyKernelQuarticCyl(1,parameters);
			match = true;
			}
			else if (strcmp(arg[2],"quartic_sphere")==0) {
			double parameters[1] = {atof(arg[3])}; // cutoff radius
			kernelFunction_ = new ATC_HardyKernelQuarticSphere(1,parameters);
			match = true;
			}
			}

			/*! \page man_hardy_on_the_fly fix_modify AtC transfer on_the_fly
			\section syntax
			fix_modify AtC transfer on_the_fly <bond \| kernel> <optional on \| off> \n
			- bond \| kernel (keyword) = specifies on-the-fly calculation of bond or kernel
			matrix elements \n
			- on \| off (keyword) = activate or discontinue on-the-fly mode \n

			\section examples
			<TT> fix_modify AtC transfer on_the_fly bond on </TT> \n
			<TT> fix_modify AtC transfer on_the_fly kernel </TT> \n
			<TT> fix_modify AtC transfer on_the_fly kernel off </TT> \n
			\section description
			Overrides normal mode of pre-calculating and storing bond pair-to-node and
			kernel atom-to-node matrices. If activated, will calculate elements of these
			matrices during repeated calls of field computations (i.e. "on-the-fly") and not store them for
			future use. \n
			on flag is optional - if omitted, on_the_fly will be activated for the specified
			matrix. Can be deactivated using off flag. \n
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			By default, on-the-fly calculation is not active (i.e. off). However, code does a memory allocation
			check to determine if it can store all needed bond and kernel matrix elements. If this allocation
			fails, on-the-fly is activated. \n
			*/
			else if (strcmp(arg[1],"on_the_fly")==0) {
			if (strcmp(arg[2],"bond")==0) {
			bondOnTheFly_ = true;
			if (narg > 3 && strcmp(arg[3],"off")==0) bondOnTheFly_ = false;
			}
			else if (strcmp(arg[2],"kernel")==0) {
			kernelOnTheFly_ = true;
			if (narg > 3 && strcmp(arg[3],"off")==0) kernelOnTheFly_ = false;
			}
			else { throw ATC_Error(0,"unsupported on_the_fly type"); }
			match = true;
			}

			/*! \page man_hardy_set fix_modify AtC set
			\section syntax
			fix_modify AtC transfer set reference_potential_energy <value>
			- value (double) : optional user specified zero point for PE
			\section examples
			<TT> fix_modify AtC transfer set reference_potential_energy </TT> \n
			<TT> fix_modify AtC transfer set reference_potential_energy -0.05 </TT> \n
			\section description
			Used to set various quantities for the post-processing algorithms.
			Currently it only
			sets the zero point for the potential energy density using
			the value provided for all nodes, or from the current
			configuration of the lattice if no value is provided
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			Defaults to lammps zero point i.e. isolated atoms
			*/
			else if (strcmp(arg[1],"set")==0) {
			if (strcmp(arg[2],"reference_potential_energy")==0) {
			if (narg > 3) {
			double value = atof(arg[3]);
			nodalRefPEvalue_ = value;
			setRefPEvalue_ = true;
			}
			else { // NOTE make this part of initialize
			setRefPE_ = true;
			}
			match = true;
			}
			} // end "set"

			/*! \page man_boundary_integral fix_modify AtC boundary_integral
			\section syntax
			fix_modify AtC transfer boundary_integral [field] faceset [name]
			- field (string) : name of hardy field
			- name (string) : name of faceset
			\section examples
			<TT> fix_modify AtC transfer boundary_integral stress faceset loop1 </TT> \n
			\section description
			Calculates a surface integral of the given field dotted with the
			outward normal of the faces and puts output in the "GLOBALS" file
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			none
			*/
			else if (strcmp(arg[1],"boundary_integral")==0) {
			hardyFieldName field;
			// if (string_to_hardy_field(arg[2],field)) { }
			// else { throw ATC_Error(0,"unsupported Hardy field"); }
			if(strcmp(arg[3],"faceset")==0) {
			string name(arg[4]);
			string field(arg[2]);
			const set< pair<int,int> > * faceSet
			= & ( (feEngine_->get_feMesh())->get_faceset(name));
			pair <string,string> pair_name(name,field);
			bndyIntegralData_[pair_name] = faceSet;
			match = true;
			}
			} // end "boundary_integral"

			/*! \page man_contour_integral fix_modify AtC contour_integral
			\section syntax
			fix_modify AtC transfer contour_integral [field] faceset [name] <axis [x \| y \| z]>
			- field (string) : name of hardy field
			- name (string) : name of faceset
			- axis (string) : x or y or z
			\section examples
			<TT> fix_modify AtC transfer contour_integral stress faceset loop1 </TT> \n
			\section description
			Calculates a surface integral of the given field dotted with the
			outward normal of the faces and puts output in the "GLOBALS" file
			\section restrictions
			Must be used with the hardy AtC transfer
			( see \ref man_fix_atc )
			\section related
			\section default
			none
			*/
			else if (strcmp(arg[1],"contour_integral")==0) {
			hardyFieldName field;
			// if (string_to_hardy_field(arg[2],field)) { }
			// else { throw ATC_Error(0,"unsupported Hardy field"); }
			if(strcmp(arg[3],"faceset")==0) {
			string name(arg[4]);
			string field(arg[2]);
			const set< pair<int,int> > * faceSet
			= & ( (feEngine_->get_feMesh())->get_faceset(name));
			pair <string,string> pair_name(name,field);
			contourIntegralData_[pair_name] = faceSet;
			match = true;
			}
			} // end "contour_integral"
			}

			// no match, call base class parser
			if (!match) {
			match = ATC_Transfer::modify(narg, arg);
			}

			return match;
			}

			//-------------------------------------------------------------------
			// called at the beginning of a timestep
			void ATC_TransferHardy::pre_init_integrate()
			{
			// output initial configuration
			if (stepCounter_ == 0 && outputFrequency_ > 0) {
			double dt = lammpsInterface_->dt();
			time_filter_pre (dt);
			compute_fields();
			// initialize filtered data
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			string name;
			hardy_field_to_string(index,name);
			filteredHardyData_[name] = hardyData_[name];
			}
			time_filter_post(dt);
			compute_boundary_integrals();
			output();
			}
			}
			//-------------------------------------------------------------------
			// called at the end of first half of a timestep
			void ATC_TransferHardy::post_init_integrate()
			{
			}

			//-------------------------------------------------------------------
			// called at the begining of second half timestep
			void ATC_TransferHardy::pre_final_integrate()
			{
			// cases to recompute transfer matrices:
			// (1) if atoms have changed processors
			bool atomSwitch = atomSwitch_;
			int localSwitch = atomSwitch_, globalSwitch = 0;
			lammpsInterface_->int_allmax(&localSwitch,&globalSwitch);
			atomSwitch = globalSwitch;
			// (2) if eulerian and we are at the atom_to_element_map reset frequency
			// which is used to track the convection of atoms
			bool eulerianReset = (atomToElementMapType_ == EULERIAN)
			&& (stepCounter_ % atomToElementMapFrequency_ ==0);

			bool needsKernel = ( kernelFunction_ && (! kernelOnTheFly_ ) );

			bool needsBond = (! bondOnTheFly_ ) &&
			(fieldFlags_(HARDY_STRESS)
			\|\| fieldFlags_(HARDY_ESHELBY_STRESS)
			\|\| fieldFlags_(HARDY_HEAT_FLUX));

			// (3) if we suspect neighbor lists changing and are
			// at the neighbor reset frequency
			int ago = lammpsInterface_->neighbor_ago();
			int interval = outputFrequency_;
			if (timeFilterManager_.filter_dynamics() && sampleFrequency_ < interval)
			{ interval = sampleFrequency_; }
			bool neighborReset = (ago <= interval);
			if (neighborResetFrequency_ > 0) {
			neighborReset = (neighborReset
			\|\| (stepCounter_ % neighborResetFrequency_ == 0) );
			}

			if (atomSwitch \|\| eulerianReset) {
			reset_nlocal();
			if (needsKernel) {
			compute_kernel_matrix();
			}
			if (useAtomicShapeFunctions_) {
			reset_shape_functions();
			set_atomic_weights();
			}
			if (needsBond) {
			compute_bond_matrix();
			}
			if (gradFlags_.has_member(true)) {
			compute_gradient_matrix();
			}
			}
			else if (neighborReset && needsBond) {
			check_pair_map(); // will recompute bond matrix if necessary
			}

			}

			//-------------------------------------------------------------------
			// called at the end of second half timestep
			void ATC_TransferHardy::post_final_integrate()
			{
			double dt = lammpsInterface_->dt();
			simTime_ += dt;
			if (dt == 0.0) simTime_ = stepCounter_;
			++stepCounter_;

			bool output_now = ( (outputFrequency_ > 0)
			&& ((stepCounter_ % outputFrequency_ == 0)) );
			bool sample_now = ( (output_now)
			\|\| (timeFilterManager_.filter_dynamics() && (stepCounter_ % sampleFrequency_ == 0)));

			// compute spatially smoothed quantities
			if (sample_now) {
			time_filter_pre (dt, output_now); // NOTE drop output_now
			compute_fields();
			time_filter_post(dt, output_now);
			compute_boundary_integrals();
			}

			// output
			if ( output_now ) output();
			if ( (outputFrequencyAtom_ > 0)
			&& ((stepCounter_ % outputFrequencyAtom_ == 0)) )
			atomic_output();

			// "1" should be outputFrequency_ or sampleFrequency_
			lammpsInterface_->atomPE_addstep(lammpsInterface_->ntimestep()+1);
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_fields(void)
			{
			if (setRefPE_) {
			compute_potential_energy(nodalRefPotentialEnergy_);
			setRefPE_ = false;
			}

			if (setRefPEvalue_) {
			nodalRefPotentialEnergy_ = nodalRefPEvalue_;
			setRefPEvalue_ = false;
			}


			// NOTE remove this
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ** v = lammpsInterface_->vatom();
			double ** x = lammpsInterface_->xatom();
			int atomIdx;

			bool needs_velocity = fieldFlags_(HARDY_VELOCITY);
			bool needs_projected_velocity = fieldFlags_(HARDY_PROJECTED_VELOCITY);
			bool needs_momentum = fieldFlags_(HARDY_MOMENTUM);
			bool needs_density = fieldFlags_(HARDY_DENSITY);
			bool needs_displacement = fieldFlags_(HARDY_DISPLACEMENT);
			bool needs_temperature = fieldFlags_(HARDY_TEMPERATURE);
			bool needs_kinetic_temperature = fieldFlags_(HARDY_KINETIC_TEMPERATURE);
			bool needs_stress = fieldFlags_(HARDY_STRESS);
			bool needs_eshelby_stress = fieldFlags_(HARDY_ESHELBY_STRESS);
			bool needs_cauchy_born_stress = fieldFlags_(HARDY_CAUCHY_BORN_STRESS);
			bool needs_heat = fieldFlags_(HARDY_HEAT_FLUX);
			bool needs_energy = fieldFlags_(HARDY_ENERGY);
			bool needs_number_density = fieldFlags_(HARDY_NUMBER_DENSITY);
			bool needs_transformed_stress = fieldFlags_(HARDY_TRANSFORMED_STRESS);
			bool needs_vacancy_concentration = fieldFlags_(HARDY_VACANCY_CONCENTRATION);
			bool needs_type_concentration = fieldFlags_(HARDY_TYPE_CONCENTRATION);

			// (1) direct quantities
			//compute : kernel norm is a _weighted_ number density
			if (needs_number_density) {
			compute_number_density(hardyData_["number_density"]);
			}
			// compute: density
			if (needs_density) {
			compute_mass_density(hardyData_["density"]);
			}
			// compute: displacement
			if (needs_displacement) {
			compute_displacement(hardyData_["displacement"],hardyData_["density"]);
			}
			// compute: momentum
			if (needs_momentum) {
			compute_momentum(hardyData_["momentum"]);
			}
			//compute: projected velocity
			if (needs_projected_velocity) {
			compute_projected_velocity(hardyData_["projected_velocity"]);
			}
			// compute: velocity
			if (needs_velocity) {
			compute_velocity(hardyData_["velocity"],
			hardyData_["density"], hardyData_["momentum"]);
			// NOTE make this only for mesh
			compute_variation_velocity(uVariationVelocity_, hardyData_["velocity"]);
			}

			// compute : temperature & (full) kinetic temperature
			if (needs_temperature) {
			compute_temperature(hardyData_["temperature"]);
			}
			if (needs_kinetic_temperature) {
			compute_kinetic_temperature(hardyData_["kinetic_temperature"]);
			}
			// compute: stress
			if (needs_stress) {
			compute_stress(hardyData_["stress"]);
			}
			// compute: heat flux
			if (needs_heat) {
			compute_heatflux(hardyData_["heat_flux"]);
			}
			// compute: energy
			if (needs_energy) {
			compute_total_energy(hardyData_["energy"]);
			}

			// (2) derived quantities
			// compute: gradients
			if (gradFlags_.has_member(true)) {
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			if (gradFlags_(index)) {
			string field;
			hardy_field_to_string(index,field);
			string grad_field = field + "_gradient";
			if (hardyData_.find(field) == hardyData_.end() ) {
			throw ATC_Error(0,"field " + field + " needs to be defined for gradient");
			}
			gradient_compute(hardyData_[field], hardyData_[grad_field]);
			}
			}
			}

			// compute:eshelby_stress
			if (fieldFlags_(HARDY_ESHELBY_STRESS)) {
			if (! (gradFlags_(DISPLACEMENT)) ) {
			throw ATC_Error(0,"displacement_gradient needs to be defined for eshelby stress");
			}
			compute_eshelby_stress(hardyData_["eshelby_stress"],
			hardyData_["energy"],hardyData_["stress"],
			hardyData_["displacement_gradient"]);
			}

			// (3) computes
			map <string,int>::const_iterator iter;
			for (iter = computes_.begin(); iter != computes_.end(); iter++) {
			string tag = iter->first;
			int projection = iter->second;
			int ncols = lammpsInterface_->compute_ncols(tag.c_str());;
			DENS_MAT atomicData(nLocal_,ncols);
			if (ncols == 1) {
			double * atomData = lammpsInterface_->compute_scalar_data(tag.c_str());
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			atomicData(i,0) = atomData[atomIdx];
			}
			}
			else {
			double ** atomData = lammpsInterface_->compute_vector_data(tag.c_str());
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			for (int k = 0; k < ncols; k++) {
			atomicData(i,k) = atomData[atomIdx][k];
			}
			}
			}
			if (projection == VOLUME_NORMALIZATION) {
			project_volume_normalized(atomicData, hardyData_[tag]);
			}
			else if (projection == NUMBER_NORMALIZATION) {
			project_count_normalized(atomicData, hardyData_[tag]);
			}
			else {
			ATC_Transfer::restrict_unscaled(atomicData, hardyData_[tag]);
			}
			}


			}// end of compute_fields routine

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_boundary_integrals(void)
			{
			if (! bndyIntegralData_.empty()) {
			map < pair<string,string>, const set< PAIR > * >::const_iterator iter;
			DENS_MAT values;
			for (iter = bndyIntegralData_.begin();
			iter != bndyIntegralData_.end(); iter++) {
			string bndyName = (iter->first).first;
			string fieldName = (iter->first).second;
			const set< PAIR > & faceSet = *(iter->second);
			const DENS_MAT & data = (hardyData_.find(fieldName)->second);
			feEngine_->field_surface_flux(data,faceSet,values);
			for (int i = 0; i < values.nRows() ; ++i ) {
			string name = bndyName + "_" + fieldName + "_"
			+ ATC_STRING::tostring(i+1);
			feEngine_->add_global(name, values(i,0));
			}
			}
			}
			if (! contourIntegralData_.empty()) {
			map < pair<string,string>, const set< PAIR > * >::const_iterator iter;
			DENS_MAT values;
			for (iter = contourIntegralData_.begin();
			iter != contourIntegralData_.end(); iter++) {
			string bndyName = (iter->first).first;
			string fieldName = (iter->first).second;
			const set< PAIR > & faceSet = *(iter->second);
			const DENS_MAT & data = (hardyData_.find(fieldName)->second);
			feEngine_->field_surface_flux(data,faceSet,values,true);
			for (int i = 0; i < values.nRows() ; ++i ) {
			string name = "contour_"+bndyName + "_" + fieldName + "_"
			+ ATC_STRING::tostring(i+1);
			feEngine_->add_global(name, values(i,0));
			}
			}
			}
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_potential_energy(DENS_MAT & nodalPE)
			{
			DENS_MAT atomicEnergy(nLocal_,1);

			// compute pair energy per atom
			double * atomPE = lammpsInterface_->atomPE_compute();

			// add up potential energy per atom
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			atomicEnergy(i,0) += atomPE[atomIdx];
			}
			project_volume_normalized(atomicEnergy, nodalPE );
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_number_density(DENS_MAT & density)
			{
			DENS_MAT atomCnt(nLocal_,1);
			atomCnt = 1;
			project_volume_normalized(atomCnt, density);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_mass_density(DENS_MAT & density)
			{
			atomicDensity_.reset(nLocal_,1);
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			// density
			atomicDensity_(i,0) = ma;
			}
			project_volume_normalized(atomicDensity_, density);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_displacement(DENS_MAT & displacement,
			const DENS_MAT & density)
			{
			atomicDisplacement_.reset(nLocal_,nsd_);
			double ** x = lammpsInterface_->xatom();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			for (int j = 0; j < nsd_; j++) {
			// displacement
			atomicDisplacement_(i,j) = x[atomIdx][j] - xref_[atomIdx][j];
			// the cheap version of unwrapping atomic positions
			if ((bool) periodicity[j]) {
			double u = atomicDisplacement_(i,j);
			if (u >= 0.5*box_length[j]) { u -= box_length[j]; }
			if (u <= -0.5*box_length[j]) { u += box_length[j]; }
			atomicDisplacement_(i,j) = u;
			}
			atomicDisplacement_(i,j) *= ma;
			}
			}
			project_volume_normalized(atomicDisplacement_, displacement);
			for (int i = 0; i < nNodes_; i++) {
			double rho_i = density(i,0);
			if (rho_i > 0.0) {
			for (int j = 0; j < nsd_; j++) {
			displacement(i,j) = displacement(i,j)/rho_i;
			}
			}
			}
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_momentum(DENS_MAT & momentum)
			{
			atomicMomentum_.reset(nLocal_,nsd_);
			double ** v = lammpsInterface_->vatom();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			// momentum
			for (int j = 0; j < nsd_; j++) {
			atomicMomentum_(i,j) = ma*(v[atomIdx][j]);
			}
			}
			project_volume_normalized(atomicMomentum_, momentum);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_projected_velocity(DENS_MAT & velocity)
			{
			atomicVelocity_.reset(nLocal_,nsd_);
			double ** v = lammpsInterface_->vatom();
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			for (int j = 0; j < nsd_; j++) {
			// velocity
			atomicVelocity_(i,j) += v[atomIdx][j];
			}
			}
			project_count_normalized(atomicVelocity_, velocity);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_velocity(DENS_MAT & velocity,
			const DENS_MAT & density,
			const DENS_MAT & momentum)
			{
			velocity.reset(nNodes_,nsd_);
			for (int i = 0; i < nNodes_; i++) {
			double rho_i = density(i,0);
			if (rho_i > 0.0) {
			for (int j = 0; j < nsd_; j++) {
			velocity(i,j) = momentum(i,j)/rho_i;
			}
			}
			}
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_variation_velocity(DENS_MAT & velocity,
			const DENS_MAT & vI)
			{
			if (nLocal_>0) {
			// interpolate nodal velocities to the atoms
			vbar_.reset(nLocal_,nsd_);
			double ** v = lammpsInterface_->vatom();
			// use of prolong assumes atom system contained within mesh
			if (useAtomicShapeFunctions_) {
			ATC_Transfer::prolong(vI,vbar_);
			}
			// compute and store variation velocities of atoms
			uVariationVelocity_.reset(nLocal_,nsd_);
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			for (int j = 0; j < nsd_; j++) {
			velocity(i,j) = v[atomIdx][j] - vbar_(i,j);
			}
			}
			}
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_temperature(DENS_MAT & temperature)
			{
			atomicTemperature_.reset(nLocal_,1);
			double kB = lammpsInterface_->kBoltzmann();
			double Tcoef = 1./(nsd_*kB);
			double ** v = lammpsInterface_->vatom();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			for (int j = 0; j < nsd_; j++) {
			// variation temperature
			atomicTemperature_(i,0) += Tcoefma(uVariationVelocity_(i,j)*uVariationVelocity_(i,j));
			}
			}
			project_count_normalized(atomicTemperature_, temperature);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_kinetic_temperature(DENS_MAT & temperature)
			{
			atomicKineticTemperature_.reset(nLocal_,1);
			double kB = lammpsInterface_->kBoltzmann();
			double Tcoef = 1./(nsd_*kB);
			double ** v = lammpsInterface_->vatom();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			for (int j = 0; j < nsd_; j++) {
			// full temperature
			atomicKineticTemperature_(i,0) += Tcoefma(v[atomIdx][j]*v[atomIdx][j]);
			}
			}
			project_count_normalized(atomicKineticTemperature_, temperature);
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_stress(DENS_MAT & stress)
			{
			// table of bond functions already calculated in initialize function
			// get conversion factor for mvv to e units
			double mvv2e = lammpsInterface_->mvv2e();
			// get conversion factor for nktV to p units
			double nktv2p = lammpsInterface_->nktv2p();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;

			// calculate kinetic energy tensor part of stress for Eulerian analysis
			DENS_MAT & vI = uVariationVelocity_;
			if (atomToElementMapType_ == EULERIAN && nLocal_>0) {
			atomicStress_.reset(nLocal_,6);
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			// convert mass to appropriate units
			ma = mvv2e*ma;
			atomicStress_(i,0) -= mavI(i,0)vI(i,0);
			atomicStress_(i,1) -= mavI(i,1)vI(i,1);
			atomicStress_(i,2) -= mavI(i,2)vI(i,2);
			atomicStress_(i,3) -= mavI(i,0)vI(i,1);
			atomicStress_(i,4) -= mavI(i,0)vI(i,2);
			atomicStress_(i,5) -= mavI(i,1)vI(i,2);
			}
			project_volume_normalized(atomicStress_, stress);
			}
			else {
			// zero stress table for Lagrangian analysis or if nLocal_ = 0
			stress.zero();
			}
			// add-in potential part of stress tensor
			int nrows = stress.nRows();
			int ncols = stress.nCols();
			DENS_MAT local_potential_hardy_stress(nrows,ncols);
			if (nLocal_>0) {
			if (bondOnTheFly_) {
			compute_potential_stress(local_potential_hardy_stress);
			}
			else {
			// compute table of force & position dyad
			compute_force_matrix();
			// calculate force part of stress tensor
			local_potential_hardy_stress = atomicBondTable_*atomicForceTable_;
			local_potential_hardy_stress *= 0.5; // NOTE: work around for matrix multiplication
			}
			}
			// global summation of potential part of stress tensor
			DENS_MAT potential_hardy_stress(nrows,ncols);
			int count = nrows*ncols;
			lammpsInterface_->allsum(local_potential_hardy_stress.get_ptr(),
			potential_hardy_stress.get_ptr(), count);
			stress += potential_hardy_stress;
			stress = nktv2p*stress;
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_heatflux(DENS_MAT & flux)
			{
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			// calculate kinetic part of heat flux
			if (atomToElementMapType_ == EULERIAN && nLocal_>0) {
			// get conversion factor for mvv to e units
			double mvv2e = lammpsInterface_->mvv2e();
			// compute pair energy per atom
			double * atomPE = lammpsInterface_->atomPE_compute();
			double atomKE, atomEnergy;
			atomicHeat_.reset(nLocal_,3);
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			// convert mass to appropriate units
			ma = mvv2e*ma;
			atomKE = 0.0;
			for (int j = 0; j < nsd_; j++) {
			atomKE += 0.5ma(uVariationVelocity_(i,j)*uVariationVelocity_(i,j));
			}
			atomEnergy = atomKE + atomPE[atomIdx];
			for (int j = 0; j < nsd_; j++) {
			atomicHeat_(i,j) += atomEnergy*uVariationVelocity_(i,j);
			}
			}
			project_volume_normalized(atomicHeat_,flux);
			}
			else {
			// zero stress table for Lagrangian analysis or if nLocal_ = 0
			flux.zero();
			}
			// add-in potential part of heat flux vector
			int nrows = flux.nRows();
			int ncols = flux.nCols();
			DENS_MAT local_hardy_heat(nrows,ncols);
			if (nLocal_>0) {
			if (bondOnTheFly_) {
			compute_potential_heatflux(local_hardy_heat);
			}
			else {
			// compute table of heat vectors
			compute_heat_matrix();
			// calculate force/potential-derivative part of heat flux
			local_hardy_heat = atomicBondTable_*atomicHeatTable_;
			}
			}
			// global summation of potential part of heat flux vector
			DENS_MAT hardy_heat(nrows,ncols);
			int count = nrows*ncols;
			lammpsInterface_->allsum(local_hardy_heat.get_ptr(),
			hardy_heat.get_ptr(), count);
			flux += hardy_heat;
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_total_energy(DENS_MAT & energy)
			{
			atomicEnergy_.reset(nLocal_,1);
			// compute pair energy per atom
			double * atomPE = lammpsInterface_->atomPE_compute();
			// get conversion factor for mvv to e units
			double mvv2e = lammpsInterface_->mvv2e();
			int * type = lammpsInterface_->atom_type();
			double * mass = lammpsInterface_->atom_mass();
			double * rmass = lammpsInterface_->atom_rmass();
			double ma;
			for (int i = 0; i < nLocal_; i++) {
			int atomIdx = internalToAtom_(i);
			if (mass) ma = mass[type[atomIdx]];
			else ma = rmass[atomIdx];
			// convert mass to appropriate units
			ma = mvv2e*ma;
			// compute kinetic energy per atom
			double atomKE = 0.0;
			for (int k = 0; k < nsd_; k++) {
			atomKE += 0.5ma(uVariationVelocity_(i,k)*uVariationVelocity_(i,k));
			}
			// add up total energy per atom
			atomicEnergy_(i,0) += atomPE[atomIdx] + atomKE;
			}
			project_volume_normalized(atomicEnergy_, energy);
			// subtract zero point energy
			energy -= nodalRefPotentialEnergy_;
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::time_filter_pre(double dt, bool output_now)
			{
			sampleCounter_++;
			string name;
			double delta_t = dt*sampleFrequency_;
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			hardy_field_to_string(index,name);
			timeFilters_(index)->apply_pre_step1(filteredHardyData_[name],
			hardyData_[name], delta_t);
			}
			// NOTE add computes_ here
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::time_filter_post(double dt, bool output_now)
			{
			sampleCounter_++;
			string name;
			double delta_t = dt*sampleFrequency_;
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			hardy_field_to_string(index,name);
			timeFilters_(index)->apply_post_step2(filteredHardyData_[name],
			hardyData_[name], delta_t);
			}
			if (rateFlags_.has_member(true)) {
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			if (rateFlags_(index)) {
			string field;
			hardy_field_to_string(index,field);
			string rate_field = field + "_rate";
			timeFilters_(index)->rate(hardyData_[rate_field],
			filteredHardyData_[field],
			hardyData_[field], delta_t);
			}
			}
			}
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			hardy_field_to_string(index,name);
			if (rateFlags_(index)) {
			string rate_field = name + "_rate";
			filteredHardyData_[rate_field] = hardyData_[rate_field];
			}
			if (gradFlags_(index)) {
			string grad_field = name + "_gradient";
			filteredHardyData_[grad_field] = hardyData_[grad_field];
			}
			}
			// NOTE add computes_ here
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::output()
			{
			if (lammpsInterface_->comm_rank() == 0) {
			// data
			OUTPUT_LIST output_data;
			for(int index=0; index < NUM_HARDY_FIELDS; ++index) {
			string name;
			hardy_field_to_string(index,name);
			if (outputFlags_(index)) {
			output_data[name] = & ( filteredHardyData_[name]);
			}
			if (rateFlags_(index)) {
			string rate_name = name + "_rate";
			output_data[rate_name] = & ( filteredHardyData_[rate_name]);
			}
			if (gradFlags_(index)) {
			string grad_name = name + "_gradient";
			output_data[grad_name] = & ( filteredHardyData_[grad_name]);
			}
			}
			map <string,int>::const_iterator iter;
			for (iter = computes_.begin(); iter != computes_.end(); iter++) {
			string tag = iter->first;
			output_data[tag] = & ( hardyData_[tag]); // NOTE should be filtered
			}

			// output
			feEngine_->write_data(simTime_, & output_data);
			}
			}
			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_kernel_matrix(void)
			{
			if (nLocal_>0) {
			set_xPointer();
			kernelShpFcn_.reset(nNodes_,nLocal_);
			DENS_VEC xI(nsd_),xa(nsd_),xaI(nsd_);
			if (lammpsInterface_->comm_rank() == 0 ) {
			cout << "ATC:: computing kernel matrix " << flush;
			}
			int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
			for (int i = 0; i < nNodes_; i++) {
			if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
			cout << "." << flush;
			}
			// Hardy point
			xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
			//xI = (feEngine_->get_feMesh())->global_coordinates(i);
			int ii = i;
			//int ii = feEngine_->get_feMesh()->map_global_to_unique(i);
			for (int j = 0; j < nLocal_; j++) {
			// atom location
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			// difference vector
			xaI = xa - xI;
			periodicity_correction(xaI);
			// compute kernel value & add it to the matrix
			double val = kernelFunction_->value(xaI);
			if (val > 0) kernelShpFcn_.add(ii,j,val);
			}
			}
			kernelShpFcn_.compress();
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "done\n";
			}
			}
			}

			//-------------------------------------------------------------------
			// count normalized
			void ATC_TransferHardy::project_count_normalized(const DENS_MAT & atomData,
			DENS_MAT & nodeData)
			{
			int ncols_atomData = atomData.nCols();
			nodeData.reset(nNodes_,ncols_atomData);
			nodeData.zero();
			// kernel functions
			if (kernelFunction_ ) {
			project_volume_normalized(atomData,nodeData);
			for (int i = 0; i < nNodes_; i++) {
			double numDens = hardyData_["number_density"](i,0);
			if (numDens > 0) {
			double invNumDens = 1.0/numDens;
			for (int j = 0; j < ncols_atomData; j++) {
			nodeData(i,j) *= invNumDens;
			}
			}
			else {
			for (int j = 0; j < ncols_atomData; j++) {
			nodeData(i,j) = 0;
			}
			}
			}
			}
			// mesh-based kernel functions
			else {
			ATC_Transfer::restrict(atomData,nodeData);
			}
			}

			//-------------------------------------------------------------------
			// volume normalized
			void ATC_TransferHardy::project_volume_normalized(const DENS_MAT & atomData,
			DENS_MAT & nodeData)
			{
			int ncols_atomData = atomData.nCols();
			DENS_MAT workNodeArray(nNodes_,ncols_atomData);
			workNodeArray.zero();
			nodeData.reset(workNodeArray.nRows(),workNodeArray.nCols());
			nodeData.zero();
			// kernel functions
			if (kernelFunction_ ) {
			if (nLocal_>0) {
			double invVol = kernelFunction_->inv_vol();
			// on the fly calculation
			if (kernelOnTheFly_) {
			set_xPointer();
			DENS_VEC xI(nsd_),xa(nsd_),xaI(nsd_);
			double val;
			for (int i = 0; i < nNodes_; i++) {
			xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
			for (int j = 0; j < nLocal_; j++) {
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			xaI = xa - xI;
			periodicity_correction(xaI);
			val = kernelFunction_->value(xaI);
			if (val > 0) {
			for (int k=0; k < ncols_atomData; k++) {
			workNodeArray(i,k) += val*atomData(j,k);
			}
			}
			}
			}
			}
			// matrix calculation
			else {
			workNodeArray = kernelShpFcn_*atomData;
			}
			workNodeArray *= invVol;
			}
			}
			// mesh-based kernel functions
			else {
			// computes nodeData = N*atomData w/ N: volume normalized shape functions
			if (nLocal_>0)
			workNodeArray = invNodeVolumes_*shpFcn_.transMat(atomData);
			}
			// accumulate across processors
			int count = workNodeArray.nRows()*workNodeArray.nCols();
			lammpsInterface_->allsum(workNodeArray.get_ptr(),nodeData.get_ptr(),count);
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::gradient_compute(const DENS_MAT & inNodeData,
			DENS_MAT & outNodeData)
			{
			int nrows = inNodeData.nRows();
			int ncols = inNodeData.nCols();
			outNodeData.reset(nrows,ncols*nsd_);
			for (int j = 0; j < nNodes_; j++) {
			for (int i = 0; i < nNodes_; i++) {
			int index = 0;
			for (int n = 0; n < ncols; n++) { //output v1,1 v1,2 v1,3 ...
			for (int m = 0; m < nsd_; m++) {
			outNodeData(j,index) += gradientTable_[m](j,i)*inNodeData(i,n);
			index++;
			}
			}
			}
			}
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_pair_map()
			{
			// neighbor lists for "internal" group
			// a: a "real" atom in the "internal" group
			// b: any neighbor to a (it should project into the mesh for mesh-based)
			int inum = lammpsInterface_->neighbor_list_inum();
			if (inum != nLocalTotal_)
			throw ATC_Error(0,"size mismatch in neighbor list");
			int *ilist = lammpsInterface_->neighbor_list_ilist();
			int *numneigh = lammpsInterface_->neighbor_list_numneigh();
			int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
			int * mask = lammpsInterface_->atom_mask();

			pairMap_.clear();
			int pair_index = 0;
			pair< int,int > pair_ij;
			for (int i = 0; i < nLocalTotal_; i++) { // NOTE why nLocal__Total__?
			int lammps_i = ilist[i];
			// filter out atoms not in the internal group
			if (mask[lammps_i] & groupbit_) {
			//int tag_i = (lammpsInterface_->atom_tag())[lammps_i];
			for (int j = 0; j < numneigh[lammps_i]; j++) {
			int lammps_j = firstneigh[lammps_i][j];
			//int tag_j = (lammpsInterface_->atom_tag())[lammps_j];
			//if (lammps_i > lammps_j) continue; // skip a > b pairs
			//if (tag_i > tag_j) continue; // skip a > b pairs
			pair_ij.first = lammps_i; // alpha
			pair_ij.second = lammps_j; // beta
			pairMap_[pair_ij] = pair_index;
			pair_index++;
			}
			}
			}
			nPairs_ = pair_index;
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_bond_matrix()
			{
			set_xPointer();
			if (lammpsInterface_->comm_rank() == 0) {
			cout << " ATC:: computing bond matrix " << flush;
			}
			// compute pair map
			compute_pair_map();

			// neighbor lists
			int *numneigh = lammpsInterface_->neighbor_list_numneigh();
			int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();

			Array<bool> latticePeriodicity(3);
			latticePeriodicity(0) = (bool) periodicity[0];
			latticePeriodicity(1) = (bool) periodicity[1];
			latticePeriodicity(2) = (bool) periodicity[2];


			double lam1,lam2;
			double bond_value;
			pair< int,int > pair_jk;
			map< pair< int,int >,int >::iterator pairMapIterator;
			int pair_index;


			atomicBondTable_.reset(nNodes_,nPairs_);

			// process differently for mesh vs translation-invariant kernels
			if (kernelFunction_ ) {
			int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
			//int heartbeatFreq = (nNodesGlobal_ <= 10 ? 1 : (int) nNodesGlobal_ / 10);
			// "normal" kernel functions
			DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
			double kernel_inv_vol = kernelFunction_->inv_vol();
			for (int i = 0; i < nNodes_; i++) {
			//for (int i = 0; i < nNodesGlobal_; i++) {
			if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
			cout << "." << flush;
			}
			int ii = i;
			//int ii = feEngine_->get_feMesh()->map_global_to_unique(i);
			// Hardy point
			xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
			//xI = (feEngine_->get_feMesh())->global_coordinates(i);
			for (pairMapIterator = pairMap_.begin();
			pairMapIterator != pairMap_.end(); pairMapIterator++){
			int lammps_j = (pairMapIterator->first).first ;
			int lammps_k = (pairMapIterator->first).second;
			xa.copy(xPointer_[lammps_j],3);
			xaI = xa - xI;
			periodicity_correction(xaI);
			xb.copy(xPointer_[lammps_k],3);
			xba = xb - xa;
			xbI = xba + xaI;
			// intersect segment with kernel support
			kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
			if (lam1 < lam2) {
			bond_value
			= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
			pair_index = pairMapIterator->second;
			atomicBondTable_.add(ii,pair_index,bond_value);
			} // if lam1 < lam2
			} // pair map
			}// end nodes loop
			atomicBondTable_.compress();
			}
			else {
			// mesh-based kernel functions
			int heartbeatFreq = (int) pairMap_.size() / 10;
			int i=0;
			int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
			Array<int> node_list(nodes_per_element);
			DENS_VEC shp(nodes_per_element);
			DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
			double **xatom = lammpsInterface_->xatom();
			for (pairMapIterator = pairMap_.begin();
			pairMapIterator != pairMap_.end(); pairMapIterator++){
			if (lammpsInterface_->comm_rank() == 0 && i++ % heartbeatFreq == 0 ) {
			cout << "." << flush;
			}
			int lammps_j = (pairMapIterator->first).first ;
			int lammps_k = (pairMapIterator->first).second;
			pair_index = pairMapIterator->second;
			xa.copy(xPointer_[lammps_j],3);
			xb.copy(xPointer_[lammps_k],3);
			lam1 = 0.0; lam2 = 1.0;
			double del_lambda = 0.5*(lam2 - lam1);
			double avg_lambda = 0.5*(lam2 + lam1);
			xab = xa - xb;
			for (int i = 0; i < line_ngauss; i++) {
			double lambda = del_lambda*line_xg[i] +avg_lambda;
			xlambda = lambda*xab + xb;
			// NOTE doesn't work for the case that fe_region \subset md_region
			int dummyEltID;
			feEngine_->shape_functions(xlambda,shp,dummyEltID,
			node_list,latticePeriodicity);
			// accumulate to nodes whose support overlaps the integration point
			for (int I = 0; I < nodes_per_element; I++) {
			// Use factor of 0.5 to account for line integration
			// domain of 0 to 1 instead of -1 to 1
			int inode = node_list(I);
			double val = invNodeVolumes_(inode,inode)shp(I)line_wg[i]*0.5;
			atomicBondTable_.add(inode,pair_index,val);
			}
			}
			}
			atomicBondTable_.compress();
			}// end if kernelFunction_
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "done\n";
			}
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_gradient_matrix()
			{
			feEngine_->compute_gradient_matrix(gradientTable_);
			}

			//-------------------------------------------------------------------
			// checks if atoms have changed neighbors
			void ATC_TransferHardy::check_pair_map()
			{
			// (1) compute fresh temporary pairMap
			// (2) check if tmpPairMap == pairMap_

			// neighbor lists
			int *numneigh = lammpsInterface_->neighbor_list_numneigh();
			int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();

			pair< int,int > pair_ij;
			map< pair< int,int >,int >::iterator pairMapIterator;

			for (int i = 0; i < nLocal_; i++) {
			int lammps_i = internalToAtom_(i);
			for (int k = 0; k < numneigh[lammps_i]; k++) {
			int j = firstneigh[lammps_i][k];
			pair_ij.first = lammps_i;
			pair_ij.second = j;
			pairMapIterator = pairMap_.find(pair_ij);
			// if pair is not found in pairMap_, recompute pair map and bond table
			// NOTE need to check if pairMap_ has a stale pair <<<<
			// construct new pair map and use "=="
			if (pairMapIterator == pairMap_.end()) {
			cout << "ATC::check_pair_map has found that pair map and bond table need to be recomputed" << endl;
			compute_bond_matrix();
			return;
			}
			}
			}
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_force_matrix()
			{
			// reset force table
			if (atomToElementMapType_ == LAGRANGIAN) {
			atomicForceTable_.reset(nPairs_,9);
			}
			else if (atomToElementMapType_ == EULERIAN) {
			atomicForceTable_.reset(nPairs_,6);
			}
			// calculate force table
			double **xatom = lammpsInterface_->xatom();
			map< pair< int,int >,int >::iterator pairMapIterator;
			for (pairMapIterator = pairMap_.begin();
			pairMapIterator != pairMap_.end(); pairMapIterator++){
			int lammps_i = (pairMapIterator->first).first ;
			int lammps_j = (pairMapIterator->first).second;
			int pair_index = pairMapIterator->second;
			double * xi = xatom[lammps_i];
			double * xj = xatom[lammps_j];
			double delx = xi[0] - xj[0];
			double dely = xi[1] - xj[1];
			double delz = xi[2] - xj[2];
			double rsq = delxdelx + delydely + delz*delz;
			double fforce = 0;
			lammpsInterface_->pair_force(lammps_i,lammps_j,rsq,fforce);
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_i][0] - xref_[lammps_j][0];
			double delY = xref_[lammps_i][1] - xref_[lammps_j][1];
			double delZ = xref_[lammps_i][2] - xref_[lammps_j][2];
			atomicForceTable_(pair_index,0)=-delxfforcedelX;
			atomicForceTable_(pair_index,1)=-delxfforcedelY;
			atomicForceTable_(pair_index,2)=-delxfforcedelZ;
			atomicForceTable_(pair_index,3)=-delyfforcedelX;
			atomicForceTable_(pair_index,4)=-delyfforcedelY;
			atomicForceTable_(pair_index,5)=-delyfforcedelZ;
			atomicForceTable_(pair_index,6)=-delzfforcedelX;
			atomicForceTable_(pair_index,7)=-delzfforcedelY;
			atomicForceTable_(pair_index,8)=-delzfforcedelZ;
			}
			else if (atomToElementMapType_ == EULERIAN) {
			atomicForceTable_(pair_index,0)=-delxdelxfforce;
			atomicForceTable_(pair_index,1)=-delydelyfforce;
			atomicForceTable_(pair_index,2)=-delzdelzfforce;
			atomicForceTable_(pair_index,3)=-delxdelyfforce;
			atomicForceTable_(pair_index,4)=-delxdelzfforce;
			atomicForceTable_(pair_index,5)=-delydelzfforce;
			}
			}
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_heat_matrix()
			{
			// reset heat flux and energy tables
			atomicHeatTable_.reset(nPairs_,3);
			// calculate force table
			map< pair< int,int >,int >::iterator pairMapIterator;
			double ** xatom = lammpsInterface_->xatom();
			double fforce;
			for (pairMapIterator = pairMap_.begin();
			pairMapIterator != pairMap_.end(); pairMapIterator++){
			int lammps_i = (pairMapIterator->first).first ;
			int lammps_j = (pairMapIterator->first).second;
			int pair_index = pairMapIterator->second;
			double * xi = xatom[lammps_i];
			double * xj = xatom[lammps_j];
			double energy_i = 0.0;
			double delx = xi[0] - xj[0];
			double dely = xi[1] - xj[1];
			double delz = xi[2] - xj[2];
			double rsq = delxdelx + delydely + delz*delz;
			double pair_pe=lammpsInterface_->pair_force(lammps_i,lammps_j,rsq,fforce);
			fforce *= 0.5;
			// This term is correct ONLY for pair potentials.
			// fforce = - (1/x_{ij})*(d(phi_i)/d(x_{ij}) + d(phi_j)/d(x_{ij}))
			// What we need for the heat flux calculation is the portion of this that belongs only to
			// atom j: - (1/x_{ij})*d(phi_j)/d(x_{ij}).
			int i = lammps_i; // NOTE is this right?
			fforce = (delxuVariationVelocity_(i,0) +
			dely*uVariationVelocity_(i,1) +
			delz*uVariationVelocity_(i,2));
			// NOTE REJ don't think this is right, where is energy_i from above?
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_i][0] - xref_[lammps_j][0];
			double delY = xref_[lammps_i][1] - xref_[lammps_j][1];
			double delZ = xref_[lammps_i][2] - xref_[lammps_j][2];
			atomicHeatTable_(pair_index,0)=delX*fforce;
			atomicHeatTable_(pair_index,1)=delY*fforce;
			atomicHeatTable_(pair_index,2)=delZ*fforce;
			}
			else if (atomToElementMapType_ == EULERIAN) {
			atomicHeatTable_(pair_index,0)=delx*fforce;
			atomicHeatTable_(pair_index,1)=dely*fforce;
			atomicHeatTable_(pair_index,2)=delz*fforce;
			}
			}
			}

			//-------------------------------------------------------------------
			// on-the-fly calculation of stress
			void ATC_TransferHardy::compute_potential_stress(DENS_MAT& stress)
			{
			set_xPointer();
			stress.zero();
			// neighbor lists
			int *numneigh = lammpsInterface_->neighbor_list_numneigh();
			int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
			double ** xatom = lammpsInterface_->xatom();
			Array<bool> latticePeriodicity(3);
			latticePeriodicity(0) = (bool) periodicity[0];
			latticePeriodicity(1) = (bool) periodicity[1];
			latticePeriodicity(2) = (bool) periodicity[2];
			double lam1,lam2;
			double bond_value;
			// process differently for mesh vs translation-invariant kernels
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "ATC:: computing potential stress: " << flush;
			}
			int heartbeatFreq = (nNodes_ <= 10 ? 1 : (int) nNodes_ / 10);
			if (kernelFunction_ ) {
			// "normal" kernel functions
			DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
			double kernel_inv_vol = kernelFunction_->inv_vol();
			for (int i = 0; i < nNodes_; i++) {
			//for (int i = 0; i < nNodesGlobal_; i++) {
			if (lammpsInterface_->comm_rank() == 0 && i % heartbeatFreq == 0 ) {
			cout << "." << flush;
			}
			// Hardy point
			xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
			//xI = (feEngine_->get_feMesh())->global_coordinates(i);
			int inode = i;
			//int inode = feEngine_->get_feMesh()->map_global_to_unique(i);
			for (int j = 0; j < nLocal_; j++) {
			// second (neighbor) atom location
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			double * xj = xatom[lammps_j];
			// difference vector
			xaI = xa - xI;
			periodicity_correction(xaI);
			for (int k = 0; k < numneigh[lammps_j]; ++k) {
			int lammps_k = firstneigh[lammps_j][k];
			// first atom location
			xb.copy(xPointer_[lammps_k],3);
			double * xk = xatom[lammps_k];
			// difference vector
			xba = xb - xa;
			xbI = xba + xaI;
			kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
			// compute virial
			if (lam1 < lam2) {
			bond_value
			= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
			double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
			double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
			double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
			double rsq = delxdelx + delydely + delz*delz;
			double fforce = 0;
			lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
			fforce *= 0.5; // dbl count
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
			double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
			double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
			stress(inode,0) +=-delxfforcedelX*bond_value;
			stress(inode,1) +=-delxfforcedelY*bond_value;
			stress(inode,2) +=-delxfforcedelZ*bond_value;
			stress(inode,3) +=-delyfforcedelX*bond_value;
			stress(inode,4) +=-delyfforcedelY*bond_value;
			stress(inode,5) +=-delyfforcedelZ*bond_value;
			stress(inode,6) +=-delzfforcedelX*bond_value;
			stress(inode,7) +=-delzfforcedelY*bond_value;
			stress(inode,8) +=-delzfforcedelZ*bond_value;
			}
			else { //if (atomToElementMapType_ == EULERIAN) {
			stress(inode,0) +=-delxdelxfforce*bond_value;
			stress(inode,1) +=-delydelyfforce*bond_value;
			stress(inode,2) +=-delzdelzfforce*bond_value;
			stress(inode,3) +=-delxdelyfforce*bond_value;
			stress(inode,4) +=-delxdelzfforce*bond_value;
			stress(inode,5) +=-delydelzfforce*bond_value;
			}
			}
			}
			}
			}
			}
			else {
			// mesh-based kernel functions
			int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
			Array<int> node_list(nodes_per_element);
			DENS_VEC shp(nodes_per_element);
			DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
			for (int j = 0; j < nLocal_; j++) {
			if (lammpsInterface_->comm_rank() == 0 && j % heartbeatFreq == 0 ) {
			cout << "." << flush;
			}
			// first atom location
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			for (int k = 0; k < numneigh[lammps_j]; ++k) {
			int lammps_k = firstneigh[lammps_j][k];
			//if (lammps_k < lammps_j) continue; // full neighbor list
			// second (neighbor) atom location
			xb.copy(xPointer_[lammps_k],3);
			double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
			double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
			double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
			double rsq = delxdelx + delydely + delz*delz;
			double fforce = 0;
			lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
			fforce *= 0.5; // 1/2 sum_ab = sum_(ab)
			double virial[9];
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
			double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
			double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
			virial[0] =-delxfforcedelX;
			virial[1] =-delxfforcedelY;
			virial[2] =-delxfforcedelZ;
			virial[3] =-delyfforcedelX;
			virial[4] =-delyfforcedelY;
			virial[5] =-delyfforcedelZ;
			virial[6] =-delzfforcedelX;
			virial[7] =-delzfforcedelY;
			virial[8] =-delzfforcedelZ;
			}
			else {//if (atomToElementMapType_ == EULERIAN) {
			virial[0] =-delxdelxfforce;
			virial[1] =-delydelyfforce;
			virial[2] =-delzdelzfforce;
			virial[3] =-delxdelyfforce;
			virial[4] =-delxdelzfforce;
			virial[5] =-delydelzfforce;
			}
			lam1 = 0.0; lam2 = 1.0;
			double del_lambda = 0.5*(lam2 - lam1);
			double avg_lambda = 0.5*(lam2 + lam1);
			xab = xa - xb;
			for (int i = 0; i < line_ngauss; i++) {
			double lambda = del_lambda*line_xg[i] +avg_lambda;
			xlambda = lambda*xab + xb;
			// NOTE doesn't work for the case that fe_region \subset md_region
			int dummyEltID;
			feEngine_->shape_functions(xlambda,shp,dummyEltID,
			node_list,latticePeriodicity);
			// accumulate to nodes whose support overlaps the integration point
			for (int I = 0; I < nodes_per_element; I++) {
			// Use factor of 0.5 to account for line integration
			// domain of 0 to 1 instead of -1 to 1
			int inode = node_list(I);
			double bond_value
			= invNodeVolumes_(inode,inode)shp(I)line_wg[i]*0.5;
			stress(inode,0) += virial[0]*bond_value;
			stress(inode,1) += virial[1]*bond_value;
			stress(inode,2) += virial[2]*bond_value;
			stress(inode,3) += virial[3]*bond_value;
			stress(inode,4) += virial[4]*bond_value;
			stress(inode,5) += virial[5]*bond_value;
			if (atomToElementMapType_ == LAGRANGIAN) {
			stress(inode,6) += virial[6]*bond_value;
			stress(inode,7) += virial[7]*bond_value;
			stress(inode,8) += virial[8]*bond_value;
			}
			}
			}
			}
			}
			}// end if kernelFunction_
			if (lammpsInterface_->comm_rank() == 0) {
			cout << "done\n" << flush;
			}
			}

			//-------------------------------------------------------------------
			// on-the-fly calculation of the heat flux
			void ATC_TransferHardy::compute_potential_heatflux(DENS_MAT& flux)
			{
			set_xPointer();
			flux.zero();
			// neighbor lists
			int *numneigh = lammpsInterface_->neighbor_list_numneigh();
			int **firstneigh = lammpsInterface_->neighbor_list_firstneigh();
			double ** xatom = lammpsInterface_->xatom();
			Array<bool> latticePeriodicity(3);
			latticePeriodicity(0) = (bool) periodicity[0];
			latticePeriodicity(1) = (bool) periodicity[1];
			latticePeriodicity(2) = (bool) periodicity[2];
			double lam1,lam2;
			double bond_value;
			// process differently for mesh vs translation-invariant kernels
			if (kernelFunction_ ) {
			// "normal" kernel functions
			DENS_VEC xa(nsd_),xI(nsd_),xaI(nsd_),xb(nsd_),xbI(nsd_),xba(nsd_);
			double kernel_inv_vol = kernelFunction_->inv_vol();
			for (int i = 0; i < nNodes_; i++) {
			//for (int i = 0; i < nNodesGlobal_; i++) {
			int inode = i;
			//int inode = feEngine_->get_feMesh()->map_global_to_unique(i);
			// Hardy point
			xI = (feEngine_->get_feMesh())->nodal_coordinates(i);
			//xI = (feEngine_->get_feMesh())->global_coordinates(i);
			for (int j = 0; j < nLocal_; j++) {
			// second (neighbor) atom location
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			double * xj = xatom[lammps_j];
			// difference vector
			xaI = xa - xI;
			periodicity_correction(xaI);
			for (int k = 0; k < numneigh[lammps_j]; ++k) {
			int lammps_k = firstneigh[lammps_j][k];
			// first atom location
			xb.copy(xPointer_[lammps_k],3);
			double * xk = xatom[lammps_k];
			// difference vector
			xba = xb - xa;
			xbI = xba + xaI;
			kernelFunction_->bond_intercepts(xaI,xbI,lam1,lam2);
			// compute virial
			if (lam1 < lam2) {
			bond_value
			= kernel_inv_vol*(kernelFunction_->bond(xaI,xbI,lam1,lam2));
			double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
			double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
			double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
			double rsq = delxdelx + delydely + delz*delz;
			double fforce = 0;
			double pair_pe = lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
			fforce *= 0.5; // dbl count
			fforce = (delxuVariationVelocity_(j,0) +
			dely*uVariationVelocity_(j,1) +
			delz*uVariationVelocity_(j,2));
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
			double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
			double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
			flux(inode,0) +=fforcedelXbond_value;
			flux(inode,1) +=fforcedelYbond_value;
			flux(inode,2) +=fforcedelZbond_value;
			}
			else { //if (atomToElementMapType_ == EULERIAN) {
			flux(inode,0) +=fforcedelxbond_value;
			flux(inode,1) +=fforcedelybond_value;
			flux(inode,2) +=fforcedelzbond_value;
			}
			}
			}
			}
			}
			}
			else {
			// mesh-based kernel functions
			int nodes_per_element = feEngine_->get_feMesh()->get_nNodesPerElement();
			Array<int> node_list(nodes_per_element);
			DENS_VEC shp(nodes_per_element);
			DENS_VEC xa(nsd_),xb(nsd_),xab(nsd_),xlambda(nsd_);
			for (int j = 0; j < nLocal_; j++) {
			// first atom location
			int lammps_j = internalToAtom_(j);
			xa.copy(xPointer_[lammps_j],3);
			for (int k = 0; k < numneigh[lammps_j]; ++k) {
			int lammps_k = firstneigh[lammps_j][k];
			//if (lammps_k < lammps_j) continue; // full neighbor list
			// second (neighbor) atom location
			xb.copy(xPointer_[lammps_k],3);
			double delx = xatom[lammps_j][0] - xatom[lammps_k][0];
			double dely = xatom[lammps_j][1] - xatom[lammps_k][1];
			double delz = xatom[lammps_j][2] - xatom[lammps_k][2];
			double rsq = delxdelx + delydely + delz*delz;
			double fforce = 0;
			double pair_pe = lammpsInterface_->pair_force(lammps_j,lammps_k,rsq,fforce);
			fforce *= 0.5; // 1/2 sum_ab = sum_(ab)
			fforce = (delxuVariationVelocity_(j,0) +
			dely*uVariationVelocity_(j,1) +
			delz*uVariationVelocity_(j,2));
			double flux_vec[3];
			if (atomToElementMapType_ == LAGRANGIAN) {
			double delX = xref_[lammps_j][0] - xref_[lammps_k][0];
			double delY = xref_[lammps_j][1] - xref_[lammps_k][1];
			double delZ = xref_[lammps_j][2] - xref_[lammps_k][2];
			flux_vec[0] =fforce*delX;
			flux_vec[1] =fforce*delY;
			flux_vec[2] =fforce*delZ;
			}
			else {//if (atomToElementMapType_ == EULERIAN) {
			flux_vec[0] =fforce*delx;
			flux_vec[1] =fforce*dely;
			flux_vec[2] =fforce*delz;
			}
			lam1 = 0.0; lam2 = 1.0;
			double del_lambda = 0.5*(lam2 - lam1);
			double avg_lambda = 0.5*(lam2 + lam1);
			xab = xa - xb;
			for (int i = 0; i < line_ngauss; i++) {
			double lambda = del_lambda*line_xg[i] +avg_lambda;
			xlambda = lambda*xab + xb;
			// NOTE doesn't work for the case that fe_region \subset md_region
			int dummyEltID;
			feEngine_->shape_functions(xlambda,shp,dummyEltID,
			node_list,latticePeriodicity);
			// accumulate to nodes whose support overlaps the integration point
			for (int I = 0; I < nodes_per_element; I++) {
			// Use factor of 0.5 to account for line integration
			// domain of 0 to 1 instead of -1 to 1
			int inode = node_list(I);
			double bond_value
			= invNodeVolumes_(inode,inode)shp(I)line_wg[i]*0.5;
			flux(inode,0) += flux_vec[0]*bond_value;
			flux(inode,1) += flux_vec[1]*bond_value;
			flux(inode,2) += flux_vec[2]*bond_value;
			}
			}
			}
			}
			}// end if kernelFunction_
			}

			//-------------------------------------------------------------------
			void ATC_TransferHardy::compute_eshelby_stress(DENS_MAT & M,
			const DENS_MAT & E, const DENS_MAT & S, const DENS_MAT & H)
			{
			// eshelby stess:M, energy:E, stress:S, displacement gradient: H
			// eshelby stress = W I - F^T.P = W I - C.S [energy]
			// symmetric if isotropic S = a_0 I + a_1 C + a_2 C^2
			M.reset(nNodes_,fieldSizes_(HARDY_ESHELBY_STRESS));
			double nktv2p = lammpsInterface_->nktv2p();
			DENS_MAT P(3,3),FT(3,3),FTP(3,3),ESH(3,3);
			for (int i = 0; i < nNodes_; i++) {
			double W = E(i,0);
			ESH = 0.0;
			ESH(0,0) = W; ESH(1,1) = W; ESH(2,2) = W;
			// NOTE make matrix function to make this easier? see VoigtOp HACK
			// copy to local
			if (atomToElementMapType_ == LAGRANGIAN) {
			// Stress notation convention:: 0:11 1:12 2:13 3:21 4:22 5:23 6:31 7:32 8:33
			P(0,0) = S(i,0); P(0,1) = S(i,1); P(0,2) = S(i,2);
			P(1,0) = S(i,3); P(1,1) = S(i,4); P(1,2) = S(i,5);
			P(2,0) = S(i,6); P(2,1) = S(i,7); P(2,2) = S(i,8);
			// NOTE can just remove \int P N dA in post-processing
			//#define H_BASED
			#ifdef H_BASED
			FT(0,0) = H(i,0); FT(1,0) = H(i,1); FT(2,0) = H(i,2);
			FT(0,1) = H(i,3); FT(1,1) = H(i,4); FT(2,1) = H(i,5);
			FT(0,2) = H(i,6); FT(1,2) = H(i,7); FT(2,2) = H(i,8);
			#else
			FT(0,0) = 1+H(i,0); FT(1,0) = H(i,1); FT(2,0) = H(i,2);
			FT(0,1) = H(i,3); FT(1,1) = 1+H(i,4); FT(2,1) = H(i,5);
			FT(0,2) = H(i,6); FT(1,2) = H(i,7); FT(2,2) = 1+H(i,8);
			#endif
			}
			else if (atomToElementMapType_ == EULERIAN) {
			// Stress notation convention:: 0:11 1:22 2:33 3:12=21 4:13=31 5:23=32
			P(0,0) = S(i,0); P(0,1) = S(i,3); P(0,2) = S(i,4);
			P(1,0) = S(i,3); P(1,1) = S(i,1); P(1,2) = S(i,5);
			P(2,0) = S(i,4); P(2,1) = S(i,5); P(2,2) = S(i,2);
			FT(0,0) = H(i,0); FT(1,0) = H(i,1); FT(2,0) = H(i,2);
			FT(0,1) = H(i,3); FT(1,1) = H(i,4); FT(2,1) = H(i,5);
			FT(0,2) = H(i,6); FT(1,2) = H(i,7); FT(2,2) = H(i,8);
			}
			FTP = (1.0/nktv2p)FTP;
			ESH -= FTP;
			if (atomToElementMapType_ == EULERIAN) {
			// For Eulerian analysis, M = F^T*(w-H^T.CauchyStress)
			DENS_MAT Q(3,3);
			// Q stores (1-H)
			Q(0,0) = 1; Q(1,1) = 1; Q(2,2) = 1;
			Q(0,1) = 0; Q(1,0) = 0;
			Q(0,2) = 0; Q(2,0) = 0;
			Q(1,2) = 0; Q(2,1) = 0;
			Q -= FT.transpose();
			DENS_MAT F; F.reset(3,3);
			F = inv(Q);
			FT = F.transpose();
			ESH = FT*ESH;
			}
			// copy to global
			M(i,0) = ESH(0,0); M(i,1) = ESH(0,1); M(i,2) = ESH(0,2);
			M(i,3) = ESH(1,0); M(i,4) = ESH(1,1); M(i,5) = ESH(1,2);
			M(i,6) = ESH(2,0); M(i,7) = ESH(2,1); M(i,8) = ESH(2,2);
			#ifdef EXTENDED_ERROR_CHECKING
			if (atomToElementMapType_ == LAGRANGIAN) {
			cout << i << " W: " << W << "\n";
			#ifdef H_BASED
			FTP.print("H^T.P");
			#else
			FTP.print("F^T.P");
			#endif
			}
			else if (atomToElementMapType_ == EULERIAN) {
			cout << i << " w: " << W << "\n";
			FTP.print("H^T.CauchyStress");
			}
			ESH.print("Eshelby stress");
			#endif
			}
			}

			//-------------------------------------------------------------------
			// correction of node-pair distances due to periodicity
			void ATC_TransferHardy::periodicity_correction(DENS_VEC& xaI)
			{
			for (int m = 0; m < nsd_; m++) {
			xaI(m) -= periodicity[m]box_length[m]rnd(xaI(m)/box_length[m]);
			}
			}

			//-------------------------------------------------------------------
			// set xPointer_ to xref or xatom depending on Lagrangian/Eulerian analysis
			void ATC_TransferHardy::set_xPointer()
			{
			xPointer_ = xref_;
			if (atomToElementMapType_ == EULERIAN) {
			xPointer_ = lammpsInterface_->xatom();
			}
			}

			//-------------------------------------------------------------------
			// check consistency of fieldFlags_
			void ATC_TransferHardy::check_fieldFlags_consistency()
			{
			if (fieldFlags_(HARDY_TRANSFORMED_STRESS)) {
			fieldFlags_(HARDY_STRESS) = true;
			fieldFlags_(HARDY_DISPLACEMENT) = true;
			}
			if (fieldFlags_(HARDY_ESHELBY_STRESS)) {
			fieldFlags_(HARDY_STRESS) = true;
			fieldFlags_(HARDY_ENERGY) = true;
			fieldFlags_(HARDY_DISPLACEMENT) = true;
			}
			if (fieldFlags_(HARDY_CAUCHY_BORN_STRESS)) {
			if (! (cauchyBornStress_) ) {
			throw ATC_Error(0,"can't compute cauchy-born stress w/o cauchy born model");
			}
			fieldFlags_(HARDY_DISPLACEMENT) = true;
			}
			if (fieldFlags_(HARDY_VACANCY_CONCENTRATION)) {
			fieldFlags_(HARDY_DISPLACEMENT) = true;
			fieldFlags_(HARDY_NUMBER_DENSITY) = true;
			}
			if (fieldFlags_(HARDY_TEMPERATURE) \|\| fieldFlags_(HARDY_HEAT_FLUX) \|\|
			fieldFlags_(HARDY_ENERGY) \|\| (fieldFlags_(HARDY_STRESS) &&
			atomToElementMapType_ == EULERIAN) ) {
			fieldFlags_(HARDY_VELOCITY) = true;
			}
			if (fieldFlags_(HARDY_VELOCITY)) {
			fieldFlags_(HARDY_DENSITY) = true;
			fieldFlags_(HARDY_MOMENTUM) = true;
			}
			if (fieldFlags_(HARDY_DISPLACEMENT)) {
			fieldFlags_(HARDY_DENSITY) = true;
			}
			if (fieldFlags_(HARDY_TEMPERATURE) \|\|
			fieldFlags_(HARDY_TYPE_CONCENTRATION)) {
			fieldFlags_(HARDY_NUMBER_DENSITY) = true;
			}

			}

			} // end namespace ATC

lammps/lib/atc/ATC_TransferHardy.cppf70511be8768lammm-devel

ATC_TransferHardy.cppView Options

lammps/lib/atc/ATC_TransferHardy.cpp
f70511be8768lammm-devel

ATC_TransferHardy.cpp
View Options