fe_base.c

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

		    std::cerr << "ERROR: Don't build an infinite element " << std::endl
			      << " with InfMapType = " << fet.inf_map << std::endl;
		    libmesh_error();
		}
	    }


	    
	  default:
	    std::cerr << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
	    libmesh_error();
	  }

      }

      


      // 3D
    case 3:
      {
	switch (fet.radial_family)
	  {
	  case INFINITE_MAP:
	    {
	      std::cerr << "ERROR: Don't build an infinite element " << std::endl
			<< " with FEFamily = " << fet.radial_family << std::endl;
	      libmesh_error();
	    }

	  case JACOBI_20_00:
	    {
  	      switch (fet.inf_map)
	        {
		  case CARTESIAN:
		    {
		      AutoPtr<FEBase> ap(new InfFE<3,JACOBI_20_00,CARTESIAN>(fet));
		      return ap;
		    }
		  default:
		    std::cerr << "ERROR: Don't build an infinite element " << std::endl
			      << " with InfMapType = " << fet.inf_map << std::endl;
		    libmesh_error();
		}
	    }

	  case JACOBI_30_00:
	    {
  	      switch (fet.inf_map)
	        {
		  case CARTESIAN:
		    {
		      AutoPtr<FEBase> ap(new InfFE<3,JACOBI_30_00,CARTESIAN>(fet));
		      return ap;
		    }
		  default:
		    std::cerr << "ERROR: Don't build an infinite element " << std::endl
			      << " with InfMapType = " << fet.inf_map << std::endl;
		    libmesh_error();
		}
	    }

	  case LEGENDRE:
	    {
  	      switch (fet.inf_map)
	        {
		  case CARTESIAN:
		    {
		      AutoPtr<FEBase> ap(new InfFE<3,LEGENDRE,CARTESIAN>(fet));
		      return ap;
		    }
		  default:
		    std::cerr << "ERROR: Don't build an infinite element " << std::endl
			      << " with InfMapType = " << fet.inf_map << std::endl;
		    libmesh_error();
		}
	    }

	  case LAGRANGE:
	    {
  	      switch (fet.inf_map)
	        {
		  case CARTESIAN:
		    {
		      AutoPtr<FEBase> ap(new InfFE<3,LAGRANGE,CARTESIAN>(fet));
		      return ap;
		    }
		  default:
		    std::cerr << "ERROR: Don't build an infinite element " << std::endl
			      << " with InfMapType = " << fet.inf_map << std::endl;
		    libmesh_error();
		}
	    }


	    
	  default:
	    std::cerr << "ERROR: Bad FEType.radial_family= " << fet.radial_family << std::endl;
	    libmesh_error();
	  }
      }

    default:
      libmesh_error();
    }

  libmesh_error();
  AutoPtr<FEBase> ap(NULL);
  return ap;
}




#endif // ifdef ENABLE_INFINITE_ELEMENTS











void FEBase::compute_shape_functions (const Elem*)
{
  //-------------------------------------------------------------------------
  // Compute the shape function values (and derivatives)
  // at the Quadrature points.  Note that the actual values
  // have already been computed via init_shape_functions

  // Start logging the shape function computation
  START_LOG("compute_shape_functions()", "FE");

  calculations_started = true;

  // If the user forgot to request anything, we'll be safe and
  // calculate everything:
#ifdef ENABLE_SECOND_DERIVATIVES
  if (!calculate_phi && !calculate_dphi && !calculate_d2phi)
    calculate_phi = calculate_dphi = calculate_d2phi = true;
#else
  if (!calculate_phi && !calculate_dphi)
    calculate_phi = calculate_dphi = true;
#endif // ENABLE_SECOND_DERIVATIVES

  // Compute the value of the derivative shape function i at quadrature point p
  switch (dim)
    {
      
    case 1:
      {
	if (calculate_dphi)
	  for (unsigned int i=0; i<dphi.size(); i++)
	    for (unsigned int p=0; p<dphi[i].size(); p++)
	      {
	        // dphi/dx    = (dphi/dxi)*(dxi/dx)
	        dphi[i][p](0) =
		  dphidx[i][p] = dphidxi[i][p]*dxidx_map[p];
	      
	        dphi[i][p](1) = dphidy[i][p] = 0.;
	        dphi[i][p](2) = dphidz[i][p] = 0.;
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<d2phi.size(); i++)
	    for (unsigned int p=0; p<d2phi[i].size(); p++)
	      {
	        d2phi[i][p](0,0) = d2phidx2[i][p] = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p];
#if DIM>1
	        d2phi[i][p](0,1) = d2phidxdy[i][p] = 
		  d2phi[i][p](1,0) = 0.;
	        d2phi[i][p](1,1) = d2phidy2[i][p] = 0.;
#if DIM>2
	        d2phi[i][p](0,2) = d2phidxdz[i][p] =
		  d2phi[i][p](2,0) = 0.;
	        d2phi[i][p](1,2) = d2phidydz[i][p] = 
		  d2phi[i][p](2,1) = 0.;
	        d2phi[i][p](2,2) = d2phidz2[i][p] = 0.;
#endif
#endif
	      }
#endif

	// All done
	break;
      }

    case 2:
      {
	if (calculate_dphi)
	  for (unsigned int i=0; i<dphi.size(); i++)
	    for (unsigned int p=0; p<dphi[i].size(); p++)
	      {
	        // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx)
	        dphi[i][p](0) =
		  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
				  dphideta[i][p]*detadx_map[p]);
	      
	        // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy)
	        dphi[i][p](1) =
		  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
				  dphideta[i][p]*detady_map[p]);
	      
	        // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz)
#if DIM == 3  
	        dphi[i][p](2) = // can only assign to the Z component if DIM==3
#endif
		dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
				dphideta[i][p]*detadz_map[p]);
	      }

#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<d2phi.size(); i++)
	    for (unsigned int p=0; p<d2phi[i].size(); p++)
	      {
	        d2phi[i][p](0,0) = d2phidx2[i][p] = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
		  2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detadx_map[p];
	        d2phi[i][p](0,1) = d2phidxdy[i][p] =
		  d2phi[i][p](1,0) = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +
		  d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
		  d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p];
	        d2phi[i][p](1,1) = d2phidy2[i][p] =
		  d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
		  2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
		  d2phideta2[i][p]*detady_map[p]*detady_map[p];
#if DIM == 3  
	        d2phi[i][p](0,2) = d2phidxdz[i][p] = 
		  d2phi[i][p](2,0) = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +
		  d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
		  d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p];
	        d2phi[i][p](1,2) = d2phidydz[i][p] = 
		  d2phi[i][p](2,1) =
		  d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +
		  d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] +
		  d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
		  d2phidxideta[i][p]*detady_map[p]*dxidz_map[p];
	        d2phi[i][p](2,2) = d2phidz2[i][p] =
		  d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
		  2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
		  d2phideta2[i][p]*detadz_map[p]*detadz_map[p];
#endif
	      }
#endif

	// All done
	break;
      }
    
    case 3:
      {
	if (calculate_dphi)
	  for (unsigned int i=0; i<dphi.size(); i++)
	    for (unsigned int p=0; p<dphi[i].size(); p++)
	      {
	        // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
	        dphi[i][p](0) =
		  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
				  dphideta[i][p]*detadx_map[p] +
				  dphidzeta[i][p]*dzetadx_map[p]);
		
	        // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
	        dphi[i][p](1) =
		  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
				  dphideta[i][p]*detady_map[p] +
				  dphidzeta[i][p]*dzetady_map[p]);
		
	        // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
	        dphi[i][p](2) =
		  dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
				  dphideta[i][p]*detadz_map[p] +
				  dphidzeta[i][p]*dzetadz_map[p]);	      
	      }

#ifdef ENABLE_SECOND_DERIVATIVES
	if (calculate_d2phi)
	  for (unsigned int i=0; i<d2phi.size(); i++)
	    for (unsigned int p=0; p<d2phi[i].size(); p++)
	      {
	        d2phi[i][p](0,0) = d2phidx2[i][p] = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidx_map[p] +
		  2*d2phidxideta[i][p]*dxidx_map[p]*detadx_map[p] +
		  2*d2phidxidzeta[i][p]*dxidx_map[p]*dzetadx_map[p] +
		  2*d2phidetadzeta[i][p]*detadx_map[p]*dzetadx_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detadx_map[p] +
		  d2phidzeta2[i][p]*dzetadx_map[p]*dzetadx_map[p];
	        d2phi[i][p](0,1) = d2phidxdy[i][p] =
		  d2phi[i][p](1,0) = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidy_map[p] +
		  d2phidxideta[i][p]*dxidx_map[p]*detady_map[p] +
		  d2phidxidzeta[i][p]*dxidx_map[p]*dzetady_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detady_map[p] +
		  d2phidxideta[i][p]*detadx_map[p]*dxidy_map[p] +
		  d2phidetadzeta[i][p]*detadx_map[p]*dzetady_map[p] +
		  d2phidzeta2[i][p]*dzetadx_map[p]*dzetady_map[p] +
		  d2phidxidzeta[i][p]*dzetadx_map[p]*dxidy_map[p] +
		  d2phidetadzeta[i][p]*dzetadx_map[p]*detady_map[p];
	        d2phi[i][p](0,2) = d2phidxdz[i][p] = 
		  d2phi[i][p](2,0) = 
		  d2phidxi2[i][p]*dxidx_map[p]*dxidz_map[p] +
		  d2phidxideta[i][p]*dxidx_map[p]*detadz_map[p] +
		  d2phidxidzeta[i][p]*dxidx_map[p]*dzetadz_map[p] +
		  d2phideta2[i][p]*detadx_map[p]*detadz_map[p] +
		  d2phidxideta[i][p]*detadx_map[p]*dxidz_map[p] +
		  d2phidetadzeta[i][p]*detadx_map[p]*dzetadz_map[p] +
		  d2phidzeta2[i][p]*dzetadx_map[p]*dzetadz_map[p] +
		  d2phidxidzeta[i][p]*dzetadx_map[p]*dxidz_map[p] +
		  d2phidetadzeta[i][p]*dzetadx_map[p]*detadz_map[p];
	        d2phi[i][p](1,1) = d2phidy2[i][p] =
		  d2phidxi2[i][p]*dxidy_map[p]*dxidy_map[p] +
		  2*d2phidxideta[i][p]*dxidy_map[p]*detady_map[p] +
		  2*d2phidxidzeta[i][p]*dxidy_map[p]*dzetady_map[p] +
		  2*d2phidetadzeta[i][p]*detady_map[p]*dzetady_map[p] +
		  d2phideta2[i][p]*detady_map[p]*detady_map[p] +
		  d2phidzeta2[i][p]*dzetady_map[p]*dzetady_map[p];
	        d2phi[i][p](1,2) = d2phidydz[i][p] = 
		  d2phi[i][p](2,1) =
		  d2phidxi2[i][p]*dxidy_map[p]*dxidz_map[p] +
		  d2phidxideta[i][p]*dxidy_map[p]*detadz_map[p] +
		  d2phidxidzeta[i][p]*dxidy_map[p]*dzetadz_map[p] +
		  d2phideta2[i][p]*detady_map[p]*detadz_map[p] +
		  d2phidxideta[i][p]*detady_map[p]*dxidz_map[p] +
		  d2phidetadzeta[i][p]*detady_map[p]*dzetadz_map[p] +
		  d2phidzeta2[i][p]*dzetady_map[p]*dzetadz_map[p] +
		  d2phidxidzeta[i][p]*dzetady_map[p]*dxidz_map[p] +
		  d2phidetadzeta[i][p]*dzetady_map[p]*detadz_map[p];
	        d2phi[i][p](2,2) = d2phidz2[i][p] =
		  d2phidxi2[i][p]*dxidz_map[p]*dxidz_map[p] +
		  2*d2phidxideta[i][p]*dxidz_map[p]*detadz_map[p] +
		  2*d2phidxidzeta[i][p]*dxidz_map[p]*dzetadz_map[p] +
		  2*d2phidetadzeta[i][p]*detadz_map[p]*dzetadz_map[p] +
		  d2phideta2[i][p]*detadz_map[p]*detadz_map[p] +
		  d2phidzeta2[i][p]*dzetadz_map[p]*dzetadz_map[p];
	      }
#endif
	// All done
	break;
      }

    default:
      {
	libmesh_error();
      }
    }
  
  // Stop logging the shape function computation
  STOP_LOG("compute_shape_functions()", "FE");
}



bool FEBase::on_reference_element(const Point& p, const ElemType t, const Real eps)
{
  libmesh_assert (eps >= 0.);
  
  const Real xi   = p(0);
  const Real eta  = p(1);
  const Real zeta = p(2);
  
  switch (t)
    {

    case EDGE2:
    case EDGE3:
    case EDGE4:
      {
	// The reference 1D element is [-1,1].
	if ((xi >= -1.-eps) &&
	    (xi <=  1.+eps))
	  return true;

	return false;
      }

      
    case TRI3:
    case TRI6:
      {
	// The reference triangle is isocoles
	// and is bound by xi=0, eta=0, and xi+eta=1.
	if ((xi  >= 0.-eps) &&
	    (eta >= 0.-eps) &&
	    ((xi + eta) <= 1.+eps))
	  return true;

	return false;
      }

      
    case QUAD4:
    case QUAD8:
    case QUAD9:
      {
	// The reference quadrilateral element is [-1,1]^2.
	if ((xi  >= -1.-eps) &&
	    (xi  <=  1.+eps) &&
	    (eta >= -1.-eps) &&
	    (eta <=  1.+eps))
	  return true;
		
	return false;
      }


    case TET4:
    case TET10:
      {
	// The reference tetrahedral is isocoles
	// and is bound by xi=0, eta=0, zeta=0,
	// and xi+eta+zeta=1.
	if ((xi   >= 0.-eps) &&