Changeset 3534


Ignore:
Timestamp:
04/14/10 11:59:07 (15 years ago)
Author:
Mathieu Morlighem
Message:

Some bug fix in prognostic2, still not working

Location:
issm/trunk/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp

    r3440 r3534  
    1616
    1717        /*Intermediary*/
    18         int i,j,k,n;
     18        int i,j;
     19        int vertex_index;
     20        int node_index;
    1921
    2022        /*DataSets: */
    … …  
    106108                        if(iomodel->my_nodes[3*i+j]){
    107109
    108                                 //Get id of the vertex on which the current node is located
    109                                 k=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
    110                                 ISSMASSERT(k>=0 && k<iomodel->numberofvertices);
     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;
    111116
    112117                                /*Add node to nodes dataset: */
    113                                 nodes->AddObject(new Node(k,i,iomodel));
     118                                nodes->AddObject(new Node(vertex_index,node_index,iomodel));
    114119
    115120                        }

  • issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp

    r3459 r3534  
    1 /*!\file: CreateParametersPrognostic2.cpp
     1/*!\file: CreateParametersPrognostic.cpp
    22 * \brief driver for creating parameters dataset, for prognostic analysis.
    33 */
    … …  
    99#include "../../objects/objects.h"
    1010#include "../../shared/shared.h"
    11 #include "../../include/macros.h"
    1211#include "../IoModel.h"
    1312
    … …  
    1716        DataSet* parameters=NULL;
    1817        int      count;
    19         int      i,j;
     18        int      i;
    2019        int      dim;
    21         int*     part;
    22 
    23         double* vx_g=NULL;
    24         double* vy_g=NULL;
    25         double* h_g=NULL;
    26         double* a_g=NULL;
    27         double* m_g=NULL;
     20        double*  vx_g=NULL;
     21        double*  vy_g=NULL;
    2822
    2923        /*recover parameters : */
    … …  
    3226        count=parameters->Size();
    3327
    34         /*Create transformation vector for DG inputs*/
    35         IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
    36         part=(int*)xcalloc(iomodel->numberofelements*3,sizeof(int));
    37         for(i=0;i<iomodel->numberofelements;i++){
    38                 for(j=0;j<3;j++){
    39                         part[3*i+j]=(int)iomodel->elements[3*i+j]-1; //Matlab to C indexing
    40                         ISSMASSERT(part[3*i+j]<iomodel->numberofvertices);
    41                 }
    42         }
    43         xfree((void**)&iomodel->elements);
     28        /*Get vx and vy: */
     29        IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
     30        IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
    4431
    45         /*Get vx: */
    46         IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
    47         vx_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
    48         if(iomodel->vx) for(i=0;i<3*iomodel->numberofelements;i++) vx_g[i]=iomodel->vx[part[i]]/iomodel->yts;
     32        vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
     33        vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double));
     34
     35        if(iomodel->vx) for(i=0;i<iomodel->numberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts;
     36        if(iomodel->vy) for(i=0;i<iomodel->numberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts;
    4937
    5038        count++;
    5139        param= new Param(count,"vx_g",DOUBLEVEC);
    52         param->SetDoubleVec(vx_g,3*iomodel->numberofelements,1);
     40        param->SetDoubleVec(vx_g,iomodel->numberofvertices,1);
    5341        parameters->AddObject(param);
    54         xfree((void**)&iomodel->vx);
    55         xfree((void**)&vx_g);
    56 
    57         /*Get vy: */
    58         IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
    59         vy_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
    60         if(iomodel->vy) for(i=0;i<3*iomodel->numberofelements;i++) vy_g[i]=iomodel->vy[part[i]]/iomodel->yts;
    61 
    6242        count++;
    6343        param= new Param(count,"vy_g",DOUBLEVEC);
    64         param->SetDoubleVec(vy_g,3*iomodel->numberofelements,1);
     44        param->SetDoubleVec(vy_g,iomodel->numberofvertices,1);
    6545        parameters->AddObject(param);
     46
     47        xfree((void**)&iomodel->vx);
    6648        xfree((void**)&iomodel->vy);
     49        xfree((void**)&vx_g);
    6750        xfree((void**)&vy_g);
    6851
    6952        /*Get thickness: */
    7053        IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
    71         h_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
    72         if(iomodel->thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=iomodel->thickness[part[i]];
    73 
     54       
    7455        count++;
    7556        param= new Param(count,"h_g",DOUBLEVEC);
    76         param->SetDoubleVec(h_g,3*iomodel->numberofelements,1);
     57        if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1);
     58        else param->SetDoubleVec(iomodel->thickness,0,0);
    7759        parameters->AddObject(param);
     60
     61        /*Free thickness: */
    7862        xfree((void**)&iomodel->thickness);
    79         xfree((void**)&h_g);
     63
     64        /*Get surface: */
     65        IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
     66       
     67        count++;
     68        param= new Param(count,"s_g",DOUBLEVEC);
     69        if(iomodel->surface) param->SetDoubleVec(iomodel->surface,iomodel->numberofvertices,1);
     70        else param->SetDoubleVec(iomodel->surface,0,0);
     71        parameters->AddObject(param);
     72
     73        /*Free surface: */
     74        xfree((void**)&iomodel->surface);
     75
     76        /*Get bed: */
     77        IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
     78       
     79        count++;
     80        param= new Param(count,"b_g",DOUBLEVEC);
     81        if(iomodel->bed) param->SetDoubleVec(iomodel->bed,iomodel->numberofvertices,1);
     82        else param->SetDoubleVec(iomodel->bed,0,0);
     83        parameters->AddObject(param);
     84
     85        /*Free bed: */
     86        xfree((void**)&iomodel->bed);
     87
     88        /*Get melting: */
     89        IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting");
     90        if(iomodel->melting) for(i=0;i<iomodel->numberofvertices;i++)iomodel->melting[i]=iomodel->melting[i]/iomodel->yts;
     91       
     92        count++;
     93        param= new Param(count,"m_g",DOUBLEVEC);
     94        if(iomodel->melting) param->SetDoubleVec(iomodel->melting,iomodel->numberofvertices,1);
     95        else param->SetDoubleVec(iomodel->melting,0,1);
     96        parameters->AddObject(param);
     97
     98        /*Free melting: */
     99        xfree((void**)&iomodel->melting);
    80100
    81101        /*Get accumulation: */
    82102        IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation");
    83         a_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
    84         if(iomodel->accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=iomodel->accumulation[part[i]]/iomodel->yts;
    85 
     103        if(iomodel->accumulation) for(i=0;i<iomodel->numberofvertices;i++)iomodel->accumulation[i]=iomodel->accumulation[i]/iomodel->yts;
     104       
    86105        count++;
    87106        param= new Param(count,"a_g",DOUBLEVEC);
    88         param->SetDoubleVec(a_g,3*iomodel->numberofelements,1);
     107        if(iomodel->accumulation) param->SetDoubleVec(iomodel->accumulation,iomodel->numberofvertices,1);
     108        else param->SetDoubleVec(iomodel->accumulation,0,0);
    89109        parameters->AddObject(param);
     110
     111        /*Free accumulation: */
    90112        xfree((void**)&iomodel->accumulation);
    91         xfree((void**)&a_g);
    92 
    93         /*Get melting: */
    94         IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting");
    95         m_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
    96         if(iomodel->melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=iomodel->melting[part[i]]/iomodel->yts;
    97 
    98         count++;
    99         param= new Param(count,"m_g",DOUBLEVEC);
    100         param->SetDoubleVec(m_g,3*iomodel->numberofelements,1);
    101         parameters->AddObject(param);
    102         xfree((void**)&iomodel->melting);
    103         xfree((void**)&m_g);
    104 
    105         /*Free partitioning vector*/
    106         xfree((void**)&part);
    107113
    108114        /*Assign output pointer: */

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

    r3332 r3534  
    2222
    2323        /*diverse: */
    24         int     numberofnodes;
     24        int numberofvertices;
    2525        int i,j,k;
    2626        int ndof;
    … …  
    3232
    3333        /*Some parameters needed: */
    34         parameters->FindParam(&numberofnodes,"numberofnodes");
     34        parameters->FindParam(&numberofvertices,"numberofvertices");
    3535
    3636        /*serialize partition vector: */
    … …  
    5151                                if (ndof!=0){ /*ok, we are dealing with a parameter that needs to be repartitioned, for ndof degrees of freedom: */
    5252
    53                                         new_parameter=(double*)xcalloc(numberofnodes*ndof,sizeof(double));
     53                                        new_parameter=(double*)xcalloc(numberofvertices*ndof,sizeof(double));
    5454                                        param->GetParameterValue(&parameter);
    5555
    5656                                        for(j=0;j<ndof;j++){
    5757
    58                                                 for(k=0;k<numberofnodes;k++){
     58                                                for(k=0;k<numberofvertices;k++){
    5959                                               
    6060                                                        new_parameter[(int)(ndof*partition[k]+j)]=parameter[ndof*k+j];
    … …  
    6666                                       
    6767                                        /*Now, replace old parameter with new parameter: */
    68                                         param->SetDoubleVec(new_parameter,ndof*numberofnodes,ndof);
     68                                        param->SetDoubleVec(new_parameter,ndof*numberofvertices,ndof);
    6969
    7070                                        /*Free ressources: */

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

    r3522 r3534  
    580580                        for(;;){
    581581
    582                                 if (node->IsOnSurface())break;
     582                                if (node->IsOnSurface()) break;
    583583
    584584                                vertexdof=node->GetVertexDof();
    … …  
    700700                        GetDofList(&doflist,&numberofdofspernode);
    701701
    702                         //initilaize node
     702                        //initialize node
    703703                        node=this;
    704704
    … …  
    707707
    708708                        //go through all nodes which sit on top of this node, until we reach the surface,
    709                         //and pltg  fieldnode in field:
     709                        //and plug  fieldnode in field:
    710710                        for(;;){
    711711

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

    r3529 r3534  
    7676
    7777        /*hooks: */
    78         for(j=0;j<3;j++){ //go recover node ids, needed to initialize the node hook.
    79                 tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab
     78        //go recover node ids, needed to initialize the node hook.
     79        if (iomodel->analysis_type==Prognostic2AnalysisEnum()){
     80                /*Discontinuous Galerkin*/
     81                tria_node_ids[0]=3*i+1;
     82                tria_node_ids[1]=3*i+2;
     83                tria_node_ids[2]=3*i+3;
     84        }
     85        else{
     86                /*Continuous Galerkin*/
     87                for(j=0;j<3;j++){
     88                        tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab
     89                }
    8090        }
    8191        tria_matice_id=i+1; //refers to the corresponding ice material object

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

    r3446 r3534  
    2828        int numberofdofspernode;
    2929        int numberofnodes;
     30        int numberofvertices;
    3031        int dofs[1]={1};
    3132
    … …  
    3940        model->FindParam(&verbose,"verbose");
    4041        model->FindParam(&numberofnodes,"numberofnodes");
     42        model->FindParam(&numberofvertices,"numberofvertices");
    4143        model->FindParam(&numberofdofspernode,"numberofdofspernode");
    4244
    … …  
    4446        vx_g=inputs->Get("vx",&dofs[0],1);
    4547        vy_g=inputs->Get("vy",&dofs[0],1);
     48        /* NOT WORKING YET....
    4649        FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
    4750        FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
    48         inputs->Add("vx_average",vx_g,1,numberofnodes);
    49         inputs->Add("vy_average",vy_g,1,numberofnodes);
     51        */
     52        inputs->Add("vx_average",vx_g,1,numberofvertices);
     53        inputs->Add("vy_average",vy_g,1,numberofvertices);
    5054       
    5155        _printf_("call computational core:\n");

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

    r3359 r3534  
    2020        displaystring(md.verbose,'\n%s',['setup inputs...']);
    2121        inputs=inputlist;
    22         inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofnodes);
    23         inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofnodes);
    24         inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes);
    25         inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes);
    26         inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.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,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofvertices);
     26        inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofvertices);
    2727        inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double');
    2828

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

    r3481 r3534  
    1414        vx_g=get(inputs,'vx',1);
    1515        vy_g=get(inputs,'vy',1);
    16         vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
    17         vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
    18         inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.numberofnodes);
    19         inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofnodes);
     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);
    2023
    2124        displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
Note: See TracChangeset for help on using the changeset viewer.