assembly_bench.C

// TYPE
// 0 - Normal
// 1 - DenseMatrix
// 2 - Eigen
// 3 - ColumnMajorMatrix
#define TYPE 0


#include "Var_evalApp.h"
#include "MooseInit.h"
#include "Moose.h"
#include "MooseApp.h"
#include "AppFactory.h"
#include "ColumnMajorMatrix.h"

// Create a performance log
PerfLog Moose::perf_log("Var_eval");

// C++ include files that we need
#include <iostream>
#include <algorithm>
#include <sstream>
#include <math.h>

// Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/perf_log.h"
#include "libmesh/boundary_info.h"
#include "libmesh/utility.h"

// For systems of equations the \p DenseSubMatrix
// and \p DenseSubVector provide convenient ways for
// assembling the element matrix and vector on a
// component-by-component basis.
#include "libmesh/dense_submatrix.h"
#include "libmesh/dense_subvector.h"

// The definition of a geometric element
#include "libmesh/elem.h"

#include "Eigen/Core"

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

// Function prototype.  This function will assemble the system
// matrix and right-hand-side.
void assemble_stokes (EquationSystems& es,
                      const std::string& system_name);

// The main program.
int main (int argc, char** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

  // Skip this 2D example if libMesh was compiled as 1D-only.
  libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");

  // Trilinos solver NaNs by default on the zero pressure block.
  // We'll skip this example for now.
  if (libMesh::default_solver_package() == TRILINOS_SOLVERS)
    {
      std::cout << "We skip example 13 when using the Trilinos solvers.\n"
                << std::endl;
      return 0;
    }

  // 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.  We instruct the mesh generator
  // to build a mesh of 8x8 \p Quad9 elements in 2D.  Building these
  // higher-order elements allows us to use higher-order
  // approximation, as in example 3.
  
  MeshTools::Generation::build_cube (mesh,
                                     25, 25, 25,
                                     0., 1.,
                                     0., 1.,
                                     0., 1.,
                                     HEX27);
  
//  mesh.read("cube.e");

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

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Creates a transient system named "Navier-Stokes"
  TransientExplicitSystem & system =
    equation_systems.add_system<TransientExplicitSystem> ("Navier-Stokes");

  // Add the variables "u" & "v" to "Navier-Stokes".  They
  // will be approximated using second-order approximation.
  for(unsigned int var=0; var<2000; var++)
  {
    std::ostringstream name;
    name << "u" <<var;
  
    system.add_variable (name.str(), FIRST, LAGRANGE);
  }

  // Give the system a pointer to the matrix assembly
  // function.
  system.attach_assemble_function (assemble_stokes);

  // Initialize the data structures for the equation system.
  equation_systems.init ();

  (*system.solution) = 4.5;

  equation_systems.update();
  
  // Prints information about the system to the screen.
  equation_systems.print_info();

  // Create a performance-logging object for this example
  PerfLog perf_log("Systems Example 2");

  // Get a reference to the Stokes system to use later.
  TransientExplicitSystem&  navier_stokes_system =
        equation_systems.get_system<TransientExplicitSystem>("Navier-Stokes");

  for(unsigned int i=0; i<5; i++)
  {
    std::cout<<"Sleeping: "<<i<<std::endl;
    sleep(1);
  }

  std::cout<<"Assembling"<<std::endl;

  Moose::perf_log.push("Assembling");
  
  for(unsigned int i=0; i<10; i++)
    equation_systems.get_system("Navier-Stokes").assemble();

  Moose::perf_log.pop("Assembling");

  std::cout<<"Finished Assembling"<<std::endl;

  // All done.
  return 0;
}

// The matrix assembly function to be called at each time step to
// prepare for the linear solve.
void assemble_stokes (EquationSystems& es,
                      const std::string& system_name)
{
  // Get a constant reference to the mesh object.
  const MeshBase& mesh = es.get_mesh();

  // The dimension that we are running
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the Stokes system object.
  TransientExplicitSystem & navier_stokes_system =
    es.get_system<TransientExplicitSystem> ("Navier-Stokes");

  FEType fe_vel_type = navier_stokes_system.variable_type(0);

  unsigned int n_vars = navier_stokes_system.n_vars();

  AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));

  QGauss qrule (dim, fe_vel_type.default_quadrature_order());

  fe_vel->attach_quadrature_rule (&qrule);

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

  const DofMap & dof_map = navier_stokes_system.get_dof_map();

//  std::vector<std::vector<dof_id_type> > dof_indices;
  std::vector<dof_id_type> dof_indices;
  
  std::vector<unsigned int> dof_indices_vars(n_vars);

  for(unsigned int i=0; i<n_vars; i++)
    dof_indices_vars[i] = i;

  NumericVector<Number> & solution = *navier_stokes_system.current_local_solution.get();

  MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();


#if TYPE==0
  std::vector<std::vector<Real> > u_vals(n_vars);
  std::vector<std::vector<RealGradient> > grad_u_vals(n_vars);
#endif

#if TYPE==1
  DenseMatrix<Real>  phi_qp_dof;
  DenseMatrix<Real>  dphi_qp_dof;
  DenseMatrix<Real>  dof_values;
  
  DenseMatrix<Real>  u;
  DenseMatrix<Real>  grad_u;
#endif
  
#if TYPE==2
  Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> phi_qp_dof;
  Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> dphi_qp_dof;
  Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> dof_values;
  Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> u;

  Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> grad_u;
#endif

  for ( ; el != end_el; ++el)
  {
    // Store a pointer to the element we are currently
    // working on.  This allows for nicer syntax later.
    const Elem* elem = *el;
      
    // Compute the element-specific data for the current
    // element.  This involves computing the location of the
    // quadrature points (q_point) and the shape functions
    // (phi, dphi) for the current element.
    fe_vel->reinit  (elem);
    
    const unsigned int n_dofs = phi.size();
    const unsigned int n_points = qrule.n_points();


#if TYPE==1 || TYPE==2
    phi_qp_dof.resize(n_points, n_dofs);
    dof_values.resize(n_dofs, n_vars);
    u.resize(n_points, n_vars);

    dphi_qp_dof.resize(LIBMESH_DIM*n_points, n_dofs);
    grad_u.resize(LIBMESH_DIM*n_points, n_vars);
#endif

#if TYPE==3
    ColumnMajorMatrix phi_qp_dof(n_points, n_dofs);
    ColumnMajorMatrix dof_values(n_dofs, n_vars);
    ColumnMajorMatrix u(n_points, n_vars);

    ColumnMajorMatrix dphi_qp_dof(LIBMESH_DIM*n_points, n_dofs);
    ColumnMajorMatrix grad_u(LIBMESH_DIM*n_points, n_vars);
#endif

    
#if TYPE==1 || TYPE==2
    // Copy out the phi values
    for (unsigned int qp=0; qp<n_points; qp++)
      for (unsigned int l=0; l<n_dofs; l++)
      {
        phi_qp_dof(qp,l) = phi[l][qp];

        for(unsigned int d=0; d<LIBMESH_DIM; d++)
          dphi_qp_dof((qp*LIBMESH_DIM)+d,l) = dphi[l][qp](d);
      } 

//    dof_map.dof_indices_for_vars_of_same_type(elem, dof_indices, dof_indices_vars);

    for(unsigned int var=0; var<n_vars; var++)
    {
      dof_map.dof_indices (elem, dof_indices, var);
      for (unsigned int l=0; l<n_dofs; l++)
      {
//        Real val = solution(dof_indices[var][l]);
        Real val = solution(dof_indices[l]);
        
        dof_values(l,var) = val;        
      }
    }
#endif    

#if TYPE==0
    for(unsigned int var=0; var<n_vars; var++)
    {
      u_vals[var].resize(n_points);
      grad_u_vals[var].resize(n_points);
      for(unsigned int qp=0; qp<n_points; qp++)
      {
        u_vals[var][qp] = 0;
        grad_u_vals[var][qp] = 0;
      }
    }

//    dof_map.dof_indices_for_vars_of_same_type(elem, dof_indices, dof_indices_vars);

    for(unsigned int var=0; var<n_vars; var++)
    {      
      dof_map.dof_indices (elem, dof_indices, var);

      for (unsigned int l=0; l<n_dofs; l++)
      {
//        Real val = solution(dof_indices[var][l]);
        Real val = solution(dof_indices[l]);
        for(unsigned int qp=0; qp<n_points; qp++)
        {
          u_vals[var][qp] += phi[l][qp]*val;
          grad_u_vals[var][qp].add_scaled (dphi[l][qp],val);
        }
      }
    }
#endif

#if TYPE==1    
    phi_qp_dof.right_multiply(dof_values);
    dphi_qp_dof.right_multiply(dof_values);
#endif

#if TYPE==2  
    u.noalias() = phi_qp_dof * dof_values;
    grad_u.noalias() = dphi_qp_dof * dof_values;
#endif

#if TYPE==3
    u = phi_qp_dof * dof_values;
    grad_u = dphi_qp_dof * dof_values;
#endif
  }
    
  return;
}