Changeset 17426

Timestamp:

03/12/14 15:55:03 (11 years ago)

Author:

jbondzio

Message:

CHG: adaptation to LSM in 3D

File:

• issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (modified) (10 diffs)

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

@@ -33 +33 @@
 void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
         int finiteelement=P1Enum;
+        if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
         ::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement);
+        iomodel->DeleteData(2,MeshVertexonbedEnum,MeshVertexonsurfaceEnum);
 }
 /*}}}*/
@@ -76 +78 @@
 ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
+
         /*Intermediaries */
-        const int dim = 2;
+        int               meshtype,dim;
         int        i,row,col,stabilization;
         bool       extrapolatebydiffusion = true;
         IssmDouble Jdet,D_scalar,h;
-        IssmDouble dlsf[dim],normal[dim];
         IssmDouble norm_dlsf;
         IssmDouble hx,hy,hz,kappa;
         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*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble     D[dim][dim];
+        IssmDouble*             dlsf   = xNew<IssmDouble>(dim);
+        IssmDouble*             normal = xNew<IssmDouble>(dim);
+        IssmDouble*    D                 = xNew<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
-        Input* lsf_slopex_input=element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-        Input* lsf_slopey_input=element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-        element->GetVerticesCoordinates(&xyz_list);
-        h = element->CharacteristicLength();
+        Input* lsf_slopex_input=basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+        Input* lsf_slopey_input=basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        h = basalelement->CharacteristicLength();
 
         /* 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);
-                GetB(B,element,xyz_list,gauss);
-                GetBprime(Bprime,element,xyz_list,gauss);
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,basalelement,xyz_list,gauss);
+                GetBprime(Bprime,basalelement,xyz_list,gauss);
                 
                 D_scalar=gauss->weight*Jdet;
@@ -117 +132 @@
                                 for(col=0;col<dim;col++)
                                         if(row==col)
-                                                D[row][col] = D_scalar;
+                                                D[row*dim+col] = D_scalar;
                                         else
-                                                D[row][col] = 0.;
+                                                D[row*dim+col] = 0.;
 
                         TripleMultiply(Bprime,dim,numnodes,1,
-                                        &D[0][0],dim,dim,0,
+                                        D,dim,dim,0,
                                         Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
@@ -143 +158 @@
                                 for(col=0;col<dim;col++)
                                         if(row==col)
-                                                D[row][col]=D_scalar*normal[row];
+                                                D[row*dim+col] = D_scalar*normal[row];
                                         else
-                                                D[row][col]=0.;
+                                                D[row*dim+col] = 0.;
                         TripleMultiply(B,dim,numnodes,1,
-                                                &D[0][0],dim,dim,0,
+                                                D,dim,dim,0,
                                                 Bprime,dim,numnodes,0,
                                                 &Ke->values[0],1);
@@ -156 +171 @@
                         else if(stabilization==1){
                                 /* Artificial Diffusion */
-                                element->ElementSizes(&hx,&hy,&hz);
+                                basalelement->ElementSizes(&hx,&hy,&hz);
                                 h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
                                 kappa=h/2.+1.e-14; 
-                                D[0][0]=D_scalar*kappa;
-                                D[0][1]=0.;
-                                D[1][0]=0.;
-                                D[1][1]=D_scalar*kappa;
+                                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[0][0],dim,dim,0,
+                                                        D,dim,dim,0,
                                                         Bprime,dim,numnodes,0,
                                                         &Ke->values[0],1);
@@ -172 +189 @@
                                 for(row=0;row<dim;row++)
                                         for(col=0;col<dim;col++)
-                                                D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+                                                D[row*dim+col]=h/(2.*1.)*normal[row]*normal[col];
 
                                 TripleMultiply(Bprime,dim,numnodes,1,
-                                                        &D[0][0],dim,dim,0,
+                                                        D,dim,dim,0,
                                                         Bprime,dim,numnodes,0,
                                                         &Ke->values[0],1);
@@ -181 +198 @@
                 }
         }/*}}}*/
+
+        //TESTING
+        _printf_("in element: " << basalelement->Id() << "\n");
+        Ke->Echo();
 
         /*Clean up and return*/
@@ -186 +207 @@
         xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(dlsf);
+        xDelete<IssmDouble>(normal);
         delete gauss;
         return Ke;
@@ -191 +214 @@
 }/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
+
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
 
         /*Intermediaries */
@@ -196 +222 @@
         
         /*Fetch number of nodes */
-        int numnodes = element->GetNumberOfNodes();
+        int numnodes = basalelement->GetNumberOfNodes();
 
         /*Initialize Element vector*/
-        ElementVector* pe = element->NewElementVector();
+        ElementVector* pe = basalelement->NewElementVector();
         for(i=0;i<numnodes;i++) 
                 pe->values[i]=0.;