Changeset 5624


Ignore:
Timestamp:
08/30/10 15:53:35 (15 years ago)
Author:
Mathieu Morlighem
Message:

renamed some Gauss points for new Gauss class

Location:
issm/trunk/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5596 r5624  
    563563
    564564        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    565         GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
     565        GaussLegendreTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
    566566
    567567        /*Retrieve all inputs we will be needing: */
    … …  
    24992499          get tria gaussian points as well as segment gaussian points. For tria gaussian
    25002500          points, order of integration is 2, because we need to integrate the product tB*D*B'
    2501           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     2501          which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    25022502          points, same deal, which yields 3 gaussian points.*/
    25032503
    … …  
    25052505        num_vert_gauss=5;
    25062506
    2507         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     2507        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    25082508
    25092509        /*Retrieve all inputs we will be needing: */
    … …  
    28302830          get tria gaussian points as well as segment gaussian points. For tria gaussian
    28312831          points, order of integration is 2, because we need to integrate the product tB*D*B'
    2832           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     2832          which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    28332833          points, same deal, which yields 3 gaussian points.*/
    28342834
    … …  
    28362836        num_vert_gauss=5;
    28372837
    2838         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     2838        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    28392839
    28402840        /*Retrieve all inputs we will be needing: */
    … …  
    30213021                get tria gaussian points as well as segment gaussian points. For tria gaussian
    30223022                points, order of integration is 2, because we need to integrate the product tB*D*B'
    3023                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3023                which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    30243024                points, same deal, which yields 3 gaussian points.*/
    30253025
    … …  
    30283028
    30293029        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3030         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     3030        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    30313031
    30323032        /*Retrieve all inputs we will be needing: */
    … …  
    30973097
    30983098                xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
    3099                 GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
     3099                GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    31003100
    31013101                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    32493249          get tria gaussian points as well as segment gaussian points. For tria gaussian
    32503250          points, order of integration is 2, because we need to integrate the product tB*D*B'
    3251           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3251          which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    32523252          points, same deal, which yields 3 gaussian points.*/
    32533253
    … …  
    32553255        num_vert_gauss=2;
    32563256
    3257         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     3257        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    32583258
    32593259        /* Start  looping on the number of gaussian points: */
    … …  
    35023502                get tria gaussian points as well as segment gaussian points. For tria gaussian
    35033503                points, order of integration is 2, because we need to integrate the product tB*D*B'
    3504                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3504                which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    35053505                points, same deal, which yields 3 gaussian points.: */
    35063506
    … …  
    35093509        num_vert_gauss=2;
    35103510
    3511         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     3511        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    35123512
    35133513        /* Start  looping on the number of gaussian points: */
    … …  
    38613861        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    38623862        num_gauss=3;
    3863         GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
     3863        gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    38643864
    38653865        /*Loop on the three segments*/
    … …  
    40344034        num_vert_gauss=3;
    40354035
    4036         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     4036        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    40374037
    40384038        /*Retrieve all inputs we will be needing: */
    … …  
    41974197                get tria gaussian points as well as segment gaussian points. For tria gaussian
    41984198                points, order of integration is 2, because we need to integrate the product tB*D*B'
    4199                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     4199                which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    42004200                points, same deal, which yields 3 gaussian points.*/
    42014201
    42024202        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4203         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     4203        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    42044204
    42054205        /*Retrieve all inputs we will be needing: */
    … …  
    42854285
    42864286                xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
    4287                 GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
     4287                GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    42884288
    42894289                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    44374437        num_vert_gauss=2;
    44384438
    4439         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     4439        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    44404440
    44414441        /*Retrieve all inputs we will be needing: */
    … …  
    46614661                get tria gaussian points as well as segment gaussian points. For tria gaussian
    46624662                points, order of integration is 2, because we need to integrate the product tB*D*B'
    4663                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     4663                which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    46644664                points, same deal, which yields 3 gaussian points.: */
    46654665
    46664666        /*Get gaussian points: */
    4667         GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     4667        gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    46684668
    46694669        /*Retrieve all inputs we will be needing: */

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r5579 r5624  
    727727
    728728        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    729         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     729        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    730730
    731731        /* Start  looping on the number of gaussian points: */
    … …  
    10811081
    10821082        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1083         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     1083        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    10841084
    10851085        /*Retrieve all inputs we will be needing: */
    … …  
    12541254
    12551255        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1256         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     1256        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    12571257
    12581258        /*Retrieve all inputs we will be needing: */
    … …  
    20032003
    20042004        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2005         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2005        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    20062006
    20072007        /* Start  looping on the number of gaussian points: */
    … …  
    21282128
    21292129        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2130         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2130        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    21312131
    21322132        /* Start  looping on the number of gaussian points: */
    … …  
    22602260
    22612261        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2262         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2262        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    22632263
    22642264        /* Start  looping on the number of gaussian points: */
    … …  
    23902390
    23912391        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2392         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2392        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    23932393
    23942394        /* Start  looping on the number of gaussian points: */
    … …  
    25222522
    25232523        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2524         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2524        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    25252525
    25262526        /* Start  looping on the number of gaussian points: */
    … …  
    26532653
    26542654        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2655         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2655        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    26562656
    26572657        /* Start  looping on the number of gaussian points: */
    … …  
    31053105
    31063106        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3107         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3107        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    31083108
    31093109        /* Start  looping on the number of gaussian points: */
    … …  
    32323232
    32333233        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3234         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3234        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    32353235
    32363236        /*Retrieve all inputs we will be needed*/
    … …  
    33963396
    33973397        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3398         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3398        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    33993399
    34003400        /* Start  looping on the number of gaussian points: */
    … …  
    35373537
    35383538        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3539         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3539        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    35403540
    35413541        /*Retrieve all inputs we will be needing: */
    … …  
    36953695
    36963696        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3697         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3697        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    36983698
    36993699        /* Start  looping on the number of gaussian points: */
    … …  
    38383838
    38393839        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3840         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3840        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    38413841
    38423842        /* Start  looping on the number of gaussian points: */
    … …  
    39793979
    39803980        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3981         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3981        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    39823982
    39833983        /* Start  looping on the number of gaussian points: */
    … …  
    41184118
    41194119        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4120         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4120        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    41214121
    41224122        /*Build normal vector to the surface:*/
    … …  
    42324232
    42334233        /* Get gaussian points and weights: */
    4234         GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4234        GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    42354235
    42364236        /* Start looping on the number of gauss  (nodes on the bedrock) */
    … …  
    43584358
    43594359        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4360         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4360        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    43614361
    43624362        /* Start  looping on the number of gaussian points: */
    … …  
    45064506
    45074507        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4508         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4508        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    45094509
    45104510        /*Retrieve all inputs we will be needed*/
    … …  
    46144614
    46154615        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4616         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4616        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    46174617
    46184618        /* Start  looping on the number of gaussian points: */
    … …  
    47054705        heatcapacity=matpar->GetHeatCapacity();
    47064706
    4707         GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4707        GaussLegendreTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    47084708
    47094709        /* Start looping on the number of gauss (nodes on the bedrock) */
    … …  
    47904790
    47914791        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4792         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4792        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    47934793
    47944794        /*retrieve inputs :*/
    … …  
    48744874
    48754875        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4876         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4876        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    48774877
    48784878        /*Retrieve all inputs we will be needing: */
    … …  
    49564956
    49574957        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4958         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4958        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    49594959
    49604960        /*Retrieve all inputs we will be needing: */
    … …  
    50475047
    50485048        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5049         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5049        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    50505050
    50515051        /*Retrieve all inputs we will be needing: */
    … …  
    51665166
    51675167        /* Get gaussian points and weights: */
    5168         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
     5168        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    51695169
    51705170        /* Start  looping on the number of gaussian points: */
    … …  
    52755275
    52765276        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5277         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5277        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    52785278
    52795279        /* Start  looping on the number of gaussian points: */
    … …  
    55115511
    55125512        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5513         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5513        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    55145514
    55155515        /* Start  looping on the number of gaussian points: */
    … …  
    57505750
    57515751        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5752         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5752        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    57535753
    57545754        /* Start  looping on the number of gaussian points: */
    … …  
    59085908
    59095909        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5910         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5910        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    59115911
    59125912        /*Retrieve all inputs we will be needing: */
    … …  
    59965996
    59975997        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5998         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     5998        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    59995999
    60006000        /*Retrieve all inputs we will be needing: */
    … …  
    60826082
    60836083        /* Get gaussian points and weights: */
    6084         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
     6084        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    60856085
    60866086        /*Retrieve all inputs we will be needing: */
    … …  
    61936193        /* Ice/ocean heat exchange flux on ice shelf base */
    61946194
    6195         GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     6195        GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    61966196
    61976197        /*Retrieve all inputs we will be needing: */
    … …  
    63256325       
    63266326        /* Ice/ocean heat exchange flux on ice shelf base */
    6327         GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     6327        GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    63286328
    63296329        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    … …  
    68446844
    68456845        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    6846         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     6846        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    68476847
    68486848        /* Start  looping on the number of gaussian points: */

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5596 r5624  
    456456        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    457457        num_gauss=3;
    458         GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
     458        gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    459459
    460460        for(ig=0;ig<num_gauss;ig++){
    … …  
    916916
    917917        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    918         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     918        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    919919       
    920920        //Recover the surface of the four nodes
    … …  
    11771177
    11781178        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1179         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1179        GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    11801180       
    11811181        //Recover the surface of the four nodes

  • issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r5578 r5624  
    455455        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    456456        num_gauss=2;
    457         GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     457        gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
    458458
    459459        /* Start  looping on the number of gaussian points: */
    … …  
    598598        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    599599        num_gauss=2;
    600         GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     600        gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
    601601
    602602        /* Start  looping on the number of gaussian points: */
    … …  
    745745        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    746746        num_gauss=2;
    747         GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     747        gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
    748748
    749749        /* Start  looping on the number of gaussian points: */

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

    r4417 r5624  
    1010
    1111/*General Gauss points*/
    12 /*FUNCTION GaussLegendre {{{1*/
    13 void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus){
     12/*FUNCTION GaussLegendreLinear {{{1*/
     13void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus){
    1414        /* Gauss-Legendre quadrature points.
    1515
    … …  
    8888        }
    8989}/*}}}1*/
    90 /*FUNCTION GaussLobatto{{{1*/
    91 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
    92         /*Gauss-Lobatto quadrature points.
    93 
    94           The recurrence coefficients for Legendre polynomials on (-1,1)
    95           are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
    96 
    97           alpha(i)=0.
    98           beta (i)=1./(4.-1./(i-1)^2))
    99 
    100           and then modified for the Gauss-Lobatto quadrature rule on (-1,1)
    101           (from the ORTHPOL subroutine LOB).
    102 
    103           For degree p, the required number of Gauss-Lobatto points is
    104           n>=(p+1)/2+1 (one more than Gauss-Legendre).*/
    105 
    106         int i;
    107         double *alpha,*beta;
    108         double left=-1.,right= 1.;
    109         double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
    110 
    111         /*p= 1, npoint= 1 (Gauss-Legendre)*/
    112         static double wgt1[]={2.000000000000000};
    113         static double xi1[]={0.000000000000000};
    114 
    115         /*p= 1, npoint= 2*/
    116         static double wgt2[]={1.000000000000000, 1.000000000000000};
    117         static double xi2[]={-1.000000000000000, 1.000000000000000};
    118 
    119         /*p= 3, npoint= 3*/
    120         static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
    121         static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
    122 
    123         /*p= 5, npoint= 4*/
    124         static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
    125         static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
    126 
    127         /*p= 7, npoint= 5*/
    128         static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
    129         static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
    130 
    131         static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
    132         static double* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
    133 
    134         static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(double),
    135                 sizeof(wgt2 )/sizeof(double),
    136                 sizeof(wgt3 )/sizeof(double),
    137                 sizeof(wgt4 )/sizeof(double),
    138                 sizeof(wgt5 )/sizeof(double)};
    139 
    140         //      _printf_("Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
    141         *pxgaus = (double *) xmalloc(ngaus*sizeof(double));
    142         *pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
    143 
    144         /*  check to see if Gauss points need to be calculated  */
    145         if (ngaus <= MAX_LINE_GLOB_PTS) {
    146 
    147                 /*  copy the points from the static arrays (noting that the pointers
    148                          could be returned directly, but then the calling function would
    149                          have to know to not free them)  */
    150                 for (i=0; i<ngaus; i++) {
    151                         (*pxgaus)[i]=xip [ngaus-1][i];
    152                         (*pxwgt )[i]=wgtp[ngaus-1][i];
    153                 }
    154         }
    155         else {
    156 
    157                 /*  calculate the Gauss points using recurrence relations  */
    158                 alpha=(double *) xmalloc(ngaus*sizeof(double));
    159                 beta =(double *) xmalloc(ngaus*sizeof(double));
    160 
    161                 /*  calculate the Legendre recurrence coefficients  */
    162                 alpha[0]=0.;
    163                 beta [0]=2.;
    164 
    165                 for (i=1; i<ngaus; i++) {
    166                         alpha[i]=0.;
    167                         beta [i]=1./(4.-1./(i*i));
    168                 }
    169 
    170                 /*  calculate the Gauss-Lobatto quadrature formula  */
    171                 for (i=0; i<ngaus-1; i++) {
    172                         pm1l=p0l;
    173                         p0l=p1l;
    174                         pm1r=p0r;
    175                         p0r=p1r;
    176                         p1l=(left -alpha[i])*p0l-beta[i]*pm1l;
    177                         p1r=(right-alpha[i])*p0r-beta[i]*pm1r;
    178                 }
    179 
    180                 /*  Normalize system to prevent underflow:
    181                          [ p1l p0l ]{ a } = {left *p1l}
    182                          [ p1r p0r ]{ b }   {right*p1r}
    183                          dividing by p1l in the first equation and p1r in the second.  */
    184 
    185                 //              det=p1l*p0r-p1r*p0l;
    186                 det=p0r/p1r-p0l/p1l;
    187                 //              alpha[ngaus-1]=(left*p1l*p0r-right*p1r*p0l)/det;
    188                 //              beta [ngaus-1]=(right-left)*p1l*p1r/det;
    189                 alpha[ngaus-1]=(left *(p0r/p1r)-right*(p0l/p1l))/det;
    190                 beta [ngaus-1]=(right          -left           )/det;
    191 
    192                 /*  calculate the Gauss points  */
    193                 GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
    194                 xfree((void **)&beta );
    195                 xfree((void **)&alpha);
    196         }
    197 
    198 }/*}}}1*/
    199 /*FUNCTION GaussRecur{{{1*/
    200 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
    201         /*Gauss quadrature points from recursion coefficients.
    202          *
    203          *The routine uses the algorithm from the ORTHPOL routine GAUSS, which
    204          *finds the eigenvalues of a tridiagonal matrix.*/
    205 
    206         /*Intermediaries*/
    207         int i,j,k,l,m,ii,mml,iter;
    208         double p,g,r,s,c,f,b;
    209         double* work;
    210 
    211         if (n==1){
    212                 zero[0]  =alpha[0];
    213                 weight[0]=beta[0];
    214                 return;
    215         }
    216 
    217         work=(double*)xmalloc(n*sizeof(double));
    218 
    219         zero[0]  =alpha[0];
    220         weight[0]=1.;
    221         work[n-1]=0.;
    222         for (i=1; i<=n-1; i++){
    223                 zero[i]=alpha[i];
    224                 work[i-1]=sqrt(beta[i]);
    225                 weight[i]=0;
    226         }
    227 
    228         for (l=0; l<=n-1; l++){
    229                 iter=0;
    230                 do {
    231 
    232                         /*  Look for a small subdiagonal element.  */
    233                         for (m=l; m<=n-1; m++) {
    234                                 if (m == n-1) break;
    235                                 if (fabs(work[m])
    236                                                         <= DBL_EPSILON*(fabs(zero[m])+fabs(zero[m+1])))
    237                                  break;
    238                         }
    239                         p=zero[l];
    240                         if (m==l) break;
    241                         ++iter;
    242 
    243                         /*  Form shift.  */
    244                         g=(zero[l+1]-p)/(2.*work[l]);
    245                         r=sqrt(g*g+1.);
    246                         //                      g=zero[m]-p+work[l]/(g+FortranSign(r,g));
    247                         g=zero[m]-p+work[l]/(g+(g>=0 ? fabs(r) : -fabs(r)));
    248                         s=1.;
    249                         c=1.;
    250                         p=0.;
    251                         mml=m-l;
    252 
    253                         /*  For i=m-1 step -1 until l do ...  */
    254                         for (ii=1; ii<=mml; ii++) {
    255                                 i=m-ii;
    256                                 f=s*work[i];
    257                                 b=c*work[i];
    258                                 if (fabs(f) >= fabs(g)) {
    259                                         c=g/f;
    260                                         r=sqrt(c*c+1.);
    261                                         work[i+1]=f*r;
    262                                         s=1./r;
    263                                         c*=s;
    264                                 }
    265                                 else {
    266                                         s=f/g;
    267                                         r=sqrt(s*s+1.);
    268                                         work[i+1]=g*r;
    269                                         c=1./r;
    270                                         s*=c;
    271                                 }
    272                                 g=zero[i+1]-p;
    273                                 r=(zero[i]-g)*s+2.*c*b;
    274                                 p=s*r;
    275                                 zero[i+1]=g+p;
    276                                 g=c*r-b;
    277 
    278                                 /*  Form first component of vector.  */
    279                                 f=weight[i+1];
    280                                 weight[i+1]=s*weight[i]+c*f;
    281                                 weight[i  ]=c*weight[i]-s*f;
    282                         }
    283                         zero[l]-=p;
    284                         work[l]=g;
    285                         work[m]=0.;
    286                 } while (iter < MAX_GAUS_ITER);
    287                 if (iter >= MAX_GAUS_ITER) {
    288                         xfree((void **)&work);
    289                         ISSMERROR("%s%i"," Max iterations exceeded for l=",MAX_GAUS_ITER);
    290                 }
    291         }
    292 
    293         /*  Order eigenvalues and eigenvectors.  */
    294         for (i=0;i<n-1;i++) {
    295                 k=i;
    296                 p=zero[i];
    297                 for (j=i+1;j<n;j++){
    298                         if (zero[j] < p){
    299                                 k=j;
    300                                 p=zero[j];
    301                         }
    302                 }
    303                 if (k > i){
    304                         p=zero[i];
    305                         zero[i]=zero[k];
    306                         zero[k]=p;
    307                         p=weight[i];
    308                         weight[i]=weight[k];
    309                         weight[k]=p;
    310                 }
    311         }
    312         for (i=0;i<n;i++){
    313                 weight[i]=beta[0]*weight[i]*weight[i];
    314         }
    315 
    316         /*Cleanup*/
    317         xfree((void **)&work);
    318 
    319 }/*}}}1*/
    320 
    321 /*Element Gauss points*/
    322 /*FUNCTION GaussQuad{{{1*/
    323 void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) {
    324         /*Gauss quadrature points for the quadrilaterial.*/
    325 
    326         /*get the gauss points using the product of two line rules  */
    327         GaussLegendre(pxgaus, pxwgt, nigaus);
    328         GaussLegendre(pegaus, pewgt, njgaus);
    329 }/*}}}1*/
    330 /*FUNCTION GaussTria{{{1*/
    331 void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
     90/*FUNCTION GaussLegendreTria{{{1*/
     91void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
    33292        /*Gauss quadrature points for the triangle.
    33393
    … …  
    13791139                sizeof(wgt20)/sizeof(double)};
    13801140
    1381         //      _printf_("GaussTria: iord=%d\n",iord);
     1141        //      _printf_("GaussLegendreTria: iord=%d\n",iord);
    13821142
    13831143        /*  check to see if Gauss points need to be calculated  */
    … …  
    14141174
    14151175                /*  get the gauss points in each direction  */
    1416                 GaussLegendre(&xgaus, &xwgt, nigaus);
     1176                GaussLegendreLinear(&xgaus, &xwgt, nigaus);
    14171177
    14181178                egaus=xgaus;
    … …  
    14391199        }
    14401200
    1441         //      _printf_("GaussTria - ngaus=%d\n",*pngaus);
     1201        //      _printf_("GaussLegendreTria - ngaus=%d\n",*pngaus);
    14421202        //      for (i=0; i<*pngaus; i++)
    14431203        //              _printf_("i=%d: l1gaus=%f,l2gaus=%f,l3gaus=%f,wgt=%f\n",
    … …  
    14461206        return;
    14471207}/*}}}1*/
    1448 /*FUNCTION GaussHexa{{{1*/
    1449 void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
    1450         /*Gauss quadrature points for the hexahedron.*/
    1451 
    1452         /*  get the gauss points using the product of three line rules  */
    1453         GaussLegendre(pxgaus, pxwgt, nigaus);
    1454         GaussLegendre(pegaus, pewgt, njgaus);
    1455         GaussLegendre(pzgaus, pzwgt, nkgaus);
    1456 }/*}}}1*/
    1457 /*FUNCTION GaussPenta{{{1*/
    1458 void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
    1459         /*Gauss quadrature points for the pentahedron.*/
    1460 
    1461         /*  get the gauss points using the product of triangular and line rules  */
    1462         GaussTria(pngaus, pl1, pl2, pl3, pwgt, iord);
    1463         GaussLegendre(pzgaus, pzwgt, nkgaus);
    1464        
    1465 }/*}}}1*/
    1466 /*FUNCTION GaussTetra{{{1*/
    1467 void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
     1208/*FUNCTION GaussLegendreTetra{{{1*/
     1209void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
    14681210        /* Gauss quadrature points for the tetrahedron.
    14691211
    … …  
    16561398                sizeof(wgt6 )/sizeof(double)};
    16571399
    1658         //      _printf_("GaussTetra: iord=%d\n",iord);
     1400        //      _printf_("GaussLegendreTetra: iord=%d\n",iord);
    16591401
    16601402        /*  check to see if Gauss points need to be calculated  */
    … …  
    16951437
    16961438                /*  get the gauss points in each direction  */
    1697                 GaussLegendre(&xgaus, &xwgt, nigaus);
     1439                GaussLegendreLinear(&xgaus, &xwgt, nigaus);
    16981440
    16991441                egaus=xgaus;
    … …  
    17291471        }
    17301472}/*}}}1*/
    1731 /*FUNCTION GaussSegment{{{1*/
    1732 void GaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss){
    1733         /* GaussSegment:
    1734          * This routine computes gaussian points on a reference segment from -1 to 1.
    1735          * The order p of integration depends on the number of gaussian points and can be computed from the polynomial degree p that needs to be integrated.
    1736          * numgauss= (p+1) /2 */
    1737 
    1738         /*WARNING: This is a copy of GaussLegendre: Merge the two!*/
    1739 
    1740         double* segment_coord=NULL;
    1741         double* gauss_weights=NULL;
     1473/*FUNCTION GaussLobatto{{{1*/
     1474void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
     1475        /*Gauss-Lobatto quadrature points.
     1476
     1477          The recurrence coefficients for Legendre polynomials on (-1,1)
     1478          are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
     1479
     1480          alpha(i)=0.
     1481          beta (i)=1./(4.-1./(i-1)^2))
     1482
     1483          and then modified for the Gauss-Lobatto quadrature rule on (-1,1)
     1484          (from the ORTHPOL subroutine LOB).
     1485
     1486          For degree p, the required number of Gauss-Lobatto points is
     1487          n>=(p+1)/2+1 (one more than Gauss-Legendre).*/
     1488
     1489        int i;
     1490        double *alpha,*beta;
     1491        double left=-1.,right= 1.;
     1492        double p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
     1493
     1494        /*p= 1, npoint= 1 (Gauss-Legendre)*/
     1495        static double wgt1[]={2.000000000000000};
     1496        static double xi1[]={0.000000000000000};
     1497
     1498        /*p= 1, npoint= 2*/
     1499        static double wgt2[]={1.000000000000000, 1.000000000000000};
     1500        static double xi2[]={-1.000000000000000, 1.000000000000000};
     1501
     1502        /*p= 3, npoint= 3*/
     1503        static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
     1504        static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
     1505
     1506        /*p= 5, npoint= 4*/
     1507        static double wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
     1508        static double xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
     1509
     1510        /*p= 7, npoint= 5*/
     1511        static double wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
     1512        static double xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
     1513
     1514        static double* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
     1515        static double* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
     1516
     1517        static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(double),
     1518                sizeof(wgt2 )/sizeof(double),
     1519                sizeof(wgt3 )/sizeof(double),
     1520                sizeof(wgt4 )/sizeof(double),
     1521                sizeof(wgt5 )/sizeof(double)};
     1522
     1523        //      _printf_("Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
     1524        *pxgaus = (double *) xmalloc(ngaus*sizeof(double));
     1525        *pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
     1526
     1527        /*  check to see if Gauss points need to be calculated  */
     1528        if (ngaus <= MAX_LINE_GLOB_PTS) {
     1529
     1530                /*  copy the points from the static arrays (noting that the pointers
     1531                         could be returned directly, but then the calling function would
     1532                         have to know to not free them)  */
     1533                for (i=0; i<ngaus; i++) {
     1534                        (*pxgaus)[i]=xip [ngaus-1][i];
     1535                        (*pxwgt )[i]=wgtp[ngaus-1][i];
     1536                }
     1537        }
     1538        else {
     1539
     1540                /*  calculate the Gauss points using recurrence relations  */
     1541                alpha=(double *) xmalloc(ngaus*sizeof(double));
     1542                beta =(double *) xmalloc(ngaus*sizeof(double));
     1543
     1544                /*  calculate the Legendre recurrence coefficients  */
     1545                alpha[0]=0.;
     1546                beta [0]=2.;
     1547
     1548                for (i=1; i<ngaus; i++) {
     1549                        alpha[i]=0.;
     1550                        beta [i]=1./(4.-1./(i*i));
     1551                }
     1552
     1553                /*  calculate the Gauss-Lobatto quadrature formula  */
     1554                for (i=0; i<ngaus-1; i++) {
     1555                        pm1l=p0l;
     1556                        p0l=p1l;
     1557                        pm1r=p0r;
     1558                        p0r=p1r;
     1559                        p1l=(left -alpha[i])*p0l-beta[i]*pm1l;
     1560                        p1r=(right-alpha[i])*p0r-beta[i]*pm1r;
     1561                }
     1562
     1563                /*  Normalize system to prevent underflow:
     1564                         [ p1l p0l ]{ a } = {left *p1l}
     1565                         [ p1r p0r ]{ b }   {right*p1r}
     1566                         dividing by p1l in the first equation and p1r in the second.  */
     1567
     1568                //              det=p1l*p0r-p1r*p0l;
     1569                det=p0r/p1r-p0l/p1l;
     1570                //              alpha[ngaus-1]=(left*p1l*p0r-right*p1r*p0l)/det;
     1571                //              beta [ngaus-1]=(right-left)*p1l*p1r/det;
     1572                alpha[ngaus-1]=(left *(p0r/p1r)-right*(p0l/p1l))/det;
     1573                beta [ngaus-1]=(right          -left           )/det;
     1574
     1575                /*  calculate the Gauss points  */
     1576                GaussRecur(*pxgaus, *pxwgt, ngaus, alpha, beta );
     1577                xfree((void **)&beta );
     1578                xfree((void **)&alpha);
     1579        }
     1580
     1581}/*}}}1*/
     1582/*FUNCTION GaussRecur{{{1*/
     1583void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
     1584        /*Gauss quadrature points from recursion coefficients.
     1585         *
     1586         *The routine uses the algorithm from the ORTHPOL routine GAUSS, which
     1587         *finds the eigenvalues of a tridiagonal matrix.*/
     1588
     1589        /*Intermediaries*/
     1590        int i,j,k,l,m,ii,mml,iter;
     1591        double p,g,r,s,c,f,b;
     1592        double* work;
     1593
     1594        if (n==1){
     1595                zero[0]  =alpha[0];
     1596                weight[0]=beta[0];
     1597                return;
     1598        }
     1599
     1600        work=(double*)xmalloc(n*sizeof(double));
     1601
     1602        zero[0]  =alpha[0];
     1603        weight[0]=1.;
     1604        work[n-1]=0.;
     1605        for (i=1; i<=n-1; i++){
     1606                zero[i]=alpha[i];
     1607                work[i-1]=sqrt(beta[i]);
     1608                weight[i]=0;
     1609        }
     1610
     1611        for (l=0; l<=n-1; l++){
     1612                iter=0;
     1613                do {
     1614
     1615                        /*  Look for a small subdiagonal element.  */
     1616                        for (m=l; m<=n-1; m++) {
     1617                                if (m == n-1) break;
     1618                                if (fabs(work[m])
     1619                                                        <= DBL_EPSILON*(fabs(zero[m])+fabs(zero[m+1])))
     1620                                 break;
     1621                        }
     1622                        p=zero[l];
     1623                        if (m==l) break;
     1624                        ++iter;
     1625
     1626                        /*  Form shift.  */
     1627                        g=(zero[l+1]-p)/(2.*work[l]);
     1628                        r=sqrt(g*g+1.);
     1629                        //                      g=zero[m]-p+work[l]/(g+FortranSign(r,g));
     1630                        g=zero[m]-p+work[l]/(g+(g>=0 ? fabs(r) : -fabs(r)));
     1631                        s=1.;
     1632                        c=1.;
     1633                        p=0.;
     1634                        mml=m-l;
     1635
     1636                        /*  For i=m-1 step -1 until l do ...  */
     1637                        for (ii=1; ii<=mml; ii++) {
     1638                                i=m-ii;
     1639                                f=s*work[i];
     1640                                b=c*work[i];
     1641                                if (fabs(f) >= fabs(g)) {
     1642                                        c=g/f;
     1643                                        r=sqrt(c*c+1.);
     1644                                        work[i+1]=f*r;
     1645                                        s=1./r;
     1646                                        c*=s;
     1647                                }
     1648                                else {
     1649                                        s=f/g;
     1650                                        r=sqrt(s*s+1.);
     1651                                        work[i+1]=g*r;
     1652                                        c=1./r;
     1653                                        s*=c;
     1654                                }
     1655                                g=zero[i+1]-p;
     1656                                r=(zero[i]-g)*s+2.*c*b;
     1657                                p=s*r;
     1658                                zero[i+1]=g+p;
     1659                                g=c*r-b;
     1660
     1661                                /*  Form first component of vector.  */
     1662                                f=weight[i+1];
     1663                                weight[i+1]=s*weight[i]+c*f;
     1664                                weight[i  ]=c*weight[i]-s*f;
     1665                        }
     1666                        zero[l]-=p;
     1667                        work[l]=g;
     1668                        work[m]=0.;
     1669                } while (iter < MAX_GAUS_ITER);
     1670                if (iter >= MAX_GAUS_ITER) {
     1671                        xfree((void **)&work);
     1672                        ISSMERROR("%s%i"," Max iterations exceeded for l=",MAX_GAUS_ITER);
     1673                }
     1674        }
     1675
     1676        /*  Order eigenvalues and eigenvectors.  */
     1677        for (i=0;i<n-1;i++) {
     1678                k=i;
     1679                p=zero[i];
     1680                for (j=i+1;j<n;j++){
     1681                        if (zero[j] < p){
     1682                                k=j;
     1683                                p=zero[j];
     1684                        }
     1685                }
     1686                if (k > i){
     1687                        p=zero[i];
     1688                        zero[i]=zero[k];
     1689                        zero[k]=p;
     1690                        p=weight[i];
     1691                        weight[i]=weight[k];
     1692                        weight[k]=p;
     1693                }
     1694        }
     1695        for (i=0;i<n;i++){
     1696                weight[i]=beta[0]*weight[i]*weight[i];
     1697        }
     1698
     1699        /*Cleanup*/
     1700        xfree((void **)&work);
     1701
     1702}/*}}}1*/
     1703
     1704/*Element Gauss points TO BE REMOVED*/
     1705/*FUNCTION gaussQuad{{{1*/
     1706void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) {
     1707        /*Gauss quadrature points for the quadrilaterial.*/
     1708
     1709        /*get the gauss points using the product of two line rules  */
     1710        GaussLegendreLinear(pxgaus, pxwgt, nigaus);
     1711        GaussLegendreLinear(pegaus, pewgt, njgaus);
     1712}/*}}}1*/
     1713/*FUNCTION gaussHexa{{{1*/
     1714void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
     1715        /*Gauss quadrature points for the hexahedron.*/
     1716
     1717        /*  get the gauss points using the product of three line rules  */
     1718        GaussLegendreLinear(pxgaus, pxwgt, nigaus);
     1719        GaussLegendreLinear(pegaus, pewgt, njgaus);
     1720        GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
     1721}/*}}}1*/
     1722/*FUNCTION gaussPenta{{{1*/
     1723void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
     1724        /*Gauss quadrature points for the pentahedron.*/
     1725
     1726        /*  get the gauss points using the product of triangular and line rules  */
     1727        GaussLegendreTria(pngaus, pl1, pl2, pl3, pwgt, iord);
     1728        GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
    17421729       
    1743        
    1744         if (num_gauss==1){
    1745                 //order=1, num_gauss=1. Can integrate polynomials of degree 0 to 1
    1746 
    1747                 gauss_weights=(double*)xmalloc(sizeof(double));
    1748                 gauss_weights[0]= 2.0;
    1749                
    1750                 segment_coord=(double*)xmalloc(sizeof(double));
    1751                 segment_coord[0]= 0.0;
    1752                
    1753         }
    1754         else if (num_gauss==2){
    1755                 //order=2, num_gauss=3. Can integrate polynomials of degree 0 to 3
    1756                
    1757                 gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
    1758                 gauss_weights[0]=1.0;
    1759                 gauss_weights[1]=1.0;
    1760        
    1761                 segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
    1762                 segment_coord[0]=-0.577350269189626;
    1763                 segment_coord[1]=0.577350269189626;
    1764        
    1765         }
    1766         else if (num_gauss==3){
    1767                 //order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
    1768 
    1769                 gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
    1770                 gauss_weights[0]=0.555555555555555556;
    1771                 gauss_weights[1]=0.88888888888889;
    1772                 gauss_weights[2]=0.555555555555555556;
    1773                
    1774                 segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
    1775                 segment_coord[0]=-0.774596669241483;
    1776                 segment_coord[1]=0.0;
    1777                 segment_coord[2]=0.774596669241483;
    1778                
    1779         }
    1780         else if (num_gauss==4){
    1781                 //order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
    1782 
    1783                 gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
    1784                 gauss_weights[0]=0.347854845137454;
    1785                 gauss_weights[1]=0.652145154862546;
    1786                 gauss_weights[2]=0.652145154862546;
    1787                 gauss_weights[3]=0.347854845137454;
    1788                
    1789                 segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
    1790                 segment_coord[0]=-0.861136311594053;
    1791                 segment_coord[1]=-0.339981043584856;
    1792                 segment_coord[2]=0.339981043584856;
    1793                 segment_coord[3]=0.861136311594053;
    1794                
    1795         }
    1796         else if (num_gauss==5){
    1797                 //order=3, num_gauss=4. Can integrate polynomials of degree 0 to 5
    1798 
    1799                 gauss_weights=(double*)xmalloc(num_gauss*sizeof(double));
    1800                 gauss_weights[0]=0.236926885056189;
    1801                 gauss_weights[1]=0.478628670499366;
    1802                 gauss_weights[2]=0.568888888888889;
    1803                 gauss_weights[3]=0.478628670499366;
    1804                 gauss_weights[4]=0.236926885056189;
    1805                
    1806                 segment_coord=(double*)xmalloc(num_gauss*sizeof(double));
    1807                 segment_coord[0]=-0.906179845938664;
    1808                 segment_coord[1]=-0.538469310105683;
    1809                 segment_coord[2]=0.0;
    1810                 segment_coord[3]=0.538469310105683;
    1811                 segment_coord[4]=0.906179845938664;
    1812                
    1813         }
    1814         else{
    1815                 _printf_("GaussSegment error message: order not supported yet");
    1816         }
    1817 
    1818         /*Assign output pointers: */
    1819         *pgauss_weights=gauss_weights;
    1820         *psegment_coord=segment_coord;
    1821         return;
    18221730}/*}}}1*/
     1731/*FUNCTION gaussSegment{{{1*/
     1732void gaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss){
     1733        /* gaussSegment*/
     1734
     1735        GaussLegendreLinear(psegment_coord,pgauss_weights,num_gauss);
     1736}/*}}}1*/

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

    r4417 r5624  
    77
    88#define    MAX_LINE_GAUS_PTS    4
    9 void GaussLegendre( double** pxgaus, double** pxwgt, int ngaus );
    10 
     9void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus );
     10#define    MAX_TRIA_SYM_ORD    20
     11void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
     12#define    MAX_TETRA_SYM_ORD    6
     13void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
    1114#define    MAX_LINE_GLOB_PTS    5
    1215void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus );
    13 
    1416#define MAX_GAUS_ITER   30
    1517void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );
    1618
    17 void GaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
    18 
    19 #define    MAX_TRIA_SYM_ORD    20
    20 void GaussTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
    21 
    22 void GaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
    23 
    24 void GaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
    25 
    26 #define    MAX_TETRA_SYM_ORD    6
    27 void GaussTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
    28 
    29 void GaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss);
     19void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
     20void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
     21void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
     22void gaussSegment(double** psegment_coord,double** pgauss_weights,int num_gauss);
    3023
    3124#endif
Note: See TracChangeset for help on using the changeset viewer.