Minimal libMesh example with parallel leak

libmesh_minimal_example.C

#include <math.h>

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/equation_systems.h"

#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"

#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"
#include "libmesh/enum_solver_package.h"

#include "libmesh/dof_map.h"

using namespace libMesh;

void assemble_poisson(EquationSystems & es,
                      const std::string & system_name);

Real exact_solution (const Real x,
                     const Real y,
                     const Real z = 0.)
{
  static const Real pi = acos(-1.);

  return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
}

int main (int argc, char ** argv)
{
  LibMeshInit init (argc, argv);

  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");

  Mesh mesh(init.comm());

  MeshTools::Generation::build_square (mesh,
                                       2, 1,
                                       -1., 1.,
                                       -1., 1.,
                                       QUAD9);

  EquationSystems equation_systems (mesh);
  equation_systems.add_system<LinearImplicitSystem> ("Poisson");
  equation_systems.get_system("Poisson").add_variable("u", SECOND);
  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
  equation_systems.init();

  equation_systems.get_system("Poisson").solve();

  return 0;
}



void assemble_poisson(EquationSystems & es,
                      const std::string &)
{
  const MeshBase & mesh = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();

  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> ("Poisson");

  const DofMap & dof_map = system.get_dof_map();

  FEType fe_type = dof_map.variable_type(0);
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
  QGauss qrule (dim, FIFTH);
  fe->attach_quadrature_rule (&qrule);

  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
  QGauss qface(dim-1, FIFTH);
  fe_face->attach_quadrature_rule (&qface);

  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<Point> & q_point = fe->get_xyz();
  const std::vector<std::vector<Real>> & phi = fe->get_phi();
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();

  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  std::vector<dof_id_type> dof_indices;

  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      dof_map.dof_indices (elem, dof_indices);

      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());

      fe->reinit (elem);

      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);

      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          for (unsigned int i=0; i != n_dofs; i++)
            for (unsigned int j=0; j != n_dofs; j++)
              {
                Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
              }
          {
            const Real x = q_point[qp](0);
            const Real y = q_point[qp](1);
            const Real eps = 1.e-3;

            const Real fxy = -(exact_solution(x, y-eps) +
                               exact_solution(x, y+eps) +
                               exact_solution(x-eps, y) +
                               exact_solution(x+eps, y) -
                               4.*exact_solution(x, y))/eps/eps;

            for (unsigned int i=0; i != n_dofs; i++)
              Fe(i) += JxW[qp]*fxy*phi[i][qp];
          }
        }

      {

        for (auto side : elem->side_index_range())
          if (elem->neighbor_ptr(side) == nullptr)
            {
              const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
              const std::vector<Real> & JxW_face = fe_face->get_JxW();
              const std::vector<Point> & qface_point = fe_face->get_xyz();

              fe_face->reinit(elem, side);


              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  const Real xf = qface_point[qp](0);
                  const Real yf = qface_point[qp](1);

                  const Real penalty = 1.e10;

                  const Real value = exact_solution(xf, yf);

                  for (unsigned int i=0; i != n_dofs; i++)
                    for (unsigned int j=0; j != n_dofs; j++)
                      Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];

                  for (unsigned int i=0; i != n_dofs; i++)
                    Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
                }
            }
      }
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      system.matrix->add_matrix (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);
    }
}

Makefile

######################################################################
#
# Template libMesh application Makefile
LIBMESH_DIR ?= /ssd/lindad/moose/libmesh/installed-mpich3.3b1-all-dbg


# include the library options determined by configure
include $(LIBMESH_DIR)/Make.common

target     := ./example-$(METHOD)


###############################################################################
# File management.  This is where the source, header, and object files are
# defined

#
# source files
srcfiles 	:= $(wildcard *.C)

#
# object files
objects		:= $(patsubst %.C, %.$(obj-suffix), $(srcfiles))
###############################################################################



.PHONY: dust clean distclean

###############################################################################
# Target:
#

all:: $(notdir $(target))

# Production rules:  how to make the target - depends on library configuration
$(notdir $(target)): $(objects)
	@echo "Linking "$@"..."
	@$(libmesh_LIBTOOL) --tag=CXX $(LIBTOOLFLAGS) --mode=link \
	  $(libmesh_CXX) $(libmesh_CXXFLAGS) $(objects) -o $@ $(libmesh_LIBS) $(libmesh_LDFLAGS) $(EXTERNAL_FLAGS)


# Useful rules.
dust:
	@echo "Deleting old output and runtime files"
	@rm -f out*.m job_output.txt output.txt* *.gmv.* *.plt.* out*.xdr* out*.xda* PI* complete

clean: dust
	@rm -f $(objects) *.$(obj-suffix) *.lo

clobber: clean 
	@rm -f $(target)

distclean: clean
	@rm -rf *.o .libs

echo:
	@echo srcfiles = $(srcfiles)
	@echo objects = $(objects)
	@echo target = $(target)

run: complete

complete: $(wildcard *.in)
#	@$(MAKE) dust
	@$(MAKE) -C $(dir $(target)) $(notdir $(target))
	@echo "***************************************************************"
	@echo "* Running App " $(notdir $(target))
	@echo "***************************************************************"
	@echo " "
	${LIBMESH_RUN} $(target) ${LIBMESH_OPTIONS} 2>&1 | tee output.txt
	@bzip2 -f output.txt
	@echo " "
	@echo "***************************************************************"
	@echo "* Done Running App " $(notdir $(target))
	@echo "***************************************************************"

###############################################################################