Changeset 5166


Ignore:
Timestamp:
08/11/10 16:21:07 (15 years ago)
Author:
seroussi
Message:

moved Horiz to MacAyeal and Pattyn when possible

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

Legend:

Unmodified
Added
Removed

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

    r5158 r5166  
    744744        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    745745        if (analysis_type==DiagnosticHorizAnalysisEnum){
    746                 CreateKMatrixDiagnosticHoriz( Kgg);
     746                int approximation;
     747                inputs->GetParameterValue(&approximation,ApproximationEnum);
     748                if(approximation==MacAyealApproximationEnum){
     749                        CreateKMatrixDiagnosticMacAyeal( Kgg);
     750                }
     751                else if(approximation==PattynApproximationEnum){
     752                        CreateKMatrixDiagnosticPattyn( Kgg);
     753                }
     754                else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    747755        }
    748756        else if (analysis_type==AdjointHorizAnalysisEnum){
    749757                /*Same as diagnostic*/
    750                 CreateKMatrixDiagnosticHoriz( Kgg);
     758                int approximation;
     759                inputs->GetParameterValue(&approximation,ApproximationEnum);
     760                if(approximation==MacAyealApproximationEnum){
     761                        CreateKMatrixDiagnosticMacAyeal( Kgg);
     762                }
     763                else if(approximation==PattynApproximationEnum){
     764                        CreateKMatrixDiagnosticPattyn( Kgg);
     765                }
     766                else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    751767        }
    752768        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    798814        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    799815        if (analysis_type==DiagnosticHorizAnalysisEnum){
    800                 CreatePVectorDiagnosticHoriz( pg);
     816                int approximation;
     817                inputs->GetParameterValue(&approximation,ApproximationEnum);
     818                if(approximation==MacAyealApproximationEnum){
     819                        CreatePVectorDiagnosticMacAyeal( pg);
     820                }
     821                else if(approximation==PattynApproximationEnum){
     822                        CreatePVectorDiagnosticPattyn( pg);
     823                }
     824                else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    801825        }
    802826        else if (analysis_type==AdjointHorizAnalysisEnum){
    803                 CreatePVectorAdjointHoriz( pg);
     827                int approximation;
     828                inputs->GetParameterValue(&approximation,ApproximationEnum);
     829                if(approximation==MacAyealApproximationEnum){
     830                        CreatePVectorAdjointMacAyeal( pg);
     831                }
     832                else if(approximation==PattynApproximationEnum){
     833                        CreatePVectorAdjointPattyn( pg);
     834                }
     835                else ISSMERROR("Approximation %s not supported yet",EnumToString(approximation));
    804836        }
    805837        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    19151947                                this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum));
    19161948                        }
    1917                         else ISSMERROR("Approximation type %i (%s) not supported yet",*(iomodel->elements_type+2*index+0),EnumToString(*(iomodel->elements_type+2*index+0)));
     1949                        else{
     1950                                ISSMERROR("Approximation type %s not supported yet",EnumToString(*(iomodel->elements_type+2*index+0)));
     1951                        }
    19181952                }
    19191953        }
    … …  
    20412075        return;
    20422076
    2043 }
    2044 /*}}}*/
    2045 /*FUNCTION Penta::CreateKMatrixDiagnosticHoriz {{{1*/
    2046 void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg){
    2047 
    2048         /* local declarations */
    2049         int             i,j;
    2050 
    2051         /* node data: */
    2052         const int    numgrids=6;
    2053         const int    numdof=2*numgrids;
    2054         double       xyz_list[numgrids][3];
    2055         int*         doflist=NULL;
    2056 
    2057         /* 3d gaussian points: */
    2058         int     num_gauss,ig;
    2059         double* first_gauss_area_coord  =  NULL;
    2060         double* second_gauss_area_coord =  NULL;
    2061         double* third_gauss_area_coord  =  NULL;
    2062         double* fourth_gauss_vert_coord  =  NULL;
    2063         double* area_gauss_weights           =  NULL;
    2064         double* vert_gauss_weights           =  NULL;
    2065         int     ig1,ig2;
    2066         double  gauss_weight1,gauss_weight2;
    2067         double  gauss_coord[4];
    2068         int     order_area_gauss;
    2069         int     num_vert_gauss;
    2070         int     num_area_gauss;
    2071         double  gauss_weight;
    2072 
    2073         /* 2d gaussian point: */
    2074         int     num_gauss2d;
    2075         double* first_gauss_area_coord2d  =  NULL;
    2076         double* second_gauss_area_coord2d =  NULL;
    2077         double* third_gauss_area_coord2d  =  NULL;
    2078         double* gauss_weights2d=NULL;
    2079         double  gauss_l1l2l3[3];
    2080 
    2081         /* material data: */
    2082         double viscosity; //viscosity
    2083         double oldviscosity; //viscosity
    2084         double newviscosity; //viscosity
    2085 
    2086         /* strain rate: */
    2087         double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2088         double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2089 
    2090         /* matrices: */
    2091         double B[5][numdof];
    2092         double Bprime[5][numdof];
    2093         double L[2][numdof];
    2094         double D[5][5]={0.0};            // material matrix, simple scalar matrix.
    2095         double D_scalar;
    2096         double DL[2][2]={0.0}; //for basal drag
    2097         double DL_scalar;
    2098 
    2099         /* local element matrices: */
    2100         double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
    2101 
    2102         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    2103         double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    2104         double Jdet;
    2105 
    2106         /*slope: */
    2107         double  slope[2]={0.0};
    2108         double  slope_magnitude;
    2109 
    2110         /*friction: */
    2111         double  alpha2_list[3];
    2112         double  alpha2;
    2113 
    2114         double MAXSLOPE=.06; // 6 %
    2115         double MOUNTAINKEXPONENT=10;
    2116 
    2117         /*parameters: */
    2118         double viscosity_overshoot;
    2119 
    2120         /*Collapsed formulation: */
    2121         Tria*  tria=NULL;
    2122 
    2123         /*inputs: */
    2124         bool onwater;
    2125         bool onbed;
    2126         bool shelf;
    2127         int  approximation;
    2128         Input* vx_input=NULL;
    2129         Input* vy_input=NULL;
    2130         Input* vxold_input=NULL;
    2131         Input* vyold_input=NULL;
    2132 
    2133         /*retrieve inputs :*/
    2134         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    2135         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    2136         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    2137         inputs->GetParameterValue(&approximation,ApproximationEnum);
    2138 
    2139         /*retrieve some parameters: */
    2140         this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    2141 
    2142         /*If on water, skip stiffness: */
    2143         if(onwater)return;
    2144 
    2145         /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
    2146           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    2147           the stiffness matrix. */
    2148 
    2149         if ((approximation==MacAyealApproximationEnum) && (onbed==0)){
    2150                 /*This element should be collapsed, but this element is not on the bedrock, therefore all its
    2151                  * dofs have already been frozen! Do nothing: */
    2152                 return;
    2153         }
    2154         else if ((approximation==MacAyealApproximationEnum) && (onbed==1)){
    2155 
    2156                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    2157                  *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
    2158 
    2159                 /*Depth Averaging B*/
    2160                 this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
    2161 
    2162                 /*Call Tria function*/
    2163                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2164                 tria->CreateKMatrix(Kgg);
    2165                 delete tria->matice; delete tria;
    2166 
    2167                 /*Delete B averaged*/
    2168                 this->matice->inputs->DeleteInput(RheologyB2dEnum);
    2169 
    2170                 return;
    2171         }
    2172         else{
    2173 
    2174                 /*Implement standard penta element: */
    2175 
    2176                 /* Get node coordinates and dof list: */
    2177                 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    2178                 GetDofList(&doflist);
    2179 
    2180                 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    2181                   get tria gaussian points as well as segment gaussian points. For tria gaussian
    2182                   points, order of integration is 2, because we need to integrate the product tB*D*B'
    2183                   which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    2184                   points, same deal, which yields 3 gaussian points.*/
    2185 
    2186                 order_area_gauss=5;
    2187                 num_vert_gauss=5;
    2188 
    2189                 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    2190 
    2191                 /*Retrieve all inputs we will be needing: */
    2192                 vx_input=inputs->GetInput(VxEnum);
    2193                 vy_input=inputs->GetInput(VyEnum);
    2194                 vxold_input=inputs->GetInput(VxOldEnum);
    2195                 vyold_input=inputs->GetInput(VyOldEnum);
    2196 
    2197                 /* Start  looping on the number of gaussian points: */
    2198                 for (ig1=0; ig1<num_area_gauss; ig1++){
    2199                         for (ig2=0; ig2<num_vert_gauss; ig2++){
    2200 
    2201                                 /*Pick up the gaussian point: */
    2202                                 gauss_weight1=*(area_gauss_weights+ig1);
    2203                                 gauss_weight2=*(vert_gauss_weights+ig2);
    2204                                 gauss_weight=gauss_weight1*gauss_weight2;
    2205 
    2206 
    2207                                 gauss_coord[0]=*(first_gauss_area_coord+ig1);
    2208                                 gauss_coord[1]=*(second_gauss_area_coord+ig1);
    2209                                 gauss_coord[2]=*(third_gauss_area_coord+ig1);
    2210                                 gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    2211 
    2212 
    2213                                 /*Get strain rate from velocity: */
    2214                                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
    2215                                 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
    2216 
    2217                                 /*Get viscosity: */
    2218                                 matice->GetViscosity3d(&viscosity, &epsilon[0]);
    2219                                 matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    2220 
    2221                                 /*Get B and Bprime matrices: */
    2222                                 GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord);
    2223                                 GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
    2224 
    2225                                 /* Get Jacobian determinant: */
    2226                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    2227 
    2228                                 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
    2229                                   onto this scalar matrix, so that we win some computational time: */
    2230 
    2231                                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    2232                                 D_scalar=newviscosity*gauss_weight*Jdet;
    2233                                 for (i=0;i<5;i++){
    2234                                         D[i][i]=D_scalar;
    2235                                 }
    2236 
    2237                                 /*  Do the triple product tB*D*Bprime: */
    2238                                 TripleMultiply( &B[0][0],5,numdof,1,
    2239                                                         &D[0][0],5,5,0,
    2240                                                         &Bprime[0][0],5,numdof,0,
    2241                                                         &Ke_gg_gaussian[0][0],0);
    2242 
    2243                                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    2244                                 for( i=0; i<numdof; i++){
    2245                                         for(j=0;j<numdof;j++){
    2246                                                 Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2247                                         }
    2248                                 }
    2249                         } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    2250                 } //for (ig1=0; ig1<num_area_gauss; ig1++)
    2251 
    2252 
    2253                 /*Add Ke_gg to global matrix Kgg: */
    2254                 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    2255 
    2256                 //Deal with 2d friction at the bedrock interface
    2257                 if((onbed && !shelf)){
    2258 
    2259                         /*Build a tria element using the 3 grids of the base of the penta. Then use
    2260                          * the tria functionality to build a friction stiffness matrix on these 3
    2261                          * grids: */
    2262 
    2263                         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2264                         tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
    2265                         delete tria->matice; delete tria;
    2266                 }
    2267 
    2268         }
    2269 
    2270         xfree((void**)&first_gauss_area_coord);
    2271         xfree((void**)&second_gauss_area_coord);
    2272         xfree((void**)&third_gauss_area_coord);
    2273         xfree((void**)&fourth_gauss_vert_coord);
    2274         xfree((void**)&area_gauss_weights);
    2275         xfree((void**)&vert_gauss_weights);
    2276         xfree((void**)&first_gauss_area_coord2d);
    2277         xfree((void**)&second_gauss_area_coord2d);
    2278         xfree((void**)&third_gauss_area_coord2d);
    2279         xfree((void**)&gauss_weights2d);
    2280         xfree((void**)&doflist);
    22812077}
    22822078/*}}}*/
    … …  
    23642160
    23652161}
     2162/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal{{{1*/
     2163void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){
     2164       
     2165        /*Collapsed formulation: */
     2166        Tria*  tria=NULL;
     2167
     2168        /*inputs: */
     2169        bool onwater;
     2170        bool onbed;
     2171
     2172        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2173        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2174
     2175        /*If on water, skip stiffness: */
     2176        if(onwater)return;
     2177
     2178        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
     2179          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
     2180          the stiffness matrix. */
     2181
     2182        if (onbed==0){
     2183                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
     2184                 * dofs have already been frozen! Do nothing: */
     2185                return;
     2186        }
     2187        else if (onbed==1){
     2188
     2189                /*This element should be collapsed into a tria element at its base. Create this tria element,
     2190                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
     2191
     2192                /*Depth Averaging B*/
     2193                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
     2194
     2195                /*Call Tria function*/
     2196                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2197                tria->CreateKMatrix(Kgg);
     2198                delete tria->matice; delete tria;
     2199
     2200                /*Delete B averaged*/
     2201                this->matice->inputs->DeleteInput(RheologyB2dEnum);
     2202
     2203                return;
     2204        }
     2205}
     2206/*}}}*/
     2207/*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/
     2208void Penta::CreateKMatrixDiagnosticPattyn( Mat Kgg){
     2209
     2210        /* local declarations */
     2211        int             i,j;
     2212
     2213        /* node data: */
     2214        const int    numgrids=6;
     2215        const int    numdof=2*numgrids;
     2216        double       xyz_list[numgrids][3];
     2217        int*         doflist=NULL;
     2218
     2219        /* 3d gaussian points: */
     2220        int     num_gauss,ig;
     2221        double* first_gauss_area_coord  =  NULL;
     2222        double* second_gauss_area_coord =  NULL;
     2223        double* third_gauss_area_coord  =  NULL;
     2224        double* fourth_gauss_vert_coord  =  NULL;
     2225        double* area_gauss_weights           =  NULL;
     2226        double* vert_gauss_weights           =  NULL;
     2227        int     ig1,ig2;
     2228        double  gauss_weight1,gauss_weight2;
     2229        double  gauss_coord[4];
     2230        int     order_area_gauss;
     2231        int     num_vert_gauss;
     2232        int     num_area_gauss;
     2233        double  gauss_weight;
     2234
     2235        /* 2d gaussian point: */
     2236        int     num_gauss2d;
     2237        double* first_gauss_area_coord2d  =  NULL;
     2238        double* second_gauss_area_coord2d =  NULL;
     2239        double* third_gauss_area_coord2d  =  NULL;
     2240        double* gauss_weights2d=NULL;
     2241        double  gauss_l1l2l3[3];
     2242
     2243        /* material data: */
     2244        double viscosity; //viscosity
     2245        double oldviscosity; //viscosity
     2246        double newviscosity; //viscosity
     2247
     2248        /* strain rate: */
     2249        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     2250        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     2251
     2252        /* matrices: */
     2253        double B[5][numdof];
     2254        double Bprime[5][numdof];
     2255        double L[2][numdof];
     2256        double D[5][5]={0.0};            // material matrix, simple scalar matrix.
     2257        double D_scalar;
     2258        double DL[2][2]={0.0}; //for basal drag
     2259        double DL_scalar;
     2260
     2261        /* local element matrices: */
     2262        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
     2263        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
     2264        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
     2265        double Jdet;
     2266
     2267        /*slope: */
     2268        double  slope[2]={0.0};
     2269        double  slope_magnitude;
     2270
     2271        /*friction: */
     2272        double  alpha2_list[3];
     2273        double  alpha2;
     2274
     2275        double MAXSLOPE=.06; // 6 %
     2276        double MOUNTAINKEXPONENT=10;
     2277
     2278        /*parameters: */
     2279        double viscosity_overshoot;
     2280
     2281        /*Collapsed formulation: */
     2282        Tria*  tria=NULL;
     2283
     2284        /*inputs: */
     2285        bool onwater;
     2286        bool onbed;
     2287        bool shelf;
     2288        Input* vx_input=NULL;
     2289        Input* vy_input=NULL;
     2290        Input* vxold_input=NULL;
     2291        Input* vyold_input=NULL;
     2292
     2293        /*retrieve inputs :*/
     2294        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     2295        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2296        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     2297
     2298        /*retrieve some parameters: */
     2299        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
     2300
     2301        /*If on water, skip stiffness: */
     2302        if(onwater)return;
     2303
     2304        /*Implement standard penta element: */
     2305
     2306        /* Get node coordinates and dof list: */
     2307        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2308        GetDofList(&doflist);
     2309
     2310        /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     2311          get tria gaussian points as well as segment gaussian points. For tria gaussian
     2312          points, order of integration is 2, because we need to integrate the product tB*D*B'
     2313          which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     2314          points, same deal, which yields 3 gaussian points.*/
     2315
     2316        order_area_gauss=5;
     2317        num_vert_gauss=5;
     2318
     2319        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     2320
     2321        /*Retrieve all inputs we will be needing: */
     2322        vx_input=inputs->GetInput(VxEnum);
     2323        vy_input=inputs->GetInput(VyEnum);
     2324        vxold_input=inputs->GetInput(VxOldEnum);
     2325        vyold_input=inputs->GetInput(VyOldEnum);
     2326
     2327        /* Start  looping on the number of gaussian points: */
     2328        for (ig1=0; ig1<num_area_gauss; ig1++){
     2329                for (ig2=0; ig2<num_vert_gauss; ig2++){
     2330
     2331                        /*Pick up the gaussian point: */
     2332                        gauss_weight1=*(area_gauss_weights+ig1);
     2333                        gauss_weight2=*(vert_gauss_weights+ig2);
     2334                        gauss_weight=gauss_weight1*gauss_weight2;
     2335
     2336
     2337                        gauss_coord[0]=*(first_gauss_area_coord+ig1);
     2338                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
     2339                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
     2340                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
     2341
     2342
     2343                        /*Get strain rate from velocity: */
     2344                        this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
     2345                        this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
     2346
     2347                        /*Get viscosity: */
     2348                        matice->GetViscosity3d(&viscosity, &epsilon[0]);
     2349                        matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
     2350
     2351                        /*Get B and Bprime matrices: */
     2352                        GetBPattyn(&B[0][0], &xyz_list[0][0], gauss_coord);
     2353                        GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss_coord);
     2354
     2355                        /* Get Jacobian determinant: */
     2356                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
     2357
     2358                        /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
     2359                          onto this scalar matrix, so that we win some computational time: */
     2360
     2361                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
     2362                        D_scalar=newviscosity*gauss_weight*Jdet;
     2363                        for (i=0;i<5;i++){
     2364                                D[i][i]=D_scalar;
     2365                        }
     2366
     2367                        /*  Do the triple product tB*D*Bprime: */
     2368                        TripleMultiply( &B[0][0],5,numdof,1,
     2369                                                &D[0][0],5,5,0,
     2370                                                &Bprime[0][0],5,numdof,0,
     2371                                                &Ke_gg_gaussian[0][0],0);
     2372
     2373                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
     2374                        for( i=0; i<numdof; i++){
     2375                                for(j=0;j<numdof;j++){
     2376                                        Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2377                                }
     2378                        }
     2379                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     2380        } //for (ig1=0; ig1<num_area_gauss; ig1++)
     2381
     2382
     2383        /*Add Ke_gg to global matrix Kgg: */
     2384        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     2385
     2386        //Deal with 2d friction at the bedrock interface
     2387        if((onbed && !shelf)){
     2388
     2389                /*Build a tria element using the 3 grids of the base of the penta. Then use
     2390                 * the tria functionality to build a friction stiffness matrix on these 3
     2391                 * grids: */
     2392
     2393                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2394                tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
     2395                delete tria->matice; delete tria;
     2396        }
     2397
     2398        xfree((void**)&first_gauss_area_coord);
     2399        xfree((void**)&second_gauss_area_coord);
     2400        xfree((void**)&third_gauss_area_coord);
     2401        xfree((void**)&fourth_gauss_vert_coord);
     2402        xfree((void**)&area_gauss_weights);
     2403        xfree((void**)&vert_gauss_weights);
     2404        xfree((void**)&first_gauss_area_coord2d);
     2405        xfree((void**)&second_gauss_area_coord2d);
     2406        xfree((void**)&third_gauss_area_coord2d);
     2407        xfree((void**)&gauss_weights2d);
     2408        xfree((void**)&doflist);
     2409}
    23662410/*}}}*/
    23672411/*FUNCTION Penta::CreateKMatrixDiagnosticStokes {{{1*/
    … …  
    31083152}
    31093153/*}}}*/
     3154/*FUNCTION Penta::CreatePVectorAdjointMacAyeal{{{1*/
     3155void  Penta::CreatePVectorAdjointMacAyeal(Vec p_g){
     3156
     3157        int i;
     3158        Tria* tria=NULL;
     3159
     3160        /*inputs: */
     3161        bool onwater;
     3162        bool onbed;
     3163
     3164        /*retrieve inputs :*/
     3165        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     3166        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     3167
     3168        /*If on water, skip: */
     3169        if(onwater) return;
     3170
     3171        /*Bail out if this element if:
     3172         * -> Non MacAyeal and not on the surface
     3173         * -> MacAyeal(2d model) and not on bed) */
     3174        if (!onbed){
     3175                return;
     3176        }
     3177        else{
     3178
     3179                /*This element should be collapsed into a tria element at its base. Create this tria element,
     3180                 * and compute pe_g*/
     3181                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     3182                tria->CreatePVectorAdjointHoriz(p_g);
     3183                delete tria->matice; delete tria;
     3184                return;
     3185        }
     3186}
     3187/*}}}*/
     3188/*FUNCTION Penta::CreatePVectorAdjointPattyn{{{1*/
     3189void  Penta::CreatePVectorAdjointPattyn(Vec p_g){
     3190
     3191        int i;
     3192        Tria* tria=NULL;
     3193
     3194        /*inputs: */
     3195        bool onwater;
     3196        bool onsurface;
     3197
     3198        /*retrieve inputs :*/
     3199        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     3200        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     3201
     3202        /*If on water, skip: */
     3203        if(onwater) return;
     3204
     3205        /*Bail out if this element if:
     3206         * -> Non MacAyeal and not on the surface
     3207         * -> MacAyeal(2d model) and not on bed) */
     3208        if (!onsurface){
     3209                return;
     3210        }
     3211        else{
     3212
     3213                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     3214                tria->CreatePVectorAdjointHoriz(p_g);
     3215                delete tria->matice; delete tria;
     3216                return;
     3217        }
     3218}
     3219/*}}}*/
    31103220/*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/
    31113221void Penta::CreatePVectorBalancedthickness( Vec pg){
    … …  
    31803290}
    31813291/*}}}*/
    3182 /*FUNCTION Penta::CreatePVectorDiagnosticHoriz {{{1*/
    3183 void Penta::CreatePVectorDiagnosticHoriz( Vec pg){
     3292/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
     3293void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
     3294
     3295        int i,j,k;
     3296       
     3297        const int numgrids=6;
     3298        const int NDOF2=2;
     3299        const int numdofs=NDOF2*numgrids;
     3300        int*      doflist=NULL;
     3301        double    pe_g[numdofs]={0.0};
     3302        double    xyz_list[numgrids][3];
     3303        double    z_list[numgrids];
     3304        double    z_segment[2];
     3305        double    slope2,constant_part;
     3306        double    gauss[numdofs][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}};
     3307        int       node0,node1;
     3308        BeamRef*  beam=NULL;
     3309
     3310        /*gaussian points: */
     3311        int      num_gauss;
     3312        double*  segment_gauss_coord=NULL;
     3313        double   gauss_coord;
     3314        double*  gauss_weights=NULL;
     3315        double   gauss_weight;
     3316        int      ig;
     3317        double   slope[2];
     3318
     3319        double   z_g;
     3320        double   rho_ice,gravity,n,B;
     3321        double   Jdet;
     3322        double   ub,vb;
     3323        double   surface,thickness;
     3324
     3325        /*flags: */
     3326        bool onwater;
     3327        bool onbed;
     3328        bool onsurface;
     3329        int  connectivity[2];
     3330
     3331        /*Inputs*/
     3332        Input* thickness_input;
     3333        Input* surface_input;
     3334        Input* slopex_input;
     3335        Input* slopey_input;
     3336
     3337        /*recover some inputs: */
     3338        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     3339
     3340        /*If on water, skip: */
     3341        if(onwater)return;
     3342
     3343        /*recover doflist: */
     3344        GetDofList(&doflist);
     3345
     3346        /*recover some inputs: */
     3347        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     3348        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     3349
     3350        /*recover parameters: */
     3351        rho_ice=matpar->GetRhoIce();
     3352        gravity=matpar->GetG();
     3353        n=matice->GetN();
     3354        B=matice->GetB();
     3355
     3356        //Initilaize beam
     3357        beam=new BeamRef();
     3358
     3359        //recover extra inputs
     3360        thickness_input=inputs->GetInput(ThicknessEnum);
     3361        surface_input=inputs->GetInput(SurfaceEnum);
     3362        slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
     3363        slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
     3364
     3365        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3366        for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
     3367
     3368        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
     3369        num_gauss=3;
     3370        GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
     3371
     3372        /*Loop on the three segments*/
     3373        for(i=0;i<3;i++){
     3374                slopex_input->GetParameterValue(&slope[0], &gauss[i][0]);
     3375                slopey_input->GetParameterValue(&slope[1], &gauss[i][0]);
     3376                surface_input->GetParameterValue(&surface, &gauss[i][0]);
     3377                thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
     3378
     3379                //compute slope2 slopex and slopey
     3380                slope2=pow(slope[0],2)+pow(slope[1],2);
     3381
     3382                //%compute constant_part
     3383                constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
     3384
     3385                z_segment[0]=z_list[i];
     3386                z_segment[1]=z_list[i+3];
     3387
     3388                node0=i;
     3389                node1=i+3;
     3390
     3391                connectivity[0]=nodes[node0]->GetConnectivity();
     3392                connectivity[1]=nodes[node1]->GetConnectivity();
     3393
     3394                /*Loop on the Gauss points: */
     3395                for(ig=0;ig<num_gauss;ig++){
     3396
     3397                        //Pick up the gaussian point and its weight:
     3398                        gauss_weight=gauss_weights[ig];
     3399                        gauss_coord=segment_gauss_coord[ig];
     3400
     3401                        beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
     3402
     3403                        //Get Jacobian determinant:
     3404                        beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
     3405
     3406                        if (onsurface){
     3407                                for(j=0;j<NDOF2;j++){
     3408                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
     3409                                }
     3410                        }
     3411                        else{//connectivity is too large, should take only half on it
     3412                                for(j=0;j<NDOF2;j++){
     3413                                        pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
     3414                                }
     3415                        }
     3416                }
     3417
     3418                //Deal with lower surface
     3419                if (onbed){
     3420
     3421                        //compute ub
     3422                        constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
     3423                        ub=constant_part*slope[0];
     3424                        vb=constant_part*slope[1];
     3425
     3426                        //Add to pe:
     3427                        pe_g[2*node0]+=ub/(double)connectivity[0];
     3428                        pe_g[2*node0+1]+=vb/(double)connectivity[0];
     3429                }
     3430        }
     3431
     3432        /*Add pe_g to global vector pg: */
     3433        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     3434
     3435        /*Clean up */
     3436        delete beam;
     3437        xfree((void**)&gauss_weights);
     3438        xfree((void**)&segment_gauss_coord);
     3439        xfree((void**)&doflist);
     3440}
     3441/*}}}*/
     3442/*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{1*/
     3443void  Penta::CreatePVectorDiagnosticMacAyeal(Vec pg){
     3444
     3445        /*Spawning: */
     3446        Tria* tria=NULL;
     3447
     3448        /*inputs: */
     3449        bool onwater;
     3450        bool onbed;
     3451
     3452        /*retrieve inputs :*/
     3453        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     3454        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     3455
     3456        /*If on water, skip load: */
     3457        if(onwater)return;
     3458
     3459        /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the
     3460          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
     3461          the load vector. */
     3462
     3463        if (onbed==0){
     3464                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
     3465                 * dofs have already been frozen! Do nothing: */
     3466                return;
     3467        }
     3468        else if (onbed==1){
     3469
     3470                /*This element should be collapsed into a tria element at its base. Create this tria element,
     3471                 *and use its CreatePVector functionality to return an elementary load vector: */
     3472                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     3473                tria->CreatePVector(pg);
     3474                delete tria->matice; delete tria;
     3475                return;
     3476        }
     3477}
     3478/*}}}*/
     3479/*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{1*/
     3480void Penta::CreatePVectorDiagnosticPattyn( Vec pg){
    31843481
    31853482        int i,j;
    … …  
    32233520        double  pe_g_gaussian[numdof];
    32243521
    3225         /*Spawning: */
    3226         Tria* tria=NULL;
    3227 
    32283522        /*inputs: */
    32293523        bool onwater;
    3230         bool onbed;
    3231         int  approximation;
    32323524        Input* surface_input=NULL;
    32333525        Input* thickness_input=NULL;
    … …  
    32353527        /*retrieve inputs :*/
    32363528        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    3237         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    3238         inputs->GetParameterValue(&approximation,ApproximationEnum);
    32393529
    32403530        /*If on water, skip load: */
    32413531        if(onwater)return;
    32423532
    3243         /*Figure out if this pentaelem is Macayeal. If so, then bailout, except if it is at the
    3244           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    3245           the load vector. */
    3246 
    3247         if ((approximation==MacAyealApproximationEnum) && (onbed==0)){
    3248                 /*This element should be collapsed, but this element is not on the bedrock, therefore all its
    3249                  * dofs have already been frozen! Do nothing: */
    3250                 return;
    3251         }
    3252         else if ((approximation==MacAyealApproximationEnum) && (onbed==1)){
    3253 
    3254                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3255                  *and use its CreatePVector functionality to return an elementary load vector: */
    3256                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    3257                 tria->CreatePVector(pg);
    3258                 delete tria->matice; delete tria;
    3259                 return;
    3260         }
    3261         else{
    3262 
    3263                 /*Implement standard penta element: */
    3264 
    3265                 /* Get node coordinates and dof list: */
    3266                 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    3267                 GetDofList(&doflist);
    3268 
    3269                 /*Get gaussian points and weights :*/
    3270                 order_area_gauss=2;
    3271                 num_vert_gauss=3;
    3272 
    3273                 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    3274 
    3275                 /*Retrieve all inputs we will be needing: */
    3276                 thickness_input=inputs->GetInput(ThicknessEnum);
    3277                 surface_input=inputs->GetInput(SurfaceEnum);
    3278 
    3279                 /* Start  looping on the number of gaussian points: */
    3280                 for (ig1=0; ig1<num_area_gauss; ig1++){
    3281                         for (ig2=0; ig2<num_vert_gauss; ig2++){
    3282 
    3283                                 /*Pick up the gaussian point: */
    3284                                 gauss_weight1=*(area_gauss_weights+ig1);
    3285                                 gauss_weight2=*(vert_gauss_weights+ig2);
    3286                                 gauss_weight=gauss_weight1*gauss_weight2;
    3287 
    3288                                 gauss_coord[0]=*(first_gauss_area_coord+ig1);
    3289                                 gauss_coord[1]=*(second_gauss_area_coord+ig1);
    3290                                 gauss_coord[2]=*(third_gauss_area_coord+ig1);
    3291                                 gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
    3292 
    3293                                 /*Compute thickness at gaussian point: */
    3294                                 thickness_input->GetParameterValue(&thickness, gauss_coord);
    3295 
    3296                                 /*Compute slope at gaussian point: */
    3297                                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord);
    3298 
    3299                                 /* Get Jacobian determinant: */
    3300                                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    3301 
    3302                                 /*Get nodal functions: */
    3303                                 GetNodalFunctionsP1(l1l6, gauss_coord);
    3304 
    3305                                 /*Compute driving stress: */
    3306                                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    3307 
    3308                                 /*Build pe_g_gaussian vector: */
    3309                                 for (i=0;i<numgrids;i++){
    3310                                         for (j=0;j<NDOF2;j++){
    3311                                                 pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l6[i];
    3312                                         }
     3533        /*Implement standard penta element: */
     3534
     3535        /* Get node coordinates and dof list: */
     3536        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     3537        GetDofList(&doflist);
     3538
     3539        /*Get gaussian points and weights :*/
     3540        order_area_gauss=2;
     3541        num_vert_gauss=3;
     3542
     3543        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
     3544
     3545        /*Retrieve all inputs we will be needing: */
     3546        thickness_input=inputs->GetInput(ThicknessEnum);
     3547        surface_input=inputs->GetInput(SurfaceEnum);
     3548
     3549        /* Start  looping on the number of gaussian points: */
     3550        for (ig1=0; ig1<num_area_gauss; ig1++){
     3551                for (ig2=0; ig2<num_vert_gauss; ig2++){
     3552
     3553                        /*Pick up the gaussian point: */
     3554                        gauss_weight1=*(area_gauss_weights+ig1);
     3555                        gauss_weight2=*(vert_gauss_weights+ig2);
     3556                        gauss_weight=gauss_weight1*gauss_weight2;
     3557
     3558                        gauss_coord[0]=*(first_gauss_area_coord+ig1);
     3559                        gauss_coord[1]=*(second_gauss_area_coord+ig1);
     3560                        gauss_coord[2]=*(third_gauss_area_coord+ig1);
     3561                        gauss_coord[3]=*(fourth_gauss_vert_coord+ig2);
     3562
     3563                        /*Compute thickness at gaussian point: */
     3564                        thickness_input->GetParameterValue(&thickness, gauss_coord);
     3565
     3566                        /*Compute slope at gaussian point: */
     3567                        surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord);
     3568
     3569                        /* Get Jacobian determinant: */
     3570                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
     3571
     3572                        /*Get nodal functions: */
     3573                        GetNodalFunctionsP1(l1l6, gauss_coord);
     3574
     3575                        /*Compute driving stress: */
     3576                        driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
     3577
     3578                        /*Build pe_g_gaussian vector: */
     3579                        for (i=0;i<numgrids;i++){
     3580                                for (j=0;j<NDOF2;j++){
     3581                                        pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l6[i];
    33133582                                }
    3314 
    3315                                 /*Add pe_g_gaussian vector to pe_g: */
    3316                                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    3317 
    3318                         } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    3319                 } //for (ig1=0; ig1<num_area_gauss; ig1++)
    3320 
    3321         } //else if ((approximation==MacAyealApproximationEnum) && (onbed==1))
     3583                        }
     3584
     3585                        /*Add pe_g_gaussian vector to pe_g: */
     3586                        for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
     3587
     3588                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     3589        } //for (ig1=0; ig1<num_area_gauss; ig1++)
    33223590
    33233591        /*Add pe_g to global vector pg: */
    … …  
    33303598        xfree((void**)&area_gauss_weights);
    33313599        xfree((void**)&vert_gauss_weights);
    3332         xfree((void**)&doflist);
    3333 }
    3334 /*}}}*/
    3335 /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
    3336 void  Penta::CreatePVectorAdjointHoriz(Vec p_g){
    3337 
    3338         int i;
    3339         Tria* tria=NULL;
    3340 
    3341         /*inputs: */
    3342         bool onwater;
    3343         bool onsurface;
    3344         bool onbed;
    3345         int  approximation;
    3346 
    3347         /*retrieve inputs :*/
    3348         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    3349         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    3350         inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    3351         inputs->GetParameterValue(&approximation,ApproximationEnum);
    3352 
    3353         /*If on water, skip: */
    3354         if(onwater) return;
    3355 
    3356         /*Bail out if this element if:
    3357          * -> Non MacAyeal and not on the surface
    3358          * -> MacAyeal(2d model) and not on bed) */
    3359         if ((approximation!=MacAyealApproximationEnum && !onsurface) || (approximation==MacAyealApproximationEnum && !onbed)){
    3360                 return;
    3361         }
    3362         else if (approximation==MacAyealApproximationEnum){
    3363 
    3364                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    3365                  * and compute pe_g*/
    3366                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    3367                 tria->CreatePVectorAdjointHoriz(p_g);
    3368                 delete tria->matice; delete tria;
    3369                 return;
    3370         }
    3371         else{
    3372 
    3373                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    3374                 tria->CreatePVectorAdjointHoriz(p_g);
    3375                 delete tria->matice; delete tria;
    3376                 return;
    3377         }
    3378 }
    3379 /*}}}*/
    3380 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
    3381 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
    3382 
    3383         int i,j,k;
    3384        
    3385         const int numgrids=6;
    3386         const int NDOF2=2;
    3387         const int numdofs=NDOF2*numgrids;
    3388         int*      doflist=NULL;
    3389         double    pe_g[numdofs]={0.0};
    3390         double    xyz_list[numgrids][3];
    3391         double    z_list[numgrids];
    3392         double    z_segment[2];
    3393         double    slope2,constant_part;
    3394         double    gauss[numdofs][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}};
    3395         int       node0,node1;
    3396         BeamRef*  beam=NULL;
    3397 
    3398         /*gaussian points: */
    3399         int      num_gauss;
    3400         double*  segment_gauss_coord=NULL;
    3401         double   gauss_coord;
    3402         double*  gauss_weights=NULL;
    3403         double   gauss_weight;
    3404         int      ig;
    3405         double   slope[2];
    3406 
    3407         double   z_g;
    3408         double   rho_ice,gravity,n,B;
    3409         double   Jdet;
    3410         double   ub,vb;
    3411         double   surface,thickness;
    3412 
    3413         /*flags: */
    3414         bool onwater;
    3415         bool onbed;
    3416         bool onsurface;
    3417         int  connectivity[2];
    3418 
    3419         /*Inputs*/
    3420         Input* thickness_input;
    3421         Input* surface_input;
    3422         Input* slopex_input;
    3423         Input* slopey_input;
    3424 
    3425         /*recover some inputs: */
    3426         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    3427 
    3428         /*If on water, skip: */
    3429         if(onwater)return;
    3430 
    3431         /*recover doflist: */
    3432         GetDofList(&doflist);
    3433 
    3434         /*recover some inputs: */
    3435         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    3436         inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    3437 
    3438         /*recover parameters: */
    3439         rho_ice=matpar->GetRhoIce();
    3440         gravity=matpar->GetG();
    3441         n=matice->GetN();
    3442         B=matice->GetB();
    3443 
    3444         //Initilaize beam
    3445         beam=new BeamRef();
    3446 
    3447         //recover extra inputs
    3448         thickness_input=inputs->GetInput(ThicknessEnum);
    3449         surface_input=inputs->GetInput(SurfaceEnum);
    3450         slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
    3451         slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
    3452 
    3453         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    3454         for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
    3455 
    3456         //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    3457         num_gauss=3;
    3458         GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    3459 
    3460         /*Loop on the three segments*/
    3461         for(i=0;i<3;i++){
    3462                 slopex_input->GetParameterValue(&slope[0], &gauss[i][0]);
    3463                 slopey_input->GetParameterValue(&slope[1], &gauss[i][0]);
    3464                 surface_input->GetParameterValue(&surface, &gauss[i][0]);
    3465                 thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
    3466 
    3467                 //compute slope2 slopex and slopey
    3468                 slope2=pow(slope[0],2)+pow(slope[1],2);
    3469 
    3470                 //%compute constant_part
    3471                 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
    3472 
    3473                 z_segment[0]=z_list[i];
    3474                 z_segment[1]=z_list[i+3];
    3475 
    3476                 node0=i;
    3477                 node1=i+3;
    3478 
    3479                 connectivity[0]=nodes[node0]->GetConnectivity();
    3480                 connectivity[1]=nodes[node1]->GetConnectivity();
    3481 
    3482                 /*Loop on the Gauss points: */
    3483                 for(ig=0;ig<num_gauss;ig++){
    3484 
    3485                         //Pick up the gaussian point and its weight:
    3486                         gauss_weight=gauss_weights[ig];
    3487                         gauss_coord=segment_gauss_coord[ig];
    3488 
    3489                         beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
    3490 
    3491                         //Get Jacobian determinant:
    3492                         beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
    3493 
    3494                         if (onsurface){
    3495                                 for(j=0;j<NDOF2;j++){
    3496                                         pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
    3497                                 }
    3498                         }
    3499                         else{//connectivity is too large, should take only half on it
    3500                                 for(j=0;j<NDOF2;j++){
    3501                                         pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
    3502                                 }
    3503                         }
    3504                 }
    3505 
    3506                 //Deal with lower surface
    3507                 if (onbed){
    3508 
    3509                         //compute ub
    3510                         constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
    3511                         ub=constant_part*slope[0];
    3512                         vb=constant_part*slope[1];
    3513 
    3514                         //Add to pe:
    3515                         pe_g[2*node0]+=ub/(double)connectivity[0];
    3516                         pe_g[2*node0+1]+=vb/(double)connectivity[0];
    3517                 }
    3518         }
    3519 
    3520         /*Add pe_g to global vector pg: */
    3521         VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
    3522 
    3523         /*Clean up */
    3524         delete beam;
    3525         xfree((void**)&gauss_weights);
    3526         xfree((void**)&segment_gauss_coord);
    35273600        xfree((void**)&doflist);
    35283601}

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

    r5096 r5166  
    115115                void      CreateKMatrixBalancedthickness(Mat Kggg);
    116116                void      CreateKMatrixBalancedvelocities(Mat Kggg);
    117                 void      CreateKMatrixDiagnosticHoriz( Mat Kgg);
    118117                void      CreateKMatrixDiagnosticHutter( Mat Kgg);
     118                void      CreateKMatrixDiagnosticMacAyeal( Mat Kgg);
     119                void      CreateKMatrixDiagnosticPattyn( Mat Kgg);
    119120                void      CreateKMatrixDiagnosticStokes( Mat Kgg);
    120121                void      CreateKMatrixDiagnosticVert( Mat Kgg);
    … …  
    125126                void      CreatePVectorBalancedthickness( Vec pg);
    126127                void      CreatePVectorBalancedvelocities( Vec pg);
    127                 void      CreatePVectorDiagnosticHoriz( Vec pg);
    128                 void      CreatePVectorAdjointHoriz( Vec pg);
     128                void      CreatePVectorAdjointMacAyeal( Vec pg);
     129                void      CreatePVectorAdjointPattyn( Vec pg);
    129130                void      CreatePVectorDiagnosticHutter( Vec pg);
     131                void      CreatePVectorDiagnosticMacAyeal( Vec pg);
     132                void      CreatePVectorDiagnosticPattyn( Vec pg);
    130133                void      CreatePVectorDiagnosticStokes( Vec pg);
    131134                void      CreatePVectorAdjointStokes( Vec pg);

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

    r5158 r5166  
    658658        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    659659        if (analysis_type==DiagnosticHorizAnalysisEnum){
    660                 CreateKMatrixDiagnosticHoriz( Kgg);
     660                CreateKMatrixDiagnosticMacAyeal( Kgg);
    661661        }
    662662        else if (analysis_type==AdjointHorizAnalysisEnum){
    663663                /*Same as diagnostic*/
    664                 CreateKMatrixDiagnosticHoriz( Kgg);
     664                CreateKMatrixDiagnosticMacAyeal( Kgg);
    665665        }
    666666        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    710710        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    711711        if (analysis_type==DiagnosticHorizAnalysisEnum){
    712                 CreatePVectorDiagnosticHoriz( pg);
     712                CreatePVectorDiagnosticMacAyeal( pg);
    713713        }
    714714        else if (analysis_type==AdjointHorizAnalysisEnum){
    … …  
    27512751}
    27522752/*}}}*/
    2753 /*FUNCTION Tria::CreateKMatrixDiagnosticHoriz {{{1*/
    2754 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg){
     2753/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
     2754void  Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg){
    27552755
    27562756        /* local declarations */
    … …  
    40994099}
    41004100/*}}}*/
    4101 /*FUNCTION Tria::CreatePVectorDiagnosticHoriz {{{1*/
    4102 void Tria::CreatePVectorDiagnosticHoriz( Vec pg){
     4101/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
     4102void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg){
    41034103
    41044104        int             i,j;

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

    r5096 r5166  
    114114                void      CreateKMatrixBalancedthickness_CG(Mat Kgg);
    115115                void      CreateKMatrixBalancedvelocities(Mat Kgg);
    116                 void      CreateKMatrixDiagnosticHoriz(Mat Kgg);
     116                void      CreateKMatrixDiagnosticMacAyeal(Mat Kgg);
    117117                void      CreateKMatrixDiagnosticHorizFriction(Mat Kgg);
    118118                void      CreateKMatrixDiagnosticHutter(Mat Kgg);
    … …  
    127127                void      CreatePVectorBalancedvelocities(Vec pg);
    128128                void      CreatePVectorDiagnosticBaseVert(Vec pg);
    129                 void      CreatePVectorDiagnosticHoriz(Vec pg);
     129                void      CreatePVectorDiagnosticMacAyeal(Vec pg);
    130130                void      CreatePVectorAdjointHoriz(Vec pg);
    131131                void      CreatePVectorAdjointStokes(Vec pg);
Note: See TracChangeset for help on using the changeset viewer.