Menu
▾
▴
Re: [Libmesh-users] Augmenting Sparsity to add Nonlocal DOFs
|
From: Charlie T. <cha...@gm...> - 2018-09-18 08:06:03
|
After looking at example 9 again, I have to ask:
Better yet: is there a way to augment the sparsity pattern without this
function? I've read that I should
"
To do this, you need to cast your matrix to a PetscMatrix, then do
petsc_matrix->mat() to get the PETSc Mat that you can pass to
MatSetOption.
PetscMatrix<Number> * NonlocalStiffness_petsc =
&(this->get_matrix("NonlocalStiffness"));
MatSetOption( NonlocalStiffness_petsc.mat() ,
MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
"
Suppose I added a sparse matrix "TestMat" to the system. How do I do this,
explicitly?
Thanks!
On Tue, Sep 18, 2018 at 3:40 AM Charlie Talbot <cha...@gm...> wrote:
> 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);
>
>
>
>
>
>
>
>
>
>
|
