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

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

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

revision 1342 by gross, Thu Nov 8 23:56:58 2007 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Finley:  */  /*   Finley: quadrature schemes */
18    
19  /**************************************************************/  /**************************************************************/
20    
21  #include "Quadrature.h"  #include "Quadrature.h"
22    
23    
24  #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]  #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
25  #define QUADWEIGHTS(_I_) quadWeights[_I_]  #define QUADWEIGHTS(_I_) quadWeights[_I_]
26    
27  /**************************************************************/  /**************************************************************/
28    
29    Finley_QuadInfo Finley_QuadInfoList[]={
30        {PointQuad, "Point", 0,  1,     Finley_Quad_getNodesPoint,      Finley_Quad_getNumNodesPoint,   Finley_Quad_MacroPoint} ,
31        {LineQuad,  "Line",  1,  2,     Finley_Quad_getNodesLine,       Finley_Quad_getNumNodesLine,    Finley_Quad_MacroLine} ,
32        {TriQuad,   "Tri",   2,  3,     Finley_Quad_getNodesTri,        Finley_Quad_getNumNodesTri,     Finley_Quad_MacroTri},
33        {RecQuad,   "Rec",   2,  4,     Finley_Quad_getNodesRec,        Finley_Quad_getNumNodesRec,     Finley_Quad_MacroRec},
34        {TetQuad,   "Tet",   3,  4,     Finley_Quad_getNodesTet,        Finley_Quad_getNumNodesTet,     Finley_Quad_MacroTet},
35        {HexQuad,   "Hex",   3,  8,     Finley_Quad_getNodesHex,        Finley_Quad_getNumNodesHex,     Finley_Quad_MacroHex},
36        {NoQuad, "NoType", 0,  1,   Finley_Quad_getNodesPoint,      Finley_Quad_getNumNodesPoint,   Finley_Quad_MacroPoint}
37    };
38    
39    Finley_QuadInfo* Finley_QuadInfo_getInfo(Finley_QuadTypeId id)
40    {
41        int ptr=0;
42        Finley_QuadInfo* out=NULL;
43        while (Finley_QuadInfoList[ptr].TypeId!=NoQuad && out==NULL) {
44           if (Finley_QuadInfoList[ptr].TypeId==id) out=&(Finley_QuadInfoList[ptr]);
45           ptr++;
46        }
47        if (out==NULL) {
48            Finley_setError(VALUE_ERROR,"Finley_QuadInfo_getInfo: canot find requested quadrature scheme.");
49        }
50        return out;
51    }
52    
53    /**************************************************************/
54    
55  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tri  */  /*   get a quadrature scheme with numQuadNodes quadrature nodes for the tri  */
56  /*   as a queezed scheme on a quad [0,1]^2 */  /*   as a queezed scheme on a quad [0,1]^2 */
57    
# Line 1027  void Finley_Quad_getNodesHex(int numQuad Line 1053  void Finley_Quad_getNodesHex(int numQuad
1053    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
1054    int numQuadNodes1d,i,j,k,l;    int numQuadNodes1d,i,j,k,l;
1055    double *quadNodes1d=NULL,*quadWeights1d=NULL;    double *quadNodes1d=NULL,*quadWeights1d=NULL;
1056    bool_t set;    bool_t set=FALSE;
1057    #define DIM 3    #define DIM 3
1058        
1059    /*  find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */    /*  find numQuadNodes1d with numQuadNodes1d**3==numQuadNodes: */
# Line 1079  void Finley_Quad_getNodesHex(int numQuad Line 1105  void Finley_Quad_getNodesHex(int numQuad
1105  /*   an error. */  /*   an error. */
1106    
1107  void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {  void Finley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
1108    if (numQuadNodes!=0) {    if (numQuadNodes==0) {
1109         Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: There is no quadrature scheme on points.");          return;
1110      } else {
1111           Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesPoint: Illegal number of quadrature nodes.");
1112    }    }
1113  }  }
1114    
# Line 1232  void Finley_Quad_getNodesLine(int numQua Line 1260  void Finley_Quad_getNodesLine(int numQua
1260          break;          break;
1261    
1262        default:        default:
1263          Finley_setError(VALUE_ERROR,"__FILE__: Negative intergration order.");          Finley_setError(VALUE_ERROR,"Finley_Quad_getNodesLine: Negative intergration order.");
1264          break;          break;
1265    }    }
1266  }  }
 /**************************************************************/  
 /*                                                            */  
 /*  the following function are used define the meshes on the surface in the xy-plane */  
1267    
 /* 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=NULL;  
     #define DIM dim  
     quadNodesOnFace=TMPMEMALLOC(numQuadNodes*(dim-1),double);  
   
     if (! Finley_checkPtr(quadNodesOnFace) ) {  
        getFaceNodes(numQuadNodes,quadNodesOnFace,quadWeights);  
        if (Finley_noError()) {  
           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;  
           }  
        }  
        TMPMEMFREE(quadNodesOnFace);  
     }  
     #undef DIM  
 }  
1268    
1269  /**************************************************************/  /**************************************************************/
1270    
# Line 1284  void Finley_Quad_makeNodesOnFace(int dim Line 1274  void Finley_Quad_makeNodesOnFace(int dim
1274    
1275    
1276  int Finley_Quad_getNumNodesPoint(int order) {  int Finley_Quad_getNumNodesPoint(int order) {
   if (order <0 ) {  
     Finley_setError(VALUE_ERROR,"Finley_Quad_getNumNodesPoint: Negative intergration order.");  
     return -1;  
   } else {  
1277      return 0;      return 0;
   }  
1278  }  }
1279    
1280  int Finley_Quad_getNumNodesLine(int order) {  int Finley_Quad_getNumNodesLine(int order) {
# Line 1387  int Finley_Quad_getNumNodesHex(int order Line 1372  int Finley_Quad_getNumNodesHex(int order
1372        return -1;        return -1;
1373    }    }
1374  }  }
1375    dim_t Finley_Quad_MacroPoint(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1376                                dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1377    {
1378        return 0;
1379        
1380    }
1381    dim_t Finley_Quad_MacroLine(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1382                                dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )
1383    {
1384        #define DIM 1
1385        dim_t s,q,i;
1386        register double x0, w;
1387        const double f=1./((double)numSubElements);
1388        
1389        if (new_len < numSubElements*numQuadNodes) {
1390            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroLine: array for new qurature scheme is too small");
1391        }
1392        for (q=0; q<numQuadNodes; ++q) {
1393                
1394            x0=quadNodes[INDEX2(0,q,DIM)];
1395            w=f*quadWeights[q];
1396            
1397            for (s=0; s<numSubElements; ++s) {
1398               new_quadWeights[INDEX2(q,s, numQuadNodes)]   =w;
1399                   new_quadNodes[INDEX3(0,q,s, DIM,numQuadNodes)]   =(x0+s)*f;
1400                       for (i=0;i<numF;i++) new_dFdv[INDEX4(i,0,q,s, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,numF, q,DIM)]*f;
1401            }
1402    
1403        }
1404        #undef DIM
1405        return numSubElements*numQuadNodes;
1406    }
1407    #define HALF 0.5
1408    #define TWO 2.
1409    dim_t Finley_Quad_MacroTri(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1410                                        dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv  )
1411    {
1412    
1413        #define DIM 2
1414        dim_t q,i;
1415        register double x0, x1, w, df0, df1;
1416        
1417        if (new_len < numSubElements*numQuadNodes) {
1418            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: array for new qurature scheme is too small");
1419            return -1;
1420        }
1421        if (numSubElements==1) {
1422            
1423            for (q=0; q<numQuadNodes; ++q) {
1424                
1425                x0=quadNodes[INDEX2(0,q,DIM)];
1426                x1=quadNodes[INDEX2(1,q,DIM)];
1427                w=quadWeights[q];
1428                
1429                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1430                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1431                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1432                            for (i=0;i<numF;i++) {
1433                                new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1434                                new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1435                            }
1436    
1437            }
1438            
1439        } else if (numSubElements==4) {
1440                    const double f = 0.25;
1441            for (q=0; q<numQuadNodes; ++q) {
1442                
1443                x0=quadNodes[INDEX2(0,q,DIM)];
1444                x1=quadNodes[INDEX2(1,q,DIM)];
1445                w=f*quadWeights[q];
1446                
1447                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1448                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1449                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1450                
1451                new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1452                    new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;
1453                    new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);
1454                
1455                new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1456                    new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);
1457                    new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;
1458                
1459                new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1460                    new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(1-x0);
1461                    new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(1-x1);
1462    
1463                            for (i=0;i<numF;i++) {
1464                                df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1465                                df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1466    
1467                                new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1468                                new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1469    
1470                                new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1471                                new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1472    
1473                                new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1474                                new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1475    
1476                                new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = -df0;
1477                                new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = -df1;
1478                            }
1479    
1480    
1481            }
1482        } else {
1483            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTri: unable to create quadrature scheme for macro element.");
1484            return -1;
1485        }
1486        #undef DIM
1487        return numSubElements*numQuadNodes;
1488    }
1489    
1490    dim_t Finley_Quad_MacroRec(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1491                            dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv   )
1492    {
1493    
1494        #define DIM 2
1495        dim_t q,i;
1496        register double x0, x1, w, df0, df1;
1497        
1498        if (new_len < numSubElements*numQuadNodes) {
1499            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: array for new qurature scheme is too small");
1500            return -1;
1501        }
1502        if (numSubElements==1) {
1503            
1504            for (q=0; q<numQuadNodes; ++q) {
1505                
1506                x0=quadNodes[INDEX2(0,q,DIM)];
1507                x1=quadNodes[INDEX2(1,q,DIM)];
1508                w=quadWeights[q];
1509                
1510                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1511                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1512                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1513                            for (i=0;i<numF;i++) {
1514                                new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0, q,numF, DIM)];
1515                                new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1, q,numF, DIM)];
1516                            }
1517    
1518            }
1519            
1520        } else if (numSubElements==4) {
1521                    const double f = 0.25;
1522            for (q=0; q<numQuadNodes; ++q) {
1523                
1524                x0=quadNodes[INDEX2(0,q,DIM)];
1525                x1=quadNodes[INDEX2(1,q,DIM)];
1526                w=f*quadWeights[q];
1527                
1528                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1529                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1530                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1531                
1532                new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1533                    new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*x0;
1534                    new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*(x1+1);
1535                
1536                new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1537                    new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*(x0+1);
1538                    new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*x1;
1539                
1540                new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1541                    new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);
1542                    new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);
1543    
1544                            for (i=0;i<numF;i++) {
1545                                df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1546                                df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1547    
1548                                new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1549                                new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1550    
1551                                new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1552                                new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1553    
1554                                new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1555                                new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1556    
1557                                new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1558                                new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1559                            }
1560    
1561            }
1562        } else {
1563            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroRec: unable to create quadrature scheme for macro element.");
1564            return -1;
1565        }
1566        #undef DIM
1567        return numSubElements*numQuadNodes;
1568    }
1569    
1570    
1571    dim_t Finley_Quad_MacroTet(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1572                                    dim_t new_len, double* new_quadNodes, double* new_quadWeights, double* new_dFdv )
1573    {
1574        #define DIM 3
1575        dim_t q, i;
1576        register double x0, x1, x2, w, df0, df1, df2;
1577        
1578        if (new_len < numSubElements*numQuadNodes) {
1579            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: array for new qurature scheme is too small");
1580            return -1;
1581        }
1582        if (numSubElements==1) {
1583            
1584            for (q=0; q<numQuadNodes; ++q) {
1585                
1586                x0=quadNodes[INDEX2(0,q,DIM)];
1587                x1=quadNodes[INDEX2(1,q,DIM)];
1588                x2=quadNodes[INDEX2(2,q,DIM)];
1589                w=quadWeights[q];
1590                
1591                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1592                        new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1593                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1594                new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;
1595    
1596                            for (i=0;i<numF;i++) {
1597                                new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q, numF, DIM)];
1598                                new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q, numF, DIM)];
1599                                new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q, numF, DIM)];
1600                            }
1601            }
1602            
1603        } else if (numSubElements==8) {
1604                    const double f = 0.125;
1605            for (q=0; q<numQuadNodes; ++q) {
1606                
1607                x0=quadNodes[INDEX2(0,q,DIM)];
1608                x1=quadNodes[INDEX2(1,q,DIM)];
1609                x2=quadNodes[INDEX2(2,q,DIM)];
1610                w=f*quadWeights[q];
1611                
1612                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1613                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1614                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1615                new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;
1616                
1617                new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1618                    new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);
1619                    new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;
1620                new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;
1621                
1622                new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1623                    new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;
1624                    new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);
1625                new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;
1626                
1627                new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1628                    new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*x0;
1629                    new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*x1;
1630                new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*(x2+1);
1631    
1632                new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;
1633                    new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*(1-x1);
1634                    new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*(x0+x1);
1635                new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*x2;
1636                
1637                new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;
1638                    new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(1-x0-x2);
1639                    new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*(1-x1);
1640                new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x0+x1);
1641                
1642                new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;
1643                    new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x2;
1644                    new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(1-x0-x2);
1645                new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(1-x1);
1646                
1647                new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;
1648                    new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+x2);
1649                    new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*x1;
1650                new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(1-x0-x1);
1651    
1652                            for (i=0;i<numF;i++) {
1653                                df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1654                                df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1655                                df2=dFdv[INDEX3(i,2,q, numF, DIM)]*TWO;
1656    
1657                                new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1658                                new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1659                                new_dFdv[INDEX4(i,2,q,0, numF,DIM,numQuadNodes)] = df2;
1660    
1661                                new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1662                                new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1663                                new_dFdv[INDEX4(i,2,q,1, numF,DIM,numQuadNodes)] = df2;
1664    
1665                                new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1666                                new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1667                                new_dFdv[INDEX4(i,2,q,2, numF,DIM,numQuadNodes)] = df2;
1668    
1669                                new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1670                                new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1671                                new_dFdv[INDEX4(i,2,q,3, numF,DIM,numQuadNodes)] = df2;
1672    
1673                                new_dFdv[INDEX4(i,0,q,4, numF,DIM,numQuadNodes)] = df0-df1;
1674                                new_dFdv[INDEX4(i,1,q,4, numF,DIM,numQuadNodes)] = df0;
1675                                new_dFdv[INDEX4(i,2,q,4, numF,DIM,numQuadNodes)] = df2;
1676    
1677                                new_dFdv[INDEX4(i,0,q,5, numF,DIM,numQuadNodes)] = -df2;
1678                                new_dFdv[INDEX4(i,1,q,5, numF,DIM,numQuadNodes)] = df0-df2-df1;
1679                                new_dFdv[INDEX4(i,2,q,5, numF,DIM,numQuadNodes)] = df0-df2;
1680    
1681                                new_dFdv[INDEX4(i,0,q,6, numF,DIM,numQuadNodes)] = -df0+df2;
1682                                new_dFdv[INDEX4(i,1,q,6, numF,DIM,numQuadNodes)] = -df0;
1683                                new_dFdv[INDEX4(i,2,q,6, numF,DIM,numQuadNodes)] = -df1;
1684    
1685                                new_dFdv[INDEX4(i,0,q,7, numF,DIM,numQuadNodes)] = df2;
1686                                new_dFdv[INDEX4(i,1,q,7, numF,DIM,numQuadNodes)] = -df0+df1+df2;
1687                                new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = -df0+df2;
1688                            }
1689    
1690            }
1691        } else {
1692            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroTet: unable to create quadrature scheme for macro element.");
1693            return -1;
1694        }
1695        #undef DIM
1696        return numSubElements*numQuadNodes;
1697    }
1698    
1699    
1700    dim_t Finley_Quad_MacroHex(dim_t numSubElements, int numQuadNodes, double* quadNodes, double* quadWeights, dim_t numF, double* dFdv,
1701                       dim_t new_len, double* new_quadNodes, double* new_quadWeights , double* new_dFdv)
1702    {
1703    
1704        #define DIM 3
1705        dim_t q, i;
1706        register double x0, x1, x2, w, df0, df1, df2;
1707        
1708        if (new_len < numSubElements*numQuadNodes) {
1709            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: array for new qurature scheme is too small");
1710            return -1;
1711        }
1712        if (numSubElements==1) {
1713            
1714            for (q=0; q<numQuadNodes; ++q) {
1715                
1716                x0=quadNodes[INDEX2(0,q,DIM)];
1717                x1=quadNodes[INDEX2(1,q,DIM)];
1718                x2=quadNodes[INDEX2(2,q,DIM)];
1719                w=quadWeights[q];
1720                
1721                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1722                        new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =x0;
1723                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =x1;
1724                new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =x2;
1725    
1726                            for (i=0;i<numF;i++) {
1727                                new_dFdv[INDEX4(i,0,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,0,q,numF, DIM)];
1728                                new_dFdv[INDEX4(i,1,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,1,q,numF, DIM)];
1729                                new_dFdv[INDEX4(i,2,q,0, numF, DIM,numQuadNodes)] = dFdv[INDEX3(i,2,q,numF, DIM)];
1730                            }
1731            }
1732            
1733        } else if (numSubElements==8) {
1734                    const double f = 0.125;
1735            for (q=0; q<numQuadNodes; ++q) {
1736                
1737                x0=quadNodes[INDEX2(0,q,DIM)];
1738                x1=quadNodes[INDEX2(1,q,DIM)];
1739                x2=quadNodes[INDEX2(2,q,DIM)];
1740                w=f*quadWeights[q];
1741                
1742                new_quadWeights[INDEX2(q,0,numQuadNodes)]       =w;
1743                    new_quadNodes[INDEX3(0,q,0,DIM,numQuadNodes)]   =HALF*x0;
1744                    new_quadNodes[INDEX3(1,q,0,DIM,numQuadNodes)]   =HALF*x1;
1745                new_quadNodes[INDEX3(2,q,0,DIM,numQuadNodes)]   =HALF*x2;
1746                
1747                new_quadWeights[INDEX2(q,1,numQuadNodes)]       =w;
1748                    new_quadNodes[INDEX3(0,q,1,DIM,numQuadNodes)]   =HALF*(x0+1);
1749                    new_quadNodes[INDEX3(1,q,1,DIM,numQuadNodes)]   =HALF*x1;
1750                new_quadNodes[INDEX3(2,q,1,DIM,numQuadNodes)]   =HALF*x2;
1751                
1752                new_quadWeights[INDEX2(q,2,numQuadNodes)]       =w;
1753                    new_quadNodes[INDEX3(0,q,2,DIM,numQuadNodes)]   =HALF*x0;
1754                    new_quadNodes[INDEX3(1,q,2,DIM,numQuadNodes)]   =HALF*(x1+1);
1755                new_quadNodes[INDEX3(2,q,2,DIM,numQuadNodes)]   =HALF*x2;
1756                
1757                new_quadWeights[INDEX2(q,3,numQuadNodes)]       =w;
1758                    new_quadNodes[INDEX3(0,q,3,DIM,numQuadNodes)]   =HALF*(x0+1);
1759                    new_quadNodes[INDEX3(1,q,3,DIM,numQuadNodes)]   =HALF*(x1+1);
1760                new_quadNodes[INDEX3(2,q,3,DIM,numQuadNodes)]   =HALF*x2;
1761                
1762                new_quadWeights[INDEX2(q,4,numQuadNodes)]       =w;
1763                    new_quadNodes[INDEX3(0,q,4,DIM,numQuadNodes)]   =HALF*x0;
1764                    new_quadNodes[INDEX3(1,q,4,DIM,numQuadNodes)]   =HALF*x1;
1765                new_quadNodes[INDEX3(2,q,4,DIM,numQuadNodes)]   =HALF*(x2+1);
1766                
1767                new_quadWeights[INDEX2(q,5,numQuadNodes)]       =w;
1768                    new_quadNodes[INDEX3(0,q,5,DIM,numQuadNodes)]   =HALF*(x0+1);
1769                    new_quadNodes[INDEX3(1,q,5,DIM,numQuadNodes)]   =HALF*x1;
1770                new_quadNodes[INDEX3(2,q,5,DIM,numQuadNodes)]   =HALF*(x2+1);
1771                
1772                new_quadWeights[INDEX2(q,6,numQuadNodes)]       =w;
1773                    new_quadNodes[INDEX3(0,q,6,DIM,numQuadNodes)]   =HALF*x0;
1774                    new_quadNodes[INDEX3(1,q,6,DIM,numQuadNodes)]   =HALF*(x1+1);
1775                new_quadNodes[INDEX3(2,q,6,DIM,numQuadNodes)]   =HALF*(x2+1);
1776                
1777                new_quadWeights[INDEX2(q,7,numQuadNodes)]       =w;
1778                new_quadNodes[INDEX3(0,q,7,DIM,numQuadNodes)]   =HALF*(x0+1);
1779                    new_quadNodes[INDEX3(1,q,7,DIM,numQuadNodes)]   =HALF*(x1+1);
1780                new_quadNodes[INDEX3(2,q,7,DIM,numQuadNodes)]   =HALF*(x2+1);
1781    
1782                            for (i=0;i<numF;i++) {
1783                                df0=dFdv[INDEX3(i,0,q, numF, DIM)]*TWO;
1784                                df1=dFdv[INDEX3(i,1,q, numF, DIM)]*TWO;
1785                                df2=dFdv[INDEX3(i,2,q, numF, DIM)]*TWO;
1786    
1787                                new_dFdv[INDEX4(i,0,q,0, numF,DIM,numQuadNodes)] = df0;
1788                                new_dFdv[INDEX4(i,1,q,0, numF,DIM,numQuadNodes)] = df1;
1789                                new_dFdv[INDEX4(i,2,q,0, numF,DIM,numQuadNodes)] = df2;
1790    
1791                                new_dFdv[INDEX4(i,0,q,1, numF,DIM,numQuadNodes)] = df0;
1792                                new_dFdv[INDEX4(i,1,q,1, numF,DIM,numQuadNodes)] = df1;
1793                                new_dFdv[INDEX4(i,2,q,1, numF,DIM,numQuadNodes)] = df2;
1794    
1795                                new_dFdv[INDEX4(i,0,q,2, numF,DIM,numQuadNodes)] = df0;
1796                                new_dFdv[INDEX4(i,1,q,2, numF,DIM,numQuadNodes)] = df1;
1797                                new_dFdv[INDEX4(i,2,q,2, numF,DIM,numQuadNodes)] = df2;
1798    
1799                                new_dFdv[INDEX4(i,0,q,3, numF,DIM,numQuadNodes)] = df0;
1800                                new_dFdv[INDEX4(i,1,q,3, numF,DIM,numQuadNodes)] = df1;
1801                                new_dFdv[INDEX4(i,2,q,3, numF,DIM,numQuadNodes)] = df2;
1802    
1803                                new_dFdv[INDEX4(i,0,q,4, numF,DIM,numQuadNodes)] = df0;
1804                                new_dFdv[INDEX4(i,1,q,4, numF,DIM,numQuadNodes)] = df1;
1805                                new_dFdv[INDEX4(i,2,q,4, numF,DIM,numQuadNodes)] = df2;
1806    
1807                                new_dFdv[INDEX4(i,0,q,5, numF,DIM,numQuadNodes)] = df0;
1808                                new_dFdv[INDEX4(i,1,q,5, numF,DIM,numQuadNodes)] = df1;
1809                                new_dFdv[INDEX4(i,2,q,5, numF,DIM,numQuadNodes)] = df2;
1810    
1811                                new_dFdv[INDEX4(i,0,q,6, numF,DIM,numQuadNodes)] = df0;
1812                                new_dFdv[INDEX4(i,1,q,6, numF,DIM,numQuadNodes)] = df1;
1813                                new_dFdv[INDEX4(i,2,q,6, numF,DIM,numQuadNodes)] = df2;
1814    
1815                                new_dFdv[INDEX4(i,0,q,7, numF,DIM,numQuadNodes)] = df0;
1816                                new_dFdv[INDEX4(i,1,q,7, numF,DIM,numQuadNodes)] = df1;
1817                                new_dFdv[INDEX4(i,2,q,7, numF,DIM,numQuadNodes)] = df2;
1818                            }
1819    
1820            }
1821        } else {
1822            Finley_setError(MEMORY_ERROR,"Finley_Quad_MacroHex: unable to create quadrature scheme for macro element.");
1823            return -1;
1824        }
1825        #undef DIM
1826        return numSubElements*numQuadNodes;
1827    }

Legend:
Removed from v.1342  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26