Changeset 5638


Ignore:
Timestamp:
08/31/10 14:26:29 (15 years ago)
Author:
Mathieu Morlighem
Message:

more GaussTria

Location:
issm/trunk/src/c/objects/Elements
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r5637 r5638  
    10071007
    10081008        /* gaussian points: */
    1009         int     num_gauss,ig;
    1010         double* first_gauss_area_coord  =  NULL;
    1011         double* second_gauss_area_coord =  NULL;
    1012         double* third_gauss_area_coord  =  NULL;
    1013         double* gauss_weights           =  NULL;
    1014         double  gauss_weight;
    1015         double  gauss_l1l2l3[3];
     1009        int     ig;
     1010        GaussTria *gauss=NULL;
    10161011
    10171012        /* parameters: */
    … …  
    10351030        double l1l2l3[3];
    10361031        double    alpha2complement_list[numvertices];                          //TO BE DELETED
    1037         double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     1032        double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    10381033
    10391034        /* strain rate: */
    … …  
    10751070        /*COMPUT alpha2_list (TO BE DELETED)*/
    10761071        for(i=0;i<numvertices;i++){
    1077                 friction->GetAlphaComplement(&alpha2complement_list[i],&gauss[i][0],VxEnum,VyEnum);
    1078         }
    1079 
    1080         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1081         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     1072                friction->GetAlphaComplement(&alpha2complement_list[i],&gaussgrids[i][0],VxEnum,VyEnum);
     1073        }
    10821074
    10831075        /*Retrieve all inputs we will be needing: */
    … …  
    10891081
    10901082        /* Start  looping on the number of gaussian points: */
    1091         for (ig=0; ig<num_gauss; ig++){
    1092                 /*Pick up the gaussian point: */
    1093                 gauss_weight=*(gauss_weights+ig);
    1094                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1095                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1096                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1083        gauss=new GaussTria(4);
     1084        for (ig=gauss->begin();ig<gauss->end();ig++){
     1085
     1086                gauss->GaussPoint(ig);
    10971087
    10981088                /*Build alpha_complement_list: */
    1099                 //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxEnum,VyEnum); // TO BE UNCOMMENTED
     1089                //if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); // TO BE UNCOMMENTED
    11001090                //else alpha_complement=0;
    1101                 TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss_l1l2l3); // TO BE DELETED
     1091                TriaRef::GetParameterValue(&alpha_complement,&alpha2complement_list[0],gauss); // TO BE DELETED
    11021092       
    11031093                /*Recover alpha_complement and k: */
    1104                 dragcoefficient_input->GetParameterValue(&drag, gauss_l1l2l3);
     1094                dragcoefficient_input->GetParameterValue(&drag, gauss);
    11051095
    11061096                /*recover lambda and mu: */
    1107                 adjointx_input->GetParameterValue(&lambda, gauss_l1l2l3);
    1108                 adjointy_input->GetParameterValue(&mu, gauss_l1l2l3);
     1097                adjointx_input->GetParameterValue(&lambda, gauss);
     1098                adjointy_input->GetParameterValue(&mu, gauss);
    11091099                       
    11101100                /*recover vx and vy: */
    1111                 inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
    1112                 inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
     1101                vx_input->GetParameterValue(&vx,gauss);
     1102                vy_input->GetParameterValue(&vy,gauss);
    11131103
    11141104                /* Get Jacobian determinant: */
    1115                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1105                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    11161106               
    11171107                /* Get nodal functions value at gaussian point:*/
    1118                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     1108                GetNodalFunctions(l1l2l3, gauss);
    11191109
    11201110                /*Get nodal functions derivatives*/
    1121                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
     1111                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
    11221112
    11231113                /*Get k derivative: dk/dx */
    1124                 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     1114                dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    11251115
    11261116                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    … …  
    11281118
    11291119                        //standard term dJ/dki
    1130                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
     1120                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i];
    11311121
    11321122                        //noise dampening d/dki(1/2*(dk/dx)^2)
    1133                         grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     1123                        grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    11341124                }
    11351125               
    … …  
    11411131        VecSetValues(gradient,numvertices,doflist1,(const double*)grade_g,ADD_VALUES);
    11421132
    1143         xfree((void**)&first_gauss_area_coord);
    1144         xfree((void**)&second_gauss_area_coord);
    1145         xfree((void**)&third_gauss_area_coord);
    1146         xfree((void**)&gauss_weights);
     1133        /*Clean up and return*/
     1134        delete gauss;
    11471135        delete friction;
    11481136
    … …  
    36723660
    36733661        /* gaussian points: */
    3674         int     num_gauss,ig;
    3675         double* first_gauss_area_coord  =  NULL;
    3676         double* second_gauss_area_coord =  NULL;
    3677         double* third_gauss_area_coord  =  NULL;
    3678         double* gauss_weights           =  NULL;
    3679         double  gauss_weight;
    3680         double  gauss_l1l2l3[3];
    3681 
     3662        int     ig;
     3663        GaussTria *gauss=NULL;
    36823664
    36833665        /* surface normal: */
    … …  
    37043686        GetDofList(&doflist);
    37053687
    3706         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3707         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    3708 
    37093688        /*Build normal vector to the surface:*/
    3710 
    37113689        x4=xyz_list[0][0];
    37123690        y4=xyz_list[0][1];
    … …  
    37373715
    37383716        /* Start  looping on the number of gaussian points: */
    3739         for (ig=0; ig<num_gauss; ig++){
    3740                 /*Pick up the gaussian point: */
    3741                 gauss_weight=*(gauss_weights+ig);
    3742                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3743                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3744                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3717        gauss=new GaussTria(2);
     3718        for (ig=gauss->begin();ig<gauss->end();ig++){
     3719
     3720                gauss->GaussPoint(ig);
    37453721
    37463722                /* Get Jacobian determinant: */
    3747                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     3723                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
    37483724
    37493725                //Get L matrix if viscous basal drag present:
    3750                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     3726                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    37513727
    37523728                /**********************Do not forget the sign**********************************/
    3753                 DL_scalar=- gauss_weight*Jdet*nz;
     3729                DL_scalar=- gauss->weight*Jdet*nz;
    37543730                /******************************************************************************/
    37553731
    … …  
    37693745        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    37703746
    3771         xfree((void**)&first_gauss_area_coord);
    3772         xfree((void**)&second_gauss_area_coord);
    3773         xfree((void**)&third_gauss_area_coord);
    3774         xfree((void**)&gauss_weights);
     3747        /*Clean up and return*/
     3748        delete gauss;
    37753749        xfree((void**)&doflist);
    37763750}
    … …  
    37943768
    37953769        /* gaussian points: */
    3796         int     num_area_gauss,ig;
    3797         double* gauss_weights  =  NULL;
    3798         double* first_gauss_area_coord  =  NULL;
    3799         double* second_gauss_area_coord =  NULL;
    3800         double* third_gauss_area_coord  =  NULL;
    3801         double  gauss_weight;
    3802         double  gauss_coord[3];
     3770        int     ig;
     3771        GaussTria *gauss=NULL;
    38033772
    38043773        /*matrices: */
    … …  
    38183787        GetDofList(&doflist);
    38193788
    3820         /* Get gaussian points and weights: */
    3821         GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    38223789
    38233790        /* Start looping on the number of gauss  (nodes on the bedrock) */
    3824         for (ig=0; ig<num_area_gauss; ig++){
    3825                 gauss_weight=*(gauss_weights+ig);
    3826                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    3827                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    3828                 gauss_coord[2]=*(third_gauss_area_coord+ig);
    3829 
    3830                 //Get the Jacobian determinant
    3831                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
     3791        gauss=new GaussTria(2);
     3792        for (ig=gauss->begin();ig<gauss->end();ig++){
     3793
     3794                gauss->GaussPoint(ig);
     3795
     3796                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
    38323797
    38333798                /*Get L matrix : */
    3834                 GetL(&L[0], &xyz_list[0][0], gauss_coord,NDOF1);
     3799                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    38353800
    38363801                /*Calculate DL on gauss point */
    3837                 D_scalar=latentheat/heatcapacity*gauss_weight*Jdet;
     3802                D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
    38383803
    38393804                /*  Do the triple product tL*D*L: */
    … …  
    38513816        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    38523817
    3853         xfree((void**)&first_gauss_area_coord);
    3854         xfree((void**)&second_gauss_area_coord);
    3855         xfree((void**)&third_gauss_area_coord);
    3856         xfree((void**)&gauss_weights);
     3818        /*Clean up and return*/
     3819        delete gauss;
    38573820        xfree((void**)&doflist);
    38583821}
    … …  
    38733836
    38743837        /* gaussian points: */
    3875         int     num_gauss,ig;
    3876         double* first_gauss_area_coord  =  NULL;
    3877         double* second_gauss_area_coord =  NULL;
    3878         double* third_gauss_area_coord  =  NULL;
    3879         double* gauss_weights           =  NULL;
    3880         double  gauss_weight;
    3881         double  gauss_l1l2l3[3];
     3838        int     ig;
     3839        GaussTria *gauss=NULL;
    38823840
    38833841        /* matrices: */
    … …  
    39333891        if(artdiff){
    39343892                //Get the Jacobian determinant
    3935                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    3936                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     3893                gauss=new GaussTria();
     3894                gauss->GaussCenter();
     3895                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     3896                delete gauss;
    39373897
    39383898                //Build K matrix (artificial diffusivity matrix)
    … …  
    39443904        }
    39453905
    3946         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3947         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    3948 
    39493906        /* Start  looping on the number of gaussian points: */
    3950         for (ig=0; ig<num_gauss; ig++){
    3951 
    3952                 /*Pick up the gaussian point: */
    3953                 gauss_weight=*(gauss_weights+ig);
    3954                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3955                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3956                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3907        gauss=new GaussTria(2);
     3908        for (ig=gauss->begin();ig<gauss->end();ig++){
     3909
     3910                gauss->GaussPoint(ig);
    39573911
    39583912                /* Get Jacobian determinant: */
    3959                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     3913                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    39603914
    39613915                /*Get L matrix: */
    3962                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
    3963 
    3964                 DL_scalar=gauss_weight*Jdettria;
     3916                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
     3917
     3918                DL_scalar=gauss->weight*Jdettria;
    39653919
    39663920                /*  Do the triple product tL*D*L: */
    … …  
    39713925
    39723926                /*Get B  and B prime matrix: */
    3973                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    3974                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     3927                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     3928                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    39753929
    39763930                //Get vx, vy and their derivatives at gauss point
    3977                 vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]);
    3978                 vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]);
    3979 
    3980                 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
    3981                 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     3931                vxaverage_input->GetParameterValue(&vx,gauss);
     3932                vyaverage_input->GetParameterValue(&vy,gauss);
     3933
     3934                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     3935                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    39823936
    39833937                dvxdx=dvx[0];
    39843938                dvydy=dvy[1];
    39853939
    3986                 DL_scalar=dt*gauss_weight*Jdettria;
     3940                DL_scalar=dt*gauss->weight*Jdettria;
    39873941
    39883942                //Create DL and DLprime matrix
    … …  
    40323986        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    40333987
    4034         xfree((void**)&first_gauss_area_coord);
    4035         xfree((void**)&second_gauss_area_coord);
    4036         xfree((void**)&third_gauss_area_coord);
    4037         xfree((void**)&gauss_weights);
     3988        /*Clean up and return*/
     3989        delete gauss;
    40383990        xfree((void**)&doflist);
    40393991}
    … …  
    40544006
    40554007        /* gaussian points: */
    4056         int     num_gauss,ig;
    4057         double* first_gauss_area_coord  =  NULL;
    4058         double* second_gauss_area_coord =  NULL;
    4059         double* third_gauss_area_coord  =  NULL;
    4060         double* gauss_weights           =  NULL;
    4061         double  gauss_weight;
    4062         double  gauss_l1l2l3[3];
     4008        int     ig;
     4009        GaussTria *gauss=NULL;
    40634010
    40644011        /* matrices: */
    … …  
    40924039        GetDofList(&doflist);
    40934040
    4094         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4095         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4096 
    40974041        /*Retrieve all inputs we will be needed*/
    40984042        if(dim==2){
    … …  
    41064050
    41074051        /* Start  looping on the number of gaussian points: */
    4108         for (ig=0; ig<num_gauss; ig++){
    4109 
    4110                 /*Pick up the gaussian point: */
    4111                 gauss_weight=*(gauss_weights+ig);
    4112                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4113                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4114                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4052        gauss=new GaussTria(2);
     4053        for (ig=gauss->begin();ig<gauss->end();ig++){
     4054
     4055                gauss->GaussPoint(ig);
    41154056
    41164057                /* Get Jacobian determinant: */
    4117                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     4058                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    41184059
    41194060                /*Get L matrix: */
    4120                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
    4121 
    4122                 DL_scalar=gauss_weight*Jdettria;
     4061                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
     4062
     4063                DL_scalar=gauss->weight*Jdettria;
    41234064
    41244065                /*  Do the triple product tL*D*L: */
    … …  
    41304071                /*Get B  and B prime matrix: */
    41314072                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    4132                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
    4133                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     4073                GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     4074                GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
    41344075
    41354076                //Get vx, vy and their derivatives at gauss point
    4136                 vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]);
    4137                 vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]);
    4138 
    4139                 DL_scalar=-dt*gauss_weight*Jdettria;
     4077                vxaverage_input->GetParameterValue(&vx,gauss);
     4078                vyaverage_input->GetParameterValue(&vy,gauss);
     4079
     4080                DL_scalar=-dt*gauss->weight*Jdettria;
    41404081
    41414082                DLprime[0][0]=DL_scalar*vx;
    … …  
    41574098        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    41584099
    4159         xfree((void**)&first_gauss_area_coord);
    4160         xfree((void**)&second_gauss_area_coord);
    4161         xfree((void**)&third_gauss_area_coord);
    4162         xfree((void**)&gauss_weights);
     4100        /*Clean up and return*/
     4101        delete gauss;
    41634102        xfree((void**)&doflist);
    41644103}
    … …  
    42314170        double heatcapacity;
    42324171
    4233         int     num_gauss,ig;
    4234         double* first_gauss_area_coord  =  NULL;
    4235         double* second_gauss_area_coord =  NULL;
    4236         double* third_gauss_area_coord  =  NULL;
    4237         double* gauss_weights           =  NULL;
    4238         double  gauss_weight;
    4239         double  gauss_coord[3];
     4172        int     ig;
     4173        GaussTria *gauss=NULL;
    42404174
    42414175        /*matrices: */
    … …  
    42644198        heatcapacity=matpar->GetHeatCapacity();
    42654199
    4266         GaussLegendreTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4267 
    42684200        /* Start looping on the number of gauss (nodes on the bedrock) */
    4269         for (ig=0; ig<num_gauss; ig++){
    4270                 gauss_weight=*(gauss_weights+ig);
    4271                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    4272                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    4273                 gauss_coord[2]=*(third_gauss_area_coord+ig);
     4201        gauss=new GaussTria(2);
     4202        for (ig=gauss->begin();ig<gauss->end();ig++){
     4203
     4204                gauss->GaussPoint(ig);
    42744205               
    42754206                //Get the Jacobian determinant
    4276                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss_coord);
     4207                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
    42774208               
    42784209                /*Get nodal functions values: */
    4279                 GetNodalFunctions(&l1l2l3[0], gauss_coord);
     4210                GetNodalFunctions(&l1l2l3[0], gauss);
    42804211                               
    42814212                /*Calculate DL on gauss point */
    4282                 D_scalar=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
     4213                D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
    42834214                if(dt){
    42844215                        D_scalar=dt*D_scalar;
    … …  
    42994230        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    43004231
    4301         xfree((void**)&first_gauss_area_coord);
    4302         xfree((void**)&second_gauss_area_coord);
    4303         xfree((void**)&third_gauss_area_coord);
    4304         xfree((void**)&gauss_weights);
     4232
     4233        /*Clean up and return*/
     4234        delete gauss;
    43054235        xfree((void**)&doflist);
    43064236}
    … …  
    43214251
    43224252        /* gaussian points: */
    4323         int     num_gauss,ig;
    4324         double* first_gauss_area_coord  =  NULL;
    4325         double* second_gauss_area_coord =  NULL;
    4326         double* third_gauss_area_coord  =  NULL;
    4327         double* gauss_weights           =  NULL;
    4328         double  gauss_weight;
    4329         double  gauss_l1l2l3[3];
     4253        int     ig;
     4254        GaussTria* gauss=NULL;
    43304255
    43314256        /* matrix */
    … …  
    43474272        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    43484273        GetDofList(&doflist);
    4349 
    4350         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4351         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    43524274
    43534275        /*retrieve inputs :*/
    … …  
    43574279       
    43584280        /* Start  looping on the number of gaussian points: */
    4359         for (ig=0; ig<num_gauss; ig++){
    4360                 /*Pick up the gaussian point: */
    4361                 gauss_weight=*(gauss_weights+ig);
    4362                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4363                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4364                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4281        gauss=new GaussTria(2);
     4282        for(ig=gauss->begin();ig<gauss->end();ig++){
     4283
     4284                gauss->GaussPoint(ig);
    43654285
    43664286                /* Get Jacobian determinant: */
    4367                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     4287                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    43684288
    43694289                /*Get L matrix: */
    4370                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     4290                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
    43714291
    43724292                /* Get accumulation, melting at gauss point */
    4373                 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);
    4374                 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);
    4375                 dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]);
     4293                accumulation_input->GetParameterValue(&accumulation_g,gauss);
     4294                melting_input->GetParameterValue(&melting_g,gauss);
     4295                dhdt_input->GetParameterValue(&dhdt_g,gauss);
    43764296
    43774297                /* Add value into pe_g: */
    4378                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g-dhdt_g)*L[i];
     4298                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
    43794299
    43804300        } // for (ig=0; ig<num_gauss; ig++)
    … …  
    43834303        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    43844304
    4385         xfree((void**)&first_gauss_area_coord);
    4386         xfree((void**)&second_gauss_area_coord);
    4387         xfree((void**)&third_gauss_area_coord);
    4388         xfree((void**)&gauss_weights);
     4305        /*Clean up and return*/
     4306        delete gauss;
    43894307        xfree((void**)&doflist);
    43904308}
    … …  
    44054323
    44064324        /* gaussian points: */
    4407         int     num_gauss,ig;
    4408         double* first_gauss_area_coord  =  NULL;
    4409         double* second_gauss_area_coord =  NULL;
    4410         double* third_gauss_area_coord  =  NULL;
    4411         double* gauss_weights           =  NULL;
    4412         double  gauss_weight;
    4413         double  gauss_l1l2l3[3];
     4325        int     ig;
     4326        GaussTria* gauss=NULL;
    44144327
    44154328        /* matrix */
    … …  
    44324345        GetDofList(&doflist);
    44334346
    4434         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4435         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4436 
    44374347        /*Retrieve all inputs we will be needing: */
    44384348        accumulation_input=inputs->GetInput(AccumulationRateEnum);
    … …  
    44414351
    44424352        /* Start  looping on the number of gaussian points: */
    4443         for (ig=0; ig<num_gauss; ig++){
    4444                 /*Pick up the gaussian point: */
    4445                 gauss_weight=*(gauss_weights+ig);
    4446                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4447                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4448                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4353        gauss=new GaussTria(2);
     4354        for(ig=gauss->begin();ig<gauss->end();ig++){
     4355
     4356                gauss->GaussPoint(ig);
    44494357
    44504358                /* Get Jacobian determinant: */
    4451                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     4359                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    44524360
    44534361                /*Get L matrix: */
    4454                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     4362                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
    44554363
    44564364                /* Get accumulation, melting and thickness at gauss point */
    4457                 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);
    4458                 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);
    4459                 dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]);
     4365                accumulation_input->GetParameterValue(&accumulation_g,gauss);
     4366                melting_input->GetParameterValue(&melting_g,gauss);
     4367                dhdt_input->GetParameterValue(&dhdt_g,gauss);
    44604368
    44614369                /* Add value into pe_g: */
    4462                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g-dhdt_g)*L[i];
     4370                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
    44634371
    44644372        } // for (ig=0; ig<num_gauss; ig++)
    … …  
    44674375        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    44684376
    4469         xfree((void**)&first_gauss_area_coord);
    4470         xfree((void**)&second_gauss_area_coord);
    4471         xfree((void**)&third_gauss_area_coord);
    4472         xfree((void**)&gauss_weights);
     4377        /*Clean up and return*/
     4378        delete gauss;
    44734379        xfree((void**)&doflist);
    44744380}
    … …  
    44894395
    44904396        /* gaussian points: */
    4491         int     num_gauss,ig;
    4492         double* first_gauss_area_coord  =  NULL;
    4493         double* second_gauss_area_coord =  NULL;
    4494         double* third_gauss_area_coord  =  NULL;
    4495         double* gauss_weights           =  NULL;
    4496         double  gauss_weight;
    4497         double  gauss_l1l2l3[3];
     4397        int     ig;
     4398        GaussTria* gauss=NULL;
    44984399
    44994400        /* matrix */
    … …  
    45144415        GetDofList(&doflist);
    45154416
    4516         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4517         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4518 
    45194417        /*Retrieve all inputs we will be needing: */
    45204418        accumulation_input=inputs->GetInput(AccumulationRateEnum);
    … …  
    45224420
    45234421        /* Start  looping on the number of gaussian points: */
    4524         for (ig=0; ig<num_gauss; ig++){
    4525                 /*Pick up the gaussian point: */
    4526                 gauss_weight=*(gauss_weights+ig);
    4527                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4528                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4529                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4422        gauss=new GaussTria(2);
     4423        for(ig=gauss->begin();ig<gauss->end();ig++){
     4424
     4425                gauss->GaussPoint(ig);
    45304426
    45314427                /* Get Jacobian determinant: */
    4532                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     4428                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    45334429
    45344430                /*Get L matrix: */
    4535                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     4431                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
    45364432
    45374433                /* Get accumulation, melting at gauss point */
    4538                 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);
    4539                 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);
     4434                accumulation_input->GetParameterValue(&accumulation_g,gauss);
     4435                melting_input->GetParameterValue(&melting_g,gauss);
    45404436
    45414437                /* Add value into pe_g: */
    4542                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
     4438                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i];
    45434439
    45444440        } // for (ig=0; ig<num_gauss; ig++)
    … …  
    45474443        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    45484444
    4549         /*Free ressources:*/
    4550         xfree((void**)&first_gauss_area_coord);
    4551         xfree((void**)&second_gauss_area_coord);
    4552         xfree((void**)&third_gauss_area_coord);
    4553         xfree((void**)&gauss_weights);
     4445        /*Clean up and return*/
     4446        delete gauss;
    45544447        xfree((void**)&doflist);
    45554448}
    … …  
    45684461
    45694462        /* gaussian points: */
    4570         int     num_gauss,ig;
    4571         double* first_gauss_area_coord  =  NULL;
    4572         double* second_gauss_area_coord =  NULL;
    4573         double* third_gauss_area_coord  =  NULL;
    4574         double* gauss_weights           =  NULL;
    4575         double  gauss_weight;
    4576         double  gauss_l1l2l3[3];
     4463        int     ig;
     4464        GaussTria* gauss=NULL;
    45774465
    45784466        /* Jacobian: */
    … …  
    46054493        GetDofList(&doflist);
    46064494
    4607         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4608         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4609 
    46104495        /*Retrieve all inputs we will be needing: */
    46114496        melting_input=inputs->GetInput(MeltingRateEnum);
    … …  
    46164501        /*For icesheets: */
    46174502        /* Start  looping on the number of gaussian points: */
    4618         for (ig=0; ig<num_gauss; ig++){
    4619 
    4620                 /*Pick up the gaussian point: */
    4621                 gauss_weight=*(gauss_weights+ig);
    4622                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4623                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4624                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4503        gauss=new GaussTria(2);
     4504        for(ig=gauss->begin();ig<gauss->end();ig++){
     4505
     4506                gauss->GaussPoint(ig);
    46254507
    46264508                /*Get melting at gaussian point: */
    4627                 melting_input->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0]);
     4509                melting_input->GetParameterValue(&meltingvalue, gauss);
    46284510
    46294511                /*Get velocity at gaussian point: */
    4630                 vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);
    4631                 vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);
     4512                vx_input->GetParameterValue(&vx, gauss);
     4513                vy_input->GetParameterValue(&vy, gauss);
    46324514
    46334515                /*Get bed slope: */
    4634                 bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     4516                bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    46354517                dbdx=slope[0];
    46364518                dbdy=slope[1];
    46374519
    46384520                /* Get Jacobian determinant: */
    4639                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4521                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
    46404522
    46414523                //Get L matrix if viscous basal drag present:
    4642                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     4524                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    46434525
    46444526                /*Build gaussian vector: */
    46454527                for(i=0;i<numvertices;i++){
    4646                         pe_g_gaussian[i]=-Jdet*gauss_weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
     4528                        pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
    46474529                }
    46484530
    … …  
    46554537        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    46564538
    4657         /*Free ressources:*/
    4658         xfree((void**)&first_gauss_area_coord);
    4659         xfree((void**)&second_gauss_area_coord);
    4660         xfree((void**)&third_gauss_area_coord);
    4661         xfree((void**)&gauss_weights);
     4539        /*Clean up and return*/
     4540        delete gauss;
    46624541        xfree((void**)&doflist);
    46634542}

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

    r5636 r5638  
    487487}
    488488/*}}}*/
     489/*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/
     490void TriaRef::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,GaussTria* gauss){
     491        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     492         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     493
     494        double x1,x2,x3,y1,y2,y3,z1,z2,z3;
     495
     496        x1=*(xyz_list+3*0+0);
     497        y1=*(xyz_list+3*0+1);
     498        z1=*(xyz_list+3*0+2);
     499        x2=*(xyz_list+3*1+0);
     500        y2=*(xyz_list+3*1+1);
     501        z2=*(xyz_list+3*1+2);
     502        x3=*(xyz_list+3*2+0);
     503        y3=*(xyz_list+3*2+1);
     504        z3=*(xyz_list+3*2+2);
     505
     506        *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
     507        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     508
     509}
     510/*}}}*/
    489511/*FUNCTION TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss) {{{1*/
    490512void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){

  • issm/trunk/src/c/objects/Elements/TriaRef.h

    r5636 r5638  
    4040                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);
    4141                void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
     42                void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss);
    4243                void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
    4344                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss);
Note: See TracChangeset for help on using the changeset viewer.