Changeset 3582


Ignore:
Timestamp:
04/20/10 16:35:48 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added discontinuous Galerkin for Balanced thickness solution

Location:
issm/trunk/src
Files:
24 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp

    r3567 r3582  
    1111
    1212void    CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
    13 
    14         int i;
    15         int count=0;
     13       
    1614        DataSet* constraints = NULL;
    17 
     15       
    1816        /*Create constraints: */
    1917        constraints = new DataSet(ConstraintsEnum);
    20 
    21         /*Fetch data: */
    22         IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
    23 
    24         count=1; //matlab indexing
    25         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    26         for (i=0;i<iomodel->numberofvertices;i++){
    27                 if(iomodel->my_vertices[i]){
    28                        
    29                         if ((int)iomodel->spcthickness[2*i]){
    30                
    31                                 constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));//we enforce first translation degree of freedom, for temperature
    32                                 count++;
    33                         }
    34                 }
    35         }
    36 
    37         /*Free data: */
    38         xfree((void**)&iomodel->spcthickness);
    39 
    40         /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    41          * datasets, it will not be redone: */
    42         constraints->Presort();
    43 
    44         cleanup_and_return:
    4518       
    4619        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp

    r3567 r3582  
    88#include "../../objects/objects.h"
    99#include "../../shared/shared.h"
     10#include "../../include/macros.h"
    1011#include "../../include/typedefs.h"
    1112#include "../../MeshPartitionx/MeshPartitionx.h"
    1213#include "../IoModel.h"
    1314
    14 void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
     15void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1516
    1617        /*Intermediary*/
    17         int i,j,k,n;
     18        int i,j;
     19        int vertex_index;
     20        int node_index;
    1821
    1922        /*DataSets: */
    20         DataSet*    elements  = NULL;
    21         DataSet*    nodes = NULL;
    22         DataSet*    vertices = NULL;
    23         DataSet*    materials = NULL;
     23        DataSet* elements  = NULL;
     24        DataSet* nodes = NULL;
     25        DataSet* vertices = NULL;
     26        DataSet* materials = NULL;
    2427
    2528        /*First create the elements, nodes and material properties: */
    … …  
    3235        Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
    3336
     37        /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
    3438        /*2d mesh: */
    3539        if (strcmp(iomodel->meshtype,"2d")==0){
    … …  
    5458                        }
    5559                }//for (i=0;i<numberofelements;i++)
    56 
     60       
    5761                /*Free data : */
    58                 xfree((void**)&iomodel->elements);
    5962                xfree((void**)&iomodel->thickness);
    6063                xfree((void**)&iomodel->surface);
    … …  
    6568        }
    6669        else{ //        if (strcmp(meshtype,"2d")==0)
    67 
    68                 /*Fetch data needed: */
    69                 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    70                 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
    71                 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
    72                 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
    73                 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
    74                 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
    75                 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
    76                 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
    77        
    78                 for (i=0;i<iomodel->numberofelements;i++){
    79                         if(iomodel->my_elements[i]){
    80                                 /*Create and add penta element to elements dataset: */
    81                                 elements->AddObject(new Penta(i,iomodel));
    82 
    83                                 /*Create and add material property to materials dataset: */
    84                                 materials->AddObject(new Matice(i,iomodel,6));
    85                         }
    86                 }//for (i=0;i<numberofelements;i++)
    87 
    88                 /*Free data: */
    89                 xfree((void**)&iomodel->elements);
    90                 xfree((void**)&iomodel->thickness);
    91                 xfree((void**)&iomodel->surface);
    92                 xfree((void**)&iomodel->bed);
    93                 xfree((void**)&iomodel->elementoniceshelf);
    94                 xfree((void**)&iomodel->elementonbed);
    95                 xfree((void**)&iomodel->elementonsurface);
    96                 xfree((void**)&iomodel->elementonwater);
    97 
     70                ISSMERROR("not implemented yet");
    9871        } //if (strcmp(meshtype,"2d")==0)
    9972
    100         /*Add new constrant material property to materials, at the end: */
     73        /*Add new constrant material property tgo materials, at the end: */
    10174        materials->AddObject(new Matpar(iomodel));
    10275
    103         /*First fetch data: */
    104         if (strcmp(iomodel->meshtype,"3d")==0){
    105                 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
    106                 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
    107         }
     76        /*Create nodes and vertices: */
    10877        IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
    10978        IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
    … …  
    11382        IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
    11483        IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
     84        IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
    11585        IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
    11686        IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
     87        if (strcmp(iomodel->meshtype,"3d")==0){
     88                IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
     89                IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
     90        }
    11791
     92        /*Build Vertices dataset*/
    11893        for (i=0;i<iomodel->numberofvertices;i++){
    11994
    … …  
    12499                        vertices->AddObject(new Vertex(i,iomodel));
    125100
    126                         /*Add node to nodes dataset: */
    127                         nodes->AddObject(new Node(i,iomodel));
     101                }
     102        }
    128103
     104        /*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
     105        for (i=0;i<iomodel->numberofelements;i++){
     106                for (j=0;j<3;j++){
     107
     108                        if(iomodel->my_nodes[3*i+j]){
     109
     110                                //Get index of the vertex on which the current node is located
     111                                vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
     112                                ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);
     113
     114                                //Compute Node index (id-1)
     115                                node_index=3*i+j;
     116
     117                                /*Add node to nodes dataset: */
     118                                nodes->AddObject(new Node(vertex_index,node_index,iomodel));
     119
     120                        }
    129121                }
    130122        }
    … …  
    139131        xfree((void**)&iomodel->gridonbed);
    140132        xfree((void**)&iomodel->gridonsurface);
     133        xfree((void**)&iomodel->gridonhutter);
    141134        xfree((void**)&iomodel->uppernodes);
    142135        xfree((void**)&iomodel->gridonicesheet);
    … …  
    146139         * datasets, it will not be redone: */
    147140        elements->Presort();
     141        nodes->Presort();
    148142        vertices->Presort();
    149         nodes->Presort();
    150143        materials->Presort();
    151144
    … …  
    157150        *pvertices=vertices;
    158151        *pmaterials=materials;
    159 
    160152}

  • issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp

    r3567 r3582  
    88#include "../../shared/shared.h"
    99#include "../../include/macros.h"
     10#include "../../include/typedefs.h"
    1011#include "../IoModel.h"
    11 
    1212
    1313void    CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1414
    15         DataSet*    loads    = NULL;
     15        /*Intermediary*/
     16        int i,j;
     17        int i1,i2;
     18        int pos1,pos2;
     19        double e1,e2;
     20
     21        /*Output*/
     22        DataSet* loads=NULL;
     23
     24        /*numericalflux intermediary data: */
     25        char numericalflux_type[NUMERICALFLUXSTRING];
     26        int  numericalflux_id;
     27        int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
     28        int  numericalflux_elem_id;
     29        double numericalflux_h[MAX_NUMERICALFLUX_NODES];
    1630
    1731        /*Create loads: */
    1832        loads   = new DataSet(LoadsEnum);
    19        
    20        
     33
     34        /*Get edges and elements*/
     35        IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
     36        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     37        IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
     38
     39        /*First load data:*/
     40        for (i=0;i<iomodel->numberofedges;i++){
     41
     42                /*Get left and right elements*/
     43                e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     44                e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
     45
     46                /*Now, if this element is not in the partition, pass: */
     47                if(!iomodel->my_elements[(int)e1]) continue;
     48
     49                /*Create load*/
     50                numericalflux_id=i+1; //Matlab indexing
     51                numericalflux_elem_id=(int)e1+1;//id is in matlab index
     52
     53                /*1: Get vertices ids*/
     54                i1=(int)iomodel->edges[4*i+0];
     55                i2=(int)iomodel->edges[4*i+1];
     56
     57                if (!isnan(e2)){
     58                        strcpy(numericalflux_type,"internal");
     59
     60                        /*Now, we must get the nodes of the 4 nodes located on the edge*/
     61
     62                        /*2: Get the column where these ids are located in the index*/
     63                        pos1=pos2=UNDEF;
     64                        for(j=0;j<3;j++){
     65                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
     66                                if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1;
     67                        }
     68                        ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
     69
     70                        /*3: We have the id of the elements and the position of the vertices in the index
     71                         * we can compute their dofs!*/
     72                        numericalflux_node_ids[0]=3*(int)e1+pos1;       //ex: 1 2 3
     73                        numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1
     74                        numericalflux_node_ids[2]=3*(int)e2+pos2;           //ex: 1 2 3
     75                        numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2
     76
     77                        numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1];
     78                        numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
     79                        numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1];
     80                        numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1];
     81                }
     82                else{
     83                        strcpy(numericalflux_type,"boundary");
     84
     85                        /*2: Get the column where these ids are located in the index*/
     86                        pos1==UNDEF;
     87                        for(j=0;j<3;j++){
     88                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
     89                        }
     90                        ISSMASSERT(pos1!=UNDEF);
     91
     92                        /*3: We have the id of the elements and the position of the vertices in the index
     93                         * we can compute their dofs!*/
     94                        numericalflux_node_ids[0]=3*(int)e1+pos1;
     95                        numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1;
     96
     97                        numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1];
     98                        numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
     99                        numericalflux_h[2]=UNDEF;
     100                        numericalflux_h[3]=UNDEF;
     101                }
     102
     103                loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h));
     104        }
     105
     106        /*Free data: */
     107        xfree((void**)&iomodel->edges);
     108        xfree((void**)&iomodel->elements);
     109        xfree((void**)&iomodel->thickness);
     110
    21111        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    22112         * datasets, it will not be redone: */
    23113        loads->Presort();
    24 
    25         cleanup_and_return:
    26114
    27115        /*Assign output pointer: */
    … …  
    29117
    30118}
    31 
    32 

  • issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp

    r3462 r3582  
    1  /* \brief driver for creating parameters dataset, for prognostic analysis.
     1/*!\file: CreateParametersPrognostic.cpp
     2 * \brief driver for creating parameters dataset, for prognostic analysis.
    23 */
    34
    … …  
    1516        int      count;
    1617        int      i;
    17         double* u_g=NULL;
     18        int      dim;
     19        double*  vx_g=NULL;
     20        double*  vy_g=NULL;
    1821
    1922        /*recover parameters : */
    … …  
    2528        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    2629        IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
    27         IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz");
    2830
    29         u_g=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double));
     31        vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
     32        vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
    3033
    31         if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+0]=iomodel->vx[i]/iomodel->yts;
    32         if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+1]=iomodel->vy[i]/iomodel->yts;
    33         if(iomodel->vz) for(i=0;i<iomodel->numberofvertices;i++)u_g[3*i+2]=iomodel->vz[i]/iomodel->yts;
     34        if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts;
     35        if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts;
    3436
    3537        count++;
    36         param= new Param(count,"u_g",DOUBLEVEC);
    37         param->SetDoubleVec(u_g,3*iomodel->numberofvertices,3);
     38        param= new Param(count,"vx_g",DOUBLEVEC);
     39        param->SetDoubleVec(vx_g,iomodel->numberofvertices,1);
    3840        parameters->AddObject(param);
    39 
     41        count++;
     42        param= new Param(count,"vy_g",DOUBLEVEC);
     43        param->SetDoubleVec(vy_g,iomodel->numberofvertices,1);
     44        parameters->AddObject(param);
    4045
    4146        xfree((void**)&iomodel->vx);
    4247        xfree((void**)&iomodel->vy);
    43         xfree((void**)&iomodel->vz);
    44         xfree((void**)&u_g);
     48        xfree((void**)&vx_g);
     49        xfree((void**)&vy_g);
     50
     51        /*Get thickness: */
     52        IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
     53
     54        count++;
     55        param= new Param(count,"h_g",DOUBLEVEC);
     56        if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1);
     57        else param->SetDoubleVec(iomodel->thickness,0,0);
     58        parameters->AddObject(param);
     59
     60        /*Free thickness: */
     61        xfree((void**)&iomodel->thickness);
     62
     63        /*Get dhdt: */
     64        IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
     65        if(iomodel->dhdt) for(i=0;i<iomodel->numberofvertices;i++)iomodel->dhdt[i]=iomodel->dhdt[i]/iomodel->yts;
     66
     67        count++;
     68        param= new Param(count,"dhdt_g",DOUBLEVEC);
     69        if(iomodel->dhdt) param->SetDoubleVec(iomodel->dhdt,iomodel->numberofvertices,1);
     70        else param->SetDoubleVec(iomodel->dhdt,0,1);
     71        parameters->AddObject(param);
     72
     73        /*Free dhdt: */
     74        xfree((void**)&iomodel->dhdt);
    4575
    4676        /*Get melting: */

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

    r3567 r3582  
    245245        count++;
    246246        param= new Param(count,"numberofnodes",DOUBLE);
    247         if (iomodel->analysis_type==Prognostic2AnalysisEnum) param->SetDouble(3*iomodel->numberofelements);
     247        if (
     248                                iomodel->analysis_type==Prognostic2AnalysisEnum ||
     249                                iomodel->analysis_type==BalancedthicknessAnalysisEnum
     250                                )
     251         param->SetDouble(3*iomodel->numberofelements);
    248252        else param->SetDouble(iomodel->numberofvertices);
    249253        parameters->AddObject(param);

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

    r3417 r3582  
    171171        /*!basal: */
    172172        iomodel->accumulation=NULL;
     173        iomodel->dhdt=NULL;
    173174       
    174175        /*parameter output: */
    … …  
    254255        xfree((void**)&iomodel->melting);
    255256        xfree((void**)&iomodel->accumulation);
     257        xfree((void**)&iomodel->dhdt);
    256258        xfree((void**)&iomodel->B);
    257259        xfree((void**)&iomodel->n);

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

    r3446 r3582  
    172172        double*  melting;
    173173        double*  accumulation;
     174        double*  dhdt;
    174175
    175176        /*parameter output: */

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

    r3570 r3582  
    2222void  Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){
    2323       
    24         if (iomodel->analysis_type==Prognostic2AnalysisEnum)
     24        if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum)
    2525                DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle);
    2626        else

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

    r3570 r3582  
    222222        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
    223223        if(!found)ISSMERROR(" could not find vy_average in inputs!");
    224         found=inputs->Recover("dt",&dt);
    225         if(!found)ISSMERROR(" could not find dt in inputs!");
     224        if (analysis_type==Prognostic2AnalysisEnum){
     225                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
     226        }
     227        else if (analysis_type==BalancedthicknessAnalysisEnum){
     228                /*No transient term is involved*/
     229                dt=1;
     230        }
     231        else{
     232                ISSMERROR("analysis_type %i not supported yet");
     233        }
    226234
    227235        /* Get node coordinates, dof list and normal vector: */
    … …  
    330338        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
    331339        if(!found)ISSMERROR(" could not find vy_average in inputs!");
    332         found=inputs->Recover("dt",&dt);
    333         if(!found)ISSMERROR(" could not find dt in inputs!");
     340        if (analysis_type==Prognostic2AnalysisEnum){
     341                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
     342        }
     343        else if (analysis_type==BalancedthicknessAnalysisEnum){
     344                /*No transient term is involved*/
     345                dt=1;
     346        }
     347        else{
     348                ISSMERROR("analysis_type %i not supported yet");
     349        }
    334350
    335351        /* Get node coordinates, dof list and normal vector: */
    … …  
    464480        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
    465481        if(!found)ISSMERROR(" could not find vy_average in inputs!");
    466         found=inputs->Recover("dt",&dt);
    467         if(!found)ISSMERROR(" could not find dt in inputs!");
     482        if (analysis_type==Prognostic2AnalysisEnum){
     483                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
     484        }
     485        else if (analysis_type==BalancedthicknessAnalysisEnum){
     486                /*No transient term is involved*/
     487                dt=1;
     488        }
     489        else{
     490                ISSMERROR("analysis_type %i not supported yet");
     491        }
    468492
    469493        /* Get node coordinates, dof list and normal vector: */

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

    r3570 r3582  
    585585                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
    586586        }
    587         else if (analysis_type==BalancedthicknessAnalysisEnum){
    588 
    589                 CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
    590         }
    591         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    592 
    593                 CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
    594         }
    595587        else if (analysis_type==ThermalAnalysisEnum){
    596588
    … …  
    604596                ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet");
    605597        }
    606 
    607 }
    608 /*}}}*/
    609 /*FUNCTION CreateKMatrixBalancedthickness {{{1*/
    610 
    611 void  Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    612 
    613         /*Collapsed formulation: */
    614         Tria*  tria=NULL;
    615 
    616         /*If on water, skip: */
    617         if(this->properties.onwater)return;
    618 
    619         /*Is this element on the bed? :*/
    620         if(!this->properties.onbed)return;
    621 
    622         /*Spawn Tria element from the base of the Penta: */
    623         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    624         tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
    625         delete tria;
    626         return;
    627 
    628 }
    629 /*}}}*/
    630 /*FUNCTION CreateKMatrixBalancedvelocities {{{1*/
    631 
    632 void  Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    633 
    634         /*Collapsed formulation: */
    635         Tria*  tria=NULL;
    636 
    637         /*If on water, skip: */
    638         if(this->properties.onwater)return;
    639 
    640         /*Is this element on the bed? :*/
    641         if(!this->properties.onbed)return;
    642 
    643         /*Spawn Tria element from the base of the Penta: */
    644         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    645         tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
    646         delete tria;
    647         return;
    648598
    649599}
    … …  
    17111661                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    17121662        }
    1713         else if (analysis_type==BalancedthicknessAnalysisEnum){
    1714 
    1715                 CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
    1716         }
    1717         else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    1718 
    1719                 CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type);
    1720         }
    17211663        else if (analysis_type==ThermalAnalysisEnum){
    17221664
    … …  
    17311673        }       
    17321674
    1733 }
    1734 /*}}}*/
    1735 /*FUNCTION CreatePVectorBalancedthickness {{{1*/
    1736 
    1737 void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    1738 
    1739         /*Collapsed formulation: */
    1740         Tria*  tria=NULL;
    1741 
    1742         /*If on water, skip: */
    1743         if(this->properties.onwater)return;
    1744 
    1745         /*Is this element on the bed? :*/
    1746         if(!this->properties.onbed)return;
    1747 
    1748         /*Spawn Tria element from the base of the Penta: */
    1749         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1750         tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
    1751         delete tria;
    1752         return;
    1753 }
    1754 /*}}}*/
    1755 /*FUNCTION CreatePVectorBalancedvelocities {{{1*/
    1756 
    1757 void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    1758 
    1759         /*Collapsed formulation: */
    1760         Tria*  tria=NULL;
    1761 
    1762         /*If on water, skip: */
    1763         if(this->properties.onwater)return;
    1764 
    1765         /*Is this element on the bed? :*/
    1766         if(!this->properties.onbed)return;
    1767 
    1768         /*Spawn Tria element from the base of the Penta: */
    1769         tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    1770         tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
    1771         delete tria;
    1772         return;
    17731675}
    17741676/*}}}*/

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

    r3529 r3582  
    110110                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    111111                void  CreatePVectorPrognostic( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    112                 void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    113                 void  CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    114                 void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    115                 void  CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    116112
    117113                void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);

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

    r3570 r3582  
    7777        /*hooks: */
    7878        //go recover node ids, needed to initialize the node hook.
    79         if (iomodel->analysis_type==Prognostic2AnalysisEnum){
     79        if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum){
    8080                /*Discontinuous Galerkin*/
    8181                tria_node_ids[0]=3*i+1;
    … …  
    633633
    634634        /* matrices: */
     635        double B[2][numgrids];
     636        double Bprime[2][numgrids];
     637        double DL[2][2]={0.0};
     638        double DLprime[2][2]={0.0};
     639        double DL_scalar;
     640        double Ke_gg[numdof][numdof]={0.0};
     641        double Ke_gg2[numdof][numdof]={0.0};
     642        double Jdettria;
     643
     644        /*input parameters for structural analysis (diagnostic): */
     645        double  vx_list[numgrids]={0.0};
     646        double  vy_list[numgrids]={0.0};
     647        double  vx,vy;
     648        int     dofs[1]={0};
     649        int     found;
     650
     651        /*dynamic objects pointed to by hooks: */
     652        Node**  nodes=NULL;
     653        Matpar* matpar=NULL;
     654        Matice* matice=NULL;
     655        Numpar* numpar=NULL;
     656
     657        ParameterInputs* inputs=NULL;
     658
     659        /*recover pointers: */
     660        inputs=(ParameterInputs*)vinputs;
     661
     662        /*recover objects from hooks: */
     663        nodes=(Node**)   hnodes.deliverp();
     664        matpar=(Matpar*)hmatpar.delivers();
     665        matice=(Matice*)hmatice.delivers();
     666        numpar=(Numpar*)hnumpar.delivers();
     667
     668        /*recover extra inputs from users, at current convergence iteration: */
     669        found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
     670        if(!found)ISSMERROR(" could not find vx_average in inputs!");
     671        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
     672        if(!found)ISSMERROR(" could not find vy_average in inputs!");
     673
     674        /* Get node coordinates and dof list: */
     675        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     676        GetDofList(&doflist[0],&numberofdofspernode);
     677
     678        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     679        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     680
     681        /* Start  looping on the number of gaussian points: */
     682        for (ig=0; ig<num_gauss; ig++){
     683                /*Pick up the gaussian point: */
     684                gauss_weight=*(gauss_weights+ig);
     685                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     686                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     687                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     688
     689                /* Get Jacobian determinant: */
     690                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     691
     692                /*Get B  and B prime matrix: */
     693                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
     694                GetB_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     695                GetBPrime_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     696
     697                //Get vx, vy and their derivatives at gauss point
     698                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     699                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     700
     701                DL_scalar=-gauss_weight*Jdettria;
     702
     703                DLprime[0][0]=DL_scalar*vx;
     704                DLprime[1][1]=DL_scalar*vy;
     705
     706                //Do the triple product tL*D*L.
     707                TripleMultiply( &B[0][0],2,numdof,1,
     708                                        &DLprime[0][0],2,2,0,
     709                                        &Bprime[0][0],2,numdof,0,
     710                                        &Ke_gg2[0][0],0);
     711
     712                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     713                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
     714
     715        } // for (ig=0; ig<num_gauss; ig++)
     716
     717        /*Add Ke_gg to global matrix Kgg: */
     718        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     719
     720cleanup_and_return:
     721        xfree((void**)&first_gauss_area_coord);
     722        xfree((void**)&second_gauss_area_coord);
     723        xfree((void**)&third_gauss_area_coord);
     724        xfree((void**)&gauss_weights);
     725
     726}
     727/*}}}*/
     728/*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/
     729void  Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     730
     731        /* local declarations */
     732        int             i,j;
     733
     734        /* node data: */
     735        const int    numgrids=3;
     736        const int    NDOF1=1;
     737        const int    numdof=NDOF1*numgrids;
     738        double       xyz_list[numgrids][3];
     739        int          doflist[numdof];
     740        int          numberofdofspernode;
     741
     742        /* gaussian points: */
     743        int     num_gauss,ig;
     744        double* first_gauss_area_coord  =  NULL;
     745        double* second_gauss_area_coord =  NULL;
     746        double* third_gauss_area_coord  =  NULL;
     747        double* gauss_weights           =  NULL;
     748        double  gauss_weight;
     749        double  gauss_l1l2l3[3];
     750
     751        /* matrices: */
    635752        double L[numgrids];
    636753        double B[2][numgrids];
    … …  
    641758        double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
    642759        double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    643         double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    644         double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    645 
     760        double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     761        double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    646762        double Jdettria;
    647763
    648764        /*input parameters for structural analysis (diagnostic): */
     765        double  surface_normal[3];
     766        double  nx,ny,norm;
    649767        double  vxvy_list[numgrids][2]={0.0};
    650768        double  vx_list[numgrids]={0.0};
    … …  
    690808        GetDofList(&doflist[0],&numberofdofspernode);
    691809
     810        /*Modify z so that it reflects the surface*/
     811        for(i=0;i<numgrids;i++) xyz_list[i][2]=this->properties.s[i];
     812
     813        /*Get normal vector to the surface*/
     814        nx=(vx_list[0]+vx_list[1]+vx_list[2])/3;
     815        ny=(vy_list[0]+vy_list[1]+vy_list[2])/3;
     816        if(nx==0 && ny==0){
     817                SurfaceNormal(&surface_normal[0],xyz_list);
     818                nx=surface_normal[0];
     819                ny=surface_normal[1];
     820        }
     821        if(nx==0 && ny==0){
     822                nx=0;
     823                ny=1;
     824        }
     825        norm=pow( pow(nx,2)+pow(ny,2) , (double).5);
     826        nx=nx/norm;
     827        ny=ny/norm;
     828
    692829        //Create Artificial diffusivity once for all if requested
    693830        if(numpar->artdiff){
    … …  
    700837                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    701838
    702                 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
    703                 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     839                K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
     840                K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
    704841        }
    705842
    … …  
    734871                DL_scalar=gauss_weight*Jdettria;
    735872
    736                 //Create DL and DLprime matrix
    737                 DL[0][0]=DL_scalar*dvxdx;
    738                 DL[1][1]=DL_scalar*dvydy;
    739 
    740                 DLprime[0][0]=DL_scalar*vx;
    741                 DLprime[1][1]=DL_scalar*vy;
     873                DLprime[0][0]=DL_scalar*nx;
     874                DLprime[1][1]=DL_scalar*ny;
    742875
    743876                //Do the triple product tL*D*L.
    744                 //Ke_gg_thickness=B'*DLprime*Bprime;
    745 
    746                 TripleMultiply( &B[0][0],2,numdof,1,
    747                                         &DL[0][0],2,2,0,
    748                                         &B[0][0],2,numdof,0,
    749                                         &Ke_gg_thickness1[0][0],0);
    750 
     877                //Ke_gg_velocities=B'*DLprime*Bprime;
    751878                TripleMultiply( &B[0][0],2,numdof,1,
    752879                                        &DLprime[0][0],2,2,0,
    753880                                        &Bprime[0][0],2,numdof,0,
    754                                         &Ke_gg_thickness2[0][0],0);
     881                                        &Ke_gg_velocities2[0][0],0);
    755882
    756883                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    757                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
    758                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
     884                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
    759885
    760886                if(numpar->artdiff){
    … …  
    822948}
    823949/*}}}*/
    824 /*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/
    825 void  Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     950/*FUNCTION Tria::CreateKMatrixDiagnosticHoriz {{{1*/
     951void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    826952
    827953        /* local declarations */
    … …  
    830956        /* node data: */
    831957        const int    numgrids=3;
    832         const int    NDOF1=1;
    833         const int    numdof=NDOF1*numgrids;
     958        const int    numdof=2*numgrids;
    834959        double       xyz_list[numgrids][3];
    835960        int          doflist[numdof];
    … …  
    845970        double  gauss_l1l2l3[3];
    846971
     972        /* material data: */
     973        double viscosity; //viscosity
     974        double newviscosity; //viscosity
     975        double oldviscosity; //viscosity
     976
     977        /* strain rate: */
     978        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     979        double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/
     980
     981        /* matrices: */
     982        double B[3][numdof];
     983        double Bprime[3][numdof];
     984        double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
     985        double D_scalar;
     986
     987        /* local element matrices: */
     988        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     989        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
     990
     991        double Jdet;
     992
     993        /*input parameters for structural analysis (diagnostic): */
     994        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     995        double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     996        double  thickness;
     997        int     dofs[2]={0,1};
     998
     999        /*dynamic objects pointed to by hooks: */
     1000        Node**  nodes=NULL;
     1001        Matpar* matpar=NULL;
     1002        Matice* matice=NULL;
     1003        Numpar* numpar=NULL;
     1004
     1005        ParameterInputs* inputs=NULL;
     1006
     1007        /*First, if we are on water, return empty matrix: */
     1008        if(this->properties.onwater) return;
     1009
     1010        /*recover pointers: */
     1011        inputs=(ParameterInputs*)vinputs;
     1012
     1013        /*recover objects from hooks: */
     1014        nodes =(Node**) hnodes.deliverp();
     1015        matpar=(Matpar*)hmatpar.delivers();
     1016        matice=(Matice*)hmatice.delivers();
     1017        numpar=(Numpar*)hnumpar.delivers();
     1018
     1019        /*recover extra inputs from users, at current convergence iteration: */
     1020        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     1021        inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     1022
     1023        /* Get node coordinates and dof list: */
     1024        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1025        GetDofList(&doflist[0],&numberofdofspernode);
     1026
     1027        /* Set Ke_gg to 0: */
     1028        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     1029
     1030        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1031        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1032
     1033        /* Start  looping on the number of gaussian points: */
     1034        for (ig=0; ig<num_gauss; ig++){
     1035                /*Pick up the gaussian point: */
     1036                gauss_weight=*(gauss_weights+ig);
     1037                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1038                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1039                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1040
     1041
     1042                /*Compute thickness at gaussian point: */
     1043                GetParameterValue(&thickness, &this->properties.h[0],gauss_l1l2l3);
     1044
     1045                /*Get strain rate from velocity: */
     1046                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
     1047                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
     1048
     1049                /*Get viscosity: */
     1050                matice->GetViscosity2d(&viscosity, &epsilon[0]);
     1051                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
     1052
     1053                /* Get Jacobian determinant: */
     1054                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1055
     1056                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
     1057                        onto this scalar matrix, so that we win some computational time: */
     1058                newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity);
     1059                D_scalar=newviscosity*thickness*gauss_weight*Jdet;
     1060
     1061                for (i=0;i<3;i++){
     1062                        D[i][i]=D_scalar;
     1063                }
     1064
     1065                /*Get B and Bprime matrices: */
     1066                GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     1067                GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     1068
     1069                /*  Do the triple product tB*D*Bprime: */
     1070                TripleMultiply( &B[0][0],3,numdof,1,
     1071                                        &D[0][0],3,3,0,
     1072                                        &Bprime[0][0],3,numdof,0,
     1073                                        &Ke_gg_gaussian[0][0],0);
     1074
     1075                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     1076                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     1077
     1078        } // for (ig=0; ig<num_gauss; ig++)
     1079
     1080        /*Add Ke_gg to global matrix Kgg: */
     1081        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1082
     1083        /*Do not forget to include friction: */
     1084        if(!this->properties.shelf){
     1085                CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
     1086        }
     1087
     1088cleanup_and_return:
     1089        xfree((void**)&first_gauss_area_coord);
     1090        xfree((void**)&second_gauss_area_coord);
     1091        xfree((void**)&third_gauss_area_coord);
     1092        xfree((void**)&gauss_weights);
     1093
     1094}
     1095/*}}}*/
     1096/*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/
     1097void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     1098
     1099
     1100        /* local declarations */
     1101        int             i,j;
     1102
     1103        /* node data: */
     1104        const int    numgrids=3;
     1105        const int    numdof=2*numgrids;
     1106        double       xyz_list[numgrids][3];
     1107        int          doflist[numdof];
     1108        int          numberofdofspernode;
     1109       
     1110        /* gaussian points: */
     1111        int     num_gauss,ig;
     1112        double* first_gauss_area_coord  =  NULL;
     1113        double* second_gauss_area_coord =  NULL;
     1114        double* third_gauss_area_coord  =  NULL;
     1115        double* gauss_weights           =  NULL;
     1116        double  gauss_weight;
     1117        double  gauss_l1l2l3[3];
     1118
     1119        /* matrices: */
     1120        double L[2][numdof];
     1121        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
     1122        double DL_scalar;
     1123
     1124        /* local element matrices: */
     1125        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1126        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
     1127       
     1128        double Jdet;
     1129       
     1130        /*slope: */
     1131        double  slope[2]={0.0,0.0};
     1132        double  slope_magnitude;
     1133
     1134        /*input parameters for structural analysis (diagnostic): */
     1135        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     1136        int     dofs[2]={0,1};
     1137
     1138        /*friction: */
     1139        double alpha2_list[numgrids]={0.0,0.0,0.0};
     1140        double alpha2;
     1141
     1142        double MAXSLOPE=.06; // 6 %
     1143        double MOUNTAINKEXPONENT=10;
     1144
     1145        /*dynamic objects pointed to by hooks: */
     1146        Node**  nodes=NULL;
     1147        Matpar* matpar=NULL;
     1148        Matice* matice=NULL;
     1149        Numpar* numpar=NULL;
     1150
     1151        ParameterInputs* inputs=NULL;
     1152
     1153        /*recover pointers: */
     1154        inputs=(ParameterInputs*)vinputs;
     1155
     1156        /*recover objects from hooks: */
     1157        nodes=(Node**)hnodes.deliverp();
     1158        matpar=(Matpar*)hmatpar.delivers();
     1159        matice=(Matice*)hmatice.delivers();
     1160        numpar=(Numpar*)hnumpar.delivers();
     1161       
     1162        /* Get node coordinates and dof list: */
     1163        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1164        GetDofList(&doflist[0],&numberofdofspernode);
     1165
     1166        /* Set Ke_gg to 0: */
     1167        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     1168
     1169        if (this->properties.shelf){
     1170                /*no friction, do nothing*/
     1171                return;
     1172        }
     1173
     1174        if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
     1175
     1176        /*recover extra inputs from users, at current convergence iteration: */
     1177        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     1178
     1179        /*Build alpha2_list used by drag stiffness matrix*/
     1180        Friction* friction=NewFriction();
     1181       
     1182        /*Initialize all fields: */
     1183        friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
     1184        strcpy(friction->element_type,"2d");
     1185       
     1186        friction->gravity=matpar->GetG();
     1187        friction->rho_ice=matpar->GetRhoIce();
     1188        friction->rho_water=matpar->GetRhoWater();
     1189        friction->K=&this->properties.k[0];
     1190        friction->bed=&this->properties.b[0];
     1191        friction->thickness=&this->properties.h[0];
     1192        friction->velocities=&vxvy_list[0][0];
     1193        friction->p=this->properties.p;
     1194        friction->q=this->properties.q;
     1195
     1196        /*Compute alpha2_list: */
     1197        FrictionGetAlpha2(&alpha2_list[0],friction);
     1198
     1199        /*Erase friction object: */
     1200        DeleteFriction(&friction);
     1201
     1202        #ifdef _DEBUGELEMENTS_
     1203        if(my_rank==RANK && id==ELID){
     1204                printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
     1205        }
     1206        #endif
     1207
     1208        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1209        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1210
     1211        #ifdef _DEBUGELEMENTS_
     1212        if(my_rank==RANK && id==ELID){
     1213                printf("   gaussian points: \n");
     1214                for(i=0;i<num_gauss;i++){
     1215                        printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
     1216                }
     1217        }
     1218        #endif
     1219
     1220        /* Start  looping on the number of gaussian points: */
     1221        for (ig=0; ig<num_gauss; ig++){
     1222                /*Pick up the gaussian point: */
     1223                gauss_weight=*(gauss_weights+ig);
     1224                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1225                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1226                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1227
     1228
     1229                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
     1230                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
     1231                GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3);
     1232                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
     1233
     1234                if (slope_magnitude>MAXSLOPE){
     1235                        alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT);
     1236                        alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT);
     1237                        alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT);
     1238                }
     1239
     1240                /* Get Jacobian determinant: */
     1241                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1242
     1243                /*Get L matrix: */
     1244                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     1245
     1246                /*Now, take care of the basal friction if there is any: */
     1247                GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
     1248
     1249                DL_scalar=alpha2*gauss_weight*Jdet;
     1250                for (i=0;i<2;i++){
     1251                        DL[i][i]=DL_scalar;
     1252                }
     1253               
     1254                /*  Do the triple producte tL*D*L: */
     1255                TripleMultiply( &L[0][0],2,numdof,1,
     1256                                        &DL[0][0],2,2,0,
     1257                                        &L[0][0],2,numdof,0,
     1258                                        &Ke_gg_gaussian[0][0],0);
     1259
     1260                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     1261
     1262        } // for (ig=0; ig<num_gauss; ig++)
     1263
     1264        /*Add Ke_gg to global matrix Kgg: */
     1265        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1266
     1267        cleanup_and_return:
     1268        xfree((void**)&first_gauss_area_coord);
     1269        xfree((void**)&second_gauss_area_coord);
     1270        xfree((void**)&third_gauss_area_coord);
     1271        xfree((void**)&gauss_weights);
     1272
     1273}       
     1274/*}}}*/
     1275/*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/
     1276void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     1277
     1278        int i,j;
     1279
     1280        /* node data: */
     1281        const int    numgrids=3;
     1282        const int    NDOF1=1;
     1283        const int    numdof=NDOF1*numgrids;
     1284        double       xyz_list[numgrids][3];
     1285        int          doflist[numdof];
     1286        int          numberofdofspernode;
     1287
     1288        /* gaussian points: */
     1289        int     num_gauss,ig;
     1290        double* first_gauss_area_coord  =  NULL;
     1291        double* second_gauss_area_coord =  NULL;
     1292        double* third_gauss_area_coord  =  NULL;
     1293        double* gauss_weights           =  NULL;
     1294        double  gauss_weight;
     1295        double  gauss_l1l2l3[3];
     1296
     1297
     1298        /* surface normal: */
     1299        double x4,y4,z4;
     1300        double x5,y5,z5;
     1301        double x6,y6,z6;
     1302        double v46[3];
     1303        double v56[3];
     1304        double normal[3];
     1305        double norm_normal;
     1306        double nz;
     1307
     1308        /*Matrices: */
     1309        double DL_scalar;
     1310        double L[3];
     1311        double Jdet;
     1312
     1313        /* local element matrices: */
     1314        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1315        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
     1316
     1317        /*dynamic objects pointed to by hooks: */
     1318        Node**  nodes=NULL;
     1319        Matpar* matpar=NULL;
     1320        Matice* matice=NULL;
     1321        Numpar* numpar=NULL;
     1322
     1323        ParameterInputs* inputs=NULL;
     1324
     1325        /*recover pointers: */
     1326        inputs=(ParameterInputs*)vinputs;
     1327
     1328        /*recover objects from hooks: */
     1329        nodes=(Node**)hnodes.deliverp();
     1330        matpar=(Matpar*)hmatpar.delivers();
     1331        matice=(Matice*)hmatice.delivers();
     1332        numpar=(Numpar*)hnumpar.delivers();
     1333
     1334        /* Get node coordinates and dof list: */
     1335        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1336        GetDofList(&doflist[0],&numberofdofspernode);
     1337
     1338        /* Set Ke_gg to 0: */
     1339        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     1340
     1341        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1342        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1343
     1344        /*Build normal vector to the surface:*/
     1345
     1346        x4=xyz_list[0][0];
     1347        y4=xyz_list[0][1];
     1348        z4=xyz_list[0][2];
     1349
     1350        x5=xyz_list[1][0];
     1351        y5=xyz_list[1][1];
     1352        z5=xyz_list[1][2];
     1353
     1354        x6=xyz_list[2][0];
     1355        y6=xyz_list[2][1];
     1356        z6=xyz_list[2][2];
     1357
     1358        v46[0]=x4-x6;
     1359        v46[1]=y4-y6;
     1360        v46[2]=z4-z6;
     1361
     1362        v56[0]=x5-x6;
     1363        v56[1]=y5-y6;
     1364        v56[2]=z5-z6;
     1365
     1366        normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6);
     1367        normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6);
     1368        normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6);
     1369
     1370        norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
     1371        nz=1.0/norm_normal*normal[2];
     1372
     1373        /* Start  looping on the number of gaussian points: */
     1374        for (ig=0; ig<num_gauss; ig++){
     1375                /*Pick up the gaussian point: */
     1376                gauss_weight=*(gauss_weights+ig);
     1377                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1378                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1379                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1380
     1381                /* Get Jacobian determinant: */
     1382                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1383
     1384                //Get L matrix if viscous basal drag present:
     1385                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     1386
     1387                /**********************Do not forget the sign**********************************/
     1388                DL_scalar=- gauss_weight*Jdet*nz;
     1389                /******************************************************************************/
     1390
     1391                /*  Do the triple producte tL*D*L: */
     1392                TripleMultiply( L,1,3,1,
     1393                                        &DL_scalar,1,1,0,
     1394                                        L,1,3,0,
     1395                                        &Ke_gg_gaussian[0][0],0);
     1396
     1397                /* Add the Ke_gg_gaussian, onto Ke_gg: */
     1398                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     1399
     1400
     1401        } //for (ig=0; ig<num_gauss; ig++)
     1402
     1403        /*Add Ke_gg to global matrix Kgg: */
     1404        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1405
     1406cleanup_and_return:
     1407        xfree((void**)&first_gauss_area_coord);
     1408        xfree((void**)&second_gauss_area_coord);
     1409        xfree((void**)&third_gauss_area_coord);
     1410        xfree((void**)&gauss_weights);
     1411}
     1412/*}}}*/
     1413/*FUNCTION Tria::CreateKMatrixMelting {{{1*/
     1414void  Tria::CreateKMatrixMelting(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     1415
     1416        /*indexing: */
     1417        int i,j;
     1418
     1419        const int  numgrids=3;
     1420        const int  NDOF1=1;
     1421        const int  numdof=numgrids*NDOF1;
     1422        int        doflist[numdof];
     1423        int        numberofdofspernode;
     1424
     1425        /*Grid data: */
     1426        double     xyz_list[numgrids][3];
     1427
     1428        /*Material constants */
     1429        double     heatcapacity,latentheat;
     1430
     1431        /* gaussian points: */
     1432        int     num_area_gauss,ig;
     1433        double* gauss_weights  =  NULL;
     1434        double* first_gauss_area_coord  =  NULL;
     1435        double* second_gauss_area_coord =  NULL;
     1436        double* third_gauss_area_coord  =  NULL;
     1437        double  gauss_weight;
     1438        double  gauss_coord[3];
     1439
     1440        /*matrices: */
     1441        double     Jdet;
     1442        double     D_scalar;
     1443        double     K_terms[numdof][numdof]={0.0};
     1444        double     L[3];
     1445        double     tLD[3];
     1446        double     Ke_gaussian[numdof][numdof]={0.0};
     1447
     1448        /*dynamic objects pointed to by hooks: */
     1449        Node**  nodes=NULL;
     1450        Matpar* matpar=NULL;
     1451        Matice* matice=NULL;
     1452        Numpar* numpar=NULL;
     1453
     1454        /*recover objects from hooks: */
     1455        nodes=(Node**)hnodes.deliverp();
     1456        matpar=(Matpar*)hmatpar.delivers();
     1457        matice=(Matice*)hmatice.delivers();
     1458        numpar=(Numpar*)hnumpar.delivers();
     1459
     1460        /*Recover constants of ice */
     1461        latentheat=matpar->GetLatentHeat();
     1462        heatcapacity=matpar->GetHeatCapacity();
     1463
     1464        /* Get node coordinates and dof list: */
     1465        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1466        GetDofList(&doflist[0],&numberofdofspernode);
     1467
     1468        /* Get gaussian points and weights: */
     1469        GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1470
     1471        /* Start looping on the number of gauss  (nodes on the bedrock) */
     1472        for (ig=0; ig<num_area_gauss; ig++){
     1473                gauss_weight=*(gauss_weights+ig);
     1474                gauss_coord[0]=*(first_gauss_area_coord+ig);
     1475                gauss_coord[1]=*(second_gauss_area_coord+ig);
     1476                gauss_coord[2]=*(third_gauss_area_coord+ig);
     1477
     1478                //Get the Jacobian determinant
     1479                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
     1480
     1481                /*Get L matrix : */
     1482                GetL(&L[0], &xyz_list[0][0], gauss_coord,NDOF1);
     1483
     1484                /*Calculate DL on gauss point */
     1485                D_scalar=latentheat/heatcapacity*gauss_weight*Jdet;
     1486
     1487                /*  Do the triple product tL*D*L: */
     1488                MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
     1489                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
     1490
     1491                for(i=0;i<numgrids;i++){
     1492                        for(j=0;j<numgrids;j++){
     1493                                K_terms[i][j]+=Ke_gaussian[i][j];
     1494                        }
     1495                }
     1496        }
     1497
     1498        /*Add Ke_gg to global matrix Kgg: */
     1499        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
     1500
     1501cleanup_and_return:
     1502        xfree((void**)&first_gauss_area_coord);
     1503        xfree((void**)&second_gauss_area_coord);
     1504        xfree((void**)&third_gauss_area_coord);
     1505        xfree((void**)&gauss_weights);
     1506
     1507}
     1508/*}}}*/
     1509/*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/
     1510void  Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     1511
     1512        /* local declarations */
     1513        int             i,j;
     1514
     1515        /* node data: */
     1516        const int    numgrids=3;
     1517        const int    NDOF1=1;
     1518        const int    numdof=NDOF1*numgrids;
     1519        double       xyz_list[numgrids][3];
     1520        int          doflist[numdof];
     1521        int          numberofdofspernode;
     1522
     1523        /* gaussian points: */
     1524        int     num_gauss,ig;
     1525        double* first_gauss_area_coord  =  NULL;
     1526        double* second_gauss_area_coord =  NULL;
     1527        double* third_gauss_area_coord  =  NULL;
     1528        double* gauss_weights           =  NULL;
     1529        double  gauss_weight;
     1530        double  gauss_l1l2l3[3];
     1531
    8471532        /* matrices: */
    8481533        double L[numgrids];
    … …  
    8521537        double DLprime[2][2]={0.0};
    8531538        double DL_scalar;
    854         double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
    855         double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    856         double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    857         double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     1539        double Ke_gg[numdof][numdof]={0.0};
     1540        double Ke_gg_gaussian[numdof][numdof]={0.0};
     1541        double Ke_gg_thickness1[numdof][numdof]={0.0};
     1542        double Ke_gg_thickness2[numdof][numdof]={0.0};
    8581543        double Jdettria;
    8591544
    8601545        /*input parameters for structural analysis (diagnostic): */
    861         double  surface_normal[3];
    862         double  nx,ny,norm;
    8631546        double  vxvy_list[numgrids][2]={0.0};
    8641547        double  vx_list[numgrids]={0.0};
    … …  
    8711554        double  K[2][2]={0.0};
    8721555        double  KDL[2][2]={0.0};
     1556        double  dt;
    8731557        int     dofs[2]={0,1};
    874         int     found=0;
     1558        int     found;
    8751559
    8761560        /*dynamic objects pointed to by hooks: */
    … …  
    9001584        }
    9011585
     1586        found=inputs->Recover("dt",&dt);
     1587        if(!found)ISSMERROR(" could not find dt in inputs!");
     1588
    9021589        /* Get node coordinates and dof list: */
    9031590        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    9041591        GetDofList(&doflist[0],&numberofdofspernode);
    9051592
    906         /*Modify z so that it reflects the surface*/
    907         for(i=0;i<numgrids;i++) xyz_list[i][2]=this->properties.s[i];
    908 
    909         /*Get normal vector to the surface*/
    910         nx=(vx_list[0]+vx_list[1]+vx_list[2])/3;
    911         ny=(vy_list[0]+vy_list[1]+vy_list[2])/3;
    912         if(nx==0 && ny==0){
    913                 SurfaceNormal(&surface_normal[0],xyz_list);
    914                 nx=surface_normal[0];
    915                 ny=surface_normal[1];
    916         }
    917         if(nx==0 && ny==0){
    918                 nx=0;
    919                 ny=1;
    920         }
    921         norm=pow( pow(nx,2)+pow(ny,2) , (double).5);
    922         nx=nx/norm;
    923         ny=ny/norm;
    924 
    9251593        //Create Artificial diffusivity once for all if requested
    926         if(numpar->artdiff){
     1594        if(numpar->artdiff==1){
    9271595                //Get the Jacobian determinant
    9281596                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    … …  
    9331601                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    9341602
    935                 K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
    936                 K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     1603                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
     1604                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
    9371605        }
    9381606
    … …  
    9421610        /* Start  looping on the number of gaussian points: */
    9431611        for (ig=0; ig<num_gauss; ig++){
     1612
    9441613                /*Pick up the gaussian point: */
    9451614                gauss_weight=*(gauss_weights+ig);
    … …  
    9511620                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    9521621
     1622                /*Get L matrix: */
     1623                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     1624
     1625                DL_scalar=gauss_weight*Jdettria;
     1626
     1627                /*  Do the triple product tL*D*L: */
     1628                TripleMultiply( &L[0],1,numdof,1,
     1629                                        &DL_scalar,1,1,0,
     1630                                        &L[0],1,numdof,0,
     1631                                        &Ke_gg_gaussian[0][0],0);
     1632
    9531633                /*Get B  and B prime matrix: */
    9541634                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    … …  
    9651645                dvydy=dvy[1];
    9661646
    967                 DL_scalar=gauss_weight*Jdettria;
    968 
    969                 DLprime[0][0]=DL_scalar*nx;
    970                 DLprime[1][1]=DL_scalar*ny;
     1647                DL_scalar=dt*gauss_weight*Jdettria;
     1648
     1649                //Create DL and DLprime matrix
     1650                DL[0][0]=DL_scalar*dvxdx;
     1651                DL[1][1]=DL_scalar*dvydy;
     1652
     1653                DLprime[0][0]=DL_scalar*vx;
     1654                DLprime[1][1]=DL_scalar*vy;
    9711655
    9721656                //Do the triple product tL*D*L.
    973                 //Ke_gg_velocities=B'*DLprime*Bprime;
     1657                //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
     1658
     1659                TripleMultiply( &B[0][0],2,numdof,1,
     1660                                        &DL[0][0],2,2,0,
     1661                                        &B[0][0],2,numdof,0,
     1662                                        &Ke_gg_thickness1[0][0],0);
     1663
    9741664                TripleMultiply( &B[0][0],2,numdof,1,
    9751665                                        &DLprime[0][0],2,2,0,
    9761666                                        &Bprime[0][0],2,numdof,0,
    977                                         &Ke_gg_velocities2[0][0],0);
     1667                                        &Ke_gg_thickness2[0][0],0);
    9781668
    9791669                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    980                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j];
    981 
    982                 if(numpar->artdiff){
     1670                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     1671                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
     1672                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
     1673
     1674                if(numpar->artdiff==1){
    9831675
    9841676                        /* Compute artificial diffusivity */
    … …  
    10441736}
    10451737/*}}}*/
    1046 /*FUNCTION Tria::CreateKMatrixDiagnosticHoriz {{{1*/
    1047 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1048 
    1049         /* local declarations */
    1050         int             i,j;
    1051 
    1052         /* node data: */
    1053         const int    numgrids=3;
    1054         const int    numdof=2*numgrids;
    1055         double       xyz_list[numgrids][3];
    1056         int          doflist[numdof];
    1057         int          numberofdofspernode;
    1058 
    1059         /* gaussian points: */
    1060         int     num_gauss,ig;
    1061         double* first_gauss_area_coord  =  NULL;
    1062         double* second_gauss_area_coord =  NULL;
    1063         double* third_gauss_area_coord  =  NULL;
    1064         double* gauss_weights           =  NULL;
    1065         double  gauss_weight;
    1066         double  gauss_l1l2l3[3];
    1067 
    1068         /* material data: */
    1069         double viscosity; //viscosity
    1070         double newviscosity; //viscosity
    1071         double oldviscosity; //viscosity
    1072 
    1073         /* strain rate: */
    1074         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    1075         double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/
    1076 
    1077         /* matrices: */
    1078         double B[3][numdof];
    1079         double Bprime[3][numdof];
    1080         double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
    1081         double D_scalar;
    1082 
    1083         /* local element matrices: */
    1084         double Ke_gg[numdof][numdof]; //local element stiffness matrix
    1085         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    1086 
    1087         double Jdet;
    1088 
    1089         /*input parameters for structural analysis (diagnostic): */
    1090         double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    1091         double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    1092         double  thickness;
    1093         int     dofs[2]={0,1};
    1094 
    1095         /*dynamic objects pointed to by hooks: */
    1096         Node**  nodes=NULL;
    1097         Matpar* matpar=NULL;
    1098         Matice* matice=NULL;
    1099         Numpar* numpar=NULL;
    1100 
    1101         ParameterInputs* inputs=NULL;
    1102 
    1103         /*First, if we are on water, return empty matrix: */
    1104         if(this->properties.onwater) return;
    1105 
    1106         /*recover pointers: */
    1107         inputs=(ParameterInputs*)vinputs;
    1108 
    1109         /*recover objects from hooks: */
    1110         nodes =(Node**) hnodes.deliverp();
    1111         matpar=(Matpar*)hmatpar.delivers();
    1112         matice=(Matice*)hmatice.delivers();
    1113         numpar=(Numpar*)hnumpar.delivers();
    1114 
    1115         /*recover extra inputs from users, at current convergence iteration: */
    1116         inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    1117         inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    1118 
    1119         /* Get node coordinates and dof list: */
    1120         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    1121         GetDofList(&doflist[0],&numberofdofspernode);
    1122 
    1123         /* Set Ke_gg to 0: */
    1124         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    1125 
    1126         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1127         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1128 
    1129         /* Start  looping on the number of gaussian points: */
    1130         for (ig=0; ig<num_gauss; ig++){
    1131                 /*Pick up the gaussian point: */
    1132                 gauss_weight=*(gauss_weights+ig);
    1133                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1134                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1135                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1136 
    1137 
    1138                 /*Compute thickness at gaussian point: */
    1139                 GetParameterValue(&thickness, &this->properties.h[0],gauss_l1l2l3);
    1140 
    1141                 /*Get strain rate from velocity: */
    1142                 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
    1143                 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
    1144 
    1145                 /*Get viscosity: */
    1146                 matice->GetViscosity2d(&viscosity, &epsilon[0]);
    1147                 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    1148 
    1149                 /* Get Jacobian determinant: */
    1150                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1151 
    1152                 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    1153                         onto this scalar matrix, so that we win some computational time: */
    1154                 newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity);
    1155                 D_scalar=newviscosity*thickness*gauss_weight*Jdet;
    1156 
    1157                 for (i=0;i<3;i++){
    1158                         D[i][i]=D_scalar;
    1159                 }
    1160 
    1161                 /*Get B and Bprime matrices: */
    1162                 GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    1163                 GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
    1164 
    1165                 /*  Do the triple product tB*D*Bprime: */
    1166                 TripleMultiply( &B[0][0],3,numdof,1,
    1167                                         &D[0][0],3,3,0,
    1168                                         &Bprime[0][0],3,numdof,0,
    1169                                         &Ke_gg_gaussian[0][0],0);
    1170 
    1171                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    1172                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1173 
    1174         } // for (ig=0; ig<num_gauss; ig++)
    1175 
    1176         /*Add Ke_gg to global matrix Kgg: */
    1177         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1178 
    1179         /*Do not forget to include friction: */
    1180         if(!this->properties.shelf){
    1181                 CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
    1182         }
    1183 
    1184 cleanup_and_return:
    1185         xfree((void**)&first_gauss_area_coord);
    1186         xfree((void**)&second_gauss_area_coord);
    1187         xfree((void**)&third_gauss_area_coord);
    1188         xfree((void**)&gauss_weights);
    1189 
    1190 }
    1191 /*}}}*/
    1192 /*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/
    1193 void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1194 
    1195 
    1196         /* local declarations */
    1197         int             i,j;
    1198 
    1199         /* node data: */
    1200         const int    numgrids=3;
    1201         const int    numdof=2*numgrids;
    1202         double       xyz_list[numgrids][3];
    1203         int          doflist[numdof];
    1204         int          numberofdofspernode;
    1205        
    1206         /* gaussian points: */
    1207         int     num_gauss,ig;
    1208         double* first_gauss_area_coord  =  NULL;
    1209         double* second_gauss_area_coord =  NULL;
    1210         double* third_gauss_area_coord  =  NULL;
    1211         double* gauss_weights           =  NULL;
    1212         double  gauss_weight;
    1213         double  gauss_l1l2l3[3];
    1214 
    1215         /* matrices: */
    1216         double L[2][numdof];
    1217         double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
    1218         double DL_scalar;
    1219 
    1220         /* local element matrices: */
    1221         double Ke_gg[numdof][numdof]; //local element stiffness matrix
    1222         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    1223        
    1224         double Jdet;
    1225        
    1226         /*slope: */
    1227         double  slope[2]={0.0,0.0};
    1228         double  slope_magnitude;
    1229 
    1230         /*input parameters for structural analysis (diagnostic): */
    1231         double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    1232         int     dofs[2]={0,1};
    1233 
    1234         /*friction: */
    1235         double alpha2_list[numgrids]={0.0,0.0,0.0};
    1236         double alpha2;
    1237 
    1238         double MAXSLOPE=.06; // 6 %
    1239         double MOUNTAINKEXPONENT=10;
    1240 
    1241         /*dynamic objects pointed to by hooks: */
    1242         Node**  nodes=NULL;
    1243         Matpar* matpar=NULL;
    1244         Matice* matice=NULL;
    1245         Numpar* numpar=NULL;
    1246 
    1247         ParameterInputs* inputs=NULL;
    1248 
    1249         /*recover pointers: */
    1250         inputs=(ParameterInputs*)vinputs;
    1251 
    1252         /*recover objects from hooks: */
    1253         nodes=(Node**)hnodes.deliverp();
    1254         matpar=(Matpar*)hmatpar.delivers();
    1255         matice=(Matice*)hmatice.delivers();
    1256         numpar=(Numpar*)hnumpar.delivers();
    1257        
    1258         /* Get node coordinates and dof list: */
    1259         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    1260         GetDofList(&doflist[0],&numberofdofspernode);
    1261 
    1262         /* Set Ke_gg to 0: */
    1263         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    1264 
    1265         if (this->properties.shelf){
    1266                 /*no friction, do nothing*/
    1267                 return;
    1268         }
    1269 
    1270         if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
    1271 
    1272         /*recover extra inputs from users, at current convergence iteration: */
    1273         inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    1274 
    1275         /*Build alpha2_list used by drag stiffness matrix*/
    1276         Friction* friction=NewFriction();
    1277        
    1278         /*Initialize all fields: */
    1279         friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    1280         strcpy(friction->element_type,"2d");
    1281        
    1282         friction->gravity=matpar->GetG();
    1283         friction->rho_ice=matpar->GetRhoIce();
    1284         friction->rho_water=matpar->GetRhoWater();
    1285         friction->K=&this->properties.k[0];
    1286         friction->bed=&this->properties.b[0];
    1287         friction->thickness=&this->properties.h[0];
    1288         friction->velocities=&vxvy_list[0][0];
    1289         friction->p=this->properties.p;
    1290         friction->q=this->properties.q;
    1291 
    1292         /*Compute alpha2_list: */
    1293         FrictionGetAlpha2(&alpha2_list[0],friction);
    1294 
    1295         /*Erase friction object: */
    1296         DeleteFriction(&friction);
    1297 
    1298         #ifdef _DEBUGELEMENTS_
    1299         if(my_rank==RANK && id==ELID){
    1300                 printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
    1301         }
    1302         #endif
    1303 
    1304         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1305         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1306 
    1307         #ifdef _DEBUGELEMENTS_
    1308         if(my_rank==RANK && id==ELID){
    1309                 printf("   gaussian points: \n");
    1310                 for(i=0;i<num_gauss;i++){
    1311                         printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
    1312                 }
    1313         }
    1314         #endif
    1315 
    1316         /* Start  looping on the number of gaussian points: */
    1317         for (ig=0; ig<num_gauss; ig++){
    1318                 /*Pick up the gaussian point: */
    1319                 gauss_weight=*(gauss_weights+ig);
    1320                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1321                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1322                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1323 
    1324 
    1325                 // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    1326                 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    1327                 GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3);
    1328                 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    1329 
    1330                 if (slope_magnitude>MAXSLOPE){
    1331                         alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT);
    1332                         alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT);
    1333                         alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT);
    1334                 }
    1335 
    1336                 /* Get Jacobian determinant: */
    1337                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1338 
    1339                 /*Get L matrix: */
    1340                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
    1341 
    1342                 /*Now, take care of the basal friction if there is any: */
    1343                 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
    1344 
    1345                 DL_scalar=alpha2*gauss_weight*Jdet;
    1346                 for (i=0;i<2;i++){
    1347                         DL[i][i]=DL_scalar;
    1348                 }
    1349                
    1350                 /*  Do the triple producte tL*D*L: */
    1351                 TripleMultiply( &L[0][0],2,numdof,1,
    1352                                         &DL[0][0],2,2,0,
    1353                                         &L[0][0],2,numdof,0,
    1354                                         &Ke_gg_gaussian[0][0],0);
    1355 
    1356                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1357 
    1358         } // for (ig=0; ig<num_gauss; ig++)
    1359 
    1360         /*Add Ke_gg to global matrix Kgg: */
    1361         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1362 
    1363         cleanup_and_return:
    1364         xfree((void**)&first_gauss_area_coord);
    1365         xfree((void**)&second_gauss_area_coord);
    1366         xfree((void**)&third_gauss_area_coord);
    1367         xfree((void**)&gauss_weights);
    1368 
    1369 }       
    1370 /*}}}*/
    1371 /*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/
    1372 void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1373 
    1374         int i,j;
    1375 
    1376         /* node data: */
    1377         const int    numgrids=3;
    1378         const int    NDOF1=1;
    1379         const int    numdof=NDOF1*numgrids;
    1380         double       xyz_list[numgrids][3];
    1381         int          doflist[numdof];
    1382         int          numberofdofspernode;
    1383 
    1384         /* gaussian points: */
    1385         int     num_gauss,ig;
    1386         double* first_gauss_area_coord  =  NULL;
    1387         double* second_gauss_area_coord =  NULL;
    1388         double* third_gauss_area_coord  =  NULL;
    1389         double* gauss_weights           =  NULL;
    1390         double  gauss_weight;
    1391         double  gauss_l1l2l3[3];
    1392 
    1393 
    1394         /* surface normal: */
    1395         double x4,y4,z4;
    1396         double x5,y5,z5;
    1397         double x6,y6,z6;
    1398         double v46[3];
    1399         double v56[3];
    1400         double normal[3];
    1401         double norm_normal;
    1402         double nz;
    1403 
    1404         /*Matrices: */
    1405         double DL_scalar;
    1406         double L[3];
    1407         double Jdet;
    1408 
    1409         /* local element matrices: */
    1410         double Ke_gg[numdof][numdof]; //local element stiffness matrix
    1411         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    1412 
    1413         /*dynamic objects pointed to by hooks: */
    1414         Node**  nodes=NULL;
    1415         Matpar* matpar=NULL;
    1416         Matice* matice=NULL;
    1417         Numpar* numpar=NULL;
    1418 
    1419         ParameterInputs* inputs=NULL;
    1420 
    1421         /*recover pointers: */
    1422         inputs=(ParameterInputs*)vinputs;
    1423 
    1424         /*recover objects from hooks: */
    1425         nodes=(Node**)hnodes.deliverp();
    1426         matpar=(Matpar*)hmatpar.delivers();
    1427         matice=(Matice*)hmatice.delivers();
    1428         numpar=(Numpar*)hnumpar.delivers();
    1429 
    1430         /* Get node coordinates and dof list: */
    1431         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    1432         GetDofList(&doflist[0],&numberofdofspernode);
    1433 
    1434         /* Set Ke_gg to 0: */
    1435         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    1436 
    1437         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1438         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1439 
    1440         /*Build normal vector to the surface:*/
    1441 
    1442         x4=xyz_list[0][0];
    1443         y4=xyz_list[0][1];
    1444         z4=xyz_list[0][2];
    1445 
    1446         x5=xyz_list[1][0];
    1447         y5=xyz_list[1][1];
    1448         z5=xyz_list[1][2];
    1449 
    1450         x6=xyz_list[2][0];
    1451         y6=xyz_list[2][1];
    1452         z6=xyz_list[2][2];
    1453 
    1454         v46[0]=x4-x6;
    1455         v46[1]=y4-y6;
    1456         v46[2]=z4-z6;
    1457 
    1458         v56[0]=x5-x6;
    1459         v56[1]=y5-y6;
    1460         v56[2]=z5-z6;
    1461 
    1462         normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6);
    1463         normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6);
    1464         normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6);
    1465 
    1466         norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
    1467         nz=1.0/norm_normal*normal[2];
    1468 
    1469         /* Start  looping on the number of gaussian points: */
    1470         for (ig=0; ig<num_gauss; ig++){
    1471                 /*Pick up the gaussian point: */
    1472                 gauss_weight=*(gauss_weights+ig);
    1473                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1474                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1475                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1476 
    1477                 /* Get Jacobian determinant: */
    1478                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1479 
    1480                 //Get L matrix if viscous basal drag present:
    1481                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
    1482 
    1483                 /**********************Do not forget the sign**********************************/
    1484                 DL_scalar=- gauss_weight*Jdet*nz;
    1485                 /******************************************************************************/
    1486 
    1487                 /*  Do the triple producte tL*D*L: */
    1488                 TripleMultiply( L,1,3,1,
    1489                                         &DL_scalar,1,1,0,
    1490                                         L,1,3,0,
    1491                                         &Ke_gg_gaussian[0][0],0);
    1492 
    1493                 /* Add the Ke_gg_gaussian, onto Ke_gg: */
    1494                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1495 
    1496 
    1497         } //for (ig=0; ig<num_gauss; ig++)
    1498 
    1499         /*Add Ke_gg to global matrix Kgg: */
    1500         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1501 
    1502 cleanup_and_return:
    1503         xfree((void**)&first_gauss_area_coord);
    1504         xfree((void**)&second_gauss_area_coord);
    1505         xfree((void**)&third_gauss_area_coord);
    1506         xfree((void**)&gauss_weights);
    1507 }
    1508 /*}}}*/
    1509 /*FUNCTION Tria::CreateKMatrixMelting {{{1*/
    1510 void  Tria::CreateKMatrixMelting(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1511 
    1512         /*indexing: */
    1513         int i,j;
    1514 
    1515         const int  numgrids=3;
    1516         const int  NDOF1=1;
    1517         const int  numdof=numgrids*NDOF1;
    1518         int        doflist[numdof];
    1519         int        numberofdofspernode;
    1520 
    1521         /*Grid data: */
    1522         double     xyz_list[numgrids][3];
    1523 
    1524         /*Material constants */
    1525         double     heatcapacity,latentheat;
    1526 
    1527         /* gaussian points: */
    1528         int     num_area_gauss,ig;
    1529         double* gauss_weights  =  NULL;
    1530         double* first_gauss_area_coord  =  NULL;
    1531         double* second_gauss_area_coord =  NULL;
    1532         double* third_gauss_area_coord  =  NULL;
    1533         double  gauss_weight;
    1534         double  gauss_coord[3];
    1535 
    1536         /*matrices: */
    1537         double     Jdet;
    1538         double     D_scalar;
    1539         double     K_terms[numdof][numdof]={0.0};
    1540         double     L[3];
    1541         double     tLD[3];
    1542         double     Ke_gaussian[numdof][numdof]={0.0};
    1543 
    1544         /*dynamic objects pointed to by hooks: */
    1545         Node**  nodes=NULL;
    1546         Matpar* matpar=NULL;
    1547         Matice* matice=NULL;
    1548         Numpar* numpar=NULL;
    1549 
    1550         /*recover objects from hooks: */
    1551         nodes=(Node**)hnodes.deliverp();
    1552         matpar=(Matpar*)hmatpar.delivers();
    1553         matice=(Matice*)hmatice.delivers();
    1554         numpar=(Numpar*)hnumpar.delivers();
    1555 
    1556         /*Recover constants of ice */
    1557         latentheat=matpar->GetLatentHeat();
    1558         heatcapacity=matpar->GetHeatCapacity();
    1559 
    1560         /* Get node coordinates and dof list: */
    1561         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    1562         GetDofList(&doflist[0],&numberofdofspernode);
    1563 
    1564         /* Get gaussian points and weights: */
    1565         GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1566 
    1567         /* Start looping on the number of gauss  (nodes on the bedrock) */
    1568         for (ig=0; ig<num_area_gauss; ig++){
    1569                 gauss_weight=*(gauss_weights+ig);
    1570                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    1571                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    1572                 gauss_coord[2]=*(third_gauss_area_coord+ig);
    1573 
    1574                 //Get the Jacobian determinant
    1575                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
    1576 
    1577                 /*Get L matrix : */
    1578                 GetL(&L[0], &xyz_list[0][0], gauss_coord,NDOF1);
    1579 
    1580                 /*Calculate DL on gauss point */
    1581                 D_scalar=latentheat/heatcapacity*gauss_weight*Jdet;
    1582 
    1583                 /*  Do the triple product tL*D*L: */
    1584                 MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
    1585                 MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
    1586 
    1587                 for(i=0;i<numgrids;i++){
    1588                         for(j=0;j<numgrids;j++){
    1589                                 K_terms[i][j]+=Ke_gaussian[i][j];
    1590                         }
    1591                 }
    1592         }
    1593 
    1594         /*Add Ke_gg to global matrix Kgg: */
    1595         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    1596 
    1597 cleanup_and_return:
    1598         xfree((void**)&first_gauss_area_coord);
    1599         xfree((void**)&second_gauss_area_coord);
    1600         xfree((void**)&third_gauss_area_coord);
    1601         xfree((void**)&gauss_weights);
    1602 
    1603 }
    1604 /*}}}*/
    1605 /*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/
    1606 void  Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1607 
    1608         /* local declarations */
    1609         int             i,j;
    1610 
    1611         /* node data: */
    1612         const int    numgrids=3;
    1613         const int    NDOF1=1;
    1614         const int    numdof=NDOF1*numgrids;
    1615         double       xyz_list[numgrids][3];
    1616         int          doflist[numdof];
    1617         int          numberofdofspernode;
    1618 
    1619         /* gaussian points: */
    1620         int     num_gauss,ig;
    1621         double* first_gauss_area_coord  =  NULL;
    1622         double* second_gauss_area_coord =  NULL;
    1623         double* third_gauss_area_coord  =  NULL;
    1624         double* gauss_weights           =  NULL;
    1625         double  gauss_weight;
    1626         double  gauss_l1l2l3[3];
    1627 
    1628         /* matrices: */
    1629         double L[numgrids];
    1630         double B[2][numgrids];
    1631         double Bprime[2][numgrids];
    1632         double DL[2][2]={0.0};
    1633         double DLprime[2][2]={0.0};
    1634         double DL_scalar;
    1635         double Ke_gg[numdof][numdof]={0.0};
    1636         double Ke_gg_gaussian[numdof][numdof]={0.0};
    1637         double Ke_gg_thickness1[numdof][numdof]={0.0};
    1638         double Ke_gg_thickness2[numdof][numdof]={0.0};
    1639         double Jdettria;
    1640 
    1641         /*input parameters for structural analysis (diagnostic): */
    1642         double  vxvy_list[numgrids][2]={0.0};
    1643         double  vx_list[numgrids]={0.0};
    1644         double  vy_list[numgrids]={0.0};
    1645         double  dvx[2];
    1646         double  dvy[2];
    1647         double  vx,vy;
    1648         double  dvxdx,dvydy;
    1649         double  v_gauss[2]={0.0};
    1650         double  K[2][2]={0.0};
    1651         double  KDL[2][2]={0.0};
    1652         double  dt;
    1653         int     dofs[2]={0,1};
    1654         int     found;
    1655 
    1656         /*dynamic objects pointed to by hooks: */
    1657         Node**  nodes=NULL;
    1658         Matpar* matpar=NULL;
    1659         Matice* matice=NULL;
    1660         Numpar* numpar=NULL;
    1661 
    1662         ParameterInputs* inputs=NULL;
    1663 
    1664         /*recover pointers: */
    1665         inputs=(ParameterInputs*)vinputs;
    1666 
    1667         /*recover objects from hooks: */
    1668         nodes=(Node**)hnodes.deliverp();
    1669         matpar=(Matpar*)hmatpar.delivers();
    1670         matice=(Matice*)hmatice.delivers();
    1671         numpar=(Numpar*)hnumpar.delivers();
    1672 
    1673         /*recover extra inputs from users, at current convergence iteration: */
    1674         found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    1675         if(!found)ISSMERROR(" could not find velocity_average  in inputs!");
    1676 
    1677         for(i=0;i<numgrids;i++){
    1678                 vx_list[i]=vxvy_list[i][0];
    1679                 vy_list[i]=vxvy_list[i][1];
    1680         }
    1681 
    1682         found=inputs->Recover("dt",&dt);
    1683         if(!found)ISSMERROR(" could not find dt in inputs!");
    1684 
    1685         /* Get node coordinates and dof list: */
    1686         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    1687         GetDofList(&doflist[0],&numberofdofspernode);
    1688 
    1689         //Create Artificial diffusivity once for all if requested
    1690         if(numpar->artdiff==1){
    1691                 //Get the Jacobian determinant
    1692                 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    1693                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    1694 
    1695                 //Build K matrix (artificial diffusivity matrix)
    1696                 v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
    1697                 v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    1698 
    1699                 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
    1700                 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
    1701         }
    1702 
    1703         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1704         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1705 
    1706         /* Start  looping on the number of gaussian points: */
    1707         for (ig=0; ig<num_gauss; ig++){
    1708 
    1709                 /*Pick up the gaussian point: */
    1710                 gauss_weight=*(gauss_weights+ig);
    1711                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1712                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1713                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1714 
    1715                 /* Get Jacobian determinant: */
    1716                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    1717 
    1718                 /*Get L matrix: */
    1719                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
    1720 
    1721                 DL_scalar=gauss_weight*Jdettria;
    1722 
    1723                 /*  Do the triple product tL*D*L: */
    1724                 TripleMultiply( &L[0],1,numdof,1,
    1725                                         &DL_scalar,1,1,0,
    1726                                         &L[0],1,numdof,0,
    1727                                         &Ke_gg_gaussian[0][0],0);
    1728 
    1729                 /*Get B  and B prime matrix: */
    1730                 GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    1731                 GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
    1732 
    1733                 //Get vx, vy and their derivatives at gauss point
    1734                 GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
    1735                 GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
    1736 
    1737                 GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
    1738                 GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
    1739 
    1740                 dvxdx=dvx[0];
    1741                 dvydy=dvy[1];
    1742 
    1743                 DL_scalar=dt*gauss_weight*Jdettria;
    1744 
    1745                 //Create DL and DLprime matrix
    1746                 DL[0][0]=DL_scalar*dvxdx;
    1747                 DL[1][1]=DL_scalar*dvydy;
    1748 
    1749                 DLprime[0][0]=DL_scalar*vx;
    1750                 DLprime[1][1]=DL_scalar*vy;
    1751 
    1752                 //Do the triple product tL*D*L.
    1753                 //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
    1754 
    1755                 TripleMultiply( &B[0][0],2,numdof,1,
    1756                                         &DL[0][0],2,2,0,
    1757                                         &B[0][0],2,numdof,0,
    1758                                         &Ke_gg_thickness1[0][0],0);
    1759 
    1760                 TripleMultiply( &B[0][0],2,numdof,1,
    1761                                         &DLprime[0][0],2,2,0,
    1762                                         &Bprime[0][0],2,numdof,0,
    1763                                         &Ke_gg_thickness2[0][0],0);
    1764 
    1765                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    1766                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1767                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
    1768                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
    1769 
    1770                 if(numpar->artdiff==1){
    1771 
    1772                         /* Compute artificial diffusivity */
    1773                         KDL[0][0]=DL_scalar*K[0][0];
    1774                         KDL[1][1]=DL_scalar*K[1][1];
    1775 
    1776                         TripleMultiply( &Bprime[0][0],2,numdof,1,
    1777                                                 &KDL[0][0],2,2,0,
    1778                                                 &Bprime[0][0],2,numdof,0,
    1779                                                 &Ke_gg_gaussian[0][0],0);
    1780 
    1781                         /* Add artificial diffusivity matrix */
    1782                         for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1783 
    1784                 }
    1785 
    1786 #ifdef _DEBUGELEMENTS_
    1787                 if(my_rank==RANK && id==ELID){
    1788                         printf("      B:\n");
    1789                         for(i=0;i<3;i++){
    1790                                 for(j=0;j<numdof;j++){
    1791                                         printf("%g ",B[i][j]);
    1792                                 }
    1793                                 printf("\n");
    1794                         }
    1795                         printf("      Bprime:\n");
    1796                         for(i=0;i<3;i++){
    1797                                 for(j=0;j<numdof;j++){
    1798                                         printf("%g ",Bprime[i][j]);
    1799                                 }
    1800                                 printf("\n");
    1801                         }
    1802                 }
    1803 #endif
    1804         } // for (ig=0; ig<num_gauss; ig++)
    1805 
    1806         /*Add Ke_gg to global matrix Kgg: */
    1807         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1808 
    1809 #ifdef _DEBUGELEMENTS_
    1810         if(my_rank==RANK && id==ELID){
    1811                 printf("      Ke_gg erms:\n");
    1812                 for( i=0; i<numdof; i++){
    1813                         for (j=0;j<numdof;j++){
    1814                                 printf("%g ",Ke_gg[i][j]);
    1815                         }
    1816                         printf("\n");
    1817                 }
    1818                 printf("      Ke_gg row_indices:\n");
    1819                 for( i=0; i<numdof; i++){
    1820                         printf("%i ",doflist[i]);
    1821                 }
    1822 
    1823         }
    1824 #endif
    1825 
    1826 cleanup_and_return:
    1827         xfree((void**)&first_gauss_area_coord);
    1828         xfree((void**)&second_gauss_area_coord);
    1829         xfree((void**)&third_gauss_area_coord);
    1830         xfree((void**)&gauss_weights);
    1831 
    1832 }
    1833 /*}}}*/
    18341738/*FUNCTION Tria::CreateKMatrixPrognostic2 {{{1*/
    18351739void  Tria::CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    … …  
    18871791
    18881792        /*recover objects from hooks: */
    1889         nodes=(Node**)hnodes.deliverp();
     1793        nodes=(Node**)   hnodes.deliverp();
    18901794        matpar=(Matpar*)hmatpar.delivers();
    18911795        matice=(Matice*)hmatice.delivers();
    … …  
    22492153        double  melting_list[numgrids]={0.0};
    22502154        double  melting_g;
    2251         double  thickness_list[numgrids]={0.0};
    2252         double  thickness_g;
     2155        double  dhdt_list[numgrids]={0.0};
     2156        double  dhdt_g;
    22532157        int     dofs[1]={0};
    22542158        int     found=0;
    … …  
    22662170
    22672171        /*recover objects from hooks: */
    2268         nodes=(Node**)hnodes.deliverp();
     2172        nodes=(Node**)   hnodes.deliverp();
    22692173        matpar=(Matpar*)hmatpar.delivers();
    22702174        matice=(Matice*)hmatice.delivers();
    … …  
    22762180        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
    22772181        if(!found)ISSMERROR(" could not find melting in inputs!");
     2182        found=inputs->Recover("dhdt",&dhdt_list[0],1,dofs,numgrids,(void**)nodes);
     2183        if(!found)ISSMERROR(" could not find dhdt in inputs!");
    22782184
    22792185        /* Get node coordinates and dof list: */
    … …  
    23012207                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
    23022208                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
     2209                GetParameterValue(&dhdt_g, &dhdt_list[0],gauss_l1l2l3);
    23032210
    23042211                /* Add value into pe_g: */
    2305                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
     2212                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g+dhdt_g)*L[i];
    23062213
    23072214        } // for (ig=0; ig<num_gauss; ig++)

  • issm/trunk/src/c/parallel/balancedthickness.cpp

    r3567 r3582  
    2525        Model* model=NULL;
    2626
    27         Vec     u_g=NULL;
    28         double* u_g_serial=NULL;
    29         double* melting_g=NULL;
    30         double* accumulation_g=NULL;
     27        double* vx_g=NULL;
     28        double* vy_g=NULL;
     29        double* m_g=NULL;
     30        double* a_g=NULL;
     31        double* h_g=NULL;
     32        double* dhdt_g=NULL;
    3133        double  dt;
    3234        double  yts;
    … …  
    6062        MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
    6163
    62         _printf_("recover , input file name and output file name:\n");
     64        _printf_("recover input file name and output file name:\n");
    6365        inputfilename=argv[2];
    6466        outputfilename=argv[3];
    … …  
    8486        _printf_("initialize inputs:\n");
    8587       
    86         model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedthicknessAnalysisEnum);
    87         model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum);
    88         model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum);
     88        model->FindParam(&vx_g,NULL,NULL,"vx_g",BalancedthicknessAnalysisEnum);
     89        model->FindParam(&vy_g,NULL,NULL,"vy_g",BalancedthicknessAnalysisEnum);
     90        model->FindParam(&m_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum);
     91        model->FindParam(&a_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum);
     92        model->FindParam(&h_g,NULL,NULL,"h_g",BalancedthicknessAnalysisEnum);
     93        model->FindParam(&dhdt_g,NULL,NULL,"dhdt_g",BalancedthicknessAnalysisEnum);
    8994        model->FindParam(&numberofnodes,"numberofnodes");
    9095       
    9196        inputs=new ParameterInputs;
    92         inputs->Add("velocity",u_g_serial,3,numberofnodes);
    93         inputs->Add("melting",melting_g,1,numberofnodes);
    94         inputs->Add("accumulation",accumulation_g,1,numberofnodes);
     97        inputs->Add("vx",vx_g,1,numberofnodes);
     98        inputs->Add("vy",vy_g,1,numberofnodes);
     99        inputs->Add("melting",m_g,1,numberofnodes);
     100        inputs->Add("accumulation",a_g,1,numberofnodes);
     101        inputs->Add("thickness",h_g,1,numberofnodes);
     102        inputs->Add("dhdt",dhdt_g,1,numberofnodes);
    95103
    96104        _printf_("initialize results:\n");
    … …  
    138146
    139147        /*Free ressources:*/
     148        xfree((void**)&vx_g);
     149        xfree((void**)&vy_g);
     150        xfree((void**)&m_g);
     151        xfree((void**)&a_g);
     152        xfree((void**)&dhdt_g);
     153        xfree((void**)&h_g);
    140154        delete processedresults;
    141155        delete results;
     156        delete model;
    142157        delete inputs;
    143158

  • issm/trunk/src/c/parallel/balancedthickness_core.cpp

    r3567 r3582  
    1919
    2020        /*intermediary: */
    21         Vec u_g=NULL;
     21        Vec vx_g=NULL;
     22        Vec vy_g=NULL;
    2223
    2324        /*solutions: */
    … …  
    2829        int numberofdofspernode;
    2930        int numberofnodes;
    30         int dofs[2]={1,1};
     31        int numberofvertices;
     32        int dofs[1]={1};
    3133
    3234        /*fem balancedthickness model: */
    … …  
    3436
    3537
    36         /*recover fem model: */
    3738        fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
    3839
    … …  
    4041        model->FindParam(&verbose,"verbose");
    4142        model->FindParam(&numberofnodes,"numberofnodes");
     43        model->FindParam(&numberofvertices,"numberofvertices");
    4244        model->FindParam(&numberofdofspernode,"numberofdofspernode");
    4345
    4446        _printf_("depth averaging velocity...\n");
    45         u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
    46         FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
    47         inputs->Add("velocity_average",u_g,2,numberofnodes);
     47        vx_g=inputs->Get("vx",&dofs[0],1);
     48        vy_g=inputs->Get("vy",&dofs[0],1);
     49        /* NOT WORKING YET....
     50        FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
     51        FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
     52        */
     53        inputs->Add("vx_average",vx_g,1,numberofvertices);
     54        inputs->Add("vy_average",vy_g,1,numberofvertices);
    4855       
    4956        _printf_("call computational core:\n");
    5057        diagnostic_core_linear(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
    5158
    52         _printf_("extrude computed thickness on all layers:\n");
    53         FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
     59        _printf_("Averaging over vertices:\n");
     60        FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
     61
     62//      _printf_("extrude computed thickness on all layers:\n");
     63//      FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
    5464
    5565        /*Plug results into output dataset: */
    … …  
    5868
    5969        /*Free ressources:*/
    60         VecFree(&u_g);
     70        VecFree(&vx_g);
     71        VecFree(&vy_g);
    6172        VecFree(&h_g);
    6273}

  • issm/trunk/src/m/classes/@model/model.m

    r3359 r3582  
    161161        md.vel_obs_raw=NaN;
    162162        md.accumulation=NaN;
     163        md.dhdt=NaN;
    163164        md.geothermalflux=NaN;
    164165        md.observed_temperature=NaN;

  • issm/trunk/src/m/classes/public/display/displayobservations.m

    r1315 r3582  
    1818fielddisplay(md,'vel_obs_raw','raw observed magnitude [m/a]');
    1919fielddisplay(md,'accumulation','surface accumulation rate [m/a]');
     20fielddisplay(md,'dhdt','surface dhdt rate [m/a]');
    2021fielddisplay(md,'observed_temperature','observed temperature [K]');
    2122fielddisplay(md,'geothermalflux','geothermal heat flux [W/m^2]');

  • issm/trunk/src/m/classes/public/ismodelselfconsistent.m

    r3357 r3582  
    340340
    341341        %VELOCITIES MELTING AND ACCUMULATION
    342         fields={'vx','vy','accumulation','melting'};
    343         checksize(md,fields,[md.numberofgrids 1]);
    344         checknan(md,fields);
    345 
    346         %SPC
    347         if any(md.spcthickness(find(md.gridonboundary))~=1),
    348                 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);
    349         end
     342        fields={'vx','vy','accumulation','melting','dhdt'};
     343        checksize(md,fields,[md.numberofgrids 1]);
     344        checknan(md,fields);
     345
    350346end
    351347

  • issm/trunk/src/m/classes/public/marshall.m

    r3468 r3582  
    8787WriteData(fid,md.melting,'Mat','melting');
    8888WriteData(fid,md.accumulation,'Mat','accumulation');
     89WriteData(fid,md.dhdt,'Mat','dhdt');
    8990
    9091%Get materials

  • issm/trunk/src/m/solutions/jpl/balancedthickness.m

    r2715 r3582  
    1212       
    1313        displaystring(md.verbose,'%s',['reading balancedthickness model data']);
    14         md.analysis_type=BalancedthicknessAnalysisEnum; models.bt=CreateFemModel(md);
     14        md.analysis_type=BalancedthicknessAnalysisEnum; models.p=CreateFemModel(md);
    1515
    1616        % figure out number of dof: just for information purposes.
    … …  
    2020        displaystring(md.verbose,'\n%s',['setup inputs...']);
    2121        inputs=inputlist;
    22         inputs=add(inputs,'velocity',models.bt.parameters.u_g,'doublevec',3,models.bt.parameters.numberofnodes);
    23         inputs=add(inputs,'melting',models.bt.parameters.m_g,'doublevec',1,models.bt.parameters.numberofnodes);
    24         inputs=add(inputs,'accumulation',models.bt.parameters.a_g,'doublevec',1,models.bt.parameters.numberofnodes);
     22        inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofvertices);
     23        inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofvertices);
     24        inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofvertices);
     25        inputs=add(inputs,'dhdt',models.p.parameters.dhdt_g,'doublevec',1,models.p.parameters.numberofvertices);
     26        inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofvertices);
     27        inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofvertices);
    2528
    2629        displaystring(md.verbose,'\n%s',['call computational core:']);

  • issm/trunk/src/m/solutions/jpl/balancedthickness_core.m

    r3533 r3582  
    66
    77        %get FE model
    8         m=models.bt;
     8        m=models.p;
    99        results.time=0;
    1010        results.step=1;
    … …  
    1212        displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']);
    1313        %Take only the first two dofs of m.parameters.u_g
    14         u_g=get(inputs,'velocity',[1 1 0 0]);
    15         u_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,u_g,'velocity');
    16         inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
     14        vx_g=get(inputs,'vx',1);
     15        vy_g=get(inputs,'vy',1);
     16
     17        %NOT WORKING YET!!!!!
     18        %vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
     19        %vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
     20
     21        inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.numberofvertices);
     22        inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofvertices);
    1723
    1824        displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
    1925        results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
    2026
     27        displaystring(m.parameters.verbose,'\n%s',['averaging over vertices']);
     28        results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
     29
    2130        displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
    22         results.h_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g,'thickness',0);
     31        %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
    2332
    2433end %end function

  • issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m

    r3487 r3582  
    1313        %system matrices
    1414        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
     15        save A K_gg p_g
    1516        [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
    1617       

  • issm/trunk/src/m/utils/BC/SetIceSheetBC.m

    r3099 r3582  
    3636        disp('      no melting specified: values set as zero');
    3737end
     38if isnan(md.dhdt),
     39        md.dhdt=zeros(md.numberofgrids,1);
     40        disp('      no dhdt specified: values set as zero');
     41end
    3842
    3943displaystring(md.verbose,'%s',['      boundary conditions for prognostic model initialization']);

  • issm/trunk/src/m/utils/BC/SetIceShelfBC.m

    r3099 r3582  
    6060        disp('      no melting specified: values set as zero');
    6161end
     62if isnan(md.dhdt),
     63        md.dhdt=zeros(md.numberofgrids,1);
     64        disp('      no dhdt specified: values set as zero');
     65end
    6266
    6367displaystring(md.verbose,'%s',['      boundary conditions for prognostic model initialization']);

  • issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

    r3099 r3582  
    7171        disp('      no melting specified: values set as zero');
    7272end
     73if isnan(md.dhdt),
     74        md.dhdt=zeros(md.numberofgrids,1);
     75        disp('      no dhdt specified: values set as zero');
     76end
    7377
    7478displaystring(md.verbose,'%s',['      boundary conditions for prognostic model initialization']);
Note: See TracChangeset for help on using the changeset viewer.