finley/src/Quadrature.c

/*
 ************************************************************
 *          Copyright 2006 by ACcESS MNRF                   *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       Primary Business: Queensland, Australia            *
 *  Licensed under the Open Software License version 3.0    *
 *     http://www.opensource.org/licenses/osl-3.0.php       *
 *                                                          *
 ************************************************************
*/

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

/*   Finley:  */

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

/*   Author: gross@access.edu.au */
/*   Version: $Id$ */

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

#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_noError()) 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_noError()) 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) {
  char error_msg[LenErrorMsg_MAX];
  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_noError()) {
        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;
    }
  }
  sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
  Finley_setError(VALUE_ERROR,error_msg);
  #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) {
  char error_msg[LenErrorMsg_MAX];
  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_noError()) {
        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;
    }
  }
  sprintf(error_msg,"__FILE__: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
  Finley_setError(VALUE_ERROR,error_msg);
  #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_setError(VALUE_ERROR,"__FILE__: 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_setError(VALUE_ERROR,"__FILE__: 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_noError()) 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_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
    return -1;
  } else { 
    return 0;
  }
}

int Finley_Quad_getNumNodesLine(int order) {
  char error_msg[LenErrorMsg_MAX];
  if (order <0 ) {
    Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");
    return -1;
  } else { 
    if (order > 2*MAX_numQuadNodesLine-1) {
      sprintf(error_msg,"__FILE__: requested integration order %d on line is too large (>%d).",
                                                           order,2*MAX_numQuadNodesLine-1);
      Finley_setError(VALUE_ERROR,error_msg);
      return -1;
    } else { 
       Finley_resetError();
       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_noError()) {
         return numQuadNodesLine*numQuadNodesLine;
      } else {
         return -1;
      }
  }
}

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

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

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