Changeset 19113


Ignore:
Timestamp:
02/17/15 15:10:59 (10 years ago)
Author:
Mathieu Morlighem
Message:

BUG: fixed problem when spcvz is on multiple column

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

    r18930 r19113  
    1010
    1111        /*Intermediary*/
    12         int count;
    13         IssmDouble yts;
    14 
    15         /*Fetch parameters: */
    16         iomodel->Constant(&yts,ConstantsYtsEnum);
     12        bool        isSIA,isSSA,isL1L2,isHO,isFS,iscoupling;
     13        int         Mz,Nz;
     14        IssmDouble *spcvz = NULL;
    1715
    1816        /*return if not 3d mesh*/
    1917        if(iomodel->domaintype!=Domain3DEnum) return;
    2018
    21         /*Fetch data: */
    22         iomodel->FetchData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
    23 
    24         /*Initialize counter*/
    25         count=0;
    26 
    27         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    28         for(int i=0;i<iomodel->numberofvertices;i++){
    29 
    30                 /*keep only this partition's nodes:*/
    31                 if(iomodel->my_vertices[i]){
    32 
    33                         if (reCast<int,IssmDouble>(iomodel->Data(FlowequationBorderFSEnum)[i])){
    34                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
    35                                 count++;
    36                         }
    37                         else if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvzEnum)[i])){
    38                                 constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,
    39                                                                 iomodel->Data(StressbalanceSpcvzEnum)[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
    40                                 count++;
    41 
    42                         }
    43                 }
    44         }
    45 
    46         /*Free data: */
    47         iomodel->DeleteData(2,StressbalanceSpcvzEnum,FlowequationBorderFSEnum);
     19        /*Fetch parameters: */
     20        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);
     21        iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
     22        iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
     23        iomodel->Constant(&isHO,FlowequationIsHOEnum);
     24        iomodel->Constant(&isFS,FlowequationIsFSEnum);
     25
     26        /*Do we have coupling*/
     27        if((isSIA?1.:0.) + (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
     28         iscoupling = true;
     29        else
     30         iscoupling = false;
     31
     32
     33        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
     34        if(!iscoupling){
     35                IoModelToConstraintsx(constraints,iomodel,StressbalanceSpcvzEnum,StressbalanceVerticalAnalysisEnum,P1Enum,0);
     36        }
     37        else{
     38                /*Fetch data: */
     39                iomodel->FetchData(1,FlowequationBorderFSEnum);
     40                /*Fetch Spc*/
     41                iomodel->FetchData(&spcvz,&Mz,&Nz,StressbalanceSpcvzEnum);
     42                if(Nz>1) _error_("not supported yet (needs to be coded)");
     43
     44                /*Initialize counter*/
     45                int count=0;
     46
     47                /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     48                for(int i=0;i<iomodel->numberofvertices;i++){
     49
     50                        /*keep only this partition's nodes:*/
     51                        if(iomodel->my_vertices[i]){
     52
     53                                if (reCast<int,IssmDouble>(iomodel->Data(FlowequationBorderFSEnum)[i])){
     54                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,0,StressbalanceVerticalAnalysisEnum)); //spc to zero as vertical velocity is done in Horiz for FS
     55                                        count++;
     56                                }
     57                                else if (!xIsNan<IssmDouble>(spcvz[i])){
     58                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,0,
     59                                                                        spcvz[i],StressbalanceVerticalAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
     60                                        count++;
     61
     62                                }
     63                        }
     64                }
     65
     66                /*Free data: */
     67                iomodel->DeleteData(1,FlowequationBorderFSEnum);
     68                iomodel->DeleteData(spcvz,StressbalanceSpcvzEnum);
     69        }
    4870
    4971}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.