Changeset 15401


Ignore:
Timestamp:
07/02/13 11:43:54 (12 years ago)
Author:
seroussi
Message:

NEW: start to work on dynamic ice front

Location:
issm/trunk-jpl/src
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15387 r15401  
    21782178                                name==FrictionCoefficientEnum ||
    21792179                                name==GLlevelsetEnum ||
     2180                                name==IcelevelsetEnum ||
    21802181                                name==GradientEnum ||
    21812182                                name==OldGradientEnum  ||

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15387 r15401  
    17331733                                name==BedEnum ||
    17341734                                name==GLlevelsetEnum ||
     1735                                name==IcelevelsetEnum ||
    17351736                                name==SurfaceSlopeXEnum ||
    17361737                                name==SurfaceSlopeYEnum ||
    … …  
    29422943ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
    29432944
     2945        /*compute all load vectors for this element*/
     2946        ElementVector* pe1=CreatePVectorDiagnosticMacAyealDrivingStress();
     2947        ElementVector* pe2=CreatePVectorDiagnosticMacAyealFront();
     2948        ElementVector* pe =new ElementVector(pe1,pe2);
     2949
     2950        /*clean-up and return*/
     2951        delete pe1;
     2952        delete pe2;
     2953        return pe;
     2954}
     2955/*}}}*/
     2956/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealDrivingStress {{{*/
     2957ElementVector* Tria::CreatePVectorDiagnosticMacAyealDrivingStress(){
     2958
    29442959        /*Intermediaries */
    29452960        int            i,j;
    … …  
    29482963        IssmDouble     xyz_list[NUMVERTICES][3];
    29492964        IssmDouble     slope[2];
     2965        IssmDouble     icefrontlevel[3];
    29502966
    29512967        /*Fetch number of nodes and dof for this finite element*/
    … …  
    29632979        Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    29642980        Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
     2981        GetInputListOnVertices(&icefrontlevel[0],BedEnum);
    29652982
    29662983        /* Start  looping on the number of gaussian points: */
    … …  
    29812998                                pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    29822999                        }
     3000                }
     3001        }
     3002
     3003        /*Transform coordinate system*/
     3004        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
     3005
     3006        /*Clean up and return*/
     3007        delete gauss;
     3008        xDelete<IssmDouble>(basis);
     3009        return pe;
     3010}
     3011/*}}}*/
     3012/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealFront {{{*/
     3013ElementVector* Tria::CreatePVectorDiagnosticMacAyealFront(){
     3014
     3015        /*Intermediaries */
     3016        int            i,j;
     3017        IssmDouble     ls[3];
     3018        IssmDouble     xyz_list[NUMVERTICES][3];
     3019
     3020        /*Fetch number of nodes and dof for this finite element*/
     3021        int numnodes = this->NumberofNodes();
     3022        int numdof   = numnodes*NDOF2;
     3023        Icefront *icefront=NULL;
     3024
     3025        return NULL;
     3026        /*Retrieve all inputs and parameters*/
     3027        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3028        GetInputListOnVertices(&ls[0],IcelevelsetEnum);
     3029        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     3030        GaussTria*     gauss  = new GaussTria(2);
     3031        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     3032
     3033        /*Create Ice Front if necessary*/
     3034        if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
     3035                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
     3036                        //icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type);
    29833037                }
    29843038        }

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r15387 r15401  
    224224                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
    225225                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
     226                ElementVector* CreatePVectorDiagnosticMacAyealDrivingStress(void);
     227                ElementVector* CreatePVectorDiagnosticMacAyealFront(void);
    226228                ElementVector* CreatePVectorDiagnosticHutter(void);
    227229                ElementMatrix* CreateJacobianDiagnosticMacayeal(void);

  • issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp

    r15298 r15401  
    122122        this->element    = NULL;
    123123        this->matpar     = NULL;
     124}
     125
     126/*}}}*/
     127/*FUNCTION Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type) {{{*/
     128Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in,int in_icefront_type,  int in_analysis_type){
     129
     130        int segment_width;
     131        int element;
     132        int numnodes;
     133        int numvertices;
     134        int dim;
     135        int numberofelements;
     136
     137        /*icefront constructor data: */
     138        int  icefront_eid;
     139        int  icefront_mparid;
     140        int  icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
     141        int  icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
     142
     143//      /*find parameters: */
     144//      iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
     145//
     146        /*First, retrieve element index and element type: */
     147        if(strcmp(element_type_in,"2d")==0){
     148                segment_width=4;
     149        }
     150        else{
     151                segment_width=6;
     152        }
     153//      element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)
     154//
     155//      /*Build ids for hook constructors: */
     156//      icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
     157//      icefront_mparid=numberofelements+1; //matlab indexing
     158//
     159        if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
     160//              icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
     161//              icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
     162//              icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
     163//              icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
     164        }
     165        else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
     166//              icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
     167//              icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
     168//              icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
     169//              icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
     170//              icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
     171//              icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
     172//              icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
     173//              icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
     174        }
     175        else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");
     176
     177        if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
     178                numnodes=4;
     179                numvertices=4;
     180        }
     181        else{
     182                numnodes=2;
     183                numvertices=2;
     184        }
     185
     186        /*Ok, we have everything to build the object: */
     187        this->id=1;
     188        this->analysis_type=in_analysis_type;
     189
     190        /*Hooks: */
     191        this->hnodes=new Hook(icefront_node_ids,numnodes);
     192        this->hvertices=new Hook(icefront_vertex_ids,numvertices);
     193        this->helement=new Hook(&icefront_eid,1);
     194        this->hmatpar=new Hook(&icefront_mparid,1);
     195
     196        //intialize  and add as many inputs per element as requested:
     197        this->inputs=inputs_in;
     198        this->inputs->AddInput(new IntInput(FillEnum,1)); //We always consider we have water, if above sea level, only air will be applied
     199        this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type));
     200
     201        //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
     202        this->parameters = NULL;
     203        this->nodes      = NULL;
     204        this->vertices   = NULL;
     205        this->element    = NULL;
     206        this->matpar     = matpar_in;
    124207}
    125208

  • issm/trunk-jpl/src/c/classes/Loads/Icefront.h

    r14975 r15401  
    4545                Icefront();
    4646                Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int analysis_type);
     47                Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type);
    4748                ~Icefront();
    4849                /*}}}*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r15375 r15401  
    9292
    9393        /*Fetch data:*/
    94         iomodel->FetchData(6,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum);
     94        iomodel->FetchData(7,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum,IcelevelsetEnum);
    9595        CreateNumberNodeToElementConnectivity(iomodel);
    9696

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r15375 r15401  
    144144        MaskVertexongroundediceEnum,
    145145        MaskVertexonwaterEnum,
     146        IcelevelsetEnum,
    146147        MaterialsBetaEnum,
    147148        MaterialsHeatcapacityEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r15375 r15401  
    152152                case MaskVertexongroundediceEnum : return "MaskVertexongroundedice";
    153153                case MaskVertexonwaterEnum : return "MaskVertexonwater";
     154                case IcelevelsetEnum : return "Icelevelset";
    154155                case MaterialsBetaEnum : return "MaterialsBeta";
    155156                case MaterialsHeatcapacityEnum : return "MaterialsHeatcapacity";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r15375 r15401  
    155155              else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum;
    156156              else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum;
     157              else if (strcmp(name,"Icelevelset")==0) return IcelevelsetEnum;
    157158              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    158159              else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
    … …  
    259260              else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum;
    260261              else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum;
    261               else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
     265              if (strcmp(name,"Thickness")==0) return ThicknessEnum;
     266              else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
    266267              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    267268              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    … …  
    382383              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    383384              else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
    384               else if (strcmp(name,"Segment")==0) return SegmentEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
     388              if (strcmp(name,"Segment")==0) return SegmentEnum;
     389              else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    389390              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
    390391              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    … …  
    505506              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    506507              else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
    507               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     511              if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
     512              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    512513              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    513514              else if (strcmp(name,"J")==0) return JEnum;

  • issm/trunk-jpl/src/m/classes/mask.m

    r15131 r15401  
    5454                        WriteData(fid,'object',obj,'fieldname','vertexongroundedice','format','DoubleMat','mattype',1);
    5555                        WriteData(fid,'object',obj,'fieldname','vertexonwater','format','DoubleMat','mattype',1);
     56                        icelevelset=ones(md.mesh.numberofvertices,1);
     57                        pos=md.diagnostic.icefront(:,1:end-1);
     58                        icelevelset(pos(:))=0;
     59                        WriteData(fid,'data',icelevelset,'format','DoubleMat','mattype',1,'enum',IcelevelsetEnum());
    5660                end % }}}
    5761        end

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r15376 r15401  
    18991899        return StringToEnum('MaskVertexonwater')[0]
    19001900
     1901def IcelevelsetEnum():
     1902        """
     1903        ICELEVELSETENUM - Enum of Icelevelset
     1904
     1905        WARNING: DO NOT MODIFY THIS FILE
     1906                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1907                                Please read src/c/shared/Enum/README for more information
     1908
     1909           Usage:
     1910              macro=IcelevelsetEnum()
     1911        """
     1912
     1913        return StringToEnum('Icelevelset')[0]
     1914
    19011915def MaterialsBetaEnum():
    19021916        """
    … …  
    78617875        """
    78627876
    7863         return 560
    7864 
     7877        return 561
     7878

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r15376 r15401  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=560;
     11macro=561;
Note: See TracChangeset for help on using the changeset viewer.