Changeset 92


Ignore:
Timestamp:
04/28/09 15:18:29 (16 years ago)
Author:
Eric.Larour
Message:

Pattyn element capability

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r46 r92  
    1414#include "../EnumDefinitions/EnumDefinitions.h"
    1515#include "../shared/shared.h"
     16#include "../include/typedefs.h"
    1617
    1718Penta::Penta(){
    … …  
    3132        for(i =0;i<6;i++){
    3233                node_ids[i] = penta_node_ids[i];
     34                node_offsets[i]=UNDEF;
     35                nodes[i]=NULL;
    3336                h[i] = penta_h[i];
    3437                s[i] = penta_s[i];
    … …  
    3942                geothermalflux[i] = penta_geothermalflux[i];
    4043        }
     44        matice=NULL;
     45        matice_offset=UNDEF;
     46        matpar=NULL;
     47        matpar_offset=UNDEF;
    4148
    4249        friction_type = penta_friction_type;
    … …  
    6572        printf("   mparid: %i\n",mparid);
    6673
    67         printf("   nodes=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]);
     74        printf("   node_ids=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]);
     75        printf("   node_offsets=[%i,%i,%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3],node_offsets[4],node_offsets[5]);
     76        printf("   matice_offset=%i\n",matice_offset);
     77        printf("   matpar_offset=%i\n",matpar_offset);
     78
    6879        printf("   h=[%i,%i,%i,%i,%i,%i]\n",h[0],h[1],h[2],h[3],h[4],h[5]);
    6980        printf("   s=[%i,%i,%i,%i,%i,%i]\n",s[0],s[1],s[2],s[3],s[4],s[5]);
    … …  
    108119        memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
    109120        memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
     121        memcpy(marshalled_dataset,&nodes,sizeof(nodes));marshalled_dataset+=sizeof(nodes);
     122        memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
     123        memcpy(marshalled_dataset,&matice,sizeof(matice));marshalled_dataset+=sizeof(matice);
     124        memcpy(marshalled_dataset,&matice_offset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset);
     125        memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
     126        memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
    110127        memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h);
    111128        memcpy(marshalled_dataset,&s,sizeof(s));marshalled_dataset+=sizeof(s);
    … …  
    133150int   Penta::MarshallSize(){
    134151
    135         return sizeof(id)+sizeof(mid)+sizeof(mparid)+sizeof(node_ids)+sizeof(h)+sizeof(s)+sizeof(b)+sizeof(k)+sizeof(friction_type)+sizeof(p)+sizeof(q)+sizeof(shelf)+sizeof(onbed)+sizeof(onsurface)+sizeof(meanvel)+sizeof(epsvel)+sizeof(collapse)+sizeof(melting)+sizeof(accumulation)+sizeof(geothermalflux)+sizeof(artdiff)+sizeof(thermal_steadystate) +sizeof(int); //sizeof(int) for enum type
     152        return sizeof(id)+
     153                sizeof(mid)+
     154                sizeof(mparid)+
     155                sizeof(node_ids)+
     156                sizeof(nodes)+
     157                sizeof(node_offsets)+
     158                sizeof(matice)+
     159                sizeof(matice_offset)+
     160                sizeof(matpar)+
     161                sizeof(matpar_offset)+
     162                sizeof(h)+
     163                sizeof(s)+
     164                sizeof(b)+
     165                sizeof(k)+
     166                sizeof(friction_type)+
     167                sizeof(p)+
     168                sizeof(q)+
     169                sizeof(shelf)+
     170                sizeof(onbed)+
     171                sizeof(onsurface)+
     172                sizeof(meanvel)+
     173                sizeof(epsvel)+
     174                sizeof(collapse)+
     175                sizeof(melting)+
     176                sizeof(accumulation)+
     177                sizeof(geothermalflux)+
     178                sizeof(artdiff)+
     179                sizeof(thermal_steadystate) +
     180                sizeof(int); //sizeof(int) for enum type
    136181}
    137182
    … …  
    142187void  Penta::Demarshall(char** pmarshalled_dataset){
    143188
     189        int i;
    144190        char* marshalled_dataset=NULL;
    145191
    … …  
    154200        memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
    155201        memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
     202        memcpy(&nodes,marshalled_dataset,sizeof(nodes));marshalled_dataset+=sizeof(nodes);
     203        memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
     204        memcpy(&matice,marshalled_dataset,sizeof(matice));marshalled_dataset+=sizeof(matice);
     205        memcpy(&matice_offset,marshalled_dataset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset);
     206        memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
     207        memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
    156208        memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
    157209        memcpy(&s,marshalled_dataset,sizeof(s));marshalled_dataset+=sizeof(s);
    … …  
    173225        memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
    174226
     227        /*nodes and materials are not pointing to correct objects anymore:*/
     228        for(i=0;i<6;i++)nodes[i]=NULL;
     229        matice=NULL;
     230        matpar=NULL;
     231
    175232        /*return: */
    176233        *pmarshalled_dataset=marshalled_dataset;
    … …  
    189246}
    190247
     248
    191249#undef __FUNCT__
    192250#define __FUNCT__ "Penta::Configure"
    193 void  Penta::Configure(void* loads,void* nodes,void* materials){
    194 
    195         throw ErrorException(__FUNCT__," not supported yet!");
    196 
    197 }
     251void  Penta::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin){
     252
     253        int i;
     254       
     255        DataSet* loadsin=NULL;
     256        DataSet* nodesin=NULL;
     257        DataSet* materialsin=NULL;
     258
     259        /*Recover pointers :*/
     260        loadsin=(DataSet*)ploadsin;
     261        nodesin=(DataSet*)pnodesin;
     262        materialsin=(DataSet*)pmaterialsin;
     263
     264        /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
     265        ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin);
     266       
     267        /*Same for materials: */
     268        ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
     269        ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
     270
     271}
     272
    198273#undef __FUNCT__
    199274#define __FUNCT__ "Penta::CreateKMatrix"
    … …  
    201276void  Penta::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
    202277
    203         throw ErrorException(__FUNCT__," not supported yet!");
    204 
     278        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     279        if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
     280
     281                CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);
     282
     283        }
     284        else{
     285                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
     286        }
     287
     288}
     289
     290
     291#undef __FUNCT__
     292#define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz"
     293void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type){
     294
     295
     296        /* local declarations */
     297        int             i,j;
     298
     299        /* node data: */
     300        const int    numgrids=6;
     301        const int    numdofs=2*numgrids;
     302        double       xyz_list[numgrids][3];
     303        int          doflist[numdofs];
     304        int          numberofdofspernode;
     305       
     306       
     307        /* 3d gaussian points: */
     308        int     num_gauss,ig;
     309        double* first_gauss_area_coord  =  NULL;
     310        double* second_gauss_area_coord =  NULL;
     311        double* third_gauss_area_coord  =  NULL;
     312        double* fourth_gauss_vert_coord  =  NULL;
     313        double* area_gauss_weights           =  NULL;
     314        double* vert_gauss_weights           =  NULL;
     315        int     ig1,ig2;
     316        double  gauss_weight1,gauss_weight2;
     317        double  gauss_l1l2l3l4[4];
     318        int     order_area_gauss;
     319        int     num_vert_gauss;
     320        int     num_area_gauss;
     321        double  gauss_weight;
     322       
     323        /* 2d gaussian point: */
     324        int     num_gauss2d;
     325        double* first_gauss_area_coord2d  =  NULL;
     326        double* second_gauss_area_coord2d =  NULL;
     327        double* third_gauss_area_coord2d  =  NULL;
     328        double* gauss_weights2d=NULL;
     329        double  gauss_l1l2l3[3];
     330
     331        /* material data: */
     332        double viscosity; //viscosity
     333
     334        /* strain rate: */
     335        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     336
     337        /* matrices: */
     338        double B[5][numdofs];
     339        double Bprime[5][numdofs];
     340        double L[2][numdofs];
     341        double D[5][5]={{ 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};              // material matrix, simple scalar matrix.
     342        double D_scalar;
     343        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
     344        double DL_scalar;
     345
     346        /* local element matrices: */
     347        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
     348        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
     349        double Ke_gg_drag_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
     350        double Jdet;
     351       
     352        /*slope: */
     353        double  slope[2]={0.0,0.0};
     354        double  slope_magnitude;
     355
     356        /*input parameters for structural analysis (diagnostic): */
     357        double* velocity_param=NULL;
     358        double  vxvy_list[numgrids][2];
     359        double  thickness;
     360
     361        /*friction: */
     362        double  alpha2_list[3];
     363        double  alpha2;
     364
     365        double MAXSLOPE=.06; // 6 %
     366        double MOUNTAINKEXPONENT=10;
     367
     368        /*Collapsed formulation: */
     369        Tria*  tria=NULL;
     370
     371        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
     372          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
     373          the stiffness matrix. */
     374
     375        if ((collapse==1) && (onbed==0)){
     376                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
     377                 * dofs have already been frozen! Do nothing: */
     378                return;
     379        }
     380        else if ((collapse==1) && (onbed==1)){
     381
     382                /*This element should be collapsed into a tria element at its base. Create this tria element,
     383                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
     384                tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     385                tria->CreateKMatrix(Kgg,inputs, analysis_type);
     386                delete tria;
     387                return;
     388        }
     389        else{
     390
     391                /*Implement standard penta element: */
     392
     393                /*recover extra inputs from users, at current convergence iteration: */
     394                velocity_param=ParameterInputsRecover(inputs,"velocity");
     395
     396                /* Get node coordinates and dof list: */
     397                GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     398                GetDofList(&doflist[0],&numberofdofspernode);
     399
     400                /* Set Ke_gg to 0: */
     401                for(i=0;i<numdofs;i++){
     402                        for(j=0;j<numdofs;j++){
     403                                Ke_gg[i][j]=0.0;
     404                        }
     405                }
     406
     407                /*Initialize vxvy_list: */
     408                for(i=0;i<numgrids;i++){
     409                        if(velocity_param){
     410                                vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
     411                                vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
     412                        }
     413                        else{
     414                                vxvy_list[i][0]=0;
     415                                vxvy_list[i][1]=0;
     416                        }
     417                }
     418                       
     419                #ifdef _DEBUGELEMENTS_
     420                if(my_rank==RANK && id==ELID){
     421                        printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK);
     422                        printf("   rho_ice: %g\n",matice->GetRhoIce());
     423                        printf("   rho_water:%g \n",matice->GetRhoWater());
     424                        printf("   gravity: %g\n",matpar->GetG());
     425                        printf("   Velocity: \n");
     426                        for (i=0;i<numgrids;i++){
     427                                printf("      grid %i  [%g,%g,%g]\n",i,vxvy_list[i][0],vxvy_list[i][1],vxvy_list[i][2]);
     428                        }
     429                        printf("   B [%g %g %g %g %g %g]\n",B_list[0],B_list[1],B_list[2],B_list[3],B_list[4],B_list[5]);
     430                        printf("   K [%g %g %g %g %g %g]\n",K_list[0],K_list[1],K_list[2],K_list[3],K_list[4],K_list[5]);
     431                        printf("   thickness [%g %g %g %g %g %g]\n",thickness_list[0],thickness_list[1],thickness_list[2],thickness_list[3],thickness_list[4],thickness_list[5]);
     432                        printf("   surface [%g %g %g %g %g %g]\n",surface_list[0],surface_list[1],surface_list[2],surface_list[3],surface_list[4],surface_list[5]);
     433                        printf("   bed [%g %g %g %g %g %g]\n",bed_list[0],bed_list[1],bed_list[2],bed_list[3],bed_list[4],bed_list[5]);
     434                        printf("   temperature_average [%g %g %g %g %g %g]\n",temperature_average_list[0],temperature_average_list[1],temperature_average_list[2],temperature_average_list[3],temperature_average_list[4],temperature_average_list[5]);
     435                }
     436                #endif
     437
     438                /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     439                  get tria gaussian points as well as segment gaussian points. For tria gaussian
     440                  points, order of integration is 2, because we need to integrate the product tB*D*B'
     441                  which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     442                  points, same deal, which yields 3 gaussian points.*/
     443
     444                order_area_gauss=5;
     445                num_vert_gauss=5;
     446
     447                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);
     448
     449                #ifdef _DEBUGGAUSS_
     450                if(my_rank==RANK && id==ELID){
     451                        printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK);
     452                        for (i=0;i<num_area_gauss;i++){
     453                                printf("   Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
     454                        }
     455                        for (i=0;i<num_vert_gauss;i++){
     456                                printf("   Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
     457                        }       
     458                }
     459                #endif
     460                /* Start  looping on the number of gaussian points: */
     461                for (ig1=0; ig1<num_area_gauss; ig1++){
     462                        for (ig2=0; ig2<num_vert_gauss; ig2++){
     463                       
     464                                /*Pick up the gaussian point: */
     465                                gauss_weight1=*(area_gauss_weights+ig1);
     466                                gauss_weight2=*(vert_gauss_weights+ig2);
     467                                gauss_weight=gauss_weight1*gauss_weight2;
     468                               
     469                               
     470                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
     471                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
     472                                gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
     473                                gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     474
     475
     476                                /*Get strain rate from velocity: */
     477                                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
     478                       
     479                                /*Get viscosity: */
     480                                matice->GetViscosity3d(&viscosity, &epsilon[0]);
     481
     482                                /*Get B and Bprime matrices: */
     483                                GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     484                                GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     485
     486                                /* Get Jacobian determinant: */
     487                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     488       
     489                                /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     490                                  onto this scalar matrix, so that we win some computational time: */
     491                                D_scalar=viscosity*gauss_weight*Jdet;
     492                                for (i=0;i<5;i++){
     493                                        D[i][i]=D_scalar;
     494                                }
     495               
     496                                /*  Do the triple product tB*D*Bprime: */
     497                                TripleMultiply( &B[0][0],5,numdofs,1,
     498                                                &D[0][0],5,5,0,
     499                                                &Bprime[0][0],5,numdofs,0,
     500                                                &Ke_gg_gaussian[0][0],0);
     501
     502                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     503                                for( i=0; i<numdofs; i++){
     504                                        for(j=0;j<numdofs;j++){
     505                                                Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     506                                        }
     507                                }
     508                        } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     509                } //for (ig1=0; ig1<num_area_gauss; ig1++)
     510               
     511
     512                /*Add Ke_gg to global matrix Kgg: */
     513                MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     514       
     515                //Deal with 2d friction at the bedrock interface
     516                if((onbed && !shelf)){
     517
     518                        /*Build a tria element using the 3 grids of the base of the penta. Then use
     519                         * the tria functionality to build a friction stiffness matrix on these 3
     520                         * grids: */
     521
     522                        tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     523                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
     524                        delete tria;
     525                }
     526
     527        }
     528               
     529        cleanup_and_return:
     530        xfree((void**)&first_gauss_area_coord);
     531        xfree((void**)&second_gauss_area_coord);
     532        xfree((void**)&third_gauss_area_coord);
     533        xfree((void**)&fourth_gauss_vert_coord);
     534        xfree((void**)&area_gauss_weights);
     535        xfree((void**)&vert_gauss_weights);
     536        xfree((void**)&first_gauss_area_coord2d);
     537        xfree((void**)&second_gauss_area_coord2d);
     538        xfree((void**)&third_gauss_area_coord2d);
     539        xfree((void**)&gauss_weights2d);
    205540}
    206541
    … …  
    208543#define __FUNCT__ "Penta::CreatePVector"
    209544void  Penta::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
    210        
    211         throw ErrorException(__FUNCT__," not supported yet!");
     545
     546        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     547        if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
     548
     549                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
     550
     551        }
     552        else{
     553                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
     554        }       
    212555
    213556}
    … …  
    215558#define __FUNCT__ "Penta::UpdateFromInputs"
    216559void  Penta::UpdateFromInputs(ParameterInputs* inputs){
    217        
    218         throw ErrorException(__FUNCT__," not supported yet!");
    219 
    220 }
    221 
     560
     561        int i;
     562        int doflist[MAXDOFSPERNODE*3];
     563        int numberofdofspernode;
     564
     565        double* thickness=NULL;
     566        double* bed=NULL;
     567        double* surface=NULL;
     568        double* drag=NULL;
     569        double* temperature_average_param=NULL;
     570        double* flow_law_param=NULL;
     571        double temperature_average;
     572        double B_average;
     573
     574
     575
     576        /*Get dof list: */
     577        GetDofList(&doflist[0],&numberofdofspernode);
     578
     579        /*Update internal data if inputs holds new values: */
     580        thickness=ParameterInputsRecover(inputs,"thickness");
     581        bed=ParameterInputsRecover(inputs,"bed");
     582        surface=ParameterInputsRecover(inputs,"surface");
     583        drag=ParameterInputsRecover(inputs,"drag");
     584
     585        for(i=0;i<5;i++){
     586                if(thickness)h[i]=thickness[doflist[i*numberofdofspernode+0]];
     587                if(bed)b[i]=bed[doflist[i*numberofdofspernode+0]];
     588                if(surface)s[i]=surface[doflist[i*numberofdofspernode+0]];
     589                if(drag)k[i]=drag[doflist[i*numberofdofspernode+0]];
     590        }
     591
     592        //Update material if necessary
     593        temperature_average_param=ParameterInputsRecover(inputs,"temperature_average");
     594        flow_law_param=ParameterInputsRecover(inputs,"B");
     595        if(temperature_average_param){
     596                for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
     597                temperature_average=temperature_average/3.0;
     598                B_average=Paterson(temperature_average);
     599                matice->SetB(B_average);
     600        }
     601
     602        //Update material if B is provided.
     603        if(flow_law_param){
     604                for(i=0;i<3;i++) B_average+=flow_law_param[doflist[i*numberofdofspernode+0]];
     605                B_average=B_average/3.0;
     606                matice->SetB(B_average);
     607        }
     608
     609}
     610               
    222611Matpar* Penta::GetMatPar(){
    223612        return matpar;
    … …  
    259648#define __FUNCT__ "Penta::Du"
    260649void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){
    261        
    262         throw ErrorException(__FUNCT__," not supported yet!");
    263 
    264 
     650
     651        int i;
     652        Tria* tria=NULL;
     653       
     654        /*Bail out if this element does not touch the surface: */
     655        if (onsurface){
     656                return;
     657        }
     658        else{
     659               
     660                tria=SpawnTria(3,4,5); //grids 0, 1 and 2 make the new tria.
     661                tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type);
     662                delete tria;
     663                return;
     664        }
    265665}
    266666
    267667#undef __FUNCT__
    268668#define __FUNCT__ "Penta::Gradj"
    269 void  Penta::Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type){
    270         throw ErrorException(__FUNCT__," not supported yet!");
    271 }
     669void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){
     670
     671        if (strcmp(control_type,"drag")==0){
     672                GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
     673        }
     674        else if (strcmp(control_type,"B")==0){
     675                GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
     676        }
     677        else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
     678}
     679
    272680#undef __FUNCT__
    273681#define __FUNCT__ "Penta::GradjDrag"
    274682void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type){
    275         throw ErrorException(__FUNCT__," not supported yet!");
    276 }
     683       
     684       
     685        Tria* tria=NULL;
     686       
     687        /*Bail out if this element does not touch the bedrock: */
     688        if (!onbed){
     689                return;
     690        }
     691        else{
     692               
     693                tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     694                tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type);
     695                delete tria;
     696                return;
     697        }
     698}
     699
    277700#undef __FUNCT__
    278701#define __FUNCT__ "Penta::GradjB"
    … …  
    283706#undef __FUNCT__
    284707#define __FUNCT__ "Penta::Misfit"
    285 double Penta::Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){
    286         throw ErrorException(__FUNCT__," not supported yet!");
    287 }
     708double Penta::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     709       
     710        double J;
     711        Tria* tria=NULL;
     712       
     713        /*Bail out if this element does not touch the surface: */
     714        if (!onsurface){
     715                return 0;
     716        }
     717        else{
     718
     719                /*use the Tria::CreateMisfit routine, on a Tria Element which is made of the 3 upper grids at the surface of the ice sheet. : */
     720                tria=SpawnTria(3,4,5);
     721                J=tria->Misfit( velocity,obs_velocity,inputs,analysis_type);
     722                delete tria;
     723                return J;
     724        }
     725}
     726               
     727#undef __FUNCT__
     728#define __FUNCT__ "Penta::SpawnTria"
     729Tria*  Penta::SpawnTria(int g0, int g1, int g2){
     730
     731        /*out of grids g0,g1 and g2 from Penta, build a tria element: */
     732        Tria* tria=NULL;
     733        double   tria_h[3];
     734        double   tria_s[3];
     735        double   tria_b[3];
     736        double   tria_c[3];
     737        double   tria_k[3];
     738       
     739        /*configuration: */
     740        int tria_node_ids[3];
     741        Node* tria_nodes[3];
     742        int tria_node_offsets[3];
     743
     744        tria_h[0]=h[g0];
     745        tria_h[1]=h[g1];
     746        tria_h[2]=h[g2];
     747
     748        tria_s[0]=s[g0];
     749        tria_s[1]=s[g1];
     750        tria_s[2]=s[g2];
     751
     752        tria_b[0]=b[g0];
     753        tria_b[1]=b[g1];
     754        tria_b[2]=b[g2];
     755
     756        tria_k[0]=k[g0];
     757        tria_k[1]=k[g1];
     758        tria_k[2]=k[g2];
     759
     760        tria_node_ids[0]=node_ids[g0];
     761        tria_node_ids[1]=node_ids[g1];
     762        tria_node_ids[2]=node_ids[g2];
     763
     764        tria_nodes[0]=nodes[g0];
     765        tria_nodes[1]=nodes[g1];
     766        tria_nodes[2]=nodes[g2];
     767
     768        tria_node_offsets[0]=node_offsets[g0];
     769        tria_node_offsets[1]=node_offsets[g1];
     770        tria_node_offsets[2]=node_offsets[g2];
     771
     772        tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, friction_type,p,q,shelf,meanvel,epsvel);
     773
     774        tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets);
     775        tria->MaticeConfiguration(matice,matice_offset);
     776        tria->MatparConfiguration(matpar,matpar_offset);
     777
     778        return tria;
     779
     780}
     781
     782
     783void  Penta::GetDofList(int* doflist,int* pnumberofdofspernode){
     784
     785        int i,j;
     786        int doflist_per_node[MAXDOFSPERNODE];
     787        int numberofdofspernode;
     788       
     789        for(i=0;i<6;i++){
     790                nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
     791                for(j=0;j<numberofdofspernode;j++){
     792                        doflist[i*numberofdofspernode+j]=doflist_per_node[j];
     793                }
     794        }
     795
     796        /*Assign output pointers:*/
     797        *pnumberofdofspernode=numberofdofspernode;
     798
     799}
     800
     801
     802#undef __FUNCT__
     803#define __FUNCT__ "Penta::GetStrainRate"
     804void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4){
     805
     806        int i;
     807        const int numgrids=6;
     808        const int NDOF2=2;
     809
     810        double B[5][NDOF2*numgrids];
     811
     812        /*Get B matrix: */
     813        GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
     814
     815        #ifdef _DEBUG_
     816        _printf_("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
     817        _printf_("              [ %lf   %lf  \n",B[1][0],B[1][1]);
     818        _printf_("              [ %lf   %lf  ]\n",B[2][0],B[2][1]);
     819        _printf_("              [ %lf   %lf  ]\n",B[3][0],B[3][1]);
     820        _printf_("              [ %lf   %lf  ]\n",B[4][0],B[4][1]);
     821       
     822        _printf_("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
     823        _printf_("              [ %lf   %lf  \n",B[1][2],B[1][3]);
     824        _printf_("              [ %lf   %lf  ]\n",B[2][2],B[2][3]);
     825        _printf_("              [ %lf   %lf  ]\n",B[3][2],B[3][3]);
     826        _printf_("              [ %lf   %lf  ]\n",B[4][2],B[4][3]);
     827       
     828        _printf_("B for grid3 : [ %lf   %lf  \n", B[0][4],B[0][5]);
     829        _printf_("              [ %lf   %lf  \n", B[1][4],B[1][5]);
     830        _printf_("              [ %lf   %lf  ]\n",B[2][4],B[2][5]);
     831        _printf_("              [ %lf   %lf  ]\n",B[3][4],B[3][5]);
     832        _printf_("              [ %lf   %lf  ]\n",B[4][4],B[4][5]);
     833       
     834        _printf_("B for grid4 : [ %lf   %lf  \n", B[0][6],B[0][7]);
     835        _printf_("              [ %lf   %lf  \n", B[1][6],B[1][7]);
     836        _printf_("              [ %lf   %lf  ]\n",B[2][6],B[2][7]);
     837        _printf_("              [ %lf   %lf  ]\n",B[3][6],B[3][7]);
     838        _printf_("              [ %lf   %lf  ]\n",B[4][6],B[4][7]);
     839                               
     840        _printf_("B for grid5 : [ %lf   %lf  \n", B[0][8],B[0][9]);
     841        _printf_("              [ %lf   %lf  \n", B[1][8],B[1][9]);
     842        _printf_("              [ %lf   %lf  ]\n",B[2][8],B[2][9]);
     843        _printf_("              [ %lf   %lf  ]\n",B[3][8],B[3][9]);
     844        _printf_("              [ %lf   %lf  ]\n",B[4][8],B[4][9]);
     845
     846        _printf_("B for grid6 : [ %lf   %lf  \n", B[0][10],B[0][11]);
     847        _printf_("              [ %lf   %lf  \n", B[1][10],B[1][11]);
     848        _printf_("              [ %lf   %lf  ]\n",B[2][10],B[2][11]);
     849        _printf_("              [ %lf   %lf  ]\n",B[3][10],B[3][11]);
     850        _printf_("              [ %lf   %lf  ]\n",B[4][10],B[4][11]);
     851
     852        for (i=0;i<numgrids;i++){
     853                _printf_("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
     854        }
     855        #endif
     856
     857        /*Multiply B by velocity, to get strain rate: */
     858        MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
     859                                      velocity,NDOF2*numgrids,1,0,
     860                                                  epsilon,0);
     861
     862}
     863
     864#undef __FUNCT__
     865#define __FUNCT__ "Penta::GetB"
     866void Penta::GetB(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     867
     868        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     869         * For grid i, Bi can be expressed in the basic coordinate system
     870         * by:
     871         *       Bi_basic=[ dh/dx          0      ]
     872         *                [   0           dh/dy   ]
     873         *                [ 1/2*dh/dy  1/2*dh/dx  ]
     874         *                [ 1/2*dh/dz      0      ]
     875         *                [  0         1/2*dh/dz  ]
     876         * where h is the interpolation function for grid i.
     877         *
     878         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
     879         */
     880       
     881        int i;
     882        const int numgrids=6;
     883        const int NDOF3=3;
     884        const int NDOF2=2;
     885
     886        double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
     887
     888        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
     889        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     890
     891        #ifdef _DEBUG_
     892        for (i=0;i<numgrids;i++){
     893                _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
     894        }
     895        #endif
     896
     897        /*Build B: */
     898        for (i=0;i<numgrids;i++){
     899                *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i];
     900                *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
     901
     902                *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
     903                *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
     904
     905                *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[1][i];
     906                *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[0][i];
     907
     908                *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
     909                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
     910               
     911                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
     912                *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
     913        }
     914
     915}
     916
     917
     918#undef __FUNCT__
     919#define __FUNCT__ "Penta::GetBPrime"
     920void Penta::GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     921
     922        /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
     923         * For grid i, Bi can be expressed in the basic coordinate system
     924         * by:
     925         *       Bi_basic=[ 2*dh/dx     dh/dy   ]
     926         *                [   dh/dx    2*dh/dy  ]
     927         *                [ dh/dy      dh/dx    ]
     928         *                [ dh/dz         0     ]
     929         *                [  0         dh/dz    ]
     930         * where h is the interpolation function for grid i.
     931         *
     932         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
     933         */
     934       
     935        int i;
     936        const int NDOF3=3;
     937        const int NDOF2=2;
     938        const int numgrids=6;
     939
     940        double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
     941
     942        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
     943        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     944
     945        #ifdef _DEBUG_
     946        for (i=0;i<numgrids;i++){
     947                _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
     948        }
     949        #endif
     950
     951        /*Build BPrime: */
     952        for (i=0;i<numgrids;i++){
     953                *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh2dh3dh4dh5dh6_basic[0][i];
     954                *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[1][i];
     955
     956                *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[0][i];
     957                *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh2dh3dh4dh5dh6_basic[1][i];
     958
     959                *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[1][i];
     960                *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[0][i];
     961
     962                *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i];
     963                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
     964               
     965                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
     966                *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i];
     967        }
     968}
     969
     970#undef __FUNCT__
     971#define __FUNCT__ "Penta::GetJacobianDeterminant"
     972void Penta::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4){
     973
     974        /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
     975         * the determinant of it: */
     976        const int NDOF3=3;
     977       
     978        double J[NDOF3][NDOF3];
     979
     980        GetJacobian(&J[0][0],xyz_list,gauss_l1l2l3l4);
     981
     982        *Jdet= J[0][0]*J[1][1]*J[2][2]-J[0][0]*J[1][2]*J[2][1]-J[1][0]*J[0][1]*J[2][2]+J[1][0]*J[0][2]*J[2][1]+J[2][0]*J[0][1]*J[1][2]-J[2][0]*J[0][2]*J[1][1];
     983        if(*Jdet<0){
     984                //printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
     985        }
     986}
     987
     988#undef __FUNCT__
     989#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic"
     990void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
     991       
     992        /*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
     993
     994       
     995        int i;
     996        const int NDOF3=3;
     997        const int numgrids=6;
     998
     999        double dh1dh2dh3dh4dh5dh6_param[NDOF3][numgrids];
     1000        double Jinv[NDOF3][NDOF3];
     1001
     1002        /*Get derivative values with respect to parametric coordinate system: */
     1003        GetNodalFunctionsDerivativesParams(&dh1dh2dh3dh4dh5dh6_param[0][0], gauss_l1l2l3l4);
     1004
     1005        /*Get Jacobian invert: */
     1006        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3l4);
     1007
     1008        /*Build dh1dh2dh3_basic:
     1009         *
     1010         * [dhi/dx]= Jinv*[dhi/dr]
     1011         * [dhi/dy]       [dhi/ds]
     1012         * [dhi/dz]       [dhi/dn]
     1013         */
     1014
     1015        for (i=0;i<numgrids;i++){
     1016                *(dh1dh2dh3dh4dh5dh6_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[0][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[0][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
     1017                *(dh1dh2dh3dh4dh5dh6_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[1][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[1][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
     1018                *(dh1dh2dh3dh4dh5dh6_basic+numgrids*2+i)=Jinv[2][0]*dh1dh2dh3dh4dh5dh6_param[0][i]+Jinv[2][1]*dh1dh2dh3dh4dh5dh6_param[1][i]+Jinv[2][2]*dh1dh2dh3dh4dh5dh6_param[2][i];
     1019        }
     1020
     1021}
     1022
     1023#undef __FUNCT__
     1024#define __FUNCT__ "Penta::GetJacobian"
     1025void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4){
     1026
     1027        const int NDOF3=3;
     1028        int i,j;
     1029       
     1030        /*The Jacobian is constant over the element, discard the gaussian points.
     1031         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1032
     1033        double A1,A2,A3; //area coordinates
     1034        double xi,eta,zi; //parametric coordinates
     1035
     1036        double x1,x2,x3,x4,x5,x6;
     1037        double y1,y2,y3,y4,y5,y6;
     1038        double z1,z2,z3,z4,z5,z6;
     1039       
     1040        double sqrt3=sqrt(3.0);
     1041       
     1042        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
     1043        A1=gauss_l1l2l3l4[0];
     1044        A2=gauss_l1l2l3l4[1];
     1045        A3=gauss_l1l2l3l4[2];
     1046
     1047        xi=A2-A1;
     1048        eta=sqrt3*A3;
     1049        zi=gauss_l1l2l3l4[3];
     1050
     1051        x1=*(xyz_list+3*0+0);
     1052        x2=*(xyz_list+3*1+0);
     1053        x3=*(xyz_list+3*2+0);
     1054        x4=*(xyz_list+3*3+0);
     1055        x5=*(xyz_list+3*4+0);
     1056        x6=*(xyz_list+3*5+0);
     1057       
     1058        y1=*(xyz_list+3*0+1);
     1059        y2=*(xyz_list+3*1+1);
     1060        y3=*(xyz_list+3*2+1);
     1061        y4=*(xyz_list+3*3+1);
     1062        y5=*(xyz_list+3*4+1);
     1063        y6=*(xyz_list+3*5+1);
     1064
     1065        z1=*(xyz_list+3*0+2);
     1066        z2=*(xyz_list+3*1+2);
     1067        z3=*(xyz_list+3*2+2);
     1068        z4=*(xyz_list+3*3+2);
     1069        z5=*(xyz_list+3*4+2);
     1070        z6=*(xyz_list+3*5+2);
     1071
     1072
     1073        *(J+NDOF3*0+0)=1.0/4.0*(x1-x2-x4+x5)*zi+1.0/4.0*(-x1+x2-x4+x5);
     1074        *(J+NDOF3*1+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+sqrt3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
     1075        *(J+NDOF3*2+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +1.0/4.0*(-x1+x5-x2+x4);
     1076
     1077        *(J+NDOF3*0+1)=1.0/4.0*(y1-y2-y4+y5)*zi+1.0/4.0*(-y1+y2-y4+y5);
     1078        *(J+NDOF3*1+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+sqrt3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
     1079        *(J+NDOF3*2+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+1.0/4.0*(y1-y2-y4+y5)*xi+1.0/4.0*(y4-y1+y5-y2);
     1080
     1081        *(J+NDOF3*0+2)=1.0/4.0*(z1-z2-z4+z5)*zi+1.0/4.0*(-z1+z2-z4+z5);
     1082        *(J+NDOF3*1+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+sqrt3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
     1083        *(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
     1084
     1085        #ifdef _DEBUG_
     1086        for(i=0;i<3;i++){
     1087                for (j=0;j<3;j++){
     1088                        printf("%lf ",*(J+NDOF3*i+j));
     1089                }
     1090                printf("\n");
     1091        }
     1092        #endif
     1093}
     1094
     1095#undef __FUNCT__
     1096#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams"
     1097void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
     1098       
     1099        /*This routine returns the values of the nodal functions derivatives  (with respect to the
     1100         * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
     1101
     1102        const int numgrids=6;
     1103        double A1,A2,A3,z;
     1104        double sqrt3=sqrt(3.0);
     1105       
     1106        A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
     1107        A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
     1108        A3=gauss_l1l2l3l4[2]; //third area coordinate value  In term of xi and eta: A3=y/sqrt(3);
     1109        z=gauss_l1l2l3l4[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
     1110
     1111
     1112        /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
     1113        *(dl1dl2dl3dl4dl5dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
     1114        *(dl1dl2dl3dl4dl5dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
     1115        *(dl1dl2dl3dl4dl5dl6+numgrids*2+0)=-1.0/2.0*A1;
     1116
     1117        /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
     1118        *(dl1dl2dl3dl4dl5dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
     1119        *(dl1dl2dl3dl4dl5dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
     1120        *(dl1dl2dl3dl4dl5dl6+numgrids*2+1)=-1.0/2.0*A2;
     1121
     1122        /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
     1123        *(dl1dl2dl3dl4dl5dl6+numgrids*0+2)=0.0;
     1124        *(dl1dl2dl3dl4dl5dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
     1125        *(dl1dl2dl3dl4dl5dl6+numgrids*2+2)=-1.0/2.0*A3;
     1126
     1127        /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
     1128        *(dl1dl2dl3dl4dl5dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
     1129        *(dl1dl2dl3dl4dl5dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
     1130        *(dl1dl2dl3dl4dl5dl6+numgrids*2+3)=1.0/2.0*A1;
     1131
     1132        /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
     1133        *(dl1dl2dl3dl4dl5dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
     1134        *(dl1dl2dl3dl4dl5dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
     1135        *(dl1dl2dl3dl4dl5dl6+numgrids*2+4)=1.0/2.0*A2;
     1136
     1137        /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
     1138        *(dl1dl2dl3dl4dl5dl6+numgrids*0+5)=0.0;
     1139        *(dl1dl2dl3dl4dl5dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
     1140        *(dl1dl2dl3dl4dl5dl6+numgrids*2+5)=1.0/2.0*A3;
     1141}
     1142
     1143#undef __FUNCT__
     1144#define __FUNCT__ "Penta::GetJacobianInvert"
     1145void Penta::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4){
     1146
     1147        double Jdet;
     1148        const int NDOF3=3;
     1149
     1150        /*Call Jacobian routine to get the jacobian:*/
     1151        GetJacobian(Jinv, xyz_list, gauss_l1l2l3l4);
     1152
     1153        /*Invert Jacobian matrix: */
     1154        MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
     1155}
     1156
     1157
     1158
     1159#undef __FUNCT__
     1160#define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz"
     1161void Penta::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type){
     1162
     1163        int i,j;
     1164
     1165        /* node data: */
     1166        const int    numgrids=6;
     1167        const int    NDOF2=2;
     1168        const int    numdofs=NDOF2*numgrids;
     1169        double       xyz_list[numgrids][3];
     1170        int          doflist[numdofs];
     1171        int          numberofdofspernode;
     1172       
     1173        /* parameters: */
     1174        double  slope[NDOF2];
     1175        double  driving_stress_baseline;
     1176        double  thickness;
     1177
     1178        /* gaussian points: */
     1179        int     num_gauss,ig;
     1180        double* first_gauss_area_coord  =  NULL;
     1181        double* second_gauss_area_coord =  NULL;
     1182        double* third_gauss_area_coord  =  NULL;
     1183        double* fourth_gauss_vert_coord  =  NULL;
     1184        double* area_gauss_weights      =  NULL;
     1185        double* vert_gauss_weights      =  NULL;
     1186        double  gauss_l1l2l3l4[4];
     1187        int     order_area_gauss;
     1188        int     num_vert_gauss;
     1189        int     num_area_gauss;
     1190        int     ig1,ig2;
     1191        double  gauss_weight1,gauss_weight2;
     1192        double  gauss_weight;
     1193
     1194        /* Jacobian: */
     1195        double Jdet;
     1196
     1197        /*nodal functions: */
     1198        double l1l2l3l4l5l6[6];
     1199
     1200        /*element vector at the gaussian points: */
     1201        double  pe_g[numdofs];
     1202        double  pe_g_gaussian[numdofs];
     1203
     1204        /*Spawning: */
     1205        Tria* tria=NULL;
     1206
     1207
     1208        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
     1209          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
     1210          the load vector. */
     1211
     1212        if ((collapse==1) && (onbed==0)){
     1213                /*This element should be collapsed, but this element is not on the bedrock, therefore all its
     1214                 * dofs have already been frozen! Do nothing: */
     1215                return;
     1216        }
     1217        else if ((collapse==1) && (onbed==1)){
     1218
     1219                /*This element should be collapsed into a tria element at its base. Create this tria element,
     1220                 *and use its CreatePVector functionality to return an elementary load vector: */
     1221                tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1222                tria->CreatePVector(pg,inputs, analysis_type);
     1223                delete tria;
     1224                return;
     1225        }
     1226        else{
     1227
     1228                /*Implement standard penta element: */
     1229
     1230                /* Get node coordinates and dof list: */
     1231                GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1232                GetDofList(&doflist[0],&numberofdofspernode);
     1233
     1234                /* Set pe_g to 0: */
     1235                for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     1236
     1237                /*Get gaussian points and weights :*/
     1238                order_area_gauss=2;
     1239                num_vert_gauss=3;
     1240
     1241                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);
     1242                #ifdef _DEBUG_
     1243                for (i=0;i<num_area_gauss;i++){
     1244                        _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
     1245                }
     1246                for (i=0;i<num_vert_gauss;i++){
     1247                        _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
     1248                }
     1249                #endif
     1250       
     1251                /* Start  looping on the number of gaussian points: */
     1252                for (ig1=0; ig1<num_area_gauss; ig1++){
     1253                        for (ig2=0; ig2<num_vert_gauss; ig2++){
     1254                       
     1255                                /*Pick up the gaussian point: */
     1256                                gauss_weight1=*(area_gauss_weights+ig1);
     1257                                gauss_weight2=*(vert_gauss_weights+ig2);
     1258                                gauss_weight=gauss_weight1*gauss_weight2;
     1259                               
     1260                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
     1261                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
     1262                                gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
     1263                                gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     1264               
     1265                                /*Compute thickness at gaussian point: */
     1266                                GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
     1267
     1268                                /*Compute slope at gaussian point: */
     1269                                GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3l4);
     1270
     1271                                /* Get Jacobian determinant: */
     1272                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     1273               
     1274                                 /*Get nodal functions: */
     1275                                GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
     1276
     1277                                /*Compute driving stress: */
     1278                                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
     1279
     1280                                /*Build pe_g_gaussian vector: */
     1281                                for (i=0;i<numgrids;i++){
     1282                                        for (j=0;j<NDOF2;j++){
     1283                                                pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3l4l5l6[i];
     1284                                        }
     1285                                }
     1286
     1287                                /*Add pe_g_gaussian vector to pe_g: */
     1288                                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     1289
     1290                        } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     1291                } //for (ig1=0; ig1<num_area_gauss; ig1++)
     1292
     1293        } //else if ((collapse==1) && (onbed==1))
     1294
     1295        /*Add pe_g to global vector pg: */
     1296        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     1297
     1298        cleanup_and_return:
     1299        xfree((void**)&first_gauss_area_coord);
     1300        xfree((void**)&second_gauss_area_coord);
     1301        xfree((void**)&third_gauss_area_coord);
     1302        xfree((void**)&fourth_gauss_vert_coord);
     1303        xfree((void**)&area_gauss_weights);
     1304        xfree((void**)&vert_gauss_weights);
     1305
     1306}
     1307
     1308#undef __FUNCT__
     1309#define __FUNCT__ "Penta::CreateKMatrix"
     1310void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
     1311       
     1312        const int numgrids=6;
     1313        double l1l2l3l4l5l6[numgrids];
     1314
     1315        GetNodalFunctions(&l1l2l3l4l5l6[0], gauss_l1l2l3l4);
     1316
     1317        *pvalue=l1l2l3l4l5l6[0]*v_list[0]+l1l2l3l4l5l6[1]*v_list[1]+l1l2l3l4l5l6[2]*v_list[2]+l1l2l3l4l5l6[3]*v_list[3]+l1l2l3l4l5l6[4]*v_list[4]+l1l2l3l4l5l6[5]*v_list[5];
     1318}
     1319
     1320#undef __FUNCT__
     1321#define __FUNCT__ "Penta::GetParameterDerivativeValue"
     1322void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
     1323                               
     1324        /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4:
     1325         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
     1326         *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
     1327         *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
     1328         *
     1329         *   p is a vector of size 3x1 already allocated.
     1330         */
     1331
     1332        const int NDOF3=3;
     1333        const int numgrids=6;
     1334        double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
     1335
     1336        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
     1337        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     1338       
     1339        *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
     1340;
     1341        *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
     1342
     1343        *(p+2)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[2][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[2][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[2][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[2][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[2][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[2][5];
     1344
     1345}
     1346
     1347#undef __FUNCT__
     1348#define __FUNCT__ "Penta::GetNodalFunctions"
     1349void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
     1350       
     1351        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     1352
     1353        l1l2l3l4l5l6[0]=gauss_l1l2l3l4[0]*(1-gauss_l1l2l3l4[3])/2.0;
     1354
     1355        l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
     1356       
     1357        l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
     1358
     1359        l1l2l3l4l5l6[3]=gauss_l1l2l3l4[0]*(1+gauss_l1l2l3l4[3])/2.0;
     1360
     1361        l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
     1362       
     1363        l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
     1364
     1365}
     1366
     1367#undef __FUNCT__
     1368#define __FUNCT__ "Penta::VelocityExtrude"
     1369void  Penta::VelocityExtrude(Vec ug,double* ug_serial){
     1370
     1371        /* node data: */
     1372        const int    numgrids=6;
     1373        const int    numdofs=2*numgrids;
     1374        int          doflist[numdofs];
     1375        int          nodedofs[2];
     1376        int          numberofdofspernode;
     1377       
     1378        Node* node=NULL;
     1379        int   i;
     1380        double velocity[2];
     1381               
     1382
     1383        if((collapse==1) && (onbed==1)){
     1384                       
     1385                GetDofList(&doflist[0],&numberofdofspernode);
     1386
     1387                /*this penta is a collapsed macayeal. For each node on the base of this penta,
     1388                 * we grab the velocity. Once we know the velocity, we follow the upper nodes,
     1389                 * inserting the same velocity value into ug, until we reach the surface: */
     1390                for(i=0;i<3;i++){
     1391       
     1392                        node=nodes[i]; //base nodes
     1393       
     1394                        /*get velocity for this base node: */
     1395                        velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
     1396                        velocity[1]=ug_serial[doflist[numberofdofspernode*i+1]];
     1397
     1398                        //go througn all nodes which sit on top of this node, until we reach the surface,
     1399                        //and plug  velocity in ug
     1400                        for(;;){
     1401
     1402                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
     1403                                VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES);
     1404                                VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES);
     1405
     1406                                if (node->IsOnSurface())break;
     1407                                /*get next node: */
     1408                                node=node->GetUpperNode();
     1409                        }
     1410                }
     1411
     1412        }
     1413
     1414}
Note: See TracChangeset for help on using the changeset viewer.