source: issm/oecreview/Archive/15392-16133/ISSM-15483-15484.diff

Last change on this file was 16134, checked in by Mathieu Morlighem, 12 years ago

Added Archive/15392-16133
File size: 12.6 KB

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp

     
    1515
    1616        /*Do not add constraints in DG, they are weakly imposed*/
    1717        if(stabilization!=3){
    18                 IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum);
     18                IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum,P1Enum);
    1919        }
    2020
    2121        /*Assign output pointer: */

  • ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp

     
    1919
    2020        if(hydrology_model!=HydrologyshreveEnum) return;
    2121
    22         IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum);
     22        IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum);
    2323
    2424        /*Assign output pointer: */
    2525        *pconstraints=constraints;

  • ../trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

     
    1010void    CreateConstraintsDiagnosticHoriz(Constraints** pconstraints, IoModel* iomodel){
    1111
    1212        /*Intermediary*/
    13         int     i,j;
    14         int     count;
    15         IssmDouble  yts;
    16         IssmDouble  g;
    17         IssmDouble  rho_ice;
    18         IssmDouble  stokesreconditioning;
    19         bool    isstokes,isl1l2,ismacayealpattyn;
    20         int     fe_ssa;
    21    bool    spcpresent=false;
    22         int Mx,Nx;
    23         int My,Ny;
    24         int Mz,Nz;
     13        int        i,j;
     14        int        count;
     15        IssmDouble yts;
     16        IssmDouble g;
     17        IssmDouble rho_ice;
     18        IssmDouble stokesreconditioning;
     19        bool       isstokes,isl1l2,ismacayealpattyn;
     20        int        fe_ssa;
     21        bool       spcpresent = false;
     22        int        Mx,Nx;
     23        int        My,Ny;
     24        int        Mz,Nz;
    2525        IssmDouble *spcvx          = NULL;
    2626        IssmDouble *spcvy          = NULL;
    2727        IssmDouble *spcvz          = NULL;
    … …  
    328328
    329329                                if(!xIsNan<IssmDouble>(spcvx[v1]) && !xIsNan<IssmDouble>(spcvx[v2])){
    330330                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
    331                                                                         1,(spcvx[v1]+spcvx[v2])/yts,DiagnosticHorizAnalysisEnum));
     331                                                                        1,(spcvx[v1]+spcvx[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
    332332                                        count++;
    333333                                }
    334334                                if(!xIsNan<IssmDouble>(spcvy[v1]) && !xIsNan<IssmDouble>(spcvy[v2])){
    335335                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
    336                                                                         2,(spcvy[v1]+spcvy[v2])/yts,DiagnosticHorizAnalysisEnum));
     336                                                                        2,(spcvy[v1]+spcvy[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
    337337                                        count++;
    338338                                }
    339339                                if (reCast<int,IssmDouble>(vertices_type[v1])==StokesApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[v1])==NoneApproximationEnum)){
    340340                                        if(!xIsNan<IssmDouble>(spcvz[v1]) && !xIsNan<IssmDouble>(spcvz[v2])){
    341341                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
    342                                                                                 3,(spcvz[v1]+spcvz[v2])/yts,DiagnosticHorizAnalysisEnum));
     342                                                                                3,(spcvz[v1]+spcvz[v2])/(2.*yts),DiagnosticHorizAnalysisEnum));
    343343                                                count++;
    344344                                        }
    345345                                }

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp

     
    1515
    1616        /*Do not add constraints in DG*/
    1717        if(stabilization!=3){
    18                 IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum);
     18                IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum);
    1919        }
    2020
    2121        /*Assign output pointer: */

  • ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp

     
    2323        iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    2424        if(!isefficientlayer) return;
    2525
    26         IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum);
     26        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum);
    2727
    2828        /*Assign output pointer: */
    2929        *pconstraints=constraints;

  • ../trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp

     
    1818        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
    1919        if(hydrology_model!=HydrologydcEnum) return;
    2020
    21         IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum);
     21        IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum);
    2222
    2323        /*Assign output pointer: */
    2424        *pconstraints=constraints;

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp

     
    1515
    1616        /*Only 3d mesh supported*/
    1717        if(iomodel->dim==3){
    18                 IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum);
     18                IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum);
    1919        }
    2020
    2121        /*Assign output pointer: */

  • ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h

     
    77#include "../../classes/classes.h"
    88
    99/* local prototypes: */
    10 void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type);
     10void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element);
    1111
    1212#endif  /* _IOMODELTOELEMENTINPUTX_H */

  • ../trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

     
    55#include "./IoModelToConstraintsx.h"
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
     8#include "../ModelProcessorx/ModelProcessorx.h"
    89
    9 void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){
     10void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element){
    1011
    1112        /*intermediary: */
    12         int     i,j;
    13         IssmDouble yts;
    14         bool        transient        = false;
     13        int         i,j,count;
     14        IssmDouble  yts;
    1515        FILE       *fid              = NULL;
    16         int         code             = 0;
    17         int         vector_layout    = 0;
     16        int         code,vector_layout;
    1817        IssmDouble *times            = NULL;
    1918        IssmDouble *values           = NULL;
    2019        bool        spcpresent       = false;
    21         int         count            = 0;
    2220
     21        /*P2 finite elements*/
     22        int   v1,v2;
     23        bool *my_edges = NULL;
     24
    2325        /*variables being fetched: */
    24         IssmDouble *IssmDoublevector = NULL;
    25         int     M,N;
     26        IssmDouble *spcdata = NULL;
     27        int         M,N;
    2628
    2729        /*Fetch parameters: */
    2830        iomodel->Constant(&yts,ConstantsYtsEnum);
    2931
    3032        /*First of, find the record for the enum, and get code  of data type: */
    3133        fid=iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum);
    32 
    3334        if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum));
    3435        if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum));
    3536
    3637        /*Fetch vector:*/
    37         iomodel->FetchData(&IssmDoublevector,&M,&N,vector_enum);
     38        iomodel->FetchData(&spcdata,&M,&N,vector_enum);
    3839
    39         /*Transient or static?:*/
     40        /*Partition edges if we are using P2 finite elements*/
     41        if(finite_element==P2Enum){
     42                EdgesPartitioning(&my_edges,iomodel);
     43        }
     44
    4045        if(M==iomodel->numberofvertices){
    41                 /*static: just create Constraints objects*/
     46                /*Static constraint*/
    4247                count=0;
    43 
    44                 /*Create Constraints from x,y,z: */
    4548                for (i=0;i<iomodel->numberofvertices;i++){
    46 
    47                         /*keep only this partition's nodes:*/
    4849                        if((iomodel->my_vertices[i])){
    49 
    50                                 if (!xIsNan<IssmDouble>(IssmDoublevector[i])){
    51 
    52                                         constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type));
     50                                if (!xIsNan<IssmDouble>(spcdata[i])){
     51                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcdata[i],analysis_type));
    5352                                        count++;
    5453                                }
    5554                        }
    5655                }
     56                if(finite_element==P2Enum){
     57                        for(i=0;i<iomodel->numberofedges;i++){
     58                                if(my_edges[i]){
     59                                        v1 = iomodel->edges[4*i+0]-1;
     60                                        v2 = iomodel->edges[4*i+1]-1;
     61                                        if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     62                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
     63                                                                                1,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     64                                                count++;
     65                                        }
     66                                }
     67                        }
     68                }
    5769        }
    5870        else if (M==(iomodel->numberofvertices+1)){
    5971                /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
    … …  
    6274
    6375                /*figure out times: */
    6476                times=xNew<IssmDouble>(N);
    65                 for(j=0;j<N;j++){
    66                         times[j]=IssmDoublevector[(M-1)*N+j]*yts;
    67                 }
     77                for(j=0;j<N;j++) times[j]=spcdata[(M-1)*N+j]*yts;
    6878
    69                 /*Create constraints from x,y,z: */
     79                /*Create constraints: */
    7080                for (i=0;i<iomodel->numberofvertices;i++){
    71 
    72                         /*keep only this partition's nodes:*/
    7381                        if((iomodel->my_vertices[i])){
    7482
    7583                                /*figure out times and values: */
    7684                                values=xNew<IssmDouble>(N);
    7785                                spcpresent=false;
    7886                                for(j=0;j<N;j++){
    79                                         values[j]=IssmDoublevector[i*N+j];
     87                                        values[j]=spcdata[i*N+j];
    8088                                        if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
    8189                                }
    8290
    … …  
    8795                                xDelete<IssmDouble>(values);
    8896                        }
    8997                }
     98                if(finite_element==P2Enum){
     99                        for(i=0;i<iomodel->numberofedges;i++){
     100                                if(my_edges[i]){
     101                                        v1 = iomodel->edges[4*i+0]-1;
     102                                        v2 = iomodel->edges[4*i+1]-1;
     103                                        values=xNew<IssmDouble>(N);
     104                                        spcpresent=false;
     105                                        for(j=0;j<N;j++){
     106                                                values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
     107                                                if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
     108                                        }
     109                                        if(spcpresent){
     110                                                constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1,
     111                                                                                N,times,values,analysis_type));
     112                                                count++;
     113                                        }
     114                                }
     115                        }
     116                }
    90117        }
    91118        else{
    92119                _error_("Size of field " << EnumToStringx(vector_enum) << " not supported");
    93120        }
    94121
    95122        /*Free ressources:*/
    96         xDelete<IssmDouble>(IssmDoublevector);
     123        xDelete<IssmDouble>(spcdata);
    97124        xDelete<IssmDouble>(times);
    98125        xDelete<IssmDouble>(values);
    99126}
Note: See TracBrowser for help on using the repository browser.