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
***************************************************************

Site Created By: libMesh Developers
Last modified: December 11 2012 18:19:52 UTC

Hosted By:
SourceForge.net Logo