#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 = n*weights*fieldsAtIPs;
      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()); }
  }

};