Changeset 17808


Ignore:
Timestamp:
04/22/14 10:51:09 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: back to normal

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

Legend:

Unmodified
Added
Removed

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

    r17801 r17808  
    8989ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
    9090
    91         _error_("NOT IMPLEMENTED YET");
    92 
     91        if(!element->IsOnBase()) return NULL;
     92        Element* basalelement = element->SpawnBasalElement();
     93
     94        /*Intermediaries */
     95        int  dim, domaintype;
     96        int i, row, col;
     97        IssmDouble kappa;
     98        IssmDouble Jdet, dt, D_scalar;
     99        IssmDouble h,hx,hy,hz;
     100        IssmDouble vel, calvingrate;
     101        IssmDouble norm_dlsf;
     102        IssmDouble* xyz_list = NULL;
     103
     104        /*Get problem dimension*/
     105        basalelement->FindParam(&domaintype,DomainTypeEnum);
     106        switch(domaintype){
     107                case Domain2DverticalEnum:   dim = 1; break;
     108                case Domain2DhorizontalEnum: dim = 2; break;
     109                case Domain3DEnum:           dim = 2; break;
     110                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     111        }
     112
     113        /*Fetch number of nodes and dof for this finite element*/
     114        int numnodes    = basalelement->GetNumberOfNodes();
     115
     116        /*Initialize Element vector and other vectors*/
     117        ElementMatrix* Ke       = basalelement->NewElementMatrix();
     118        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
     119        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
     120        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
     121        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
     122        IssmDouble*    v        = xNew<IssmDouble>(dim);
     123        IssmDouble*    w        = xNew<IssmDouble>(dim);
     124        IssmDouble*    c        = xNew<IssmDouble>(dim);
     125        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
     126
     127        /*Retrieve all inputs and parameters*/
     128        basalelement->GetVerticesCoordinates(&xyz_list);
     129        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     130        Input* vx_input=NULL;
     131        Input* vy_input=NULL;
     132        if(domaintype==Domain2DhorizontalEnum){
     133                vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     134                vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     135        }
     136        else{
     137                if(dim==1){
     138                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     139                }
     140                if(dim==2){
     141                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
     142                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
     143                }
     144        }
     145
     146        Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     147        Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     148        Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
     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                D_scalar=gauss->weight*Jdet;
     157
     158                /* Transient */
     159                if(dt!=0.){
     160                        basalelement->NodalFunctions(basis,gauss);
     161                        TripleMultiply(basis,numnodes,1,0,
     162                                                &D_scalar,1,1,0,
     163                                                basis,1,numnodes,0,
     164                                                &Ke->values[0],1);
     165                        D_scalar*=dt;
     166                }
     167
     168                /* Advection */
     169                GetB(B,basalelement,xyz_list,gauss);
     170                GetBprime(Bprime,basalelement,xyz_list,gauss);
     171                vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
     172                vy_input->GetInputValue(&v[1],gauss);
     173                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     174                lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     175                calvingrate_input->GetInputValue(&calvingrate,gauss);
     176
     177                norm_dlsf=0.;
     178                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     179                norm_dlsf=sqrt(norm_dlsf);
     180
     181                if(norm_dlsf>1.e-10)
     182                        for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
     183                else
     184                        for(i=0;i<dim;i++) c[i]=0.;
     185               
     186                for(i=0;i<dim;i++) w[i]=v[i]-c[i];
     187
     188                for(row=0;row<dim;row++)
     189                        for(col=0;col<dim;col++)
     190                                if(row==col)
     191                                        D[row*dim+col]=D_scalar*w[row];
     192                                else
     193                                   D[row*dim+col]=0.;
     194
     195                TripleMultiply(B,dim,numnodes,1,
     196                                        D,dim,dim,0,
     197                                        Bprime,dim,numnodes,0,
     198                                        &Ke->values[0],1);
     199
     200                /* Stabilization */
     201                int stabilization=2;
     202                vel=0.;
     203                for(i=0;i<dim;i++) vel+=v[i]*v[i];
     204                vel=sqrt(vel)+1.e-14;
     205                switch(stabilization){
     206                        case 0:
     207                                // no stabilization, do nothing
     208                                break;
     209                        case 1:
     210                                /* Artificial Diffusion */
     211                                basalelement->ElementSizes(&hx,&hy,&hz);
     212                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
     213                                kappa=h*vel/2.;
     214                                for(row=0;row<dim;row++)
     215                                        for(col=0;col<dim;col++)
     216                                        if(row==col)
     217                                                D[row*dim+col]=D_scalar*kappa;
     218                                        else
     219                                                D[row*dim+col]=0.;
     220
     221                                TripleMultiply(Bprime,dim,numnodes,1,
     222                                                        D,dim,dim,0,
     223                                                        Bprime,dim,numnodes,0,
     224                                                        &Ke->values[0],1);
     225                                break; 
     226                        case 2:
     227                                /* Streamline Upwinding */
     228                                basalelement->ElementSizes(&hx,&hy,&hz);
     229                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
     230                                for(row=0;row<dim;row++)
     231                                        for(col=0;col<dim;col++)
     232                                                D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
     233
     234                                TripleMultiply(Bprime,dim,numnodes,1,
     235                                                        D,dim,dim,0,
     236                                                        Bprime,dim,numnodes,0,
     237                                                        &Ke->values[0],1);
     238                                break;
     239                        default:
     240                                _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
     241                }
     242        }
     243
     244        /*Clean up and return*/
     245        xDelete<IssmDouble>(xyz_list);
     246        xDelete<IssmDouble>(basis);
     247        xDelete<IssmDouble>(B);
     248        xDelete<IssmDouble>(D);
     249        xDelete<IssmDouble>(Bprime);
     250        xDelete<IssmDouble>(v);
     251        xDelete<IssmDouble>(w);
     252        xDelete<IssmDouble>(c);
     253        xDelete<IssmDouble>(dlsf);
     254        delete gauss;
     255        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     256        return Ke;
    93257}/*}}}*/
    94258ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.