void
DiracKernel::updateCaches(const Elem* old_elem,
                          const Elem* new_elem,
                          Point p,
                          unsigned id)
{
  // Update the point cache.  Remove old cached data, only cache
  // new_elem if it is non-NULL and local.
  _point_cache.erase(id);
  if (new_elem && (new_elem->processor_id() == processor_id()))
    _point_cache[id] = std::make_pair(new_elem, p);

  // Update the reverse cache
  //
  // First, remove the Point from the old_elem's vector
  reverse_cache_t::iterator it = _reverse_point_cache.find(old_elem);
  if (it != _reverse_point_cache.end())
  {
    reverse_cache_t::mapped_type & points = it->second;
    {
      reverse_cache_t::mapped_type::iterator
        points_it = points.begin(),
        points_end = points.end();

      for (; points_it != points_end; ++points_it)
      {
        // If the point matches, remove it from the vector of points
        if (p.relative_fuzzy_equals(points_it->first))
        {
          // Vector erasure.  It can be slow but these vectors are
          // generally very short.  It also invalidates existing
          // iterators, so we need to break out of the loop and
          // not use them any more.
          points.erase(points_it);
          break;
        }
      }

      // If the points vector is now empty, remove old_elem from the reverse cache entirely
      if (points.empty())
        _reverse_point_cache.erase(old_elem);
    }
  }

  // Next, if new_elem is not NULL and local, add the point to the new_elem's vector
  if (new_elem && (new_elem->processor_id() == processor_id()))
  {
    reverse_cache_t::mapped_type & points = _reverse_point_cache[new_elem];
    points.push_back(std::make_pair(p, id));
  }
}

void
ContactSlipDamper::timestepSetup()
{
  GeometricSearchData & displaced_geom_search_data = _displaced_problem->geomSearchData();
  std::map<std::pair<unsigned int, unsigned int>, PenetrationLocator *> * penetration_locators = &displaced_geom_search_data._penetration_locators;

  for (pl_iterator plit = penetration_locators->begin(); plit != penetration_locators->end(); ++plit)
  {
    PenetrationLocator & pen_loc = *plit->second;

    if (operateOnThisInteraction(pen_loc))
    {
      std::vector<dof_id_type> & slave_nodes = pen_loc._nearest_node._slave_nodes;

      for (unsigned int i = 0; i < slave_nodes.size(); i++)
      {
        dof_id_type slave_node_num = slave_nodes[i];

        if (pen_loc._penetration_info[slave_node_num])
        {
          PenetrationInfo & info = *pen_loc._penetration_info[slave_node_num];
          const Node * node = info._node;

          if (node->processor_id() == processor_id())
//              && info.isCaptured())                  //TODO maybe just set this everywhere?
            info._slip_reversed = false;
        }
      }
    }
  }
}

void
NodalNormalsEvaluator::execute()
{
  if (_current_node->processor_id() == processor_id())
  {
    if (_current_node->n_dofs(_aux.number(), _fe_problem.getVariable(_tid, "nodal_normal_x").number()) > 0)
    {
      Threads::spin_mutex::scoped_lock lock(nodal_normals_evaluator_mutex);

      dof_id_type dof_x = _current_node->dof_number(_aux.number(), _fe_problem.getVariable(_tid, "nodal_normal_x").number(), 0);
      dof_id_type dof_y = _current_node->dof_number(_aux.number(), _fe_problem.getVariable(_tid, "nodal_normal_y").number(), 0);
      dof_id_type dof_z = _current_node->dof_number(_aux.number(), _fe_problem.getVariable(_tid, "nodal_normal_z").number(), 0);

      NumericVector<Number> & sln = _aux.solution();
      Real nx = sln(dof_x);
      Real ny = sln(dof_y);
      Real nz = sln(dof_z);

      Real n = std::sqrt((nx * nx) + (ny * ny) + (nz * nz));
      if (std::abs(n) >= 1e-13)
      {
        // divide by n only if it is not close to zero to avoid NaNs
        sln.set(dof_x, nx / n);
        sln.set(dof_y, ny / n);
        sln.set(dof_z, nz / n);
      }
    }
  }
}

Real
FindValueOnLine::getValueAtPoint(const Point & p)
{
  const Elem * elem = (*_pl)(p);

  if (elem)
  {
    Real value;

    if (elem->processor_id() == processor_id())
    {
      // element is local
      _point_vec[0] = p;
      _subproblem.reinitElemPhys(elem, _point_vec, 0);
      value = _coupled_var->sln()[0];
    }

    // broadcast value
    _communicator.broadcast(value, elem->processor_id());
    return value;
  }
  else
  {
    // there is no element
    mooseError("No element found at the current search point. Please make sure the sampling line stays inside the mesh completely.");
  }
}

void
MaterialTensorOnLine::finalize()
{

    _communicator.set_union(_dist);
    _communicator.set_union(_value);

   if (processor_id() == 0)
   {
     std::vector<std::pair<Real, Real> > table;
//     std::map<unsigned int, Real>::iterator id(_dist.begin());
//     std::map<unsigned int, Real>::iterator is(_value.begin());
     std::map<std::pair<unsigned int, unsigned int>, Real>::iterator id(_dist.begin());
     std::map<std::pair<unsigned int, unsigned int>, Real>::iterator is(_value.begin());

     while (id != _dist.end() && is != _value.end())
     {
       table.push_back(std::make_pair(id->second,is->second));
       id++;
       is++;
     }

     std::sort(table.begin(),table.end());

     _output_file << "time " << _fe_problem.time() << std::endl;
     for (std::vector<std::pair<Real, Real> >::iterator it = table.begin();
          it != table.end();
          ++it)
     {
       _output_file << it->first << "  " << it->second << std::endl;
     }
   }

}

ExodusTimeSequenceStepper::ExodusTimeSequenceStepper(const InputParameters & parameters) :
    TimeSequenceStepperBase(parameters),
    _mesh_file(getParam<MeshFileName>("mesh"))
{
  // Read the Exodus file on processor 0
  std::vector<Real> times;
  if (processor_id() == 0)
  {
    // Check that the required file exists
    MooseUtils::checkFileReadable(_mesh_file);

    // dummy mesh
    ReplicatedMesh mesh(_communicator);

    ExodusII_IO exodusII_io(mesh);
    exodusII_io.read(_mesh_file);
    times = exodusII_io.get_time_steps();
  }

  // distribute timestep list
  unsigned int num_steps = times.size();
  _communicator.broadcast(num_steps);
  times.resize(num_steps);
  _communicator.broadcast(times);

  setupSequence(times);
}

void
CSV::output()
{
  // Call the base class output (populates tables)
  TableOutput::output();

  // Print the table containing all the data to a file
  if (!_all_data_table.empty() && processor_id() == 0)
    _all_data_table.printCSV(filename(), 1, _align);

  // Output each VectorPostprocessor's data to a file
  for (std::map<std::string, FormattedTable>::iterator it = _vector_postprocessor_tables.begin(); it != _vector_postprocessor_tables.end(); ++it)
  {
    std::ostringstream output;
    output << _file_base << "_" << it->first;

    output << "_"
           << std::setw(_padding)
           << std::setprecision(0)
           << std::setfill('0')
           << std::right
           << timeStep();

    output << ".csv";

    if (_set_delimiter)
      it->second.setDelimiter(_delimiter);
    it->second.setPrecision(_precision);
    it->second.printCSV(output.str(), 1, _align);
  }
}

Real
ExceptionKernel::computeQpResidual()
{
  // We need a thread lock here so that we don't introduce a race condition when inspecting or
  // changing the static variable
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mutex);

  if (_when == WhenType::INITIAL_CONDITION)
    throw MooseException("MooseException thrown during initial condition computation");

  // Make sure we have called computeQpResidual enough times to
  // guarantee that we are in the middle of a linear solve, to verify
  // that we can throw an exception at that point.
  // Also, only throw on one rank if the user requests it
  else if (_when == WhenType::RESIDUAL && !_res_has_thrown && time_to_throw() &&
           (_rank == DofObject::invalid_processor_id || _rank == processor_id()))
  {
    // The residual has now thrown
    _res_has_thrown = true;

    if (_should_throw)
      throw MooseException("MooseException thrown during residual calculation");
    else
      mooseError("Intentional error triggered during residual calculation");
  }
  else
    return 0.;
}

void
SlopeLimitingBase::deserialize(std::vector<std::string> & serialized_buffers)
{
  // The input string stream used for deserialization
  std::istringstream iss;

  mooseAssert(serialized_buffers.size() == _app.n_processors(),
              "Unexpected size of serialized_buffers: " << serialized_buffers.size());

  for (auto rank = decltype(_app.n_processors())(0); rank < serialized_buffers.size(); ++rank)
  {
    if (rank == processor_id())
      continue;

    iss.str(serialized_buffers[rank]); // populate the stream with a new buffer
    iss.clear();                       // reset the string stream state

    // Load the communicated data into all of the other processors' slots

    unsigned int size = 0;
    iss.read((char *)&size, sizeof(size));

    for (unsigned int i = 0; i < size; i++)
    {
      dof_id_type key;
      loadHelper(iss, key, this);

      std::vector<RealGradient> value;
      loadHelper(iss, value, this);

      // merge the data we received from other procs
      _lslope.insert(std::pair<dof_id_type, std::vector<RealGradient>>(key, value));
    }
  }
}

void
FeatureFloodCount::visitNodalNeighbors(const Node * node, int current_idx, FeatureData * feature, bool recurse)
{
  mooseAssert(node, "Node is NULL");

  std::vector<const Node *> neighbors;
  MeshTools::find_nodal_neighbors(_mesh.getMesh(), *node, _nodes_to_elem_map, neighbors);

  // Loop over all nodal neighbors
  for (unsigned int i = 0; i < neighbors.size(); ++i)
  {
    const Node * neighbor_node = neighbors[i];
    processor_id_type my_proc_id = processor_id();

    // Only recurse on nodes this processor can see
    if (_mesh.isSemiLocal(const_cast<Node *>(neighbor_node)))
    {
      if (neighbor_node->processor_id() != my_proc_id)
        feature->_ghosted_ids.insert(neighbor_node->id());

      // Premark Halo values
      feature->_halo_ids.insert(neighbor_node->id());

      if (recurse)
        flood(neighbors[i], current_idx, feature);
    }
  }
}

void
PointSamplerBase::finalize()
{
  // Save off for speed
  unsigned int pid = processor_id();

  /*
   * Figure out which processor is actually going "claim" each point.
   * If multiple processors found the point and computed values what happens is that
   * maxloc will give us the smallest PID in max_id
   */
  std::vector<unsigned int> max_id(_found_points.size());

  _communicator.maxloc(_found_points, max_id);

  for (unsigned int i=0; i<max_id.size(); i++)
  {
    // Only do this check on the proc zero because it's the same on every processor
    // _found_points should contain all 1's at this point (ie every point was found by a proc)
    if (pid == 0 && !_found_points[i])
      mooseError("In " << name() << ", sample point not found: " << _points[i]);

    if (max_id[i] == pid)
      SamplerBase::addSample(_points[i], _ids[i], _values[i]);
  }

  SamplerBase::finalize();
}

Real
SolutionTimeAdaptiveDT::computeDT()
{
    //Ratio grew so switch direction
    if (_sol_time_vs_dt > _old_sol_time_vs_dt && _sol_time_vs_dt > _older_sol_time_vs_dt)
    {
        _direction *= -1;

        // Make sure we take at least two steps in this new direction
        _old_sol_time_vs_dt = std::numeric_limits<Real>::max();
        _older_sol_time_vs_dt = std::numeric_limits<Real>::max();
    }

//  if (_t_step > 1)
    Real local_dt =  _dt + _dt * _percent_change*_direction;

    if ((_adapt_log) && (processor_id() == 0))
    {
        Real out_dt = getCurrentDT();
        if (out_dt > _dt_max)
            out_dt = _dt_max;

        _adaptive_log<<"***Time step: "<<_t_step<<", time = "<<_time+out_dt<<std::endl;
        _adaptive_log<<"Cur DT: "<<out_dt<<std::endl;
        _adaptive_log<<"Older Ratio: "<<_older_sol_time_vs_dt<<std::endl;
        _adaptive_log<<"Old Ratio: "<<_old_sol_time_vs_dt<<std::endl;
        _adaptive_log<<"New Ratio: "<<_sol_time_vs_dt<<std::endl;
    }

    return local_dt;
}

void MeshASMPartitioning::DoPartition( const unsigned *block_size, vector < vector< unsigned > > &block_elements,
                                         vector <unsigned> &block_type_range){

    unsigned iproc = processor_id();
    unsigned ElemOffset    = _mesh._elementOffset[iproc];
    unsigned ElemOffsetp1  = _mesh._elementOffset[iproc + 1];
    unsigned OwnedElements = ElemOffsetp1 - ElemOffset;

    unsigned counter[3] = {0, 0, 0};
    unsigned flag_block[3] = {4, 3, 2};
    
    for (unsigned iel = ElemOffset; iel < ElemOffsetp1; iel++) {
      unsigned flag_mat   = _mesh.GetElementMaterial(iel);

      if (flag_mat == flag_block[0]) {
        counter[0]++;
      }
      else if (flag_mat == flag_block[1]) {
        counter[1]++;
      }
    }
    counter[2] = OwnedElements - counter[0] - counter[1];

    block_type_range.resize(3);

  

    unsigned block_start = 0;
    unsigned iMaterial = 0;
    while (iMaterial < 3) {
      if (counter[iMaterial] != 0 ) { //material of this type is there
        unsigned reminder = counter[iMaterial] % block_size[iMaterial];
        unsigned blocks = (0 == reminder) ? counter[iMaterial] / block_size[iMaterial] : counter[iMaterial] / block_size[iMaterial] + 1 ;
        block_elements.resize(block_start + blocks);

        for (int i = 0; i < blocks; i++) {
          block_elements[block_start + i].resize(block_size[iMaterial]);
        }
        if  (0 != reminder ) {
          block_elements[block_start + (blocks - 1u)].resize(reminder);
        }

        unsigned counter = 0;
        for (unsigned iel = ElemOffset; iel < ElemOffsetp1; iel++) {
          unsigned flag_mat   = _mesh.GetElementMaterial(iel);
          if ( flag_block[iMaterial] == flag_mat ) {
            block_elements[ block_start + (counter / block_size[iMaterial]) ][ counter % block_size[iMaterial] ] = iel;
            counter++;
          }
        }
        block_type_range[iMaterial] = block_start + blocks;
        block_start += blocks;
      } 
      else {
        block_type_range[iMaterial] = block_start;
      }
      iMaterial++;
    }
  }

void GCTaskThread::run() {
  // Set up the thread for stack overflow support
  this->record_stack_base_and_size();
  this->initialize_thread_local_storage();
  // Bind yourself to your processor.
  if (processor_id() != GCTaskManager::sentinel_worker()) {
    if (TraceGCTaskThread) {
      tty->print_cr("GCTaskThread::run: "
                    "  binding to processor %u", processor_id());
    }
    if (!os::bind_to_processor(processor_id())) {
      DEBUG_ONLY(
        warning("Couldn't bind GCTaskThread %u to processor %u",
                      which(), processor_id());
      )
    }
  }

void
SlaveConstraint::addPoints()
{
    _point_to_info.clear();

    std::set<dof_id_type> & has_penetrated = _penetration_locator._has_penetrated;

    std::map<dof_id_type, PenetrationInfo *>::iterator
    it  = _penetration_locator._penetration_info.begin(),
    end = _penetration_locator._penetration_info.end();
    for (; it!=end; ++it)
    {
        PenetrationInfo * pinfo = it->second;

        if (!pinfo)
        {
            continue;
        }

        dof_id_type slave_node_num = it->first;

        const Node * node = pinfo->_node;

        std::set<dof_id_type>::iterator hpit = has_penetrated.find(slave_node_num);
        if (hpit != has_penetrated.end() && node->processor_id() == processor_id())
        {
            // Find an element that is connected to this node that and that is also on this processor

            std::vector<dof_id_type> & connected_elems = _mesh.nodeToElemMap()[slave_node_num];

            Elem * elem = NULL;

            for (unsigned int i=0; i<connected_elems.size() && !elem; ++i)
            {
                Elem * cur_elem = _mesh.elem(connected_elems[i]);
                if (cur_elem->processor_id() == processor_id())
                    elem = cur_elem;
            }

            mooseAssert(elem, "Couldn't find an element on this processor that is attached to the slave node!");

            addPoint(elem, *node);
            _point_to_info[*node] = pinfo;
        }
    }
}

MaterialTensorOnLine::~MaterialTensorOnLine()
{
  if (_stream_open && processor_id() == 0)
  {
    _output_file.close();
    _stream_open = false;
  }
}

void
SlaveConstraint::addPoints()
{
  _point_to_info.clear();

  std::map<dof_id_type, PenetrationInfo *>::iterator
      it = _penetration_locator._penetration_info.begin(),
      end = _penetration_locator._penetration_info.end();

  const auto & node_to_elem_map = _mesh.nodeToElemMap();
  for (; it != end; ++it)
  {
    PenetrationInfo * pinfo = it->second;

    // Skip this pinfo if there are no DOFs on this node.
    if (!pinfo || pinfo->_node->n_comp(_sys.number(), _vars[_component]) < 1)
      continue;

    dof_id_type slave_node_num = it->first;
    const Node * node = pinfo->_node;

    if (pinfo->isCaptured() && node->processor_id() == processor_id())
    {
      // Find an element that is connected to this node that and that is also on this processor
      auto node_to_elem_pair = node_to_elem_map.find(slave_node_num);
      mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing node in node to elem map");
      const std::vector<dof_id_type> & connected_elems = node_to_elem_pair->second;

      Elem * elem = NULL;

      for (unsigned int i = 0; i < connected_elems.size() && !elem; ++i)
      {
        Elem * cur_elem = _mesh.elemPtr(connected_elems[i]);
        if (cur_elem->processor_id() == processor_id())
          elem = cur_elem;
      }

      mooseAssert(elem,
                  "Couldn't find an element on this processor that is attached to the slave node!");

      addPoint(elem, *node);
      _point_to_info[*node] = pinfo;
    }
  }
}

void
SplitTester::execute()
{
  if (!_fe_problem.mesh().isDistributedMesh())
    mooseError("Not using DistributedMesh but should be");
  auto & m = _fe_problem.mesh().getMesh();
  if (m.n_elem_on_proc(processor_id()) == m.n_elem())
    mooseError("Elements not shared properly between distributed mesh procs");
}

//----------------------------------------------------------------------------------------------------------------
void MeshASMPartitioning::DoPartition( const unsigned *block_size, vector < vector< unsigned > > &block_elements,
					 vector <unsigned> &block_type_range){

  unsigned iproc=processor_id();
  unsigned ElemOffset    = _mesh.IS_Mts2Gmt_elem_offset[iproc];
  unsigned ElemOffsetp1  = _mesh.IS_Mts2Gmt_elem_offset[iproc+1];
  unsigned OwnedElements = ElemOffsetp1 - ElemOffset;

  unsigned counter[2]={0,0};
  for (unsigned iel_mts = ElemOffset; iel_mts < ElemOffsetp1; iel_mts++) {
    unsigned kel        = _mesh.IS_Mts2Gmt_elem[iel_mts]; 
    unsigned flag_mat   = _mesh.el->GetElementMaterial(kel);

    if(2 == flag_mat){
      counter[1]++;
    } 
  }
  counter[0]=OwnedElements - counter[1];
  
  block_type_range.resize(2);
 
  unsigned flag_block[2]={4,2};
  
  unsigned block_start=0;
  unsigned iblock=0;
  while(iblock < 2){
    if(counter[iblock] !=0 ){
      unsigned reminder = counter[iblock] % block_size[iblock];
      unsigned blocks = (0 == reminder)? counter[iblock]/block_size[iblock] : counter[iblock]/block_size[iblock] + 1 ;
      block_elements.resize(block_start+blocks);
  
      for(int i = 0; i < blocks; i++) 
	block_elements[block_start + i].resize(block_size[iblock]);
      if  (0 != reminder ){
	block_elements[block_start + (blocks-1u)].resize(reminder); 
      }
      
      unsigned counter=0;
      for (unsigned iel_mts = ElemOffset; iel_mts < ElemOffsetp1; iel_mts++) {
	unsigned kel        = _mesh.IS_Mts2Gmt_elem[iel_mts]; 
 	unsigned flag_mat   = _mesh.el->GetElementMaterial(kel);
	if( flag_block[iblock] == flag_mat ){
	  block_elements[ block_start + (counter / block_size[iblock]) ][ counter % block_size[iblock] ]=iel_mts;
	  counter++;
	}	 
      }
      block_type_range[iblock]=block_start+blocks;
      block_start += blocks;
    }
    else{
      block_type_range[iblock]=block_start;
    }
    iblock++;
  }

} 

void
CSV::output()
{
  // Call the base class output (populates tables)
  TableOutput::output();

  // Print the table containing all the data to a file
  if (!_all_data_table.empty() && processor_id() == 0)
    _all_data_table.printCSV(filename());
}