fem2d_poisson.cpp

c++ for poisson

//****************************************************************************80*
//
//  Purpose:
//
//    NODES_WRITE writes the nodes to a file.
//
//  Modified:
//
//    22 March 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NODE_NUM, the number of nodes.
//
//    Input, double NODE_XY[2*NODE_NUM], the X and Y coordinates of nodes.
//
//    Input, char *OUTPUT_FILENAME, the name of the file
//    in which the data should be stored.
//
{
  int node;
  ofstream output;
  double x;
  double y;

  output.open ( output_filename );

  if ( !output )
  {
    cout << "\n";
    cout << "NODES_WRITE - Warning!\n";
    cout << "  Could not write the node file.\n";
    return;
  }

  for ( node = 0; node < node_num; node++ )
  {
    x = node_xy[0+node*2];
    y = node_xy[1+node*2];

    output << setw(8)  << x << "  "
           << setw(8)  << y << "\n";
  }

  output.close ( );

  return;
}
//****************************************************************************80

void qbf ( double x, double y, int element, int inode, double node_xy[], 
  int element_node[], int element_num, int nnodes, 
  int node_num, double *b, double *dbdx, double *dbdy )

//****************************************************************************80
//
//  Purpose:
//
//    QBF evaluates the quadratic basis functions.
//
//  Discussion:
//
//    This routine assumes that the "midpoint" nodes are, in fact,
//    exactly the average of the two extreme nodes.  This is NOT true
//    for a general quadratic triangular element.
//
//    Assuming this property of the midpoint nodes makes it easy to
//    determine the values of (R,S) in the reference element that
//    correspond to (X,Y) in the physical element.
//
//    Once we know the (R,S) coordinates, it's easy to evaluate the
//    basis functions and derivatives.
//
//  The physical element T6:
//
//    In this picture, we don't mean to suggest that the bottom of
//    the physical triangle is horizontal.  However, we do assume that
//    each of the sides is a straight line, and that the intermediate
//    points are exactly halfway on each side.
//
//    |
//    |
//    |        3
//    |       / \
//    |      /   \
//    Y     6     5
//    |    /       \
//    |   /         \
//    |  1-----4-----2
//    |
//    +--------X-------->
//
//  Reference element T6:
//
//    In this picture of the reference element, we really do assume
//    that one side is vertical, one horizontal, of length 1.
//
//    |
//    |
//    1  3
//    |  |\
//    |  | \
//    S  6  5
//    |  |   \
//    |  |    \
//    0  1--4--2
//    |
//    +--0--R--1-------->
//
//  Modified:
//
//    18 June 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the (global) coordinates of the point
//    at which the basis function is to be evaluated.
//
//    Input, int ELEMENT, the index of the element which contains the point.
//
//    Input, int INODE, the local index (between 1 and 6) that
//    specifies which basis function is to be evaluated.
//
//    Input, double NODE_XY[2*NODE_NUM], the nodes.
//
//    Input, int ELEMENT_NODE[NNODES*ELEMENT_NUM];
//    ELEMENT_NODE(I,J) is the global index of local node I in element J.
//
//    Input, int ELEMENT_NUM, the number of elements.
//
//    Input, int NNODES, the number of nodes used to form one element.
//
//    Input, int NODE_NUM, the number of nodes.
//
//    Output, double *B, *DBDX, *DBDY, the value of the basis function
//    and its X and Y derivatives at (X,Y).
//
{
  double dbdr;
  double dbds;
  double det;
  double drdx;
  double drdy;
  double dsdx;
  double dsdy;
  int i;
  double r;
  double s;
  double xn[6];
  double yn[6];

  for ( i = 0; i < 6; i++ )
  {
    xn[i] = node_xy[0+(element_node[i+(element-1)*nnodes]-1)*2];
    yn[i] = node_xy[1+(element_node[i+(element-1)*nnodes]-1)*2];
  }
//
//  Determine the (R,S) coordinates corresponding to (X,Y).
//
//  What is happening here is that we are solving the linear system:
//
//    ( X2-X1  X3-X1 ) * ( R ) = ( X - X1 )
//    ( Y2-Y1  Y3-Y1 )   ( S )   ( Y - Y1 )
//
//  by computing the inverse of the coefficient matrix and multiplying
//  it by the right hand side to get R and S.
//
//  The values of dRdX, dRdY, dSdX and dSdY are easily from the formulas
//  for R and S.
//
  det =   ( xn[1] - xn[0] ) * ( yn[2] - yn[0] ) 
        - ( xn[2] - xn[0] ) * ( yn[1] - yn[0] );

  r = ( ( yn[2] - yn[0] ) * ( x     - xn[0] ) 
      + ( xn[0] - xn[2] ) * ( y     - yn[0] ) ) / det;

  drdx = ( yn[2] - yn[0] ) / det;
  drdy = ( xn[0] - xn[2] ) / det;

  s = ( ( yn[0] - yn[1] ) * ( x     - xn[0] ) 
      + ( xn[1] - xn[0] ) * ( y     - yn[0] ) ) / det;

  dsdx = ( yn[0] - yn[1] ) / det;
  dsdy = ( xn[1] - xn[0] ) / det;
//
//  The basis functions can now be evaluated in terms of the
//  reference coordinates R and S.  It's also easy to determine
//  the values of the derivatives with respect to R and S.
//
  if ( inode == 1 )
  {
    *b   =   2.0E+00 *     ( 1.0E+00 - r - s ) * ( 0.5E+00 - r - s );
    dbdr = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s;
    dbds = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s;
  }
  else if ( inode == 2 )
  {
    *b   =   2.0E+00 * r * ( r - 0.5E+00 );
    dbdr = - 1.0E+00 + 4.0E+00 * r;
    dbds =   0.0E+00;
  }
  else if ( inode == 3 )
  {
    *b   =   2.0E+00 * s * ( s - 0.5E+00 );
    dbdr =   0.0E+00;
    dbds = - 1.0E+00               + 4.0E+00 * s;
  }
  else if ( inode == 4 )
  {
    *b   =   4.0E+00 * r * ( 1.0E+00 - r - s );
    dbdr =   4.0E+00 - 8.0E+00 * r - 4.0E+00 * s;
    dbds =           - 4.0E+00 * r;
  }
  else if ( inode == 5 )
  {
    *b   =   4.0E+00 * r * s;
    dbdr =                           4.0E+00 * s;
    dbds =             4.0E+00 * r;
  }
  else if ( inode == 6 )
  {
    *b   =   4.0E+00 * s * ( 1.0E+00 - r - s );
    dbdr =                         - 4.0E+00 * s;
    dbds =   4.0E+00 - 4.0E+00 * r - 8.0E+00 * s;
  }
  else
  {
    cout << "\n";
    cout << "QBF - Fatal error!\n";
    cout << "  Request for local basis function INODE = " << inode << "\n";
    exit ( 1 );
  }
//
//  We need to convert the derivative information from (R(X,Y),S(X,Y))
//  to (X,Y) using the chain rule.
//
  *dbdx = dbdr * drdx + dbds * dsdx;
  *dbdy = dbdr * drdy + dbds * dsdy;

  return;
}
//****************************************************************************80

void quad_a ( double node_xy[], int element_node[], 
  int element_num, int node_num, int nnodes, double wq[], double xq[], 
  double yq[] )

//****************************************************************************80
//
//  Purpose:
//
//    QUAD_A sets the quadrature rule for assembly.
//
//  Modified:
//
//    07 April 2004
//
//  Parameters:
//
//    Input, double NODE_XY[2*NODE_NUM], the nodes.
//
//    Input, int ELEMENT_NODE[NNODES*ELEMENT_NUM]; 
//    ELEMENT_NODE(I,J) is the global index of local node I in element J.
//
//    Input, int ELEMENT_NUM, the number of elements.
//
//    Input, int NODE_NUM, the number of nodes.
//
//    Input, int NNODES, the number of nodes used to form one element.
//
//    Output, double WQ[3], quadrature weights.
//
//    Output, double XQ[3*ELEMENT_NUM], YQ[3*ELEMENT_NUM], the  
//    coordinates of the quadrature points in each element.
//
{
  int element;
  int ip1;
  int ip2;
  int ip3;
  double x1;
  double x2;
  double x3;
  double y1;
  double y2;
  double y3;

  wq[0] = 1.0E+00 / 3.0E+00;
  wq[1] = wq[0];
  wq[2] = wq[0];

  for ( element = 1; element <= element_num; element++ )
  {
    ip1 = element_node[0+(element-1)*nnodes];
    ip2 = element_node[1+(element-1)*nnodes];
    ip3 = element_node[2+(element-1)*nnodes];

    x1 = node_xy[0+(ip1-1)*2];
    x2 = node_xy[0+(ip2-1)*2];
    x3 = node_xy[0+(ip3-1)*2];

    y1 = node_xy[1+(ip1-1)*2];
    y2 = node_xy[1+(ip2-1)*2];
    y3 = node_xy[1+(ip3-1)*2];

    xq[0+(element-1)*3] = 0.5E+00 * ( x1 + x2 );
    xq[1+(element-1)*3] = 0.5E+00 * ( x2 + x3 );
    xq[2+(element-1)*3] = 0.5E+00 * ( x1 + x3 );

    yq[0+(element-1)*3] = 0.5E+00 * ( y1 + y2 );
    yq[1+(element-1)*3] = 0.5E+00 * ( y2 + y3 );
    yq[2+(element-1)*3] = 0.5E+00 * ( y1 + y3 );
  }

  return;
}
//****************************************************************************80*

void quad_e ( double node_xy[], int element_node[], 
  int element, int element_num, int nnodes, int node_num, int nqe, 
  double wqe[], double xqe[], double yqe[] )

//****************************************************************************80*
//
//  Purpose:
//
//    QUAD_E sets a quadrature rule for the error calculation.
//
//  Modified:
//
//    07 April 2004
//
//  Parameters:
//
//    Input, double NODE_XY[2*NODE_NUM], the X and Y coordinates of nodes.
//
//    Input, int ELEMENT_NODE[NNODES*ELEMENT_NUM]; ELEMENT_NODE(I,J) is the global 
//    index of local node I in element J.
//
//    Input, int ELEMENT, the index of the element for which the quadrature
//    points are to be computed.
//
//    Input, int ELEMENT_NUM, the number of elements.
//
//    Input, int NNODES, the number of nodes used to form one element.
//
//    Input, int NODE_NUM, the number of nodes.
//
//    Input, int NQE, the number of points in the quadrature rule.
//    This is actually fixed at 13.
//
//    Output, double WQE[NQE], the quadrature weights.
//
//    Output, double XQE[NQE], YQE[NQE], the X and Y coordinates
//    of the quadrature points.
//
{
  int i;
  int ii;
  int iii;
  int ip1;
  int ip2;
  int ip3;
  double x1;
  double x2;
  double x3;
  double y1;
  double y2;
  double y3;
  double z1;
  double z2;
  double z3;
  double z4;
  double z5;
  double z6;
  double z7;
//
  for ( i = 1; i <= 3; i++ )
  {
    wqe[i-1] = 0.175615257433204E+00;
    ii = i + 3;
    wqe[ii-1] = 0.053347235608839E+00;
    ii = i + 6;
    iii = ii + 3;
    wqe[ii-1] = 0.077113760890257E+00;
    wqe[iii-1] = wqe[ii-1];
  }

  wqe[13-1] = -0.14957004446767E+00;

  z1 = 0.479308067841923E+00;
  z2 = 0.260345966079038E+00;
  z3 = 0.869739794195568E+00;
  z4 = 0.065130102902216E+00;
  z5 = 0.638444188569809E+00;
  z6 = 0.312865496004875E+00;
  z7 = 0.048690315425316E+00;
 
  ip1 = element_node[0+(element-1)*nnodes];
  ip2 = element_node[1+(element-1)*nnodes];
  ip3 = element_node[2+(element-1)*nnodes];

  x1 = node_xy[0+(ip1-1)*2];
  x2 = node_xy[0+(ip2-1)*2];
  x3 = node_xy[0+(ip3-1)*2];

  y1 = node_xy[1+(ip1-1)*2];
  y2 = node_xy[1+(ip2-1)*2];
  y3 = node_xy[1+(ip3-1)*2];

  xqe[ 1-1] = z1 * x1 + z2 * x2 + z2 * x3;
  yqe[ 1-1] = z1 * y1 + z2 * y2 + z2 * y3;
  xqe[ 2-1] = z2 * x1 + z1 * x2 + z2 * x3;
  yqe[ 2-1] = z2 * y1 + z1 * y2 + z2 * y3;
  xqe[ 3-1] = z2 * x1 + z2 * x2 + z1 * x3;
  yqe[ 3-1] = z2 * y1 + z2 * y2 + z1 * y3;
  xqe[ 4-1] = z3 * x1 + z4 * x2 + z4 * x3;
  yqe[ 4-1] = z3 * y1 + z4 * y2 + z4 * y3;
  xqe[ 5-1] = z4 * x1 + z3 * x2 + z4 * x3;
  yqe[ 5-1] = z4 * y1 + z3 * y2 + z4 * y3;
  xqe[ 6-1] = z4 * x1 + z4 * x2 + z3 * x3;
  yqe[ 6-1] = z4 * y1 + z4 * y2 + z3 * y3;
  xqe[ 7-1] = z5 * x1 + z6 * x2 + z7 * x3;
  yqe[ 7-1] = z5 * y1 + z6 * y2 + z7 * y3;
  xqe[ 8-1] = z5 * x1 + z7 * x2 + z6 * x3;
  yqe[ 8-1] = z5 * y1 + z7 * y2 + z6 * y3;
  xqe[ 9-1] = z6 * x1 + z5 * x2 + z7 * x3;
  yqe[ 9-1] = z6 * y1 + z5 * y2 + z7 * y3;
  xqe[10-1] = z6 * x1 + z7 * x2 + z5 * x3;
  yqe[10-1] = z6 * y1 + z7 * y2 + z5 * y3;
  xqe[11-1] = z7 * x1 + z5 * x2 + z6 * x3;
  yqe[11-1] = z7 * y1 + z5 * y2 + z6 * y3;
  xqe[12-1] = z7 * x1 + z6 * x2 + z5 * x3;
  yqe[12-1] = z7 * y1 + z6 * y2 + z5 * y3;
  xqe[13-1] = ( x1 + x2 + x3 ) / 3.0;
  yqe[13-1] = ( y1 + y2 + y3 ) / 3.0;
 
  return;
}
//****************************************************************************80

double r8_huge ( void )

//****************************************************************************80
//
//  Purpose:
//
//    R8_HUGE returns a "huge" R8.
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_HUGE, a "huge" R8.
//
{
  return ( double ) HUGE_VAL;
}
//****************************************************************************80

double r8_max ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_MAX returns the maximum of two R8's.
//
//  Modified:
//
//    10 January 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the quantities to compare.
//
//    Output, double R8_MAX, the maximum of X and Y.
//
{
  if ( y < x )
  {
    return x;
  } 
  else
  {
    return y;
  }
}
//****************************************************************************80