Changeset 5637


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

more GaussTria

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

Legend:

Unmodified
Added
Removed

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

    r5636 r5637  
    32473247
    32483248        /* gaussian points: */
    3249         int     num_gauss,ig;
    3250         double *first_gauss_area_coord  = NULL;
    3251         double *second_gauss_area_coord = NULL;
    3252         double *third_gauss_area_coord  = NULL;
    3253         double *gauss_weights           = NULL;
    3254         double  gauss_weight;
    3255         double  gauss_l1l2l3[3];
     3249        int     ig;
     3250        GaussTria *gauss=NULL;
    32563251
    32573252        /* matrices: */
    … …  
    32743269        double    alpha2;
    32753270        double    alpha2_list[numvertices];                                       //TO BE DELETED
    3276         double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     3271        double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    32773272
    32783273        double MAXSLOPE=.06; // 6 %
    … …  
    33143309        /*COMPUT alpha2_list (TO BE DELETED)*/
    33153310        for(i=0;i<numvertices;i++){
    3316                 friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
    3317         }
    3318 
    3319         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3320         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3311                friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
     3312        }
    33213313
    33223314        /* Start  looping on the number of gaussian points: */
    3323         for (ig=0; ig<num_gauss; ig++){
    3324                 /*Pick up the gaussian point: */
    3325                 gauss_weight=*(gauss_weights+ig);
    3326                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3327                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3328                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3315        gauss=new GaussTria(2);
     3316        for (ig=gauss->begin();ig<gauss->end();ig++){
     3317
     3318                gauss->GaussPoint(ig);
    33293319
    33303320                /*Friction: */
    3331                 // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3332                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
     3321                // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
     3322                TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
    33333323
    33343324                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    33353325                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    3336                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     3326                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    33373327                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    33383328
    … …  
    33423332
    33433333                /* Get Jacobian determinant: */
    3344                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     3334                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    33453335
    33463336                /*Get L matrix: */
    3347                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     3337                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    33483338
    33493339               
    3350                 DL_scalar=alpha2*gauss_weight*Jdet;
     3340                DL_scalar=alpha2*gauss->weight*Jdet;
    33513341                for (i=0;i<2;i++){
    33523342                        DL[i][i]=DL_scalar;
    … …  
    33673357        MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES);
    33683358
    3369         xfree((void**)&first_gauss_area_coord);
    3370         xfree((void**)&second_gauss_area_coord);
    3371         xfree((void**)&third_gauss_area_coord);
    3372         xfree((void**)&gauss_weights);
     3359        /*Clean up and return*/
     3360        delete gauss;
     3361        delete friction;
    33733362        xfree((void**)&doflistm);
    33743363        xfree((void**)&doflistp);
    3375         delete friction;
    33763364}       
    33773365/*}}}*/
    … …  
    33913379
    33923380        /* gaussian points: */
    3393         int     num_gauss,ig;
    3394         double *first_gauss_area_coord  = NULL;
    3395         double *second_gauss_area_coord = NULL;
    3396         double *third_gauss_area_coord  = NULL;
    3397         double *gauss_weights           = NULL;
    3398         double  gauss_weight;
    3399         double  gauss_l1l2l3[3];
     3381        int     ig;
     3382        GaussTria *gauss=NULL;
    34003383
    34013384        /* matrices: */
    … …  
    34183401        double    alpha2;
    34193402        double    alpha2_list[numvertices];                                       //TO BE DELETED
    3420         double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     3403        double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    34213404
    34223405        double MAXSLOPE=.06; // 6 %
    … …  
    34573440        /*COMPUT alpha2_list (TO BE DELETED)*/
    34583441        for(i=0;i<numvertices;i++){
    3459                 friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
    3460         }
    3461 
    3462         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3463         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3442                friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
     3443        }
    34643444
    34653445        /* Start  looping on the number of gaussian points: */
    3466         for (ig=0; ig<num_gauss; ig++){
    3467                 /*Pick up the gaussian point: */
    3468                 gauss_weight=*(gauss_weights+ig);
    3469                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3470                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3471                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3446        gauss=new GaussTria(2);
     3447        for (ig=gauss->begin();ig<gauss->end();ig++){
     3448
     3449                gauss->GaussPoint(ig);
    34723450
    34733451                /*Friction: */
    3474                 // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3475                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
     3452                // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
     3453                TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
    34763454
    34773455                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    34783456                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    3479                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     3457                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    34803458                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    34813459
    … …  
    34853463
    34863464                /* Get Jacobian determinant: */
    3487                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     3465                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    34883466
    34893467                /*Get L matrix: */
    3490                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     3468                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    34913469
    34923470               
    3493                 DL_scalar=alpha2*gauss_weight*Jdet;
     3471                DL_scalar=alpha2*gauss->weight*Jdet;
    34943472                for (i=0;i<2;i++){
    34953473                        DL[i][i]=DL_scalar;
    … …  
    35093487        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    35103488
    3511         xfree((void**)&first_gauss_area_coord);
    3512         xfree((void**)&second_gauss_area_coord);
    3513         xfree((void**)&third_gauss_area_coord);
    3514         xfree((void**)&gauss_weights);
     3489        /*Clean up and return*/
     3490        delete gauss;
     3491        delete friction;
    35153492        xfree((void**)&doflist);
    3516         delete friction;
    35173493}       
    35183494/*}}}*/
    … …  
    35323508
    35333509        /* gaussian points: */
    3534         int     num_gauss,ig;
    3535         double *first_gauss_area_coord  = NULL;
    3536         double *second_gauss_area_coord = NULL;
    3537         double *third_gauss_area_coord  = NULL;
    3538         double *gauss_weights           = NULL;
    3539         double  gauss_weight;
    3540         double  gauss_l1l2l3[3];
     3510        int     ig;
     3511        GaussTria *gauss=NULL;
    35413512
    35423513        /* matrices: */
    … …  
    35593530        double    alpha2;
    35603531        double    alpha2_list[numvertices];                                       //TO BE DELETED
    3561         double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     3532        double    gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    35623533
    35633534        double MAXSLOPE=.06; // 6 %
    … …  
    35983569        /*COMPUT alpha2_list (TO BE DELETED)*/
    35993570        for(i=0;i<numvertices;i++){
    3600                 friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
    3601         }
    3602 
    3603         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3604         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     3571                friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum);
     3572        }
    36053573
    36063574        /* Start  looping on the number of gaussian points: */
    3607         for (ig=0; ig<num_gauss; ig++){
    3608                 /*Pick up the gaussian point: */
    3609                 gauss_weight=*(gauss_weights+ig);
    3610                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3611                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3612                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3575        gauss=new GaussTria(2);
     3576        for (ig=gauss->begin();ig<gauss->end();ig++){
     3577
     3578                gauss->GaussPoint(ig);
    36133579
    36143580                /*Friction: */
    3615                 // friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    3616                 TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
     3581                // friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
     3582                TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss); // TO BE DELETED
    36173583
    36183584                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    36193585                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    3620                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     3586                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    36213587                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    36223588
    … …  
    36263592
    36273593                /* Get Jacobian determinant: */
    3628                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     3594                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    36293595
    36303596                /*Get L matrix: */
    3631                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     3597                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    36323598
    36333599               
    3634                 DL_scalar=alpha2*gauss_weight*Jdet;
     3600                DL_scalar=alpha2*gauss->weight*Jdet;
    36353601                for (i=0;i<2;i++){
    36363602                        DL[i][i]=DL_scalar;
    … …  
    36503616        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    36513617
    3652         xfree((void**)&first_gauss_area_coord);
    3653         xfree((void**)&second_gauss_area_coord);
    3654         xfree((void**)&third_gauss_area_coord);
    3655         xfree((void**)&gauss_weights);
     3618        /*Clean up and return*/
     3619        delete gauss;
     3620        delete friction;
    36563621        xfree((void**)&doflist);
    3657         delete friction;
    36583622}       
    36593623/*}}}*/

  • issm/trunk/src/c/objects/Gauss/GaussTria.cpp

    r5630 r5637  
    102102}
    103103/*}}}*/
     104/*FUNCTION GaussTria::GaussCenter{{{1*/
     105void GaussTria::GaussCenter(void){
     106
     107        /*update static arrays*/
     108        coord1=ONETHIRD;
     109        coord2=ONETHIRD;
     110        coord3=ONETHIRD;
     111
     112}
     113/*}}}*/
    104114/*FUNCTION GaussTria::GaussPoint{{{1*/
    105115void GaussTria::GaussPoint(int ig){

  • issm/trunk/src/c/objects/Gauss/GaussTria.h

    r5630 r5637  
    4040                void GaussPoint(int ig);
    4141                void GaussVertex(int iv);
     42                void GaussCenter(void);
    4243};
    4344#endif  /* _GAUSSTRIA_H_ */
Note: See TracChangeset for help on using the changeset viewer.