ex8.c

来自「一个用来实现偏微分方程中网格的计算库」· C语言 代码 · 共 660 行 · 第 1/2 页

      // Write nodal results to file.  The results can then
      // be viewed with e.g. gnuplot (run gnuplot and type
      // 'plot "pressure_node.res" with lines' in the command line)
      res_out << t_time << "\t"
              << global_displacement[dof_no]
              << std::endl;
    }
#endif
  
  // All done.  
  return 0;
}

// This function assembles the system matrix and right-hand-side
// for our wave equation.
void assemble_wave(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 == "Wave");

  // Get a constant reference to the mesh object.
  const MeshBase& mesh = es.get_mesh();

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

  // Copy the speed of sound and fluid density
  // to a local variable.
  const Real speed = es.parameters.get<Real>("speed");
  const Real rho   = es.parameters.get<Real>("fluid density");

  // Get a reference to our system, as before.
  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = t_system.get_dof_map().variable_type(0);

  // In here, we will add the element matrices to the
  // @e additional matrices "stiffness_mass" and "damping"
  // and the additional vector "force", not to the members 
  // "matrix" and "rhs".  Therefore, get writable
  // references to them.
  SparseMatrix<Number>&   stiffness = t_system.get_matrix("stiffness");
  SparseMatrix<Number>&   damping   = t_system.get_matrix("damping");
  SparseMatrix<Number>&   mass      = t_system.get_matrix("mass");

  NumericVector<Number>&  force     = t_system.get_vector("force");

  // 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     = *t_system.matrix;
  DenseMatrix<Number>    zero_matrix;

  // Build a Finite Element object of the specified type.  Since the
  // \p FEBase::build() member dynamically creates memory we will
  // store the object as an \p 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 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);

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

  // A reference to the \p DofMap object for this system.  The \p DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.
  const DofMap& dof_map = t_system.get_dof_map();

  // The element mass, damping and stiffness matrices
  // and the element contribution to the rhs.
  DenseMatrix<Number>   Ke, Ce, 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_elem_iterator           el (mesh.elements_begin());
//   const const_elem_iterator end_el (mesh.elements_end());

  MeshBase::const_element_iterator       el     = mesh.elements_begin();
  const MeshBase::const_element_iterator end_el = mesh.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);

      // 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 the element matrices and rhs before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.  Note that this will be the case if the
      // element type is different (i.e. the last element was HEX8
      // and now have a PRISM6).
      {
        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);
      }

      // Now loop over the quadrature points.  This handles
      // the numeric integration.
      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).
          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]
                           *1./(speed*speed);
              } // end of the matrix summation loop          
        } // end of quadrature point loop

      // Now compute the contribution to the element matrix and the
      // right-hand-side vector if the current element lies on the
      // boundary. 
      {
        // In this example no natural boundary conditions will
        // be considered.  The code is left here so it can easily
        // be extended.
        // 
        // don't do this for any side
        for (unsigned int side=0; side<elem->n_sides(); side++)
          if (!true)
            // if (elem->neighbor(side) == NULL)
            {
              // 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);

              // Here we consider a normal acceleration acc_n=1 applied to
              // the whole boundary of our mesh.
              const Real acc_n_value = 1.0;
              
              // Loop over the face quadrature points for integration.
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  // Right-hand-side contribution due to prescribed
                  // normal acceleration.
                  for (unsigned int i=0; i<phi_face.size(); i++)
                    {
                      Fe(i) += acc_n_value*rho
                        *phi_face[i][qp]*JxW_face[qp];
                    }
                } // end face quadrature point loop          
            } // end if (elem->neighbor(side) == NULL)
        
        // In this example the Dirichlet boundary conditions will be 
        // imposed via panalty method after the
        // system is assembled.
        
      } // end boundary condition section          

      // 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);
      
      force.add_vector     (Fe, 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);
    
    } // end of element loop
  
  // All done!
  return;
}

// This function applies the initial conditions
void apply_initial(EquationSystems& es,
                   const std::string& system_name)
{
  // Get a reference to our system, as before
  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);
  
  // Numeric vectors for the pressure, velocity and acceleration
  // values.
  NumericVector<Number>&  pres_vec       = t_system.get_vector("displacement");
  NumericVector<Number>&  vel_vec        = t_system.get_vector("velocity");
  NumericVector<Number>&  acc_vec        = t_system.get_vector("acceleration");

  // Assume our fluid to be at rest, which would
  // also be the default conditions in class NewmarkSystem,
  // but let us do it explicetly here.
  pres_vec.zero();
  vel_vec.zero();
  acc_vec.zero();
}

// This function applies the Dirichlet boundary conditions
void fill_dirichlet_bc(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 == "Wave");

  // Get a reference to our system, as before.
  NewmarkSystem & t_system = es.get_system<NewmarkSystem> (system_name);

  // Get writable references to the overall matrix and vector.
  SparseMatrix<Number>&  matrix = *t_system.matrix;
  NumericVector<Number>& rhs    = *t_system.rhs;

  // Get a constant reference to the mesh object.
  const MeshBase& mesh = es.get_mesh();

  // Get \p libMesh's  \f$ \pi \f$
  const Real pi = libMesh::pi;

  // Ask the \p EquationSystems flag whether
  // we should do this also for the matrix
  const bool do_for_matrix =
    es.parameters.get<bool>("Newmark set BC for Matrix");

  // Ge the current time from \p EquationSystems
  const Real current_time = es.parameters.get<Real>("time");

  // Number of nodes in the mesh.
  unsigned int n_nodes = mesh.n_nodes();

  for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
    {
      // Get a reference to the current node.
      const Node& curr_node = mesh.node(n_cnt);

      // Check if Dirichlet BCs should be applied to this node.
      // Use the \p TOLERANCE from \p mesh_common.h as tolerance.
      // Here a pressure value is applied if the z-coord.
      // is equal to 4, which corresponds to one end of the
      // pipe-mesh in this directory.
      const Real z_coo = 4.;

      if (fabs(curr_node(2)-z_coo) < TOLERANCE)
        {
          // The global number of the respective degree of freedom.
          unsigned int dn = curr_node.dof_number(0,0,0);

          // The penalty parameter.
          const Real penalty = 1.e10;

          // Here we apply sinusoidal pressure values for 0<t<0.002
          // at one end of the pipe-mesh.
          Real p_value;
          if (current_time < .002 )
            p_value = sin(2*pi*current_time/.002);
          else
            p_value = .0;

          // Now add the contributions to the matrix and the rhs.
          rhs.add(dn, p_value*penalty);

          // Add the panalty parameter to the global matrix
          // if desired.
          if (do_for_matrix)
            matrix.add(dn, dn, penalty);
        }
    } // loop n_cnt
}