Changeset 3588


Ignore:
Timestamp:
04/21/10 09:27:00 (15 years ago)
Author:
Mathieu Morlighem
Message:

Readded initial Balancedthickness and moved DG to Balancedthickness2

Location:
issm/trunk/src/c
Files:
7 added
19 edited

Legend:

Unmodified
Added
Removed

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

    r3446 r3588  
    1818        extern int my_rank;
    1919       
    20         //_printf_("      Configuring elements...\n");
     20        _printf_("      Configuring elements...\n");
    2121        elements->Configure(elements,loads,nodes,vertices,materials,parameters);
    22         //_printf_("      Configuring loads...\n");
     22        _printf_("      Configuring loads...\n");
    2323        loads->Configure(elements,loads,nodes,vertices,materials,parameters);
    24         //_printf_("      Configuring nodes...\n");
     24        _printf_("      Configuring nodes...\n");
    2525        nodes->Configure(elements,loads,nodes,vertices,materials,parameters);
    26         //_printf_("      Configuring parameters...\n");
     26        _printf_("      Configuring parameters...\n");
    2727        parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
    2828

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

    r3567 r3588  
    4848        Prognostic2AnalysisEnum,
    4949        BalancedthicknessAnalysisEnum,
     50        Balancedthickness2AnalysisEnum,
    5051        BalancedvelocitiesAnalysisEnum,
    5152        //melting

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

    r3529 r3588  
    245245                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
    246246                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
     247                                        ./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\
     248                                        ./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\
     249                                        ./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\
     250                                        ./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\
    247251                                        ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
    248252                                        ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
    … …  
    645649                                        ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\
    646650                                        ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\
     651                                        ./ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp\
     652                                        ./ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp\
     653                                        ./ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp\
     654                                        ./ModelProcessorx/Balancedthickness2/CreateParametersBalancedthickness2.cpp\
    647655                                        ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\
    648656                                        ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\
    … …  
    761769                                        ./parallel/prognostic2_core.cpp\
    762770                                        ./parallel/balancedthickness_core.cpp\
     771                                        ./parallel/balancedthickness2_core.cpp\
    763772                                        ./parallel/balancedvelocities_core.cpp\
    764773                                        ./parallel/slopecompute_core.cpp\
    … …  
    835844bin_PROGRAMS =
    836845else
    837 bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe
     846bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe prognostic2.exe balancedthickness.exe balancedthickness2.exe balancedvelocities.exe transient.exe steadystate.exe slopecompute.exe
    838847endif
    839848
    … …  
    858867balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    859868
     869balancedthickness2_exe_SOURCES = parallel/balancedthickness2.cpp
     870balancedthickness2_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     871
    860872balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp
    861873balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_

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

    r3582 r3588  
    1111
    1212void    CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
    13        
     13
     14        int i;
     15        int count=0;
    1416        DataSet* constraints = NULL;
    15        
     17
    1618        /*Create constraints: */
    1719        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:
    1845       
    1946        /*Assign output pointer: */

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

    r3582 r3588  
    88#include "../../objects/objects.h"
    99#include "../../shared/shared.h"
    10 #include "../../include/macros.h"
    1110#include "../../include/typedefs.h"
    1211#include "../../MeshPartitionx/MeshPartitionx.h"
    1312#include "../IoModel.h"
    1413
    15 void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
     14void    CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1615
    1716        /*Intermediary*/
    18         int i,j;
    19         int vertex_index;
    20         int node_index;
     17        int i,j,k,n;
    2118
    2219        /*DataSets: */
    23         DataSet* elements  = NULL;
    24         DataSet* nodes = NULL;
    25         DataSet* vertices = NULL;
    26         DataSet* materials = NULL;
     20        DataSet*    elements  = NULL;
     21        DataSet*    nodes = NULL;
     22        DataSet*    vertices = NULL;
     23        DataSet*    materials = NULL;
    2724
    2825        /*First create the elements, nodes and material properties: */
    … …  
    3532        Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
    3633
    37         /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
    3834        /*2d mesh: */
    3935        if (strcmp(iomodel->meshtype,"2d")==0){
    … …  
    5854                        }
    5955                }//for (i=0;i<numberofelements;i++)
    60        
     56
    6157                /*Free data : */
     58                xfree((void**)&iomodel->elements);
    6259                xfree((void**)&iomodel->thickness);
    6360                xfree((void**)&iomodel->surface);
    … …  
    6865        }
    6966        else{ //        if (strcmp(meshtype,"2d")==0)
    70                 ISSMERROR("not implemented yet");
     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
    7198        } //if (strcmp(meshtype,"2d")==0)
    7299
    73         /*Add new constrant material property tgo materials, at the end: */
     100        /*Add new constrant material property to materials, at the end: */
    74101        materials->AddObject(new Matpar(iomodel));
    75102
    76         /*Create nodes and vertices: */
     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        }
    77108        IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
    78109        IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
    … …  
    82113        IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
    83114        IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
    84         IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
    85115        IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
    86116        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         }
    91117
    92         /*Build Vertices dataset*/
    93118        for (i=0;i<iomodel->numberofvertices;i++){
    94119
    … …  
    99124                        vertices->AddObject(new Vertex(i,iomodel));
    100125
    101                 }
    102         }
     126                        /*Add node to nodes dataset: */
     127                        nodes->AddObject(new Node(i,iomodel));
    103128
    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                         }
    121129                }
    122130        }
    … …  
    131139        xfree((void**)&iomodel->gridonbed);
    132140        xfree((void**)&iomodel->gridonsurface);
    133         xfree((void**)&iomodel->gridonhutter);
    134141        xfree((void**)&iomodel->uppernodes);
    135142        xfree((void**)&iomodel->gridonicesheet);
    … …  
    139146         * datasets, it will not be redone: */
    140147        elements->Presort();
     148        vertices->Presort();
    141149        nodes->Presort();
    142         vertices->Presort();
    143150        materials->Presort();
    144151
    … …  
    150157        *pvertices=vertices;
    151158        *pmaterials=materials;
     159
    152160}

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

    r3582 r3588  
    88#include "../../shared/shared.h"
    99#include "../../include/macros.h"
    10 #include "../../include/typedefs.h"
    1110#include "../IoModel.h"
     11
    1212
    1313void    CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1414
    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];
     15        DataSet*    loads    = NULL;
    3016
    3117        /*Create loads: */
    3218        loads   = new DataSet(LoadsEnum);
    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 
     19       
     20       
    11121        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    11222         * datasets, it will not be redone: */
    11323        loads->Presort();
     24
     25        cleanup_and_return:
    11426
    11527        /*Assign output pointer: */
    … …  
    11729
    11830}
     31
     32

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

    r3582 r3588  
    1 /*!\file: CreateParametersPrognostic.cpp
    2  * \brief driver for creating parameters dataset, for prognostic analysis.
     1 /* \brief driver for creating parameters dataset, for prognostic analysis.
    32 */
    43
    … …  
    1615        int      count;
    1716        int      i;
    18         int      dim;
    19         double*  vx_g=NULL;
    20         double*  vy_g=NULL;
     17        double* u_g=NULL;
    2118
    2219        /*recover parameters : */
    … …  
    2825        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    2926        IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
     27        IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz");
    3028
    31         vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
    32         vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
     29        u_g=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double));
    3330
    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;
     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;
    3634
    3735        count++;
    38         param= new Param(count,"vx_g",DOUBLEVEC);
    39         param->SetDoubleVec(vx_g,iomodel->numberofvertices,1);
     36        param= new Param(count,"u_g",DOUBLEVEC);
     37        param->SetDoubleVec(u_g,3*iomodel->numberofvertices,3);
    4038        parameters->AddObject(param);
    41         count++;
    42         param= new Param(count,"vy_g",DOUBLEVEC);
    43         param->SetDoubleVec(vy_g,iomodel->numberofvertices,1);
    44         parameters->AddObject(param);
     39
    4540
    4641        xfree((void**)&iomodel->vx);
    4742        xfree((void**)&iomodel->vy);
    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);
     43        xfree((void**)&iomodel->vz);
     44        xfree((void**)&u_g);
    7545
    7646        /*Get melting: */

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

    r3570 r3588  
    100100
    101101        }
     102        else if (iomodel->analysis_type==Balancedthickness2AnalysisEnum){
     103
     104                CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
     105                CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle);
     106                CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle);
     107                CreateParametersBalancedthickness2(pparameters,iomodel,iomodel_handle);
     108
     109        }
    102110        else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum){
    103111

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

    r3582 r3588  
    247247        if (
    248248                                iomodel->analysis_type==Prognostic2AnalysisEnum ||
    249                                 iomodel->analysis_type==BalancedthicknessAnalysisEnum
     249                                iomodel->analysis_type==Balancedthickness2AnalysisEnum
    250250                                )
    251251         param->SetDouble(3*iomodel->numberofelements);

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

    r3582 r3588  
    263263        void  CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
    264264
     265        /*balancedthickness2:*/
     266        void    CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
     267        void    CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
     268        void  CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     269        void  CreateParametersBalancedthickness2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
     270
    265271        /*balancedvelocities:*/
    266272        void    CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);

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

    r3582 r3588  
    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 || iomodel->analysis_type==BalancedthicknessAnalysisEnum)
     24        if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum)
    2525                DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle);
    2626        else

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

    r3582 r3588  
    225225                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
    226226        }
    227         else if (analysis_type==BalancedthicknessAnalysisEnum){
     227        else if (analysis_type==Balancedthickness2AnalysisEnum){
    228228                /*No transient term is involved*/
    229229                dt=1;
    … …  
    256256                GetParameterValue(&vy, &vy_list[0],gauss_coord);
    257257                UdotN=vx*normal[0]+vy*normal[1];
     258                if (fabs(UdotN)<1.0e-9) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN);
     259                //if (fabs(UdotN)>0 && fabs(UdotN)< 1.0e-7) UdotN= 1.0e-7;
     260                //if (fabs(UdotN)<0 && fabs(UdotN)>-1.0e-7) UdotN=-1.0e-7;
    258261
    259262                /*Get L and B: */
    … …  
    341344                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
    342345        }
    343         else if (analysis_type==BalancedthicknessAnalysisEnum){
     346        else if (analysis_type==Balancedthickness2AnalysisEnum){
    344347                /*No transient term is involved*/
    345348                dt=1;
    … …  
    483486                if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!");
    484487        }
    485         else if (analysis_type==BalancedthicknessAnalysisEnum){
     488        else if (analysis_type==Balancedthickness2AnalysisEnum){
    486489                /*No transient term is involved*/
    487490                dt=1;

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

    r3582 r3588  
    7777        /*hooks: */
    7878        //go recover node ids, needed to initialize the node hook.
    79         if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum){
     79        if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){
    8080                /*Discontinuous Galerkin*/
    8181                tria_node_ids[0]=3*i+1;
    … …  
    598598                CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
    599599        }
     600        else if (analysis_type==Balancedthickness2AnalysisEnum){
     601
     602                CreateKMatrixBalancedthickness2( Kgg,inputs,analysis_type,sub_analysis_type);
     603        }
    600604        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    601605
    … …  
    611615/*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/
    612616void  Tria::CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     617
     618        /* local declarations */
     619        int             i,j;
     620
     621        /* node data: */
     622        const int    numgrids=3;
     623        const int    NDOF1=1;
     624        const int    numdof=NDOF1*numgrids;
     625        double       xyz_list[numgrids][3];
     626        int          doflist[numdof];
     627        int          numberofdofspernode;
     628
     629        /* gaussian points: */
     630        int     num_gauss,ig;
     631        double* first_gauss_area_coord  =  NULL;
     632        double* second_gauss_area_coord =  NULL;
     633        double* third_gauss_area_coord  =  NULL;
     634        double* gauss_weights           =  NULL;
     635        double  gauss_weight;
     636        double  gauss_l1l2l3[3];
     637
     638        /* matrices: */
     639        double L[numgrids];
     640        double B[2][numgrids];
     641        double Bprime[2][numgrids];
     642        double DL[2][2]={0.0};
     643        double DLprime[2][2]={0.0};
     644        double DL_scalar;
     645        double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
     646        double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     647        double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     648        double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     649
     650        double Jdettria;
     651
     652        /*input parameters for structural analysis (diagnostic): */
     653        double  vxvy_list[numgrids][2]={0.0};
     654        double  vx_list[numgrids]={0.0};
     655        double  vy_list[numgrids]={0.0};
     656        double  dvx[2];
     657        double  dvy[2];
     658        double  vx,vy;
     659        double  dvxdx,dvydy;
     660        double  v_gauss[2]={0.0};
     661        double  K[2][2]={0.0};
     662        double  KDL[2][2]={0.0};
     663        int     dofs[2]={0,1};
     664        int     found=0;
     665
     666        /*dynamic objects pointed to by hooks: */
     667        Node**  nodes=NULL;
     668        Matpar* matpar=NULL;
     669        Matice* matice=NULL;
     670        Numpar* numpar=NULL;
     671
     672        ParameterInputs* inputs=NULL;
     673
     674        /*recover pointers: */
     675        inputs=(ParameterInputs*)vinputs;
     676
     677        /*recover objects from hooks: */
     678        nodes=(Node**)hnodes.deliverp();
     679        matpar=(Matpar*)hmatpar.delivers();
     680        matice=(Matice*)hmatice.delivers();
     681        numpar=(Numpar*)hnumpar.delivers();
     682
     683        /*recover extra inputs from users, at current convergence iteration: */
     684        found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     685        if(!found)ISSMERROR(" could not find velocity_average  in inputs!");
     686
     687        for(i=0;i<numgrids;i++){
     688                vx_list[i]=vxvy_list[i][0];
     689                vy_list[i]=vxvy_list[i][1];
     690        }
     691
     692        /* Get node coordinates and dof list: */
     693        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     694        GetDofList(&doflist[0],&numberofdofspernode);
     695
     696        //Create Artificial diffusivity once for all if requested
     697        if(numpar->artdiff){
     698                //Get the Jacobian determinant
     699                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
     700                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     701
     702                //Build K matrix (artificial diffusivity matrix)
     703                v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     704                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     705
     706                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
     707                K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
     708        }
     709
     710        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     711        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     712
     713        /* Start  looping on the number of gaussian points: */
     714        for (ig=0; ig<num_gauss; ig++){
     715                /*Pick up the gaussian point: */
     716                gauss_weight=*(gauss_weights+ig);
     717                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     718                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     719                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     720
     721                /* Get Jacobian determinant: */
     722                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     723
     724                /*Get B  and B prime matrix: */
     725                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     726                GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     727
     728                //Get vx, vy and their derivatives at gauss point
     729                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     730                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     731
     732                GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
     733                GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
     734
     735                dvxdx=dvx[0];
     736                dvydy=dvy[1];
     737
     738                DL_scalar=gauss_weight*Jdettria;
     739
     740                //Create DL and DLprime matrix
     741                DL[0][0]=DL_scalar*dvxdx;
     742                DL[1][1]=DL_scalar*dvydy;
     743
     744                DLprime[0][0]=DL_scalar*vx;
     745                DLprime[1][1]=DL_scalar*vy;
     746
     747                //Do the triple product tL*D*L.
     748                //Ke_gg_thickness=B'*DLprime*Bprime;
     749
     750                TripleMultiply( &B[0][0],2,numdof,1,
     751                                        &DL[0][0],2,2,0,
     752                                        &B[0][0],2,numdof,0,
     753                                        &Ke_gg_thickness1[0][0],0);
     754
     755                TripleMultiply( &B[0][0],2,numdof,1,
     756                                        &DLprime[0][0],2,2,0,
     757                                        &Bprime[0][0],2,numdof,0,
     758                                        &Ke_gg_thickness2[0][0],0);
     759
     760                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     761                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
     762                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
     763
     764                if(numpar->artdiff){
     765
     766                        /* Compute artificial diffusivity */
     767                        KDL[0][0]=DL_scalar*K[0][0];
     768                        KDL[1][1]=DL_scalar*K[1][1];
     769
     770                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     771                                                &KDL[0][0],2,2,0,
     772                                                &Bprime[0][0],2,numdof,0,
     773                                                &Ke_gg_gaussian[0][0],0);
     774
     775                        /* Add artificial diffusivity matrix */
     776                        for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     777
     778                }
     779
     780#ifdef _DEBUGELEMENTS_
     781                if(my_rank==RANK && id==ELID){
     782                        printf("      B:\n");
     783                        for(i=0;i<3;i++){
     784                                for(j=0;j<numdof;j++){
     785                                        printf("%g ",B[i][j]);
     786                                }
     787                                printf("\n");
     788                        }
     789                        printf("      Bprime:\n");
     790                        for(i=0;i<3;i++){
     791                                for(j=0;j<numdof;j++){
     792                                        printf("%g ",Bprime[i][j]);
     793                                }
     794                                printf("\n");
     795                        }
     796                }
     797#endif
     798        } // for (ig=0; ig<num_gauss; ig++)
     799
     800        /*Add Ke_gg to global matrix Kgg: */
     801        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     802
     803#ifdef _DEBUGELEMENTS_
     804        if(my_rank==RANK && id==ELID){
     805                printf("      Ke_gg erms:\n");
     806                for( i=0; i<numdof; i++){
     807                        for (j=0;j<numdof;j++){
     808                                printf("%g ",Ke_gg[i][j]);
     809                        }
     810                        printf("\n");
     811                }
     812                printf("      Ke_gg row_indices:\n");
     813                for( i=0; i<numdof; i++){
     814                        printf("%i ",doflist[i]);
     815                }
     816
     817        }
     818#endif
     819
     820cleanup_and_return:
     821        xfree((void**)&first_gauss_area_coord);
     822        xfree((void**)&second_gauss_area_coord);
     823        xfree((void**)&third_gauss_area_coord);
     824        xfree((void**)&gauss_weights);
     825
     826}
     827/*}}}*/
     828/*FUNCTION Tria::CreateKMatrixBalancedthickness2 {{{1*/
     829void  Tria::CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    613830
    614831        /* local declarations */
    … …  
    21092326                CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type);
    21102327        }
     2328        else if (analysis_type==Balancedthickness2AnalysisEnum){
     2329
     2330                CreatePVectorBalancedthickness2( pg,inputs,analysis_type,sub_analysis_type);
     2331        }
    21112332        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    21122333
    … …  
    21212342/*FUNCTION Tria::CreatePVectorBalancedthickness {{{1*/
    21222343void  Tria::CreatePVectorBalancedthickness(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
     2344
     2345
     2346        /* local declarations */
     2347        int             i,j;
     2348
     2349        /* node data: */
     2350        const int    numgrids=3;
     2351        const int    NDOF1=1;
     2352        const int    numdof=NDOF1*numgrids;
     2353        double       xyz_list[numgrids][3];
     2354        int          doflist[numdof];
     2355        int          numberofdofspernode;
     2356
     2357        /* gaussian points: */
     2358        int     num_gauss,ig;
     2359        double* first_gauss_area_coord  =  NULL;
     2360        double* second_gauss_area_coord =  NULL;
     2361        double* third_gauss_area_coord  =  NULL;
     2362        double* gauss_weights           =  NULL;
     2363        double  gauss_weight;
     2364        double  gauss_l1l2l3[3];
     2365
     2366        /* matrix */
     2367        double pe_g[numgrids]={0.0};
     2368        double L[numgrids];
     2369        double Jdettria;
     2370
     2371        /*input parameters for structural analysis (diagnostic): */
     2372        double  accumulation_list[numgrids]={0.0};
     2373        double  accumulation_g;
     2374        double  melting_list[numgrids]={0.0};
     2375        double  melting_g;
     2376        double  thickness_list[numgrids]={0.0};
     2377        double  thickness_g;
     2378        int     dofs[1]={0};
     2379        int     found=0;
     2380
     2381        /*dynamic objects pointed to by hooks: */
     2382        Node**  nodes=NULL;
     2383        Matpar* matpar=NULL;
     2384        Matice* matice=NULL;
     2385        Numpar* numpar=NULL;
     2386
     2387        ParameterInputs* inputs=NULL;
     2388
     2389        /*recover pointers: */
     2390        inputs=(ParameterInputs*)vinputs;
     2391
     2392        /*recover objects from hooks: */
     2393        nodes=(Node**)hnodes.deliverp();
     2394        matpar=(Matpar*)hmatpar.delivers();
     2395        matice=(Matice*)hmatice.delivers();
     2396        numpar=(Numpar*)hnumpar.delivers();
     2397
     2398        /*recover extra inputs from users, at current convergence iteration: */
     2399        found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
     2400        if(!found)ISSMERROR(" could not find accumulation in inputs!");
     2401        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
     2402        if(!found)ISSMERROR(" could not find melting in inputs!");
     2403
     2404        /* Get node coordinates and dof list: */
     2405        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     2406        GetDofList(&doflist[0],&numberofdofspernode);
     2407
     2408        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2409        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2410
     2411        /* Start  looping on the number of gaussian points: */
     2412        for (ig=0; ig<num_gauss; ig++){
     2413                /*Pick up the gaussian point: */
     2414                gauss_weight=*(gauss_weights+ig);
     2415                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2416                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2417                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2418
     2419                /* Get Jacobian determinant: */
     2420                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     2421
     2422                /*Get L matrix: */
     2423                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     2424
     2425                /* Get accumulation, melting and thickness at gauss point */
     2426                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
     2427                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
     2428
     2429                /* Add value into pe_g: */
     2430                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(accumulation_g-melting_g)*L[i];
     2431
     2432        } // for (ig=0; ig<num_gauss; ig++)
     2433
     2434        /*Add pe_g to global matrix Kgg: */
     2435        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     2436
     2437cleanup_and_return:
     2438        xfree((void**)&first_gauss_area_coord);
     2439        xfree((void**)&second_gauss_area_coord);
     2440        xfree((void**)&third_gauss_area_coord);
     2441        xfree((void**)&gauss_weights);
     2442
     2443}
     2444/*}}}*/
     2445/*FUNCTION Tria::CreatePVectorBalancedthickness2 {{{1*/
     2446void  Tria::CreatePVectorBalancedthickness2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
    21232447
    21242448

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

    r3529 r3588  
    117117                void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    118118                void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
     119                void  CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
     120                void  CreatePVectorBalancedthickness2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
    119121                void  CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    120122                void  CreatePVectorBalancedvelocities(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);

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

    r3567 r3588  
    135135        if(analysis_type==BalancedthicknessAnalysisEnum){
    136136                fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
     137        }
     138        if(analysis_type==Balancedthickness2AnalysisEnum){
     139                fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
    137140        }
    138141        if(analysis_type==BalancedvelocitiesAnalysisEnum){

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

    r3582 r3588  
    2525        Model* model=NULL;
    2626
    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;
     27        Vec     u_g=NULL;
     28        double* u_g_serial=NULL;
     29        double* melting_g=NULL;
     30        double* accumulation_g=NULL;
    3331        double  dt;
    3432        double  yts;
    … …  
    6260        MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
    6361
    64         _printf_("recover input file name and output file name:\n");
     62        _printf_("recover , input file name and output file name:\n");
    6563        inputfilename=argv[2];
    6664        outputfilename=argv[3];
    … …  
    8684        _printf_("initialize inputs:\n");
    8785       
    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);
     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);
    9489        model->FindParam(&numberofnodes,"numberofnodes");
    9590       
    9691        inputs=new ParameterInputs;
    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);
     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);
    10395
    10496        _printf_("initialize results:\n");
    … …  
    146138
    147139        /*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);
    154140        delete processedresults;
    155141        delete results;
    156         delete model;
    157142        delete inputs;
    158143

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

    r3582 r3588  
    1919
    2020        /*intermediary: */
    21         Vec vx_g=NULL;
    22         Vec vy_g=NULL;
     21        Vec u_g=NULL;
    2322
    2423        /*solutions: */
    … …  
    2928        int numberofdofspernode;
    3029        int numberofnodes;
    31         int numberofvertices;
    32         int dofs[1]={1};
     30        int dofs[2]={1,1};
    3331
    3432        /*fem balancedthickness model: */
    … …  
    3634
    3735
     36        /*recover fem model: */
    3837        fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
    3938
    … …  
    4140        model->FindParam(&verbose,"verbose");
    4241        model->FindParam(&numberofnodes,"numberofnodes");
    43         model->FindParam(&numberofvertices,"numberofvertices");
    4442        model->FindParam(&numberofdofspernode,"numberofdofspernode");
    4543
    4644        _printf_("depth averaging velocity...\n");
    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);
     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);
    5548       
    5649        _printf_("call computational core:\n");
    5750        diagnostic_core_linear(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
    5851
    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);
     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);
    6454
    6555        /*Plug results into output dataset: */
    … …  
    6858
    6959        /*Free ressources:*/
    70         VecFree(&vx_g);
    71         VecFree(&vy_g);
     60        VecFree(&u_g);
    7261        VecFree(&h_g);
    7362}

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

    r3373 r3588  
    1919void prognostic2_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2020void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs);
     21void balancedthickness2_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2122void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs);
    2223void slopecompute_core(DataSet* results,Model* model, ParameterInputs* inputs);

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

    r3570 r3588  
    5959                numdofs=1;
    6060        }
     61        else if (analysis_type==Balancedthickness2AnalysisEnum){
     62                numdofs=1;
     63        }
    6164        else if (analysis_type==BalancedvelocitiesAnalysisEnum){
    6265                numdofs=1;
Note: See TracChangeset for help on using the changeset viewer.