Changeset 16832

Timestamp:

11/19/13 14:25:55 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added mass transport

Location:

issm/trunk-jpl/src/c

Files:

• analyses/BalancethicknessAnalysis.cpp (modified) (2 diffs)
• analyses/MasstransportAnalysis.cpp (modified) (4 diffs)
• analyses/MasstransportAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)

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

@@ -136 +136 @@
         /*Intermediaries */
         IssmDouble  dhdt,mb,ms,Jdet;
-        IssmDouble* xyz_list;
+        IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -177 +177 @@
         /*Intermediaries */
         IssmDouble  dhdt,mb,ms,Jdet;
-        IssmDouble* xyz_list;
+        IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/

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

@@ -213 +213 @@
 /*Finite Element Analysis*/
 ElementMatrix* MasstransportAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
+
+        int meshtype;
+        element->FindParam(&meshtype,MeshTypeEnum);
+
+        switch(element->FiniteElement()){
+                case P1Enum: case P2Enum:
+                        return CreateKMatrixCG(basalelement);
+                case P1DGEnum:
+                        return CreateKMatrixDG(basalelement);
+                default:
+                        _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
+        }
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateKMatrixCG(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        int        stabilization;
+        int        meshtype;
+        IssmDouble Jdet,D_scalar,dt,h;
+        IssmDouble vel,vx,vy,dvxdx,dvydy;
+        IssmDouble dvx[2],dvy[2];
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        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];
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        element->FindParam(&meshtype,MeshTypeEnum);
+        element->FindParam(&stabilization,MasstransportStabilizationEnum);
+        Input* vxaverage_input=NULL;
+        Input* vyaverage_input=NULL;
+        if(meshtype==Mesh2DhorizontalEnum){
+                vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+                vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
+        }
+        else{
+                vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+                vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+        }
+        h = element->CharacteristicLength();
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                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);
+
+                D_scalar=gauss->weight*Jdet;
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D_scalar,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+
+                GetB(B,element,xyz_list,gauss);
+                GetBprime(Bprime,element,xyz_list,gauss);
+
+                dvxdx=dvx[0];
+                dvydy=dvy[1];
+                D_scalar=dt*gauss->weight*Jdet;
+
+                D[0][0]=D_scalar*dvxdx;
+                D[0][1]=0.;
+                D[1][0]=0.;
+                D[1][1]=D_scalar*dvydy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        B,2,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,
+                                        &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;
+                }
+                else if(stabilization==1){
+                        /*SSA*/
+                        vxaverage_input->GetInputAverage(&vx);
+                        vyaverage_input->GetInputAverage(&vy);
+                        D[0][0]=h/2.0*fabs(vx);
+                        D[0][1]=0.;
+                        D[1][0]=0.;
+                        D[1][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,
+                                                &Ke->values[0],1);
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        return Ke;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateKMatrixDG(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        int        meshtype;
+        IssmDouble Jdet,D_scalar,dt,vx,vy;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        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];
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        element->FindParam(&meshtype,MeshTypeEnum);
+        Input* vxaverage_input=NULL;
+        Input* vyaverage_input=NULL;
+        if(meshtype==Mesh2DhorizontalEnum){
+                vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+                vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
+        }
+        else{
+                vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+                vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+        }
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                vxaverage_input->GetInputValue(&vx,gauss);
+                vyaverage_input->GetInputValue(&vy,gauss);
+
+                D_scalar=gauss->weight*Jdet;
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D_scalar,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+
+                /*WARNING: B and Bprime are inverted compared to CG*/
+                GetB(Bprime,element,xyz_list,gauss);
+                GetBprime(B,element,xyz_list,gauss);
+
+                D_scalar = - dt*gauss->weight*Jdet;
+                D[0][0]  = D_scalar*vx;
+                D[0][1]  = 0.;
+                D[1][0]  = 0.;
+                D[1][1]  = D_scalar*vy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        Bprime,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        return Ke;
 }/*}}}*/
 ElementVector* MasstransportAnalysis::CreatePVector(Element* element){/*{{{*/
@@ -235 +438 @@
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,mb_correction=0.,thickness;
-        IssmDouble* xyz_list;
+        IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -281 +484 @@
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,mb_correction=0.,thickness;
-        IssmDouble* xyz_list;
+        IssmDouble* xyz_list = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -321 +524 @@
         delete gauss;
         return pe;
+}/*}}}*/
+void MasstransportAnalysis::GetB(IssmDouble* B,Element* element,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
+         * by: 
+         *       Bi=[ N ]
+         *          [ N ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        element->NodalFunctions(basis,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[numnodes*0+i] = basis[i];
+                B[numnodes*1+i] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
+}/*}}}*/
+void MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,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
+         * by: 
+         *       Bi_prime=[ dN/dx ]
+         *                [ dN/dy ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*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];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -22 +22 @@
                 /*Finite element Analysis*/
                 ElementMatrix* CreateKMatrix(Element* element);
+                ElementMatrix* CreateKMatrixCG(Element* element);
+                ElementMatrix* CreateKMatrixDG(Element* element);
                 ElementVector* CreatePVector(Element* element);
                 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 GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -37 +37 @@
                 virtual        ~Element(){};
 
+                virtual IssmDouble CharacteristicLength(void)=0;
                 virtual void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
                 virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3011 +3011 @@
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
         this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
+        if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
         this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
         this->material->inputs->DeleteInput(DamageDbarEnum);
+        this->inputs->DeleteInput(VxAverageEnum);
+        this->inputs->DeleteInput(VyAverageEnum);
 
         return tria;

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -69 +69 @@
                 /*Element virtual functions definitions: {{{*/
                 void   BasalFrictionCreateInput(void);
+                IssmDouble CharacteristicLength(void){_error_("not implemented yet");};
                 void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
                 void   ComputeStrainRate(Vector<IssmDouble>* eps);

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -70 +70 @@
                 void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
                 void        AddMaterialInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
+                IssmDouble  CharacteristicLength(void){_error_("not implemented yet");};
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
                 void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -196 +196 @@
         _assert_(this->material);
         this->material->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+/*FUNCTION Tria::AddInput{{{*/
+IssmDouble Tria::CharacteristicLength(void){
+
+        return sqrt(2*this->GetArea());
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -67 +67 @@
                 /*}}}*/
                 /*Element virtual functions definitions: {{{*/
+                IssmDouble  CharacteristicLength(void);
                 void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
                 void        ComputeStrainRate(Vector<IssmDouble>* eps);