/[escript]/trunk/finley/src/Quadrature.cpp
ViewVC logotype

Diff of /trunk/finley/src/Quadrature.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4491 by jfenwick, Tue Apr 2 04:46:45 2013 UTC revision 4492 by caltinay, Tue Jul 2 01:44:11 2013 UTC
# Line 14  Line 14 
14  *****************************************************************************/  *****************************************************************************/
15    
16    
17  /************************************************************************************/  /****************************************************************************
18    
19  /*   Finley: quadrature schemes */    Finley: quadrature schemes
20    
21  /************************************************************************************/  *****************************************************************************/
22    
23  #include "Quadrature.h"  #include "Quadrature.h"
24  #include "esysUtils/index.h"  #include "esysUtils/index.h"
# Line 28  Line 28 
28  #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]  #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
29  #define QUADWEIGHTS(_I_) quadWeights[_I_]  #define QUADWEIGHTS(_I_) quadWeights[_I_]
30    
31  /************************************************************************************/  namespace finley {
32    
33  Finley_QuadInfo Finley_QuadInfoList[]={  QuadInfo QuadInfoList[]={
34      {PointQuad, "Point", 0,  1,     Finley_Quad_getNodesPoint,      Finley_Quad_getNumNodesPoint,   Finley_Quad_MacroPoint} ,      {PointQuad, "Point", 0,  1,     Quad_getNodesPoint,      Quad_getNumNodesPoint,   Quad_MacroPoint} ,
35          {LineQuad,  "Line",  1,  2, Finley_Quad_getNodesLine,       Finley_Quad_getNumNodesLine,    Finley_Quad_MacroLine} ,          {LineQuad,  "Line",  1,  2, Quad_getNodesLine,       Quad_getNumNodesLine,    Quad_MacroLine} ,
36      {TriQuad,   "Tri",   2,  3,     Finley_Quad_getNodesTri,           Finley_Quad_getNumNodesTri,  Finley_Quad_MacroTri},      {TriQuad,   "Tri",   2,  3,     Quad_getNodesTri,           Quad_getNumNodesTri,  Quad_MacroTri},
37      {RecQuad,   "Rec",   2,  4,     Finley_Quad_getNodesRec,        Finley_Quad_getNumNodesRec,     Finley_Quad_MacroRec},      {RecQuad,   "Rec",   2,  4,     Quad_getNodesRec,        Quad_getNumNodesRec,     Quad_MacroRec},
38      {TetQuad,   "Tet",   3,  4,     Finley_Quad_getNodesTet,        Finley_Quad_getNumNodesTet,     Finley_Quad_MacroTet},      {TetQuad,   "Tet",   3,  4,     Quad_getNodesTet,        Quad_getNumNodesTet,     Quad_MacroTet},
39      {HexQuad,   "Hex",   3,  8,     Finley_Quad_getNodesHex,        Finley_Quad_getNumNodesHex,     Finley_Quad_MacroHex},      {HexQuad,   "Hex",   3,  8,     Quad_getNodesHex,        Quad_getNumNodesHex,     Quad_MacroHex},
40      {NoQuad, "NoType", 0,  1,   Finley_Quad_getNodesPoint,      Finley_Quad_getNumNodesPoint,   Finley_Quad_MacroPoint}      {NoQuad, "NoType", 0,  1,   Quad_getNodesPoint,      Quad_getNumNodesPoint,   Quad_MacroPoint}
41  };  };
42    
43  Finley_QuadInfo* Finley_QuadInfo_getInfo(Finley_QuadTypeId id)  QuadInfo* QuadInfo_getInfo(QuadTypeId id)
44  {  {
45      int ptr=0;      int ptr=0;
46      Finley_QuadInfo* out=NULL;      QuadInfo* out=NULL;
47      while (Finley_QuadInfoList[ptr].TypeId!=NoQuad && out==NULL) {      while (QuadInfoList[ptr].TypeId!=NoQuad && out==NULL) {
48         if (Finley_QuadInfoList[ptr].TypeId==id) out=&(Finley_QuadInfoList[ptr]);         if (QuadInfoList[ptr].TypeId==id) out=&(QuadInfoList[ptr]);
49         ptr++;         ptr++;
50      }      }
51      if (out==NULL) {      if (out==NULL) {
52          Finley_setError(VALUE_ERROR,"Finley_QuadInfo_getInfo: cannot find requested quadrature scheme.");          Finley_setError(VALUE_ERROR,"QuadInfo_getInfo: cannot find requested quadrature scheme.");
53      }      }
54      return out;      return out;
55  }  }
# Line 59  Finley_QuadInfo* Finley_QuadInfo_getInfo Line 59  Finley_QuadInfo* Finley_QuadInfo_getInfo
59  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tri  */  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tri  */
60  /*   as a squeezed scheme on a quad [0,1]^2 */  /*   as a squeezed scheme on a quad [0,1]^2 */
61    
62  void Finley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
63    int i;    int i;
64    double Q1,Q2,a,b,c,d,e,f,g,u,v,w;    double Q1,Q2,a,b,c,d,e,f,g,u,v,w;
65    #define DIM 2    #define DIM 2
# Line 350  void Finley_Quad_getNodesTri(int numQuad Line 350  void Finley_Quad_getNodesTri(int numQuad
350    } else {    } else {
351            
352      /*  get scheme on [0.1]^2 */      /*  get scheme on [0.1]^2 */
353      Finley_Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);      Quad_getNodesRec(numQuadNodes,quadNodes,quadWeights);
354      if (! Finley_noError()) return;      if (! Finley_noError()) return;
355            
356      /*  squeeze it: */      /*  squeeze it: */
# Line 373  void Finley_Quad_getNodesTri(int numQuad Line 373  void Finley_Quad_getNodesTri(int numQuad
373  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tet */  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
374  /*   as a squeezed scheme on a hex [0,1]^3 */  /*   as a squeezed scheme on a hex [0,1]^3 */
375    
376  void Finley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
377    int i;    int i;
378    double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;    double Q1,Q2,Q3,JA11,JA12,JA13,JA21,JA22,JA23,JA31,JA32,JA33,DET;
379    double a,b,c,d,e,f,g,h;    double a,b,c,d,e,f,g,h;
# Line 969  void Finley_Quad_getNodesTet(int numQuad Line 969  void Finley_Quad_getNodesTet(int numQuad
969            
970      /*  get scheme on [0.1]^3 */      /*  get scheme on [0.1]^3 */
971            
972      Finley_Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);      Quad_getNodesHex(numQuadNodes,quadNodes,quadWeights);
973      if (! Finley_noError()) return;      if (! Finley_noError()) return;
974            
975      /*  squeeze it: */      /*  squeeze it: */
# Line 1002  void Finley_Quad_getNodesTet(int numQuad Line 1002  void Finley_Quad_getNodesTet(int numQuad
1002  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the quad [0.1]^2 */
1003  /*   as a X-product of a 1D scheme. */  /*   as a X-product of a 1D scheme. */
1004    
1005  void Finley_Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesRec(int numQuadNodes,double* quadNodes,double* quadWeights) {
1006    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1007    int numQuadNodes1d,i,j,l;    int numQuadNodes1d,i,j,l;
1008    double *quadNodes1d=NULL,*quadWeights1d=NULL;    double *quadNodes1d=NULL,*quadWeights1d=NULL;
# Line 1019  void Finley_Quad_getNodesRec(int numQuad Line 1019  void Finley_Quad_getNodesRec(int numQuad
1019                
1020           /*  get 1D scheme: */           /*  get 1D scheme: */
1021                    
1022           Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);           Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1023                
1024           /*  make 2D scheme: */           /*  make 2D scheme: */
1025                
# Line 1039  void Finley_Quad_getNodesRec(int numQuad Line 1039  void Finley_Quad_getNodesRec(int numQuad
1039         }         }
1040       }       }
1041       if (!set) {       if (!set) {
1042           sprintf(error_msg,"Finley_Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);           sprintf(error_msg,"Quad_getNodesRec: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1043           Finley_setError(VALUE_ERROR,error_msg);           Finley_setError(VALUE_ERROR,error_msg);
1044       }       }
1045       delete[] quadNodes1d;       delete[] quadNodes1d;
# Line 1053  void Finley_Quad_getNodesRec(int numQuad Line 1053  void Finley_Quad_getNodesRec(int numQuad
1053  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the hex [0.1]^3 */
1054  /*   as a X-product of a 1D scheme. */  /*   as a X-product of a 1D scheme. */
1055    
1056  void Finley_Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesHex(int numQuadNodes,double* quadNodes,double* quadWeights) {
1057    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1058    int numQuadNodes1d,i,j,k,l;    int numQuadNodes1d,i,j,k,l;
1059    double *quadNodes1d=NULL,*quadWeights1d=NULL;    double *quadNodes1d=NULL,*quadWeights1d=NULL;
# Line 1070  void Finley_Quad_getNodesHex(int numQuad Line 1070  void Finley_Quad_getNodesHex(int numQuad
1070                
1071           /*  get 1D scheme: */           /*  get 1D scheme: */
1072                
1073           Finley_Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);           Quad_getNodesLine(numQuadNodes1d,quadNodes1d,quadWeights1d);
1074                
1075           /*  make 3D scheme: */           /*  make 3D scheme: */
1076                
# Line 1093  void Finley_Quad_getNodesHex(int numQuad Line 1093  void Finley_Quad_getNodesHex(int numQuad
1093         }         }
1094       }       }
1095       if (!set) {       if (!set) {
1096            sprintf(error_msg,"Finley_Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);            sprintf(error_msg,"Quad_getNodesHex: Illegal number of quadrature nodes %d on hexahedron.",numQuadNodes);
1097            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1098       }       }
1099       delete[] quadNodes1d;       delete[] quadNodes1d;
# Line 1108  void Finley_Quad_getNodesHex(int numQuad Line 1108  void Finley_Quad_getNodesHex(int numQuad
1108  /*   in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */  /*   in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
1109  /*   an error. */  /*   an error. */
1110    
1111  void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
1112    if (numQuadNodes>=0) {    if (numQuadNodes>=0) {
1113          return;          return;
1114    } else {    } else {
1115         Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: Illegal number of quadrature nodes.");         Finley_setError(VALUE_ERROR,"Quad_getNodesPoint: Illegal number of quadrature nodes.");
1116    }    }
1117  }  }
1118    
# Line 1121  void Finley_Quad_getNodesPoint(int numQu Line 1121  void Finley_Quad_getNodesPoint(int numQu
1121  /*   get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */  /*   get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
1122  /*   The nodes and weights are set from a table. */  /*   The nodes and weights are set from a table. */
1123    
1124  void Finley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
1125    switch(numQuadNodes) {    switch(numQuadNodes) {
1126        case 1:        case 1:
1127          quadNodes[0]=0.5;          quadNodes[0]=0.5;
# Line 1264  void Finley_Quad_getNodesLine(int numQua Line 1264  void Finley_Quad_getNodesLine(int numQua
1264          break;          break;
1265    
1266        default:        default:
1267          Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesLine: Negative integration order.");          Finley_setError(VALUE_ERROR,"Quad_getNodesLine: Negative integration order.");
1268          break;          break;
1269    }    }
1270  }  }
# Line 1272  void Finley_Quad_getNodesLine(int numQua Line 1272  void Finley_Quad_getNodesLine(int numQua
1272    
1273  /************************************************************************************/  /************************************************************************************/
1274    
1275  /*    The following functions Finley_Quad_getNumNodes* return the number of quadrature points needed to */  /*    The following functions Quad_getNumNodes* return the number of quadrature points needed to */
1276  /*    achieve a certain accuracy. Notice that for Tet and Tri the order is increased */  /*    achieve a certain accuracy. Notice that for Tet and Tri the order is increased */
1277  /*    to consider the accuracy reduction through the construction process.  */  /*    to consider the accuracy reduction through the construction process.  */
1278    
1279    
1280  int Finley_Quad_getNumNodesPoint(int order) {  int Quad_getNumNodesPoint(int order) {
1281      return 1;      return 1;
1282  }  }
1283    
1284  int Finley_Quad_getNumNodesLine(int order) {  int Quad_getNumNodesLine(int order) {
1285    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1286    if (order <0 ) {    if (order <0 ) {
1287      Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesLine: Negative integration order.");      Finley_setError(VALUE_ERROR,"Quad_getNumNodesLine: Negative integration order.");
1288      return -1;      return -1;
1289    } else {    } else {
1290      if (order > 2*MAX_numQuadNodesLine-1) {      if (order > 2*MAX_numQuadNodesLine-1) {
1291        sprintf(error_msg,"Finley_Quad_getNumNodesLine: requested integration order %d on line is too large (>%d).",        sprintf(error_msg,"Quad_getNumNodesLine: requested integration order %d on line is too large (>%d).",
1292                                                             order,2*MAX_numQuadNodesLine-1);                                                             order,2*MAX_numQuadNodesLine-1);
1293        Finley_setError(VALUE_ERROR,error_msg);        Finley_setError(VALUE_ERROR,error_msg);
1294        return -1;        return -1;
# Line 1299  int Finley_Quad_getNumNodesLine(int orde Line 1299  int Finley_Quad_getNumNodesLine(int orde
1299    }    }
1300  }  }
1301    
1302  int Finley_Quad_getNumNodesTri(int order) {  int Quad_getNumNodesTri(int order) {
1303    int numQuadNodesLine;    int numQuadNodesLine;
1304    if (order<=1) {    if (order<=1) {
1305        return 1;        return 1;
# Line 1320  int Finley_Quad_getNumNodesTri(int order Line 1320  int Finley_Quad_getNumNodesTri(int order
1320    } else if (order<=9){    } else if (order<=9){
1321       return 19;       return 19;
1322    } else {    } else {
1323        numQuadNodesLine=Finley_Quad_getNumNodesLine(order+1);        numQuadNodesLine=Quad_getNumNodesLine(order+1);
1324        if (Finley_noError()) {        if (Finley_noError()) {
1325           return numQuadNodesLine*numQuadNodesLine;           return numQuadNodesLine*numQuadNodesLine;
1326        } else {        } else {
# Line 1329  int Finley_Quad_getNumNodesTri(int order Line 1329  int Finley_Quad_getNumNodesTri(int order
1329    }    }
1330  }  }
1331    
1332  int Finley_Quad_getNumNodesRec(int order) {  int Quad_getNumNodesRec(int order) {
1333    int numQuadNodesLine;    int numQuadNodesLine;
1334    numQuadNodesLine=Finley_Quad_getNumNodesLine(order);    numQuadNodesLine=Quad_getNumNodesLine(order);
1335    if (Finley_noError()) {    if (Finley_noError()) {
1336        return numQuadNodesLine*numQuadNodesLine;        return numQuadNodesLine*numQuadNodesLine;
1337    } else {    } else {
# Line 1339  int Finley_Quad_getNumNodesRec(int order Line 1339  int Finley_Quad_getNumNodesRec(int order
1339    }    }
1340  }  }
1341    
1342  int Finley_Quad_getNumNodesTet(int order) {  int Quad_getNumNodesTet(int order) {
1343    int numQuadNodesLine;    int numQuadNodesLine;
1344    if (order<=1) {    if (order<=1) {
1345        return 1;        return 1;
# Line 1358  int Finley_Quad_getNumNodesTet(int order Line 1358  int Finley_Quad_getNumNodesTet(int order
1358    } else if (order<=8){    } else if (order<=8){
1359        return 45;        return 45;
1360    } else {    } else {
1361       numQuadNodesLine=Finley_Quad_getNumNodesLine(order+2);       numQuadNodesLine=Quad_getNumNodesLine(order+2);
1362       if (Finley_noError()) {       if (Finley_noError()) {
1363           return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;           return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1364       } else {       } else {
# Line 1367  int Finley_Quad_getNumNodesTet(int order Line 1367  int Finley_Quad_getNumNodesTet(int order
1367    }    }
1368  }  }
1369    
1370  int Finley_Quad_getNumNodesHex(int order) {  int Quad_getNumNodesHex(int order) {
1371    int numQuadNodesLine;    int numQuadNodesLine;
1372    numQuadNodesLine=Finley_Quad_getNumNodesLine(order);    numQuadNodesLine=Quad_getNumNodesLine(order);
1373    if (Finley_noError()) {    if (Finley_noError()) {
1374        return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;        return numQuadNodesLine*numQuadNodesLine*numQuadNodesLine;
1375    } else {    } else {
1376        return -1;        return -1;
1377    }    }
1378  }  }
1379  dim_t Finley_Quad_MacroPoint(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroPoint(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1380                              dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )                              dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1381  {  {
1382      return 0;      return 0;
1383            
1384  }  }
1385  dim_t Finley_Quad_MacroLine(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroLine(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1386                              dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )                              dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )
1387  {  {
1388      #define DIM 1      #define DIM 1
1389      dim_t s,q,i;      dim_t s,q,i;
1390      register double x0, w;      register double x0, w;
1391      const double f=1./((double)numSubElements);      const double f=1./((double)numSubElements);
1392            
1393      if (new_len < numSubElements*numQuadNodes) {      if (new_len < numSubElements*numQuadNodes) {
1394          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroLine: array for new quadrature scheme is too small");          Finley_setError(MEMORY_ERROR,"Quad_MacroLine: array for new quadrature scheme is too small");
1395      }      }
1396      for (q=0; q<numQuadNodes; ++q) {      for (q=0; q<numQuadNodes; ++q) {
1397                            
1398          x0=quadNodes[INDEX2(0,q,DIM)];          x0=quadNodes[INDEX2(0,q,DIM)];
1399          w=f*quadWeights[q];          w=f*quadWeights[q];
1400                    
1401          for (s=0; s<numSubElements; ++s) {          for (s=0; s<numSubElements; ++s) {
1402             new_quadWeights[INDEX2(q,s, numQuadNodes)]   =w;             new_quadWeights[INDEX2(q,s, numQuadNodes)]   =w;
1403                 new_quadNodes[INDEX3(0,q,s, DIM,numQuadNodes)]   =(x0+s)*f;                 new_quadNodes[INDEX3(0,q,s, DIM,numQuadNodes)]   =(x0+s)*f;
1404                     for (i=0;i<numF;i++) new_dFdv[INDEX4(i,0,q,s, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,numF, q,DIM)]*f;                     for (i=0;i<numF;i++) new_dFdv[INDEX4(i,0,q,s, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,numF, q,DIM)]*f;
1405          }          }
1406    
1407      }      }
1408      #undef DIM      #undef DIM
1409      return numSubElements*numQuadNodes;      return numSubElements*numQuadNodes;
1410  }  }
1411  #define HALF 0.5  #define HALF 0.5
1412  #define TWO 2.  #define TWO 2.
1413  dim_t Finley_Quad_MacroTri(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroTri(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1414                                      dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )                                      dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )
1415  {  {
1416    
1417      #define DIM 2      #define DIM 2
1418      dim_t q,i;      dim_t q,i;
1419      register double x0, x1, w, df0, df1;      register double x0, x1, w, df0, df1;
1420            
1421      if (new_len < numSubElements*numQuadNodes) {      if (new_len < numSubElements*numQuadNodes) {
1422          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: array for new quadrature scheme is too small");          Finley_setError(MEMORY_ERROR,"Quad_MacroTri: array for new quadrature scheme is too small");
1423          return -1;          return -1;
1424      }      }
1425      if (numSubElements==1) {      if (numSubElements==1) {
1426                    
1427          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1428                            
1429              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1430              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1431              w=quadWeights[q];              w=quadWeights[q];
1432                            
1433              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1434                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1435                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1436                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1437                              new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];                              new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1438                              new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];                              new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1439                          }                          }
1440    
1441          }          }
1442                    
1443      } else if (numSubElements==4) {      } else if (numSubElements==4) {
1444                  const double f = 0.25;                  const double f = 0.25;
1445          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1446                            
1447              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1448              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1449              w=f*quadWeights[q];              w=f*quadWeights[q];
1450                            
1451              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1452                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1453                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1454                            
1455              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1456                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;
1457                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);
1458                            
1459              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1460                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);
1461                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;
1462                            
1463              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1464                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(1-x0);                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(1-x0);
1465                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(1-x1);                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(1-x1);
1466    
1467                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1468                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
# Line 1482  dim_t Finley_Quad_MacroTri(dim_t numSubE Line 1482  dim_t Finley_Quad_MacroTri(dim_t numSubE
1482                          }                          }
1483    
1484    
1485          }          }
1486      } else {      } else {
1487          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: unable to create quadrature scheme for macro element.");          Finley_setError(MEMORY_ERROR,"Quad_MacroTri: unable to create quadrature scheme for macro element.");
1488          return -1;          return -1;
1489      }      }
1490      #undef DIM      #undef DIM
1491      return numSubElements*numQuadNodes;      return numSubElements*numQuadNodes;
1492  }  }
1493    
1494  dim_t Finley_Quad_MacroRec(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroRec(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1495                          dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv   )                          dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv   )
1496  {  {
1497    
1498      #define DIM 2      #define DIM 2
1499      dim_t q,i;      dim_t q,i;
1500      register double x0, x1, w, df0, df1;      register double x0, x1, w, df0, df1;
1501            
1502      if (new_len < numSubElements*numQuadNodes) {      if (new_len < numSubElements*numQuadNodes) {
1503          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: array for new quadrature scheme is too small");          Finley_setError(MEMORY_ERROR,"Quad_MacroRec: array for new quadrature scheme is too small");
1504          return -1;          return -1;
1505      }      }
1506      if (numSubElements==1) {      if (numSubElements==1) {
1507                    
1508          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1509                            
1510              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1511              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1512              w=quadWeights[q];              w=quadWeights[q];
1513                            
1514              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1515                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1516                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1517                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1518                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0, q,numF, DIM)];                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0, q,numF, DIM)];
1519                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1, q,numF, DIM)];                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1, q,numF, DIM)];
1520                          }                          }
1521    
1522          }          }
1523                    
1524      } else if (numSubElements==4) {      } else if (numSubElements==4) {
1525                  const double f = 0.25;                  const double f = 0.25;
1526          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1527                            
1528              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1529              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1530              w=f*quadWeights[q];              w=f*quadWeights[q];
1531                            
1532              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1533                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1534                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1535                            
1536              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1537                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;
1538                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);
1539                            
1540              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1541                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);
1542                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;
1543                            
1544              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1545                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);
1546                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);
1547    
1548                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1549                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
# Line 1562  dim_t Finley_Quad_MacroRec(dim_t numSubE Line 1562  dim_t Finley_Quad_MacroRec(dim_t numSubE
1562                              new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;                              new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1563                          }                          }
1564    
1565          }          }
1566      } else {      } else {
1567          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: unable to create quadrature scheme for macro element.");          Finley_setError(MEMORY_ERROR,"Quad_MacroRec: unable to create quadrature scheme for macro element.");
1568          return -1;          return -1;
1569      }      }
1570      #undef DIM      #undef DIM
1571      return numSubElements*numQuadNodes;      return numSubElements*numQuadNodes;
1572  }  }
1573    
1574    
1575  dim_t Finley_Quad_MacroTet(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroTet(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1576                                  dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )                                  dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1577  {  {
1578      #define DIM 3      #define DIM 3
1579      dim_t q, i;      dim_t q, i;
1580      register double x0, x1, x2, w, df0, df1, df2;      register double x0, x1, x2, w, df0, df1, df2;
1581            
1582      if (new_len < numSubElements*numQuadNodes) {      if (new_len < numSubElements*numQuadNodes) {
1583          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: array for new quadrature scheme is too small");          Finley_setError(MEMORY_ERROR,"Quad_MacroTet: array for new quadrature scheme is too small");
1584          return -1;          return -1;
1585      }      }
1586      if (numSubElements==1) {      if (numSubElements==1) {
1587                    
1588          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1589                            
1590              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1591              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1592              x2=quadNodes[INDEX2(2,q,DIM)];              x2=quadNodes[INDEX2(2,q,DIM)];
1593              w=quadWeights[q];              w=quadWeights[q];
1594                            
1595              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1596                      new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;                      new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1597                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1598              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;
1599    
1600                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1601                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1602                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1603                              new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q, numF, DIM)];                              new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q, numF, DIM)];
1604                          }                          }
1605          }          }
1606                    
1607      } else if (numSubElements==8) {      } else if (numSubElements==8) {
1608                  const double f = 0.125;                  const double f = 0.125;
1609          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1610                            
1611              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1612              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1613              x2=quadNodes[INDEX2(2,q,DIM)];              x2=quadNodes[INDEX2(2,q,DIM)];
1614              w=f*quadWeights[q];              w=f*quadWeights[q];
1615                            
1616              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1617                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1618                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1619              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;
1620                            
1621              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1622                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);
1623                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;
1624              new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;
1625                            
1626              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1627                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;
1628                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);
1629              new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;
1630                            
1631              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1632                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*x0;
1633                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*x1;
1634              new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*(x2+1);              new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*(x2+1);
1635    
1636              new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;
1637                  new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*(1-x1);                  new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*(1-x1);
1638                  new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*(x0+x1);                  new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*(x0+x1);
1639              new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*x2;
1640                            
1641              new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;
1642                  new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(1-x0-x2);                  new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(1-x0-x2);
1643                  new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*(1-x1);                  new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*(1-x1);
1644              new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x0+x1);              new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x0+x1);
1645                            
1646              new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;
1647                  new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x2;                  new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x2;
1648                  new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(1-x0-x2);                  new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(1-x0-x2);
1649              new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(1-x1);              new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(1-x1);
1650                            
1651              new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;
1652                  new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+x2);                  new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+x2);
1653                  new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*x1;
1654              new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(1-x0-x1);              new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(1-x0-x1);
1655    
1656                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1657                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
# Line 1691  dim_t Finley_Quad_MacroTet(dim_t numSubE Line 1691  dim_t Finley_Quad_MacroTet(dim_t numSubE
1691                              new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = -df0+df2;                              new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = -df0+df2;
1692                          }                          }
1693    
1694          }          }
1695      } else {      } else {
1696          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: unable to create quadrature scheme for macro element.");          Finley_setError(MEMORY_ERROR,"Quad_MacroTet: unable to create quadrature scheme for macro element.");
1697          return -1;          return -1;
1698      }      }
1699      #undef DIM      #undef DIM
1700      return numSubElements*numQuadNodes;      return numSubElements*numQuadNodes;
1701  }  }
1702    
1703    
1704  dim_t Finley_Quad_MacroHex(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,  dim_t Quad_MacroHex(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1705                     dim_t new_len, double* new_quadNodes, double* new_quadWeights , double* new_dFdv)                     dim_t new_len, double* new_quadNodes, double* new_quadWeights , double* new_dFdv)
1706  {  {
1707    
1708      #define DIM 3      #define DIM 3
1709      dim_t q, i;      dim_t q, i;
1710      register double x0, x1, x2, w, df0, df1, df2;      register double x0, x1, x2, w, df0, df1, df2;
1711            
1712      if (new_len < numSubElements*numQuadNodes) {      if (new_len < numSubElements*numQuadNodes) {
1713          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: array for new quadrature scheme is too small");          Finley_setError(MEMORY_ERROR,"Quad_MacroHex: array for new quadrature scheme is too small");
1714          return -1;          return -1;
1715      }      }
1716      if (numSubElements==1) {      if (numSubElements==1) {
1717                    
1718          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1719                            
1720              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1721              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1722              x2=quadNodes[INDEX2(2,q,DIM)];              x2=quadNodes[INDEX2(2,q,DIM)];
1723              w=quadWeights[q];              w=quadWeights[q];
1724                            
1725              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1726                      new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;                      new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1727                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1728              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;
1729    
1730                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1731                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q,numF, DIM)];                              new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q,numF, DIM)];
1732                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q,numF, DIM)];                              new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q,numF, DIM)];
1733                              new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q,numF, DIM)];                              new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q,numF, DIM)];
1734                          }                          }
1735          }          }
1736                    
1737      } else if (numSubElements==8) {      } else if (numSubElements==8) {
1738                  const double f = 0.125;                  const double f = 0.125;
1739          for (q=0; q<numQuadNodes; ++q) {          for (q=0; q<numQuadNodes; ++q) {
1740                            
1741              x0=quadNodes[INDEX2(0,q,DIM)];              x0=quadNodes[INDEX2(0,q,DIM)];
1742              x1=quadNodes[INDEX2(1,q,DIM)];              x1=quadNodes[INDEX2(1,q,DIM)];
1743              x2=quadNodes[INDEX2(2,q,DIM)];              x2=quadNodes[INDEX2(2,q,DIM)];
1744              w=f*quadWeights[q];              w=f*quadWeights[q];
1745                            
1746              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1747                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1748                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1749              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;
1750                            
1751              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1752                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);
1753                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;
1754              new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;
1755                            
1756              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1757                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;
1758                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);
1759              new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;
1760                            
1761              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1762                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);
1763                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);
1764              new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*x2;              new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*x2;
1765                            
1766              new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;
1767                  new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*x0;
1768                  new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*x1;
1769              new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*(x2+1);              new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*(x2+1);
1770                            
1771              new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;
1772                  new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(x0+1);                  new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(x0+1);
1773                  new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*x1;                  new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*x1;
1774              new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x2+1);              new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x2+1);
1775                            
1776              new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;
1777                  new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x0;                  new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x0;
1778                  new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(x1+1);
1779              new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(x2+1);              new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(x2+1);
1780                            
1781              new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;              new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;
1782              new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+1);              new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+1);
1783                  new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*(x1+1);                  new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*(x1+1);
1784              new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(x2+1);              new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(x2+1);
1785    
1786                          for (i=0;i<numF;i++) {                          for (i=0;i<numF;i++) {
1787                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;                              df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
# Line 1821  dim_t Finley_Quad_MacroHex(dim_t numSubE Line 1821  dim_t Finley_Quad_MacroHex(dim_t numSubE
1821                              new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = df2;                              new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = df2;
1822                          }                          }
1823    
1824          }          }
1825      } else {      } else {
1826          Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: unable to create quadrature scheme for macro element.");          Finley_setError(MEMORY_ERROR,"Quad_MacroHex: unable to create quadrature scheme for macro element.");
1827          return -1;          return -1;
1828      }      }
1829      #undef DIM      #undef DIM
1830      return numSubElements*numQuadNodes;      return numSubElements*numQuadNodes;
1831  }  }
1832    
1833    } // namespace finley
1834    

Legend:
Removed from v.4491  
changed lines
  Added in v.4492

  ViewVC Help
Powered by ViewVC 1.1.26