Changeset 17331

Timestamp:

02/21/14 09:21:47 (11 years ago)

Author:

seroussi

Message:

NEW: added mass transport for flowline models

Location:

issm/trunk-jpl/src/c/analyses

Files:

• MasstransportAnalysis.cpp (modified) (12 diffs)
• MasstransportAnalysis.h (modified) (1 diff)

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

@@ -261 +261 @@
         /*Intermediaries */
         int        stabilization;
-        int        meshtype;
+        int        meshtype,dim;
         IssmDouble Jdet,D_scalar,dt,h;
         IssmDouble vel,vx,vy,dvxdx,dvydy;
@@ -267 +267 @@
         IssmDouble* xyz_list = NULL;
 
+        /*Get problem dimension*/
+        element->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();
@@ -273 +282 @@
         ElementMatrix* Ke     = element->NewElementMatrix();
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
-        IssmDouble     D[2][2]={0.};
+        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
@@ -290 +299 @@
         else{
                 vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
-                vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+                if(dim==2) vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+
         }
         h = element->CharacteristicLength();
@@ -303 +313 @@
 
                 vxaverage_input->GetInputValue(&vx,gauss);
-                vyaverage_input->GetInputValue(&vy,gauss);
                 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
-                vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                if(dim==2){
+                        vyaverage_input->GetInputValue(&vy,gauss);
+                        vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                }
 
                 D_scalar=gauss->weight*Jdet;
@@ -313 +325 @@
                                         &Ke->values[0],1);
 
-                GetB(B,element,xyz_list,gauss);
-                GetBprime(Bprime,element,xyz_list,gauss);
+                GetB(B,element,dim,xyz_list,gauss);
+                GetBprime(Bprime,element,dim,xyz_list,gauss);
 
                 dvxdx=dvx[0];
-                dvydy=dvy[1];
+                if(dim==2) dvydy=dvy[1];
                 D_scalar=dt*gauss->weight*Jdet;
 
-                D[0][0]=D_scalar*dvxdx;
-                D[1][1]=D_scalar*dvydy;
-                TripleMultiply(B,2,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        B,2,numnodes,0,
+                D[0*dim+0]=D_scalar*dvxdx;
+                if(dim==2) D[1*dim+1]=D_scalar*dvydy;
+
+                TripleMultiply(B,dim,numnodes,1,
+                                        D,dim,dim,0,
+                                        B,dim,numnodes,0,
                                         &Ke->values[0],1);
 
-                D[0][0]=D_scalar*vx;
-                D[1][1]=D_scalar*vy;
-                TripleMultiply(B,2,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        Bprime,2,numnodes,0,
+                D[0*dim+0]=D_scalar*vx;
+                if(dim==2) D[1*dim+1]=D_scalar*vy;
+
+                TripleMultiply(B,dim,numnodes,1,
+                                        D,dim,dim,0,
+                                        Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
 
                 if(stabilization==2){
-                        /*Streamline upwinding*/
-                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
-                        D[0][0]=h/(2*vel)*vx*vx;
-                        D[1][0]=h/(2*vel)*vy*vx;
-                        D[0][1]=h/(2*vel)*vx*vy;
-                        D[1][1]=h/(2*vel)*vy*vy;
+                        if(dim==1){
+                                vel=fabs(vx)+1.e-8;
+                                D[0]=h/(2*vel)*vx*vx;
+                        }
+                        else{
+                                /*Streamline upwinding*/
+                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                                D[0*dim+0]=h/(2*vel)*vx*vx;
+                                D[1*dim+0]=h/(2*vel)*vy*vx;
+                                D[0*dim+1]=h/(2*vel)*vx*vy;
+                                D[1*dim+1]=h/(2*vel)*vy*vy;
+                        }
                 }
                 else if(stabilization==1){
                         /*SSA*/
                         vxaverage_input->GetInputAverage(&vx);
-                        vyaverage_input->GetInputAverage(&vy);
-                        D[0][0]=h/2.0*fabs(vx);
-                        D[1][1]=h/2.0*fabs(vy);
+                        if(dim==2) vyaverage_input->GetInputAverage(&vy);
+                        D[0*dim+0]=h/2.0*fabs(vx);
+                        if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
                 }
                 if(stabilization==1 || stabilization==2){
-                        D[0][0]=D_scalar*D[0][0];
-                        D[1][0]=D_scalar*D[1][0];
-                        D[0][1]=D_scalar*D[0][1];
-                        D[1][1]=D_scalar*D[1][1];
-                        TripleMultiply(Bprime,2,numnodes,1,
-                                                &D[0][0],2,2,0,
-                                                Bprime,2,numnodes,0,
+                        if(dim==1) D[0]=D_scalar*D[0];
+                        else{
+                                D[0*dim+0]=D_scalar*D[0*dim+0];
+                                D[1*dim+0]=D_scalar*D[1*dim+0];
+                                D[0*dim+1]=D_scalar*D[0*dim+1];
+                                D[1*dim+1]=D_scalar*D[1*dim+1];
+                        }
+
+                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                D,dim,dim,0,
+                                                Bprime,dim,numnodes,0,
                                                 &Ke->values[0],1);
                 }
@@ -366 +390 @@
         xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(D);
         delete gauss;
         return Ke;
@@ -422 +447 @@
 
                 /*WARNING: B and Bprime are inverted compared to CG*/
-                GetB(Bprime,element,xyz_list,gauss);
-                GetBprime(B,element,xyz_list,gauss);
+                GetB(Bprime,element,2,xyz_list,gauss);
+                GetBprime(B,element,2,xyz_list,gauss);
 
                 D_scalar = - dt*gauss->weight*Jdet;
@@ -568 +593 @@
         return pe;
 }/*}}}*/
-void MasstransportAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void MasstransportAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
          * For node i, Bi can be expressed in the actual coordinate system
@@ -588 +613 @@
         /*Build B: */
         for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = basis[i];
-                B[numnodes*1+i] = basis[i];
+                for(int j=0;j<dim;j++){
+                        B[numnodes*j+i] = basis[i];
+                }
         }
 
@@ -595 +621 @@
         xDelete<IssmDouble>(basis);
 }/*}}}*/
-void MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
          * For node i, Bi' can be expressed in the actual coordinate system
@@ -610 +636 @@
 
         /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B': */
         for(int i=0;i<numnodes;i++){
-                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+                for(int j=0;j<dim;j++){
+                        Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
+                }
         }
 

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h

@@ -30 +30 @@
                 ElementVector* CreatePVectorCG(Element* element);
                 ElementVector* CreatePVectorDG(Element* element);
-                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-                void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+                void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+                void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);