Changeset 5677


Ignore:
Timestamp:
09/03/10 14:02:48 (15 years ago)
Author:
seroussi
Message:

removed ugly gauss in KMatrixStokes

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

Legend:

Unmodified
Added
Removed

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

    r5674 r5677  
    27462746
    27472747        /* gaussian points: */
    2748         int     num_area_gauss;
    2749         int     igarea,igvert;
    2750         double* first_gauss_area_coord  =  NULL;
    2751         double* second_gauss_area_coord =  NULL;
    2752         double* third_gauss_area_coord  =  NULL;
    2753         double* vert_gauss_coord = NULL;
    2754         double* area_gauss_weights  =  NULL;
    2755         double* vert_gauss_weights  =  NULL;
    2756         double  gaussgrids[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
     2748        int     ig;
     2749        GaussPenta *gauss=NULL;
     2750        GaussTria  *gauss_tria=NULL;
     2751
    27572752
    27582753        /* specific gaussian point: */
    … …  
    28062801                points, same deal, which yields 3 gaussian points.*/
    28072802
    2808         area_order=5;
    2809         num_vert_gauss=5;
    2810 
    2811         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2812         gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    2813 
    28142803        /*Retrieve all inputs we will be needing: */
    28152804        vx_input=inputs->GetInput(VxEnum);
    … …  
    28182807
    28192808        /* Start  looping on the number of gaussian points: */
    2820         for (igarea=0; igarea<num_area_gauss; igarea++){
    2821                 for (igvert=0; igvert<num_vert_gauss; igvert++){
    2822                         /*Pick up the gaussian point: */
    2823                         area_gauss_weight=*(area_gauss_weights+igarea);
    2824                         vert_gauss_weight=*(vert_gauss_weights+igvert);
    2825                         gauss_weight=area_gauss_weight*vert_gauss_weight;
    2826                         gauss_coord[0]=*(first_gauss_area_coord+igarea);
    2827                         gauss_coord[1]=*(second_gauss_area_coord+igarea);
    2828                         gauss_coord[2]=*(third_gauss_area_coord+igarea);
    2829                         gauss_coord[3]=*(vert_gauss_coord+igvert);
    2830 
    2831                         /*Compute strain rate: */
    2832                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    2833 
    2834                         /*Get viscosity: */
    2835                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    2836 
    2837                         /*Get B and Bprime matrices: */
    2838                         GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord);
    2839                         GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord);
    2840 
    2841                         /* Get Jacobian determinant: */
    2842                         GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]);
    2843 
    2844                         /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    2845                          * onto this scalar matrix, so that we win some computational time: */
    2846                         D_scalar=gauss_weight*Jdet;
    2847                         for (i=0;i<6;i++){
    2848                                 D[i][i]=D_scalar*2*viscosity;
    2849                         }
    2850                         for (i=6;i<8;i++){
    2851                                 D[i][i]=-D_scalar*stokesreconditioning;
    2852                         }
    2853 
    2854                         /*  Do the triple product tB*D*Bprime: */
    2855                         TripleMultiply( &B[0][0],8,27,1,
    2856                                                 &D[0][0],8,8,0,
    2857                                                 &B_prime[0][0],8,27,0,
    2858                                                 &Ke_gaussian[0][0],0);
    2859 
    2860                         /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
    2861                         for(i=0;i<27;i++){
    2862                                 for(j=0;j<27;j++){
    2863                                         Ke_temp[i][j]+=Ke_gaussian[i][j];
    2864                                 }
    2865                         }
    2866                 }
     2809        /* Start  looping on the number of gaussian points: */
     2810        gauss=new GaussPenta(5,5);
     2811        for (ig=gauss->begin();ig<gauss->end();ig++){
     2812
     2813                gauss->GaussPoint(ig);
     2814
     2815                /*Compute strain rate: */
     2816                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     2817
     2818                /*Get viscosity: */
     2819                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     2820
     2821                /*Get B and Bprime matrices: */
     2822                GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
     2823                GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss);
     2824
     2825                /* Get Jacobian determinant: */
     2826                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     2827
     2828                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     2829                 * onto this scalar matrix, so that we win some computational time: */
     2830                D_scalar=gauss->weight*Jdet;
     2831                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
     2832                for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     2833
     2834                /*  Do the triple product tB*D*Bprime: */
     2835                TripleMultiply( &B[0][0],8,27,1,
     2836                                        &D[0][0],8,8,0,
     2837                                        &B_prime[0][0],8,27,0,
     2838                                        &Ke_gaussian[0][0],0);
     2839
     2840                /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
     2841                for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    28672842        }
    28682843
    … …  
    28782853                }
    28792854
    2880                 xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
    2881                 GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    2882 
    2883                 /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    2884                 for (igarea=0; igarea<num_area_gauss; igarea++){
    2885                         gauss_weight=*(area_gauss_weights+igarea);
    2886                         gauss_coord[0]=*(first_gauss_area_coord+igarea);
    2887                         gauss_coord[1]=*(second_gauss_area_coord+igarea);
    2888                         gauss_coord[2]=*(third_gauss_area_coord+igarea);
    2889                         gauss_coord[3]=-1;
    2890 
    2891                         gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
    2892                         gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    2893                         gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
     2855                /* Start  looping on the number of gaussian points: */
     2856                delete gauss; //gauss of previous loop
     2857                gauss=new GaussPenta();
     2858                gauss->GaussFaceTria(0,1,2,2);
     2859                gauss_tria=new GaussTria();
     2860                for (ig=gauss->begin();ig<gauss->end();ig++){
     2861
     2862                        gauss->GaussPoint(ig);
     2863                        gauss->SynchronizeGaussTria(gauss_tria);
    28942864
    28952865                        /*Get the Jacobian determinant */
    2896                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
     2866                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);
    28972867
    28982868                        /*Get L matrix if viscous basal drag present: */
    2899                         GetLStokes(&LStokes[0][0],  gauss_coord);
    2900                         GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord);
     2869                        GetLStokes(&LStokes[0][0], gauss);
     2870                        GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
    29012871
    29022872                        /*Compute strain rate: */
    2903                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
     2873                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    29042874
    29052875                        /*Get viscosity at last iteration: */
    … …  
    29102880
    29112881                        /*Calculate DL on gauss point */
    2912                         friction->GetAlpha2(&alpha2_gauss, gauss_coord,VxEnum,VyEnum,VzEnum);
    2913 
    2914                         DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d;
    2915                         DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d;
    2916                         DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
    2917                         DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
    2918                         DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
    2919                         DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
    2920                         DLStokes[6][6]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0];
    2921                         DLStokes[7][7]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1];
    2922                         DLStokes[8][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[2];
    2923                         DLStokes[9][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0;
    2924                         DLStokes[10][10]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0;
    2925                         DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0];
    2926                         DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1];
    2927                         DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2];
     2882                        friction->GetAlpha2(&alpha2_gauss, gauss_tria,VxEnum,VyEnum,VzEnum);
     2883
     2884                        DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     2885                        DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     2886                        DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     2887                        DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     2888                        DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     2889                        DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     2890                        DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
     2891                        DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
     2892                        DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2];
     2893                        DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0;
     2894                        DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0;
     2895                        DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
     2896                        DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
     2897                        DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2];
    29282898
    29292899                        /*  Do the triple product tL*D*L: */
    … …  
    29332903                                                &Ke_drag_gaussian[0][0],0);
    29342904
    2935                         for(i=0;i<numdof2d;i++){
    2936                                 for(j=0;j<numdof2d;j++){
    2937                                         Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
    2938                                 }
    2939                         }
     2905                        for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
    29402906                }
    29412907       
    … …  
    29522918
    29532919        /*Free ressources:*/
    2954         xfree((void**)&first_gauss_area_coord);
    2955         xfree((void**)&second_gauss_area_coord);
    2956         xfree((void**)&third_gauss_area_coord);
    2957         xfree((void**)&area_gauss_weights);
    2958         xfree((void**)&vert_gauss_coord);
    2959         xfree((void**)&vert_gauss_weights);
     2920        delete gauss;
     2921        delete gauss_tria;
    29602922        xfree((void**)&doflist);
    29612923
    … …  
    39553917
    39563918                        /*Get L matrix : */
    3957                         GetNodalFunctionsP1(l1l6, gauss);
     3919                        GetNodalFunctionsP1(l1l6, gauss_coord);
    39583920
    39593921                        /*Get water_pressure at gaussian point */
    … …  
    39653927                        for(i=0;i<NUMVERTICES2D;i++){
    39663928                                for(j=0;j<3;j++){
    3967                                         Pe_temp[i*NDOF4+j]+=water_pressure*gauss_weight*Jdet2d*L[i]*bed_normal[j];
     3929                                        Pe_temp[i*NDOF4+j]+=water_pressure*gauss_weight*Jdet2d*l1l6[i]*bed_normal[j];
    39683930                                }
    39693931                        }

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

    r5670 r5677  
    14061406}
    14071407/*}}}*/
     1408/*FUNCTION PentaRef::GetLStokes {{{1*/
     1409void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){
     1410        /*
     1411         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     1412         * For grid i, Li can be expressed in the actual coordinate system
     1413         * by:
     1414         *       Li=[ h    0    0   0]
     1415         *                    [ 0    h    0   0]
     1416         *                    [ 0    0    h   0]
     1417         *                    [ 0    0    h   0]
     1418         *                    [ h    0    0   0]
     1419         *                    [ 0    h    0   0]
     1420         *                    [ h    0    0   0]
     1421         *                    [ 0    h    0   0]
     1422         *                    [ 0    0    h   0]
     1423         *                    [ 0    0    h   0]
     1424         *                    [ 0    0    h   0]
     1425         *                    [ h    0    0   0]
     1426         *                    [ 0    h    0   0]
     1427         *                    [ 0    0    h   0]
     1428         * where h is the interpolation function for grid i.
     1429         */
     1430
     1431        int i;
     1432        const int numgrids2d=3;
     1433        int num_dof=4;
     1434
     1435        double l1l2l3[numgrids2d];
     1436
     1437
     1438        /*Get l1l2l3 in actual coordinate system: */
     1439        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     1440        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     1441        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     1442
     1443        /*Build LStokes: */
     1444        for (i=0;i<3;i++){
     1445                *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
     1446                *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
     1447                *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
     1448                *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
     1449                *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
     1450                *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
     1451                *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
     1452                *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
     1453                *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
     1454                *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
     1455                *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
     1456                *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
     1457                *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
     1458                *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
     1459                *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
     1460                *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
     1461                *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
     1462                *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
     1463                *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
     1464                *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
     1465                *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
     1466                *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
     1467                *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
     1468                *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
     1469                *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
     1470                *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
     1471                *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
     1472                *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
     1473                *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
     1474                *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
     1475                *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
     1476                *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
     1477                *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
     1478                *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
     1479                *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
     1480                *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
     1481                *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
     1482                *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
     1483                *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
     1484                *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
     1485                *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
     1486                *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
     1487                *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
     1488                *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
     1489                *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
     1490                *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
     1491                *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
     1492                *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
     1493                *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
     1494                *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
     1495                *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
     1496                *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
     1497                *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
     1498                *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
     1499                *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
     1500                *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
     1501
     1502        }
     1503}
     1504/*}}}*/
    14081505/*FUNCTION PentaRef::GetLprimeStokes {{{1*/
    1409 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss){
     1506void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
    14101507
    14111508        /*
    … …  
    14371534
    14381535        /*Get l1l2l3 in actual coordinate system: */
    1439         l1l2l3[0]=gauss_tria[0];
    1440         l1l2l3[1]=gauss_tria[1];
    1441         l1l2l3[2]=gauss_tria[2];
     1536        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     1537        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     1538        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    14421539
    14431540        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);

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

    r5670 r5677  
    7070                void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
    7171                void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss);
    72                 //void GetLStokes(double* LStokes, double* gauss_tria);
    73                 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss);
    74                 void GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss);
     72                void GetLStokes(double* LStokes, GaussPenta* gauss);
     73                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
     74                void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss);
    7575                //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    7676                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss);
Note: See TracChangeset for help on using the changeset viewer.