Changeset 17172


Ignore:
Timestamp:
01/26/14 13:35:28 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: trying to reimplement free surfaces (testing soon)

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

Legend:

Unmodified
Added
Removed

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

    r17005 r17172  
    101101}/*}}}*/
    102102ElementMatrix* FreeSurfaceBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
    103         _error_("not implemented yet");
     103
     104        /*Intermediaries*/
     105        int         meshtype,dim,stabilization;
     106        Element*    basalelement = NULL;
     107        IssmDouble *xyz_list  = NULL;
     108        IssmDouble  Jdet,D_scalar,dt,h;
     109        IssmDouble  vel,vx,vy;
     110
     111        /*Get basal element*/
     112        element->FindParam(&meshtype,MeshTypeEnum);
     113        switch(meshtype){
     114                case Mesh2DhorizontalEnum:
     115                        basalelement = element;
     116                        dim = 2;
     117                        break;
     118                case Mesh2DverticalEnum:
     119                        if(!element->IsOnBed()) return NULL;
     120                        basalelement = element->SpawnBasalElement();
     121                        dim = 1;
     122                        break;
     123                case Mesh3DEnum:
     124                        if(!element->IsOnBed()) return NULL;
     125                        basalelement = element->SpawnBasalElement();
     126                        dim = 2;
     127                        break;
     128                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     129        }
     130
     131        /*Fetch number of nodes and dof for this finite element*/
     132        int numnodes = basalelement->GetNumberOfNodes();
     133
     134        /*Initialize Element vector*/
     135        ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
     136        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     137        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
     138        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     139        IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
     140
     141        /*Retrieve all inputs and parameters*/
     142        basalelement->GetVerticesCoordinates(&xyz_list);
     143        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     144        basalelement->FindParam(&stabilization,MasstransportStabilizationEnum);
     145        Input* vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     146        Input* vy_input=NULL;
     147        if(dim>1) basalelement->GetInput(VyEnum); _assert_(vy_input);
     148        h = basalelement->CharacteristicLength();
     149
     150        /* Start  looping on the number of gaussian points: */
     151        Gauss* gauss=basalelement->NewGauss(2);
     152        for(int ig=gauss->begin();ig<gauss->end();ig++){
     153                gauss->GaussPoint(ig);
     154
     155                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     156                basalelement->NodalFunctions(basis,gauss);
     157
     158                vx_input->GetInputValue(&vx,gauss);
     159                if(dim==2) vy_input->GetInputValue(&vy,gauss);
     160
     161                D_scalar=gauss->weight*Jdet;
     162                TripleMultiply(basis,1,numnodes,1,
     163                                        &D_scalar,1,1,0,
     164                                        basis,1,numnodes,0,
     165                                        &Ke->values[0],1);
     166
     167                GetB(B,element,dim,xyz_list,gauss);
     168                GetBprime(Bprime,element,dim,xyz_list,gauss);
     169
     170                D_scalar=dt*gauss->weight*Jdet;
     171
     172                for(int i=0;i<dim*dim;i++) D[i]=0.;
     173                D[0] = D_scalar*vx;
     174                if(dim==2) D[1*dim+1] = D_scalar*vy;
     175                TripleMultiply(B,dim,numnodes,1,
     176                                        D,dim,dim,0,
     177                                        B,dim,numnodes,0,
     178                                        &Ke->values[0],1);
     179
     180                D[0*dim+0]=D_scalar*vx;
     181                if(dim==2) D[1*dim+1]=D_scalar*vy;
     182                TripleMultiply(B,dim,numnodes,1,
     183                                        D,dim,dim,0,
     184                                        Bprime,dim,numnodes,0,
     185                                        &Ke->values[0],1);
     186
     187                if(stabilization==2){
     188                        /*Streamline upwinding*/
     189                        if(dim==1){
     190                         vel=fabs(vx)+1.e-8;
     191                         D[0] = h/(2.*vel)*vx;
     192                        }
     193                        else{
     194                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
     195                         D[0*dim+0]=h/(2*vel)*vx*vx;
     196                         D[1*dim+0]=h/(2*vel)*vy*vx;
     197                         D[0*dim+1]=h/(2*vel)*vx*vy;
     198                         D[1*dim+1]=h/(2*vel)*vy*vy;
     199                        }
     200                }
     201                else if(stabilization==1){
     202                        /*SSA*/
     203                        if(dim==1){
     204                                vx_input->GetInputAverage(&vx);
     205                                D[0]=h/2.*fabs(vx);
     206                        }
     207                        else{
     208                                vx_input->GetInputAverage(&vx);
     209                                vy_input->GetInputAverage(&vy);
     210                                D[0*dim+0]=h/2.0*fabs(vx);
     211                                D[1*dim+1]=h/2.0*fabs(vy);
     212                        }
     213                }
     214                if(stabilization==1 || stabilization==2){
     215                        for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
     216                        TripleMultiply(Bprime,dim,numnodes,1,
     217                                                D,dim,dim,0,
     218                                                Bprime,dim,numnodes,0,
     219                                                &Ke->values[0],1);
     220                }
     221        }
     222
     223        /*Clean up and return*/
     224        xDelete<IssmDouble>(xyz_list);
     225        xDelete<IssmDouble>(basis);
     226        xDelete<IssmDouble>(B);
     227        xDelete<IssmDouble>(Bprime);
     228        xDelete<IssmDouble>(D);
     229        delete gauss;
     230        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     231        return Ke;
    104232}/*}}}*/
    105233ElementVector* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/
    106 _error_("not implemented yet");
     234        /*Intermediaries*/
     235        int         meshtype,dim;
     236        IssmDouble  Jdet,dt;
     237        IssmDouble  mb,mb_correction,bed,vz;
     238        Element*    basalelement = NULL;
     239        IssmDouble *xyz_list  = NULL;
     240
     241        /*Get basal element*/
     242        element->FindParam(&meshtype,MeshTypeEnum);
     243        switch(meshtype){
     244                case Mesh2DhorizontalEnum:
     245                        basalelement = element;
     246                        dim = 2;
     247                        break;
     248                case Mesh2DverticalEnum:
     249                        if(!element->IsOnBed()) return NULL;
     250                        basalelement = element->SpawnBasalElement();
     251                        dim = 1;
     252                        break;
     253                case Mesh3DEnum:
     254                        if(!element->IsOnBed()) return NULL;
     255                        basalelement = element->SpawnBasalElement();
     256                        dim = 2;
     257                        break;
     258                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     259        }
     260
     261        /*Fetch number of nodes and dof for this finite element*/
     262        int numnodes = basalelement->GetNumberOfNodes();
     263
     264        /*Initialize Element vector and other vectors*/
     265        ElementVector* pe    = basalelement->NewElementVector();
     266        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     267
     268        /*Retrieve all inputs and parameters*/
     269        basalelement->GetVerticesCoordinates(&xyz_list);
     270        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     271        Input* mb_input            = basalelement->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
     272        Input* mb_correction_input = basalelement->GetInput(BasalforcingsMeltingRateCorrectionEnum);
     273        Input* bed_input           = basalelement->GetInput(BedEnum);                        _assert_(bed_input);
     274        Input* vz_input      = NULL;
     275        switch(dim){
     276                case 1: vz_input = basalelement->GetInput(VyEnum); _assert_(vz_input); break;
     277                case 2: vz_input = basalelement->GetInput(VzEnum); _assert_(vz_input); break;
     278                default: _error_("not implemented");
     279        }
     280
     281        /*Initialize mb_correction to 0, do not forget!:*/
     282        /* Start  looping on the number of gaussian points: */
     283        Gauss* gauss=basalelement->NewGauss(2);
     284        for(int ig=gauss->begin();ig<gauss->end();ig++){
     285                gauss->GaussPoint(ig);
     286
     287                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     288                basalelement->NodalFunctions(basis,gauss);
     289
     290                vz_input->GetInputValue(&vz,gauss);
     291                mb_input->GetInputValue(&mb,gauss);
     292                bed_input->GetInputValue(&bed,gauss);
     293                if(mb_correction_input)
     294                 mb_correction_input->GetInputValue(&mb_correction,gauss);
     295                else
     296                 mb_correction=0.;
     297
     298                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
     299        }
     300
     301        /*Clean up and return*/
     302        xDelete<IssmDouble>(xyz_list);
     303        xDelete<IssmDouble>(basis);
     304        delete gauss;
     305        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     306        return pe;
     307
     308}/*}}}*/
     309void FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     310        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     311         * For node i, Bi can be expressed in the actual coordinate system
     312         * by:
     313         *       Bi=[ N ]
     314         *          [ N ]
     315         * where N is the finiteelement function for node i.
     316         *
     317         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
     318         */
     319
     320        /*Fetch number of nodes for this finite element*/
     321        int numnodes = element->GetNumberOfNodes();
     322
     323        /*Get nodal functions*/
     324        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     325        element->NodalFunctions(basis,gauss);
     326
     327        /*Build B: */
     328        for(int i=0;i<numnodes;i++){
     329                for(int j=0;j<dim;i++){
     330                        B[numnodes*j+i] = basis[i];
     331                }
     332        }
     333
     334        /*Clean-up*/
     335        xDelete<IssmDouble>(basis);
     336}/*}}}*/
     337void FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     338        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     339         * For node i, Bi' can be expressed in the actual coordinate system
     340         * by:
     341         *       Bi_prime=[ dN/dx ]
     342         *                [ dN/dy ]
     343         * where N is the finiteelement function for node i.
     344         *
     345         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     346         */
     347
     348        /*Fetch number of nodes for this finite element*/
     349        int numnodes = element->GetNumberOfNodes();
     350
     351        /*Get nodal functions derivatives*/
     352        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     353        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     354
     355        /*Build B': */
     356        for(int i=0;i<numnodes;i++){
     357                for(int j=0;j<dim;j++){
     358                        Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
     359                }
     360        }
     361
     362        /*Clean-up*/
     363        xDelete<IssmDouble>(dbasis);
     364
    107365}/*}}}*/
    108366void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h

    r17005 r17172  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
     28                void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     29                void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    2830                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2931                void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

    r17083 r17172  
    110110}/*}}}*/
    111111ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
    112         _error_("not implemented yet");
     112
     113        /*Intermediaries*/
     114        int         meshtype,dim,stabilization;
     115        Element*    topelement = NULL;
     116        IssmDouble *xyz_list  = NULL;
     117        IssmDouble  Jdet,D_scalar,dt,h;
     118        IssmDouble  vel,vx,vy;
     119
     120        /*Get top element*/
     121        element->FindParam(&meshtype,MeshTypeEnum);
     122        switch(meshtype){
     123                case Mesh2DhorizontalEnum:
     124                        topelement = element;
     125                        dim = 2;
     126                        break;
     127                case Mesh2DverticalEnum:
     128                        if(!element->IsOnSurface()) return NULL;
     129                        topelement = element->SpawnTopElement();
     130                        dim = 1;
     131                        break;
     132                case Mesh3DEnum:
     133                        if(!element->IsOnSurface()) return NULL;
     134                        topelement = element->SpawnTopElement();
     135                        dim = 2;
     136                        break;
     137                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     138        }
     139
     140        /*Fetch number of nodes and dof for this finite element*/
     141        int numnodes = topelement->GetNumberOfNodes();
     142
     143        /*Initialize Element vector*/
     144        ElementMatrix* Ke     = topelement->NewElementMatrix(NoneApproximationEnum);
     145        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     146        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
     147        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     148        IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
     149
     150        /*Retrieve all inputs and parameters*/
     151        topelement->GetVerticesCoordinates(&xyz_list);
     152        topelement->FindParam(&dt,TimesteppingTimeStepEnum);
     153        topelement->FindParam(&stabilization,MasstransportStabilizationEnum);
     154        Input* vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
     155        Input* vy_input=NULL;
     156        if(dim>1) topelement->GetInput(VyEnum); _assert_(vy_input);
     157        h = topelement->CharacteristicLength();
     158
     159        /* Start  looping on the number of gaussian points: */
     160        Gauss* gauss=topelement->NewGauss(2);
     161        for(int ig=gauss->begin();ig<gauss->end();ig++){
     162                gauss->GaussPoint(ig);
     163
     164                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     165                topelement->NodalFunctions(basis,gauss);
     166
     167                vx_input->GetInputValue(&vx,gauss);
     168                if(dim==2) vy_input->GetInputValue(&vy,gauss);
     169
     170                D_scalar=gauss->weight*Jdet;
     171                TripleMultiply(basis,1,numnodes,1,
     172                                        &D_scalar,1,1,0,
     173                                        basis,1,numnodes,0,
     174                                        &Ke->values[0],1);
     175
     176                GetB(B,element,dim,xyz_list,gauss);
     177                GetBprime(Bprime,element,dim,xyz_list,gauss);
     178
     179                D_scalar=dt*gauss->weight*Jdet;
     180
     181                for(int i=0;i<dim*dim;i++) D[i]=0.;
     182                D[0] = D_scalar*vx;
     183                if(dim==2) D[1*dim+1] = D_scalar*vy;
     184                TripleMultiply(B,dim,numnodes,1,
     185                                        D,dim,dim,0,
     186                                        B,dim,numnodes,0,
     187                                        &Ke->values[0],1);
     188
     189                D[0*dim+0]=D_scalar*vx;
     190                if(dim==2) D[1*dim+1]=D_scalar*vy;
     191                TripleMultiply(B,dim,numnodes,1,
     192                                        D,dim,dim,0,
     193                                        Bprime,dim,numnodes,0,
     194                                        &Ke->values[0],1);
     195
     196                if(stabilization==2){
     197                        /*Streamline upwinding*/
     198                        if(dim==1){
     199                                vel=fabs(vx)+1.e-8;
     200                                D[0] = h/(2.*vel)*vx;
     201                        }
     202                        else{
     203                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
     204                                D[0*dim+0]=h/(2*vel)*vx*vx;
     205                                D[1*dim+0]=h/(2*vel)*vy*vx;
     206                                D[0*dim+1]=h/(2*vel)*vx*vy;
     207                                D[1*dim+1]=h/(2*vel)*vy*vy;
     208                        }
     209                }
     210                else if(stabilization==1){
     211                        /*SSA*/
     212                        if(dim==1){
     213                                vx_input->GetInputAverage(&vx);
     214                                D[0]=h/2.*fabs(vx);
     215                        }
     216                        else{
     217                                vx_input->GetInputAverage(&vx);
     218                                vy_input->GetInputAverage(&vy);
     219                                D[0*dim+0]=h/2.0*fabs(vx);
     220                                D[1*dim+1]=h/2.0*fabs(vy);
     221                        }
     222                }
     223                if(stabilization==1 || stabilization==2){
     224                        for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
     225                        TripleMultiply(Bprime,dim,numnodes,1,
     226                                                D,dim,dim,0,
     227                                                Bprime,dim,numnodes,0,
     228                                                &Ke->values[0],1);
     229                }
     230        }
     231
     232        /*Clean up and return*/
     233        xDelete<IssmDouble>(xyz_list);
     234        xDelete<IssmDouble>(basis);
     235        xDelete<IssmDouble>(B);
     236        xDelete<IssmDouble>(Bprime);
     237        xDelete<IssmDouble>(D);
     238        delete gauss;
     239        if(meshtype!=Mesh2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
     240        return Ke;
    113241}/*}}}*/
    114242ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
    115 _error_("not implemented yet");
     243        /*Intermediaries*/
     244        int         meshtype,dim;
     245        IssmDouble  Jdet,dt;
     246        IssmDouble  ms,surface,vz;
     247        Element*    topelement = NULL;
     248        IssmDouble *xyz_list  = NULL;
     249
     250        /*Get top element*/
     251        element->FindParam(&meshtype,MeshTypeEnum);
     252        switch(meshtype){
     253                case Mesh2DhorizontalEnum:
     254                        topelement = element;
     255                        dim = 2;
     256                        break;
     257                case Mesh2DverticalEnum:
     258                        if(!element->IsOnSurface()) return NULL;
     259                        topelement = element->SpawnTopElement();
     260                        dim = 1;
     261                        break;
     262                case Mesh3DEnum:
     263                        if(!element->IsOnSurface()) return NULL;
     264                        topelement = element->SpawnTopElement();
     265                        dim = 2;
     266                        break;
     267                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     268        }
     269
     270        /*Fetch number of nodes and dof for this finite element*/
     271        int numnodes = topelement->GetNumberOfNodes();
     272
     273        /*Initialize Element vector and other vectors*/
     274        ElementVector* pe    = topelement->NewElementVector();
     275        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     276
     277        /*Retrieve all inputs and parameters*/
     278        topelement->GetVerticesCoordinates(&xyz_list);
     279        topelement->FindParam(&dt,TimesteppingTimeStepEnum);
     280        Input* ms_input      = topelement->GetInput(SurfaceforcingsMassBalanceEnum);  _assert_(ms_input);
     281        Input* surface_input = topelement->GetInput(SurfaceEnum);                     _assert_(surface_input);
     282        Input* vz_input      = NULL;
     283        switch(dim){
     284                case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break;
     285                case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break;
     286                default: _error_("not implemented");
     287        }
     288
     289        /*Initialize mb_correction to 0, do not forget!:*/
     290        /* Start  looping on the number of gaussian points: */
     291        Gauss* gauss=topelement->NewGauss(2);
     292        for(int ig=gauss->begin();ig<gauss->end();ig++){
     293                gauss->GaussPoint(ig);
     294
     295                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     296                topelement->NodalFunctions(basis,gauss);
     297
     298                ms_input->GetInputValue(&ms,gauss);
     299                vz_input->GetInputValue(&vz,gauss);
     300                surface_input->GetInputValue(&surface,gauss);
     301
     302                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
     303        }
     304
     305        /*Clean up and return*/
     306        xDelete<IssmDouble>(xyz_list);
     307        xDelete<IssmDouble>(basis);
     308        delete gauss;
     309        if(meshtype!=Mesh2DhorizontalEnum){topelement->DeleteMaterials(); delete topelement;};
     310        return pe;
     311
     312}/*}}}*/
     313void FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     314        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     315         * For node i, Bi can be expressed in the actual coordinate system
     316         * by:
     317         *       Bi=[ N ]
     318         *          [ N ]
     319         * where N is the finiteelement function for node i.
     320         *
     321         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
     322         */
     323
     324        /*Fetch number of nodes for this finite element*/
     325        int numnodes = element->GetNumberOfNodes();
     326
     327        /*Get nodal functions*/
     328        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     329        element->NodalFunctions(basis,gauss);
     330
     331        /*Build B: */
     332        for(int i=0;i<numnodes;i++){
     333                for(int j=0;j<dim;i++){
     334                        B[numnodes*j+i] = basis[i];
     335                }
     336        }
     337
     338        /*Clean-up*/
     339        xDelete<IssmDouble>(basis);
     340}/*}}}*/
     341void FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     342        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     343         * For node i, Bi' can be expressed in the actual coordinate system
     344         * by:
     345         *       Bi_prime=[ dN/dx ]
     346         *                [ dN/dy ]
     347         * where N is the finiteelement function for node i.
     348         *
     349         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     350         */
     351
     352        /*Fetch number of nodes for this finite element*/
     353        int numnodes = element->GetNumberOfNodes();
     354
     355        /*Get nodal functions derivatives*/
     356        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     357        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     358
     359        /*Build B': */
     360        for(int i=0;i<numnodes;i++){
     361                for(int j=0;j<dim;j++){
     362                        Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
     363                }
     364        }
     365
     366        /*Clean-up*/
     367        xDelete<IssmDouble>(dbasis);
     368
    116369}/*}}}*/
    117370void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h

    r17005 r17172  
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
     28                void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     29                void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    2830                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2931                void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

    r17118 r17172  
    208208                virtual void   SmbGradients(void)=0;
    209209           virtual Element*   SpawnBasalElement(void)=0;
     210                virtual Element*   SpawnTopElement(void)=0;
    210211                virtual void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
    211212                virtual void   ResetCoordinateSystem()=0;

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

    r17137 r17172  
    26052605        this->inputs->DeleteInput(VxAverageEnum);
    26062606        this->inputs->DeleteInput(VyAverageEnum);
     2607
     2608        return tria;
     2609}
     2610/*}}}*/
     2611/*FUNCTION Penta::SpawnTopElement{{{*/
     2612Element*  Penta::SpawnTopElement(void){
     2613
     2614        _assert_(this->IsOnSurface());
     2615
     2616        Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    26072617
    26082618        return tria;

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

    r17118 r17172  
    121121                void   SetTemporaryElementType(int element_type_in);
    122122           Element* SpawnBasalElement(void);
     123                Element* SpawnTopElement(void);
    123124                IssmDouble SurfaceArea(void);
    124125                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);

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

    r17118 r17172  
    121121                int         NumberofNodesPressure(void){_error_("not implemented yet");};
    122122           Element*    SpawnBasalElement(void){_error_("not implemented yet");};
     123                Element*    SpawnTopElement(void){_error_("not implemented yet");};
    123124                IssmDouble  StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
    124125                IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};

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

    r17162 r17172  
    22612261}
    22622262/*}}}*/
     2263/*FUNCTION Tria::SpawnTopElement{{{*/
     2264Element*  Tria::SpawnTopElement(void){
     2265
     2266        int index1,index2;
     2267        int meshtype;
     2268
     2269        this->parameters->FindParam(&meshtype,MeshTypeEnum);
     2270        switch(meshtype){
     2271                case Mesh2DhorizontalEnum:
     2272                        return this;
     2273                case Mesh2DverticalEnum:
     2274                        _assert_(HasEdgeOnSurface());
     2275                        this->EdgeOnSurfaceIndices(&index1,&index2);
     2276                        return SpawnSeg(index1,index2);
     2277                default:
     2278                        _error_("not implemented yet");
     2279        }
     2280}
     2281/*}}}*/
    22632282/*FUNCTION Tria::SurfaceArea {{{*/
    22642283IssmDouble Tria::SurfaceArea(void){

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

    r17158 r17172  
    119119                void          SmbGradients();
    120120           Element*    SpawnBasalElement(void);
     121                Element*    SpawnTopElement(void);
    121122                int         VelocityInterpolation();
    122123                IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/shared/Exceptions/Exceptions.cpp

    r17107 r17172  
    4545
    4646        file_line= what_line;
     47        this->Report();
    4748
    4849}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.