DG_EXO_test.cpp

// C++ include files that we need
#include <iostream>

// Basic include files needed for the mesh functionality.
#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/exodusII_io.h"

// Define the Finite Element object.
#include "libmesh/fe.h"

// Define Gauss quadrature rules.
#include "libmesh/quadrature_gauss.h"

// Define useful datatypes for finite element
// matrix and vector components.
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"

// Define the DofMap, which handles degree of freedom
// indexing.
#include "libmesh/dof_map.h"

// Bring in everything from the libMesh namespace
using namespace libMesh;


int main (int argc, char ** argv)
{
  // Initialize libraries, like in example 2.
  LibMeshInit init (argc, argv);
  // Create a mesh, with dimension to be overridden later, distributed
  // across the default MPI communicator.
  Mesh mesh(init.comm());
  // Use the MeshTools::Generation mesh generator to create a uniform
  // 2D grid on the square [-1,1]^2.
  MeshTools::Generation::build_square (mesh,
                                        1,   1,
                                       -1., 1.,
                                       -1., 1.,
                                       TRI3);

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);
  // Initialize the data structures for the equation system.
  equation_systems.add_system<LinearImplicitSystem> ("DG");

  // Adds the variable "u" to "DG".  "u"
  // will be approximated using first-order approximation.
  equation_systems.get_system("DG").add_variable("u", FIRST, L2_LAGRANGE);
  equation_systems.init();
  // Prints information about the system to the screen.
  equation_systems.print_info();

  //INITIALIZE u VECTOR:
  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  std::vector<libMesh::dof_id_type> dof_indices;
  libMesh::DenseVector<libMesh::Number> Fe;
  const libMesh::DofMap & dof_map = equation_systems.get_system("DG").get_dof_map();
  libMesh::FEType fe_disp = dof_map.variable_type(0);
  const unsigned int dim = mesh.mesh_dimension();

  for (; el != end_el; ++el)
  {
      const libMesh::Elem * elem = *el;
      auto elID = elem->id();

      dof_map.dof_indices(elem, dof_indices);
      Fe.resize( dof_indices.size() );
      Fe(0) = dof_indices[0];
      Fe(1) = dof_indices[1];
      Fe(2) = dof_indices[2];

      equation_systems.get_system("DG").solution->insert(Fe, dof_indices);
  }
  equation_systems.get_system("DG").solution->close();
  equation_systems.get_system("DG").solution->print(std::cout);
  
  ExodusII_IO exp(mesh);
  exp.write_discontinuous_exodusII("dg.e", equation_systems);
  exp.append(true);
  // COMMENTING THE LINE BELOW THE DATA ARE EXPORTED CORRECTLY
  exp.write_timestep("dg.e", equation_systems, 1, 3.0);
  return 0;
}