ex15.c

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

{
  // Second derivatives to be returned.
  Tensor hessu;
  
  // x and y coordinates in space
  const Real x = p(0);
  const Real y = p(1);

  hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x)*(y-y*y)*(y-y*y);
  hessu(0,1) = 256.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y);
  hessu(1,1) = 256.*2.*(x-x*x)*(x-x*x)*(1.-6.*y+6.*y*y);

  // Hessians are always symmetric
  hessu(1,0) = hessu(0,1);
  return hessu;
}



Number forcing_function_2D(const Point& p)
{
  // x and y coordinates in space
  const Real x = p(0);
  const Real y = p(1);

  // Equals laplacian(laplacian(u))
  return 256. * 8. * (3.*((y-y*y)*(y-y*y)+(x-x*x)*(x-x*x))
         + (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y));
}



Number exact_3D_solution(const Point& p,
                         const Parameters&,  // parameters, not needed
                         const std::string&, // sys_name, not needed
                         const std::string&) // unk_name, not needed
{
  // xyz coordinates in space
  const Real x = p(0);
  const Real y = p(1);
  const Real z = p(2);
  
  // analytic solution value
  return 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
}


Gradient exact_3D_derivative(const Point& p,
                             const Parameters&,  // parameters, not needed
                             const std::string&, // sys_name, not needed
                             const std::string&) // unk_name, not needed
{
  // First derivatives to be returned.
  Gradient gradu;
  
  // xyz coordinates in space
  const Real x = p(0);
  const Real y = p(1);
  const Real z = p(2);

  gradu(0) = 4096.*2.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
  gradu(1) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
  gradu(2) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);

  return gradu;
}


// We now define the hessian of the exact solution
Tensor exact_3D_hessian(const Point& p,
                        const Parameters&,  // parameters, not needed
                        const std::string&, // sys_name, not needed
                        const std::string&) // unk_name, not needed
{
  // Second derivatives to be returned.
  Tensor hessu;
  
  // xyz coordinates in space
  const Real x = p(0);
  const Real y = p(1);
  const Real z = p(2);

  hessu(0,0) = 4096.*(2.-12.*x+12.*x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
  hessu(0,1) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
  hessu(0,2) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
  hessu(1,1) = 4096.*(x-x*x)*(x-x*x)*(2.-12.*y+12.*y*y)*(z-z*z)*(z-z*z);
  hessu(1,2) = 4096.*4.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(1.-2.*z);
  hessu(2,2) = 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(2.-12.*z+12.*z*z);

  // Hessians are always symmetric
  hessu(1,0) = hessu(0,1);
  hessu(2,0) = hessu(0,2);
  hessu(2,1) = hessu(1,2);

  return hessu;
}



Number forcing_function_3D(const Point& p)
{
  // xyz coordinates in space
  const Real x = p(0);
  const Real y = p(1);
  const Real z = p(2);

  // Equals laplacian(laplacian(u))
  return 4096. * 8. * (3.*((y-y*y)*(y-y*y)*(x-x*x)*(x-x*x) +
                           (z-z*z)*(z-z*z)*(x-x*x)*(x-x*x) +
                           (z-z*z)*(z-z*z)*(y-y*y)*(y-y*y)) +
         (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y)*(z-z*z)*(z-z*z) +
         (1.-6.*x+6.*x*x)*(1.-6.*z+6.*z*z)*(y-y*y)*(y-y*y) +
         (1.-6.*y+6.*y*y)*(1.-6.*z+6.*z*z)*(x-x*x)*(x-x*x));
}



// We now define the matrix assembly function for the
// Biharmonic 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_biharmonic(EquationSystems& es,
                      const std::string& system_name)
{
#ifdef ENABLE_AMR
#ifdef ENABLE_SECOND_DERIVATIVES

  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert (system_name == "Biharmonic");

  // 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>("Biharmonic");
  
  // 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));
  
  // Quadrature rule for numerical integration.
  // With 2D triangles, the Clough quadrature rule puts a Gaussian
  // quadrature rule on each of the 3 subelements
  AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(dim));

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

  // Declare a special finite element object for
  // boundary integration.
  AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
              
  // Boundary integration requires another quadraure rule,
  // with dimensionality one less than the dimensionality
  // of the element.
  // In 1D, the Clough and Gauss quadrature rules are identical.
  AutoPtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));

  // Tell the finte element object to use our
  // quadrature rule.
  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();

  // The physical XY locations of the quadrature points on the element.
  // These might be useful for evaluating spatially varying material
  // properties at the quadrature points.
  const std::vector<Point>& q_point = fe->get_xyz();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real> >& phi = fe->get_phi();

  // The element shape function second derivatives evaluated at the
  // quadrature points.  Note that for the simple biharmonic, shape
  // function first derivatives are unnecessary.
  const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi();

  // For efficiency we will compute shape function laplacians n times,
  // not n^2
  std::vector<Real> shape_laplacian;

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

  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.
      Ke.resize (dof_indices.size(),
                 dof_indices.size());

      Fe.resize (dof_indices.size());

      // Make sure there is enough room in this cache
      shape_laplacian.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 laplacians of the test funcions
      // (i) against laplacians of the trial functions (j).
      //
      // This step is why we need the Clough-Tocher elements -
      // these C1 differentiable elements have square-integrable
      // second derivatives.
      //
      // 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<phi.size(); i++)
            {
              shape_laplacian[i] = d2phi[i][qp](0,0)+d2phi[i][qp](1,1);
              if (dim == 3)
                 shape_laplacian[i] += d2phi[i][qp](2,2);
            }
          for (unsigned int i=0; i<phi.size(); i++)
            for (unsigned int j=0; j<phi.size(); j++)
              Ke(i,j) += JxW[qp]*
                         shape_laplacian[i]*shape_laplacian[j];
        }

      // 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.  Note that this is a fourth-order
      // problem: Dirichlet boundary conditions include *both*
      // boundary values and boundary normal fluxes.
      {
        // Start logging the boundary condition computation
        perf_log.push ("BCs");

        // The penalty values, for solution boundary trace and flux.  
        const Real penalty = 1e10;
        const Real penalty2 = 1e10;

        // 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)
            {
              // The value of the shape functions at the quadrature
              // points.
              const std::vector<std::vector<Real> >&  phi_face =
                              fe_face->get_phi();

              // The value of the shape function derivatives at the
              // quadrature points.
              const std::vector<std::vector<RealGradient> >& dphi_face =
                              fe_face->get_dphi();

              // The Jacobian * Quadrature Weight at the quadrature
              // points on the face.
              const std::vector<Real>& JxW_face = fe_face->get_JxW();
                                                                               
              // The XYZ locations (in physical space) of the
              // quadrature points on the face.  This is where
              // we will interpolate the boundary value function.
              const std::vector<Point>& qface_point = fe_face->get_xyz();

              const std::vector<Point>& face_normals =
                              fe_face->get_normals();

              // Compute the shape function values on the element
              // face.
              fe_face->reinit(elem, s);
                                                                                
              // Loop over the face quagrature points for integration.
              for (unsigned int qp=0; qp<qface->n_points(); qp++)
                {
                  // The boundary value.
                  Number value = exact_solution(qface_point[qp],
                                                es.parameters, "null",
                                                "void");
                  Gradient flux = exact_2D_derivative(qface_point[qp],
                                                      es.parameters,
                                                      "null", "void");

                  // Matrix contribution of the L2 projection.
                  // Note that the basis function values are
                  // integrated against test function values while
                  // basis fluxes are integrated against test function
                  // fluxes.
                  for (unsigned int i=0; i<phi_face.size(); i++)
                    for (unsigned int j=0; j<phi_face.size(); j++)
                      Ke(i,j) += JxW_face[qp] *
                                 (penalty * phi_face[i][qp] *
                                  phi_face[j][qp] + penalty2
                                  * (dphi_face[i][qp] *
                                  face_normals[qp]) *
                                  (dphi_face[j][qp] *
                                   face_normals[qp]));

                  // Right-hand-side contribution of the L2
                  // projection.
                  for (unsigned int i=0; i<phi_face.size(); i++)
                    Fe(i) += JxW_face[qp] *
                                    (penalty * value * phi_face[i][qp]
                                     + penalty2 * 
                                     (flux * face_normals[qp])
                                    * (dphi_face[i][qp]
                                       * face_normals[qp]));

                }
            } 
        
        // Stop logging the boundary condition computation
        perf_log.pop ("BCs");
      } 

      for (unsigned int qp=0; qp<qrule->n_points(); qp++)
        for (unsigned int i=0; i<phi.size(); i++)
          Fe(i) += JxW[qp]*phi[i][qp]*forcing_function(q_point[qp]);

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

      // Stop 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?

#else

#endif // #ifdef ENABLE_SECOND_DERIVATIVES
#endif // #ifdef ENABLE_AMR
}