ex14.c

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

            }
        
          // 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 ();
        }
    }            
  
  // Write out the solution
  // After solving the system write the solution
  // to a GMV-formatted plot file.
  GMVIO (mesh).write_equation_systems ("lshaped.gmv",
                                       equation_systems);

  // Close up the output file.
  out << "];" << std::endl;
  out << "hold on" << std::endl;
  out << "plot(e(:,1), e(:,2), 'bo-');" << std::endl;
  out << "plot(e(:,1), e(:,3), 'ro-');" << std::endl;
  //    out << "set(gca,'XScale', 'Log');" << std::endl;
  //    out << "set(gca,'YScale', 'Log');" << std::endl;
  out << "xlabel('dofs');" << std::endl;
  out << "title('" << approx_name << " elements');" << std::endl;
  out << "legend('L2-error', 'H1-error');" << std::endl;
  //     out << "disp('L2-error linear fit');" << std::endl;
  //     out << "polyfit(log10(e(:,1)), log10(e(:,2)), 1)" << std::endl;
  //     out << "disp('H1-error linear fit');" << std::endl;
  //     out << "polyfit(log10(e(:,1)), log10(e(:,3)), 1)" << std::endl;
#endif // #ifndef ENABLE_AMR
  
  // All done.  
  return 0;
}




// We now define the exact solution, being careful
// to obtain an angle from atan2 in the correct
// quadrant.
Number exact_solution(const Point& p,
                      const Parameters&,  // parameters, not needed
                      const std::string&, // sys_name, not needed
                      const std::string&) // unk_name, not needed
{
  const Real x = p(0);
  const Real y = (dim > 1) ? p(1) : 0.;
  
  if (singularity)
    {
      // The exact solution to the singular problem,
      // u_exact = r^(2/3)*sin(2*theta/3).
      Real theta = atan2(y,x);

      // Make sure 0 <= theta <= 2*pi
      if (theta < 0)
        theta += 2. * libMesh::pi;

      // Make the 3D solution similar
      const Real z = (dim > 2) ? p(2) : 0;
                  
      return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
    }
  else
    {
      // The exact solution to a nonsingular problem,
      // good for testing ideal convergence rates
      const Real z = (dim > 2) ? p(2) : 0;

      return cos(x) * exp(y) * (1. - z);
    }
}





// We now define the gradient of the exact solution, again being careful
// to obtain an angle from atan2 in the correct
// quadrant.
Gradient exact_derivative(const Point& p,
                          const Parameters&,  // parameters, not needed
                          const std::string&, // sys_name, not needed
                          const std::string&) // unk_name, not needed
{
  // Gradient value to be returned.
  Gradient gradu;
  
  // x and y coordinates in space
  const Real x = p(0);
  const Real y = dim > 1 ? p(1) : 0.;

  if (singularity)
    {
      // We can't compute the gradient at x=0, it is not defined.
      libmesh_assert (x != 0.);

      // For convenience...
      const Real tt = 2./3.;
      const Real ot = 1./3.;
  
      // The value of the radius, squared
      const Real r2 = x*x + y*y;

      // The boundary value, given by the exact solution,
      // u_exact = r^(2/3)*sin(2*theta/3).
      Real theta = atan2(y,x);
  
      // Make sure 0 <= theta <= 2*pi
      if (theta < 0)
        theta += 2. * libMesh::pi;

      // du/dx
      gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;

      // du/dy
      if (dim > 1)
        gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;

      if (dim > 2)
        gradu(2) = 1.;
    }
  else
    {
      const Real z = (dim > 2) ? p(2) : 0;

      gradu(0) = -sin(x) * exp(y) * (1. - z);
      if (dim > 1)
        gradu(1) = cos(x) * exp(y) * (1. - z);
      if (dim > 2)
        gradu(2) = -cos(x) * exp(y);
    }

  return gradu;
}






// We now define the matrix assembly function for the
// Laplace system.  We need to first compute element
// matrices and right-hand sides, and then take into
// account the boundary conditions, which will be handled
// via a penalty method.
void assemble_laplace(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 == "Laplace");


  // Declare a performance log.  Give it a descriptive
  // string to identify what part of the code we are
  // logging, since there may be many PerfLogs in an
  // application.
  PerfLog perf_log ("Matrix Assembly",false);
  
    // 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 LinearImplicitSystem we are solving
  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Laplace");
  
  // 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();

  // Get a constant reference to the Finite Element type
  // for the first (and only) variable in the system.
  FEType fe_type = dof_map.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));
  
  // Quadrature rules for numerical integration.
  AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(dim));
  AutoPtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));

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

  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  // We begin 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 physical XY locations of the quadrature points on the element.
  // These might be useful for evaluating spatially varying material
  // properties or forcing functions at the quadrature points.
  const std::vector<Point>& q_point = fe->get_xyz();

  // The element shape functions evaluated at the quadrature points.
  // For this simple problem we usually only need them on element
  // boundaries.
  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();

  // 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". More detail is in example 3.
  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;

  // Now we will loop over all the elements in the mesh.  We will
  // compute the element matrix and right-hand-side contribution.  See
  // example 3 for a discussion of the element iterators.  Here we use
  // the \p const_active_local_elem_iterator to indicate we only want
  // to loop over elements that are assigned to the local processor
  // which are "active" in the sense of AMR.  This allows each
  // processor to compute its components of the global matrix for
  // active elements while ignoring parent elements which have been
  // refined.
  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 shape function initialization.
      // This is done through a simple function call with
      // the name of the event to log.
      perf_log.push("elem init");      

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

      // Stop logging the shape function initialization.
      // If you forget to stop logging an event the PerfLog
      // object will probably catch the error and abort.
      perf_log.pop("elem init");      

      // Now we will build the element matrix.  This involves
      // a double loop to integrate the test funcions (i) against
      // the trial functions (j).
      //
      // Now start logging the element matrix computation
      perf_log.push ("Ke");

      for (unsigned int qp=0; qp<qrule->n_points(); qp++)
        for (unsigned int i=0; i<dphi.size(); i++)
          for (unsigned int j=0; j<dphi.size(); j++)
            Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);

      // We need a forcing function to make the 1D case interesting
      if (dim == 1)
        for (unsigned int qp=0; qp<qrule->n_points(); qp++)
          {
            Real x = q_point[qp](0);
            Real f = singularity ? sqrt(3.)/9.*pow(-x, -4./3.) :
                                   cos(x);
            for (unsigned int i=0; i<dphi.size(); ++i)
              Fe(i) += JxW[qp]*phi[i][qp]*f;
          }

      // Stop logging the matrix computation
      perf_log.pop ("Ke");


      // 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.
      //
      // This approach adds the L2 projection of the boundary
      // data in penalty form to the weak statement.  This is
      // a more generic approach for applying Dirichlet BCs
      // which is applicable to non-Lagrange finite element
      // discretizations.
      {
        // Start logging the boundary condition computation
        perf_log.push ("BCs");

        // 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],
                                                       es.parameters,
                                                       "null",
                                                       "void");

                  // 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];
                }
            } 
        
        // Stop logging the boundary condition computation
        perf_log.pop ("BCs");
      } 
      

      // 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.
      // Start logging the insertion of the local (element)
      // matrix and vector into the global matrix and vector
      perf_log.push ("matrix insertion");

      dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
      system.matrix->add_matrix (Ke, dof_indices);
      system.rhs->add_vector    (Fe, dof_indices);

      // Start logging the insertion of the local (element)
      // matrix and vector into the global matrix and vector
      perf_log.pop ("matrix insertion");
    }

  // That's it.  We don't need to do anything else to the
  // PerfLog.  When it goes out of scope (at this function return)
  // it will print its log to the screen. Pretty easy, huh?
#endif // #ifdef ENABLE_AMR
}