Changeset 128


Ignore:
Timestamp:
04/29/09 14:23:19 (16 years ago)
Author:
Eric.Larour
Message:

Vertical velocity solution

Location:
issm/trunk
Files:
5 added
1 deleted
20 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp

    r102 r128  
    2424                return DiagnosticVertAnalysisEnum();
    2525        }
    26         else if (strcmp(analysis_type,"diagnostic_basevert")==0){
    27                 return DiagnosticBaseVertAnalysisEnum();
    28         }
    2926        else if (strcmp(analysis_type,"prognostic")==0){
    3027                return PrognosticAnalysisEnum();

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

    r8 r128  
    2020int DiagnosticHorizAnalysisEnum(void){  return          23; }
    2121int DiagnosticVertAnalysisEnum(void){   return          24; }
    22 int DiagnosticBaseVertAnalysisEnum(void){ return        25; }
    2322int PrognosticAnalysisEnum(void){       return          26; }
    2423int ThermalSteadyAnalysisEnum(void)     {   return          27; }

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r1 r128  
    3333int DiagnosticHorizAnalysisEnum(void);
    3434int DiagnosticVertAnalysisEnum(void);
    35 int DiagnosticBaseVertAnalysisEnum(void);
    3635int PrognosticAnalysisEnum(void);
    3736int ThermalSteadyAnalysisEnum(void);

  • issm/trunk/src/c/Makefile.am

    r118 r128  
    163163                                        ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
    164164                                        ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
     165                                        ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
     166                                        ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
     167                                        ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
    165168                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    166169                                        ./Dofx/Dofx.h\
    … …  
    378381                                        ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
    379382                                        ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
     383                                        ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
     384                                        ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
     385                                        ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
    380386                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    381387                                        ./Dofx/Dofx.h\

  • issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp

    r117 r128  
    2121                CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle);
    2222        }
    23         /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
    24                 CreateConstraintsDiagnosticBaseVert(pconstraints,model,model_handle);
    25         }
    2623        else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
    2724                CreateConstraintsDiagnosticVert(pconstraints,model,model_handle);
    28         }
     25        }/*
    2926        else if (strcmp(model->analysis_type,"melting")==0){
    3027                CreateConstraintsMelting(pconstraints,model,model_handle);

  • issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp

    r117 r128  
    2424       
    2525        }
    26         /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
     26        else if ((strcmp(model->analysis_type,"diagnostic_vert")==0)){
    2727
    28                 CreateElementsNodesAndMaterialsDataSetsDiagnosticBaseVert(pelements,pnodes,pmaterials, model,model->analysis_type);
    29 
     28                CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle);
     29       
    3030        }
    31         else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
    32 
    33                 CreateElementsNodesAndMaterialsDataSetsDiagnosticVert(pelements,pnodes,pmaterials, model,model->analysis_type);
    34         }
     31        /*
    3532        else if (strcmp(model->analysis_type,"melting")==0){
    3633
    … …  
    4643        }*/
    4744        else{
    48                 throw ErrorException(__FUNCT__," analysis_type not supported yet!");
     45                throw ErrorException(__FUNCT__,exprintf("%s%s%s"," analysis_type: ",model->analysis_type," not supported yet!"));
    4946        }
    5047}

  • issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp

    r117 r128  
    2222                CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
    2323        }
    24         /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){
    25 
    26                 CreateLoadsDiagnosticBaseVert(ploads,model,model_handle);
    27         }
    2824        else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
    2925
    3026                CreateLoadsDiagnosticVert(ploads,model,model_handle);
    31         }
     27        }/*
    3228        else if (strcmp(model->analysis_type,"melting")==0){
    3329

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp

    r117 r128  
    5959        double tria_b[3];
    6060        double tria_k[3];
     61        double tria_melting[3];
     62        double tria_accumulation[3];
    6163        int    tria_friction_type;
    6264        double tria_p;
    … …  
    278280
    279281                        /*Create tria element using its constructor:*/
    280                         tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
     282                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
    281283
    282284                        /*Add tria element to elements dataset: */

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r117 r128  
    170170        /*Create of fem datasets: specialised drivers: */
    171171       
    172         /*diagnostic:*/
     172        /*diagnostic horizontal*/
    173173        void    CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
    174174        void    CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    175175        void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle);
     176
     177        /*diagnostic vertical*/
     178        void    CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
     179        void    CreateConstraintsDiagnosticVert(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
     180        void    CreateLoadsDiagnosticVert(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    176181       
    177182        /*control:*/

  • issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp

    r1 r128  
    3232        loads->Configure(elements, loads, nodes, materials);
    3333
    34 
    3534        /*Get size of matrix: */
    3635        gsize=nodes->NumberOfDofs();
    … …  
    3938        if(kflag)Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode);
    4039        if(pflag)pg=NewVec(gsize);
    41        
     40
    4241        /*Fill stiffness matrix and right hand side vector, from elements: */
    4342        if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type);
    … …  
    5958        }
    6059       
    61 
    6260        /*Assign output pointers: */
    6361        *pKgg=Kgg;

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

    r111 r128  
    287287
    288288        }
     289        else if ((analysis_type==DiagnosticVertAnalysisEnum())){
     290
     291                CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type);
     292
     293        }
    289294        else{
    290295                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    … …  
    546551
    547552#undef __FUNCT__
     553#define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
     554void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type){
     555
     556        /* local declarations */
     557        int             i,j;
     558
     559        /* node data: */
     560        const int    numgrids=6;
     561        const int    NDOF1=1;
     562        const int    numdofs=NDOF1*numgrids;
     563        double       xyz_list[numgrids][3];
     564        int          doflist[numdofs];
     565        int          numberofdofspernode;
     566
     567        /* 3d gaussian points: */
     568        int     num_gauss,ig;
     569        double* first_gauss_area_coord  =  NULL;
     570        double* second_gauss_area_coord =  NULL;
     571        double* third_gauss_area_coord  =  NULL;
     572        double* fourth_gauss_vert_coord  =  NULL;
     573        double* area_gauss_weights           =  NULL;
     574        double* vert_gauss_weights           =  NULL;
     575        int     ig1,ig2;
     576        double  gauss_weight1,gauss_weight2;
     577        double  gauss_l1l2l3l4[4];
     578        int     order_area_gauss;
     579        int     num_vert_gauss;
     580        int     num_area_gauss;
     581        double  gauss_weight;
     582
     583        /* matrices: */
     584        double  Ke_gg[numdofs][numdofs];
     585        double  Ke_gg_gaussian[numdofs][numdofs];
     586        double  Jdet;
     587        double  B[NDOF1][numgrids];
     588        double  Bprime[NDOF1][numgrids];
     589        double  DL_scalar;
     590       
     591        /*Collapsed formulation: */
     592        Tria*  tria=NULL;
     593
     594        /*If this element is on the bed rock, spawn a Tria element out of the bed triangle, and use it
     595         * to create the stifness matrix for the base velocity: */
     596        if(onbed){
     597                tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 for the base
     598                tria->CreateKMatrixDiagnosticBaseVert(Kgg,inputs, analysis_type);
     599                delete tria;
     600        }
     601
     602        /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness
     603         * matrix: */
     604        if(onsurface){
     605                tria=SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
     606                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type);
     607                delete tria;
     608        }
     609
     610        /*Now, onto the formulation for the vertical velocity: */
     611
     612        /* Get node coordinates and dof list: */
     613        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     614        GetDofList(&doflist[0],&numberofdofspernode);
     615
     616        /* Set Ke_gg to 0: */
     617        for(i=0;i<numdofs;i++){
     618                for(j=0;j<numdofs;j++){
     619                        Ke_gg[i][j]=0.0;
     620                }
     621        }
     622
     623        /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
     624          get tria gaussian points as well as segment gaussian points. For tria gaussian
     625          points, order of integration is 2, because we need to integrate the product tB*D*B'
     626          which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     627          points, same deal, which yields 3 gaussian points.*/
     628
     629        order_area_gauss=2;
     630        num_vert_gauss=2;
     631
     632        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);
     633        #ifdef _DEBUG_
     634        for (i=0;i<num_area_gauss;i++){
     635                _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));
     636        }
     637        for (i=0;i<num_vert_gauss;i++){
     638                _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
     639        }
     640        #endif
     641       
     642        /* Start  looping on the number of gaussian points: */
     643        for (ig1=0; ig1<num_area_gauss; ig1++){
     644                for (ig2=0; ig2<num_vert_gauss; ig2++){
     645               
     646                        /*Pick up the gaussian point: */
     647                        gauss_weight1=*(area_gauss_weights+ig1);
     648                        gauss_weight2=*(vert_gauss_weights+ig2);
     649                        gauss_weight=gauss_weight1*gauss_weight2;
     650                       
     651                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
     652                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
     653                        gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
     654                        gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     655
     656                        /*Get B and Bprime matrices: */
     657                        GetB_vert(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     658                        GetBPrime_vert(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4);
     659
     660                        /* Get Jacobian determinant: */
     661                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     662
     663                        DL_scalar=gauss_weight*Jdet;
     664
     665                        /*  Do the triple product tB*D*Bprime: */
     666                        TripleMultiply( &B[0][0],1,numgrids,1,
     667                                        &DL_scalar,1,1,0,
     668                                        &Bprime[0][0],1,numgrids,0,
     669                                        &Ke_gg_gaussian[0][0],0);
     670
     671                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     672                        for( i=0; i<numdofs; i++){
     673                                for(j=0;j<numdofs;j++){
     674                                        Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     675                                }
     676                        }       
     677                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     678        } //for (ig1=0; ig1<num_area_gauss; ig1++)
     679
     680        /*Add Ke_gg to global matrix Kgg: */
     681        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     682       
     683        cleanup_and_return:
     684        xfree((void**)&first_gauss_area_coord);
     685        xfree((void**)&second_gauss_area_coord);
     686        xfree((void**)&third_gauss_area_coord);
     687        xfree((void**)&fourth_gauss_vert_coord);
     688        xfree((void**)&area_gauss_weights);
     689        xfree((void**)&vert_gauss_weights);
     690}
     691
     692#undef __FUNCT__
    548693#define __FUNCT__ "Penta::CreatePVector"
    549694void  Penta::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
    … …  
    553698
    554699                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
    555 
     700        }
     701        else if ((analysis_type==DiagnosticVertAnalysisEnum())){
     702
     703                CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
    556704        }
    557705        else{
    … …  
    741889        double   tria_c[3];
    742890        double   tria_k[3];
     891        double   tria_melting[3];
     892        double   tria_accumulation[3];
    743893       
    744894        /*configuration: */
    … …  
    763913        tria_k[2]=k[g2];
    764914
     915        tria_melting[0]=melting[g0];
     916        tria_melting[1]=melting[g1];
     917        tria_melting[2]=melting[g2];
     918       
     919        tria_accumulation[0]=accumulation[g0];
     920        tria_accumulation[1]=accumulation[g1];
     921        tria_accumulation[2]=accumulation[g2];
     922
    765923        tria_node_ids[0]=node_ids[g0];
    766924        tria_node_ids[1]=node_ids[g1];
    … …  
    775933        tria_node_offsets[2]=node_offsets[g2];
    776934
    777         tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
     935        tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
    778936
    779937        tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets);
    … …  
    14181576
    14191577}
     1578
     1579#undef __FUNCT__
     1580#define __FUNCT__ "Penta:GetB_vert"
     1581void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     1582
     1583
     1584        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
     1585        where hi is the interpolation function for grid i.*/
     1586
     1587        int i;
     1588        const int NDOF3=3;
     1589        const int numgrids=6;
     1590        double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids];
     1591
     1592        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
     1593        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
     1594
     1595        #ifdef _DEBUG_
     1596        for (i=0;i<numgrids;i++){
     1597                _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]);
     1598        }
     1599        #endif
     1600
     1601        /*Build B: */
     1602        for (i=0;i<numgrids;i++){
     1603                B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
     1604        }
     1605       
     1606}
     1607
     1608
     1609#undef __FUNCT__
     1610#define __FUNCT__ "Penta:GetBPrime_vert"
     1611void Penta::GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){
     1612
     1613        // Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i
     1614
     1615        int i;
     1616                               
     1617        GetNodalFunctions(B, gauss_l1l2l3l4);
     1618
     1619}
     1620
     1621#undef __FUNCT__
     1622#define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
     1623void  Penta::CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type){
     1624
     1625        int i;
     1626
     1627
     1628        /* node data: */
     1629        const int    numgrids=6;
     1630        const int    NDOF1=1;
     1631        const int    numdofs=NDOF1*numgrids;
     1632        double       xyz_list[numgrids][3];
     1633        int          doflist[numdofs];
     1634        int          numberofdofspernode;
     1635       
     1636        /* gaussian points: */
     1637        int     num_gauss,ig;
     1638        double* first_gauss_area_coord  =  NULL;
     1639        double* second_gauss_area_coord =  NULL;
     1640        double* third_gauss_area_coord  =  NULL;
     1641        double* fourth_gauss_vert_coord  =  NULL;
     1642        double* area_gauss_weights      =  NULL;
     1643        double* vert_gauss_weights      =  NULL;
     1644        double  gauss_l1l2l3l4[4];
     1645        int     order_area_gauss;
     1646        int     num_vert_gauss;
     1647        int     num_area_gauss;
     1648        int     ig1,ig2;
     1649        double  gauss_weight1,gauss_weight2;
     1650        double  gauss_weight;
     1651
     1652        /* Jacobian: */
     1653        double Jdet;
     1654
     1655        /*element vector at the gaussian points: */
     1656        double  pe_g[numdofs];
     1657        double  pe_g_gaussian[numdofs];
     1658        double l1l2l3l4l5l6[6];
     1659
     1660        /*Spawning: */
     1661        Tria* tria=NULL;
     1662
     1663        /*input parameters for structural analysis (diagnostic): */
     1664        double* velocity_param=NULL;
     1665        double vx_list[numgrids];
     1666        double vy_list[numgrids];
     1667        double du[3];
     1668        double dv[3];
     1669        double dudx,dvdy;
     1670
     1671        /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the
     1672         *diagnostic base vertical stifness: */
     1673        tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
     1674        tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
     1675        delete tria;
     1676
     1677        /*Now, handle the standard penta element: */
     1678        /* Get node coordinates and dof list: */
     1679        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1680        GetDofList(&doflist[0],&numberofdofspernode);
     1681
     1682        /* Set pe_g to 0: */
     1683        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     1684
     1685        /* recover input parameters: */
     1686        velocity_param=ParameterInputsRecover(inputs,"velocity");
     1687
     1688        if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
     1689
     1690        /*Initialize vx_list and vy_list: */
     1691        for(i=0;i<numgrids;i++){
     1692                /*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity
     1693                 *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
     1694                vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0];
     1695                vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
     1696        }
     1697
     1698        /*Get gaussian points and weights :*/
     1699        order_area_gauss=2;
     1700        num_vert_gauss=2;
     1701
     1702        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);
     1703        #ifdef _DEBUG_
     1704        for (i=0;i<num_area_gauss;i++){
     1705                _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));
     1706        }
     1707        for (i=0;i<num_vert_gauss;i++){
     1708                _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
     1709        }
     1710        #endif
     1711
     1712        /* Start  looping on the number of gaussian points: */
     1713        for (ig1=0; ig1<num_area_gauss; ig1++){
     1714                for (ig2=0; ig2<num_vert_gauss; ig2++){
     1715               
     1716                        /*Pick up the gaussian point: */
     1717                        gauss_weight1=*(area_gauss_weights+ig1);
     1718                        gauss_weight2=*(vert_gauss_weights+ig2);
     1719                        gauss_weight=gauss_weight1*gauss_weight2;
     1720                       
     1721                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
     1722                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
     1723                        gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
     1724                        gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
     1725       
     1726                        /*Get velocity derivative, with respect to x and y: */
     1727                        GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
     1728                        GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
     1729                        dudx=du[0];
     1730                        dvdy=dv[1];
     1731                       
     1732
     1733                        /* Get Jacobian determinant: */
     1734                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
     1735                        #ifdef _DEBUG_
     1736                        _printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet);
     1737                        #endif
     1738               
     1739                         /*Get nodal functions: */
     1740                        GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
     1741
     1742                        /*Build pe_g_gaussian vector: */
     1743                        for (i=0;i<numgrids;i++){
     1744                                pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss_weight*l1l2l3l4l5l6[i];
     1745                        }
     1746
     1747                        /*Add pe_g_gaussian vector to pe_g: */
     1748                        for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     1749
     1750                } //for (ig2=0; ig2<num_vert_gauss; ig2++)
     1751        } //for (ig1=0; ig1<num_area_gauss; ig1++)
     1752
     1753        /*Add pe_g to global vector pg: */
     1754        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     1755
     1756        cleanup_and_return:
     1757        xfree((void**)&first_gauss_area_coord);
     1758        xfree((void**)&second_gauss_area_coord);
     1759        xfree((void**)&third_gauss_area_coord);
     1760        xfree((void**)&fourth_gauss_vert_coord);
     1761        xfree((void**)&area_gauss_weights);
     1762        xfree((void**)&vert_gauss_weights);
     1763
     1764
     1765}

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

    r111 r128  
    3737                double b[6];
    3838                double k[6];
     39                double melting[6];
     40                double accumulation[6];
    3941                int    friction_type;
    4042                double p;
    … …  
    4648                double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    4749                int    collapse;
    48                 double melting[6];
    49                 double accumulation[6];
    5050                double geothermalflux[6];
    5151                int    artdiff;
    … …  
    7373                void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    7474                void  CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type);
     75                void  CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type);
    7576                void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
    7677                void  UpdateFromInputs(ParameterInputs* inputs);
    … …  
    9495                void  GetB(double* pB, double* xyz_list, double* gauss_l1l2l3l4);
    9596                void  GetBPrime(double* B, double* xyz_list, double* gauss_l1l2l3l4);
     97                void  GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
     98                void  GetBPrime_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4);
    9699                void  GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3l4);
    97100                void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4);
    … …  
    100103                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
    101104                void  CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type);
     105                void  CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type);
    102106                void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
    103107                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);

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

    r98 r128  
    3232}
    3333
    34 Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],
    35                                 int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,double tria_viscosity_overshoot){
     34Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3],
     35                                double tria_accumulation[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
     36                                double tria_viscosity_overshoot){
    3637       
    3738        int i;
    … …  
    4849                b[i]=tria_b[i];
    4950                k[i]=tria_k[i];
     51                melting[i]=tria_melting[i];
     52                accumulation[i]=tria_accumulation[i];
    5053        }
    5154        matice=NULL;
    … …  
    8285        printf("   b=[%g,%g,%g]\n",b[0],b[1],b[2]);
    8386        printf("   k=[%g,%g,%g]\n",k[0],k[1],k[2]);
     87        printf("   melting=[%g,%g,%g]\n",melting[0],melting[1],melting[2]);
     88        printf("   accumulation=[%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2]);
    8489        printf("   friction_type: %i\n",friction_type);
    8590        printf("   p: %g\n",p);
    … …  
    129134        memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b);
    130135        memcpy(marshalled_dataset,&k,sizeof(k));marshalled_dataset+=sizeof(k);
     136        memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting);
     137        memcpy(marshalled_dataset,&accumulation,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
    131138        memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
    132139        memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
    … …  
    157164                +sizeof(b)
    158165                +sizeof(k)
     166                +sizeof(melting)
     167                +sizeof(accumulation)
    159168                +sizeof(friction_type)
    160169                +sizeof(onbed)
    … …  
    197206        memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
    198207        memcpy(&k,marshalled_dataset,sizeof(k));marshalled_dataset+=sizeof(k);
     208        memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting);
     209        memcpy(&accumulation,marshalled_dataset,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation);
    199210        memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
    200211        memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
    … …  
    20422053}
    20432054
     2055
     2056
     2057#undef __FUNCT__
     2058#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
     2059void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     2060
     2061        int i,j;
     2062
     2063        /* node data: */
     2064        const int    numgrids=3;
     2065        const int    NDOF1=1;
     2066        const int    numdofs=NDOF1*numgrids;
     2067        double       xyz_list[numgrids][3];
     2068        int          doflist[numdofs];
     2069        int          numberofdofspernode;
     2070       
     2071        /* gaussian points: */
     2072        int     num_gauss,ig;
     2073        double* first_gauss_area_coord  =  NULL;
     2074        double* second_gauss_area_coord =  NULL;
     2075        double* third_gauss_area_coord  =  NULL;
     2076        double* gauss_weights           =  NULL;
     2077        double  gauss_weight;
     2078        double  gauss_l1l2l3[3];
     2079
     2080
     2081        /* matrices: */
     2082        double L[1][3];
     2083        double DL_scalar;
     2084
     2085        double Ke_gg[3][3]; //stiffness matrix
     2086        double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
     2087        double Jdet;
     2088
     2089        /* Get node coordinates and dof list: */
     2090        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2091        GetDofList(&doflist[0],&numberofdofspernode);
     2092
     2093        /* Set Ke_gg to 0: */
     2094        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     2095
     2096        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2097        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2098
     2099        /* Start  looping on the number of gaussian points: */
     2100        for (ig=0; ig<num_gauss; ig++){
     2101                /*Pick up the gaussian point: */
     2102                gauss_weight=*(gauss_weights+ig);
     2103                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2104                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2105                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2106
     2107                /* Get Jacobian determinant: */
     2108                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2109       
     2110                //Get L matrix if viscous basal drag present:
     2111                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2112
     2113                DL_scalar=gauss_weight*Jdet;
     2114               
     2115                /*  Do the triple producte tL*D*L: */
     2116                TripleMultiply( &L[0][0],1,3,1,
     2117                                        &DL_scalar,1,1,0,
     2118                                        &L[0][0],1,3,0,
     2119                                        &Ke_gg_gaussian[0][0],0);
     2120
     2121                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     2122                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2123        } //for (ig=0; ig<num_gauss; ig++
     2124
     2125        /*Add Ke_gg to global matrix Kgg: */
     2126        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     2127               
     2128        cleanup_and_return:
     2129        xfree((void**)&first_gauss_area_coord);
     2130        xfree((void**)&second_gauss_area_coord);
     2131        xfree((void**)&third_gauss_area_coord);
     2132        xfree((void**)&gauss_weights);
     2133}
     2134
     2135#undef __FUNCT__
     2136#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
     2137void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     2138
     2139        int i,j;
     2140       
     2141        /* node data: */
     2142        const int    numgrids=3;
     2143        const int    NDOF1=1;
     2144        const int    numdofs=NDOF1*numgrids;
     2145        double       xyz_list[numgrids][3];
     2146        int          doflist[numdofs];
     2147        int          numberofdofspernode;
     2148
     2149        /* gaussian points: */
     2150        int     num_gauss,ig;
     2151        double* first_gauss_area_coord  =  NULL;
     2152        double* second_gauss_area_coord =  NULL;
     2153        double* third_gauss_area_coord  =  NULL;
     2154        double* gauss_weights           =  NULL;
     2155        double  gauss_weight;
     2156        double  gauss_l1l2l3[3];
     2157
     2158       
     2159        /* surface normal: */
     2160        double x4,y4,z4;
     2161        double x5,y5,z5;
     2162        double x6,y6,z6;
     2163        double v46[3];
     2164        double v56[3];
     2165        double normal[3];
     2166        double norm_normal;
     2167        double nz;
     2168
     2169        /*Matrices: */
     2170        double DL_scalar;
     2171        double L[3];
     2172        double Jdet;
     2173
     2174        /* local element matrices: */
     2175        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
     2176        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
     2177       
     2178        /* Get node coordinates and dof list: */
     2179        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2180        GetDofList(&doflist[0],&numberofdofspernode);
     2181
     2182        /* Set Ke_gg to 0: */
     2183        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     2184
     2185        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2186        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2187
     2188        /*Build normal vector to the surface:*/
     2189       
     2190        x4=xyz_list[0][0];
     2191        y4=xyz_list[0][1];
     2192        z4=xyz_list[0][2];
     2193
     2194        x5=xyz_list[1][0];
     2195        y5=xyz_list[1][1];
     2196        z5=xyz_list[1][2];
     2197
     2198        x6=xyz_list[2][0];
     2199        y6=xyz_list[2][1];
     2200        z6=xyz_list[2][2];
     2201
     2202        v46[0]=x4-x6;
     2203        v46[1]=y4-y6;
     2204        v46[2]=z4-z6;
     2205
     2206        v56[0]=x5-x6;
     2207        v56[1]=y5-y6;
     2208        v56[2]=z5-z6;
     2209
     2210        normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6);
     2211        normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6);
     2212        normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6);
     2213
     2214        norm_normal=sqrt(pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2));
     2215        nz=1.0/norm_normal*normal[2];
     2216
     2217        /* Start  looping on the number of gaussian points: */
     2218        for (ig=0; ig<num_gauss; ig++){
     2219                /*Pick up the gaussian point: */
     2220                gauss_weight=*(gauss_weights+ig);
     2221                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2222                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2223                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2224
     2225                /* Get Jacobian determinant: */
     2226                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2227       
     2228                //Get L matrix if viscous basal drag present:
     2229                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2230
     2231                /**********************Do not forget the sign**********************************/
     2232                DL_scalar=- gauss_weight*Jdet*nz;
     2233                /******************************************************************************/
     2234       
     2235                /*  Do the triple producte tL*D*L: */
     2236                TripleMultiply( L,1,3,1,
     2237                                        &DL_scalar,1,1,0,
     2238                                        L,1,3,0,
     2239                                        &Ke_gg_gaussian[0][0],0);
     2240
     2241                /* Add the Ke_gg_gaussian, onto Ke_gg: */
     2242                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2243
     2244               
     2245        } //for (ig=0; ig<num_gauss; ig++)
     2246
     2247        /*Add Ke_gg to global matrix Kgg: */
     2248        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     2249       
     2250        cleanup_and_return:
     2251        xfree((void**)&first_gauss_area_coord);
     2252        xfree((void**)&second_gauss_area_coord);
     2253        xfree((void**)&third_gauss_area_coord);
     2254        xfree((void**)&gauss_weights);
     2255}
     2256               
     2257
     2258#undef __FUNCT__
     2259#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
     2260void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type){
     2261
     2262        int             i,j;
     2263
     2264        /* node data: */
     2265        const int    numgrids=3;
     2266        const int    NDOF1=1;
     2267        const int    numdofs=NDOF1*numgrids;
     2268        double       xyz_list[numgrids][3];
     2269        int          doflist[numdofs];
     2270        int          numberofdofspernode;
     2271       
     2272        /* gaussian points: */
     2273        int     num_gauss,ig;
     2274        double* first_gauss_area_coord  =  NULL;
     2275        double* second_gauss_area_coord =  NULL;
     2276        double* third_gauss_area_coord  =  NULL;
     2277        double* gauss_weights           =  NULL;
     2278        double  gauss_weight;
     2279        double  gauss_l1l2l3[3];
     2280
     2281        /* Jacobian: */
     2282        double Jdet;
     2283
     2284        /*nodal functions: */
     2285        double l1l2l3[3];
     2286
     2287        /*element vector at the gaussian points: */
     2288        double  pe_g[numdofs];
     2289        double  pe_g_gaussian[numdofs];
     2290
     2291        /* matrices: */
     2292        double L[numgrids];
     2293
     2294        /*input parameters for structural analysis (diagnostic): */
     2295        double* velocity_param=NULL;
     2296        double  vx_list[numgrids];
     2297        double  vy_list[numgrids];
     2298        double  vx,vy;
     2299        double  meltingvalue;
     2300        double  slope[2];
     2301        double  dbdx,dbdy;
     2302
     2303        /* recover input parameters: */
     2304        velocity_param=ParameterInputsRecover(inputs,"velocity");
     2305        if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
     2306
     2307        /* Get node coordinates and dof list: */
     2308        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2309        GetDofList(&doflist[0],&numberofdofspernode);
     2310
     2311        /* Set pe_g to 0: */
     2312        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
     2313
     2314        /*Initialize vx_list and vy_list: */
     2315        for(i=0;i<numgrids;i++){
     2316                if(velocity_param){
     2317                        /*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity
     2318                         *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
     2319                        vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0];
     2320                        vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
     2321                }
     2322                else{
     2323                        vx_list[i]=0;
     2324                        vy_list[i]=0;
     2325                }
     2326        }
     2327               
     2328        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2329        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2330
     2331        /*For icesheets: */
     2332        /* Start  looping on the number of gaussian points: */
     2333        for (ig=0; ig<num_gauss; ig++){
     2334
     2335                /*Pick up the gaussian point: */
     2336                gauss_weight=*(gauss_weights+ig);
     2337                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2338                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2339                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2340
     2341                /*Get melting at gaussian point: */
     2342                GetParameterValue(&meltingvalue, &melting[0],gauss_l1l2l3);
     2343
     2344                /*Get velocity at gaussian point: */
     2345                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     2346                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     2347
     2348                /*Get bed slope: */
     2349                GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
     2350                dbdx=slope[0];
     2351                dbdy=slope[1];
     2352
     2353                /* Get Jacobian determinant: */
     2354                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2355               
     2356                //Get L matrix if viscous basal drag present:
     2357                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2358
     2359               
     2360                /*Build gaussian vector: */
     2361                for(i=0;i<numgrids;i++){
     2362                        pe_g_gaussian[i]=Jdet*gauss_weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
     2363                }
     2364       
     2365                /*Add pe_g_gaussian vector to pe_g: */
     2366                for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];
     2367       
     2368        }
     2369
     2370        /*Add pe_g to global vector pg: */
     2371        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
     2372
     2373        cleanup_and_return:
     2374        xfree((void**)&first_gauss_area_coord);
     2375        xfree((void**)&second_gauss_area_coord);
     2376        xfree((void**)&third_gauss_area_coord);
     2377        xfree((void**)&gauss_weights);
     2378
     2379}

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

    r91 r128  
    3636                double b[3];
    3737                double k[3];
     38                double melting[3];
     39                double accumulation[3];
    3840                int    friction_type;
    3941                double p;
    … …  
    4850
    4951                Tria();
    50                 Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],
     52                Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],double melting[3],double accumulation[3],
    5153                                int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot);
    5254                ~Tria();
    … …  
    6870                void  CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    6971                void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type);
     72                void  CreateKMatrixDiagnosticBaseVert(Mat Kgg,ParameterInputs* inputs,int analysis_type);
     73                void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    7074                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    7175                void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
    … …  
    8791
    8892                void  CreatePVectorDiagnosticHoriz(Vec pg,ParameterInputs* inputs,int analysis_type);
     93                void  CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type);
    8994                Matpar* GetMatPar();
    9095                int   GetShelf();

  • issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp

    r117 r128  
    1919                numdofs=2;
    2020        }
     21        else if (analysis_type==DiagnosticVertAnalysisEnum()){
     22                numdofs=1;
     23        }
    2124        else if (analysis_type==ControlAnalysisEnum()){
    2225                numdofs=2;

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r104 r128  
    44        t1=clock;
    55
    6         analysis_dh='diagnostic_horiz';
    7         analysis_dbv='diagnostic_basevert';
    8         analysis_dv='diagnostic_vert';
    9 
    106        %Build all models requested for diagnostic simulation
    117        %remark, partitions are all identical.
    12         md.analysis_type=analysis_dh;
    13         m_dh=CreateFemModel(md);
     8        md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
    149
    1510        if strcmpi(md.type,'3d'),
    16                 md.analysis_type=analysis_dbv;
    17                 m_dbv=CreateFemModel(md);
    18                 md.analysis_type=analysis_dv;
    19                 m_dv=CreateFemModel(md);
     11                md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
    2012        end
    2113
    … …  
    3123       
    3224                %extrude velocities for collapsed penta elements
    33                 u_g=VelocityExtrude(md,u_g);
     25                disp(sprintf('\n%s',['extruding horizontal velocities...']));
     26                u_g=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
    3427
    35                 %Compute depth averaged velocity and add it to inputs
    36                 velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);
     28                disp(sprintf('\n%s',['computing vertical velocities...']));
     29                u_g_vert=diagnostic_core_linear(m_dv,struct('velocity',u_g));
    3730
    38                 %swith u_g and velocity_average to cluster partitioning
    39                 u_g=SwitchPartitioning(u_g,'cluster',part,tpart,[1 2]);
    40                 velocity_average=SwitchPartitioning(velocity_average,'cluster',part,tpart,[1 2]);
    41                
    42                 disp(sprintf('\n%s',['computing basal vertical velocities...']));
    43                 SetUset(m_dbv);
    44                 u_g_basevert=diagnostic_core_linear(m_dbv,m_dbv.parameters,struct('velocity',u_g,'velocity_average',velocity_average),analysis_dbv);
    45                
    46                 %use u_g_basevert to constrain vertical velocity
    47                 SetUset(m_dv);
    48                 m_dv.y_g=u_g_basevert;
    49 
    50                 %reduce constraining vector to s-set (the one we solve on)
    51                 m_dv.y_s = Reducevectorg( m_dv.y_g);
    52 
    53                 %run core linear to solve for vertical velocity
    54                 disp(sprintf('\n%s',['computing vertical velocities...']));
    55                 u_g_vert=diagnostic_core_linear(m_dv,m_dv.parameters,struct('velocity',u_g),analysis_dv);
    56                
    57                 %load results onto model: carefull, u_g and u_g_vert are in cluster partitioning
    58                 md.vx=u_g(indx)*md.yts;
    59                 md.vy=u_g(indy)*md.yts;
    60                 md.vz=u_g_vert(indz)*md.yts;
     31                %load results onto model:
     32                md.vx=u_g(1:2:end)*md.yts;
     33                md.vy=u_g(2:2:end)*md.yts;
     34                md.vz=u_g_vert*md.yts;
    6135                md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
    6236               
    … …  
    6539        else
    6640               
    67                 %Build partitioning vectors to recover solution
    68                 indx=m_dh.part(1:2:end);
    69                 indy=m_dh.part(2:2:end);
    70 
    7141                %load results onto model
    72                 md.vx=u_g(indx)*md.yts;
    73                 md.vy=u_g(indy)*md.yts;
    74                 md.vz=zeros(md.numberofgrids,1);
     42                md.vx=u_g(1:2:end)*md.yts;
     43                md.vy=u_g(2:2:end)*md.yts;
     44                md.vz=zeros(m_dv.nodesets.gsize,1);
    7545                md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
    7646       

  • issm/trunk/src/mex/Makefile.am

    r9 r128  
    3838                                TriMeshProcessRifts\
    3939                                TriMeshRefine\
    40                                 UpdateFromInputs
     40                                UpdateFromInputs\
     41                                VelocityExtrude
    4142
    4243endif
    … …  
    166167                          UpdateFromInputs/UpdateFromInputs.h
    167168
     169VelocityExtrude_SOURCES = VelocityExtrude/VelocityExtrude.cpp\
     170                          VelocityExtrude/VelocityExtrude.h
     171

  • issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.cpp

    r82 r128  
    3535        /*write output : */
    3636        WriteData(UGOUT,ug,0,0,"Vector",NULL);
    37 
    3837        /*Free ressources: */
    3938        delete elements;
    … …  
    4241        delete materials;
    4342        VecFree(&ug);
    44 
     43       
    4544        /*end module: */
    4645        MODULEEND();
     46
    4747}
    4848

  • issm/trunk/src/mex/VelocityExtrude/VelocityExtrude.h

    r82 r128  
    1818/* serial input macros: */
    1919#define ELEMENTS (mxArray*)prhs[0]
    20 #define LOADS (mxArray*)prhs[1]
    21 #define NODES (mxArray*)prhs[2]
     20#define NODES (mxArray*)prhs[1]
     21#define LOADS (mxArray*)prhs[2]
    2222#define MATERIALS (mxArray*)prhs[3]
    2323#define UG (mxArray*)prhs[4]

  • issm/trunk/todo

    r53 r128  
    11Create sorting routine for ids.
    2 Update all inputs. Look at that in more details, ie with respect to control methods. When do we do update of inputs?
    3 Handle input update for Nodes and Materials.
     2Update all inputs. Look at that in more details, ie with respect to control methods. When do we do update of inputs?  Handle input update for Nodes and Materials.
    43Define default cluster in /etc/cielo.rc
     4VelocityExtrude: debug
     5Add vx and vy as input parameters to diagnostic.
Note: See TracChangeset for help on using the changeset viewer.