Changeset 12424


Ignore:
Timestamp:
06/15/12 14:36:38 (13 years ago)
Author:
utke
Message:

type renames

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

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/objects/Gauss/GaussTria.cpp

    r12365 r12424  
    4242
    4343        /*Intermediaties*/
    44         double *seg_coords  = NULL;
    45         double *seg_weights = NULL;
     44        IssmPDouble *seg_coords  = NULL;
     45        IssmPDouble *seg_weights = NULL;
    4646        int     i,index3;
    4747
     
    5151
    5252        /*Allocate GaussTria fields*/
    53         coords1=(double*)xmalloc(numgauss*sizeof(double));
    54         coords2=(double*)xmalloc(numgauss*sizeof(double));
    55         coords3=(double*)xmalloc(numgauss*sizeof(double));
    56         weights=(double*)xmalloc(numgauss*sizeof(double));
     53        coords1=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
     54        coords2=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
     55        coords3=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
     56        weights=xNew<IssmDouble>(numgauss*sizeof(IssmDouble));
    5757
    5858        /*Reverse index1 and 2 if necessary*/
     
    7171        else if (index1==0 && index2==2){
    7272                for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
    73                 for(i=0;i<numgauss;i++) coords2[i]=0;
     73                for(i=0;i<numgauss;i++) coords2[i]= 0 ;
    7474                for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
    7575                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
     
    9797/*FUNCTION GaussTria::~GaussTria(){{{*/
    9898GaussTria::~GaussTria(){
    99         xfree((void**)&weights);
    100         xfree((void**)&coords1);
    101         xfree((void**)&coords2);
    102         xfree((void**)&coords3);
     99        xDelete<IssmDouble>(weights);
     100        xDelete<IssmDouble>(coords1);
     101        xDelete<IssmDouble>(coords2);
     102        xDelete<IssmDouble>(coords3);
    103103}
    104104/*}}}*/
     
    199199/*}}}*/
    200200/*FUNCTION GaussTria::GaussFromCoords{{{*/
    201 void GaussTria::GaussFromCoords(double x,double y,double* xyz_list){
     201void GaussTria::GaussFromCoords(IssmDouble x,IssmDouble y,IssmDouble* xyz_list){
    202202
    203203        /*Intermediaries*/
    204         double    area = 0;
    205         double    x1,y1,x2,y2,x3,y3;
     204        IssmDouble    area = 0;
     205        IssmDouble    x1,y1,x2,y2,x3,y3;
    206206
    207207        /*in debugging mode: check that the default constructor has been called*/
  • issm/trunk-jpl/src/c/objects/Gauss/GaussTria.h

    r12365 r12424  
    1515        private:
    1616                int numgauss;
    17                 double* weights;
    18                 double* coords1;
    19                 double* coords2;
    20                 double* coords3;
     17                IssmDouble* weights;
     18                IssmDouble* coords1;
     19                IssmDouble* coords2;
     20                IssmDouble* coords3;
    2121
    2222        public:
    23                 double weight;
    24                 double coord1;
    25                 double coord2;
    26                 double coord3;
     23                IssmDouble weight;
     24                IssmDouble coord1;
     25                IssmDouble coord2;
     26                IssmDouble coord3;
    2727               
    2828        public:
     
    3838                int  end(void);
    3939                void Echo(void);
    40                 void GaussFromCoords(double x1,double y1,double* xyz_list);
     40                void GaussFromCoords(IssmDouble x1,IssmDouble y1,IssmDouble* xyz_list);
    4141                void GaussPoint(int ig);
    4242                void GaussVertex(int iv);
  • issm/trunk-jpl/src/c/objects/Loads/Friction.cpp

    r12365 r12424  
    5555}
    5656/*}}}*/
    57 /*FUNCTION Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
    58 void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     57/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/
     58void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
    5959
    6060        /*This routine calculates the basal friction coefficient
     
    6262
    6363        /*diverse: */
    64         double  r,s;
    65         double  drag_p, drag_q;
    66         double  gravity,rho_ice,rho_water;
    67         double  Neff;
    68         double  thickness,bed;
    69         double  vx,vy,vz,vmag;
    70         double  drag_coefficient;
    71         double  alpha2;
     64        IssmDouble  r,s;
     65        IssmDouble  drag_p, drag_q;
     66        IssmDouble  gravity,rho_ice,rho_water;
     67        IssmDouble  Neff;
     68        IssmDouble  thickness,bed;
     69        IssmDouble  vx,vy,vz,vmag;
     70        IssmDouble  drag_coefficient;
     71        IssmDouble  alpha2;
    7272
    7373        /*Recover parameters: */
     
    119119}
    120120/*}}}*/
    121 /*FUNCTION Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
    122 void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
     121/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/
     122void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    123123
    124124        /*This routine calculates the basal friction coefficient
     
    126126
    127127        /*diverse: */
    128         double  r,s;
    129         double  drag_p, drag_q;
    130         double  gravity,rho_ice,rho_water;
    131         double  Neff;
    132         double  thickness,bed;
    133         double  vx,vy,vz,vmag;
    134         double  drag_coefficient;
    135         double  alpha2;
     128        IssmDouble  r,s;
     129        IssmDouble  drag_p, drag_q;
     130        IssmDouble  gravity,rho_ice,rho_water;
     131        IssmDouble  Neff;
     132        IssmDouble  thickness,bed;
     133        IssmDouble  vx,vy,vz,vmag;
     134        IssmDouble  drag_coefficient;
     135        IssmDouble  alpha2;
    136136
    137137        /*Recover parameters: */
     
    183183}
    184184/*}}}*/
    185 /*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
    186 void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     185/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/
     186void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
    187187
    188188        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     
    192192        /*diverse: */
    193193        int     i;
    194         double  r,s;
    195         double  vx,vy,vz,vmag;
    196         double  drag_p,drag_q;
    197         double  Neff;
    198         double  drag_coefficient;
    199         double  bed,thickness;
    200         double  gravity,rho_ice,rho_water;
    201         double  alpha_complement;
     194        IssmDouble  r,s;
     195        IssmDouble  vx,vy,vz,vmag;
     196        IssmDouble  drag_p,drag_q;
     197        IssmDouble  Neff;
     198        IssmDouble  drag_coefficient;
     199        IssmDouble  bed,thickness;
     200        IssmDouble  gravity,rho_ice,rho_water;
     201        IssmDouble  alpha_complement;
    202202
    203203        /*Recover parameters: */
     
    249249}
    250250/*}}}*/
    251 /*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
    252 void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
     251/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/
     252void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    253253
    254254        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     
    258258        /*diverse: */
    259259        int     i;
    260         double  r,s;
    261         double  vx,vy,vz,vmag;
    262         double  drag_p,drag_q;
    263         double  Neff;
    264         double  drag_coefficient;
    265         double  bed,thickness;
    266         double  gravity,rho_ice,rho_water;
    267         double  alpha_complement;
     260        IssmDouble  r,s;
     261        IssmDouble  vx,vy,vz,vmag;
     262        IssmDouble  drag_p,drag_q;
     263        IssmDouble  Neff;
     264        IssmDouble  drag_coefficient;
     265        IssmDouble  bed,thickness;
     266        IssmDouble  gravity,rho_ice,rho_water;
     267        IssmDouble  alpha_complement;
    268268
    269269        /*Recover parameters: */
     
    316316/*}}}*/
    317317/*FUNCTION Friction::GetInputValue{{{*/
    318 void Friction::GetInputValue(double* pvalue,GaussTria* gauss,int enum_type){
     318void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){
    319319
    320320        Input* input=inputs->GetInput(enum_type);
     
    325325/*}}}*/
    326326/*FUNCTION Friction::GetInputValue{{{*/
    327 void Friction::GetInputValue(double* pvalue,GaussPenta* gauss,int enum_type){
     327void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){
    328328
    329329        Input* input=inputs->GetInput(enum_type);
  • issm/trunk-jpl/src/c/objects/Loads/Friction.h

    r12365 r12424  
    2727       
    2828                void  Echo(void);
    29                 void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    30                 void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    31                 void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    32                 void  GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    33                 void  GetInputValue(double* pvalue,GaussTria* gauss,int enum_type);
    34                 void  GetInputValue(double* pvalue,GaussPenta* gauss,int enum_type);
     29                void  GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
     30                void  GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
     31                void  GetAlphaComplement(IssmDouble* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
     32                void  GetAlphaComplement(IssmDouble* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
     33                void  GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type);
     34                void  GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type);
    3535
    3636};
  • issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp

    r12365 r12424  
    11/*  Gauss point structures and prototypes  */
    22
     3#include "../../include/include.h"
    34#include "./GaussPoints.h"
    45#include "../Alloc/alloc.h"
    5 #include "../../include/include.h"
    66#include "../../io/io.h"
    77#include "../Exceptions/exceptions.h"
     
    1111/*General Gauss points*/
    1212/*FUNCTION GaussLegendreLinear {{{*/
    13 void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus){
     13void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus){
    1414        /* Gauss-Legendre quadrature points.
    1515
     
    2525        /*Intermediaries*/
    2626        int i;
    27         double *alpha,*beta;
     27        IssmPDouble *alpha,*beta;
    2828
    2929        /*p= 1, npoint= 1*/
    30         static double wgt1[]={2.000000000000000};
    31         static double xi1[]={0.000000000000000};
     30        static IssmPDouble wgt1[]={2.000000000000000};
     31        static IssmPDouble xi1[]={0.000000000000000};
    3232
    3333        /*p= 3, npoint= 2*/
    34         static double wgt2[]={1.000000000000000, 1.000000000000000};
    35         static double xi2[]={-0.577350269189626, 0.577350269189626};
     34        static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
     35        static IssmPDouble xi2[]={-0.577350269189626, 0.577350269189626};
    3636
    3737        /*p= 5, npoint= 3*/
    38         static double wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};
    39         static double xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};
     38        static IssmPDouble wgt3[]={0.555555555555556, 0.888888888888889, 0.555555555555556};
     39        static IssmPDouble xi3[]={-0.774596669241483, 0.000000000000000, 0.774596669241483};
    4040
    4141        /*p= 7, npoint= 4*/
    42         static double wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
    43         static double xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};
    44 
    45         static double* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 };
    46         static double* xip [MAX_LINE_GAUS_PTS]={xi1  ,xi2  ,xi3  ,xi4  };
    47 
    48         static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof(double),
    49                 sizeof(wgt2 )/sizeof(double),
    50                 sizeof(wgt3 )/sizeof(double),
    51                 sizeof(wgt4 )/sizeof(double)};
     42        static IssmPDouble wgt4[]={0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454};
     43        static IssmPDouble xi4[]={-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053};
     44
     45        static IssmPDouble* wgtp[MAX_LINE_GAUS_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 };
     46        static IssmPDouble* xip [MAX_LINE_GAUS_PTS]={xi1  ,xi2  ,xi3  ,xi4  };
     47
     48        static int np[MAX_LINE_GAUS_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble),
     49                sizeof(wgt2 )/sizeof(IssmPDouble),
     50                sizeof(wgt3 )/sizeof(IssmPDouble),
     51                sizeof(wgt4 )/sizeof(IssmPDouble)};
    5252
    5353        //      _printf_(true,"Gauss-Legendre recurrence coefficients ngaus=%d\n",ngaus);
    54         *pxgaus = (double *) xmalloc(ngaus*sizeof(double));
    55         *pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
     54        *pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
     55        *pxwgt  = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
    5656
    5757        /*  check to see if Gauss points need to be calculated  */
     
    7070
    7171                /*  calculate the Gauss points using recurrence relations  */
    72                 alpha=(double *) xmalloc(ngaus*sizeof(double));
    73                 beta =(double *) xmalloc(ngaus*sizeof(double));
     72                alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
     73                beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
    7474
    7575                /*  calculate the Legendre recurrence coefficients  */
     
    8989}/*}}}*/
    9090/*FUNCTION GaussLegendreTria{{{*/
    91 void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord ) {
     91void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord ) {
    9292        /*Gauss quadrature points for the triangle.
    9393
     
    9898        /*Intermediaries*/
    9999        int i,j,ipt,nigaus;
    100         double xi,eta;
    101         double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
     100        IssmPDouble xi,eta;
     101        IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt;
    102102
    103103        /*Hardcoded Gauss points declaration*/
    104104        /*p= 1, npoint= 1{{{*/
    105         static double wgt1[]={
     105        static IssmPDouble wgt1[]={
    106106                1.732050807568877};
    107         static double l11[]={
     107        static IssmPDouble l11[]={
    108108                0.333333333333333};
    109         static double l21[]={
     109        static IssmPDouble l21[]={
    110110                0.333333333333333};
    111         static double l31[]={
     111        static IssmPDouble l31[]={
    112112                0.333333333333333};
    113113        /*}}}*/
    114114        /*p= 2, npoint= 3  {{{*/
    115         static double wgt2[]={
     115        static IssmPDouble wgt2[]={
    116116                0.577350269189625, 0.577350269189625, 0.577350269189625};
    117         static double l12[]={
     117        static IssmPDouble l12[]={
    118118                0.666666666666667, 0.166666666666667, 0.166666666666667};
    119         static double l22[]={
     119        static IssmPDouble l22[]={
    120120                0.166666666666667, 0.666666666666667, 0.166666666666667};
    121         static double l32[]={
     121        static IssmPDouble l32[]={
    122122                0.166666666666667, 0.166666666666667, 0.666666666666667};
    123123        /*}}}*/
    124124        /*p= 3, npoint= 4  {{{*/
    125         static double wgt3[]={
     125        static IssmPDouble wgt3[]={
    126126                -0.974278579257493, 0.902109795608790, 0.902109795608790,
    127127                0.902109795608790};
    128         static double l13[]={
     128        static IssmPDouble l13[]={
    129129                0.333333333333333, 0.600000000000000, 0.200000000000000,
    130130                0.200000000000000};
    131         static double l23[]={
     131        static IssmPDouble l23[]={
    132132                0.333333333333333, 0.200000000000000, 0.600000000000000,
    133133                0.200000000000000};
    134         static double l33[]={
     134        static IssmPDouble l33[]={
    135135                0.333333333333333, 0.200000000000000, 0.200000000000000,
    136136                0.600000000000000};
    137137        /*}}}*/
    138138        /*p= 4, npoint= 6  {{{*/
    139         static double wgt4[]={
     139        static IssmPDouble wgt4[]={
    140140                0.386908262797819, 0.386908262797819, 0.386908262797819,
    141141                0.190442006391807, 0.190442006391807, 0.190442006391807};
    142         static double l14[]={
     142        static IssmPDouble l14[]={
    143143                0.108103018168070, 0.445948490915965, 0.445948490915965,
    144144                0.816847572980459, 0.091576213509771, 0.091576213509771};
    145         static double l24[]={
     145        static IssmPDouble l24[]={
    146146                0.445948490915965, 0.108103018168070, 0.445948490915965,
    147147                0.091576213509771, 0.816847572980459, 0.091576213509771};
    148         static double l34[]={
     148        static IssmPDouble l34[]={
    149149                0.445948490915965, 0.445948490915965, 0.108103018168070,
    150150                0.091576213509771, 0.091576213509771, 0.816847572980459};
    151151        /*}}}*/
    152152        /*p= 5, npoint= 7  {{{*/
    153         static double wgt5[]={
     153        static IssmPDouble wgt5[]={
    154154                0.389711431702997, 0.229313399254729, 0.229313399254729,
    155155                0.229313399254729, 0.218133059367230, 0.218133059367230,
    156156                0.218133059367230};
    157         static double l15[]={
     157        static IssmPDouble l15[]={
    158158                0.333333333333333, 0.059715871789770, 0.470142064105115,
    159159                0.470142064105115, 0.797426985353087, 0.101286507323456,
    160160                0.101286507323456};
    161         static double l25[]={
     161        static IssmPDouble l25[]={
    162162                0.333333333333333, 0.470142064105115, 0.059715871789770,
    163163                0.470142064105115, 0.101286507323456, 0.797426985353087,
    164164                0.101286507323456};
    165         static double l35[]={
     165        static IssmPDouble l35[]={
    166166                0.333333333333333, 0.470142064105115, 0.470142064105115,
    167167                0.059715871789770, 0.101286507323456, 0.101286507323456,
     
    169169        /*}}}*/
    170170        /*p= 6, npoint=12  {{{*/
    171         static double wgt6[]={
     171        static IssmPDouble wgt6[]={
    172172                0.202279763184836, 0.202279763184836, 0.202279763184836,
    173173                0.088065961139281, 0.088065961139281, 0.088065961139281,
    174174                0.143502272432755, 0.143502272432755, 0.143502272432755,
    175175                0.143502272432755, 0.143502272432755, 0.143502272432755};
    176         static double l16[]={
     176        static IssmPDouble l16[]={
    177177                0.501426509658179, 0.249286745170910, 0.249286745170910,
    178178                0.873821971016996, 0.063089014491502, 0.063089014491502,
    179179                0.053145049844817, 0.053145049844817, 0.310352451033784,
    180180                0.636502499121399, 0.310352451033784, 0.636502499121399};
    181         static double l26[]={
     181        static IssmPDouble l26[]={
    182182                0.249286745170910, 0.501426509658179, 0.249286745170910,
    183183                0.063089014491502, 0.873821971016996, 0.063089014491502,
    184184                0.310352451033784, 0.636502499121399, 0.053145049844817,
    185185                0.053145049844817, 0.636502499121399, 0.310352451033784};
    186         static double l36[]={
     186        static IssmPDouble l36[]={
    187187                0.249286745170910, 0.249286745170910, 0.501426509658179,
    188188                0.063089014491502, 0.063089014491502, 0.873821971016996,
     
    191191        /*}}}*/
    192192        /*p= 7, npoint=13  {{{*/
    193         static double wgt7[]={
     193        static IssmPDouble wgt7[]={
    194194                -0.259062916308362, 0.304174548458604, 0.304174548458604,
    195195                0.304174548458604, 0.092400122517855, 0.092400122517855,
     
    197197                0.133564951824643, 0.133564951824643, 0.133564951824643,
    198198                0.133564951824643};
    199         static double l17[]={
     199        static IssmPDouble l17[]={
    200200                0.333333333333333, 0.479308067841920, 0.260345966079040,
    201201                0.260345966079040, 0.869739794195568, 0.065130102902216,
     
    203203                0.312865496004874, 0.638444188569810, 0.312865496004874,
    204204                0.638444188569810};
    205         static double l27[]={
     205        static IssmPDouble l27[]={
    206206                0.333333333333333, 0.260345966079040, 0.479308067841920,
    207207                0.260345966079040, 0.065130102902216, 0.869739794195568,
     
    209209                0.048690315425316, 0.048690315425316, 0.638444188569810,
    210210                0.312865496004874};
    211         static double l37[]={
     211        static IssmPDouble l37[]={
    212212                0.333333333333333, 0.260345966079040, 0.260345966079040,
    213213                0.479308067841920, 0.065130102902216, 0.065130102902216,
     
    217217        /*}}}*/
    218218        /*p= 8, npoint=16  {{{*/
    219         static double wgt8[]={
     219        static IssmPDouble wgt8[]={
    220220                0.249961964823104, 0.164703541925695, 0.164703541925695,
    221221                0.164703541925695, 0.178777729989794, 0.178777729989794,
     
    224224                0.047164287656184, 0.047164287656184, 0.047164287656184,
    225225                0.047164287656184};
    226         static double l18[]={
     226        static IssmPDouble l18[]={
    227227                0.333333333333333, 0.081414823414554, 0.459292588292723,
    228228                0.459292588292723, 0.658861384496480, 0.170569307751760,
     
    231231                0.263112829634638, 0.728492392955404, 0.263112829634638,
    232232                0.728492392955404};
    233         static double l28[]={
     233        static IssmPDouble l28[]={
    234234                0.333333333333333, 0.459292588292723, 0.081414823414554,
    235235                0.459292588292723, 0.170569307751760, 0.658861384496480,
     
    238238                0.008394777409958, 0.008394777409958, 0.728492392955404,
    239239                0.263112829634638};
    240         static double l38[]={
     240        static IssmPDouble l38[]={
    241241                0.333333333333333, 0.459292588292723, 0.459292588292723,
    242242                0.081414823414554, 0.170569307751760, 0.170569307751760,
     
    247247        /*}}}*/
    248248        /*p= 9, npoint=19  {{{*/
    249         static double wgt9[]={
     249        static IssmPDouble wgt9[]={
    250250                0.168244134395468, 0.054273292833345, 0.054273292833345,
    251251                0.054273292833345, 0.134801255248419, 0.134801255248419,
     
    255255                0.074969289332873, 0.074969289332873, 0.074969289332873,
    256256                0.074969289332873};
    257         static double l19[]={
     257        static IssmPDouble l19[]={
    258258                0.333333333333333, 0.020634961602525, 0.489682519198738,
    259259                0.489682519198738, 0.125820817014127, 0.437089591492937,
     
    263263                0.221962989160766, 0.741198598784498, 0.221962989160766,
    264264                0.741198598784498};
    265         static double l29[]={
     265        static IssmPDouble l29[]={
    266266                0.333333333333333, 0.489682519198738, 0.020634961602525,
    267267                0.489682519198738, 0.437089591492937, 0.125820817014127,
     
    271271                0.036838412054736, 0.036838412054736, 0.741198598784498,
    272272                0.221962989160766};
    273         static double l39[]={
     273        static IssmPDouble l39[]={
    274274                0.333333333333333, 0.489682519198738, 0.489682519198738,
    275275                0.020634961602525, 0.437089591492937, 0.437089591492937,
     
    281281        /*}}}*/
    282282        /*p=10, npoint=25  {{{*/
    283         static double wgt10[]={
     283        static IssmPDouble wgt10[]={
    284284                0.157301373584232, 0.063611224790829, 0.063611224790829,
    285285                0.063611224790829, 0.078498377595183, 0.078498377595183,
     
    291291                0.016318805873179, 0.016318805873179, 0.016318805873179,
    292292                0.016318805873179};
    293         static double l110[]={
     293        static IssmPDouble l110[]={
    294294                0.333333333333333, 0.028844733232685, 0.485577633383657,
    295295                0.485577633383657, 0.781036849029926, 0.109481575485037,
     
    301301                0.066803251012200, 0.923655933587500, 0.066803251012200,
    302302                0.923655933587500};
    303         static double l210[]={
     303        static IssmPDouble l210[]={
    304304                0.333333333333333, 0.485577633383657, 0.028844733232685,
    305305                0.485577633383657, 0.109481575485037, 0.781036849029926,
     
    311311                0.009540815400299, 0.009540815400299, 0.923655933587500,
    312312                0.066803251012200};
    313         static double l310[]={
     313        static IssmPDouble l310[]={
    314314                0.333333333333333, 0.485577633383657, 0.485577633383657,
    315315                0.028844733232685, 0.109481575485037, 0.109481575485037,
     
    323323        /*}}}*/
    324324        /*p=11, npoint=27  {{{*/
    325         static double wgt11[]={
     325        static IssmPDouble wgt11[]={
    326326                0.001605622060698, 0.001605622060698, 0.001605622060698,
    327327                0.133626914252765, 0.133626914252765, 0.133626914252765,
     
    333333                0.035866718600836, 0.035866718600836, 0.035866718600836,
    334334                0.035866718600836, 0.035866718600836, 0.035866718600836};
    335         static double l111[]={
     335        static IssmPDouble l111[]={
    336336                -0.069222096541517, 0.534611048270758, 0.534611048270758,
    337337                0.202061394068290, 0.398969302965855, 0.398969302965855,
     
    343343                0.021022016536166, 0.021022016536166, 0.171488980304042,
    344344                0.807489003159792, 0.171488980304042, 0.807489003159792};
    345         static double l211[]={
     345        static IssmPDouble l211[]={
    346346                0.534611048270758,-0.069222096541517, 0.534611048270758,
    347347                0.398969302965855, 0.202061394068290, 0.398969302965855,
     
    353353                0.171488980304042, 0.807489003159792, 0.021022016536166,
    354354                0.021022016536166, 0.807489003159792, 0.171488980304042};
    355         static double l311[]={
     355        static IssmPDouble l311[]={
    356356                0.534611048270758, 0.534611048270758,-0.069222096541517,
    357357                0.398969302965855, 0.398969302965855, 0.202061394068290,
     
    365365        /*}}}*/
    366366        /*p=12, npoint=33  {{{*/
    367         static double wgt12[]={
     367        static IssmPDouble wgt12[]={
    368368                0.044567514407799, 0.044567514407799, 0.044567514407799,
    369369                0.075677707051848, 0.075677707051848, 0.075677707051848,
     
    377377                0.029992592075802, 0.029992592075802, 0.029992592075802,
    378378                0.029992592075802, 0.029992592075802, 0.029992592075802};
    379         static double l112[]={
     379        static IssmPDouble l112[]={
    380380                0.023565220452390, 0.488217389773805, 0.488217389773805,
    381381                0.120551215411079, 0.439724392294460, 0.439724392294460,
     
    389389                0.025734050548330, 0.025734050548330, 0.116251915907597,
    390390                0.858014033544073, 0.116251915907597, 0.858014033544073};
    391         static double l212[]={
     391        static IssmPDouble l212[]={
    392392                0.488217389773805, 0.023565220452390, 0.488217389773805,
    393393                0.439724392294460, 0.120551215411079, 0.439724392294460,
     
    401401                0.116251915907597, 0.858014033544073, 0.025734050548330,
    402402                0.025734050548330, 0.858014033544073, 0.116251915907597};
    403         static double l312[]={
     403        static IssmPDouble l312[]={
    404404                0.488217389773805, 0.488217389773805, 0.023565220452390,
    405405                0.439724392294460, 0.439724392294460, 0.120551215411079,
     
    415415        /*}}}*/
    416416        /*  p=13, npoint=37  {{{*/
    417         static double wgt13[]={
     417        static IssmPDouble wgt13[]={
    418418                0.090968907790622, 0.019537784619314, 0.019537784619314,
    419419                0.019537784619314, 0.054427130356344, 0.054427130356344,
     
    429429                0.026884523429480, 0.026884523429480, 0.026884523429480,
    430430                0.026884523429480};
    431         static double l113[]={
     431        static IssmPDouble l113[]={
    432432                0.333333333333333, 0.009903630120591, 0.495048184939705,
    433433                0.495048184939705, 0.062566729780852, 0.468716635109574,
     
    443443                0.126357385491669, 0.851409537834241, 0.126357385491669,
    444444                0.851409537834241};
    445         static double l213[]={
     445        static IssmPDouble l213[]={
    446446                0.333333333333333, 0.495048184939705, 0.009903630120591,
    447447                0.495048184939705, 0.468716635109574, 0.062566729780852,
     
    457457                0.022233076674090, 0.022233076674090, 0.851409537834241,
    458458                0.126357385491669};
    459         static double l313[]={
     459        static IssmPDouble l313[]={
    460460                0.333333333333333, 0.495048184939705, 0.495048184939705,
    461461                0.009903630120591, 0.468716635109574, 0.468716635109574,
     
    473473        /*}}}*/
    474474        /*p=14, npoint=42{{{*/
    475         static double wgt14[]={
     475        static IssmPDouble wgt14[]={
    476476                0.037903474783419, 0.037903474783419, 0.037903474783419,
    477477                0.056791094234956, 0.056791094234956, 0.056791094234956,
     
    488488                0.008677970905831, 0.008677970905831, 0.008677970905831,
    489489                0.008677970905831, 0.008677970905831, 0.008677970905831};
    490         static double l114[]={
     490        static IssmPDouble l114[]={
    491491                0.022072179275643, 0.488963910362179, 0.488963910362179,
    492492                0.164710561319092, 0.417644719340454, 0.417644719340454,
     
    503503                0.001268330932872, 0.001268330932872, 0.118974497696957,
    504504                0.879757171370171, 0.118974497696957, 0.879757171370171};
    505         static double l214[]={
     505        static IssmPDouble l214[]={
    506506                0.488963910362179, 0.022072179275643, 0.488963910362179,
    507507                0.417644719340454, 0.164710561319092, 0.417644719340454,
     
    518518                0.118974497696957, 0.879757171370171, 0.001268330932872,
    519519                0.001268330932872, 0.879757171370171, 0.118974497696957};
    520         static double l314[]={
     520        static IssmPDouble l314[]={
    521521                0.488963910362179, 0.488963910362179, 0.022072179275643,
    522522                0.417644719340454, 0.417644719340454, 0.164710561319092,
     
    535535        /*}}}*/
    536536        /*p=15, npoint=48{{{*/
    537         static double wgt15[]={
     537        static IssmPDouble wgt15[]={
    538538                0.003320126005206, 0.003320126005206, 0.003320126005206,
    539539                0.076641563419124, 0.076641563419124, 0.076641563419124,
     
    552552                0.013291658531346, 0.013291658531346, 0.013291658531346,
    553553                0.013291658531346, 0.013291658531346, 0.013291658531346};
    554         static double l115[]={
     554        static IssmPDouble l115[]={
    555555                -0.013945833716486, 0.506972916858243, 0.506972916858243,
    556556                0.137187291433955, 0.431406354283023, 0.431406354283023,
     
    569569                0.012459809331199, 0.012459809331199, 0.103575616576386,
    570570                0.883964574092416, 0.103575616576386, 0.883964574092416};
    571         static double l215[]={
     571        static IssmPDouble l215[]={
    572572                0.506972916858243,-0.013945833716486, 0.506972916858243,
    573573                0.431406354283023, 0.137187291433955, 0.431406354283023,
     
    586586                0.103575616576386, 0.883964574092416, 0.012459809331199,
    587587                0.012459809331199, 0.883964574092416, 0.103575616576386};
    588         static double l315[]={
     588        static IssmPDouble l315[]={
    589589                0.506972916858243, 0.506972916858243,-0.013945833716486,
    590590                0.431406354283023, 0.431406354283023, 0.137187291433955,
     
    605605        /*}}}*/
    606606        /*p=16, npoint=52  {{{*/
    607         static double wgt16[]={
     607        static IssmPDouble wgt16[]={
    608608                0.081191089584902, 0.011095307165226, 0.011095307165226,
    609609                0.011095307165226, 0.072244353151393, 0.072244353151393,
     
    624624                0.011864642509229, 0.011864642509229, 0.011864642509229,
    625625                0.011864642509229};
    626         static double l116[]={
     626        static IssmPDouble l116[]={
    627627                0.333333333333333, 0.005238916103123, 0.497380541948438,
    628628                0.497380541948438, 0.173061122901295, 0.413469438549352,
     
    643643                0.085283615682657, 0.900399064086661, 0.085283615682657,
    644644                0.900399064086661};
    645         static double l216[]={
     645        static IssmPDouble l216[]={
    646646                0.333333333333333, 0.497380541948438, 0.005238916103123,
    647647                0.497380541948438, 0.413469438549352, 0.173061122901295,
     
    662662                0.014317320230681, 0.014317320230681, 0.900399064086661,
    663663                0.085283615682657};
    664         static double l316[]={
     664        static IssmPDouble l316[]={
    665665                0.333333333333333, 0.497380541948438, 0.497380541948438,
    666666                0.005238916103123, 0.413469438549352, 0.413469438549352,
     
    683683        /*}}}*/
    684684        /*p=17, npoint=61{{{*/
    685         static double wgt17[]={
     685        static IssmPDouble wgt17[]={
    686686                0.057914928034477, 0.008822054327014, 0.008822054327014,
    687687                0.008822054327014, 0.025410682752829, 0.025410682752829,
     
    705705                0.011545213295771, 0.011545213295771, 0.011545213295771,
    706706                0.011545213295771};
    707         static double l117[]={
     707        static IssmPDouble l117[]={
    708708                0.333333333333333, 0.005658918886452, 0.497170540556774,
    709709                0.497170540556774, 0.035647354750751, 0.482176322624625,
     
    727727                0.080711313679564, 0.904625504095608, 0.080711313679564,
    728728                0.904625504095608};
    729         static double l217[]={
     729        static IssmPDouble l217[]={
    730730                0.333333333333333, 0.497170540556774, 0.005658918886452,
    731731                0.497170540556774, 0.482176322624625, 0.035647354750751,
     
    749749                0.014663182224828, 0.014663182224828, 0.904625504095608,
    750750                0.080711313679564};
    751         static double l317[]={
     751        static IssmPDouble l317[]={
    752752                0.333333333333333, 0.497170540556774, 0.497170540556774,
    753753                0.005658918886452, 0.482176322624625, 0.482176322624625,
     
    774774        /*  p=18, npoint=70  {{{*/
    775775
    776         static double wgt18[]={
     776        static IssmPDouble wgt18[]={
    777777                0.053364381350150, 0.015713921277179, 0.015713921277179,
    778778                0.015713921277179, 0.032495554156279, 0.032495554156279,
     
    799799                0.000079999375178, 0.000079999375178, 0.000079999375178,
    800800                0.000079999375178};
    801         static double l118[]={
     801        static IssmPDouble l118[]={
    802802                0.333333333333333, 0.013310382738157, 0.493344808630921,
    803803                0.493344808630921, 0.061578811516086, 0.469210594241957,
     
    824824                0.020874755282586, 1.014347260005363, 0.020874755282586,
    825825                1.014347260005363};
    826         static double l218[]={
     826        static IssmPDouble l218[]={
    827827                0.333333333333333, 0.493344808630921, 0.013310382738157,
    828828                0.493344808630921, 0.469210594241957, 0.061578811516086,
     
    849849                -0.035222015287949,-0.035222015287949, 1.014347260005363,
    850850                0.020874755282586};
    851         static double l318[]={
     851        static IssmPDouble l318[]={
    852852                0.333333333333333, 0.493344808630921, 0.493344808630921,
    853853                0.013310382738157, 0.469210594241957, 0.469210594241957,
     
    877877        /*p=19, npoint=73  {{{*/
    878878
    879         static double wgt19[]={
     879        static IssmPDouble wgt19[]={
    880880                0.056995437856306, 0.017893352515055, 0.017893352515055,
    881881                0.017893352515055, 0.038775849701151, 0.038775849701151,
     
    903903                0.006581669842530, 0.006581669842530, 0.006581669842530,
    904904                0.006581669842530};
    905         static double l119[]={
     905        static IssmPDouble l119[]={
    906906                0.333333333333333, 0.020780025853987, 0.489609987073006,
    907907                0.489609987073006, 0.090926214604215, 0.454536892697893,
     
    929929                0.065494628082938, 0.924344252620784, 0.065494628082938,
    930930                0.924344252620784};
    931         static double l219[]={
     931        static IssmPDouble l219[]={
    932932                0.333333333333333, 0.489609987073006, 0.020780025853987,
    933933                0.489609987073006, 0.454536892697893, 0.090926214604215,
     
    955955                0.010161119296278, 0.010161119296278, 0.924344252620784,
    956956                0.065494628082938};
    957         static double l319[]={
     957        static IssmPDouble l319[]={
    958958                0.333333333333333, 0.489609987073006, 0.489609987073006,
    959959                0.020780025853987, 0.454536892697893, 0.454536892697893,
     
    983983        /*}}}*/
    984984        /*p=20, npoint=79 {{{*/
    985         static double wgt20[]={
     985        static IssmPDouble wgt20[]={
    986986                0.057256499746719, 0.001501721280705, 0.001501721280705,
    987987                0.001501721280705, 0.020195803723819, 0.020195803723819,
     
    10111011                0.006190192638113, 0.006190192638113, 0.006190192638113,
    10121012                0.006190192638113};
    1013         static double l120[]={
     1013        static IssmPDouble l120[]={
    10141014                0.333333333333333,-0.001900928704400, 0.500950464352200,
    10151015                0.500950464352200, 0.023574084130543, 0.488212957934729,
     
    10391039                0.059696109149007, 0.929756171556853, 0.059696109149007,
    10401040                0.929756171556853};
    1041         static double l220[]={
     1041        static IssmPDouble l220[]={
    10421042                0.333333333333333, 0.500950464352200,-0.001900928704400,
    10431043                0.500950464352200, 0.488212957934729, 0.023574084130543,
     
    10671067                0.010547719294141, 0.010547719294141, 0.929756171556853,
    10681068                0.059696109149007};
    1069         static double l320[]={
     1069        static IssmPDouble l320[]={
    10701070                0.333333333333333, 0.500950464352200, 0.500950464352200,
    10711071                -0.001900928704400, 0.488212957934729, 0.488212957934729,
     
    10971097        /*}}}*/
    10981098
    1099         static double* wgtp[MAX_TRIA_SYM_ORD]={
     1099        static IssmPDouble* wgtp[MAX_TRIA_SYM_ORD]={
    11001100                wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 ,
    11011101                wgt6 ,wgt7 ,wgt8 ,wgt9 ,wgt10,
    11021102                wgt11,wgt12,wgt13,wgt14,wgt15,
    11031103                wgt16,wgt17,wgt18,wgt19,wgt20};
    1104         static double* l1p [MAX_TRIA_SYM_ORD]={
     1104        static IssmPDouble* l1p [MAX_TRIA_SYM_ORD]={
    11051105                l11  ,l12  ,l13  ,l14  ,l15  ,
    11061106                l16  ,l17  ,l18  ,l19  ,l110 ,
    11071107                l111 ,l112 ,l113 ,l114 ,l115 ,
    11081108                l116 ,l117 ,l118 ,l119 ,l120 };
    1109         static double* l2p [MAX_TRIA_SYM_ORD]={
     1109        static IssmPDouble* l2p [MAX_TRIA_SYM_ORD]={
    11101110                l21  ,l22  ,l23  ,l24  ,l25  ,
    11111111                l26  ,l27  ,l28  ,l29  ,l210 ,
    11121112                l211 ,l212 ,l213 ,l214 ,l215 ,
    11131113                l216 ,l217 ,l218 ,l219 ,l220 };
    1114         static double* l3p [MAX_TRIA_SYM_ORD]={
     1114        static IssmPDouble* l3p [MAX_TRIA_SYM_ORD]={
    11151115                l31  ,l32  ,l33  ,l34  ,l35  ,
    11161116                l36  ,l37  ,l38  ,l39  ,l310 ,
     
    11181118                l316 ,l317 ,l318 ,l319 ,l320 };
    11191119
    1120         static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(double),
    1121                 sizeof(wgt2 )/sizeof(double),
    1122                 sizeof(wgt3 )/sizeof(double),
    1123                 sizeof(wgt4 )/sizeof(double),
    1124                 sizeof(wgt5 )/sizeof(double),
    1125                 sizeof(wgt6 )/sizeof(double),
    1126                 sizeof(wgt7 )/sizeof(double),
    1127                 sizeof(wgt8 )/sizeof(double),
    1128                 sizeof(wgt9 )/sizeof(double),
    1129                 sizeof(wgt10)/sizeof(double),
    1130                 sizeof(wgt11)/sizeof(double),
    1131                 sizeof(wgt12)/sizeof(double),
    1132                 sizeof(wgt13)/sizeof(double),
    1133                 sizeof(wgt14)/sizeof(double),
    1134                 sizeof(wgt15)/sizeof(double),
    1135                 sizeof(wgt16)/sizeof(double),
    1136                 sizeof(wgt17)/sizeof(double),
    1137                 sizeof(wgt18)/sizeof(double),
    1138                 sizeof(wgt19)/sizeof(double),
    1139                 sizeof(wgt20)/sizeof(double)};
     1120        static int np[MAX_TRIA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
     1121                sizeof(wgt2 )/sizeof(IssmPDouble),
     1122                sizeof(wgt3 )/sizeof(IssmPDouble),
     1123                sizeof(wgt4 )/sizeof(IssmPDouble),
     1124                sizeof(wgt5 )/sizeof(IssmPDouble),
     1125                sizeof(wgt6 )/sizeof(IssmPDouble),
     1126                sizeof(wgt7 )/sizeof(IssmPDouble),
     1127                sizeof(wgt8 )/sizeof(IssmPDouble),
     1128                sizeof(wgt9 )/sizeof(IssmPDouble),
     1129                sizeof(wgt10)/sizeof(IssmPDouble),
     1130                sizeof(wgt11)/sizeof(IssmPDouble),
     1131                sizeof(wgt12)/sizeof(IssmPDouble),
     1132                sizeof(wgt13)/sizeof(IssmPDouble),
     1133                sizeof(wgt14)/sizeof(IssmPDouble),
     1134                sizeof(wgt15)/sizeof(IssmPDouble),
     1135                sizeof(wgt16)/sizeof(IssmPDouble),
     1136                sizeof(wgt17)/sizeof(IssmPDouble),
     1137                sizeof(wgt18)/sizeof(IssmPDouble),
     1138                sizeof(wgt19)/sizeof(IssmPDouble),
     1139                sizeof(wgt20)/sizeof(IssmPDouble)};
    11401140
    11411141        //      _printf_(true,"GaussLegendreTria: iord=%d\n",iord);
     
    11501150                *pngaus=np[iord-1];
    11511151
    1152                 *pl1  = (double *) xmalloc(*pngaus*sizeof(double));
    1153                 *pl2  = (double *) xmalloc(*pngaus*sizeof(double));
    1154                 *pl3  = (double *) xmalloc(*pngaus*sizeof(double));
    1155                 *pwgt = (double *) xmalloc(*pngaus*sizeof(double));
     1152                *pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1153                *pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1154                *pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1155                *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
    11561156
    11571157                for (i=0; i<*pngaus; i++) {
     
    11681168                *pngaus=nigaus*nigaus;
    11691169
    1170                 *pl1  = (double *) xmalloc(*pngaus*sizeof(double));
    1171                 *pl2  = (double *) xmalloc(*pngaus*sizeof(double));
    1172                 *pl3  = (double *) xmalloc(*pngaus*sizeof(double));
    1173                 *pwgt = (double *) xmalloc(*pngaus*sizeof(double));
     1170                *pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1171                *pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1172                *pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1173                *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
    11741174
    11751175                /*  get the gauss points in each direction  */
     
    12071207}/*}}}*/
    12081208/*FUNCTION GaussLegendreTetra{{{*/
    1209 void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord ) {
     1209void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord ) {
    12101210        /* Gauss quadrature points for the tetrahedron.
    12111211
     
    12201220        /*Intermediaries*/
    12211221        int i,j,k,ipt,nigaus;
    1222         double xi,eta,zeta;
    1223         double *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
     1222        IssmPDouble xi,eta,zeta;
     1223        IssmPDouble *xgaus=NULL,*xwgt=NULL,*egaus,*ewgt,*zgaus,*zwgt;
    12241224
    12251225        /*Hardcoded Gauss points definition*/
    12261226        /*p= 1, npoint= 1  {{{*/
    1227         static double wgt1[]={
     1227        static IssmPDouble wgt1[]={
    12281228                1.000000000000000};
    1229         static double l11[]={
     1229        static IssmPDouble l11[]={
    12301230                0.250000000000000};
    1231         static double l21[]={
     1231        static IssmPDouble l21[]={
    12321232                0.250000000000000};
    1233         static double l31[]={
     1233        static IssmPDouble l31[]={
    12341234                0.250000000000000};
    1235         static double l41[]={
     1235        static IssmPDouble l41[]={
    12361236                0.250000000000000};
    12371237        /*}}}*/
    12381238        /*p= 2, npoint= 4  {{{*/
    12391239
    1240         static double wgt2[]={
     1240        static IssmPDouble wgt2[]={
    12411241                0.250000000000000, 0.250000000000000, 0.250000000000000,
    12421242                0.250000000000000};
    1243         static double l12[]={
     1243        static IssmPDouble l12[]={
    12441244                0.585410196624969, 0.138196601125011, 0.138196601125011,
    12451245                0.138196601125011};
    1246         static double l22[]={
     1246        static IssmPDouble l22[]={
    12471247                0.138196601125011, 0.585410196624969, 0.138196601125011,
    12481248                0.138196601125011};
    1249         static double l32[]={
     1249        static IssmPDouble l32[]={
    12501250                0.138196601125011, 0.138196601125011, 0.585410196624969,
    12511251                0.138196601125011};
    1252         static double l42[]={
     1252        static IssmPDouble l42[]={
    12531253                0.138196601125011, 0.138196601125011, 0.138196601125011,
    12541254                0.585410196624969};
    12551255        /*}}}*/
    12561256        /*p= 3, npoint= 5  {{{*/
    1257         static double wgt3[]={
     1257        static IssmPDouble wgt3[]={
    12581258                -0.800000000000000, 0.450000000000000, 0.450000000000000,
    12591259                0.450000000000000, 0.450000000000000};
    1260         static double l13[]={
     1260        static IssmPDouble l13[]={
    12611261                0.250000000000000, 0.500000000000000, 0.166666666666667,
    12621262                0.166666666666667, 0.166666666666667};
    1263         static double l23[]={
     1263        static IssmPDouble l23[]={
    12641264                0.250000000000000, 0.166666666666667, 0.500000000000000,
    12651265                0.166666666666667, 0.166666666666667};
    1266         static double l33[]={
     1266        static IssmPDouble l33[]={
    12671267                0.250000000000000, 0.166666666666667, 0.166666666666667,
    12681268                0.500000000000000, 0.166666666666667};
    1269         static double l43[]={
     1269        static IssmPDouble l43[]={
    12701270                0.250000000000000, 0.166666666666667, 0.166666666666667,
    12711271                0.166666666666667, 0.500000000000000};
     
    12731273        /*p= 4, npoint=11  {{{*/
    12741274
    1275         static double wgt4[]={
     1275        static IssmPDouble wgt4[]={
    12761276                -0.013155555555556, 0.007622222222222, 0.007622222222222,
    12771277                0.007622222222222, 0.007622222222222, 0.024888888888889,
    12781278                0.024888888888889, 0.024888888888889, 0.024888888888889,
    12791279                0.024888888888889, 0.024888888888889};
    1280         static double l14[]={
     1280        static IssmPDouble l14[]={
    12811281                0.250000000000000, 0.785714285714286, 0.071428571428571,
    12821282                0.071428571428571, 0.071428571428571, 0.399403576166799,
    12831283                0.399403576166799, 0.399403576166799, 0.100596423833201,
    12841284                0.100596423833201, 0.100596423833201};
    1285         static double l24[]={
     1285        static IssmPDouble l24[]={
    12861286                0.250000000000000, 0.071428571428571, 0.785714285714286,
    12871287                0.071428571428571, 0.071428571428571, 0.399403576166799,
    12881288                0.100596423833201, 0.100596423833201, 0.399403576166799,
    12891289                0.399403576166799, 0.100596423833201};
    1290         static double l34[]={
     1290        static IssmPDouble l34[]={
    12911291                0.250000000000000, 0.071428571428571, 0.071428571428571,
    12921292                0.785714285714286, 0.071428571428571, 0.100596423833201,
    12931293                0.399403576166799, 0.100596423833201, 0.399403576166799,
    12941294                0.100596423833201, 0.399403576166799};
    1295         static double l44[]={
     1295        static IssmPDouble l44[]={
    12961296                0.250000000000000, 0.071428571428571, 0.071428571428571,
    12971297                0.071428571428571, 0.785714285714286, 0.100596423833201,
     
    13011301        /*p= 5, npoint=15  {{{*/
    13021302
    1303         static double wgt5[]={
     1303        static IssmPDouble wgt5[]={
    13041304                0.030283678097089, 0.006026785714286, 0.006026785714286,
    13051305                0.006026785714286, 0.006026785714286, 0.011645249086029,
     
    13071307                0.010949141561386, 0.010949141561386, 0.010949141561386,
    13081308                0.010949141561386, 0.010949141561386, 0.010949141561386};
    1309         static double l15[]={
     1309        static IssmPDouble l15[]={
    13101310                0.250000000000000, 0.000000000000000, 0.333333333333333,
    13111311                0.333333333333333, 0.333333333333333, 0.727272727272727,
     
    13131313                0.066550153573664, 0.066550153573664, 0.066550153573664,
    13141314                0.433449846426336, 0.433449846426336, 0.433449846426336};
    1315         static double l25[]={
     1315        static IssmPDouble l25[]={
    13161316                0.250000000000000, 0.333333333333333, 0.000000000000000,
    13171317                0.333333333333333, 0.333333333333333, 0.090909090909091,
     
    13191319                0.066550153573664, 0.433449846426336, 0.433449846426336,
    13201320                0.066550153573664, 0.066550153573664, 0.433449846426336};
    1321         static double l35[]={
     1321        static IssmPDouble l35[]={
    13221322                0.250000000000000, 0.333333333333333, 0.333333333333333,
    13231323                0.000000000000000, 0.333333333333333, 0.090909090909091,
     
    13251325                0.433449846426336, 0.066550153573664, 0.433449846426336,
    13261326                0.066550153573664, 0.433449846426336, 0.066550153573664};
    1327         static double l45[]={
     1327        static IssmPDouble l45[]={
    13281328                0.250000000000000, 0.333333333333333, 0.333333333333333,
    13291329                0.333333333333333, 0.000000000000000, 0.090909090909091,
     
    13341334        /*p= 6, npoint=24  {{{*/
    13351335
    1336         static double wgt6[]={
     1336        static IssmPDouble wgt6[]={
    13371337                0.006653791709695, 0.006653791709695, 0.006653791709695,
    13381338                0.006653791709695, 0.001679535175887, 0.001679535175887,
     
    13431343                0.008035714285714, 0.008035714285714, 0.008035714285714,
    13441344                0.008035714285714, 0.008035714285714, 0.008035714285714};
    1345         static double l16[]={
     1345        static IssmPDouble l16[]={
    13461346                0.356191386222545, 0.214602871259152, 0.214602871259152,
    13471347                0.214602871259152, 0.877978124396166, 0.040673958534611,
     
    13531353                0.269672331458316, 0.603005664791649, 0.269672331458316,
    13541354                0.603005664791649, 0.269672331458316, 0.603005664791649};
    1355         static double l26[]={
     1355        static IssmPDouble l26[]={
    13561356                0.214602871259152, 0.356191386222545, 0.214602871259152,
    13571357                0.214602871259152, 0.040673958534611, 0.877978124396166,
     
    13631363                0.063661001875018, 0.063661001875018, 0.063661001875018,
    13641364                0.063661001875018, 0.603005664791649, 0.269672331458316};
    1365         static double l36[]={
     1365        static IssmPDouble l36[]={
    13661366                0.214602871259152, 0.214602871259152, 0.356191386222545,
    13671367                0.214602871259152, 0.040673958534611, 0.040673958534611,
     
    13731373                0.063661001875018, 0.063661001875018, 0.603005664791649,
    13741374                0.269672331458316, 0.063661001875018, 0.063661001875018};
    1375         static double l46[]={
     1375        static IssmPDouble l46[]={
    13761376                0.214602871259152, 0.214602871259152, 0.214602871259152,
    13771377                0.356191386222545, 0.040673958534611, 0.040673958534611,
     
    13851385        /*}}}*/
    13861386
    1387         static double* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};
    1388         static double* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };
    1389         static double* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 };
    1390         static double* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };
    1391         static double* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };
    1392 
    1393         static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(double),
    1394                 sizeof(wgt2 )/sizeof(double),
    1395                 sizeof(wgt3 )/sizeof(double),
    1396                 sizeof(wgt4 )/sizeof(double),
    1397                 sizeof(wgt5 )/sizeof(double),
    1398                 sizeof(wgt6 )/sizeof(double)};
     1387        static IssmPDouble* wgtp[MAX_TETRA_SYM_ORD]={wgt1,wgt2,wgt3,wgt4,wgt5,wgt6};
     1388        static IssmPDouble* l1p [MAX_TETRA_SYM_ORD]={l11 ,l12 ,l13 ,l14 ,l15 ,l16 };
     1389        static IssmPDouble* l2p [MAX_TETRA_SYM_ORD]={l21 ,l22 ,l32 ,l24 ,l25 ,l26 };
     1390        static IssmPDouble* l3p [MAX_TETRA_SYM_ORD]={l31 ,l32 ,l33 ,l34 ,l35 ,l36 };
     1391        static IssmPDouble* l4p [MAX_TETRA_SYM_ORD]={l41 ,l42 ,l43 ,l44 ,l45 ,l46 };
     1392
     1393        static int np[MAX_TETRA_SYM_ORD]={sizeof(wgt1 )/sizeof(IssmPDouble),
     1394                sizeof(wgt2 )/sizeof(IssmPDouble),
     1395                sizeof(wgt3 )/sizeof(IssmPDouble),
     1396                sizeof(wgt4 )/sizeof(IssmPDouble),
     1397                sizeof(wgt5 )/sizeof(IssmPDouble),
     1398                sizeof(wgt6 )/sizeof(IssmPDouble)};
    13991399
    14001400        //      _printf_(true,"GaussLegendreTetra: iord=%d\n",iord);
     
    14101410                *pngaus=np[iord-1];
    14111411
    1412                 *pl1  = (double *) xmalloc(*pngaus*sizeof(double));
    1413                 *pl2  = (double *) xmalloc(*pngaus*sizeof(double));
    1414                 *pl3  = (double *) xmalloc(*pngaus*sizeof(double));
    1415                 *pl4  = (double *) xmalloc(*pngaus*sizeof(double));
    1416                 *pwgt = (double *) xmalloc(*pngaus*sizeof(double));
     1412                *pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1413                *pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1414                *pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1415                *pl4  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1416                *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
    14171417
    14181418                for (i=0; i<*pngaus; i++) {
     
    14301430                *pngaus=nigaus*nigaus*nigaus;
    14311431
    1432                 *pl1  = (double *) xmalloc(*pngaus*sizeof(double));
    1433                 *pl2  = (double *) xmalloc(*pngaus*sizeof(double));
    1434                 *pl3  = (double *) xmalloc(*pngaus*sizeof(double));
    1435                 *pl4  = (double *) xmalloc(*pngaus*sizeof(double));
    1436                 *pwgt = (double *) xmalloc(*pngaus*sizeof(double));
     1432                *pl1  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1433                *pl2  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1434                *pl3  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1435                *pl4  = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
     1436                *pwgt = (IssmPDouble *) xmalloc(*pngaus*sizeof(IssmPDouble));
    14371437
    14381438                /*  get the gauss points in each direction  */
     
    14721472}/*}}}*/
    14731473/*FUNCTION GaussLobatto{{{*/
    1474 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus ) {
     1474void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus ) {
    14751475        /*Gauss-Lobatto quadrature points.
    14761476
     
    14881488
    14891489        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;
     1490        IssmPDouble *alpha,*beta;
     1491        IssmPDouble left=-1.,right= 1.;
     1492        IssmPDouble p0l=0.,p0r=0.,p1l=1.,p1r=1.,pm1l,pm1r,det;
    14931493
    14941494        /*p= 1, npoint= 1 (Gauss-Legendre)*/
    1495         static double wgt1[]={2.000000000000000};
    1496         static double xi1[]={0.000000000000000};
     1495        static IssmPDouble wgt1[]={2.000000000000000};
     1496        static IssmPDouble xi1[]={0.000000000000000};
    14971497
    14981498        /*p= 1, npoint= 2*/
    1499         static double wgt2[]={1.000000000000000, 1.000000000000000};
    1500         static double xi2[]={-1.000000000000000, 1.000000000000000};
     1499        static IssmPDouble wgt2[]={1.000000000000000, 1.000000000000000};
     1500        static IssmPDouble xi2[]={-1.000000000000000, 1.000000000000000};
    15011501
    15021502        /*p= 3, npoint= 3*/
    1503         static double wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
    1504         static double xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
     1503        static IssmPDouble wgt3[]={0.333333333333333, 1.333333333333333, 0.333333333333333};
     1504        static IssmPDouble xi3[]={-1.000000000000000, 0.000000000000000, 1.000000000000000};
    15051505
    15061506        /*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};
     1507        static IssmPDouble wgt4[]={0.166666666666667, 0.833333333333333, 0.833333333333333, 0.166666666666667};
     1508        static IssmPDouble xi4[]={-1.000000000000000,-0.447213595499958, 0.447213595499958, 1.000000000000000};
    15091509
    15101510        /*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)};
     1511        static IssmPDouble wgt5[]={0.100000000000000, 0.544444444444444, 0.711111111111111, 0.544444444444444, 0.100000000000000};
     1512        static IssmPDouble xi5[]={-1.000000000000000,-0.654653670707977, 0.000000000000000, 0.654653670707977, 1.000000000000000};
     1513
     1514        static IssmPDouble* wgtp[MAX_LINE_GLOB_PTS]={wgt1 ,wgt2 ,wgt3 ,wgt4 ,wgt5 };
     1515        static IssmPDouble* xip [MAX_LINE_GLOB_PTS]={xi1  ,xi2  ,xi3  ,xi4  ,xi5  };
     1516
     1517        static int np[MAX_LINE_GLOB_PTS]={sizeof(wgt1 )/sizeof(IssmPDouble),
     1518                sizeof(wgt2 )/sizeof(IssmPDouble),
     1519                sizeof(wgt3 )/sizeof(IssmPDouble),
     1520                sizeof(wgt4 )/sizeof(IssmPDouble),
     1521                sizeof(wgt5 )/sizeof(IssmPDouble)};
    15221522
    15231523        //      _printf_(true,"Gauss-Lobatto recurrence coefficients ngaus=%d\n",ngaus);
    1524         *pxgaus = (double *) xmalloc(ngaus*sizeof(double));
    1525         *pxwgt  = (double *) xmalloc(ngaus*sizeof(double));
     1524        *pxgaus = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
     1525        *pxwgt  = (IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
    15261526
    15271527        /*  check to see if Gauss points need to be calculated  */
     
    15391539
    15401540                /*  calculate the Gauss points using recurrence relations  */
    1541                 alpha=(double *) xmalloc(ngaus*sizeof(double));
    1542                 beta =(double *) xmalloc(ngaus*sizeof(double));
     1541                alpha=(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
     1542                beta =(IssmPDouble *) xmalloc(ngaus*sizeof(IssmPDouble));
    15431543
    15441544                /*  calculate the Legendre recurrence coefficients  */
     
    15811581}/*}}}*/
    15821582/*FUNCTION GaussRecur{{{*/
    1583 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta ) {
     1583void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta ) {
    15841584        /*Gauss quadrature points from recursion coefficients.
    15851585         *
     
    15891589        /*Intermediaries*/
    15901590        int i,j,k,l,m,ii,mml,iter;
    1591         double p,g,r,s,c,f,b;
    1592         double* work;
     1591        IssmPDouble p,g,r,s,c,f,b;
     1592        IssmPDouble* work;
    15931593
    15941594        if (n==1){
     
    15981598        }
    15991599
    1600         work=(double*)xmalloc(n*sizeof(double));
     1600        work=(IssmPDouble*)xmalloc(n*sizeof(IssmPDouble));
    16011601
    16021602        zero[0]  =alpha[0];
     
    17041704/*Element Gauss points TO BE REMOVED*/
    17051705/*FUNCTION gaussQuad{{{*/
    1706 void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus ) {
     1706void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ) {
    17071707        /*Gauss quadrature points for the quadrilaterial.*/
    17081708
     
    17121712}/*}}}*/
    17131713/*FUNCTION gaussHexa{{{*/
    1714 void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus ) {
     1714void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ) {
    17151715        /*Gauss quadrature points for the hexahedron.*/
    17161716
  • issm/trunk-jpl/src/c/shared/Numerics/GaussPoints.h

    r5808 r12424  
    77
    88#define MAX_LINE_GAUS_PTS    4
    9 void GaussLegendreLinear( double** pxgaus, double** pxwgt, int ngaus );
     9void GaussLegendreLinear( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus );
    1010#define MAX_TRIA_SYM_ORD    20
    11 void GaussLegendreTria( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, int iord );
     11void GaussLegendreTria( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pwgt, int iord );
    1212#define MAX_TETRA_SYM_ORD    6
    13 void GaussLegendreTetra( int* pngaus, double** pl1, double** pl2, double** pl3, double** pl4, double** pwgt, int iord );
     13void GaussLegendreTetra( int* pngaus, IssmPDouble** pl1, IssmPDouble** pl2, IssmPDouble** pl3, IssmPDouble** pl4, IssmPDouble** pwgt, int iord );
    1414#define MAX_LINE_GLOB_PTS    5
    15 void GaussLobatto( double** pxgaus, double** pxwgt, int ngaus );
     15void GaussLobatto( IssmPDouble** pxgaus, IssmPDouble** pxwgt, int ngaus );
    1616#define MAX_GAUS_ITER   30
    17 void GaussRecur( double* zero, double* weight, int n, double* alpha, double* beta );
     17void GaussRecur( IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta );
    1818
    19 void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
    20 void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
     19void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus );
     20void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus );
    2121
    2222#endif
Note: See TracChangeset for help on using the changeset viewer.