ex6.c

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

  // A 2nd order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, SECOND);
  
  // Tell the finite element object to use our quadrature rule.   
  fe->attach_quadrature_rule (&qrule);
  
  // Due to its internal structure, the infinite element handles 
  // quadrature rules differently.  It takes the quadrature
  // rule which has been initialized for the FE object, but
  // creates suitable quadrature rules by @e itself.  The user
  // need not worry about this.   
  inf_fe->attach_quadrature_rule (&qrule);
  
  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke",  "Ce", "Me", and "Fe" for the stiffness, damping
  // and mass matrices, and the load vector.  Note that in 
  // Acoustics, these descriptors though do @e not match the 
  // true physical meaning of the projectors.  The final 
  // overall system, however, resembles the conventional 
  // notation again.   
  DenseMatrix<Number> Ke;
  DenseMatrix<Number> Ce;
  DenseMatrix<Number> Me;
  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.
  //   const_local_elem_iterator           el (mesh.elements_begin());
  //   const const_local_elem_iterator end_el (mesh.elements_end());

  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)
    {      
      // 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);

      
      // The mesh contains both finite and infinite elements.  These
      // elements are handled through different classes, namely
      // \p FE and \p InfFE, respectively.  However, since both
      // are derived from \p FEBase, they share the same interface,
      // and overall burden of coding is @e greatly reduced through
      // using a pointer, which is adjusted appropriately to the
      // current element type.       
      FEBase* cfe=NULL;
      
      // This here is almost the only place where we need to
      // distinguish between finite and infinite elements.
      // For faster computation, however, different approaches
      // may be feasible.
      //
      // Up to now, we do not know what kind of element we
      // have.  Aske the element of what type it is:        
      if (elem->infinite())
        {           
          // We have an infinite element.  Let \p cfe point
          // to our \p InfFE object.  This is handled through
          // an AutoPtr.  Through the \p AutoPtr::get() we "borrow"
          // the pointer, while the \p  AutoPtr \p inf_fe is
          // still in charge of memory management.           
          cfe = inf_fe.get(); 
        }
      else
        {
          // This is a conventional finite element.  Let \p fe handle it.           
            cfe = fe.get();
          
          // Boundary conditions.
          // Here we just zero the rhs-vector. For natural boundary 
          // conditions check e.g. previous examples.           
          {              
            // Zero the RHS for this element.                
            Fe.resize (dof_indices.size());
            
            system.rhs->add_vector (Fe, dof_indices);
          } // end boundary condition section             
        } // else ( if (elem->infinite())) )

      // This is slightly different from the Poisson solver:
      // Since the finite element object may change, we have to
      // initialize the constant references to the data fields
      // each time again, when a new element is processed.
      //
      // The element Jacobian * quadrature weight at each integration point.          
      const std::vector<Real>& JxW = cfe->get_JxW();
      
      // The element shape functions evaluated at the quadrature points.       
      const std::vector<std::vector<Real> >& phi = cfe->get_phi();
      
      // The element shape function gradients evaluated at the quadrature
      // points.       
      const std::vector<std::vector<RealGradient> >& dphi = cfe->get_dphi();

      // The infinite elements need more data fields than conventional FE.  
      // These are the gradients of the phase term \p dphase, an additional 
      // radial weight for the test functions \p Sobolev_weight, and its
      // gradient.
      // 
      // Note that these data fields are also initialized appropriately by
      // the \p FE method, so that the weak form (below) is valid for @e both
      // finite and infinite elements.       
      const std::vector<RealGradient>& dphase  = cfe->get_dphase();
      const std::vector<Real>&         weight  = cfe->get_Sobolev_weight();
      const std::vector<RealGradient>& dweight = cfe->get_Sobolev_dweight();

      // Now this is all independent of whether we use an \p FE
      // or an \p InfFE.  Nice, hm? ;-)
      //
      // Compute the element-specific data, as described
      // in previous examples.       
      cfe->reinit (elem);
      
      // Zero the element matrices.  Boundary conditions were already
      // processed in the \p FE-only section, see above.       
      Ke.resize (dof_indices.size(), dof_indices.size());
      Ce.resize (dof_indices.size(), dof_indices.size());
      Me.resize (dof_indices.size(), dof_indices.size());
      
      // The total number of quadrature points for infinite elements
      // @e has to be determined in a different way, compared to
      // conventional finite elements.  This type of access is also
      // valid for finite elements, so this can safely be used
      // anytime, instead of asking the quadrature rule, as
      // seen in previous examples.       
      unsigned int max_qp = cfe->n_quadrature_points();
      
      // Loop over the quadrature points.        
      for (unsigned int qp=0; qp<max_qp; qp++)
        {          
          // Similar to the modified access to the number of quadrature 
          // points, the number of shape functions may also be obtained
          // in a different manner.  This offers the great advantage
          // of being valid for both finite and infinite elements.           
          const unsigned int n_sf = cfe->n_shape_functions();

          // Now we will build the element matrices.  Since the infinite
          // elements are based on a Petrov-Galerkin scheme, the
          // resulting system matrices are non-symmetric. The additional
          // weight, described before, is part of the trial space.
          //
          // For the finite elements, though, these matrices are symmetric
          // just as we know them, since the additional fields \p dphase,
          // \p weight, and \p dweight are initialized appropriately.
          //
          // test functions:    weight[qp]*phi[i][qp]
          // trial functions:   phi[j][qp]
          // phase term:        phase[qp]
          // 
          // derivatives are similar, but note that these are of type
          // Point, not of type Real.           
          for (unsigned int i=0; i<n_sf; i++)
            for (unsigned int j=0; j<n_sf; j++)
              {
                //         (ndt*Ht + nHt*d) * nH 
                Ke(i,j) +=
                  (                            //    (                         
                   (                           //      (                       
                    dweight[qp] * phi[i][qp]   //        Point * Real  = Point 
                    +                          //        +                     
                    dphi[i][qp] * weight[qp]   //        Point * Real  = Point 
                    ) * dphi[j][qp]            //      )       * Point = Real  
                   ) * JxW[qp];                //    )         * Real  = Real  

                // (d*Ht*nmut*nH - ndt*nmu*Ht*H - d*nHt*nmu*H)
                Ce(i,j) +=
                  (                                //    (                         
                   (dphase[qp] * dphi[j][qp])      //      (Point * Point) = Real  
                   * weight[qp] * phi[i][qp]       //      * Real * Real   = Real  
                   -                               //      -                       
                   (dweight[qp] * dphase[qp])      //      (Point * Point) = Real  
                   * phi[i][qp] * phi[j][qp]       //      * Real * Real   = Real  
                   -                               //      -                       
                   (dphi[i][qp] * dphase[qp])      //      (Point * Point) = Real  
                   * weight[qp] * phi[j][qp]       //      * Real * Real   = Real  
                   ) * JxW[qp];                    //    )         * Real  = Real  
                
                // (d*Ht*H * (1 - nmut*nmu))
                Me(i,j) +=
                  (                                       //    (                                  
                   (1. - (dphase[qp] * dphase[qp]))       //      (Real  - (Point * Point)) = Real 
                   * phi[i][qp] * phi[j][qp] * weight[qp] //      * Real *  Real  * Real    = Real 
                   ) * JxW[qp];                           //    ) * Real                    = Real 

              } // end of the matrix summation loop
        } // end of quadrature point loop

      // The element matrices are now built for this element.  
      // Collect them in Ke, and then add them to the global matrix.  
      // The \p SparseMatrix::add_matrix() member does this for us.
      Ke.add(1./speed        , Ce);
      Ke.add(1./(speed*speed), Me);

      system.matrix->add_matrix (Ke, dof_indices);
    } // end of element loop

  // Note that we have not applied any boundary conditions so far.
  // Here we apply a unit load at the node located at (0,0,0).
  {
    // Iterate over local nodes
    MeshBase::const_node_iterator           nd = mesh.local_nodes_begin();
    const MeshBase::const_node_iterator nd_end = mesh.local_nodes_end();
    
    for (; nd != nd_end; ++nd)
      {        
        // Get a reference to the current node.
        const Node& node = **nd;
        
        // Check the location of the current node.
        if (fabs(node(0)) < TOLERANCE &&
            fabs(node(1)) < TOLERANCE &&
            fabs(node(2)) < TOLERANCE)
          {
            // The global number of the respective degree of freedom.
            unsigned int dn = node.dof_number(0,0,0);

            system.rhs->add (dn, 1.);
          }
      }
  }

#else

  // dummy assert 
  libmesh_assert(es.get_mesh().mesh_dimension() != 1);

#endif //ifdef ENABLE_INFINITE_ELEMENTS
  
  // All done!   
  return;
}