ex10.c

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

          // Assemble & solve the linear system
          system.solve();
          
          // Possibly refine the mesh
          if (r_step+1 != max_r_steps)
            {
              std::cout << "  Refining the mesh..." << std::endl;

              // The \p ErrorVector is a particular \p StatisticsVector
              // for computing error information on a finite element mesh.
              ErrorVector error;

              // The \p ErrorEstimator class interrogates a finite element
              // solution and assigns to each element a positive error value.
              // This value is used for deciding which elements to refine
              // and which to coarsen.
              //ErrorEstimator* error_estimator = new KellyErrorEstimator;
              KellyErrorEstimator error_estimator;
              
              // Compute the error for each active element using the provided
              // \p flux_jump indicator.  Note in general you will need to
              // provide an error estimator specifically designed for your
              // application.
              error_estimator.estimate_error (system,
                                              error);
              
              // This takes the error in \p error and decides which elements
              // will be coarsened or refined.  Any element within 20% of the
              // maximum error on any element will be refined, and any
              // element within 7% of the minimum error on any element might
              // be coarsened. Note that the elements flagged for refinement
              // will be refined, but those flagged for coarsening _might_ be
              // coarsened.
              mesh_refinement.refine_fraction() = 0.80;
              mesh_refinement.coarsen_fraction() = 0.07;
              mesh_refinement.max_h_level() = 5;
              mesh_refinement.flag_elements_by_error_fraction (error);
              
              // This call actually refines and coarsens the flagged
              // elements.
              mesh_refinement.refine_and_coarsen_elements();
              
              // This call reinitializes the \p EquationSystems object for
              // the newly refined mesh.  One of the steps in the
              // reinitialization is projecting the \p solution,
              // \p old_solution, etc... vectors from the old mesh to
              // the current one.
              equation_systems.reinit ();
            }            
        }
        
      // Output evey 10 timesteps to file.
      if ( (t_step+1)%10 == 0)
        {
          OStringStream file_name;

          file_name << "out.gmv.";
          OSSRealzeroright(file_name,3,0,t_step+1);

          GMVIO(mesh).write_equation_systems (file_name.str(),
                                              equation_systems);
        }
    }

  if(!read_solution)
    {
      mesh.write("saved_mesh.xda");
      equation_systems.write("saved_solution.xda", libMeshEnums::WRITE);
    }
#endif // #ifndef ENABLE_AMR
  
  return 0;
}

// Here we define the initialization routine for the
// Convection-Diffusion system.  This routine is
// responsible for applying the initial conditions to
// the system.
void init_cd (EquationSystems& es,
              const std::string& system_name)
{
  // It is a good idea to make sure we are initializing
  // the proper system.
  libmesh_assert (system_name == "Convection-Diffusion");

  // Get a reference to the Convection-Diffusion system object.
  TransientLinearImplicitSystem & system =
    es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");

  // Project initial conditions at time 0
  es.parameters.set<Real> ("time") = 0;
  
  system.project_solution(exact_value, NULL, es.parameters);
}



// This function defines the assembly routine which
// will be called at each time step.  It is responsible
// for computing the proper matrix entries for the
// element stiffness matrices and right-hand sides.
void assemble_cd (EquationSystems& es,
                  const std::string& system_name)
{
#ifdef ENABLE_AMR
  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert (system_name == "Convection-Diffusion");
  
  // 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();
  
  // Get a reference to the Convection-Diffusion system object.
  TransientLinearImplicitSystem & system =
    es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
  
  // Get the Finite Element type for the first (and only) 
  // variable in the system.
  FEType fe_type = system.variable_type(0);
  
  // 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));
  AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  
  // A Gauss quadrature rule for numerical integration.
  // Let the \p FEType object decide what order rule is appropriate.
  QGauss qrule (dim,   fe_type.default_quadrature_order());
  QGauss qface (dim-1, fe_type.default_quadrature_order());

  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule      (&qrule);
  fe_face->attach_quadrature_rule (&qface);

  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.  We will start
  // with the element Jacobian * quadrature weight at each integration point.
  const std::vector<Real>& JxW      = fe->get_JxW();
  const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real> >& phi = fe->get_phi();
  const std::vector<std::vector<Real> >& psi = fe_face->get_phi();

  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();

  // The XY locations of the quadrature points used for face integration
  const std::vector<Point>& qface_points = fe_face->get_xyz();
    
  // 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.  We will talk more about the \p DofMap
  // in future examples.
  const DofMap& dof_map = system.get_dof_map();
  
  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe".
  DenseMatrix<Number> Ke;
  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;

  // Here we extract the velocity & parameters that we put in the
  // EquationSystems object.
  const RealVectorValue velocity =
    es.parameters.get<RealVectorValue> ("velocity");

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

  // Now we will loop over all the elements in the mesh that
  // live on the local processor. We will compute the element
  // matrix and right-hand-side contribution.  Since the mesh
  // will be refined we want to only consider the ACTIVE elements,
  // hence we use a variant of the \p active_elem_iterator.
  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)
    {    
      // 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 matrix and right-hand side 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 a
      // triangle, now we are on a quadrilateral).
      Ke.resize (dof_indices.size(),
                 dof_indices.size());

      Fe.resize (dof_indices.size());
      
      // Now we will build the element matrix and right-hand-side.
      // Constructing the RHS requires the solution and its
      // gradient from the previous timestep.  This myst be
      // calculated at each quadrature point by summing the
      // solution degree-of-freedom values by the appropriate
      // weight functions.
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          // Values to hold the old solution & its gradient.
          Number   u_old = 0.;
          Gradient grad_u_old;
          
          // Compute the old solution & its gradient.
          for (unsigned int l=0; l<phi.size(); l++)
            {
              u_old      += phi[l][qp]*system.old_solution  (dof_indices[l]);
              
              // This will work,
              // grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]);
              // but we can do it without creating a temporary like this:
              grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
            }
          
          // Now compute the element matrix and RHS contributions.
          for (unsigned int i=0; i<phi.size(); i++)
            {
              // The RHS contribution
              Fe(i) += JxW[qp]*(
                                // Mass matrix term
                                u_old*phi[i][qp] + 
                                -.5*dt*(
                                        // Convection term
                                        // (grad_u_old may be complex, so the
                                        // order here is important!)
                                        (grad_u_old*velocity)*phi[i][qp] +
                                        
                                        // Diffusion term
                                        0.01*(grad_u_old*dphi[i][qp]))     
                                );
              
              for (unsigned int j=0; j<phi.size(); j++)
                {
                  // The matrix contribution
                  Ke(i,j) += JxW[qp]*(
                                      // Mass-matrix
                                      phi[i][qp]*phi[j][qp] + 
                                      .5*dt*(
                                             // Convection term
                                             (velocity*dphi[j][qp])*phi[i][qp] +
                                             // Diffusion term
                                             0.01*(dphi[i][qp]*dphi[j][qp]))      
                                      );
                }
            } 
        } 

      // At this point the interior element integration has
      // been completed.  However, we have not yet addressed
      // boundary conditions.  For this example we will only
      // consider simple Dirichlet boundary conditions imposed
      // via the penalty method. 
      //        
      // 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.
      {
        // The penalty value.  
        const Real penalty = 1.e10;

        // 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 s=0; s<elem->n_sides(); s++)
          if (elem->neighbor(s) == NULL)
            {
              fe_face->reinit(elem,s);
              
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  const Number value = exact_solution (qface_points[qp](0),
                                                       qface_points[qp](1),
                                                       time);
                                                       
                  // RHS contribution
                  for (unsigned int i=0; i<psi.size(); i++)
                    Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];

                  // Matrix contribution
                  for (unsigned int i=0; i<psi.size(); i++)
                    for (unsigned int j=0; j<psi.size(); j++)
                      Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
                }
            } 
      } 

      
      // We have now built the element matrix and RHS vector in terms
      // of the element degrees of freedom.  However, it is possible
      // that some of the element DOFs are constrained to enforce
      // solution continuity, i.e. they are not really "free".  We need
      // to constrain those DOFs in terms of non-constrained DOFs to
      // ensure a continuous solution.  The
      // \p DofMap::constrain_element_matrix_and_vector() method does
      // just that.
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
      
      // The element matrix and right-hand-side are now built
      // for this element.  Add them to the global matrix and
      // right-hand-side vector.  The \p PetscMatrix::add_matrix()
      // and \p PetscVector::add_vector() members do this for us.
      system.matrix->add_matrix (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);
      
    }
  // Finished computing the sytem matrix and right-hand side.
#endif // #ifdef ENABLE_AMR
}