ex7.c

一个用来实现偏微分方程中网格的计算库

  // 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<FEBase>.  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.local_elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.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");
          }

     
      // 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 // USE_COMPLEX_NUMBERS