Changeset 5636


Ignore:
Timestamp:
08/31/10 13:44:43 (15 years ago)
Author:
Mathieu Morlighem
Message:

more GaussTria

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

Legend:

Unmodified
Added
Removed

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

    r5635 r5636  
    19311931
    19321932        /* gaussian points: */
    1933         int     num_gauss,ig;
    1934         double* first_gauss_area_coord  =  NULL;
    1935         double* second_gauss_area_coord =  NULL;
    1936         double* third_gauss_area_coord  =  NULL;
    1937         double* gauss_weights           =  NULL;
    1938         double  gauss_weight;
    1939         double  gauss_l1l2l3[3];
     1933        int     ig;
     1934        GaussTria *gauss=NULL;
    19401935
    19411936        /* parameters: */
    … …  
    20082003        }
    20092004
    2010         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2011         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2012 
    20132005        /* Start  looping on the number of gaussian points: */
    2014         for (ig=0; ig<num_gauss; ig++){
    2015                 /*Pick up the gaussian point: */
    2016                 gauss_weight=*(gauss_weights+ig);
    2017                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2018                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2019                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    2020 
    2021                 /* Get Jacobian determinant: */
    2022                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2023 
    2024                 /*Compute misfit at gaussian point: */
    2025                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    2026 
    2027                 /*compute SurfaceRelVelMisfit*/
    2028                 Jelem+=misfit*Jdet*gauss_weight;
    2029         }
    2030 
    2031         xfree((void**)&first_gauss_area_coord);
    2032         xfree((void**)&second_gauss_area_coord);
    2033         xfree((void**)&third_gauss_area_coord);
    2034         xfree((void**)&gauss_weights);
    2035 
    2036         /*Return: */
     2006        gauss=new GaussTria(2);
     2007        for (ig=gauss->begin();ig<gauss->end();ig++){
     2008
     2009                gauss->GaussPoint(ig);
     2010
     2011                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2012                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
     2013                Jelem+=misfit*Jdet*gauss->weight;
     2014        }
     2015
     2016        /*clean up and Return: */
     2017        delete gauss;
    20372018        return Jelem;
    20382019}
    … …  
    20622043
    20632044        /* gaussian points: */
    2064         int     num_gauss,ig;
    2065         double* first_gauss_area_coord  =  NULL;
    2066         double* second_gauss_area_coord =  NULL;
    2067         double* third_gauss_area_coord  =  NULL;
    2068         double* gauss_weights           =  NULL;
    2069         double  gauss_weight;
    2070         double  gauss_l1l2l3[3];
     2045        int     ig;
     2046        GaussTria *gauss=NULL;
    20712047
    20722048        /* parameters: */
    … …  
    20832059        /*inputs: */
    20842060        bool onwater;
    2085        
    20862061        double  meanvel, epsvel;
    20872062
    … …  
    21372112        }
    21382113
    2139         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2140         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2141 
    21422114        /* Start  looping on the number of gaussian points: */
    2143         for (ig=0; ig<num_gauss; ig++){
    2144                 /*Pick up the gaussian point: */
    2145                 gauss_weight=*(gauss_weights+ig);
    2146                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2147                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2148                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    2149 
    2150                 /* Get Jacobian determinant: */
    2151                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2152 
    2153                 /*Compute misfit at gaussian point: */
    2154                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    2155 
    2156                 /*compute SurfaceLogVelMisfit*/
    2157                 Jelem+=misfit*Jdet*gauss_weight;
    2158         }
    2159 
    2160         xfree((void**)&first_gauss_area_coord);
    2161         xfree((void**)&second_gauss_area_coord);
    2162         xfree((void**)&third_gauss_area_coord);
    2163         xfree((void**)&gauss_weights);
    2164                
    2165         /*Return: */
     2115        gauss=new GaussTria(2);
     2116        for (ig=gauss->begin();ig<gauss->end();ig++){
     2117
     2118                gauss->GaussPoint(ig);
     2119
     2120                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2121                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
     2122                Jelem+=misfit*Jdet*gauss->weight;
     2123        }
     2124
     2125        /*clean-up and Return: */
     2126        delete gauss;
    21662127        return Jelem;
    21672128}
    … …  
    21912152
    21922153        /* gaussian points: */
    2193         int     num_gauss,ig;
    2194         double* first_gauss_area_coord  =  NULL;
    2195         double* second_gauss_area_coord =  NULL;
    2196         double* third_gauss_area_coord  =  NULL;
    2197         double* gauss_weights           =  NULL;
    2198         double  gauss_weight;
    2199         double  gauss_l1l2l3[3];
     2154        int     ig;
     2155        GaussTria *gauss=NULL;
    22002156
    22012157        /* parameters: */
    … …  
    22682224        }
    22692225
    2270         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2271         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2272 
    22732226        /* Start  looping on the number of gaussian points: */
    2274         for (ig=0; ig<num_gauss; ig++){
    2275                 /*Pick up the gaussian point: */
    2276                 gauss_weight=*(gauss_weights+ig);
    2277                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2278                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2279                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    2280 
    2281                 /* Get Jacobian determinant: */
    2282                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2283 
    2284                 /*Compute misfit at gaussian point: */
    2285                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    2286 
    2287                 /*compute SurfaceLogVxVyMisfit*/
    2288                 Jelem+=misfit*Jdet*gauss_weight;
    2289         }
    2290 
    2291         xfree((void**)&first_gauss_area_coord);
    2292         xfree((void**)&second_gauss_area_coord);
    2293         xfree((void**)&third_gauss_area_coord);
    2294         xfree((void**)&gauss_weights);
    2295                
    2296         /*Return: */
     2227        gauss=new GaussTria(2);
     2228        for (ig=gauss->begin();ig<gauss->end();ig++){
     2229
     2230                gauss->GaussPoint(ig);
     2231
     2232                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2233                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
     2234                Jelem+=misfit*Jdet*gauss->weight;
     2235        }
     2236
     2237        /*clean-up and Return: */
     2238        delete gauss;
    22972239        return Jelem;
    22982240}
    … …  
    23222264
    23232265        /* gaussian points: */
    2324         int     num_gauss,ig;
    2325         double* first_gauss_area_coord  =  NULL;
    2326         double* second_gauss_area_coord =  NULL;
    2327         double* third_gauss_area_coord  =  NULL;
    2328         double* gauss_weights           =  NULL;
    2329         double  gauss_weight;
    2330         double  gauss_l1l2l3[3];
     2266        int     ig;
     2267        GaussTria *gauss=NULL;
    23312268
    23322269        /* parameters: */
    … …  
    23982335        }
    23992336
    2400         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2401         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2402 
    24032337        /* Start  looping on the number of gaussian points: */
    2404         for (ig=0; ig<num_gauss; ig++){
    2405                 /*Pick up the gaussian point: */
    2406                 gauss_weight=*(gauss_weights+ig);
    2407                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2408                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2409                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    2410 
    2411                 /* Get Jacobian determinant: */
    2412                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2413 
    2414                 /*Compute misfit at gaussian point: */
    2415                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    2416 
    2417                 /*compute SurfaceAverageVelMisfit*/
    2418                 Jelem+=misfit*Jdet*gauss_weight;
    2419         }
    2420 
    2421         xfree((void**)&first_gauss_area_coord);
    2422         xfree((void**)&second_gauss_area_coord);
    2423         xfree((void**)&third_gauss_area_coord);
    2424         xfree((void**)&gauss_weights);
    2425 
    2426         /*Return: */
     2338        gauss=new GaussTria(2);
     2339        for (ig=gauss->begin();ig<gauss->end();ig++){
     2340
     2341                gauss->GaussPoint(ig);
     2342
     2343                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2344                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
     2345                Jelem+=misfit*Jdet*gauss->weight;
     2346        }
     2347
     2348        /*clean-up and Return: */
     2349                delete gauss;
    24272350        return Jelem;
    24282351}
    … …  
    24342357        int row;
    24352358        int vertices_ids[3];
    2436 
    24372359
    24382360        /*recover pointer: */
    … …  
    27792701
    27802702        /* gaussian points: */
    2781         int     num_gauss,ig;
    2782         double* first_gauss_area_coord  =  NULL;
    2783         double* second_gauss_area_coord =  NULL;
    2784         double* third_gauss_area_coord  =  NULL;
    2785         double* gauss_weights           =  NULL;
    2786         double  gauss_weight;
    2787         double  gauss_l1l2l3[3];
     2703        int     ig;
     2704        GaussTria *gauss=NULL;
    27882705
    27892706        /* matrices: */
    … …  
    28392756        if(artdiff){
    28402757                //Get the Jacobian determinant
    2841                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    2842                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     2758                gauss=new GaussTria();
     2759                gauss->GaussCenter();
     2760                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     2761                delete gauss;
    28432762
    28442763                //Build K matrix (artificial diffusivity matrix)
    … …  
    28502769        }
    28512770
    2852         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2853         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2854 
    28552771        /* Start  looping on the number of gaussian points: */
    2856         for (ig=0; ig<num_gauss; ig++){
    2857                 /*Pick up the gaussian point: */
    2858                 gauss_weight=*(gauss_weights+ig);
    2859                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2860                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2861                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2772        gauss=new GaussTria(2);
     2773        for (ig=gauss->begin();ig<gauss->end();ig++){
     2774
     2775                gauss->GaussPoint(ig);
    28622776
    28632777                /* Get Jacobian determinant: */
    2864                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     2778                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    28652779
    28662780                /*Get B  and B prime matrix: */
    2867                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    2868                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     2781                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     2782                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    28692783
    28702784                //Get vx, vy and their derivatives at gauss point
    2871                 vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);
    2872                 vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);
    2873 
    2874                 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
    2875                 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     2785                vxaverage_input->GetParameterValue(&vx,gauss);
     2786                vyaverage_input->GetParameterValue(&vy,gauss);
     2787
     2788                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     2789                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    28762790
    28772791                dvxdx=dvx[0];
    28782792                dvydy=dvy[1];
    28792793
    2880                 DL_scalar=gauss_weight*Jdettria;
     2794                DL_scalar=gauss->weight*Jdettria;
    28812795
    28822796                //Create DL and DLprime matrix
    … …  
    29242838        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    29252839
    2926         xfree((void**)&first_gauss_area_coord);
    2927         xfree((void**)&second_gauss_area_coord);
    2928         xfree((void**)&third_gauss_area_coord);
    2929         xfree((void**)&gauss_weights);
     2840        /*Clean up*/
     2841        delete gauss;
    29302842        xfree((void**)&doflist);
    29312843}
    … …  
    29452857
    29462858        /* gaussian points: */
    2947         int     num_gauss,ig;
    2948         double* first_gauss_area_coord  =  NULL;
    2949         double* second_gauss_area_coord =  NULL;
    2950         double* third_gauss_area_coord  =  NULL;
    2951         double* gauss_weights           =  NULL;
    2952         double  gauss_weight;
    2953         double  gauss_l1l2l3[3];
     2859        int     ig;
     2860        GaussTria *gauss=NULL;
    29542861
    29552862        /* matrices: */
    … …  
    29772884        this->parameters->FindParam(&dim,DimEnum);
    29782885
    2979         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    2980         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    2981 
    29822886        /*Retrieve all inputs we will be needed*/
    29832887        if(dim==2){
    … …  
    29872891
    29882892        /* Start  looping on the number of gaussian points: */
    2989         for (ig=0; ig<num_gauss; ig++){
    2990                 /*Pick up the gaussian point: */
    2991                 gauss_weight=*(gauss_weights+ig);
    2992                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    2993                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    2994                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2893        gauss=new GaussTria(2);
     2894        for (ig=gauss->begin();ig<gauss->end();ig++){
     2895
     2896                gauss->GaussPoint(ig);
    29952897
    29962898                /* Get Jacobian determinant: */
    2997                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     2899                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    29982900
    29992901                /*Get B  and B prime matrix: */
    30002902                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    3001                 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
    3002                 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     2903                GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     2904                GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
    30032905
    30042906                //Get vx, vy and their derivatives at gauss point
    3005                 vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);
    3006                 vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);
    3007 
    3008                 DL_scalar=-gauss_weight*Jdettria;
     2907                vx_input->GetParameterValue(&vx,gauss);
     2908                vy_input->GetParameterValue(&vy,gauss);
     2909
     2910                DL_scalar=-gauss->weight*Jdettria;
    30092911
    30102912                DLprime[0][0]=DL_scalar*vx;
    … …  
    30252927        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    30262928
    3027         xfree((void**)&first_gauss_area_coord);
    3028         xfree((void**)&second_gauss_area_coord);
    3029         xfree((void**)&third_gauss_area_coord);
    3030         xfree((void**)&gauss_weights);
     2929        /*Clean up and return*/
     2930        delete gauss;
    30312931        xfree((void**)&doflist);
    30322932}
    … …  
    30462946
    30472947        /* gaussian points: */
    3048         int     num_gauss,ig;
    3049         double* first_gauss_area_coord  =  NULL;
    3050         double* second_gauss_area_coord =  NULL;
    3051         double* third_gauss_area_coord  =  NULL;
    3052         double* gauss_weights           =  NULL;
    3053         double  gauss_weight;
    3054         double  gauss_l1l2l3[3];
     2948        int     ig;
     2949        GaussTria *gauss=NULL;
    30552950
    30562951        /* matrices: */
    … …  
    31293024        if(artdiff){
    31303025                //Get the Jacobian determinant
    3131                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    3132                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     3026                gauss=new GaussTria();
     3027                gauss->GaussCenter();
     3028                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     3029                delete gauss;
    31333030
    31343031                //Build K matrix (artificial diffusivity matrix)
    … …  
    31403037        }
    31413038
    3142         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3143         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    3144 
    31453039        /* Start  looping on the number of gaussian points: */
    3146         for (ig=0; ig<num_gauss; ig++){
    3147                 /*Pick up the gaussian point: */
    3148                 gauss_weight=*(gauss_weights+ig);
    3149                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3150                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3151                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     3040        gauss=new GaussTria(2);
     3041        for (ig=gauss->begin();ig<gauss->end();ig++){
     3042
     3043                gauss->GaussPoint(ig);
    31523044
    31533045                /* Get Jacobian determinant: */
    3154                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     3046                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    31553047
    31563048                /*Get B  and B prime matrix: */
    3157                 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    3158                 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     3049                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     3050                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    31593051
    31603052                //Get vx, vy and their derivatives at gauss point
    3161                 vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]);
    3162                 vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]);
    3163 
    3164                 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
    3165                 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     3053                vxaverage_input->GetParameterValue(&vx,gauss);
     3054                vyaverage_input->GetParameterValue(&vy,gauss);
     3055
     3056                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     3057                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
    31663058
    31673059                dvxdx=dvx[0];
    31683060                dvydy=dvy[1];
    31693061
    3170                 DL_scalar=gauss_weight*Jdettria;
     3062                DL_scalar=gauss->weight*Jdettria;
    31713063
    31723064                DLprime[0][0]=DL_scalar*nx;
    … …  
    32043096        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    32053097
    3206         xfree((void**)&first_gauss_area_coord);
    3207         xfree((void**)&second_gauss_area_coord);
    3208         xfree((void**)&third_gauss_area_coord);
    3209         xfree((void**)&gauss_weights);
     3098        /*Clean up and return*/
     3099        delete gauss;
    32103100        xfree((void**)&doflist);
    32113101}
    … …  
    32243114
    32253115        /* gaussian points: */
    3226         int     num_gauss,ig;
    3227         double* first_gauss_area_coord  =  NULL;
    3228         double* second_gauss_area_coord =  NULL;
    3229         double* third_gauss_area_coord  =  NULL;
    3230         double* gauss_weights           =  NULL;
    3231         double  gauss_weight;
    3232         double  gauss_l1l2l3[3];
     3116        int     ig;
     3117        GaussTria *gauss=NULL;
    32333118
    32343119        /* material data: */
    … …  
    32813166        GetDofList(&doflist,MacAyealApproximationEnum);
    32823167
    3283         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    3284         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    3285 
    32863168        /*Retrieve all inputs we will be needing: */
    32873169        thickness_input=inputs->GetInput(ThicknessEnum);
    … …  
    32923174
    32933175        /* Start  looping on the number of gaussian points: */
    3294         for (ig=0; ig<num_gauss; ig++){
    3295                 /*Pick up the gaussian point: */
    3296                 gauss_weight=*(gauss_weights+ig);
    3297                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    3298                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    3299                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    3300 
     3176        gauss=new GaussTria(2);
     3177        for (ig=gauss->begin();ig<gauss->end();ig++){
     3178
     3179                gauss->GaussPoint(ig);
    33013180
    33023181                /*Compute thickness at gaussian point: */
    3303                 thickness_input->GetParameterValue(&thickness, gauss_l1l2l3);
     3182                thickness_input->GetParameterValue(&thickness, gauss);
    33043183
    33053184                /*Get strain rate from velocity: */
    3306                 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input);
    3307                 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,vxold_input,vyold_input);
     3185                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3186                this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    33083187
    33093188                /*Get viscosity: */
    … …  
    33123191
    33133192                /* Get Jacobian determinant: */
    3314                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     3193                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    33153194
    33163195                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    33173196                        onto this scalar matrix, so that we win some computational time: */
    33183197                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    3319                 D_scalar=2*newviscosity*thickness*gauss_weight*Jdet;
     3198                D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
    33203199
    33213200                for (i=0;i<3;i++){
    … …  
    33243203
    33253204                /*Get B and Bprime matrices: */
    3326                 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    3327                 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     3205                GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
     3206                GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
    33283207
    33293208                /*  Do the triple product tB*D*Bprime: */
    … …  
    33463225        }
    33473226
    3348         xfree((void**)&first_gauss_area_coord);
    3349         xfree((void**)&second_gauss_area_coord);
    3350         xfree((void**)&third_gauss_area_coord);
    3351         xfree((void**)&gauss_weights);
     3227
     3228        /*Clean up and return*/
     3229        delete gauss;
    33523230        xfree((void**)&doflist);
    33533231}
    … …  
    69496827/*}}}1*/
    69506828/*FUNCTION Tria::SurfaceNormal{{{1*/
    6951 
    69526829void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
    69536830
    … …  
    69726849        *(surface_normal+1)=normal[1]/normal_norm;
    69736850        *(surface_normal+2)=normal[2]/normal_norm;
    6974 
    6975 }
    6976 /*}}}*/
     6851}
     6852/*}}}*/

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

    r5631 r5636  
    149149}
    150150/*}}}*/
     151/*FUNCTION TriaRef::GetBPrognostic{{{1*/
     152void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss){
     153        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     154         * For grid i, Bi can be expressed in the actual coordinate system
     155         * by:
     156         *       Bi=[ h ]
     157         *          [ h ]
     158         * where h is the interpolation function for grid i.
     159         *
     160         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
     161         */
     162
     163        int i;
     164        const int NDOF1=1;
     165        const int numgrids=3;
     166
     167        double l1l2l3[numgrids];
     168
     169
     170        /*Get dh1dh2dh3 in actual coordinate system: */
     171        GetNodalFunctions(&l1l2l3[0],gauss);
     172
     173        /*Build B_prog: */
     174        for (i=0;i<numgrids;i++){
     175                *(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
     176                *(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
     177        }
     178}
     179/*}}}*/
    151180/*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
    152181void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss){
    … …  
    184213}
    185214/*}}}*/
     215/*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
     216void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){
     217
     218        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     219         * For grid i, Bi' can be expressed in the actual coordinate system
     220         * by:
     221         *       Bi_prime=[ 2*dh/dx    dh/dy ]
     222         *                [   dh/dx  2*dh/dy ]
     223         *                [   dh/dy    dh/dx ]
     224         * where h is the interpolation function for grid i.
     225         *
     226         * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
     227         */
     228
     229        int i;
     230        const int NDOF2=2;
     231        const int numgrids=3;
     232
     233        /*Same thing in the actual coordinate system: */
     234        double dh1dh3[NDOF2][numgrids];
     235
     236        /*Get dh1dh2dh3 in actual coordinates system : */
     237        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     238
     239        /*Build B': */
     240        for (i=0;i<numgrids;i++){
     241                *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i];
     242                *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i];
     243                *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i];
     244                *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];
     245                *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i];
     246                *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i];
     247        }
     248}
     249/*}}}*/
    186250/*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
    187251void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss){
     252        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     253         * For grid i, Bi' can be expressed in the actual coordinate system
     254         * by:
     255         *       Bi_prime=[ dh/dx ]
     256         *                [ dh/dy ]
     257         * where h is the interpolation function for grid i.
     258         *
     259         * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
     260         */
     261
     262        int i;
     263        const int NDOF1=1;
     264        const int NDOF2=2;
     265        const int numgrids=3;
     266
     267        /*Same thing in the actual coordinate system: */
     268        double dh1dh3[NDOF2][numgrids];
     269
     270        /*Get dh1dh2dh3 in actual coordinates system : */
     271        GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
     272
     273        /*Build B': */
     274        for (i=0;i<numgrids;i++){
     275                *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i];
     276                *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i];
     277        }
     278}
     279/*}}}*/
     280/*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
     281void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss){
    188282        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    189283         * For grid i, Bi' can be expressed in the actual coordinate system

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

    r5631 r5636  
    2828                void GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss);
    2929                void GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss);
     30                void GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss);
    3031                void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);
     32                void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss);
    3133                void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);
     34                void GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss);
    3235                void GetL(double* L, double* xyz_list,double* gauss,int numdof);
    3336                void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
Note: See TracChangeset for help on using the changeset viewer.