Changeset 17427

Timestamp:

03/12/14 16:10:59 (11 years ago)

Author:

jbondzio

Message:

CHG: adaptation of LSM to 3D

File:

• issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (modified) (12 diffs)

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

@@ -29 +29 @@
                 }
         }
-
+        
         iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
         iomodel->FetchDataToInput(elements,VxEnum);
@@ -35 +35 @@
         iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
         
+        if(iomodel->meshtype==Mesh3DEnum){
+                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
+        }
 }
 /*}}}*/
 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
         int finiteelement=P1Enum;
+        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
         ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
+        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }
 /*}}}*/
@@ -92 +97 @@
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
+
         /*Intermediaries */
-        int  dim = 2; // solve for LSF in horizontal plane only
+        int  dim, meshtype; // solve for LSF in horizontal plane only
         int i, row, col;
         IssmDouble kappa;
@@ -102 +110 @@
         IssmDouble* xyz_list = NULL;
 
+        /*Get problem dimension*/
+        basalelement->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum:   dim = 1; break;
+                case Mesh2DhorizontalEnum: dim = 2; break;
+                case Mesh3DEnum:           dim = 2; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
         /*Fetch number of nodes and dof for this finite element*/
-        int numnodes    = element->GetNumberOfNodes();
+        int numnodes    = basalelement->GetNumberOfNodes();
 
         /*Initialize Element vector and other vectors*/
-        ElementMatrix* Ke       = element->NewElementMatrix();
+        ElementMatrix* Ke       = basalelement->NewElementMatrix();
         IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
         IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
@@ -117 +134 @@
 
         /*Retrieve all inputs and parameters*/
-        element->GetVerticesCoordinates(&xyz_list);
-        element->FindParam(&dt,TimesteppingTimeStepEnum);
-        Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
-        Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
-        Input* lsf_slopex_input  = element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-        Input* lsf_slopey_input  = element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-        Input* calvingrate_input  = element->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+        Input* vx_input  = NULL;
+        Input* vy_input  = NULL;
+        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);
+
+        if(meshtype==Mesh2DhorizontalEnum){
+                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);
+                }
+        }
         
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
+        Gauss* gauss=basalelement->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
 
-                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 D_scalar=gauss->weight*Jdet;
 
                 /* Transient */
                 if(dt!=0.){
-                        element->NodalFunctions(basis,gauss);
+                        basalelement->NodalFunctions(basis,gauss);
                         TripleMultiply(basis,numnodes,1,0,
                                                 &D_scalar,1,1,0,
@@ -144 +175 @@
 
                 /* Advection */
-                GetB(B,element,xyz_list,gauss); 
-                GetBprime(Bprime,element,xyz_list,gauss); 
+                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); 
@@ -186 +217 @@
                         case 1:
                                 /* Artificial Diffusion */
-                                element->ElementSizes(&hx,&hy,&hz);
+                                basalelement->ElementSizes(&hx,&hy,&hz);
                                 h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
                                 kappa=h*vel/2.;
@@ -203 +234 @@
                         case 2:
                                 /* Streamline Upwinding */
-                                element->ElementSizes(&hx,&hy,&hz);
+                                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++) 
@@ -229 +260 @@
         xDelete<IssmDouble>(c);
         xDelete<IssmDouble>(dlsf);
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         delete gauss;
         return Ke;
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
+        
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
 
         /*Intermediaries */
-        const int dim = 2;
-        int i, ig, k;
+        int i, ig, meshtype;
         IssmDouble  Jdet,dt;
         IssmDouble  lsf;
@@ -242 +276 @@
         
         /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
+        int numnodes = basalelement->GetNumberOfNodes();
 
         /*Initialize Element vector*/
-        ElementVector* pe = element->NewElementVector();
-        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        ElementVector* pe = basalelement->NewElementVector();
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
         
         if(dt!=0.){
@@ -253 +287 @@
 
                 /*Retrieve all inputs and parameters*/
-                element->GetVerticesCoordinates(&xyz_list);
-                Input* levelset_input     = element->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
+                basalelement->GetVerticesCoordinates(&xyz_list);
+                Input* levelset_input     = basalelement->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
 
                 /* Start  looping on the number of gaussian points: */
-                Gauss* gauss=element->NewGauss(2);
+                Gauss* gauss=basalelement->NewGauss(2);
                 for(ig=gauss->begin();ig<gauss->end();ig++){
                         gauss->GaussPoint(ig);
 
-                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                        element->NodalFunctions(basis,gauss);
+                        basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                        basalelement->NodalFunctions(basis,gauss);
 
                         /* old function value */
@@ -277 +311 @@
                 for(i=0;i<numnodes;i++) 
                         pe->values[i]=0.; 
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         return pe;
 }/*}}}*/