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);
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
