Changeset 99


Ignore:
Timestamp:
04/28/09 15:26:17 (16 years ago)
Author:
Eric.Larour
Message:

New cross and norm for vectors.
Simplified GaussPoints.

Location:
issm/trunk/src/c/shared/Numerics
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/shared/Numerics/GaussPoints.cpp

    r1 r99  
    4141#undef __FUNCT__
    4242#define __FUNCT__ "GaussLegendre"
    43 int GaussLegendre( double** pxgaus,
    44                                    double** pxwgt,
    45                                    int ngaus )
    46 {
    47         int i,ierr=0;
     43void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus ) {
     44       
     45        int i;
    4846        double *alpha,*beta;
    4947
    … …  
    124122/*  calculate the Gauss points  */
    125123
    126                 if ((ierr=GaussRecur(*pxgaus,
    127                                                          *pxwgt,
    128                                                          ngaus,
    129                                                          alpha,
    130                                                          beta )))
    131                         _printf_("%s -- Error %d from GaussRecur for ngaus=%d.\n",
    132                                          __FUNCT__,ierr,ngaus);
    133 
     124                GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
    134125                xfree((void **)&beta );
    135126                xfree((void **)&alpha);
    136127        }
    137128
    138 //      for (i=0; i<ngaus; i++)
    139 //              _printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
    140 
    141         return(ierr);
    142129}
    143130
    … …  
    179166#define __FUNCT__ "GaussLobatto"
    180167
    181 int GaussLobatto( double** pxgaus,
    182                                   double** pxwgt,
    183                                   int ngaus )
    184 {
    185         int i,ierr=0;
     168void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
     169
     170        int i;
    186171        double *alpha,*beta;
    187172        double left=-1.,right= 1.;
    … …  
    297282/*  calculate the Gauss points  */
    298283
    299                 if ((ierr=GaussRecur(*pxgaus,
    300                                                          *pxwgt,
    301                                                          ngaus,
    302                                                          alpha,
    303                                                          beta )))
    304                         _printf_("%s -- Error %d from GaussRecur for ngaus=%d.\n",
    305                                          __FUNCT__,ierr,ngaus);
    306 
     284                GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
    307285                xfree((void **)&beta );
    308286                xfree((void **)&alpha);
    309287        }
    310288
    311 //      for (i=0; i<ngaus; i++)
    312 //              _printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
    313 
    314         return(ierr);
    315289}
    316290
    … …  
    339313#define __FUNCT__ "GaussRecur"
    340314
    341 int GaussRecur( double* zero,
    342                                 double* weight,
    343                                 int n,
    344                                 double* alpha,
    345                                 double* beta )
    346 {
     315void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
    347316        int i,j,k,l,m,ii,mml,iter;
    348317        double p,g,r,s,c,f,b;
    349318        double* work;
    350319
    351 //      _printf_("Gauss abscissas and weights n=%d\n",n);
    352 
    353320        if ( n == 1 ) {
    354321                zero[0]  =alpha[0];
    355322                weight[0]=beta[0];
    356                 return 1;
     323                return ;
    357324        }
    358325
    … …  
    463430        xfree((void **)&work);
    464431
    465         return 1;
    466432}
    467433
    … …  
    492458#define __FUNCT__ "GaussQuad"
    493459
    494 int GaussQuad( double** pxgaus,
    495                            double** pxwgt,
    496                            double** pegaus,
    497                            double** pewgt,
    498                            int nigaus,
    499                            int njgaus )
    500 {
    501         int i,ierr=0;
    502 
    503 //      _printf_("GaussQuad: nigaus=%d,njgaus=%d\n",
    504 //                       nigaus,njgaus);
     460void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) {
    505461
    506462/*  get the gauss points using the product of two line rules  */
    507463
    508         if ((ierr=GaussLegendre(pxgaus,
    509                                                         pxwgt,
    510                                                         nigaus))) {
    511                 _printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
    512                                  __FUNCT__,ierr,nigaus);
    513                 return(ierr);
    514         }
    515 //      for (i=0; i<nigaus; i++)
    516 //              _printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
    517 
    518         if ((ierr=GaussLegendre(pegaus,
    519                                                         pewgt,
    520                                                         njgaus))) {
    521                 _printf_("%s -- Error %d from GaussLegendre for njgaus=%d.\n",
    522                                  __FUNCT__,ierr,njgaus);
    523                 return(ierr);
    524         }
    525 //      for (i=0; i<njgaus; i++)
    526 //              _printf_("j=%d: egaus=%f,ewgt=%f\n",i,(*pegaus)[i],(*pewgt)[i]);
    527 
    528         return 1;
     464        GaussLegendre(pxgaus, pxwgt, nigaus);
     465
     466        GaussLegendre(pegaus, pewgt, njgaus);
    529467}
    530468
    … …  
    559497#undef __FUNCT__
    560498#define __FUNCT__ "GaussTria"
    561 int GaussTria( int* pngaus,
    562                            double** pl1,
    563                            double** pl2,
    564                            double** pl3,
    565                            double** pwgt,
    566                            int iord )
    567 {
    568         int i,j,ipt,nigaus,ierr=0;
     499void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
     500       
     501        int i,j,ipt,nigaus;
    569502        double xi,eta;
    570503        double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
    … …  
    16591592/*  get the gauss points in each direction  */
    16601593
    1661                 if ((ierr=GaussLegendre(&xgaus,
    1662                                                                 &xwgt,
    1663                                                                 nigaus))) {
    1664                         _printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
    1665                                          __FUNCT__,ierr,nigaus);
    1666                         return(ierr);
    1667                 }
     1594                GaussLegendre(&xgaus, &xwgt, nigaus);
    16681595
    16691596                egaus=xgaus;
    … …  
    16971624//                               i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pwgt)[i]);
    16981625
    1699         return 1;
     1626        return ;
    17001627}
    17011628
    … …  
    17301657#undef __FUNCT__
    17311658#define __FUNCT__ "GaussHexa"
    1732 int GaussHexa( double** pxgaus,
    1733                            double** pxwgt,
    1734                            double** pegaus,
    1735                            double** pewgt,
    1736                            double** pzgaus,
    1737                            double** pzwgt,
    1738                            int nigaus,
    1739                            int njgaus,
    1740                            int nkgaus )
    1741 {
    1742         int i,ierr=0;
    1743 
    1744 //      _printf_("GaussHexa: nigaus=%d,njgaus=%d,nkgaus=%d\n",
    1745 //                       nigaus,njgaus,nkgaus);
    1746 
    1747 /*  get the gauss points using the product of three line rules  */
    1748 
    1749         if ((ierr=GaussLegendre(pxgaus,
    1750                                                         pxwgt,
    1751                                                         nigaus))) {
    1752                 _printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
    1753                                  __FUNCT__,ierr,nigaus);
    1754                 return(ierr);
    1755         }
    1756 //      for (i=0; i<nigaus; i++)
    1757 //              _printf_("i=%d: xgaus=%f,xwgt=%f\n",i,(*pxgaus)[i],(*pxwgt)[i]);
    1758 
    1759         if ((ierr=GaussLegendre(pegaus,
    1760                                                         pewgt,
    1761                                                         njgaus))) {
    1762                 _printf_("%s -- Error %d from GaussLegendre for njgaus=%d.\n",
    1763                                  __FUNCT__,ierr,njgaus);
    1764                 return(ierr);
    1765         }
    1766 //      for (i=0; i<njgaus; i++)
    1767 //              _printf_("j=%d: egaus=%f,ewgt=%f\n",i,(*pegaus)[i],(*pewgt)[i]);
    1768 
    1769         if ((ierr=GaussLegendre(pzgaus,
    1770                                                         pzwgt,
    1771                                                         nkgaus))) {
    1772                 _printf_("%s -- Error %d from GaussLegendre for nkgaus=%d.\n",
    1773                                  __FUNCT__,ierr,nkgaus);
    1774                 return(ierr);
    1775         }
    1776 //      for (i=0; i<nkgaus; i++)
    1777 //              _printf_("k=%d: zgaus=%f,zwgt=%f\n",i,(*pzgaus)[i],(*pzwgt)[i]);
    1778 
    1779         return 1;
     1659void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
     1660
     1661        /*  get the gauss points using the product of three line rules  */
     1662
     1663        GaussLegendre(pxgaus, pxwgt, nigaus);
     1664
     1665        GaussLegendre(pegaus, pewgt, njgaus);
     1666
     1667        GaussLegendre(pzgaus, pzwgt, nkgaus);
    17801668}
    17811669
    … …  
    18111699#define __FUNCT__ "GaussPenta"
    18121700
    1813 int GaussPenta( int* pngaus,
    1814                                 double** pl1,
    1815                                 double** pl2,
    1816                                 double** pl3,
    1817                                 double** pwgt,
    1818                                 double** pzgaus,
    1819                                 double** pzwgt,
    1820                                 int iord,
    1821                                 int nkgaus )
    1822 {
    1823         int i,ierr=0;
    1824 
    1825 //      _printf_("GaussPenta: iord=%d,nkgaus=%d\n",
    1826 //                       iord,nkgaus);
     1701void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
    18271702
    18281703/*  get the gauss points using the product of triangular and line rules  */
    18291704
    1830         if ((ierr=GaussTria(pngaus,
    1831                                                 pl1,
    1832                                                 pl2,
    1833                                                 pl3,
    1834                                                 pwgt,
    1835                                                 iord))) {
    1836                 _printf_("%s -- Error %d from GaussTria for iord=%d.\n",
    1837                                  __FUNCT__,ierr,iord);
    1838                 return(ierr);
    1839         }
    1840 //      for (i=0; i<*pngaus; i++)
    1841 //              _printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,wgt=%f\n",
    1842 //                               i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pwgt)[i]);
    1843 
    1844         if ((ierr=GaussLegendre(pzgaus,
    1845                                                         pzwgt,
    1846                                                         nkgaus))) {
    1847                 _printf_("%s -- Error %d from GaussLegendre for nkgaus=%d.\n",
    1848                                  __FUNCT__,ierr,nkgaus);
    1849                 return(ierr);
    1850         }
    1851 //      for (i=0; i<nkgaus; i++)
    1852 //              _printf_("k=%d: zgaus=%f,zwgt=%f\n",i,(*pzgaus)[i],(*pzwgt)[i]);
    1853 
    1854         return 1;
     1705        GaussTria(pngaus, pl1, pl2, pl3, pwgt, iord);
     1706
     1707        GaussLegendre(pzgaus, pzwgt, nkgaus);
     1708       
    18551709}
    18561710
    … …  
    18911745#define __FUNCT__ "GaussTetra"
    18921746
    1893 int GaussTetra( int* pngaus,
    1894                                 double** pl1,
    1895                                 double** pl2,
    1896                                 double** pl3,
    1897                                 double** pl4,
    1898                                 double** pwgt,
    1899                                 int iord )
    1900 {
    1901         int i,j,k,ipt,nigaus,ierr=0;
     1747void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
     1748        int i,j,k,ipt,nigaus;
    19021749        double xi,eta,zeta;
    19031750        double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
    … …  
    21261973/*  get the gauss points in each direction  */
    21271974
    2128                 if ((ierr=GaussLegendre(&xgaus,
    2129                                                                 &xwgt,
    2130                                                                 nigaus))) {
    2131                         _printf_("%s -- Error %d from GaussLegendre for nigaus=%d.\n",
    2132                                          __FUNCT__,ierr,nigaus);
    2133                         return(ierr);
    2134                 }
     1975                GaussLegendre(&xgaus, &xwgt, nigaus);
    21351976
    21361977                egaus=xgaus;
    … …  
    21682009        }
    21692010
    2170 //      _printf_("GaussTetra - ngaus=%d\n",*pngaus);
    2171 //      for (i=0; i<*pngaus; i++)
    2172 //              _printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,l4gaus=%f,wgt=%f\n",
    2173 //                               i,(*pl1 )[i],(*pl2 )[i],(*pl3 )[i],(*pl4 )[i],(*pwgt)[i]);
    2174 
    2175         return 1;
    21762011}
    21772012

  • issm/trunk/src/c/shared/Numerics/GaussPoints.h

    r1 r99  
    88
    99#define    MAX_LINE_GAUS_PTS    4
    10 
    11 int GaussLegendre( double** pxgaus,
    12                                    double** pxwgt,
    13                                    int ngaus );
     10void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus );
    1411
    1512#define    MAX_LINE_GLOB_PTS    5
    16 
    17 int GaussLobatto( double** pxgaus,
    18                                   double** pxwgt,
    19                                   int ngaus );
     13void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus );
    2014
    2115#define MAX_GAUS_ITER   30
     16void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );
    2217
    23 int GaussRecur( double* zero,
    24                                 double* weight,
    25                                 int n,
    26                                 double* alpha,
    27                                 double* beta );
    28 
    29 int GaussQuad( double** pxgaus,
    30                            double** pxwgt,
    31                            double** pegaus,
    32                            double** pewgt,
    33                            int nigaus,
    34                            int njgaus );
     18void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
    3519
    3620#define    MAX_TRIA_SYM_ORD    20
     21void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
    3722
    38 int GaussTria( int* pngaus,
    39                            double** pl1,
    40                            double** pl2,
    41                            double** pl3,
    42                            double** pwgt,
    43                            int iord );
     23void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
    4424
    45 int GaussHexa( double** pxgaus,
    46                            double** pxwgt,
    47                            double** pegaus,
    48                            double** pewgt,
    49                            double** pzgaus,
    50                            double** pzwgt,
    51                            int nigaus,
    52                            int njgaus,
    53                            int nkgaus );
    54 
    55 int GaussPenta( int* pngaus,
    56                                 double** pl1,
    57                                 double** pl2,
    58                                 double** pl3,
    59                                 double** pwgt,
    60                                 double** pzgaus,
    61                                 double** pzwgt,
    62                                 int iord,
    63                                 int nkgaus );
     25void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
    6426
    6527#define    MAX_TETRA_SYM_ORD    6
     28void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
    6629
    67 int GaussTetra( int* pngaus,
    68                                 double** pl1,
    69                                 double** pl2,
    70                                 double** pl3,
    71                                 double** pl4,
    72                                 double** pwgt,
    73                                 int iord );
    7430void GaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss);
    7531

  • issm/trunk/src/c/shared/Numerics/numerics.h

    r8 r99  
    1313double OptFunc(double scalar, OptArgs* optargs);
    1414void BrentSearch(double* psearch_scalar,double* pJ,OptPars* optpars,double (*f)(double,OptArgs*), OptArgs* optargs);
     15void cross(double* result,double* vector1,double* vector2);
     16double norm(double* vector);
    1517
    1618#endif //ifndef _NUMERICS_H_
Note: See TracChangeset for help on using the changeset viewer.