src/finley/Quadrature.c

/* $Id$ */

/**************************************************************/

/*   Finley:  */

/**************************************************************/

/*   Copyrights by ACcESS Australia 2003 */
/*   Author: gross@access.edu.au */
/*   Version: $Id$ */

/**************************************************************/

#include "Common.h"
#include "Finley.h"
#include "Quadrature.h"

#define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
#define QUADWEIGHTS(_I_) quadWeights[_I_]

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes for the tri  */
/*   as a queezed scheme on a quad [0,1]^2 */

void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
  int i;
  double Q1,Q2;
  #define DIM 2
  
  /*  the easy case: */
  
  if (numQuadNodes==1) {
    QUADNODES(0,0)=1./3.;
    QUADNODES(1,0)=1./3.;
    QUADWEIGHTS(0)= .5;
  } else {
    
    /*  get scheme on [0.1]^2 */
    
    Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
    if (Finley_ErrorCode!=NO_ERROR) return;
    
    /*  squeeze it: */
    
    for (i=0;i<numQuadNodes;i++) {
        Q1=QUADNODES(0,i);
        Q2=QUADNODES(1,i);
        QUADWEIGHTS(i)=QUADWEIGHTS(i)*(1.-(1./2.)*(Q1+Q2));
        QUADNODES(0,i)=Q1*(1.-(1./2.)*Q2);
        QUADNODES(1,i)=Q2*(1.-(1./2.)*Q1);
    }
  }
  #undef DIM
}

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
/*   as a queezed scheme on a hex [0,1]^3 */

void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
  int i;
  double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
  #define DIM 3
  
  /*  the easy case: */
  
  if (numQuadNodes==1) {
    QUADNODES(0,0)= .25;
    QUADNODES(1,0)= .25;
    QUADNODES(2,0)= .25;
    QUADWEIGHTS(0)=1./6.;
  } else {
    
    /*  get scheme on [0.1]^3 */
    
    Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
    if (Finley_ErrorCode!=NO_ERROR) return;
    
    /*  squeeze it: */
    
    for (i=0;i<numQuadNodes;i++) {
      Q1=QUADNODES(0,i);
      Q2=QUADNODES(1,i);
      Q3=QUADNODES(2,i);

      JA11= (1./3.)*Q2*Q3-(1./2.)*(Q2+Q3) +1.;
      JA12= (1./3.)*Q1*Q3-(1./2.)*Q1;
      JA13= (1./3.)*Q1*Q2-(1./2.)*Q1;
      JA21= (1./3.)*Q2*Q3-(1./2.)*Q2;
      JA22= (1./3.)*Q1*Q3-(1./2.)*(Q1+Q3) +1.;
      JA23= (1./3.)*Q1*Q2-(1./2.)*Q2;
      JA31= (1./3.)*Q2*Q3-(1./2.)*Q3;
      JA32= (1./3.)*Q1*Q3-(1./2.)*Q3;
      JA33= (1./3.)*Q1*Q2-(1./2.)*(Q1+Q2) +1.;
      DET=JA11*JA22*JA33+JA12*JA23*JA31+JA13*JA21*JA32-JA13*JA22*JA31-JA11*JA23*JA32-JA12*JA21*JA33;
      quadWeights[i]=quadWeights[i]*ABS(DET);
      QUADNODES(0,i)=Q1*((1./3.)*Q2*Q3-(1./2.)*(Q2+Q3)+1.);
      QUADNODES(1,i)=Q2*((1./3.)*Q1*Q3-(1./2.)*(Q1+Q3)+1.);
      QUADNODES(2,i)=Q3*((1./3.)*Q1*Q2-(1./2.)*(Q1+Q2)+1.);
    }
  }
  #undef DIM
}

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
/*   as a X-product of a 1D scheme. */

void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
  int numQuadNodes1d,i,j,l;
  double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
  #define DIM 2
  
  /*  find numQuadNodes1d with numQuadNodes1d**2==numQuadNodes: */
  
  for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
    if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
      
      /*  get 1D scheme: */
      
      Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
      
      /*  make 2D scheme: */
      
      if (Finley_ErrorCode==NO_ERROR) {
        l=0;
        for (i=0;i<numQuadNodes1d;i++) {
          for (j=0;j<numQuadNodes1d;j++) {
            QUADNODES(0,l)=quadNodes1d[i];
            QUADNODES(1,l)=quadNodes1d[j];
            QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
            l++;
          }
        }
      }
      return;
    }
  }
  Finley_ErrorCode=VALUE_ERROR;
  sprintf(Finley_ErrorMsg,"Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
  #undef DIM
}

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
/*   as a X-product of a 1D scheme. */

void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
  int numQuadNodes1d,i,j,k,l;
  double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
  #define DIM 3
  
  /*  find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
  
  for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
    if (numQuadNodes1d*numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
      
      /*  get 1D scheme: */
      
      Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
      
      /*  make 3D scheme: */
      
      if (Finley_ErrorCode==NO_ERROR) {
        l=0;
        for (i=0;i<numQuadNodes1d;i++) {
          for (j=0;j<numQuadNodes1d;j++) {
            for (k=0;k<numQuadNodes1d;k++) {
              QUADNODES(0,l)=quadNodes1d[i];
              QUADNODES(1,l)=quadNodes1d[j];
              QUADNODES(2,l)=quadNodes1d[k];
              QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j]*quadWeights1d[k];
              l++;
            }
          }
        }
      }
      
      return;
    }
  }
  Finley_ErrorCode=VALUE_ERROR;
  sprintf(Finley_ErrorMsg,"Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
  #undef DIM
}

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
/*   in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
/*   an error. */

void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
  if (numQuadNodes!=0) {
       Finley_ErrorCode=VALUE_ERROR;
       sprintf(Finley_ErrorMsg,"There is no quadrature scheme on points.");
  }
}

/**************************************************************/

/*   get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
/*   The nodes and weights are set from a table. */

void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
  switch(numQuadNodes) {
      case 1:
        quadNodes[0]=0.5;
        quadWeights[0]=1.;
        break;

      case 2:
        quadNodes[0]=(1.-.577350269189626)/2.;
        quadNodes[1]=(1.+.577350269189626)/2.;
        quadWeights[0]=.5;
        quadWeights[1]=.5;
        break;

      case 3:
        quadNodes[0]=(1.-.774596669241483)/2.;
        quadNodes[1]=.5;
        quadNodes[2]=(1.+.774596669241483)/2.;
        quadWeights[0]=5./18.;
        quadWeights[1]=4./ 9.;
        quadWeights[2]=5./18.;
        break;

      case 4:
        quadNodes[0]=(1.-.861136311594053)/2.;
        quadNodes[1]=(1.-.339981043584856)/2.;
        quadNodes[2]=(1.+.339981043584856)/2.;
        quadNodes[3]=(1.+.861136311594053)/2.;
        quadWeights[0]=.347854845137454/2.;
        quadWeights[1]=.652145154862546/2.;
        quadWeights[2]=.652145154862546/2.;
        quadWeights[3]=.347854845137454/2.;
        break;

      case 5:
        quadNodes[0]=(1.-.906179845938664)/2.;
        quadNodes[1]=(1.-.538469310105683)/2.;
        quadNodes[2]= .5;
        quadNodes[3]=(1.+.538469310105683)/2.;
        quadNodes[4]=(1.+.906179845938664)/2.;
        quadWeights[0]=.236926885056189/2.;
        quadWeights[1]=.478628670499366/2.;
        quadWeights[2]=.568888888888889/2.;
        quadWeights[3]=.478628670499366/2.;
        quadWeights[4]=.236926885056189/2.;
        break;

      case 6:
        quadNodes[0]=(1.-.932469514203152)/2.;
        quadNodes[1]=(1.-.661209386466265)/2.;
        quadNodes[2]=(1.-.238619186083197)/2.;
        quadNodes[3]=(1.+.238619186083197)/2.;
        quadNodes[4]=(1.+.661209386466265)/2.;
        quadNodes[5]=(1.+.932469514203152)/2.;
        quadWeights[0]=.171324492379170/2.;
        quadWeights[1]=.360761573048139/2.;
        quadWeights[2]=.467913934572691/2.;
        quadWeights[3]=.467913934572691/2.;
        quadWeights[4]=.360761573048139/2.;
        quadWeights[5]=.171324492379170/2.;
        break;

      case 7:
        quadNodes[0]=(1.-.949107912342759)/2.;
        quadNodes[1]=(1.-.741531185599394)/2.;
        quadNodes[2]=(1.-.405845151377397)/2.;
        quadNodes[3]=0.5;
        quadNodes[4]=(1.+.405845151377397)/2.;
        quadNodes[5]=(1.+.741531185599394)/2.;
        quadNodes[6]=(1.+.949107912342759)/2.;
        quadWeights[0]= .129484966168870/2.;
        quadWeights[1]= .279705391489277/2.;
        quadWeights[2]= .381830050505119/2.;
        quadWeights[3]= .417959183673469/2.;
        quadWeights[4]= .381830050505119/2.;
        quadWeights[5]= .279705391489277/2.;
        quadWeights[6]= .129484966168870/2.;
        break;

      case 8:
        quadNodes[0]=(1.-.960289856497536)/2.;
        quadNodes[1]=(1.-.796666477413627)/2.;
        quadNodes[2]=(1.-.525532409916329)/2.;
        quadNodes[3]=(1.-.183434642495650)/2.;
        quadNodes[4]=(1.+.183434642495650)/2.;
        quadNodes[5]=(1.+.525532409916329)/2.;
        quadNodes[6]=(1.+.796666477413627)/2.;
        quadNodes[7]=(1.+.960289856497536)/2.;
        quadWeights[0]= .101228536290376/2.;
        quadWeights[1]= .222381034453374/2.;
        quadWeights[2]= .313706645877887/2.;
        quadWeights[3]= .362683783378362/2.;
        quadWeights[4]= .362683783378362/2.;
        quadWeights[5]= .313706645877887/2.;
        quadWeights[6]= .222381034453374/2.;
        quadWeights[7]= .101228536290376/2.;
        break;

      case 9:
        quadNodes[0]=(1.-.968160239507626)/2.;
        quadNodes[1]=(1.-.836031107326636)/2.;
        quadNodes[2]=(1.-.613371432700590)/2.;
        quadNodes[3]=(1.-.324253423403809)/2.;
        quadNodes[4]= .5;
        quadNodes[5]=(1.+.324253423403809)/2.;
        quadNodes[6]=(1.+.613371432700590)/2.;
        quadNodes[7]=(1.+.836031107326636)/2.;
        quadNodes[8]=(1.+.968160239507626)/2.;
        quadWeights[0]= .081274388361574/2.;
        quadWeights[1]= .180648160694857/2.;
        quadWeights[2]= .260610696402935/2.;
        quadWeights[3]= .312347077040003/2.;
        quadWeights[4]= .330239355001260/2.;
        quadWeights[5]= .312347077040003/2.;
        quadWeights[6]= .260610696402935/2.;
        quadWeights[7]= .180648160694857/2.;
        quadWeights[8]= .081274388361574/2.;
        break;

      case 10:
        quadNodes[0]=(1.-.973906528517172)/2.;
        quadNodes[1]=(1.-.865063366688985)/2.;
        quadNodes[2]=(1.-.679409568299024)/2.;
        quadNodes[3]=(1.-.433395394129247)/2.;
        quadNodes[4]=(1.-.148874338981631)/2.;
        quadNodes[5]=(1.+.148874338981631)/2.;
        quadNodes[6]=(1.+.433395394129247)/2.;
        quadNodes[7]=(1.+.679409568299024)/2.;
        quadNodes[8]=(1.+.865063366688985)/2.;
        quadNodes[9]=(1.+.973906528517172)/2.;
        quadWeights[0]= .066671344308688/2.;
        quadWeights[1]= .149451349150581/2.;
        quadWeights[2]= .219086362515982/2.;
        quadWeights[3]= .269266719309996/2.;
        quadWeights[4]= .295524224714753/2.;
        quadWeights[5]= .295524224714753/2.;
        quadWeights[6]= .269266719309996/2.;
        quadWeights[7]= .219086362515982/2.;
        quadWeights[8]= .149451349150581/2.;
        quadWeights[9]= .066671344308688/2.;
        break;

      default:
        Finley_ErrorCode=VALUE_ERROR;
        sprintf(Finley_ErrorMsg,"Negative intergration order.");
        break;
  }
}
/**************************************************************/
/*                                                            */
/*  the following function are used define the meshes on the surface in the xy-plane */

/* triangle surface on a tetrahedron */
void Finley_Quad_getNodesTriOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
       Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesTri);
}
/* rectangular surface on a hexahedron */
void Finley_Quad_getNodesRecOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
       Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesRec);
}
/* line surface on a triangle or rectangle */
void Finley_Quad_getNodesLineOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
       Finley_Quad_makeNodesOnFace(2,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesLine);
}
/* point surface on a line */
void Finley_Quad_getNodesPointOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
     Finley_Quad_makeNodesOnFace(1,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesPoint);
}

void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
    int q,i;
    double quadNodesOnFace[numQuadNodes*(dim-1)];

    #define DIM dim
    getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
    if (Finley_ErrorCode!=NO_ERROR) return;
    
    for (q=0;q<numQuadNodes;q++) {
       for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
       QUADNODES(dim-1,q)=0;
    }
    #undef DIM
}

/**************************************************************/

/*    The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
/*    achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
/*    to consider the accuracy reduction through the construction process.  */


int Finley_Quad_getNumNodesPoint(int order) {
  if (order <0 ) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Negative intergration order.");
    return -1;
  } else { 
    return 0;
  }
}

int Finley_Quad_getNumNodesLine(int order) {
  if (order <0 ) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Negative intergration order.");
    return -1;
  } else { 
    if (order > 2*MAX_numQuadNodesLine-1) {
      Finley_ErrorCode=VALUE_ERROR;
      sprintf(Finley_ErrorMsg,"requested integration order %d on line is too large (>%d).",
                                                           order,2*MAX_numQuadNodesLine-1);
      return -1;
    } else { 
       Finley_ErrorCode=NO_ERROR;
       return order/2+1;
    }
  }
}

int Finley_Quad_getNumNodesTri(int order) {
  int numQuadNodesLine;
  if (order==1) {
      return 1;
  } else {
      numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
      if (Finley_ErrorCode==NO_ERROR) {
         return numQuadNodesLine*numQuadNodesLine;
      } else {
         return -1;
      }
  }
}

int Finley_Quad_getNumNodesRec(int order) {
  int numQuadNodesLine;
  numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
  if (Finley_ErrorCode==NO_ERROR) {
      return numQuadNodesLine*numQuadNodesLine;
  } else {
      return -1;
  }
}

int Finley_Quad_getNumNodesTet(int order) {
  int numQuadNodesLine;
  numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
  if (Finley_ErrorCode==NO_ERROR) {
      return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
  } else {
      return -1;
  }
}

int Finley_Quad_getNumNodesHex(int order) {
  int numQuadNodesLine;
  numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
  if (Finley_ErrorCode==NO_ERROR) {
      return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
  } else {
      return -1;
  }
}
/* 
* $Log$
* Revision 1.1  2004/10/26 06:53:57  jgs
* Initial revision
*
* Revision 1.2  2004/08/03 04:49:06  gross
* bug in Quadrature.c fixed
*
* Revision 1.1.1.1  2004/06/24 04:00:40  johng
* Initial version of eys using boost-python.
*
*
*/
1	jgs	82	/* $Id$ */
2
3			/**************************************************************/
4
5			/* Finley: */
6
7			/**************************************************************/
8
9			/* Copyrights by ACcESS Australia 2003 */
10			/* Author: gross@access.edu.au */
11			/* Version: $Id$ */
12
13			/**************************************************************/
14
15			#include "Common.h"
16			#include "Finley.h"
17			#include "Quadrature.h"
18
19			#define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
20			#define QUADWEIGHTS(_I_) quadWeights[_I_]
21
22			/**************************************************************/
23
24			/* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */
25			/* as a queezed scheme on a quad [0,1]^2 */
26
27			void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
28			int i;
29			double Q1,Q2;
30			#define DIM 2
31
32			/* the easy case: */
33
34			if (numQuadNodes==1) {
35			QUADNODES(0,0)=1./3.;
36			QUADNODES(1,0)=1./3.;
37			QUADWEIGHTS(0)= .5;
38			} else {
39
40			/* get scheme on [0.1]^2 */
41
42			Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
43			if (Finley_ErrorCode!=NO_ERROR) return;
44
45			/* squeeze it: */
46
47			for (i=0;i<numQuadNodes;i++) {
48			Q1=QUADNODES(0,i);
49			Q2=QUADNODES(1,i);
50			QUADWEIGHTS(i)=QUADWEIGHTS(i)(1.-(1./2.)(Q1+Q2));
51			QUADNODES(0,i)=Q1(1.-(1./2.)Q2);
52			QUADNODES(1,i)=Q2(1.-(1./2.)Q1);
53			}
54			}
55			#undef DIM
56			}
57
58			/**************************************************************/
59
60			/* get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
61			/* as a queezed scheme on a hex [0,1]^3 */
62
63			void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
64			int i;
65			double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
66			#define DIM 3
67
68			/* the easy case: */
69
70			if (numQuadNodes==1) {
71			QUADNODES(0,0)= .25;
72			QUADNODES(1,0)= .25;
73			QUADNODES(2,0)= .25;
74			QUADWEIGHTS(0)=1./6.;
75			} else {
76
77			/* get scheme on [0.1]^3 */
78
79			Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
80			if (Finley_ErrorCode!=NO_ERROR) return;
81
82			/* squeeze it: */
83
84			for (i=0;i<numQuadNodes;i++) {
85			Q1=QUADNODES(0,i);
86			Q2=QUADNODES(1,i);
87			Q3=QUADNODES(2,i);
88
89			JA11= (1./3.)Q2Q3-(1./2.)*(Q2+Q3) +1.;
90			JA12= (1./3.)Q1Q3-(1./2.)*Q1;
91			JA13= (1./3.)Q1Q2-(1./2.)*Q1;
92			JA21= (1./3.)Q2Q3-(1./2.)*Q2;
93			JA22= (1./3.)Q1Q3-(1./2.)*(Q1+Q3) +1.;
94			JA23= (1./3.)Q1Q2-(1./2.)*Q2;
95			JA31= (1./3.)Q2Q3-(1./2.)*Q3;
96			JA32= (1./3.)Q1Q3-(1./2.)*Q3;
97			JA33= (1./3.)Q1Q2-(1./2.)*(Q1+Q2) +1.;
98			DET=JA11JA22JA33+JA12JA23JA31+JA13JA21JA32-JA13JA22JA31-JA11JA23JA32-JA12JA21JA33;
99			quadWeights[i]=quadWeights[i]*ABS(DET);
100			QUADNODES(0,i)=Q1((1./3.)Q2Q3-(1./2.)(Q2+Q3)+1.);
101			QUADNODES(1,i)=Q2((1./3.)Q1Q3-(1./2.)(Q1+Q3)+1.);
102			QUADNODES(2,i)=Q3((1./3.)Q1Q2-(1./2.)(Q1+Q2)+1.);
103			}
104			}
105			#undef DIM
106			}
107
108			/**************************************************************/
109
110			/* get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
111			/* as a X-product of a 1D scheme. */
112
113			void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
114			int numQuadNodes1d,i,j,l;
115			double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
116			#define DIM 2
117
118			/* find numQuadNodes1d with numQuadNodes1d*2==numQuadNodes: /
119
120			for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
121			if (numQuadNodes1d*numQuadNodes1d==numQuadNodes) {
122
123			/* get 1D scheme: */
124
125			Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
126
127			/* make 2D scheme: */
128
129			if (Finley_ErrorCode==NO_ERROR) {
130			l=0;
131			for (i=0;i<numQuadNodes1d;i++) {
132			for (j=0;j<numQuadNodes1d;j++) {
133			QUADNODES(0,l)=quadNodes1d[i];
134			QUADNODES(1,l)=quadNodes1d[j];
135			QUADWEIGHTS(l)=quadWeights1d[i]*quadWeights1d[j];
136			l++;
137			}
138			}
139			}
140			return;
141			}
142			}
143			Finley_ErrorCode=VALUE_ERROR;
144			sprintf(Finley_ErrorMsg,"Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
145			#undef DIM
146			}
147
148			/**************************************************************/
149
150			/* get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
151			/* as a X-product of a 1D scheme. */
152
153			void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
154			int numQuadNodes1d,i,j,k,l;
155			double quadNodes1d[numQuadNodes],quadWeights1d[numQuadNodes];
156			#define DIM 3
157
158			/* find numQuadNodes1d with numQuadNodes1d*3==numQuadNodes: /
159
160			for (numQuadNodes1d=1;numQuadNodes1d<=MAX_numQuadNodesLine;numQuadNodes1d++) {
161			if (numQuadNodes1dnumQuadNodes1dnumQuadNodes1d==numQuadNodes) {
162
163			/* get 1D scheme: */
164
165			Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
166
167			/* make 3D scheme: */
168
169			if (Finley_ErrorCode==NO_ERROR) {
170			l=0;
171			for (i=0;i<numQuadNodes1d;i++) {
172			for (j=0;j<numQuadNodes1d;j++) {
173			for (k=0;k<numQuadNodes1d;k++) {
174			QUADNODES(0,l)=quadNodes1d[i];
175			QUADNODES(1,l)=quadNodes1d[j];
176			QUADNODES(2,l)=quadNodes1d[k];
177			QUADWEIGHTS(l)=quadWeights1d[i]quadWeights1d[j]quadWeights1d[k];
178			l++;
179			}
180			}
181			}
182			}
183
184			return;
185			}
186			}
187			Finley_ErrorCode=VALUE_ERROR;
188			sprintf(Finley_ErrorMsg,"Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
189			#undef DIM
190			}
191
192			/**************************************************************/
193
194			/* get a quadrature scheme with numQuadNodes quadrature nodes for a point. As there */
195			/* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
196			/* an error. */
197
198			void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
199			if (numQuadNodes!=0) {
200			Finley_ErrorCode=VALUE_ERROR;
201			sprintf(Finley_ErrorMsg,"There is no quadrature scheme on points.");
202			}
203			}
204
205			/**************************************************************/
206
207			/* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
208			/* The nodes and weights are set from a table. */
209
210			void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
211			switch(numQuadNodes) {
212			case 1:
213			quadNodes[0]=0.5;
214			quadWeights[0]=1.;
215			break;
216
217			case 2:
218			quadNodes[0]=(1.-.577350269189626)/2.;
219			quadNodes[1]=(1.+.577350269189626)/2.;
220			quadWeights[0]=.5;
221			quadWeights[1]=.5;
222			break;
223
224			case 3:
225			quadNodes[0]=(1.-.774596669241483)/2.;
226			quadNodes[1]=.5;
227			quadNodes[2]=(1.+.774596669241483)/2.;
228			quadWeights[0]=5./18.;
229			quadWeights[1]=4./ 9.;
230			quadWeights[2]=5./18.;
231			break;
232
233			case 4:
234			quadNodes[0]=(1.-.861136311594053)/2.;
235			quadNodes[1]=(1.-.339981043584856)/2.;
236			quadNodes[2]=(1.+.339981043584856)/2.;
237			quadNodes[3]=(1.+.861136311594053)/2.;
238			quadWeights[0]=.347854845137454/2.;
239			quadWeights[1]=.652145154862546/2.;
240			quadWeights[2]=.652145154862546/2.;
241			quadWeights[3]=.347854845137454/2.;
242			break;
243
244			case 5:
245			quadNodes[0]=(1.-.906179845938664)/2.;
246			quadNodes[1]=(1.-.538469310105683)/2.;
247			quadNodes[2]= .5;
248			quadNodes[3]=(1.+.538469310105683)/2.;
249			quadNodes[4]=(1.+.906179845938664)/2.;
250			quadWeights[0]=.236926885056189/2.;
251			quadWeights[1]=.478628670499366/2.;
252			quadWeights[2]=.568888888888889/2.;
253			quadWeights[3]=.478628670499366/2.;
254			quadWeights[4]=.236926885056189/2.;
255			break;
256
257			case 6:
258			quadNodes[0]=(1.-.932469514203152)/2.;
259			quadNodes[1]=(1.-.661209386466265)/2.;
260			quadNodes[2]=(1.-.238619186083197)/2.;
261			quadNodes[3]=(1.+.238619186083197)/2.;
262			quadNodes[4]=(1.+.661209386466265)/2.;
263			quadNodes[5]=(1.+.932469514203152)/2.;
264			quadWeights[0]=.171324492379170/2.;
265			quadWeights[1]=.360761573048139/2.;
266			quadWeights[2]=.467913934572691/2.;
267			quadWeights[3]=.467913934572691/2.;
268			quadWeights[4]=.360761573048139/2.;
269			quadWeights[5]=.171324492379170/2.;
270			break;
271
272			case 7:
273			quadNodes[0]=(1.-.949107912342759)/2.;
274			quadNodes[1]=(1.-.741531185599394)/2.;
275			quadNodes[2]=(1.-.405845151377397)/2.;
276			quadNodes[3]=0.5;
277			quadNodes[4]=(1.+.405845151377397)/2.;
278			quadNodes[5]=(1.+.741531185599394)/2.;
279			quadNodes[6]=(1.+.949107912342759)/2.;
280			quadWeights[0]= .129484966168870/2.;
281			quadWeights[1]= .279705391489277/2.;
282			quadWeights[2]= .381830050505119/2.;
283			quadWeights[3]= .417959183673469/2.;
284			quadWeights[4]= .381830050505119/2.;
285			quadWeights[5]= .279705391489277/2.;
286			quadWeights[6]= .129484966168870/2.;
287			break;
288
289			case 8:
290			quadNodes[0]=(1.-.960289856497536)/2.;
291			quadNodes[1]=(1.-.796666477413627)/2.;
292			quadNodes[2]=(1.-.525532409916329)/2.;
293			quadNodes[3]=(1.-.183434642495650)/2.;
294			quadNodes[4]=(1.+.183434642495650)/2.;
295			quadNodes[5]=(1.+.525532409916329)/2.;
296			quadNodes[6]=(1.+.796666477413627)/2.;
297			quadNodes[7]=(1.+.960289856497536)/2.;
298			quadWeights[0]= .101228536290376/2.;
299			quadWeights[1]= .222381034453374/2.;
300			quadWeights[2]= .313706645877887/2.;
301			quadWeights[3]= .362683783378362/2.;
302			quadWeights[4]= .362683783378362/2.;
303			quadWeights[5]= .313706645877887/2.;
304			quadWeights[6]= .222381034453374/2.;
305			quadWeights[7]= .101228536290376/2.;
306			break;
307
308			case 9:
309			quadNodes[0]=(1.-.968160239507626)/2.;
310			quadNodes[1]=(1.-.836031107326636)/2.;
311			quadNodes[2]=(1.-.613371432700590)/2.;
312			quadNodes[3]=(1.-.324253423403809)/2.;
313			quadNodes[4]= .5;
314			quadNodes[5]=(1.+.324253423403809)/2.;
315			quadNodes[6]=(1.+.613371432700590)/2.;
316			quadNodes[7]=(1.+.836031107326636)/2.;
317			quadNodes[8]=(1.+.968160239507626)/2.;
318			quadWeights[0]= .081274388361574/2.;
319			quadWeights[1]= .180648160694857/2.;
320			quadWeights[2]= .260610696402935/2.;
321			quadWeights[3]= .312347077040003/2.;
322			quadWeights[4]= .330239355001260/2.;
323			quadWeights[5]= .312347077040003/2.;
324			quadWeights[6]= .260610696402935/2.;
325			quadWeights[7]= .180648160694857/2.;
326			quadWeights[8]= .081274388361574/2.;
327			break;
328
329			case 10:
330			quadNodes[0]=(1.-.973906528517172)/2.;
331			quadNodes[1]=(1.-.865063366688985)/2.;
332			quadNodes[2]=(1.-.679409568299024)/2.;
333			quadNodes[3]=(1.-.433395394129247)/2.;
334			quadNodes[4]=(1.-.148874338981631)/2.;
335			quadNodes[5]=(1.+.148874338981631)/2.;
336			quadNodes[6]=(1.+.433395394129247)/2.;
337			quadNodes[7]=(1.+.679409568299024)/2.;
338			quadNodes[8]=(1.+.865063366688985)/2.;
339			quadNodes[9]=(1.+.973906528517172)/2.;
340			quadWeights[0]= .066671344308688/2.;
341			quadWeights[1]= .149451349150581/2.;
342			quadWeights[2]= .219086362515982/2.;
343			quadWeights[3]= .269266719309996/2.;
344			quadWeights[4]= .295524224714753/2.;
345			quadWeights[5]= .295524224714753/2.;
346			quadWeights[6]= .269266719309996/2.;
347			quadWeights[7]= .219086362515982/2.;
348			quadWeights[8]= .149451349150581/2.;
349			quadWeights[9]= .066671344308688/2.;
350			break;
351
352			default:
353			Finley_ErrorCode=VALUE_ERROR;
354			sprintf(Finley_ErrorMsg,"Negative intergration order.");
355			break;
356			}
357			}
358			/**************************************************************/
359			/* */
360			/* the following function are used define the meshes on the surface in the xy-plane */
361
362			/* triangle surface on a tetrahedron */
363			void Finley_Quad_getNodesTriOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
364			Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesTri);
365			}
366			/* rectangular surface on a hexahedron */
367			void Finley_Quad_getNodesRecOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
368			Finley_Quad_makeNodesOnFace(3,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesRec);
369			}
370			/* line surface on a triangle or rectangle */
371			void Finley_Quad_getNodesLineOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
372			Finley_Quad_makeNodesOnFace(2,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesLine);
373			}
374			/* point surface on a line */
375			void Finley_Quad_getNodesPointOnFace(int numQuadNodes,double* quadNodes,double* quadWeights) {
376			Finley_Quad_makeNodesOnFace(1,numQuadNodes,quadNodes,quadWeights,Finley_Quad_getNodesPoint);
377			}
378
379			void Finley_Quad_makeNodesOnFace(int dim, int numQuadNodes,double* quadNodes,double* quadWeights, Finley_Quad_getNodes getFaceNodes) {
380			int q,i;
381			double quadNodesOnFace[numQuadNodes*(dim-1)];
382
383			#define DIM dim
384			getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);
385			if (Finley_ErrorCode!=NO_ERROR) return;
386
387			for (q=0;q<numQuadNodes;q++) {
388			for (i=0;i<dim-1;i++) QUADNODES(i,q)=quadNodesOnFace[INDEX2(i,q,dim-1)];
389			QUADNODES(dim-1,q)=0;
390			}
391			#undef DIM
392			}
393
394			/**************************************************************/
395
396			/* The following functions Finley_Quad_getNumNodes* return the nmber of quadrature points needed to */
397			/* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
398			/* to consider the accuracy reduction through the construction process. */
399
400
401			int Finley_Quad_getNumNodesPoint(int order) {
402			if (order <0 ) {
403			Finley_ErrorCode=VALUE_ERROR;
404			sprintf(Finley_ErrorMsg,"Negative intergration order.");
405			return -1;
406			} else {
407			return 0;
408			}
409			}
410
411			int Finley_Quad_getNumNodesLine(int order) {
412			if (order <0 ) {
413			Finley_ErrorCode=VALUE_ERROR;
414			sprintf(Finley_ErrorMsg,"Negative intergration order.");
415			return -1;
416			} else {
417			if (order > 2*MAX_numQuadNodesLine-1) {
418			Finley_ErrorCode=VALUE_ERROR;
419			sprintf(Finley_ErrorMsg,"requested integration order %d on line is too large (>%d).",
420			order,2*MAX_numQuadNodesLine-1);
421			return -1;
422			} else {
423			Finley_ErrorCode=NO_ERROR;
424			return order/2+1;
425			}
426			}
427			}
428
429			int Finley_Quad_getNumNodesTri(int order) {
430			int numQuadNodesLine;
431			if (order==1) {
432			return 1;
433			} else {
434			numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);
435			if (Finley_ErrorCode==NO_ERROR) {
436			return numQuadNodesLine*numQuadNodesLine;
437			} else {
438			return -1;
439			}
440			}
441			}
442
443			int Finley_Quad_getNumNodesRec(int order) {
444			int numQuadNodesLine;
445			numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
446			if (Finley_ErrorCode==NO_ERROR) {
447			return numQuadNodesLine*numQuadNodesLine;
448			} else {
449			return -1;
450			}
451			}
452
453			int Finley_Quad_getNumNodesTet(int order) {
454			int numQuadNodesLine;
455			numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);
456			if (Finley_ErrorCode==NO_ERROR) {
457			return numQuadNodesLinenumQuadNodesLinenumQuadNodesLine;
458			} else {
459			return -1;
460			}
461			}
462
463			int Finley_Quad_getNumNodesHex(int order) {
464			int numQuadNodesLine;
465			numQuadNodesLine=Finley_Quad_getNumNodesLine(order);
466			if (Finley_ErrorCode==NO_ERROR) {
467			return numQuadNodesLinenumQuadNodesLinenumQuadNodesLine;
468			} else {
469			return -1;
470			}
471			}
472			/*
473			* $Log$
474			* Revision 1.1 2004/10/26 06:53:57 jgs
475			* Initial revision
476			*
477			* Revision 1.2 2004/08/03 04:49:06 gross
478			* bug in Quadrature.c fixed
479			*
480			* Revision 1.1.1.1 2004/06/24 04:00:40 johng
481			* Initial version of eys using boost-python.
482			*
483			*
484			*/
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision