fe.c

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

	   
	    if (Dim == 3)	     
	      dphidzeta[i].resize (n_qp);
	  }
#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  {
            d2phi[i].resize     (n_qp);
            d2phidx2[i].resize  (n_qp);
            d2phidxdy[i].resize (n_qp);
            d2phidxdz[i].resize (n_qp);
            d2phidy2[i].resize  (n_qp);
            d2phidydz[i].resize (n_qp);
            d2phidz2[i].resize  (n_qp);
            d2phidxi2[i].resize (n_qp);
            if (Dim > 1)
              {
                d2phidxideta[i].resize (n_qp);
                d2phideta2[i].resize   (n_qp);
              }
            if (Dim > 2)
              {
                d2phidxidzeta[i].resize  (n_qp);
                d2phidetadzeta[i].resize (n_qp);
                d2phidzeta2[i].resize    (n_qp);
              }
	  }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
      }
       
    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
      {
	phi_map[i].resize         (n_qp);
	dphidxi_map[i].resize     (n_qp);
#ifdef ENABLE_SECOND_DERIVATIVES
	d2phidxi2_map[i].resize     (n_qp);
	if (Dim > 1)
          {
            d2phidxideta_map[i].resize (n_qp);
            d2phideta2_map[i].resize (n_qp);
          }
        if (Dim > 2)
          {
            d2phidxidzeta_map[i].resize  (n_qp);
            d2phidetadzeta_map[i].resize (n_qp);
            d2phidzeta2_map[i].resize    (n_qp);
          }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	   
	if (Dim > 1)
	  dphideta_map[i].resize  (n_qp);
	   
	if (Dim == 3)
	  dphidzeta_map[i].resize (n_qp);
      }
  }


      
#ifdef ENABLE_INFINITE_ELEMENTS
  //------------------------------------------------------------
  // Initialize the data fields, which should only be used for infinite 
  // elements, to some sensible values, so that using a FE with the
  // variational formulation of an InfFE, correct element matrices are
  // returned

 {
    weight.resize  (n_qp);
    dweight.resize (n_qp);
    dphase.resize  (n_qp);
    
    for (unsigned int p=0; p<n_qp; p++)
      {
        weight[p] = 1.;
	dweight[p].zero();
	dphase[p].zero();
      }

 }
#endif // ifdef ENABLE_INFINITE_ELEMENTS

  // Optimize for the affine elements case:
  bool has_affine_map = elem->has_affine_map();
  
  switch (Dim)
    {

      //------------------------------------------------------------
      // 1D
    case 1:
      {
	// Compute the value of the approximation shape function i at quadrature point p
	if (calculate_phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      phi[i][p]      = FE<Dim,T>::shape       (elem, fe_type.order, i,    qp[p]);
	if (calculate_dphi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      dphidxi[i][p]  = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 0, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
              d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 0, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	
	// Compute the value of the mapping shape function i at quadrature point p
	// (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
	    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
	        phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
	        dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	        for (unsigned int p=1; p<n_qp; p++)
                  {
	            phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
	            dphidxi_map[i][p]  = dphidxi_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    d2phidxi2_map[i][p] = d2phidxi2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
	  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        phi_map[i][p]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
	        dphidxi_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	      }
		
	break;
      }


      
      //------------------------------------------------------------
      // 2D
    case 2:
      {
	// Compute the value of the approximation shape function i at quadrature point p
	if (calculate_phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      phi[i][p]      = FE<Dim,T>::shape       (elem, fe_type.order, i,    qp[p]);
	if (calculate_dphi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        dphidxi[i][p]  = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 0, qp[p]);
	        dphideta[i][p] = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 1, qp[p]);
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
                d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 0, qp[p]);
                d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 1, qp[p]);
                d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 2, qp[p]);
	      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	
	// Compute the value of the mapping shape function i at quadrature point p
	// (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
	    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
	        phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
	        dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
	        dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
                d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
                d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	        for (unsigned int p=1; p<n_qp; p++)
                  {
	            phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
	            dphidxi_map[i][p]  = dphidxi_map[i][0];
	            dphideta_map[i][p] = dphideta_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    d2phidxi2_map[i][p] = d2phidxi2_map[i][0];
                    d2phidxideta_map[i][p] = d2phidxideta_map[i][0];
                    d2phideta2_map[i][p] = d2phideta2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
	  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        phi_map[i][p]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
	        dphidxi_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
	        dphideta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
                d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
                d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	      }
			
       	break;
      }


      
      //------------------------------------------------------------
      // 3D
    case 3:
      {
	// Compute the value of the approximation shape function i at quadrature point p
	if (calculate_phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      phi[i][p]       = FE<Dim,T>::shape       (elem, fe_type.order, i,    qp[p]);
	if (calculate_dphi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        dphidxi[i][p]   = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 0, qp[p]);
	        dphideta[i][p]  = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 1, qp[p]);
	        dphidzeta[i][p] = FE<Dim,T>::shape_deriv (elem, fe_type.order, i, 2, qp[p]);
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<n_approx_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
                d2phidxi2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 0, qp[p]);
                d2phidxideta[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 1, qp[p]);
                d2phideta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 2, qp[p]);
                d2phidxidzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 3, qp[p]);
                d2phidetadzeta[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 4, qp[p]);
                d2phidzeta2[i][p] = FE<Dim,T>::shape_second_deriv (elem, fe_type.order, i, 5, qp[p]);
	      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	
	// Compute the value of the mapping shape function i at quadrature point p
	// (Lagrange shape functions are used for mapping)
        if (has_affine_map)
          {
	    for (unsigned int i=0; i<n_mapping_shape_functions; i++)
              {
	        phi_map[i][0]      = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[0]);
	        dphidxi_map[i][0]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
	        dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
	        dphidzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
                d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
                d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
                d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[0]);
                d2phidetadzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[0]);
                d2phidzeta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[0]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	        for (unsigned int p=1; p<n_qp; p++)
                  {
	            phi_map[i][p]      = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i,    qp[p]);
	            dphidxi_map[i][p]  = dphidxi_map[i][0];
	            dphideta_map[i][p] = dphideta_map[i][0];
	            dphidzeta_map[i][p] = dphidzeta_map[i][0];
#ifdef ENABLE_SECOND_DERIVATIVES
                    d2phidxi2_map[i][p] = d2phidxi2_map[i][0];
                    d2phidxideta_map[i][p] = d2phidxideta_map[i][0];
                    d2phideta2_map[i][p] = d2phideta2_map[i][0];
                    d2phidxideta_map[i][p] = d2phidxideta_map[i][0];
                    d2phidetadzeta_map[i][p] = d2phidetadzeta_map[i][0];
                    d2phidzeta2_map[i][p] = d2phidzeta2_map[i][0];
#endif // ifdef ENABLE_SECOND_DERIVATIVES
                  }
              }
          }
        else
	  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
	    for (unsigned int p=0; p<n_qp; p++)
	      {
	        phi_map[i][p]       = FE<Dim,LAGRANGE>::shape       (mapping_elem_type, mapping_order, i,    qp[p]);
	        dphidxi_map[i][p]   = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
	        dphideta_map[i][p]  = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
	        dphidzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
#ifdef ENABLE_SECOND_DERIVATIVES
                d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
                d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
                d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
                d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[p]);
                d2phidetadzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[p]);
                d2phidzeta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[p]);
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	      }
			
	break;
      }


    default:
      libmesh_error();
    }
  
  // Stop logging the shape function initialization
  STOP_LOG("init_shape_functions()", "FE");
}



    
#ifdef ENABLE_INFINITE_ELEMENTS

template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::init_base_shape_functions(const std::vector<Point>& qp,
					  const Elem* e)
{ 
  // I don't understand infinite elements well enough to risk
  // calculating too little.  :-(  RHS
  calculate_phi = calculate_dphi = calculate_d2phi = true;

  elem_type = e->type(); 
  init_shape_functions(qp, e); 
}

#endif // ENABLE_INFINITE_ELEMENTS
    


//--------------------------------------------------------------
// Explicit instantiations using macro from fe_macro.h

INSTANTIATE_FE(1);

INSTANTIATE_FE(2);

INSTANTIATE_FE(3);