[Libmesh-users] Augmenting Sparsity to add Nonlocal DOFs

[Charlie T. <cha...@gm...>]

Hi, I modified the files in miscellaneous example 9 to augment the sparsity
pattern to include elements found using Patch, and I can't seem to get this
to work, what am I missing? Thanks!

#ifndef AUGMENT_SPARSITY_NONLOCAL_H
#define AUGMENT_SPARSITY_NONLOCAL_H

#include "libmesh/ghosting_functor.h"
#include "libmesh/mesh_base.h"
#include "libmesh/patch.h"
using libMesh::Elem;
using libMesh::GhostingFunctor;
using libMesh::MeshBase;
using libMesh::boundary_id_type;
using libMesh::processor_id_type;

// Convenient typedef for a map for (element, side id) --> element neighbor
typedef std::map<std::pair<const Elem *, unsigned char>, const Elem *>
ElementSideMap;

// And the inverse map, but without sides in the key
typedef std::map<const Elem *, const Elem *> ElementMap;

class AugmentSparsityNonlocal : public GhostingFunctor
{
private:

  /**
   * The Mesh we're calculating on
   */
  MeshBase & _mesh;

  /**
   * A map from (lower element ID, side ID) to matching upper element ID.
Here "lower"
   * and "upper" refer to the lower and upper (wrt +z direction) sides of
the crack in
   * our mesh.
   */
  ElementSideMap _lower_to_upper;

  /**
   * The inverse (ignoring sides) of the above map.
   */
  ElementMap _upper_to_lower;

  /**
   * Boundary IDs for the lower and upper faces of the "crack" in the mesh.
   */
  boundary_id_type _crack_boundary_lower, _crack_boundary_upper;

  /**
   * Make sure we've been initialized before use
   */
  bool _initialized;

public:
  unsigned int _target_patch_size;

  /**
   * Constructor.
   */
  AugmentSparsityNonlocal(MeshBase & mesh,
    unsigned int target_patch_size);
                             // boundary_id_type crack_boundary_lower,
                             // boundary_id_type crack_boundary_upper);

  /**
   * @return a const reference to the lower-to-upper element ID map.
   */
  const ElementSideMap & get_lower_to_upper() const;

  /**
   * User-defined function to augment the sparsity pattern.
   */
  virtual void operator() (const MeshBase::const_element_iterator &
range_begin,
                           const MeshBase::const_element_iterator &
range_end,
                           processor_id_type p,
                           map_type & coupled_elements) override;

  /**
   * Rebuild the cached _lower_to_upper map whenever our Mesh has
   * changed.
   */
  virtual void mesh_reinit () override;

  /**
   * Update the cached _lower_to_upper map whenever our Mesh has been
   * redistributed.  We'll be lazy and just recalculate from scratch.
   */
  virtual void redistribute () override
  { this->mesh_reinit(); }


};

#endif












// local includes
#include "augment_sparsity_nonlocal.h"

// libMesh includes
#include "libmesh/elem.h"

using namespace libMesh;

AugmentSparsityNonlocal::AugmentSparsityNonlocal(MeshBase & mesh,
  unsigned int target_patch_size):
                                                       // boundary_id_type
crack_boundary_lower,
                                                       // boundary_id_type
crack_boundary_upper) :
  _mesh(mesh),
  // _crack_boundary_lower(crack_boundary_lower),
  // _crack_boundary_upper(crack_boundary_upper),
  _initialized(false)
{
  this->mesh_reinit();
  this->_target_patch_size = target_patch_size;

}


const ElementSideMap & AugmentSparsityNonlocal::get_lower_to_upper() const
{
  libmesh_assert(this->_initialized);
  return _lower_to_upper;
}

void AugmentSparsityNonlocal::mesh_reinit ()
{
  this->_initialized = true;

  // Loop over all elements (not just local elements) to make sure we find
  // "neighbor" elements on opposite sides of the crack.

  // Map from (elem, side) to centroid
  std::map<std::pair<const Elem *, unsigned char>, Point> lower_centroids;
  std::map<std::pair<const Elem *, unsigned char>, Point> upper_centroids;

  // for (const auto & elem : _mesh.active_element_ptr_range())
  //   for (auto side : elem->side_index_range())
  //     if (elem->neighbor_ptr(side) == nullptr)
  //       {
  //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
_crack_boundary_lower))
  //           {
  //             std::unique_ptr<const Elem> side_elem =
elem->build_side_ptr(side);

  //             lower_centroids[std::make_pair(elem, side)] =
side_elem->centroid();
  //           }

  //         if (_mesh.get_boundary_info().has_boundary_id(elem, side,
_crack_boundary_upper))
  //           {
  //             std::unique_ptr<const Elem> side_elem =
elem->build_side_ptr(side);

  //             upper_centroids[std::make_pair(elem, side)] =
side_elem->centroid();
  //           }
  //       }

  // If we're doing a reinit on a distributed mesh then we may not see
  // all the centroids, or even a matching number of centroids.
  // std::size_t n_lower_centroids = lower_centroids.size();
  // std::size_t n_upper_centroids = upper_centroids.size();
  // libmesh_assert(n_lower_centroids == n_upper_centroids);

  // Clear _lower_to_upper. This map will be used for matrix assembly later
on.
  _lower_to_upper.clear();

  // Clear _upper_to_lower. This map will be used for element ghosting
  // on distributed meshes, communication send_list construction in
  // parallel, and sparsity calculations
  _upper_to_lower.clear();

  // We do an N^2 search to find elements with matching centroids. This
could be optimized,
  // e.g. by first sorting the centroids based on their (x,y,z) location.
  {
    std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it
 = lower_centroids.begin();
    std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
it_end = lower_centroids.end();
    for ( ; it != it_end; ++it)
      {
        Point lower_centroid = it->second;

        // find closest centroid in upper_centroids
        Real min_distance = std::numeric_limits<Real>::max();

        std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it     = upper_centroids.begin();
        std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it_end = upper_centroids.end();

        for ( ; inner_it != inner_it_end; ++inner_it)
          {
            Point upper_centroid = inner_it->second;

            Real distance = (upper_centroid - lower_centroid).norm();
            if (distance < min_distance)
              {
                min_distance = distance;
                _lower_to_upper[it->first] = inner_it->first.first;
              }
          }

        // For pairs with local elements, we should have found a
        // matching pair by now.
        const Elem * elem     = it->first.first;
        const Elem * neighbor = _lower_to_upper[it->first];
        if (min_distance < TOLERANCE)
          {
            // fill up the inverse map
            _upper_to_lower[neighbor] = elem;
          }
        else
          {
            libmesh_assert_not_equal_to(elem->processor_id(),
_mesh.processor_id());
            // This must have a false positive; a remote element would
            // have been closer.
            _lower_to_upper.erase(it->first);
          }
      }

    // Let's make sure we didn't miss any upper elements either
// #ifndef NDEBUG
//     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it     = upper_centroids.begin();
//     std::map<std::pair<const Elem *, unsigned char>, Point>::iterator
inner_it_end = upper_centroids.end();

//     for ( ; inner_it != inner_it_end; ++inner_it)
//       {
//         const Elem * neighbor = inner_it->first.first;
//         if (neighbor->processor_id() != _mesh.processor_id())
//           continue;
//         ElementMap::const_iterator utl_it =
//           _upper_to_lower.find(neighbor);
//         libmesh_assert(utl_it != _upper_to_lower.end());
//       }
// #endif
  }
}

void AugmentSparsityNonlocal::operator()
  (const MeshBase::const_element_iterator & range_begin,
   const MeshBase::const_element_iterator & range_end,
   processor_id_type p,
   map_type & coupled_elements)
{
  libmesh_assert(this->_initialized);

  const CouplingMatrix * const null_mat = nullptr;

  for (const auto & elem : as_range(range_begin, range_end))
    {
      if (elem->processor_id() != p)
        coupled_elements.insert (std::make_pair(elem, null_mat));

      for (auto side : elem->side_index_range())
        if (elem->neighbor_ptr(side) == nullptr)
          {
            ElementSideMap::const_iterator ltu_it =
              _lower_to_upper.find(std::make_pair(elem, side));
            if (ltu_it != _lower_to_upper.end())
              {
                const Elem * neighbor = ltu_it->second;
                if (neighbor->processor_id() != p)
                  coupled_elements.insert (std::make_pair(neighbor,
null_mat));
              }
          }

      ElementMap::const_iterator utl_it =
        _upper_to_lower.find(elem);
      if (utl_it != _upper_to_lower.end())
        {
          const Elem * neighbor = utl_it->second;
          if (neighbor->processor_id() != p)
            coupled_elements.insert (std::make_pair(neighbor, null_mat));
        }

    }


/*
Nonlocal augment
*/
  // const CouplingMatrix * const null_mat = nullptr;
  libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
  // #ifndef this->_target_patch_size
  // const unsigned int _target_patch_size=4;
  // #endif
/* build element patches */
  for (const auto & elem : as_range(range_begin, range_end))
    {
      libMesh::Patch patch(_mesh.processor_id());
      patch.build_around_element(elem, _target_patch_size, patchtype);
      std::vector<const Elem *> patchvec;
      patchvec.reserve(patch.size());
      Patch::const_iterator        patch_it  = patch.begin();
      const Patch::const_iterator  patch_end = patch.end();
      for (; patch_it != patch_end; ++patch_it)
      {
        const Elem * elem2 = *patch_it;
        coupled_elements.insert (std::make_pair(elem2, null_mat));
          // patchvec.push_back(elem2);
      }
    }
}



// MatSetOption(this->get_matrix("NonlocalStiffness"),
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
//     DoFRenumbering::Cuthill_McKee (dof_handler);
//     hanging_node_constraints.clear();
//
 DoFTools::make_hanging_node_constraints(dof_handler,hanging_node_constraints);

//     hanging_node_constraints.close();

//     DynamicSparsityPattern dsp(dof_handler.n_dofs(),
dof_handler.n_dofs());
//     DoFTools::make_sparsity_pattern(dof_handler, dsp);
//     hanging_node_constraints.condense(dsp);
//     sparsity_pattern.copy_from(dsp);

//     stiffness_matrix.reinit(sparsity_pattern);
//     mass_matrix.reinit(sparsity_pattern);