Changeset 17801

Timestamp:

04/22/14 08:55:17 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing Levelset code before sync

File:

• issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)

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

@@ -89 +89 @@
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
-        if(!element->IsOnBase()) return NULL;
-        Element* basalelement = element->SpawnBasalElement();
-
-        /*Intermediaries */
-        int  dim, domaintype;
-        int i, row, col;
-        IssmDouble kappa;
-        IssmDouble Jdet, dt, D_scalar;
-        IssmDouble h,hx,hy,hz;
-        IssmDouble vel, calvingrate;
-        IssmDouble norm_dlsf;
-        IssmDouble* xyz_list = NULL;
-
-        /*Get problem dimension*/
-        basalelement->FindParam(&domaintype,DomainTypeEnum);
-        switch(domaintype){
-                case Domain2DverticalEnum:   dim = 1; break;
-                case Domain2DhorizontalEnum: dim = 2; break;
-                case Domain3DEnum:           dim = 2; break;
-                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-        }
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int numnodes    = basalelement->GetNumberOfNodes();
-
-        /*Initialize Element vector and other vectors*/
-        ElementMatrix* Ke       = basalelement->NewElementMatrix();
-        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
-        IssmDouble*    v        = xNew<IssmDouble>(dim);
-        IssmDouble*    w        = xNew<IssmDouble>(dim);
-        IssmDouble*    c        = xNew<IssmDouble>(dim);
-        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
-
-        /*Retrieve all inputs and parameters*/
-        basalelement->GetVerticesCoordinates(&xyz_list);
-        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-        Input* vx_input=NULL;
-        Input* vy_input=NULL;
-        if(domaintype==Domain2DhorizontalEnum){
-                vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-                vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-        }
-        else{
-                if(dim==1){
-                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-                }
-                if(dim==2){
-                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
-                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-                }
-        }
-
-        Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-        Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-        Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
-        
-        /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=basalelement->NewGauss(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-                gauss->GaussPoint(ig);
-
-                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                D_scalar=gauss->weight*Jdet;
-
-                /* Transient */
-                if(dt!=0.){
-                        basalelement->NodalFunctions(basis,gauss);
-                        TripleMultiply(basis,numnodes,1,0,
-                                                &D_scalar,1,1,0,
-                                                basis,1,numnodes,0,
-                                                &Ke->values[0],1);
-                        D_scalar*=dt;
-                }
-
-                /* Advection */
-                GetB(B,basalelement,xyz_list,gauss); 
-                GetBprime(Bprime,basalelement,xyz_list,gauss); 
-                vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
-                vy_input->GetInputValue(&v[1],gauss); 
-                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-                lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-                calvingrate_input->GetInputValue(&calvingrate,gauss);
-
-                norm_dlsf=0.;
-                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-                norm_dlsf=sqrt(norm_dlsf);
-
-                if(norm_dlsf>1.e-10)
-                        for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
-                else
-                        for(i=0;i<dim;i++) c[i]=0.;
-                
-                for(i=0;i<dim;i++) w[i]=v[i]-c[i];
-
-                for(row=0;row<dim;row++)
-                        for(col=0;col<dim;col++)
-                                if(row==col)
-                                        D[row*dim+col]=D_scalar*w[row];
-                                else
-                                   D[row*dim+col]=0.;
-
-                TripleMultiply(B,dim,numnodes,1,
-                                        D,dim,dim,0,
-                                        Bprime,dim,numnodes,0,
-                                        &Ke->values[0],1);
-
-                /* Stabilization */
-                int stabilization=2;
-                vel=0.;
-                for(i=0;i<dim;i++) vel+=v[i]*v[i];
-                vel=sqrt(vel)+1.e-14;
-                switch(stabilization){
-                        case 0:
-                                // no stabilization, do nothing
-                                break;
-                        case 1:
-                                /* Artificial Diffusion */
-                                basalelement->ElementSizes(&hx,&hy,&hz);
-                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
-                                kappa=h*vel/2.;
-                                for(row=0;row<dim;row++)
-                                        for(col=0;col<dim;col++)
-                                        if(row==col)
-                                                D[row*dim+col]=D_scalar*kappa;
-                                        else
-                                                D[row*dim+col]=0.;
-
-                                TripleMultiply(Bprime,dim,numnodes,1,
-                                                        D,dim,dim,0,
-                                                        Bprime,dim,numnodes,0,
-                                                        &Ke->values[0],1);
-                                break;  
-                        case 2:
-                                /* Streamline Upwinding */
-                                basalelement->ElementSizes(&hx,&hy,&hz);
-                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
-                                for(row=0;row<dim;row++) 
-                                        for(col=0;col<dim;col++) 
-                                                D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
-
-                                TripleMultiply(Bprime,dim,numnodes,1,
-                                                        D,dim,dim,0,
-                                                        Bprime,dim,numnodes,0,
-                                                        &Ke->values[0],1);
-                                break;
-                        default:
-                                _error_("unknown type of stabilization in LevelsetAnalysis.cpp");
-                }
-        }
-
-        /*Clean up and return*/
-        xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(D);
-        xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(v);
-        xDelete<IssmDouble>(w);
-        xDelete<IssmDouble>(c);
-        xDelete<IssmDouble>(dlsf);
-        delete gauss;
-        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-        return Ke;
+        _error_("NOT IMPLEMENTED YET");
+
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/