In particular, here is a slightly modified assemble function (from a
miscellaneous example for IPDG). For some reason, this does not work
void assemble_nonlocal(EquationSystems & es,
const std::string & libmesh_dbg_var(system_name))
{
libMesh::out << " assembling nonlocal stiffness... ";libMesh::out.flush();
// // _NonlocalStiffness->zero();
// EquationSystems & es = this->get_equation_systems();
// const MeshBase & mesh = es.get_mesh();
// const unsigned int dim = mesh.mesh_dimension();
const MeshBase & mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
// Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem & nonlocal_system =
es.get_system<LinearImplicitSystem> ("NonlocalSystem");
// Get some parameters that we need during assembly
const Real penalty = es.parameters.get<Real> ("penalty");
const unsigned int target_patch_size = es.parameters.get<unsigned int>
("target_patch_size");
const Real kernelparam = es.parameters.get<Real> ("kernelparam");
const Real horizon = es.parameters.get<Real> ("horizon");
const Order _quad_order = es.parameters.get<Order> ("_quad_order");
std::string refinement_type = es.parameters.get<std::string>
("refinement");
// A reference to the DofMap object for this system. The DofMap
// object handles the index translation from node and element numbers
// to degree of freedom numbers. We will talk more about the DofMap
const DofMap & dof_map = nonlocal_system.get_dof_map();
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = nonlocal_system.variable_type(0);
// Build a Finite Element object of the specified type. Since the
// FEBase::build() member dynamically creates memory we will
// store the object as a std::unique_ptr<FEBase>. This can be thought
// of as a pointer that will clean up after itself.
std::unique_ptr<FEBase> fe_elem (FEBase::build(dim, fe_type));
std::unique_ptr<FEBase> fe_neighbor (FEBase::build(dim, fe_type));
std::unique_ptr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
std::unique_ptr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
// Quadrature rules for numerical integration.
#ifdef _quad_type
// QGauss qrule (dim, QORDER);
// std::unique_ptr<QBase> qrule(QBase::build(_quad_type, dim,
_quad_order));
std::unique_ptr<QBase> qrule(QBase::build(_quad_type, dim, _quad_order));
std::unique_ptr<QBase> qrule_neigbor(QBase::build(_quad_type, dim,
_quad_order));
#else
std::unique_ptr<QBase> qrule(QBase::build(0, dim,
fe_type.default_quadrature_order()));
std::unique_ptr<QBase> qrule_neighbor(QBase::build(0, dim,
fe_type.default_quadrature_order()));
#endif
// fe->attach_quadrature_rule (&qrule);
fe_elem->attach_quadrature_rule (qrule.get());
fe_neighbor->attach_quadrature_rule (qrule_neighbor.get());
#ifdef _quad_type
std::unique_ptr<QBase> qface(QBase::build(_quad_type, dim, _quad_order));
#else
std::unique_ptr<QBase> qface(QBase::build(0, dim,
fe_type.default_quadrature_order()));
#endif
fe_elem_face->attach_quadrature_rule(qrule.get());
fe_neighbor_face->attach_quadrature_rule(qface.get());
const std::vector<std::vector<Real>> & phi_elem = fe_elem->get_phi();
const std::vector<std::vector<RealGradient>> & dphi_elem =
fe_elem->get_dphi();
const std::vector<Real> & JxW_elem = fe_elem->get_JxW();
const std::vector<Point> & qpoint_elem = fe_elem->get_xyz();
const std::vector<std::vector<Real>> & phi_neighbor =
fe_neighbor->get_phi();
const std::vector<std::vector<RealGradient>> & dphi2 =
fe_neighbor->get_dphi();
const std::vector<Real> & JxW_neighbor = fe_neighbor->get_JxW();
const std::vector<Point> & qpoint_neighbor = fe_neighbor->get_xyz();
const std::vector<Point> & qpoint_elem_global = fe_elem->get_xyz();
const std::vector<Point> & qpoint_neighbor_global =
fe_neighbor->get_xyz();
// Data for surface integrals on the element boundary
const std::vector<std::vector<Real>> & phi_face =
fe_elem_face->get_phi();
const std::vector<std::vector<RealGradient>> & dphi_face =
fe_elem_face->get_dphi();
const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
const std::vector<Point> & qface_normals = fe_elem_face->get_normals();
const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
// // Data for surface integrals on the neighbor boundary
const std::vector<std::vector<Real>> & phi_neighbor_face =
fe_neighbor_face->get_phi();
const std::vector<std::vector<RealGradient>> & dphi_neighbor_face =
fe_neighbor_face->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
// Data structures to contain the element and neighbor boundary matrix
// contribution. This matrices will do the coupling between the dofs of
// the element and those of his neighbors.
// Ken: matrix coupling elem and neighbor dofs
DenseMatrix<Number> Kne;
DenseMatrix<Number> Ken;
DenseMatrix<Number> Kee;
DenseMatrix<Number> Knn;
DenseMatrix<Number> zero_matrix, zero_matrix1, zero_matrix2;
std::vector<dof_id_type> dof_indices;
libMesh::Patch::PMF patchtype = &Patch::add_face_neighbors;
libMesh::Patch patch(mesh.processor_id());
// for (const auto& elem : mesh.active_local_element_ptr_range())
for (const auto& elem : mesh.active_element_ptr_range())
{
patch.build_around_element(elem, target_patch_size, patchtype);
fe_elem->reinit(elem);
dof_map.dof_indices (elem, dof_indices);
const unsigned int n_dofs = dof_indices.size();
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
Kee.resize (n_dofs, n_dofs);
/* reposition */
// Kee.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs,
n_u_dofs);
libMesh::out << " iterate over elements in patch ";
libMesh::out.flush();
for (auto neighbor : patch)
{
fe_neighbor->reinit(neighbor);
std::vector<dof_id_type> neighbor_dof_indices;
dof_map.dof_indices (neighbor, neighbor_dof_indices);
const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
Ken.resize (n_dofs, n_neighbor_dofs);
Kne.resize (n_neighbor_dofs, n_dofs);
Ken.resize (n_neighbor_dofs, n_neighbor_dofs);
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
{
// this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
dof_indices, dof_indices);
// this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
neighbor_dof_indices, neighbor_dof_indices);
// this->get_matrix("NonlocalStiffness").add_matrix(zero_matrix,
neighbor_dof_indices, dof_indices);
libMesh::out << " fill with zeros 1 "; libMesh::out.flush();
zero_matrix1.resize (n_dofs, n_neighbor_dofs);
nonlocal_system.matrix->add_matrix(zero_matrix1, dof_indices,
neighbor_dof_indices);
//
nonlocal_system.matrix->get_matrix("NonlocalStiffness").close();
libMesh::out << " fill with zeros 2 "; libMesh::out.flush();
zero_matrix2.resize (n_neighbor_dofs,n_dofs);
nonlocal_system.matrix->add_matrix(zero_matrix2,
neighbor_dof_indices, dof_indices);
//
nonlocal_system.matrix->get_matrix("NonlocalStiffness").close();
for (unsigned int qp2=0; qp2 < qrule->n_points(); qp2++)
{
Real kernelval = kernel(qpoint_elem_global[qp],
qpoint_neighbor_global[qp2], kernelparam);
libMesh::out << " kernelval " << kernelval;
libMesh::out.flush();
for (unsigned int i=0; i<n_dofs; i++)
for (unsigned int j=0; j<n_neighbor_dofs; j++)
Ken(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_elem[j][qp]);
for (unsigned int i=0; i<n_neighbor_dofs; i++)
for (unsigned int j=0; j<n_dofs; j++)
Kne(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_elem[j][qp]);
for (unsigned int i=0; i<n_dofs; i++)
for (unsigned int j=0; j<n_dofs; j++)
Kee(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_elem[i][qp2] * phi_elem[j][qp]);
for (unsigned int i=0; i<n_neighbor_dofs; i++)
for (unsigned int j=0; j<n_neighbor_dofs; j++)
Knn(i,j) += JxW_neighbor[qp2] * JxW_elem[qp] *
kernelval * (phi_neighbor[i][qp2] * phi_neighbor[j][qp]);
}
// We consider Dirichlet bc imposed via the interior penalty
method
for (auto side : elem->side_index_range())
{
if (elem->neighbor_ptr(side) == nullptr)
{
// Pointer to the element face
fe_elem_face->reinit(elem, side);
std::unique_ptr<const Elem> elem_side
(elem->build_side_ptr(side));
// h element dimension to compute the interior penalty
penalty parameter
const unsigned int elem_b_order = static_cast<unsigned
int> (fe_elem_face->get_order());
const double h_elem =
elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
for (unsigned int qp=0; qp<qface->n_points(); qp++)
{
Number bc_value = exact_solution(qface_points[qp],
es.parameters, "null", "void");
for (unsigned int i=0; i<n_dofs; i++)
{
// Matrix contribution
for (unsigned int j=0; j<n_dofs; j++)
{
// stability
Ke(i,j) += JxW_face[qp] * penalty/h_elem *
phi_face[i][qp] * phi_face[j][qp];
// consistency
Ke(i,j) -=
JxW_face[qp] *
(phi_face[i][qp] *
(dphi_face[j][qp]*qface_normals[qp]) +
phi_face[j][qp] *
(dphi_face[i][qp]*qface_normals[qp]));
}
// RHS contributions
// stability
Fe(i) += JxW_face[qp] * bc_value *
penalty/h_elem * phi_face[i][qp];
// consistency
Fe(i) -= JxW_face[qp] * dphi_face[i][qp] *
(bc_value*qface_normals[qp]);
}
}
}
}
libMesh::out << " constrain nonlocal stiffness...
";libMesh::out.flush();
dof_map.constrain_element_matrix(Kee, dof_indices,
dof_indices);
dof_map.constrain_element_matrix(Ken, dof_indices,
neighbor_dof_indices);
dof_map.constrain_element_matrix(Kne, neighbor_dof_indices,
dof_indices);
dof_map.constrain_element_matrix(Knn, neighbor_dof_indices,
neighbor_dof_indices);
libMesh::out << " add nonlocal stiffness...
";libMesh::out.flush();
nonlocal_system.matrix->add_matrix(Kee, dof_indices,
dof_indices);
nonlocal_system.matrix->add_matrix(Ken, dof_indices,
neighbor_dof_indices);
nonlocal_system.matrix->add_matrix(Kne, neighbor_dof_indices,
dof_indices);
nonlocal_system.matrix->add_matrix(Knn, neighbor_dof_indices,
neighbor_dof_indices);
}//qp
}//patch
dof_map.constrain_element_vector(Fe, dof_indices);
nonlocal_system.rhs->add_vector(Fe, dof_indices);
// this->_NonlocalStiffness->add_matrix(Kee, dof_indices);
}//elem
libMesh::out << " assembled nonlocal stiffness... ";libMesh::out.flush();
}//fcn
On Tue, Sep 18, 2018 at 4:05 AM Charlie Talbot <chazti...@gmail.com> wrote:
> 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 <chazti...@gmail.com>
> 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);
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
