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
20
21  /************************************************************************************/  *****************************************************************************/
22
24  #include "esysUtils/index.h"  #include "esysUtils/index.h"
# Line 28  Line 28
30
31  /************************************************************************************/  namespace finley {
32
41  };  };
42
44  {  {
45      int ptr=0;      int ptr=0;
49         ptr++;         ptr++;
50      }      }
51      if (out==NULL) {      if (out==NULL) {
53      }      }
54      return out;      return out;
55  }  }
60  /*   as a squeezed scheme on a quad [0,1]^2 */  /*   as a squeezed scheme on a quad [0,1]^2 */
61
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
350    } else {    } else {
351
352      /*  get scheme on [0.1]^2 */      /*  get scheme on [0.1]^2 */
354      if (! Finley_noError()) return;      if (! Finley_noError()) return;
355
356      /*  squeeze it: */      /*  squeeze it: */
374  /*   as a squeezed scheme on a hex [0,1]^3 */  /*   as a squeezed scheme on a hex [0,1]^3 */
375
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;
969
970      /*  get scheme on [0.1]^3 */      /*  get scheme on [0.1]^3 */
971
973      if (! Finley_noError()) return;      if (! Finley_noError()) return;
974
975      /*  squeeze it: */      /*  squeeze it: */
1003  /*   as a X-product of a 1D scheme. */  /*   as a X-product of a 1D scheme. */
1004
1006    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1019
1020           /*  get 1D scheme: */           /*  get 1D scheme: */
1021
1023
1024           /*  make 2D scheme: */           /*  make 2D scheme: */
1025
1039         }         }
1040       }       }
1041       if (!set) {       if (!set) {
1043           Finley_setError(VALUE_ERROR,error_msg);           Finley_setError(VALUE_ERROR,error_msg);
1044       }       }
1054  /*   as a X-product of a 1D scheme. */  /*   as a X-product of a 1D scheme. */
1055
1057    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1070
1071           /*  get 1D scheme: */           /*  get 1D scheme: */
1072
1074
1075           /*  make 3D scheme: */           /*  make 3D scheme: */
1076
1093         }         }
1094       }       }
1095       if (!set) {       if (!set) {
1097            Finley_setError(VALUE_ERROR,error_msg);            Finley_setError(VALUE_ERROR,error_msg);
1098       }       }
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
1113          return;          return;
1114    } else {    } else {
1116    }    }
1117  }  }
1118
1122  /*   The nodes and weights are set from a table. */  /*   The nodes and weights are set from a table. */
1123
1126        case 1:        case 1:
1264          break;          break;
1265
1266        default:        default:
1268          break;          break;
1269    }    }
1270  }  }
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
1281      return 1;      return 1;
1282  }  }
1283
1285    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1286    if (order <0 ) {    if (order <0 ) {
1288      return -1;      return -1;
1289    } else {    } else {
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).",
1293        Finley_setError(VALUE_ERROR,error_msg);        Finley_setError(VALUE_ERROR,error_msg);
1294        return -1;        return -1;
1299    }    }
1300  }  }
1301
1304    if (order<=1) {    if (order<=1) {
1305        return 1;        return 1;
1320    } else if (order<=9){    } else if (order<=9){
1321       return 19;       return 19;
1322    } else {    } else {
1324        if (Finley_noError()) {        if (Finley_noError()) {
1326        } else {        } else {
1329    }    }
1330  }  }
1331
1335    if (Finley_noError()) {    if (Finley_noError()) {
1337    } else {    } else {
1339    }    }
1340  }  }
1341
1344    if (order<=1) {    if (order<=1) {
1345        return 1;        return 1;
1358    } else if (order<=8){    } else if (order<=8){
1359        return 45;        return 45;
1360    } else {    } else {
1362       if (Finley_noError()) {       if (Finley_noError()) {
1364       } else {       } else {
1367    }    }
1368  }  }
1369
1373    if (Finley_noError()) {    if (Finley_noError()) {
1375    } else {    } else {
1376        return -1;        return -1;
1377    }    }
1378  }  }
1381  {  {
1382      return 0;      return 0;
1383
1384  }  }
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
1395      }      }
1397
1400
1401          for (s=0; s<numSubElements; ++s) {          for (s=0; s<numSubElements; ++s) {
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
1410  }  }
1411  #define HALF 0.5  #define HALF 0.5
1412  #define TWO 2.  #define TWO 2.
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
1423          return -1;          return -1;
1424      }      }
1425      if (numSubElements==1) {      if (numSubElements==1) {
1426
1428
1432
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;
1446
1450
1454
1458
1462
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;
1482                          }                          }
1483
1484
1485          }          }
1486      } else {      } else {
1488          return -1;          return -1;
1489      }      }
1490      #undef DIM      #undef DIM
1492  }  }
1493
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
1504          return -1;          return -1;
1505      }      }
1506      if (numSubElements==1) {      if (numSubElements==1) {
1507
1509
1513
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;
1527
1531
1535
1539
1543
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;
1563                          }                          }
1564
1565          }          }
1566      } else {      } else {
1568          return -1;          return -1;
1569      }      }
1570      #undef DIM      #undef DIM
1572  }  }
1573
1574
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
1584          return -1;          return -1;
1585      }      }
1586      if (numSubElements==1) {      if (numSubElements==1) {
1587
1589
1594
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;
1610
1615
1620
1625
1630
1635
1640
1645
1650
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;
1692                          }                          }
1693
1694          }          }
1695      } else {      } else {
1697          return -1;          return -1;
1698      }      }
1699      #undef DIM      #undef DIM
1701  }  }
1702
1703
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
1714          return -1;          return -1;
1715      }      }
1716      if (numSubElements==1) {      if (numSubElements==1) {
1717
1719
1724
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;
1740
1745
1750
1755
1760
1765
1770
1775
1780
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;
1822                          }                          }
1823
1824          }          }
1825      } else {      } else {