Changeset 16434


Ignore:
Timestamp:
10/16/13 14:32:53 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: preparing Free surfaces flowband model

Location:
issm/trunk-jpl/src/c
Files:
19 edited

Legend:

Unmodified
Added
Removed

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

    r16384 r16434  
    140140/*}}}*/
    141141
     142/*FUNCTION Seg::GetSize{{{*/
     143IssmDouble Seg::GetSize(void){
     144
     145        IssmDouble xyz_list[NUMVERTICES][3];
     146        IssmDouble x1,y1,x2,y2;
     147
     148        /*Get xyz list: */
     149        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     150        x1=xyz_list[0][0]; y1=xyz_list[0][1];
     151        x2=xyz_list[1][0]; y2=xyz_list[1][1];
     152
     153        return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
     154}
     155/*}}}*/
     156/*FUNCTION Seg::CreateKMatrixFreeSurfaceTop {{{*/
     157ElementMatrix* Seg::CreateKMatrixFreeSurfaceTop(void){
     158
     159        /*Intermediaries */
     160        int        stabilization;
     161        IssmDouble Jdet,D_scalar,dt,h;
     162        IssmDouble vx,vel;
     163        IssmDouble xyz_list[NUMVERTICES][3];
     164
     165        /*Fetch number of nodes for this finite element*/
     166        int numnodes = this->NumberofNodes();
     167
     168        /*Initialize Element matrix and vectors*/
     169        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     170        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     171        IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
     172        IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
     173
     174        /*Retrieve all inputs and parameters*/
     175        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     176        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     177        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
     178        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     179        h=this->GetSize();
     180
     181        /* Start  looping on the number of gaussian points: */
     182        GaussSeg *gauss=new GaussSeg(2);
     183        for(int ig=gauss->begin();ig<gauss->end();ig++){
     184
     185                gauss->GaussPoint(ig);
     186
     187                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     188                GetNodalFunctions(basis,gauss);
     189
     190                vx_input->GetInputValue(&vx,gauss);
     191
     192                D_scalar=gauss->weight*Jdet;
     193
     194                TripleMultiply(basis,1,numnodes,1,
     195                                        &D_scalar,1,1,0,
     196                                        basis,1,numnodes,0,
     197                                        &Ke->values[0],1);
     198
     199                GetNodalFunctions(B,gauss);
     200                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
     201
     202                D_scalar=dt*gauss->weight*Jdet*vx;
     203                TripleMultiply(B,1,numnodes,1,
     204                                        &D_scalar,1,1,0,
     205                                        Bprime,1,numnodes,0,
     206                                        &Ke->values[0],1);
     207
     208                if(stabilization==2){
     209                        /*Streamline upwinding*/
     210                        vel=fabs(vx)+1.e-8;
     211                        D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
     212                }
     213                else if(stabilization==1){
     214                        /*SSA*/
     215                        vx_input->GetInputAverage(&vx);
     216                        D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
     217                }
     218                if(stabilization==1 || stabilization==2){
     219                        TripleMultiply(Bprime,1,numnodes,1,
     220                                                &D_scalar,1,1,0,
     221                                                Bprime,1,numnodes,0,
     222                                                &Ke->values[0],1);
     223                }
     224        }
     225
     226        /*Clean up and return*/
     227        xDelete<IssmDouble>(basis);
     228        xDelete<IssmDouble>(B);
     229        xDelete<IssmDouble>(Bprime);
     230        delete gauss;
     231        return Ke;
     232}
     233/*}}}*/
     234/*FUNCTION Seg::CreateKMatrixFreeSurfaceBase {{{*/
     235ElementMatrix* Seg::CreateKMatrixFreeSurfaceBase(void){
     236
     237        /*Intermediaries */
     238        int        stabilization;
     239        IssmDouble Jdet,D_scalar,dt,h;
     240        IssmDouble vx,vel;
     241        IssmDouble xyz_list[NUMVERTICES][3];
     242
     243        /*Fetch number of nodes for this finite element*/
     244        int numnodes = this->NumberofNodes();
     245
     246        /*Initialize Element matrix and vectors*/
     247        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     248        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     249        IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
     250        IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
     251
     252        /*Retrieve all inputs and parameters*/
     253        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     254        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     255        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
     256        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     257        h=this->GetSize();
     258
     259        /* Start  looping on the number of gaussian points: */
     260        GaussSeg *gauss=new GaussSeg(2);
     261        for(int ig=gauss->begin();ig<gauss->end();ig++){
     262
     263                gauss->GaussPoint(ig);
     264
     265                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     266                GetNodalFunctions(basis,gauss);
     267
     268                vx_input->GetInputValue(&vx,gauss);
     269
     270                D_scalar=gauss->weight*Jdet;
     271
     272                TripleMultiply(basis,1,numnodes,1,
     273                                        &D_scalar,1,1,0,
     274                                        basis,1,numnodes,0,
     275                                        &Ke->values[0],1);
     276
     277                GetNodalFunctions(B,gauss);
     278                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
     279
     280                D_scalar=dt*gauss->weight*Jdet*vx;
     281                TripleMultiply(B,1,numnodes,1,
     282                                        &D_scalar,1,1,0,
     283                                        Bprime,1,numnodes,0,
     284                                        &Ke->values[0],1);
     285
     286                if(stabilization==2){
     287                        /*Streamline upwinding*/
     288                        vel=fabs(vx)+1.e-8;
     289                        D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
     290                }
     291                else if(stabilization==1){
     292                        /*SSA*/
     293                        vx_input->GetInputAverage(&vx);
     294                        D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
     295                }
     296                if(stabilization==1 || stabilization==2){
     297                        TripleMultiply(Bprime,1,numnodes,1,
     298                                                &D_scalar,1,1,0,
     299                                                Bprime,1,numnodes,0,
     300                                                &Ke->values[0],1);
     301                }
     302        }
     303
     304        /*Clean up and return*/
     305        xDelete<IssmDouble>(basis);
     306        xDelete<IssmDouble>(B);
     307        xDelete<IssmDouble>(Bprime);
     308        delete gauss;
     309        return Ke;
     310}
     311/*}}}*/
    142312/*FUNCTION Seg::CreateMassMatrix {{{*/
    143313ElementMatrix* Seg::CreateMassMatrix(void){
    … …  
    229399}
    230400/*}}}*/
     401/*FUNCTION Seg::CreatePVectorFreeSurfaceTop {{{*/
     402ElementVector* Seg::CreatePVectorFreeSurfaceTop(void){
     403
     404        /*Intermediaries */
     405        IssmDouble Jdet,dt;
     406        IssmDouble ms,surface,vy;
     407        IssmDouble xyz_list[NUMVERTICES][3];
     408
     409        /*Fetch number of nodes and dof for this finite element*/
     410        int numnodes = this->NumberofNodes();
     411
     412        /*Initialize Element vector and other vectors*/
     413        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     414        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     415
     416        /*Retrieve all inputs and parameters*/
     417        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     418        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     419        Input* vy_input     = inputs->GetInput(VyEnum);                         _assert_(vy_input);
     420        Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
     421        Input* surface_input= inputs->GetInput(SurfaceEnum);                    _assert_(surface_input);
     422
     423        /*Initialize mb_correction to 0, do not forget!:*/
     424        /* Start  looping on the number of gaussian points: */
     425        GaussSeg* gauss=new GaussSeg(2);
     426        for(int ig=gauss->begin();ig<gauss->end();ig++){
     427
     428                gauss->GaussPoint(ig);
     429
     430                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     431                GetNodalFunctions(basis,gauss);
     432
     433                vy_input->GetInputValue(&vy,gauss);
     434                ms_input->GetInputValue(&ms,gauss);
     435                surface_input->GetInputValue(&surface,gauss);
     436
     437                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vy)*basis[i];
     438        }
     439
     440        /*Clean up and return*/
     441        xDelete<IssmDouble>(basis);
     442        delete gauss;
     443        return pe;
     444}
     445/*}}}*/
     446/*FUNCTION Seg::CreatePVectorFreeSurfaceBase {{{*/
     447ElementVector* Seg::CreatePVectorFreeSurfaceBase(void){
     448
     449        /*Intermediaries */
     450        IssmDouble Jdet,dt;
     451        IssmDouble mb,mb_correction,bed,vy;
     452        IssmDouble xyz_list[NUMVERTICES][3];
     453
     454        /*Fetch number of nodes and dof for this finite element*/
     455        int numnodes = this->NumberofNodes();
     456
     457        /*Initialize Element vector and other vectors*/
     458        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     459        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     460
     461        /*Retrieve all inputs and parameters*/
     462        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     463        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     464        Input* vy_input            = inputs->GetInput(VyEnum);                         _assert_(vy_input);
     465        Input* mb_input            = inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     466        Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     467        Input* bed_input           = inputs->GetInput(BedEnum);                        _assert_(bed_input);
     468
     469        /*Initialize mb_correction to 0, do not forget!:*/
     470        /* Start  looping on the number of gaussian points: */
     471        GaussSeg* gauss=new GaussSeg(2);
     472        for(int ig=gauss->begin();ig<gauss->end();ig++){
     473
     474                gauss->GaussPoint(ig);
     475
     476                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     477                GetNodalFunctions(basis,gauss);
     478
     479                vy_input->GetInputValue(&vy,gauss);
     480                mb_input->GetInputValue(&mb,gauss);
     481                bed_input->GetInputValue(&bed,gauss);
     482                if(mb_correction_input)
     483                 mb_correction_input->GetInputValue(&mb_correction,gauss);
     484                else
     485                 mb_correction=0.;
     486
     487                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vy)*basis[i];
     488        }
     489
     490        /*Clean up and return*/
     491        xDelete<IssmDouble>(basis);
     492        delete gauss;
     493        return pe;
     494}
     495/*}}}*/

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

    r16388 r16434  
    192192                /*Seg specific routines:*/
    193193                ElementMatrix* CreateMassMatrix(void);
     194                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
     195                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
     196                ElementVector* CreatePVectorFreeSurfaceTop(void);
     197                ElementVector* CreatePVectorFreeSurfaceBase(void);
     198                IssmDouble     GetSize(void);
    194199};
    195200#endif  /* _SEG_H */

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

    r16386 r16434  
    4848
    4949/*Reference Element numerics*/
     50/*FUNCTION SegRef::GetBprimeMasstransport{{{*/
     51void SegRef::GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussSeg* gauss){
     52        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     53         * For node i, Bi' can be expressed in the actual coordinate system
     54         * by:
     55         *       Bi_prime=[ dN/dx ]
     56         * where N is the finiteelement function for node i.
     57         *
     58         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     59         */
     60
     61        /*Fetch number of nodes for this finite element*/
     62        int numnodes = this->NumberofNodes();
     63
     64        /*Get nodal functions derivatives*/
     65        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     66        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     67
     68        /*Build B': */
     69        for(int i=0;i<numnodes;i++){
     70                Bprime[i] = dbasis[i];
     71        }
     72
     73        /*Clean-up*/
     74        xDelete<IssmDouble>(dbasis);
     75}
     76/*}}}*/
    5077/*FUNCTION SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){{{*/
    5178void SegRef::GetNodalFunctions(IssmDouble* basis,GaussSeg* gauss){

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

    r16382 r16434  
    2222                /*Management*/
    2323                void SetElementType(int type,int type_counter);
     24                void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussSeg* gauss);
    2425                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussSeg* gauss);
    2526                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussSeg* gauss);

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

    r16433 r16434  
    208208        /*retreive parameters: */
    209209        int analysis_type;
     210        int meshtype;
    210211        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    211212
    … …  
    246247                #endif
    247248                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    248                         int meshtype;
    249249                        parameters->FindParam(&meshtype,MeshTypeEnum);
    250250                        if(meshtype==Mesh2DverticalEnum){
    … …  
    258258                case MasstransportAnalysisEnum:
    259259                        return CreateKMatrixMasstransport();
     260                        break;
     261                case FreeSurfaceTopAnalysisEnum:
     262                        parameters->FindParam(&meshtype,MeshTypeEnum);
     263                        switch(meshtype){
     264                                case Mesh2DverticalEnum:
     265                                        return CreateKMatrixFreeSurfaceTop1D();
     266                                case Mesh3DEnum:
     267                                        return CreateKMatrixFreeSurfaceTop();
     268                                default:
     269                                        _error_("Mesh not supported yet");
     270                        }
     271                        break;
     272                case FreeSurfaceBaseAnalysisEnum:
     273                        parameters->FindParam(&meshtype,MeshTypeEnum);
     274                        switch(meshtype){
     275                                case Mesh2DverticalEnum:
     276                                        return CreateKMatrixFreeSurfaceBase1D();
     277                                case Mesh3DEnum:
     278                                        return CreateKMatrixFreeSurfaceBase();
     279                                default:
     280                                        _error_("Mesh not supported yet");
     281                        }
    260282                        break;
    261283                #endif
    … …  
    396418
    397419        /*retrive parameters: */
     420        int meshtype;
    398421        int analysis_type;
    399422        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    436459                #endif
    437460                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    438                         int meshtype;
    439461                        parameters->FindParam(&meshtype,MeshTypeEnum);
    440462                        if(meshtype==Mesh2DverticalEnum){
    … …  
    448470                case MasstransportAnalysisEnum:
    449471                        return CreatePVectorMasstransport();
     472                        break;
     473                case FreeSurfaceTopAnalysisEnum:
     474                        parameters->FindParam(&meshtype,MeshTypeEnum);
     475                        switch(meshtype){
     476                                case Mesh2DverticalEnum:
     477                                        return CreatePVectorFreeSurfaceTop1D();
     478                                case Mesh3DEnum:
     479                                        return CreatePVectorFreeSurfaceTop();
     480                                default:
     481                                        _error_("Mesh not supported yet");
     482                        }
     483                        break;
     484                case FreeSurfaceBaseAnalysisEnum:
     485                        parameters->FindParam(&meshtype,MeshTypeEnum);
     486                        switch(meshtype){
     487                                case Mesh2DverticalEnum:
     488                                        return CreatePVectorFreeSurfaceBase1D();
     489                                case Mesh3DEnum:
     490                                        return CreatePVectorFreeSurfaceBase();
     491                                default:
     492                                        _error_("Mesh not supported yet");
     493                        }
    450494                        break;
    451495                #endif
    … …  
    833877IssmDouble Tria::GetArea(void){
    834878
    835         IssmDouble area=0;
    836879        IssmDouble xyz_list[NUMVERTICES][3];
    837880        IssmDouble x1,y1,x2,y2,x3,y3;
    … …  
    18221865                        InputUpdateFromSolutionMasstransport(solution);
    18231866                        break;
     1867                case FreeSurfaceTopAnalysisEnum:
     1868                        InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
     1869                        break;
     1870                case FreeSurfaceBaseAnalysisEnum:
     1871                        InputUpdateFromSolutionOneDof(solution,BedEnum);
     1872                        break;
    18241873                default:
    18251874                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    … …  
    21632212}
    21642213/*}}}*/
     2214/*FUNCTION Tria::HasEdgeOnSurface {{{*/
     2215bool Tria::HasEdgeOnSurface(){
     2216
     2217        IssmDouble values[NUMVERTICES];
     2218        IssmDouble sum;
     2219
     2220        /*Retrieve all inputs and parameters*/
     2221        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     2222        sum = values[0]+values[1]+values[2];
     2223
     2224        _assert_(sum==0. || sum==1. || sum==2.);
     2225
     2226        if(sum==3.)  _error_("Two edges on surface not supported yet...");
     2227
     2228        if(sum>1.){
     2229                return true;
     2230        }
     2231        else{
     2232                return false;
     2233        }
     2234}
     2235/*}}}*/
    21652236/*FUNCTION Tria::EdgeOnBedIndices{{{*/
    21662237void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
    … …  
    21842255}
    21852256/*}}}*/
     2257/*FUNCTION Tria::EdgeOnSurfaceIndices{{{*/
     2258void Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){
     2259
     2260        IssmDouble values[NUMVERTICES];
     2261        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     2262
     2263        /*Retrieve all inputs and parameters*/
     2264        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     2265
     2266        for(int i=0;i<3;i++){
     2267                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     2268                        *pindex1 = indices[i][0];
     2269                        *pindex2 = indices[i][1];
     2270                        return;
     2271                }
     2272        }
     2273
     2274        _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     2275        _error_("Could not find 2 vertices on surface");
     2276}
     2277/*}}}*/
    21862278/*FUNCTION Tria::EdgeOnBedIndex{{{*/
    21872279int Tria::EdgeOnBedIndex(void){
    … …  
    22012293        _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
    22022294        _error_("Could not find 2 vertices on bed");
     2295}
     2296/*}}}*/
     2297/*FUNCTION Tria::EdgeOnSurfaceIndex{{{*/
     2298int Tria::EdgeOnSurfaceIndex(void){
     2299
     2300        IssmDouble values[NUMVERTICES];
     2301        int        indices[3][2] = {{1,2},{2,0},{0,1}};
     2302
     2303        /*Retrieve all inputs and parameters*/
     2304        GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
     2305
     2306        for(int i=0;i<3;i++){
     2307                if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
     2308                        return i;
     2309                }
     2310        }
     2311
     2312        _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
     2313        _error_("Could not find 2 vertices on surface");
    22032314}
    22042315/*}}}*/
    … …  
    74507561}
    74517562/*}}}*/
     7563/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop1D {{{*/
     7564ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop1D(void){
     7565
     7566        if(!HasEdgeOnSurface()) return NULL;
     7567
     7568        int index1,index2;
     7569        this->EdgeOnSurfaceIndices(&index1,&index2);
     7570
     7571        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     7572        ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceTop();
     7573        delete seg->material; delete seg;
     7574
     7575        /*clean up and return*/
     7576        return Ke;
     7577}
     7578/*}}}*/
    74527579/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/
    74537580ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){
    … …  
    75477674}
    75487675/*}}}*/
     7676/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase1D {{{*/
     7677ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase1D(void){
     7678
     7679        if(!HasEdgeOnBed()) return NULL;
     7680
     7681        int index1,index2;
     7682        this->EdgeOnBedIndices(&index1,&index2);
     7683
     7684        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     7685        ElementMatrix* Ke=seg->CreateKMatrixFreeSurfaceBase();
     7686        delete seg->material; delete seg;
     7687
     7688        /*clean up and return*/
     7689        return Ke;
     7690}
     7691/*}}}*/
    75497692/*FUNCTION Tria::CreatePVectorMasstransport{{{*/
    75507693ElementVector* Tria::CreatePVectorMasstransport(void){
    … …  
    76997842}
    77007843/*}}}*/
     7844/*FUNCTION Tria::CreatePVectorFreeSurfaceTop1D {{{*/
     7845ElementVector* Tria::CreatePVectorFreeSurfaceTop1D(void){
     7846
     7847        if(!HasEdgeOnSurface()) return NULL;
     7848
     7849        int index1,index2;
     7850        this->EdgeOnSurfaceIndices(&index1,&index2);
     7851
     7852        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     7853        ElementVector* pe=seg->CreatePVectorFreeSurfaceTop();
     7854        delete seg->material; delete seg;
     7855
     7856        /*clean up and return*/
     7857        return pe;
     7858}
     7859/*}}}*/
    77017860/*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/
    77027861ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){
    … …  
    77467905        xDelete<IssmDouble>(basis);
    77477906        delete gauss;
     7907        return pe;
     7908}
     7909/*}}}*/
     7910/*FUNCTION Tria::CreatePVectorFreeSurfaceBase1D {{{*/
     7911ElementVector* Tria::CreatePVectorFreeSurfaceBase1D(void){
     7912
     7913        if(!HasEdgeOnBed()) return NULL;
     7914
     7915        int index1,index2;
     7916        this->EdgeOnBedIndices(&index1,&index2);
     7917
     7918        Seg* seg=(Seg*)SpawnSeg(index1,index2);
     7919        ElementVector* pe=seg->CreatePVectorFreeSurfaceBase();
     7920        delete seg->material; delete seg;
     7921
     7922        /*clean up and return*/
    77487923        return pe;
    77497924}

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

    r16388 r16434  
    8686                bool        IsOnBed();
    8787                bool        HasEdgeOnBed();
     88                bool        HasEdgeOnSurface();
     89                void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
    8890                void        EdgeOnBedIndices(int* pindex1,int* pindex);
    8991                int         EdgeOnBedIndex();
     92                int         EdgeOnSurfaceIndex();
    9093                bool        IsFloating();
    9194                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
    … …  
    191194                ElementMatrix* CreateKMatrixMasstransport_DG(void);
    192195                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
     196                ElementMatrix* CreateKMatrixFreeSurfaceTop1D(void);
    193197                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
     198                ElementMatrix* CreateKMatrixFreeSurfaceBase1D(void);
    194199                ElementMatrix* CreateMassMatrix(void);
    195200                ElementMatrix* CreateBasalMassMatrix(void);
    … …  
    205210                ElementVector* CreatePVectorMasstransport_DG(void);
    206211                ElementVector* CreatePVectorFreeSurfaceTop(void);
     212                ElementVector* CreatePVectorFreeSurfaceTop1D(void);
    207213                ElementVector* CreatePVectorFreeSurfaceBase(void);
     214                ElementVector* CreatePVectorFreeSurfaceBase1D(void);
    208215                ElementVector* CreatePVectorSlope(void);
    209216                ElementVector* CreateBasalPVectorSlope(void);

  • issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h

    r16382 r16434  
    4747                void GetInputValue(int* pvalue);
    4848                void GetInputValue(IssmDouble* pvalue);
     49                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    4950                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    5051                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r16382 r16434  
    5252                void GetInputValue(int* pvalue);
    5353                void GetInputValue(IssmDouble* pvalue);
     54                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    5455                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    5556                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h

    r16382 r16434  
    4949                void GetInputValue(int* pvalue){_error_("not implemented yet");};
    5050                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
     51                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");};
    5253                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h

    r16382 r16434  
    4646                void GetInputValue(int* pvalue);
    4747                void GetInputValue(IssmDouble* pvalue);
     48                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    4849                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    4950                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/Input.h

    r16382 r16434  
    2929                virtual void GetInputValue(int* pvalue)=0;
    3030                virtual void GetInputValue(IssmDouble* pvalue)=0;
     31                virtual void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss)=0;
    3132                virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss)=0;
    3233                virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0;

  • issm/trunk-jpl/src/c/classes/Inputs/IntInput.h

    r16382 r16434  
    4747                void GetInputValue(int* pvalue);
    4848                void GetInputValue(IssmDouble* pvalue);
     49                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    4950                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    5051                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r16382 r16434  
    4646                void GetInputValue(int* pvalue){_error_("not implemented yet");};
    4747                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
     48                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    4849                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");};
    4950                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp

    r16382 r16434  
    8585
    8686/*Object functions*/
     87/*FUNCTION SegInput::GetInputAverage{{{*/
     88void SegInput::GetInputAverage(IssmDouble* pvalue){
     89
     90        int        numnodes  = this->NumberofNodes();
     91        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
     92        IssmDouble value     = 0.;
     93
     94        for(int i=0;i<numnodes;i++) value+=values[i];
     95        value = value/numnodesd;
     96
     97        *pvalue=value;
     98}
     99/*}}}*/
    87100/*FUNCTION SegInput::GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){{{*/
    88101void SegInput::GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.h

    r16382 r16434  
    5757                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
    5858                void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
    59                 void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");};
     59                void GetInputAverage(IssmDouble* pvalue);
    6060                void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
    6161                void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

    r16382 r16434  
    5050                void GetInputValue(int* pvalue){_error_("not implemented yet");};
    5151                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");};
     52                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    5253                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    5354                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r16382 r16434  
    9797}
    9898/*}}}*/
    99 /*FUNCTION SegInput::SpawnSegInput{{{*/
     99/*FUNCTION TriaInput::SpawnSegInput{{{*/
    100100Input* TriaInput::SpawnSegInput(int index1,int index2){
    101101

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r16382 r16434  
    4747                void GetInputValue(int* pvalue){_error_("not implemented yet");}
    4848                void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}
     49                void GetInputValue(IssmDouble* pvalue,GaussSeg* gauss){_error_("not implemented yet");};
    4950                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
    5051                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error_("not implemented yet");};

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

    r16432 r16434  
    3030        iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.);
    3131        iomodel->FetchDataToInput(elements,VxEnum);
    32         iomodel->FetchDataToInput(elements,VyEnum);
     32        iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
    3333        if(iomodel->meshtype==Mesh3DEnum){
    3434                iomodel->FetchDataToInput(elements,VzEnum);
Note: See TracChangeset for help on using the changeset viewer.