Miscellaneous Example 2 - Complex Numbers and the "FrequencySystem"

This is the seventh example program. It builds on the previous example programs, introduces complex numbers and the FrequencySystem class to solve a simple Helmholtz equation grad(p)*grad(p)+(omega/c)^2*p=0, for multiple frequencies rather efficiently.

The FrequencySystem class offers two solution styles, namely to solve large systems, or to solve moderately-sized systems fast, for multiple frequencies. The latter approach is implemented here.

This example uses an L--shaped mesh and nodal boundary data given in the files lshape.un and lshape_data.unv

For this example the library has to be compiled with complex numbers enabled.

C++ include files that we need

        #include <iostream>
        #include <algorithm>
        #include <stdio.h>

Basic include files needed for overall functionality.

        #include "libmesh.h"
        #include "libmesh_logging.h"
        #include "mesh.h"
        #include "mesh_generation.h"
        #include "exodusII_io.h"
        #include "equation_systems.h"
        #include "elem.h"

Include FrequencySystem. Compared to GeneralSystem, this class offers added functionality for the solution of frequency-dependent systems.

        #include "frequency_system.h"

Define the Finite Element object.

        #include "fe.h"

Define Gauss quadrature rules.

        #include "quadrature_gauss.h"

Define useful datatypes for finite element matrix and vector components.

        #include "dense_matrix.h"
        #include "dense_vector.h"

Define matrix and vector data types for the global equation system. These are base classes, from which specific implementations, like the PETSc or LASPACK implementations, are derived.

        #include "sparse_matrix.h"
        #include "numeric_vector.h"

Define the DofMap, which handles degree of freedom indexing.

        #include "dof_map.h"

Defines the MeshData class, which allows you to store data about the mesh when reading in files, etc.

        #include "mesh_data.h"

Bring in everything from the libMesh namespace

        using namespace libMesh;

This problem is only defined on complex-valued fields, for which libMesh must be configured with Number == complex.

        #ifdef LIBMESH_USE_COMPLEX_NUMBERS

Function prototype. This is the function that will assemble the mass, damping and stiffness matrices. It will not form an overall system matrix ready for solution.

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

Function prototype. This is the function that will combine the previously-assembled mass, damping and stiffness matrices to the overall matrix, which then renders ready for solution.

        void add_M_C_K_helmholtz(EquationSystems& es,
                                 const std::string& system_name);
        #endif

Begin the main program. Note that this example only works correctly if complex numbers have been enabled in the library. In order to link against the complex PETSc libraries, you must have built PETSc with the same C++ compiler that you used to build libMesh. This is so that the name mangling will be the same for the routines in both libraries.

        int main (int argc, char** argv)
        {

Initialize libraries, like in example 2.

          LibMeshInit init (argc, argv);

This example is designed for complex numbers.

        #ifndef LIBMESH_USE_COMPLEX_NUMBERS
          libmesh_example_assert(false, "--enable-complex");
        #else

Check for proper usage.

          if (argc < 3)
            {
              if (libMesh::processor_id() == 0)
                std::cerr << "Usage: " << argv[0] << " -f [frequency]"
                          << std::endl;
              
              libmesh_error();
            }
        
          if (libMesh::n_processors() > 1)
            {
              if (libMesh::processor_id() == 0)
                {
                  std::cerr << "ERROR: Skipping example 7. " << std::endl;
                  std::cerr << "MeshData objects currently only work in serial." << std::endl;
                }
              return 0;
            }

Tell the user what we are doing.

          else 
            {
              std::cout << "Running " << argv[0];
              
              for (int i=1; i<argc; i++)
                std::cout << " " << argv[i];
              
              std::cout << std::endl << std::endl;
            }

Get the frequency from argv[2] as a float, currently, solve for 1/3rd, 2/3rd and 1/1th of the given frequency

          const Real frequency_in = atof(argv[2]);
          const unsigned int n_frequencies = 3;

Skip this 2D example if libMesh was compiled as 1D-only.

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

Create a mesh.

          Mesh mesh;

Create a corresponding MeshData and activate it. For more information on this object cf. example 12.

          MeshData mesh_data(mesh);
          mesh_data.activate();

Read the mesh file. Here the file lshape.unv contains an L--shaped domain in .unv format.

          mesh.read("lshape.unv", &mesh_data);

Print information about the mesh to the screen.

          mesh.print_info();

The load on the boundary of the domain is stored in the .unv formated mesh data file lshape_data.unv. At this, the data is given as complex valued normal velocities.

          mesh_data.read("lshape_data.unv");

Print information about the mesh to the screen.

          mesh_data.print_info();

Create an equation systems object, which now handles a frequency system, as opposed to previous examples. Also pass a MeshData pointer so the data can be accessed in the matrix and rhs assembly.

          EquationSystems equation_systems (mesh, &mesh_data);

Create a FrequencySystem named "Helmholtz" & store a reference to it.

          FrequencySystem & f_system =      
            equation_systems.add_system<FrequencySystem> ("Helmholtz");

Add the variable "p" to "Helmholtz". "p" will be approximated using second-order approximation.

          f_system.add_variable("p", SECOND);

Tell the frequency system about the two user-provided functions. In other circumstances, at least the solve function has to be attached.

          f_system.attach_assemble_function (assemble_helmholtz);
          f_system.attach_solve_function    (add_M_C_K_helmholtz);

To enable the fast solution scheme, additional global matrices and one global vector, all appropriately sized, have to be added. The system object takes care of the appropriate size, but the user should better fill explicitly the sparsity structure of the overall matrix, so that the fast matrix addition method can be used, as will be shown later.

          f_system.add_matrix ("stiffness");
          f_system.add_matrix ("damping");
          f_system.add_matrix ("mass");
          f_system.add_vector ("rhs");

Communicates the frequencies to the system. Note that the frequency system stores the frequencies as parameters in the equation systems object, so that our assemble and solve functions may directly access them. Will solve for 1/3rd, 2/3rd and 1/1th of the given frequency

          f_system.set_frequencies_by_steps (frequency_in/n_frequencies,
                                             frequency_in,
                                             n_frequencies);

Use the parameters of the equation systems object to tell the frequency system about the wave velocity and fluid density. The frequency system provides default values, but these may be overridden, as shown here.

          equation_systems.parameters.set<Real> ("wave speed") = 1.;
          equation_systems.parameters.set<Real> ("rho")        = 1.;

Initialize the data structures for the equation system. Always prior to this, the frequencies have to be communicated to the system.

          equation_systems.init ();

Prints information about the system to the screen.

          equation_systems.print_info ();
        
          for (unsigned int n=0; n < n_frequencies; n++)
            {

Solve the system "Helmholtz" for the n-th frequency. Since we attached an assemble() function to the system, the mass, damping and stiffness contributions will only be assembled once. Then, the system is solved for the given frequencies. Note that solve() may also solve the system only for specific frequencies.

              f_system.solve (n,n);

After solving the system, write the solution to an ExodusII-formatted plot file, for every frequency.

        #ifdef LIBMESH_HAVE_EXODUS_API
              char buf[14];
              sprintf (buf, "out%04d.exd", n);
        
              ExodusII_IO(mesh).write_equation_systems (buf,
                                                  equation_systems);
        #endif
            }

Alternatively, the whole EquationSystems object can be written to disk. By default, the additional vectors are also saved.

          equation_systems.write ("eqn_sys.dat", libMeshEnums::WRITE);

All done.

          return 0;
        
        #endif 
        }
        
        
        #ifdef LIBMESH_USE_COMPLEX_NUMBERS

Here we define the matrix assembly routine for the Helmholtz system. This function will be called to form the stiffness matrix and right-hand side.

        void assemble_helmholtz(EquationSystems& es,
                                const std::string& system_name)
        {

It is a good idea to make sure we are assembling the proper system.

          libmesh_assert (system_name == "Helmholtz");

Get a constant reference to the mesh object.

          const MeshBase& mesh = es.get_mesh();

The dimension that we are in

          const unsigned int dim = mesh.mesh_dimension();

Get a reference to our system, as before

          FrequencySystem & f_system =
            es.get_system<FrequencySystem> (system_name);

A const 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.

          const DofMap& dof_map = f_system.get_dof_map();

Get a constant reference to the Finite Element type for the first (and only) variable in the system.

          const FEType& fe_type = dof_map.variable_type(0);

For the admittance boundary condition, get the fluid density

          const Real rho = es.parameters.get<Real>("rho");

In here, we will add the element matrices to the additional matrices "stiffness_mass", "damping", and the additional vector "rhs", not to the members "matrix" and "rhs". Therefore, get writable references to them

          SparseMatrix<Number>&   stiffness      = f_system.get_matrix("stiffness");
          SparseMatrix<Number>&   damping        = f_system.get_matrix("damping");
          SparseMatrix<Number>&   mass           = f_system.get_matrix("mass");
          NumericVector<Number>&  freq_indep_rhs = f_system.get_vector("rhs");

Some solver packages (PETSc) are especially picky about allocating sparsity structure and truly assigning values to this structure. Namely, matrix additions, as performed later, exhibit acceptable performance only for identical sparsity structures. Therefore, explicitly zero the values in the collective matrix, so that matrix additions encounter identical sparsity structures.

          SparseMatrix<Number>&  matrix           = *f_system.matrix;

------------------------------------------------------------------ Finite Element related stuff

Build a Finite Element object of the specified type. Since the FEBase::build() member dynamically creates memory we will store the object as an AutoPtr. This can be thought of as a pointer that will clean up after itself.

          AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));

A 5th order Gauss quadrature rule for numerical integration.

          QGauss qrule (dim, FIFTH);

Tell the finite element object to use our quadrature rule.

          fe->attach_quadrature_rule (&qrule);

The element Jacobian// quadrature weight at each integration point.

          const std::vector<Real>& JxW = fe->get_JxW();

The element shape functions evaluated at the quadrature points.

          const std::vector<std::vector<Real> >& phi = fe->get_phi();

The element shape function gradients evaluated at the quadrature points.

          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();

Now this is slightly different from example 4. We will not add directly to the overall (PETSc/LASPACK) matrix, but to the additional matrices "stiffness_mass" and "damping". The same holds for the right-hand-side vector Fe, which we will later on store in the additional vector "rhs". The zero_matrix is used to explicitly induce the same sparsity structure in the overall matrix. see later on. (At least) the mass, and stiffness matrices, however, are inherently real. Therefore, store these as one complex matrix. This will definitely save memory.

          DenseMatrix<Number> Ke, Ce, Me, zero_matrix;
          DenseVector<Number> Fe;

This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.

          std::vector<unsigned int> dof_indices;

Now we will loop over all the elements in the mesh. We will compute the element matrix and right-hand-side contribution.

          MeshBase::const_element_iterator           el = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
          
          for ( ; el != end_el; ++el)
            {

Start logging the element initialization.

              START_LOG("elem init","assemble_helmholtz");

Store a pointer to the element we are currently working on. This allows for nicer syntax later.

              const Elem* elem = *el;

Get the degree of freedom indices for the current element. These define where in the global matrix and right-hand-side this element will contribute to.

              dof_map.dof_indices (elem, dof_indices);

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->reinit (elem);

Zero & resize the element matrix and right-hand side before summing them, with different element types in the mesh this is quite necessary.

              {
                const unsigned int n_dof_indices = dof_indices.size();
        
                Ke.resize          (n_dof_indices, n_dof_indices);
                Ce.resize          (n_dof_indices, n_dof_indices);
                Me.resize          (n_dof_indices, n_dof_indices);
                zero_matrix.resize (n_dof_indices, n_dof_indices);
                Fe.resize          (n_dof_indices);
              }

Stop logging the element initialization.

              STOP_LOG("elem init","assemble_helmholtz");

Now loop over the quadrature points. This handles the numeric integration.

              START_LOG("stiffness & mass","assemble_helmholtz");
        
              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {

Now we will build the element matrix. This involves a double loop to integrate the test funcions (i) against the trial functions (j). Note the braces on the rhs of Ke(i,j): these are quite necessary to finally compute Real*(Point*Point) = Real, and not something else...

                  for (unsigned int i=0; i<phi.size(); i++)
                    for (unsigned int j=0; j<phi.size(); j++)
                      {
                        Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                        Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
                      }          
                }
        
              STOP_LOG("stiffness & mass","assemble_helmholtz");

Now compute the contribution to the element matrix (due to mixed boundary conditions) if the current element lies on the boundary.

The following loops over the sides of the element. If the element has no neighbor on a side then that side MUST live on a boundary of the domain.

              for (unsigned int side=0; side<elem->n_sides(); side++)
                if (elem->neighbor(side) == NULL)
                  {
                    START_LOG("damping","assemble_helmholtz");

Declare a special finite element object for boundary integration.

                    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));

Boundary integration requires one quadraure rule, with dimensionality one less than the dimensionality of the element.

                    QGauss qface(dim-1, SECOND);

Tell the finte element object to use our quadrature rule.

                    fe_face->attach_quadrature_rule (&qface);

The value of the shape functions at the quadrature points.

                    const std::vector<std::vector<Real> >&  phi_face =
                      fe_face->get_phi();

The Jacobian// Quadrature Weight at the quadrature points on the face.

                    const std::vector<Real>& JxW_face = fe_face->get_JxW();

Compute the shape function values on the element face.

                    fe_face->reinit(elem, side);

For the Robin BCs consider a normal admittance an=1 at some parts of the bounfdary

                    const Real an_value = 1.;

Loop over the face quadrature points for integration.

                    for (unsigned int qp=0; qp<qface.n_points(); qp++)
                      {

Element matrix contributrion due to precribed admittance boundary conditions.

                        for (unsigned int i=0; i<phi_face.size(); i++)
                          for (unsigned int j=0; j<phi_face.size(); j++)
                            Ce(i,j) += rho*an_value*JxW_face[qp]*phi_face[i][qp]*phi_face[j][qp];
                      }
        
                    STOP_LOG("damping","assemble_helmholtz");
                  }

If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations by uncommenting the following lines: std::vector dof_indicesC = dof_indices; std::vector dof_indicesM = dof_indices; dof_map.constrain_element_matrix (Ke, dof_indices); dof_map.constrain_element_matrix (Ce, dof_indicesC); dof_map.constrain_element_matrix (Me, dof_indicesM);

Finally, simply add the contributions to the additional matrices and vector.

              stiffness.add_matrix      (Ke, dof_indices);
              damping.add_matrix        (Ce, dof_indices);
              mass.add_matrix           (Me, dof_indices);

For the overall matrix, explicitly zero the entries where we added values in the other ones, so that we have identical sparsity footprints.

              matrix.add_matrix(zero_matrix, dof_indices);
        
            } // loop el

It now remains to compute the rhs. Here, we simply get the normal velocities values on the boundary from the mesh data.

          {
            START_LOG("rhs","assemble_helmholtz");

get a reference to the mesh data.

            const MeshData& mesh_data = es.get_mesh_data();

We will now loop over all nodes. In case nodal data for a certain node is available in the MeshData, we simply adopt the corresponding value for the rhs vector. Note that normal data was given in the mesh data file, i.e. one value per node

            libmesh_assert(mesh_data.n_val_per_node() == 1);
        
            MeshBase::const_node_iterator       node_it  = mesh.nodes_begin();
            const MeshBase::const_node_iterator node_end = mesh.nodes_end();
            
            for ( ; node_it != node_end; ++node_it)
              {

the current node pointer

                const Node* node = *node_it;

check if the mesh data has data for the current node and do for all components

                if (mesh_data.has_data(node))
                  for (unsigned int comp=0; comp<node->n_comp(0,0); comp++)
                    {

the dof number

                      unsigned int dn = node->dof_number(0,0,comp);

set the nodal value

                      freq_indep_rhs.set(dn, mesh_data.get_data(node)[0]);
                    }
              }
         
        
            STOP_LOG("rhs","assemble_helmholtz");
          }

All done!

We now define the function which will combine the previously-assembled mass, stiffness, and damping matrices into a single system matrix.

        void add_M_C_K_helmholtz(EquationSystems& es,
                                 const std::string& system_name)
        {
          START_LOG("init phase","add_M_C_K_helmholtz");

It is a good idea to make sure we are assembling the proper system.

          libmesh_assert (system_name == "Helmholtz");

Get a reference to our system, as before

          FrequencySystem & f_system =
            es.get_system<FrequencySystem> (system_name);

Get the frequency, fluid density, and speed of sound for which we should currently solve

          const Real frequency = es.parameters.get<Real> ("current frequency");
          const Real rho       = es.parameters.get<Real> ("rho");
          const Real speed     = es.parameters.get<Real> ("wave speed");

Compute angular frequency omega and wave number k

          const Real omega = 2.0*libMesh::pi*frequency;
          const Real k     = omega / speed;

Get writable references to the overall matrix and vector, where the frequency-dependent system is to be collected

          SparseMatrix<Number>&  matrix          = *f_system.matrix;
          NumericVector<Number>& rhs             = *f_system.rhs;

Get writable references to the frequency-independent matrices and rhs, though we only need to extract values. This write access is necessary, since solver packages have to close the data structure before they can extract values for computation.

          SparseMatrix<Number>&   stiffness      = f_system.get_matrix("stiffness");
          SparseMatrix<Number>&   damping        = f_system.get_matrix("damping");
          SparseMatrix<Number>&   mass           = f_system.get_matrix("mass");
          NumericVector<Number>&  freq_indep_rhs = f_system.get_vector("rhs");

form the scaling values for the coming matrix and vector axpy's

          const Number scale_stiffness (  1., 0.   );
          const Number scale_damping   (  0., omega);
          const Number scale_mass      (-k*k, 0.   );
          const Number scale_rhs       (  0., -(rho*omega));

Now simply add the matrices together, store the result in matrix and rhs. Clear them first.

          matrix.close(); matrix.zero ();
          rhs.close();    rhs.zero    ();

The matrices from which values are added to another matrix have to be closed. The add() method does take care of that, but let us do it explicitly.

          stiffness.close ();
          damping.close   ();
          mass.close      ();
        
          STOP_LOG("init phase","add_M_C_K_helmholtz");
        
          START_LOG("global matrix & vector additions","add_M_C_K_helmholtz");

add the stiffness and mass with the proper frequency to the overall system. For this to work properly, matrix has to be not only initialized, but filled with the identical sparsity structure as the matrix added to it, otherwise solver packages like PETSc crash.

Note that we have to add the mass and stiffness contributions one at a time; otherwise, the real part of matrix would be fine, but the imaginary part cluttered with unwanted products.

          matrix.add (scale_stiffness, stiffness);
          matrix.add (scale_mass,      mass);
          matrix.add (scale_damping,   damping);
          rhs.add    (scale_rhs,       freq_indep_rhs);
        
          STOP_LOG("global matrix & vector additions","add_M_C_K_helmholtz");

The "matrix" and "rhs" are now ready for solution

        }
        
        #endif // LIBMESH_USE_COMPLEX_NUMBERS

The program without comments:

 
   
  #include <iostream>
  #include <algorithm>
  #include <stdio.h>
  
  #include "libmesh.h"
  #include "libmesh_logging.h"
  #include "mesh.h"
  #include "mesh_generation.h"
  #include "exodusII_io.h"
  #include "equation_systems.h"
  #include "elem.h"
  
  #include "frequency_system.h"
  
  #include "fe.h"
  
  #include "quadrature_gauss.h"
  
  #include "dense_matrix.h"
  #include "dense_vector.h"
  
  #include "sparse_matrix.h"
  #include "numeric_vector.h"
  
  #include "dof_map.h"
  
  #include "mesh_data.h"
  
  using namespace libMesh;
  
  
  #ifdef LIBMESH_USE_COMPLEX_NUMBERS
  void assemble_helmholtz(EquationSystems& es,
                          const std::string& system_name);
  
  void add_M_C_K_helmholtz(EquationSystems& es,
                           const std::string& system_name);
  #endif
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
    
  #ifndef LIBMESH_USE_COMPLEX_NUMBERS
    libmesh_example_assert(false, "--enable-complex");
  #else
    
    if (argc < 3)
      {
        if (libMesh::processor_id() == 0)
          std::cerr << "Usage: " << argv[0] << " -f [frequency]"
                    << std::endl;
        
        libmesh_error();
      }
  
    if (libMesh::n_processors() > 1)
      {
        if (libMesh::processor_id() == 0)
          {
            std::cerr << "ERROR: Skipping example 7. " << std::endl;
            std::cerr << "MeshData objects currently only work in serial." << std::endl;
          }
        return 0;
      }
    
    else 
      {
        std::cout << "Running " << argv[0];
        
        for (int i=1; i<argc; i++)
          std::cout << " " << argv[i];
        
        std::cout << std::endl << std::endl;
      }
    
    const Real frequency_in = atof(argv[2]);
    const unsigned int n_frequencies = 3;        
    
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
    
    Mesh mesh;
  
    MeshData mesh_data(mesh);
    mesh_data.activate();
  
    mesh.read("lshape.unv", &mesh_data);
  
    mesh.print_info();    
  
    mesh_data.read("lshape_data.unv");
  
    mesh_data.print_info();
  
    EquationSystems equation_systems (mesh, &mesh_data);
    
    FrequencySystem & f_system =      
      equation_systems.add_system<FrequencySystem> ("Helmholtz");
    
    f_system.add_variable("p", SECOND);
    
    f_system.attach_assemble_function (assemble_helmholtz);
    f_system.attach_solve_function    (add_M_C_K_helmholtz);
    
    f_system.add_matrix ("stiffness");
    f_system.add_matrix ("damping");
    f_system.add_matrix ("mass");
    f_system.add_vector ("rhs");
    
    f_system.set_frequencies_by_steps (frequency_in/n_frequencies,
                                       frequency_in,
                                       n_frequencies);
    
    equation_systems.parameters.set<Real> ("wave speed") = 1.;
    equation_systems.parameters.set<Real> ("rho")        = 1.;
    
    equation_systems.init ();
    
    equation_systems.print_info ();
  
    for (unsigned int n=0; n < n_frequencies; n++)
      {
        f_system.solve (n,n);
        
  #ifdef LIBMESH_HAVE_EXODUS_API
        char buf[14];
        sprintf (buf, "out%04d.exd", n);
  
        ExodusII_IO(mesh).write_equation_systems (buf,
                                            equation_systems);
  #endif
      }
    
  
    equation_systems.write ("eqn_sys.dat", libMeshEnums::WRITE);
    
    return 0;
  
  #endif 
  }
  
  
  #ifdef LIBMESH_USE_COMPLEX_NUMBERS
  void assemble_helmholtz(EquationSystems& es,
                          const std::string& system_name)
  {
      
    libmesh_assert (system_name == "Helmholtz");
    
    const MeshBase& mesh = es.get_mesh();
    
    const unsigned int dim = mesh.mesh_dimension();
    
    FrequencySystem & f_system =
      es.get_system<FrequencySystem> (system_name);
    
    const DofMap& dof_map = f_system.get_dof_map();
    
    const FEType& fe_type = dof_map.variable_type(0);
  
    const Real rho = es.parameters.get<Real>("rho");
    
    SparseMatrix<Number>&   stiffness      = f_system.get_matrix("stiffness");
    SparseMatrix<Number>&   damping        = f_system.get_matrix("damping");
    SparseMatrix<Number>&   mass           = f_system.get_matrix("mass");
    NumericVector<Number>&  freq_indep_rhs = f_system.get_vector("rhs");
    
    SparseMatrix<Number>&  matrix           = *f_system.matrix;
    
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    
    QGauss qrule (dim, FIFTH);
  
    fe->attach_quadrature_rule (&qrule);
    
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
    
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke, Ce, Me, zero_matrix;
    DenseVector<Number> Fe;
    
    std::vector<unsigned int> dof_indices;
  
  
    MeshBase::const_element_iterator           el = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
    
    for ( ; el != end_el; ++el)
      {
        START_LOG("elem init","assemble_helmholtz");
        
        const Elem* elem = *el;
        
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
        
        {
          const unsigned int n_dof_indices = dof_indices.size();
  
          Ke.resize          (n_dof_indices, n_dof_indices);
          Ce.resize          (n_dof_indices, n_dof_indices);
          Me.resize          (n_dof_indices, n_dof_indices);
          zero_matrix.resize (n_dof_indices, n_dof_indices);
          Fe.resize          (n_dof_indices);
        }
        
        STOP_LOG("elem init","assemble_helmholtz");
  
        START_LOG("stiffness & mass","assemble_helmholtz");
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            for (unsigned int i=0; i<phi.size(); i++)
              for (unsigned int j=0; j<phi.size(); j++)
                {
                  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                  Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
                }          
          }
  
        STOP_LOG("stiffness & mass","assemble_helmholtz");
  
           
        for (unsigned int side=0; side<elem->n_sides(); side++)
          if (elem->neighbor(side) == NULL)
            {
              START_LOG("damping","assemble_helmholtz");
                
              AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
                
              QGauss qface(dim-1, SECOND);
                
              fe_face->attach_quadrature_rule (&qface);
                
              const std::vector<std::vector<Real> >&  phi_face =
                fe_face->get_phi();
                
              const std::vector<Real>& JxW_face = fe_face->get_JxW();
                
              fe_face->reinit(elem, side);
                
              const Real an_value = 1.;
                
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  
                  for (unsigned int i=0; i<phi_face.size(); i++)
                    for (unsigned int j=0; j<phi_face.size(); j++)
                      Ce(i,j) += rho*an_value*JxW_face[qp]*phi_face[i][qp]*phi_face[j][qp];
                }
  
              STOP_LOG("damping","assemble_helmholtz");
            }
  
  
       
        stiffness.add_matrix      (Ke, dof_indices);
        damping.add_matrix        (Ce, dof_indices);
        mass.add_matrix           (Me, dof_indices);
        
        matrix.add_matrix(zero_matrix, dof_indices);
  
      } // loop el
  
  
    {
      START_LOG("rhs","assemble_helmholtz");
      
      const MeshData& mesh_data = es.get_mesh_data();
  
      libmesh_assert(mesh_data.n_val_per_node() == 1);
  
      MeshBase::const_node_iterator       node_it  = mesh.nodes_begin();
      const MeshBase::const_node_iterator node_end = mesh.nodes_end();
      
      for ( ; node_it != node_end; ++node_it)
        {
          const Node* node = *node_it;
  
          if (mesh_data.has_data(node))
            for (unsigned int comp=0; comp<node->n_comp(0,0); comp++)
              {
                unsigned int dn = node->dof_number(0,0,comp);
  
                freq_indep_rhs.set(dn, mesh_data.get_data(node)[0]);
              }
        }
   
  
      STOP_LOG("rhs","assemble_helmholtz");
    }
  
  }
  
  
  void add_M_C_K_helmholtz(EquationSystems& es,
                           const std::string& system_name)
  {
    START_LOG("init phase","add_M_C_K_helmholtz");
    
    libmesh_assert (system_name == "Helmholtz");
    
    FrequencySystem & f_system =
      es.get_system<FrequencySystem> (system_name);
    
    const Real frequency = es.parameters.get<Real> ("current frequency");
    const Real rho       = es.parameters.get<Real> ("rho");
    const Real speed     = es.parameters.get<Real> ("wave speed");
    
    const Real omega = 2.0*libMesh::pi*frequency;
    const Real k     = omega / speed;
    
    SparseMatrix<Number>&  matrix          = *f_system.matrix;
    NumericVector<Number>& rhs             = *f_system.rhs;
    
    SparseMatrix<Number>&   stiffness      = f_system.get_matrix("stiffness");
    SparseMatrix<Number>&   damping        = f_system.get_matrix("damping");
    SparseMatrix<Number>&   mass           = f_system.get_matrix("mass");
    NumericVector<Number>&  freq_indep_rhs = f_system.get_vector("rhs");
    
    const Number scale_stiffness (  1., 0.   );
    const Number scale_damping   (  0., omega);
    const Number scale_mass      (-k*k, 0.   );
    const Number scale_rhs       (  0., -(rho*omega));
    
    matrix.close(); matrix.zero ();
    rhs.close();    rhs.zero    ();
    
    stiffness.close ();
    damping.close   ();
    mass.close      ();
  
    STOP_LOG("init phase","add_M_C_K_helmholtz");
  
    START_LOG("global matrix & vector additions","add_M_C_K_helmholtz");
    
    matrix.add (scale_stiffness, stiffness);
    matrix.add (scale_mass,      mass);
    matrix.add (scale_damping,   damping);
    rhs.add    (scale_rhs,       freq_indep_rhs);
  
    STOP_LOG("global matrix & vector additions","add_M_C_K_helmholtz");
    
  }
  
  #endif // LIBMESH_USE_COMPLEX_NUMBERS

The console output of the program:

Linking miscellaneous_ex2-opt... *************************************************************** * Running Example mpirun -np 6 ./miscellaneous_ex2-opt -f .5 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary *************************************************************** Assertion `false' failed. Configuring libMesh with --enable-complex may be required to run this code. ************************************************************************************************************************ *** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document *** ************************************************************************************************************************ ---------------------------------------------- PETSc Performance Summary: ---------------------------------------------- ./miscellaneous_ex2-opt on a intel-11. named daedalus with 6 processors, by roystgnr Fri Aug 24 15:21:46 2012 Using Petsc Release Version 3.1.0, Patch 5, Mon Sep 27 11:51:54 CDT 2010 Max Max/Min Avg Total Time (sec): 1.817e-03 2.14314 1.187e-03 Objects: 0.000e+00 0.00000 0.000e+00 Flops: 0.000e+00 0.00000 0.000e+00 0.000e+00 Flops/sec: 0.000e+00 0.00000 0.000e+00 0.000e+00 MPI Messages: 0.000e+00 0.00000 0.000e+00 0.000e+00 MPI Message Lengths: 0.000e+00 0.00000 0.000e+00 0.000e+00 MPI Reductions: 0.000e+00 0.00000 Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract) e.g., VecAXPY() for real vectors of length N --> 2N flops and VecAXPY() for complex vectors of length N --> 8N flops Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions -- Avg %Total Avg %Total counts %Total Avg %Total counts %Total 0: Main Stage: 1.1434e-03 96.4% 0.0000e+00 0.0% 0.000e+00 0.0% 0.000e+00 0.0% 0.000e+00 0.0% ------------------------------------------------------------------------------------------------------------------------ See the 'Profiling' chapter of the users' manual for details on interpreting output. Phase summary info: Count: number of times phase was executed Time and Flops: Max - maximum over all processors Ratio - ratio of maximum to minimum over all processors Mess: number of messages sent Avg. len: average message length Reduct: number of global reductions Global: entire computation Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop(). %T - percent time in this phase %F - percent flops in this phase %M - percent messages in this phase %L - percent message lengths in this phase %R - percent reductions in this phase Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors) ------------------------------------------------------------------------------------------------------------------------ Event Count Time (sec) Flops --- Global --- --- Stage --- Total Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %F %M %L %R %T %F %M %L %R Mflop/s ------------------------------------------------------------------------------------------------------------------------ --- Event Stage 0: Main Stage ------------------------------------------------------------------------------------------------------------------------ Memory usage is given in bytes: Object Type Creations Destructions Memory Descendants' Mem. Reports information only for process 0. --- Event Stage 0: Main Stage ======================================================================================================================== Average time to get PetscTime(): 0 Average time for MPI_Barrier(): 3.06129e-05 Average time for zero size MPI_Send(): 7.51416e-05 #PETSc Option Table entries: -f .5 -ksp_right_pc -log_summary -pc_type bjacobi -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -sub_pc_type ilu #End of PETSc Option Table entries Compiled without FORTRAN kernels Compiled with full precision matrices (default) sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 Configure run at: Sat May 19 03:47:23 2012 Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared=1 --with-shared-libraries=1 --with-mpi-dir=/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid --with-mumps=true --download-mumps=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t/libmkl_intel_lp64.so,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t/libmkl_sequential.so,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t/libmkl_core.so]" --with-lapack-lib=/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t/libmkl_solver_lp64_sequential.a ----------------------------------------- Libraries compiled on Sat May 19 03:47:23 CDT 2012 on daedalus Machine characteristics: Linux daedalus 2.6.32-34-generic #76-Ubuntu SMP Tue Aug 30 17:05:01 UTC 2011 x86_64 GNU/Linux Using PETSc directory: /org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5 Using PETSc arch: intel-11.1-lucid-mpich2-1.4.1-cxx-opt ----------------------------------------- Using C compiler: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/bin/mpicxx -O3 -fPIC Using Fortran compiler: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/bin/mpif90 -fPIC -O3 ----------------------------------------- Using include paths: -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/include -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/include -I/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/include -I/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/include ------------------------------------------ Using C linker: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/bin/mpicxx -O3 Using Fortran linker: /org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/bin/mpif90 -fPIC -O3 Using libraries: -Wl,-rpath,/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/lib -L/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/lib -L/org/centers/pecos/LIBRARIES/PETSC3/petsc-3.1-p5/intel-11.1-lucid-mpich2-1.4.1-cxx-opt/lib -lHYPRE -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lscalapack -lblacs -lsuperlu_dist_2.4 -lparmetis -lmetis -lsuperlu_4.0 -Wl,-rpath,/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t -L/org/centers/pecos/LIBRARIES/MKL/mkl-10.0.3.020-intel-11.1-lucid/lib/em64t -lmkl_solver_lp64_sequential -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -ldl -Wl,-rpath,/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/lib -L/org/centers/pecos/LIBRARIES/MPICH2/mpich2-1.4.1-intel-11.1-lucid/lib -lmpich -lopa -lmpl -lrt -lpthread -Wl,-rpath,/opt/intel/Compiler/11.1/073/lib/intel64 -L/opt/intel/Compiler/11.1/073/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-linux-gnu/4.4.3 -L/usr/lib/gcc/x86_64-linux-gnu/4.4.3 -limf -lsvml -lipgo -ldecimal -lgcc_s -lirc -lirc_s -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -lstdc++ -lmpichcxx -lstdc++ -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lgcc_s -lirc -lirc_s -ldl ------------------------------------------ *************************************************************** * Done Running Example mpirun -np 6 ./miscellaneous_ex2-opt -f .5 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary ***************************************************************