fe_boundary.c

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

        // Calculate x at the point
	libmesh_assert (psi_map.size() == 1);
        // In the unlikely event we have multiple quadrature
        // points, they'll be in the same place
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    xyz[p].zero();
	    xyz[p].add_scaled          (side->point(0), psi_map[0][p]);
            normals[p] = normals[0];
	    JxW[p] = 1.0*qw[p];
          }

	// done computing the map
	break;
      }
      
    case 2:
      {
	// A 2D finite element living in either 2D or 3D space.
	// This means the boundary is a 1D finite element, i.e.
	// and EDGE2 or EDGE3.
	// Resize the vectors to hold data at the quadrature points
	{  
	  xyz.resize(n_qp);
	  dxyzdxi_map.resize(n_qp);
	  d2xyzdxi2_map.resize(n_qp);
	  tangents.resize(n_qp);
	  normals.resize(n_qp);
	  curvatures.resize(n_qp);
	  
	  JxW.resize(n_qp);
	}
	
	// Clear the entities that will be summed
	// Compute the tangent & normal at the quadrature point
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    tangents[p].resize(DIM-1); // 1 Tangent in 2D, 2 in 3D
	    xyz[p].zero();
	    dxyzdxi_map[p].zero();
	    d2xyzdxi2_map[p].zero();
	  }
	
	// compute x, dxdxi at the quadrature points    
	for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
	  {
	    const Point& side_point = side->point(i);
	    
	    for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
	      {	  
		xyz[p].add_scaled          (side_point, psi_map[i][p]);
		dxyzdxi_map[p].add_scaled  (side_point, dpsidxi_map[i][p]);
		d2xyzdxi2_map[p].add_scaled(side_point, d2psidxi2_map[i][p]);
	      }
	  }

	// Compute the tangent & normal at the quadrature point
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    // The first tangent comes from just the edge's Jacobian
	    tangents[p][0] = dxyzdxi_map[p].unit();
	    
#if DIM == 2
	    // For a 2D element living in 2D, the normal is given directly
	    // from the entries in the edge Jacobian.
	    normals[p] = (Point(dxyzdxi_map[p](1), -dxyzdxi_map[p](0), 0.)).unit();
	    
#elif DIM == 3
	    // For a 2D element living in 3D, there is a second tangent.
	    // For the second tangent, we need to refer to the full
	    // element's (not just the edge's) Jacobian.
	    const Elem *elem = side->parent();
	    libmesh_assert (elem != NULL);

	    // Inverse map xyz[p] to a reference point on the parent...
	    Point reference_point = FE<2,LAGRANGE>::inverse_map(elem, xyz[p]);
	    
	    // Get dxyz/dxi and dxyz/deta from the parent map.
	    Point dx_dxi  = FE<2,LAGRANGE>::map_xi (elem, reference_point);
	    Point dx_deta = FE<2,LAGRANGE>::map_eta(elem, reference_point);

	    // The second tangent vector is formed by crossing these vectors.
	    tangents[p][1] = dx_dxi.cross(dx_deta).unit();

	    // Finally, the normal in this case is given by crossing these
	    // two tangents.
	    normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
#endif 
	    

	    // The curvature is computed via the familiar Frenet formula:
	    // curvature = [d^2(x) / d (xi)^2] dot [normal]
	    // For a reference, see:
	    // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
	    //
	    // Note: The sign convention here is different from the
	    // 3D case.  Concave-upward curves (smiles) have a positive
	    // curvature.  Concave-downward curves (frowns) have a
	    // negative curvature.  Be sure to take that into account!
	    const Real numerator   = d2xyzdxi2_map[p] * normals[p];
	    const Real denominator = dxyzdxi_map[p].size_sq();
	    libmesh_assert (denominator != 0);
	    curvatures[p] = numerator / denominator;
	  }
	
	// compute the jacobian at the quadrature points
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    const Real jac = dxyzdxi_map[p].size();
	    
	    libmesh_assert (jac > 0.);
	    
	    JxW[p] = jac*qw[p];
	  }
	
	// done computing the map
	break;
      }


      
    case 3:
      {
	// A 3D finite element living in 3D space.
	// Resize the vectors to hold data at the quadrature points
	{  
	  xyz.resize(n_qp);
	  dxyzdxi_map.resize(n_qp);
	  dxyzdeta_map.resize(n_qp);
	  d2xyzdxi2_map.resize(n_qp);
	  d2xyzdxideta_map.resize(n_qp);
	  d2xyzdeta2_map.resize(n_qp);
	  tangents.resize(n_qp);
	  normals.resize(n_qp);
	  curvatures.resize(n_qp);

	  JxW.resize(n_qp);
	}
    
	// Clear the entities that will be summed
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    tangents[p].resize(DIM-1); // 1 Tangent in 2D, 2 in 3D
	    xyz[p].zero();
	    dxyzdxi_map[p].zero();
	    dxyzdeta_map[p].zero();
	    d2xyzdxi2_map[p].zero();
	    d2xyzdxideta_map[p].zero();
	    d2xyzdeta2_map[p].zero();
	  }
	
	// compute x, dxdxi at the quadrature points    
	for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
	  {
	    const Point& side_point = side->point(i);
	    
	    for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
	      {
		xyz[p].add_scaled         (side_point, psi_map[i][p]);
		dxyzdxi_map[p].add_scaled (side_point, dpsidxi_map[i][p]);
		dxyzdeta_map[p].add_scaled(side_point, dpsideta_map[i][p]);
		d2xyzdxi2_map[p].add_scaled   (side_point, d2psidxi2_map[i][p]);
		d2xyzdxideta_map[p].add_scaled(side_point, d2psidxideta_map[i][p]);
		d2xyzdeta2_map[p].add_scaled  (side_point, d2psideta2_map[i][p]);
	      }
	  }

	// Compute the tangents, normal, and curvature at the quadrature point
	for (unsigned int p=0; p<n_qp; p++)
	  {	    
	    const Point n  = dxyzdxi_map[p].cross(dxyzdeta_map[p]);
	    normals[p]     = n.unit();
	    tangents[p][0] = dxyzdxi_map[p].unit();
	    tangents[p][1] = n.cross(dxyzdxi_map[p]).unit();
	    
	    // Compute curvature using the typical nomenclature
	    // of the first and second fundamental forms.
	    // For reference, see:
	    // 1) http://mathworld.wolfram.com/MeanCurvature.html
	    //    (note -- they are using inward normal)
	    // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
	    const Real L  = -d2xyzdxi2_map[p]    * normals[p];
	    const Real M  = -d2xyzdxideta_map[p] * normals[p];
	    const Real N  = -d2xyzdeta2_map[p]   * normals[p];
	    const Real E  =  dxyzdxi_map[p].size_sq();
	    const Real F  =  dxyzdxi_map[p]      * dxyzdeta_map[p];
	    const Real G  =  dxyzdeta_map[p].size_sq();
	    
	    const Real numerator   = E*N -2.*F*M + G*L;
	    const Real denominator = E*G - F*F;
	    libmesh_assert (denominator != 0.);
	    curvatures[p] = 0.5*numerator/denominator;
	  }  
    	
	// compute the jacobian at the quadrature points, see
	// http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
	for (unsigned int p=0; p<n_qp; p++)
	  {
	    const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
			      dydxi_map(p)*dydxi_map(p) +
			      dzdxi_map(p)*dzdxi_map(p));
	    
	    const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
			      dydxi_map(p)*dydeta_map(p) +
			      dzdxi_map(p)*dzdeta_map(p));
	    
	    const Real g21 = g12;
	    
	    const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
			      dydeta_map(p)*dydeta_map(p) +
			      dzdeta_map(p)*dzdeta_map(p));
	    
	    
	    const Real jac = std::sqrt(g11*g22 - g12*g21);
	    
	    libmesh_assert (jac > 0.);

	    JxW[p] = jac*qw[p];
	  }
	
	// done computing the map
	break;
      }


    default:
      libmesh_error();
      
    }
  STOP_LOG("compute_face_map()", "FE");
}




void FEBase::compute_edge_map(const std::vector<Real>& qw,
			      const Elem* edge)
{
  libmesh_assert (edge != NULL);

  if (dim == 2)
    {
      // A 2D finite element living in either 2D or 3D space.
      // The edges here are the sides of the element, so the
      // (misnamed) compute_face_map function does what we want
      FEBase::compute_face_map(qw, edge);
      return;
    }

  libmesh_assert (dim == 3);  // 1D is unnecessary and currently unsupported

  START_LOG("compute_edge_map()", "FE");

  // The number of quadrature points.
  const unsigned int n_qp = qw.size();
  
  // Resize the vectors to hold data at the quadrature points
  xyz.resize(n_qp);
  dxyzdxi_map.resize(n_qp);
  dxyzdeta_map.resize(n_qp);
  d2xyzdxi2_map.resize(n_qp);
  d2xyzdxideta_map.resize(n_qp);
  d2xyzdeta2_map.resize(n_qp);
  tangents.resize(n_qp);
  normals.resize(n_qp);
  curvatures.resize(n_qp);

  JxW.resize(n_qp);
    
  // Clear the entities that will be summed
  for (unsigned int p=0; p<n_qp; p++)
    {
      tangents[p].resize(1);
      xyz[p].zero();
      dxyzdxi_map[p].zero();
      dxyzdeta_map[p].zero();
      d2xyzdxi2_map[p].zero();
      d2xyzdxideta_map[p].zero();
      d2xyzdeta2_map[p].zero();
    }

  // compute x, dxdxi at the quadrature points    
  for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
    {
      const Point& edge_point = edge->point(i);
      
      for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
        {
	  xyz[p].add_scaled             (edge_point, psi_map[i][p]);
	  dxyzdxi_map[p].add_scaled     (edge_point, dpsidxi_map[i][p]);
	  d2xyzdxi2_map[p].add_scaled   (edge_point, d2psidxi2_map[i][p]);
        }
    }

  // Compute the tangents at the quadrature point
  // FIXME: normals (plural!) and curvatures are uncalculated
  for (unsigned int p=0; p<n_qp; p++)
    {    
      const Point n  = dxyzdxi_map[p].cross(dxyzdeta_map[p]);
      tangents[p][0] = dxyzdxi_map[p].unit();

      // compute the jacobian at the quadrature points
      const Real jac = std::sqrt(dxdxi_map(p)*dxdxi_map(p) +
				 dydxi_map(p)*dydxi_map(p) +
				 dzdxi_map(p)*dzdxi_map(p));
	    
      libmesh_assert (jac > 0.);

      JxW[p] = jac*qw[p];
    }

  STOP_LOG("compute_edge_map()", "FE");
}




//--------------------------------------------------------------
// Explicit instantiations
template void FE<1,LAGRANGE>::reinit(Elem const*, unsigned int, Real);
template void FE<1,HIERARCHIC>::reinit(Elem const*, unsigned int, Real);
template void FE<1,CLOUGH>::reinit(Elem const*, unsigned int, Real);
template void FE<1,HERMITE>::reinit(Elem const*, unsigned int, Real);
template void FE<1,MONOMIAL>::reinit(Elem const*, unsigned int, Real);
#ifdef ENABLE_HIGHER_ORDER_SHAPES
template void FE<1,BERNSTEIN>::reinit(Elem const*, unsigned int, Real);
template void FE<1,SZABAB>::reinit(Elem const*, unsigned int, Real);
#endif
template void FE<1,XYZ>::reinit(Elem const*, unsigned int, Real);

template void FE<2,LAGRANGE>::reinit(Elem const*, unsigned int, Real);
template void FE<2,LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<2,HIERARCHIC>::reinit(Elem const*, unsigned int, Real);
template void FE<2,HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<2,CLOUGH>::reinit(Elem const*, unsigned int, Real);
template void FE<2,CLOUGH>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<2,HERMITE>::reinit(Elem const*, unsigned int, Real);
template void FE<2,HERMITE>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<2,MONOMIAL>::reinit(Elem const*, unsigned int, Real);
template void FE<2,MONOMIAL>::edge_reinit(Elem const*, unsigned int, Real);
#ifdef ENABLE_HIGHER_ORDER_SHAPES
template void FE<2,BERNSTEIN>::reinit(Elem const*, unsigned int, Real);
template void FE<2,BERNSTEIN>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<2,SZABAB>::reinit(Elem const*, unsigned int, Real);
template void FE<2,SZABAB>::edge_reinit(Elem const*, unsigned int, Real);
#endif
template void FE<2,XYZ>::reinit(Elem const*, unsigned int, Real);
template void FE<2,XYZ>::edge_reinit(Elem const*, unsigned int, Real);

// Intel 9.1 complained it needed this in devel mode.
template void FE<2,XYZ>::init_face_shape_functions(const std::vector<Point>&, const Elem*);

template void FE<3,LAGRANGE>::reinit(Elem const*, unsigned int, Real);
template void FE<3,LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<3,HIERARCHIC>::reinit(Elem const*, unsigned int, Real);
template void FE<3,HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<3,CLOUGH>::reinit(Elem const*, unsigned int, Real);
template void FE<3,CLOUGH>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<3,HERMITE>::reinit(Elem const*, unsigned int, Real);
template void FE<3,HERMITE>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<3,MONOMIAL>::reinit(Elem const*, unsigned int, Real);
template void FE<3,MONOMIAL>::edge_reinit(Elem const*, unsigned int, Real);
#ifdef ENABLE_HIGHER_ORDER_SHAPES
template void FE<3,BERNSTEIN>::reinit(Elem const*, unsigned int, Real);
template void FE<3,BERNSTEIN>::edge_reinit(Elem const*, unsigned int, Real);
template void FE<3,SZABAB>::reinit(Elem const*, unsigned int, Real);
template void FE<3,SZABAB>::edge_reinit(Elem const*, unsigned int, Real);
#endif
template void FE<3,XYZ>::reinit(Elem const*, unsigned int, Real);
template void FE<3,XYZ>::edge_reinit(Elem const*, unsigned int, Real);

// Intel 9.1 complained it needed this in devel mode.
template void FE<3,XYZ>::init_face_shape_functions(const std::vector<Point>&, const Elem*);