Diffusion lammps (lammm-master-deprecated) Browse

FE_Engine.cpp
6354d994c21alammm-master-deprecated

FE_Engine.cpp
View Options

			#include "FE_Engine.h"
			#include "ATC_Transfer.h"
			#include "FE_Element.h"
			#include "XT_Function.h"
			#include "LammpsInterface.h"
			#include "PhysicsModel.h"

			#include <stdio.h>
			#include <map>

			using namespace std;

			namespace ATC{

			//-----------------------------------------------------------------
			FE_Engine::FE_Engine(ATC_Transfer * atcTransfer)
			: atcTransfer_(atcTransfer),
			feMesh_(NULL),
			initialized_(false),
			amendedMeshData_(false),
			outputManager_()
			{
			}

			FE_Engine::~FE_Engine()
			{
			if (feMesh_) delete feMesh_;
			if (connectivity_ && amendedMeshData_) delete connectivity_;
			if (nodeMap_ && amendedMeshData_) delete nodeMap_;
			if (coordinates_ && amendedMeshData_) delete coordinates_;
			}

			//-----------------------------------------------------------------
			void FE_Engine::initialize()
			{
			if (! feMesh_)
			throw ATC_Error(0,"FE_Engine has no mesh");

			if (! initialized_ ) {
			// mesh data
			nodeMap_ = &(feMesh_->global_to_unique_map());
			connectivity_ = &(feMesh_->connectivity());
			coordinates_ = &(feMesh_->nodal_coordinates());

			// set up work spaces
			nNodesPerElement_ = feMesh_->get_nNodesPerElement();
			nIPsPerElement_ = feMesh_->get_nIPsPerElement();
			nIPsPerFace_ = feMesh_->get_nIPsPerFace();
			nSD_ = feMesh_->get_nSpatialDimensions();
			nElems_ = feMesh_->get_nElements();

			// remove specified elements
			if (nullElements_.size() > 0) delete_elements(nullElements_);

			initialized_ = true;
			}
			}

			//-----------------------------------------------------------------
			bool FE_Engine::modify(int narg, char **arg)
			{
			bool match = false;
			/*! \page man_fem_mesh fix_modify AtC fem create mesh
			\section syntax
			fix_modify AtC fem create mesh <nx> <ny> <nz> <region-id>
			<f\|p> <f\|p> <f\|p>
			- nx ny nz = number of elements in x, y, z
			- region-id = id of region that is to be meshed
			- f p p = perioidicity flags for x, y, z
			\section examples
			<TT> fix_modify AtC fem create mesh 10 1 1 feRegion p p p </TT>
			\section description
			Creates a uniform mesh in a rectangular region
			\section restrictions
			creates only uniform rectangular grids in a rectangular region
			\section related
			\section default
			none
			*/
			if (strcmp(arg[0],"fem")==0) {
			// create mesh
			if (strcmp(arg[1],"create")==0) {
			if (strcmp(arg[2],"mesh")==0) {
			int nx = atoi(arg[3]);
			int ny = atoi(arg[4]);
			int nz = atoi(arg[5]);
			int xperiodic = 0;
			if (strcmp(arg[7],"p")==0) xperiodic = 1;
			int yperiodic = 0;
			if (strcmp(arg[8],"p")==0) yperiodic = 1;
			int zperiodic = 0;
			if (strcmp(arg[9],"p")==0) zperiodic = 1;
			if (feMesh_) throw ATC_Error(0,"FE Engine already has a mesh");
			create_mesh(nx, ny, nz, arg[6], xperiodic, yperiodic, zperiodic );

			// construct prescribed data manager NOTE move to ATC_Transfer?
			atcTransfer_->initialize_mesh_data();

			match = true;
			}
			}
			}
			/*! \page man_mesh_delete_elements fix_modify AtC mesh delete_elements
			\section syntax
			fix_modify AtC mesh delete_elements <element_set>
			- <element_set> = name of an element set
			\section examples
			<TT> fix_modify AtC delete_elements gap </TT>
			\section description
			Deletes a group of elements from the mesh.
			\section restrictions
			\section related
			\section default
			none
			*/
			else if ( (strcmp(arg[0],"mesh")==0)
			&& (strcmp(arg[1],"delete_elements")==0) ) {
			string esetName = arg[2];
			set<int> elemSet = feMesh_->get_elementset(esetName);
			nullElements_.insert(elemSet.begin(), elemSet.end());
			match = true;
			}
			// FE_Mesh
			else {
			match = feMesh_->modify(narg,arg);
			}
			return match;
			}

			//-----------------------------------------------------------------
			void FE_Engine::finish()
			{
			// Nothing to do
			}

			//-----------------------------------------------------------------
			// write geometry
			//-----------------------------------------------------------------
			void FE_Engine::initialize_output(int rank,
			string outputPrefix, OutputType otype)
			{
			outputManager_.initialize(outputPrefix, otype);
			if (! feMesh_) throw ATC_Error(0,"output needs mesh");
			if (! initialized_) initialize();
			if (rank == 0) outputManager_.write_geometry(*coordinates_, connectivity_);
			}

			//-----------------------------------------------------------------
			// write geometry
			//-----------------------------------------------------------------
			void FE_Engine::write_geometry(void)
			{
			outputManager_.write_geometry(*coordinates_, connectivity_);
			}

			// -------------------------------------------------------------
			// write data
			// -------------------------------------------------------------
			void FE_Engine::write_data(double time, FIELDS &soln, OUTPUT_LIST *data)
			{
			outputManager_.write_data(time, &soln, data, nodeMap_->get_data());
			}

			// -------------------------------------------------------------
			// write data
			// -------------------------------------------------------------
			void FE_Engine::write_data(double time, OUTPUT_LIST *data)
			{
			outputManager_.write_data(time, data, nodeMap_->get_data());
			}

			// -------------------------------------------------------------
			// amend mesh for deleted elements
			// -------------------------------------------------------------
			void FE_Engine::delete_elements(const set<int> & elementList)
			{
			int nsd = feMesh_->get_nSpatialDimensions();
			set<int> elementsNew;
			feMesh_->elementset_complement(elementList,elementsNew);
			int nElementsNew = elementsNew.size();
			set<int> newToOld;
			map<int,int> oldToNewMap;
			feMesh_->elementset_to_nodeset(elementsNew,newToOld);
			int nNodesNew = newToOld.size();
			set<int>::const_iterator itr;


			// coordinates & node map (from nodes to data)
			const DENS_MAT &coordinates = feMesh_->nodal_coordinates();
			const Array<int> & node_map = feMesh_->global_to_unique_map();
			if (coordinates_ && amendedMeshData_) delete coordinates_;
			coordinates_ = new DENS_MAT(nsd,nNodesNew);
			DENS_MAT & coor = * (const_cast<DENS_MAT *> ( coordinates_));
			if (nodeMap_ && amendedMeshData_) delete nodeMap_;
			nodeMap_ = new Array<int> (nNodesNew);
			Array<int> & nmap = * (const_cast<Array<int> *> ( nodeMap_));
			int k = 0, i = 0;
			for (itr = newToOld.begin(); itr != newToOld.end(); itr++) {
			int node = *itr;
			oldToNewMap[node]=i++;
			nmap(k) = node_map(node);
			for(int j = 0; j < nsd; j++) {
			coor(j,k) = coordinates(j,node);
			}
			k++;
			}
			// connectivity
			const Array2D<int> & connectivity = feMesh_->connectivity();
			int nNodesPerElement = connectivity.nRows();
			int nElements = connectivity.nCols();
			if (connectivity_ && amendedMeshData_) delete connectivity_;
			connectivity_ = new Array2D<int>(nNodesPerElement,nElementsNew);
			Array2D<int> & conn = * (const_cast<Array2D<int> *> ( connectivity_));
			k = 0;
			for (itr = elementsNew.begin(); itr != elementsNew.end(); itr++) {
			int ielem = *itr;
			for(int j = 0; j < nNodesPerElement; j++) {
			int old_node = connectivity(j,ielem);
			map<int,int>::iterator map_itr = oldToNewMap.find(old_node);
			if (map_itr == oldToNewMap.end()) {
			cout << "map failure " << old_node << "\n";
			}
			int node = map_itr->second;
			conn(j,k) = node;
			}
			k++;
			}
			amendedMeshData_ = true;
			}

			// -------------------------------------------------------------
			// compute dimensionless stiffness matrix using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_stiffness_matrix(SPAR_MAT &matrix) const
			{
			//int nfields = field_mask.get_length();
			int nNodes = feMesh_->get_nNodesUnique();
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();

			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			vector<DENS_MAT> dN(nsd);
			dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));
			DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);


			// assemble consistent mass (nnodes X nnodes)
			matrix.reset(nNodes,nNodes);
			DENS_MAT elementMassMatrix(nNodesPerElement,nNodesPerElement);
			for (int ielem = 0; ielem < nelems; ++ielem)
			{
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, dN, weights);
			// perform quadrature
			elementMassMatrix = dN[0].transMat(weights*dN[0]);
			for (int i = 1; i < nsd; ++i) {
			elementMassMatrix += dN[i].transMat(weights*dN[i]);
			}
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			int inode = conn(i);
			for (int j = 0; j < nNodesPerElement; ++j)
			{
			int jnode = conn(j);
			matrix.add(inode, jnode, elementMassMatrix(i,j));
			}
			}
			}
			}
			// -------------------------------------------------------------
			// compute dimensionless consistent mass using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_mass_matrix(SPAR_MAT &mass_matrix) const
			{
			int nNodes = feMesh_->get_nNodesUnique();
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();

			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// assemble consistent mass (nnodes X nnodes)
			mass_matrix.reset(nNodes,nNodes);
			DENS_MAT elementMassMatrix(nNodesPerElement,nNodesPerElement);
			for (int ielem = 0; ielem < nelems; ++ielem)
			{
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, weights);
			// perform quadrature
			elementMassMatrix = N.transMat(weights*N);
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			int inode = conn(i);
			for (int j = 0; j < nNodesPerElement; ++j)
			{
			int jnode = conn(j);
			mass_matrix.add(inode, jnode, elementMassMatrix(i,j));
			}
			}
			}
			}

			// -------------------------------------------------------------
			// compute consistent mass using given quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_mass_matrix(const DIAG_MAT &weights,
			const SPAR_MAT &N,
			SPAR_MAT &mass_matrix) const
			{
			int nn = N.nCols();
			int nips = N.nRows();

			DENS_MAT tmp_mass_matrix_local(nn,nn), tmp_mass_matrix(nn,nn);
			if (nips>0) { tmp_mass_matrix_local = N.transMat(weights*N); }
			// share information between processors
			int count = nn*nn;
			LammpsInterface::instance()->allsum(
			tmp_mass_matrix_local.get_ptr(),
			tmp_mass_matrix.get_ptr(), count);

			// create sparse from dense
			mass_matrix.reset(tmp_mass_matrix);

			}

			// -------------------------------------------------------------
			// compute dimensionless lumped mass using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_lumped_mass_matrix(DiagonalMatrix<double> & M) const
			{
			// initialize
			int nNodes = feMesh_->get_nNodesUnique();
			M.reset(nNodes,nNodes);

			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();

			DENS_MAT N(nNodesPerElement,nIPsPerElement);
			DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// assemble lumped diagonal mass (nnodes X nnodes)
			int inode = -1;
			DENS_VEC weightVec(nIPsPerElement);
			DENS_VEC Nvec(nNodesPerElement);
			for (int ielem = 0; ielem < nelems; ++ielem)
			{
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, weights);
			int nips = weights.nRows(); // different than nIPsPerElement?
			weightVec = 1.;
			weightVec = weights*weightVec;
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			Nvec = N*weightVec;
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			inode = conn(i);
			M(inode,inode) += Nvec(i);
			}
			}
			}

			// -------------------------------------------------------------
			// compute physical lumped mass using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_lumped_mass_matrix(
			const Array<FieldName>& field_mask,
			const FIELDS & fields,
			const PhysicsModel * physicsModel,
			const Array<int> &elementMaterials,
			map<FieldName, DIAG_MAT>& M, // mass matrix
			const Array<bool> *element_mask) const
			{
			int nfields = field_mask.get_length();
			// initialize
			int nNodes = feMesh_->get_nNodesUnique();
			for (int j = 0; j < nfields; ++j) {
			M[field_mask(j)].reset(nNodes,nNodes);
			}

			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();

			DENS_MAT N(nNodesPerElement,nIPsPerElement);
			DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);
			FIELDS fieldsAtIPs, localElementFields;
			DENS_MAT Nmat;

			// assemble lumped diagonal mass (nnodes X nnodes)
			for (int ielem = 0; ielem < nelems; ++ielem)
			{
			// if element is masked, skip it
			if (element_mask && !(*element_mask)(ielem)) continue;
			// material id
			int imat = elementMaterials(ielem);
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, weights);
			int nips = weights.nRows();
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);

			// interpolate fields at IPs
			FIELDS::const_iterator field;
			for (field = fields.begin(); field != fields.end(); field++)
			{
			// field values at all nodes
			const DENS_MAT &vI = field->second;
			// field values at element nodes
			DENS_MAT &vIe = localElementFields[field->first];
			// field values at integration points -> to be computed
			DENS_MAT &vP = fieldsAtIPs[field->first];

			int numFieldDOF = vI.nCols();
			vIe.reset(nNodesPerElement, numFieldDOF);
			// gather local field
			for (int i = 0; i < nNodesPerElement; i++) {
			for (int j = 0; j < numFieldDOF; j++) {
			vIe(i,j) = vI(conn(i),j);
			}
			}

			// interpolate field at integration points
			vP = N*vIe;
			}

			// compute energy/momentum/mass densities
			FIELDS capacities;
			physicsModel->M_integrand(field_mask, fieldsAtIPs, capacities, imat);

			// integrate & assemble
			for (int j = 0; j < nfields; ++j) {
			FieldName thisFieldName = field_mask(j);
			Nmat = N.transMat(weights*capacities[thisFieldName]);
			for (int i = 0; i < nNodesPerElement; ++i) {
			int inode = conn(i);
			M[thisFieldName](inode,inode) += Nmat(i,0);// assume all dof same
			}
			}
			}
			}

			// -------------------------------------------------------------
			// compute physical lumped mass using given quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_lumped_mass_matrix(
			const Array<FieldName> &field_mask,
			const FIELDS &fields,
			const PhysicsModel * physicsModel,
			const Array<set<int> > & pointMaterialGroups,
			const DIAG_MAT &weights, // NOTE use these as a mask
			const SPAR_MAT &N,
			map<FieldName, DIAG_MAT> &M) const // mass matrices
			{
			int nips = weights.nCols();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nfields = field_mask.get_length();
			int nNodes = feMesh_->get_nNodesUnique();
			// initialize
			FieldName thisFieldName;
			map<FieldName, DIAG_MAT> M_local;
			for (int j = 0; j < nfields; ++j) {
			thisFieldName = field_mask(j);
			M_local[thisFieldName].reset(nNodes,nNodes);
			M[thisFieldName].reset(nNodes,nNodes);
			}

			if (nips>0) {
			// compute fields at given ips
			FIELDS fieldsAtIPs;
			for (int j = 0; j < nfields; ++j) {
			thisFieldName = field_mask(j);
			const FIELD & field = (fields.find(thisFieldName))->second;
			int numFieldDOF = field.nCols();
			fieldsAtIPs[thisFieldName].reset(nips,numFieldDOF);
			fieldsAtIPs[thisFieldName] = N*field;
			}

			// treat single material point sets specially
			int nMatls = pointMaterialGroups.size();
			int atomMatls = 0;
			for (int imat = 0; imat < nMatls; imat++) {
			const set<int> & indices = pointMaterialGroups(imat);
			if ( indices.size() > 0) atomMatls++;
			}
			bool singleMaterial = ( atomMatls == 1 );

			// setup data structures
			FIELDS capacities;
			// evaluate physics model
			if (singleMaterial) {
			physicsModel->M_integrand(field_mask, fieldsAtIPs, capacities);
			}
			else {
			throw ATC_Error(0,"unimplemented function in FE_Engine::compute_mass_matrix");
			}

			// integrate & assemble
			for (int j = 0; j < nfields; ++j) {
			FieldName thisFieldName = field_mask(j);
			M_local[thisFieldName].reset( // assume all columns same
			column(N.transMat(weights*capacities[thisFieldName]),0) );
			}
			}

			// Share information between processors
			for (int j = 0; j < nfields; ++j) {
			FieldName thisFieldName = field_mask(j);
			int count = M_local[thisFieldName].size();
			LammpsInterface::instance()->allsum(
			M_local[thisFieldName].get_ptr(),
			M[thisFieldName].get_ptr(), count);
			}
			}

			//-----------------------------------------------------------------
			// compute assembled average gradient evaluated at the nodes
			//-----------------------------------------------------------------
			void FE_Engine::compute_gradient_matrix(
			GRAD_SHPFCN & grad_matrix) const
			{
			int nNodesUnique = feMesh_->get_nNodesUnique();
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();

			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			vector< DENS_MAT > dN(nsd);
			vector< DENS_MAT > tmp_grad_matrix(nsd);
			for (int i = 0; i < nsd; i++) {
			dN[i].reset(nNodesPerElement,nNodesPerElement);
			tmp_grad_matrix[i].reset(nNodesUnique,nNodesUnique);
			}

			DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			Array<int> count(nNodesUnique); count = 0;
			feMesh_->set_quadrature(FE_Element::NODAL_QUADRATURE);
			for (int ielem = 0; ielem < nelems; ++ielem) {
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, dN, weights);
			int nips = weights.nRows();
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);

			for (int j = 0; j < nIPsPerElement; ++j) {
			int jnode = conn(j); // NOTE assumes ip order is node order
			count(jnode) += 1;
			for (int i = 0; i < nNodesPerElement; ++i) {
			int inode = conn(i);
			for (int k = 0; k < nsd; ++k) {
			tmp_grad_matrix[k](jnode,inode) += dN[k](j,i);
			}
			}
			}
			}
			feMesh_->set_quadrature(FE_Element::GAUSSIAN_QUADRATURE); // reset to default
			for (int inode = 0; inode < nNodesUnique; ++inode) {
			//cout << inode << " count " << count(inode) << "\n";
			for (int jnode = 0; jnode < nNodesUnique; ++jnode) {
			for (int k = 0; k < nsd; ++k) {
			tmp_grad_matrix[k](jnode,inode) /= count(jnode);
			}
			}
			}
			// compact dense matrices
			grad_matrix.resize(nsd);
			for (int k = 0; k < nsd; ++k) {
			grad_matrix[k].reset(tmp_grad_matrix[k]);
			}
			}

			// -------------------------------------------------------------
			// compute energy per node using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_energy(const Array<FieldName> &mask,
			const FIELDS &fields,
			const PhysicsModel * physicsModel,
			const Array<int> & elementMaterials,
			FIELDS &energy,
			const Array<bool> *element_mask) const
			{
			//NOTE move to initialization and make workspace
			const int nNodesPerElement = feMesh_->get_nNodesPerElement();
			const int nIPsPerElement = feMesh_->get_nIPsPerElement();
			const int nsd = feMesh_->get_nSpatialDimensions();
			const int nelems = feMesh_->get_nElements();
			DENS_MAT Nmat;

			//NOTE move to initialization and make workspace
			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			vector<DENS_MAT> dN(nsd);
			dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));

			FIELDS fieldsAtIPs, localElementFields;
			GRAD_FIELDS grad_fieldsAtIPs;
			FIELDS::const_iterator field;
			FieldName thisFieldName;
			int numFieldDOF;

			for (int n = 0; n < mask.get_length(); n++) {
			FieldName thisFieldName = mask(n);
			const FIELD & field = (fields.find(thisFieldName))->second;
			energy[thisFieldName].reset(field.nRows(), field.nCols());
			}

			DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			int inode = -1;
			for (int ielem=0; ielem < nelems; ielem++)
			{
			// if element is masked, skip it
			if (element_mask && !(*element_mask)(ielem)) continue;
			// material id
			int imat = elementMaterials(ielem);
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, dN, weights);
			int nips = weights.nRows();
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			// compute fields and gradients of fields ips of this element

			for (int n = 0; n < mask.get_length(); n++) {
			FieldName thisFieldName = mask(n);
			FIELDS::const_iterator field = (fields.find(thisFieldName));
			// field values at all nodes
			const DENS_MAT &vI = field->second;
			// field values at element nodes
			DENS_MAT &vIe = localElementFields[field->first];
			// field values at integration points -> to be computed
			DENS_MAT &vP = fieldsAtIPs[field->first];
			// gradients of field at integration points -> to be computed
			vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];

			numFieldDOF = vI.nCols();
			vIe.reset(nNodesPerElement, numFieldDOF);
			// gather local field
			for (int i = 0; i < nNodesPerElement; i++)
			for (int j = 0; j < numFieldDOF; j++)
			vIe(i,j) = vI(conn(i),j);

			// interpolate field at integration points
			vP = N*vIe;
			dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
			for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
			}

			FIELDS energies;
			physicsModel->E_integrand(mask, fieldsAtIPs, grad_fieldsAtIPs, energies, imat);
			// assemble
			for (int n = 0; n < mask.get_length(); n++) {
			FieldName thisFieldName = mask(n);
			FIELDS::const_iterator field = (fields.find(thisFieldName));
			numFieldDOF = (field->second).nCols();
			Nmat = N.transMat(weights*energies[thisFieldName]);
			for (int i = 0; i < nNodesPerElement; ++i) {
			inode = conn(i);
			for (int k = 0; k < numFieldDOF; ++k) {
			energy[thisFieldName](inode,k) += Nmat(i,k);
			}
			}
			}
			} // element loop
			} // function call

			// -------------------------------------------------------------
			// compute rhs using native quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhs_mask,
			const FIELDS &fields,
			const PhysicsModel * physicsModel,
			const Array<int> & elementMaterials,
			FIELDS &rhs,
			const Array<bool> *element_mask) const
			{
			//NOTE move to initialization and make workspace
			const int nNodesPerElement = feMesh_->get_nNodesPerElement();
			const int nIPsPerElement = feMesh_->get_nIPsPerElement();
			const int nsd = feMesh_->get_nSpatialDimensions();
			const int nelems = feMesh_->get_nElements();
			DENS_MAT Nmat;

			//NOTE move to initialization and make workspace
			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			vector<DENS_MAT> dN(nsd);
			dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));

			FIELDS fieldsAtIPs, localElementFields;
			GRAD_FIELDS grad_fieldsAtIPs;
			FIELDS::const_iterator field;
			FieldName thisFieldName;
			int numFieldDOF;

			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX) \|\| rhs_mask(thisFieldName,SOURCE)) {
			int nrows = (field->second).nRows();
			int ncols = (field->second).nCols();
			rhs[thisFieldName].reset(nrows,ncols);
			}
			}

			DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			int inode = -1;
			for (int ielem=0; ielem < nelems; ielem++)
			{
			// if element is masked, skip it
			if (element_mask && !(*element_mask)(ielem)) continue;
			// material id
			int imat = elementMaterials(ielem);
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, dN, weights);
			int nips = weights.nRows();
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			// compute fields and gradients of fields ips of this element

			for (field = fields.begin(); field != fields.end(); field++)
			{
			// field values at all nodes
			const DENS_MAT &vI = field->second;
			// field values at element nodes
			DENS_MAT &vIe = localElementFields[field->first];
			// field values at integration points -> to be computed
			DENS_MAT &vP = fieldsAtIPs[field->first];
			// gradients of field at integration points -> to be computed
			vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];

			numFieldDOF = vI.nCols();
			vIe.reset(nNodesPerElement, numFieldDOF);
			// gather local field
			for (int i = 0; i < nNodesPerElement; i++)
			for (int j = 0; j < numFieldDOF; j++)
			vIe(i,j) = vI(conn(i),j);

			// interpolate field at integration points
			vP = N*vIe;
			dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
			for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
			}

			// NOTE move scaling of N by weights up here

			// evaluate physics model
			// N_fluxes is a set of fields
			if(physicsModel->has_N_integrand()) {
			FIELDS N_fluxes;
			physicsModel->N_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, N_fluxes, imat);
			// assemble
			for (field = fields.begin(); field != fields.end(); field++)
			{
			thisFieldName = field->first;
			// NOTE why is this line here?
			// SOURCE should refer only to adding extrinsic coupling sources
			if (!rhs_mask(thisFieldName,SOURCE)) continue;
			numFieldDOF = (field->second).nCols();
			Nmat = N.transMat(weights*N_fluxes[thisFieldName]);
			for (int i = 0; i < nNodesPerElement; ++i) {
			inode = conn(i);
			for (int k = 0; k < numFieldDOF; ++k) {
			rhs[thisFieldName](inode,k) += Nmat(i,k);
			}
			}
			}
			}

			// evaluate Physics model
			// B_fluxes is a set of field gradients
			if (physicsModel->has_B_integrand()) {
			GRAD_FIELDS B_fluxes;
			physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes, imat);

			// assemble
			for (field = fields.begin(); field != fields.end(); field++)
			{
			thisFieldName = field->first;
			numFieldDOF = (field->second).nCols();
			if (!rhs_mask(thisFieldName,FLUX)) continue;

			for (int j = 0; j < nsd; j++)
			{
			Nmat = dN[j].transMat(weights*B_fluxes[thisFieldName][j]);
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			inode = conn(i);
			for (int k = 0; k < numFieldDOF; ++k)
			{
			rhs[thisFieldName](inode,k) += Nmat(i,k);
			} // scatter k
			} // scatter i
			} // loop over nsd
			} //loop on fields
			} // if B_integrand
			} // element loop
			} // function call

			// -------------------------------------------------------------
			// compute rhs using given quadrature
			// -------------------------------------------------------------
			void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhs_mask,
			const FIELDS &fields,
			const PhysicsModel *physicsModel,
			const Array<set<int> >&pointMaterialGroups,
			const DIAG_MAT &weights,
			const SPAR_MAT &N,
			const GRAD_SHPFCN &dN,
			FIELDS &rhs) const
			{
			int nips = weights.nCols();
			int nNodesUnique = feMesh_->get_nNodesUnique();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();

			FieldName thisFieldName;
			FIELDS::const_iterator field;
			FIELDS rhs_local;
			for (field = rhs.begin(); field != rhs.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX) \|\| rhs_mask(thisFieldName,SOURCE)) {
			int nrows = (field->second).nRows();
			int ncols = (field->second).nCols();
			rhs [thisFieldName].reset(nrows,ncols);
			rhs_local[thisFieldName].reset(nrows,ncols);
			}
			}

			if (nips>0) {
			// compute fields and gradients of fields at given ips
			GRAD_FIELDS grad_fieldsAtIPs;
			FIELDS fieldsAtIPs;
			int numFieldDOF;
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			const FIELD & field = (fields.find(thisFieldName))->second;
			numFieldDOF = field.nCols();
			grad_fieldsAtIPs[thisFieldName].assign(nsd,DENS_MAT(nips,numFieldDOF));
			fieldsAtIPs[thisFieldName].reset(nips,numFieldDOF);
			fieldsAtIPs[thisFieldName] = N*field;
			for (int j = 0; j < nsd; ++j) {
			grad_fieldsAtIPs[thisFieldName][j] = dN[j]*field;
			}
			}

			// treat single material point sets specially
			int nMatls = pointMaterialGroups.size();
			int atomMatls = 0;
			for (int imat = 0; imat < nMatls; imat++) {
			const set<int> & indices = pointMaterialGroups(imat);
			if ( indices.size() > 0) atomMatls++;
			}
			bool singleMaterial = ( atomMatls == 1 );

			// setup data structures
			FIELDS N_fluxes;
			GRAD_FIELDS B_fluxes;
			if(physicsModel->has_B_integrand()) {
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			if ( rhs_mask(index,FLUX) ) {
			B_fluxes[thisFieldName].reserve(nsd);
			}
			}
			}
			// evaluate physics model
			if (singleMaterial)
			{
			if(physicsModel->has_N_integrand()) {
			physicsModel->N_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs,
			N_fluxes);
			}
			if(physicsModel->has_B_integrand()) {
			physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs,
			B_fluxes);
			}
			}
			else
			{
			// NOTE this copy in fields/copy out fluxes is extremely slow
			// NOTE alternates:
			// 1) permanent workspace with per-row mapped clones per matl
			// from caller/atc
			// 2) per point calls to N/B_integrand
			// 3) collect internalToAtom into materials and use mem-cont clones
			// what about dof for matrices and data storage: clone _matrix_

			// for each material group:
			// set up storage
			FIELDS groupN_fluxes, fieldsGroup;
			GRAD_FIELDS groupB_fluxes, grad_fieldsGroup;
			if(physicsModel->has_B_integrand()) {
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			int ndof = (field->second).nCols();
			grad_fieldsGroup[thisFieldName].assign(nsd,DENS_MAT(nips,ndof));
			N_fluxes[thisFieldName].reset(nips,ndof);
			B_fluxes[thisFieldName].assign(nsd,DENS_MAT(nips,ndof));
			if ( rhs_mask(index,FLUX) ) {
			groupB_fluxes[thisFieldName].reserve(nsd);
			}
			}
			}
			// copy fields
			for ( int imat = 0; imat < pointMaterialGroups.get_length(); imat++)
			{
			const set<int> & indices = pointMaterialGroups(imat);
			int npts = indices.size();
			int i = 0;
			for (set<int>::const_iterator iter=indices.begin();
			iter != indices.end(); iter++, i++ ) {
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int ndof = (field->second).nCols();
			fieldsGroup[thisFieldName].reset(npts,ndof);
			for (int j = 0; j < nsd; ++j) {
			(grad_fieldsGroup[thisFieldName][j]).reset(npts,ndof);
			}
			for (int dof = 0; dof < ndof; ++dof) {
			fieldsGroup[thisFieldName](i,dof)
			= fieldsAtIPs[thisFieldName](*iter,dof);
			for (int j = 0; j < nsd; ++j) {
			grad_fieldsGroup[thisFieldName][j](i,dof)
			= grad_fieldsAtIPs[thisFieldName][j](*iter,dof);
			}
			}
			}
			}
			// calculate forces, & assemble
			if(physicsModel->has_N_integrand()) {
			physicsModel->N_integrand(rhs_mask, fieldsGroup, grad_fieldsGroup,
			groupN_fluxes, imat);
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			int ndof = (field->second).nCols();
			if ( rhs_mask(index,SOURCE) ) {
			int i = 0;
			for (set<int>::const_iterator iter=indices.begin();
			iter != indices.end(); iter++, i++ ) {
			for (int dof = 0; dof < ndof; ++dof) {
			N_fluxes[thisFieldName](*iter,dof)
			= groupN_fluxes[thisFieldName](i,dof);
			}
			}
			}
			}
			}
			// calculate forces, & assemble
			if(physicsModel->has_B_integrand()) {
			physicsModel->B_integrand(rhs_mask, fieldsGroup, grad_fieldsGroup,
			groupB_fluxes, imat);
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			int ndof = (field->second).nCols();
			if ( rhs_mask(index,FLUX) ) {
			int i = 0;
			for (set<int>::const_iterator iter=indices.begin();
			iter != indices.end(); iter++, i++ ) {
			for (int dof = 0; dof < ndof; ++dof) {
			for (int j = 0; j < nsd; ++j) {
			B_fluxes[thisFieldName][j](*iter,dof)
			= groupB_fluxes[thisFieldName][j](i,dof);
			}
			}
			}
			}
			}
			}
			}
			}

			// assemble
			for (field = fields.begin(); field != fields.end(); field++)
			{
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			if(physicsModel->has_N_integrand()) {
			if ( rhs_mask(index,SOURCE) ) {
			rhs_local[thisFieldName]
			+= N.transMat(weights*N_fluxes[thisFieldName]);
			}
			}
			if(physicsModel->has_B_integrand()) {
			if ( rhs_mask(index,FLUX) ) {
			for (int j = 0; j < nsd; ++j) {
			rhs_local[thisFieldName]
			+= dN[j].transMat(weights*B_fluxes[thisFieldName][j]);
			}
			}
			}
			}
			}

			// Share information between processors
			for (field = rhs.begin(); field != rhs.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX) \|\| rhs_mask(thisFieldName,SOURCE)) {
			int count = rhs_local[thisFieldName].size();
			//rhs_local[thisFieldName].print("RHS LCL");
			LammpsInterface::instance()->allsum(
			rhs_local[thisFieldName].get_ptr(),
			rhs[thisFieldName].get_ptr(), count);
			}
			}
			}

			// -------------------------------------------------------------
			// compute flux for post processing
			// -------------------------------------------------------------
			void FE_Engine::compute_flux(const Array2D<bool> &rhs_mask,
			const FIELDS &fields,
			const PhysicsModel * physicsModel,
			const Array<int> & elementMaterials,
			GRAD_FIELDS &flux,
			const Array<bool> *element_mask) const
			{
			//NOTE move to initialization and make workspace
			const int nNodesPerElement = feMesh_->get_nNodesPerElement();
			const int nIPsPerElement = feMesh_->get_nIPsPerElement();
			const int nsd = feMesh_->get_nSpatialDimensions();
			const int nelems = feMesh_->get_nElements();
			DENS_MAT Nmat;

			//NOTE move to initialization and make workspace
			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			vector<DENS_MAT> dN(nsd);
			dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));

			FIELDS fieldsAtIPs, localElementFields;
			GRAD_FIELDS grad_fieldsAtIPs;
			FIELDS::const_iterator field;
			FieldName thisFieldName;
			int numFieldDOF;

			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX)) {
			int nrows = (field->second).nRows();
			int ncols = (field->second).nCols();
			flux[thisFieldName].assign(nsd, DENS_MAT(nrows,ncols));
			}
			}

			DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			int inode = -1;
			for (int ielem=0; ielem < nelems; ielem++)
			{
			// if element is masked, skip it
			if (element_mask && !(*element_mask)(ielem)) continue;
			// material id
			int imat = elementMaterials(ielem);
			// evaluate shape functions
			feMesh_->shape_function(ielem, N, dN, weights);
			int nips = weights.nRows();
			// get connectivity
			feMesh_->element_connectivity_unique(ielem, conn);
			// compute fields and gradients of fields ips of this element

			for (field = fields.begin(); field != fields.end(); field++)
			{
			// field values at all nodes
			const DENS_MAT &vI = field->second;
			// field values at element nodes
			DENS_MAT &vIe = localElementFields[field->first];
			// field values at integration points -> to be computed
			DENS_MAT &vP = fieldsAtIPs[field->first];
			// gradients of field at integration points -> to be computed
			vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];

			numFieldDOF = vI.nCols();
			vIe.reset(nNodesPerElement, numFieldDOF);
			// gather local field
			for (int i = 0; i < nNodesPerElement; i++)
			for (int j = 0; j < numFieldDOF; j++)
			vIe(i,j) = vI(conn(i),j);

			// interpolate field at integration points
			vP = N*vIe;
			dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
			for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
			}

			// NOTE move scaling of N by weights up here

			// evaluate Physics model
			// B_fluxes is a set of field gradients
			if (physicsModel->has_B_integrand()) {
			GRAD_FIELDS B_fluxes;
			physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes, imat);

			// assemble
			for (field = fields.begin(); field != fields.end(); field++)
			{
			thisFieldName = field->first;
			numFieldDOF = (field->second).nCols();
			if (!rhs_mask(thisFieldName,FLUX)) continue;

			for (int j = 0; j < nsd; j++)
			{
			Nmat = N.transMat(weights*B_fluxes[thisFieldName][j]);
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			inode = conn(i);
			for (int k = 0; k < numFieldDOF; ++k)
			{
			flux[thisFieldName][j](inode,k) += Nmat(i,k);
			} // scatter k
			} // scatter i
			} // loop over nsd
			} //loop on fields
			} // if B_integrand
			} // element loop
			} // function call

			//-----------------------------------------------------------------
			// boundary flux using native quadrature
			//-----------------------------------------------------------------
			void FE_Engine::compute_boundary_flux(
			const Array2D<bool> & rhs_mask,
			const FIELDS & fields,
			const PhysicsModel * physicsModel,
			const Array<int> & elementMaterials,
			const set< pair <int,int> > & faceSet,
			FIELDS & rhs) const
			{
			//NOTE move to initialization and make workspace
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerFace = feMesh_->get_nIPsPerFace();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();
			DENS_MAT Nmat(nNodesPerElement,3);

			FieldName thisFieldName;
			FIELDS::const_iterator field;
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX)) {
			int nrows = (field->second).nRows();
			int ncols = (field->second).nCols();
			rhs [thisFieldName].reset(nrows,ncols);
			}
			}

			//NOTE move to initialization and make workspace
			// sizing working arrays
			DENS_MAT N(nIPsPerFace,nNodesPerElement);
			vector< DENS_MAT > dN(nsd), nN(nsd);
			for (int i = 0; i < nsd; i++) {
			dN[i].reset(nIPsPerFace,nNodesPerElement);
			nN[i].reset(nIPsPerFace,nNodesPerElement);
			}

			GRAD_FIELDS grad_fieldsAtIPs;
			FIELDS fieldsAtIPs;
			FIELDS localElementFields;
			int numFieldDOF;
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			numFieldDOF = field->second.nCols();
			grad_fieldsAtIPs[thisFieldName].reserve(nsd);
			for (int i = 0; i < nsd; ++i) {
			grad_fieldsAtIPs[thisFieldName].push_back(DENS_MAT(nIPsPerFace,numFieldDOF));
			}
			fieldsAtIPs[thisFieldName].reset(nIPsPerFace,numFieldDOF);
			localElementFields[thisFieldName].reset(nNodesPerElement,numFieldDOF);
			}

			DiagonalMatrix<double> weights(nIPsPerFace,nIPsPerFace);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			int inode = -1;
			set< pair <int,int> >::iterator iter;
			for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
			// evaluate shape functions at ips
			feMesh_->face_shape_function(*iter, N, dN, nN, weights);
			int nips = weights.nRows();
			// get connectivity
			int ielem = iter->first;
			int imat = elementMaterials(ielem);
			feMesh_->element_connectivity_unique(ielem, conn);

			// interpolate fields and gradients of fields ips of this element
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			const DenseMatrix<double> * thisField
			= (const DENS_MAT *) &(field->second);
			numFieldDOF = thisField->nCols();
			//NOTE an hopefully unnecessary copy (should alias)
			for (int i = 0; i < nNodesPerElement; i++) {
			for (int j = 0; j < numFieldDOF; j++) {
			localElementFields[thisFieldName](i,j) = (*thisField)(conn(i),j);
			}
			}
			// ips X dof = ips X nodes * nodes X dof
			fieldsAtIPs[thisFieldName] = N*localElementFields[thisFieldName];
			for (int j = 0; j < nsd; ++j) {
			grad_fieldsAtIPs[thisFieldName][j] = dN[j]*localElementFields[thisFieldName];
			}
			}

			// Evaluate- physics model
			// do nothing for N_integrand
			// nN : precomputed and held by ATC_Transfer
			if(physicsModel->has_B_integrand()) {
			map<FieldName, vector< DENS_MAT > > B_fluxes;
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			if ( rhs_mask(thisFieldName,FLUX) ) {
			B_fluxes[thisFieldName].reserve(nsd);
			}
			}
			physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes,imat);
			// assemble
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			int index = (int) thisFieldName;
			const DenseMatrix<double> * thisField
			= (const DENS_MAT *) &(field->second);
			if ( rhs_mask(index,FLUX) ) {
			numFieldDOF = thisField->nCols();
			Nmat.reset(nNodesPerElement,numFieldDOF);
			for (int j = 0; j < nsd; j++) {
			// nodes X dof = nodes X ips * ips X dof ???
			Nmat = nN[j].transMat(weights*B_fluxes[thisFieldName][j]);
			for (int i = 0; i < nNodesPerElement; ++i) {
			inode = conn(i);
			for (int k = 0; k < numFieldDOF; ++k) {
			rhs[thisFieldName](inode,k) += Nmat(i,k);
			}
			}
			}
			}
			}
			}

			}
			}

			// -------------------------------------------------------------
			// compute boundary flux using given quadrature and interpolation
			// -------------------------------------------------------------
			void FE_Engine::compute_boundary_flux(
			const Array2D<bool> & rhs_mask,
			const FIELDS & fields,
			const PhysicsModel * physicsModel,
			const Array< int > & elementMaterials,
			const Array< set<int> > & pointMaterialGroups,
			const DIAG_MAT & weights,
			const SPAR_MAT & N,
			const GRAD_SHPFCN & dN,
			const DIAG_MAT & flux_mask,
			FIELDS & flux) const
			{
			int nNodesUnique = feMesh_->get_nNodesUnique();
			FieldName thisFieldName;

			map<FieldName, DENS_MAT > rhs;
			map<FieldName, DENS_MAT > rhsSubset;
			map<FieldName, DENS_MAT > rhsSubset_local;

			FIELDS::const_iterator field;
			for (field = fields.begin(); field != fields.end(); field++) {
			thisFieldName = field->first;
			if (rhs_mask(thisFieldName,FLUX)) {
			int nrows = (field->second).nRows();
			int ncols = (field->second).nCols();
			rhs [thisFieldName].reset(nrows,ncols);
			rhsSubset[thisFieldName].reset(nrows,ncols);
			}
			}

			// R_I = - int_Omega Delta N_I .q dV
			compute_rhs_vector(rhs_mask, fields, physicsModel, elementMaterials, rhs);
			// R_I^md = - int_Omega^md Delta N_I .q dV
			compute_rhs_vector(rhs_mask, fields, physicsModel, pointMaterialGroups,
			weights, N, dN, rhsSubset);

			// flux_I = 1/Delta V_I V_I^md R_I + R_I^md
			// where : Delta V_I = int_Omega N_I dV
			for (int n = 0; n < NUM_FIELDS; ++n) {
			FieldName thisFieldName = (FieldName) n;
			if ( rhs_mask(thisFieldName,FLUX) ) {
			flux[thisFieldName]
			= rhsSubset[thisFieldName] - flux_mask*rhs[thisFieldName];
			}
			}
			}

			//-----------------------------------------------------------------
			// prescribed boundary flux using native quadrature
			// integrate \int_delV N_I s(x,t).n dA
			//-----------------------------------------------------------------
			void FE_Engine::add_fluxes(const Array<bool> &fieldMask,
			const double time,
			const SURFACE_SOURCE & sourceFunctions,
			FIELDS &nodalSources) const
			{
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nNodesPerFace = feMesh_->get_nNodesPerFace();
			int nIPsPerFace = feMesh_->get_nIPsPerFace();
			int nsd = feMesh_->get_nSpatialDimensions();
			FieldName thisField;
			int nFieldDOF = 1;
			XT_Function * f = NULL;

			//NOTE move arrays to initialization and make workspace
			// sizing working arrays
			vector<DENS_MAT> dN, nN;
			dN.assign(nsd, DENS_MAT(nIPsPerFace, nNodesPerElement));
			nN.assign(nsd, DENS_MAT(nIPsPerFace, nNodesPerElement));

			DENS_MAT N(nIPsPerFace,nNodesPerElement);
			DENS_MAT xCoords(nsd,nNodesPerElement);
			DENS_MAT xAtIPs(nsd,nIPsPerFace);
			DENS_MAT Nmat(nNodesPerElement,nsd);
			DIAG_MAT weights(nIPsPerFace,nIPsPerFace);
			Array<int> conn(nNodesPerElement);
			DENS_MAT faceSource(nIPsPerFace,nFieldDOF);

			// element wise assembly
			SURFACE_SOURCE::const_iterator src_iter;
			for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
			{
			FieldName thisFieldName = src_iter->first;
			if (!fieldMask((int)thisFieldName)) continue;

			typedef map<PAIR,Array<XT_Function*> > FSET;
			const FSET fset = (const FSET )&(src_iter->second);
			FSET::const_iterator fset_iter;
			for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
			{
			const PAIR &face = fset_iter->first;
			const int elem = face.first;
			const Array <XT_Function*> &fs = fset_iter->second;
			feMesh_->element_connectivity_unique(elem, conn);
			// evaluate location at ips
			feMesh_->face_shape_function(face, N, dN, nN, weights);
			feMesh_->element_coordinates(elem, xCoords);
			MultAB(xCoords,N,xAtIPs,0,1); //xAtIPs = xCoords*(N.transpose());

			// interpolate prescribed flux at ips of this element
			// NOTE this seems totally redundant
			FSET::const_iterator face_iter = fset->find(face);
			if (face_iter == fset->end())
			{
			cout << "face not found" << std::endl;
			// NOTE: do something there
			}
			// NOTE why can't I know the size earlier??
			// NOTE only normal flux here
			nFieldDOF = (face_iter->second).size();
			faceSource.reset(nIPsPerFace,nFieldDOF);
			for (int ip = 0; ip < nIPsPerFace; ++ip)
			{
			for (int idof = 0; idof<nFieldDOF; ++idof)
			{
			f = fs(idof);
			if (!f) continue;
			faceSource(ip,idof) = f->f(column(xAtIPs,ip).get_ptr(),time);
			}
			}
			// assemble
			Nmat = N.transMat(weights*faceSource);
			for (int i = 0; i < nNodesPerElement; ++i)
			{
			int inode = conn(i);
			for (int idof = 0; idof < nFieldDOF; ++idof)
			{
			nodalSources[thisFieldName](inode,idof) += Nmat(i,idof);
			} // end assemble nFieldDOF
			} // end assemble nNodesPerElement
			} // end fset loop
			} // field loop
			}

			//-----------------------------------------------------------------
			// prescribed volume flux using native quadrature
			// integrate \int_V N_I s(x,t) dV
			//-----------------------------------------------------------------
			void FE_Engine::add_sources(const Array<bool> &fieldMask,
			const double time,
			const VOLUME_SOURCE &sourceFunctions,
			FIELDS &nodalSources) const
			{
			int nNodes = feMesh_->get_nNodes();
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerElement = feMesh_->get_nIPsPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();
			int nFieldDOF = 1;
			XT_Function * f = NULL;

			//NOTE move arrays to initialization and make workspace
			// sizing working arrays
			DENS_MAT elemSource(nIPsPerElement,nFieldDOF);
			DENS_MAT N(nIPsPerElement,nNodesPerElement);
			DENS_MAT xCoords(nsd,nNodesPerElement);
			DENS_MAT xAtIPs(nsd,nIPsPerElement);
			DENS_MAT Nwg(nNodesPerElement,1);
			DENS_MAT Nmat(nNodesPerElement,nsd);
			DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
			Array<int> conn(nNodesPerElement);

			// element wise assembly
			for (int ielem = 0; ielem < nelems; ++ielem) {
			feMesh_->element_connectivity_unique(ielem, conn);

			// evaluate location at ips
			feMesh_->shape_function(ielem, N, weights);
			feMesh_->element_coordinates(ielem, xCoords);
			xAtIPs =xCoords*(N.transpose());

			VOLUME_SOURCE ::const_iterator src_iter;
			for (src_iter = sourceFunctions.begin();
			src_iter != sourceFunctions.end(); src_iter++) {
			FieldName thisField = src_iter->first;
			int index = (int) thisField;
			if ( fieldMask(index) ) {
			const Array2D<XT_Function > thisSource
			= (const Array2D<XT_Function > ) &(src_iter->second);
			nFieldDOF = thisSource->nCols();
			elemSource.reset(nIPsPerElement,nFieldDOF);

			// interpolate prescribed flux at ips of this element
			for (int ip = 0; ip < nIPsPerElement; ++ip) {
			for (int idof = 0; idof < nFieldDOF; ++idof) {
			f = (*thisSource)(ielem,idof);
			if (f) {
			elemSource(ip,idof) = f->f(column(xAtIPs,ip).get_ptr(),time);
			}
			}
			}

			// assemble
			Nmat.reset(nNodesPerElement,nFieldDOF);
			Nmat = N.transMat(weights*elemSource);
			for (int i = 0; i < nNodesPerElement; ++i) {
			int inode = conn(i);
			for (int idof = 0; idof < nFieldDOF; ++idof) {
			nodalSources[thisField](inode,idof) += Nmat(i,idof);
			}
			}
			}
			}
			}
			}

			//-----------------------------------------------------------------
			// boundary integral of a nodal field
			//-----------------------------------------------------------------
			void FE_Engine::field_surface_flux(
			const DENS_MAT & field,
			const set< PAIR > & faceSet,
			DENS_MAT & values,
			const bool contour,
			const int axis) const
			{
			int nNodesPerElement = feMesh_->get_nNodesPerElement();
			int nIPsPerFace = feMesh_->get_nIPsPerFace();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nelems = feMesh_->get_nElements();
			int numFieldDOF = field.nCols();

			double a[3] = {0,0,0};
			a[axis] = 1;

			// sizing working arrays
			DENS_MAT N(nIPsPerFace,nNodesPerElement);
			DENS_MAT n(nsd,nIPsPerFace);
			DENS_MAT fieldsAtIPs(nIPsPerFace,numFieldDOF);
			DENS_MAT localElementFields(nNodesPerElement,numFieldDOF);
			DiagonalMatrix<double> weights(nIPsPerFace,nIPsPerFace);
			Array<int> conn(nNodesPerElement);

			DENS_MAT Nmat(nsd,numFieldDOF);
			DENS_MAT integrals(numFieldDOF,nsd);

			// element wise assembly
			set< pair <int,int> >::iterator iter;
			for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
			// evaluate shape functions at ips
			feMesh_->face_shape_function(*iter, N, n, weights);
			// cross n with axis to get tangent
			if (contour) {
			double t[3];
			for (int i = 0; i < nIPsPerFace; i++) {
			t[0] = a[1]n(2,i) - a[2]n(1,i);
			t[1] = a[2]n(0,i) - a[0]n(2,i);
			t[2] = a[0]n(1,i) - a[1]n(0,i);
			n(0,i) = t[0];
			n(1,i) = t[1];
			n(2,i) = t[2];
			}
			}
			int nips = weights.nRows();
			// get connectivity
			int ielem = iter->first;
			feMesh_->element_connectivity_unique(ielem, conn);

			// interpolate fields and gradients of fields ips of this element
			for (int i = 0; i < nNodesPerElement; i++) {
			for (int j = 0; j < numFieldDOF; j++) {
			localElementFields(i,j) = field(conn(i),j);
			}
			}
			// ips X dof = ips X nodes * nodes X dof
			fieldsAtIPs = N*localElementFields;

			// assemble : integral(k,j) = sum_ip n(j,ip) wg(ip,ip) field(ip,k)
			Nmat = nweightsfieldsAtIPs;
			for (int j = 0; j < nsd; j++) {
			for (int k = 0; k < numFieldDOF; ++k) {
			integrals(k,j) += Nmat(j,k);
			}
			}
			}
			// (S.n)_1 = S_1j n_j = S_11 n_1 + S_12 n_2 + S_13 n_3
			// (S.n)_2 = S_2j n_j = S_21 n_1 + S_22 n_2 + S_23 n_3
			// (S.n)_3 = S_3j n_j = S_31 n_1 + S_32 n_2 + S_33 n_3
			if (numFieldDOF == 9) { // tensor
			values.reset(nsd,1);
			values(0,0) = integrals(0,0)+integrals(1,1)+integrals(2,2);
			values(1,0) = integrals(3,0)+integrals(4,1)+integrals(5,2);
			values(2,0) = integrals(6,0)+integrals(7,1)+integrals(8,2);
			}
			else if (numFieldDOF == 6) { // sym tensor
			values.reset(nsd,1);
			values(0,0) = integrals(0,0)+integrals(3,1)+integrals(4,2);
			values(1,0) = integrals(3,0)+integrals(1,1)+integrals(5,2);
			values(2,0) = integrals(4,1)+integrals(5,1)+integrals(2,2);
			}
			// (v.n) = v_j n_j
			else if (numFieldDOF == 3) { // vector
			values.reset(1,1);
			values(0,0) = integrals(0,0)+integrals(1,1)+integrals(2,2);
			}
			// s n = s n_j e_j
			else if (numFieldDOF == 1) { // scalar
			values.reset(nsd,1);
			values(0,0) = integrals(0,0);
			values(1,0) = integrals(0,1);
			values(2,0) = integrals(0,2);
			}
			else {
			string msg = "FE_Engine::field_surface_flux unsupported field width: ";
			msg += tostring(numFieldDOF);
			throw ATC_Error(0,msg);
			}
			}


			//-----------------------------------------------------------------
			// evaluate shape functions at given points
			//-----------------------------------------------------------------
			void FE_Engine::evaluate_shape_functions(
			const MATRIX & pt_coords,
			SPAR_MAT & N,
			Array<int> & pointToEltMap) const
			{
			// Get shape function and derivatives at atomic locations
			int npe = feMesh_->get_nNodesPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nnodes = feMesh_->get_nNodesUnique();
			int npts = pt_coords.nCols();

			pointToEltMap.reset(npts);

			// loop over point (atom) coordinates
			DENS_VEC x(nsd);
			Array<int> node_index(npe);
			DENS_VEC shp(npe);

			N.reset(npts,nnodes);
			for (int i = 0; i < npts; ++i) {
			for (int k = 0; k < nsd; ++k) {
			x(k) = pt_coords(k,i);
			}
			int eltID;
			feMesh_->shape_functions(x,shp,eltID,node_index);
			for (int j = 0; j < npe; ++j) {
			int jnode = node_index(j);
			N.add(i,jnode,shp(j));
			}
			pointToEltMap(i) = eltID;
			}
			N.compress();
			}

			//-----------------------------------------------------------------
			// evaluate shape functions and their spatial derivatives at given points
			//-----------------------------------------------------------------
			void FE_Engine::evaluate_shape_functions(
			const MATRIX & pt_coords,
			SPAR_MAT & N,
			GRAD_SHPFCN & dN,
			Array<int> & pointToEltMap) const
			{
			// Get shape function and derivatives at atomic locations
			int npe = feMesh_->get_nNodesPerElement();
			int nsd = feMesh_->get_nSpatialDimensions();
			int nnodes = feMesh_->get_nNodesUnique();
			int npts = pt_coords.nCols();

			pointToEltMap.reset(npts);

			// loop over point (atom) coordinates
			DENS_VEC x(nsd);
			Array<int> node_index(npe);
			DENS_VEC shp(npe);
			DENS_MAT dshp(nsd,npe);

			for (int k = 0; k < nsd; ++k) {
			dN[k].reset(npts,nnodes);
			}

			N.reset(npts,nnodes);
			for (int i = 0; i < npts; ++i) {
			for (int k = 0; k < nsd; ++k) {
			x(k) = pt_coords(k,i);
			}
			int eltID;
			feMesh_->shape_functions(x,shp,dshp,eltID,node_index);
			for (int j = 0; j < npe; ++j) {
			int jnode = node_index(j);
			N.add(i,jnode,shp(j));
			for (int k = 0; k < nsd; ++k) {
			dN[k].add(i,jnode,dshp(k,j));
			}
			}
			pointToEltMap(i) = eltID;
			}
			N.compress();
			for (int k = 0; k < nsd; ++k) {
			dN[k].compress();
			}
			}

			//-----------------------------------------------------------------
			// evaluate shape functions at given points
			//-----------------------------------------------------------------
			void FE_Engine::evaluate_shape_functions(const VECTOR & x,
			Array<int>& node_index,
			DENS_VEC& shp,
			DENS_MAT& dshp,
			int & eltID) const
			{
			// Get shape function and derivatives at a specific point
			int nsd = feMesh_->get_nSpatialDimensions();
			int npe = feMesh_->get_nNodesPerElement();

			dshp.reset(nsd,npe);

			feMesh_->shape_functions(x,shp,dshp,eltID,node_index);

			}

			//-----------------------------------------------------------------
			// create a uniform rectangular mesh on a rectangular region
			//-----------------------------------------------------------------
			void FE_Engine::create_mesh(int nx, int ny, int nz, char * regionName,
			int xperiodic, int yperiodic, int zperiodic)
			{
			double xmin, xmax, ymin, ymax, zmin, zmax;
			double xscale, yscale, zscale;
			double boxxlo, boxxhi, boxylo, boxyhi, boxzlo, boxzhi;

			// check to see if region exists and is it a box, and if so get the bounds
			bool isBox;
			isBox = atcTransfer_->get_region_bounds(regionName,
			xmin, xmax,
			ymin, ymax,
			zmin, zmax,
			xscale,
			yscale,
			zscale);
			if (!isBox) throw ATC_Error(0,"Region for FE mesh is not a box");

			feMesh_ = new FE_Uniform3DMesh(nx,ny,nz,
			xmin,xmax,
			ymin,ymax,
			zmin,zmax,
			xscale,
			yscale,
			zscale,
			xperiodic,
			yperiodic,
			zperiodic);
			if (LammpsInterface::instance()->comm_rank()==0) {
			printf(" ATC:: created FEM Mesh with %i Global Nodes, %i Unique Nodes, and %i Elements\n",
			feMesh_->get_nNodes(),feMesh_->get_nNodesUnique(),feMesh_->get_nElements()); }
			}

			};

lammps/lib/atc/FE_Engine.cpp6354d994c21alammm-master-deprecated

FE_Engine.cppView Options

lammps/lib/atc/FE_Engine.cpp
6354d994c21alammm-master-deprecated

FE_Engine.cpp
View Options