Changeset 16809

Timestamp:

11/16/13 20:32:45 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: done with Mass conservation CG and DG

Location:

issm/trunk-jpl/src/c

Files:

• analyses/MasstransportAnalysis.cpp (modified) (1 diff)
• analyses/MasstransportAnalysis.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (2 diffs)

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

@@ -216 +216 @@
 }/*}}}*/
 ElementVector* MasstransportAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+        if(!element->IsOnBed()) return NULL;
+        Element* basalelement = element->SpawnBasalElement();
+
+        switch(element->FiniteElement()){
+                case P1Enum: case P2Enum:
+                        return CreatePVectorCG(basalelement);
+                case P1DGEnum:
+                        return CreatePVectorDG(basalelement);
+                default:
+                        _error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
+        }
+
+}/*}}}*/
+ElementVector* MasstransportAnalysis::CreatePVectorCG(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble  Jdet,dt;
+        IssmDouble  ms,mb,mb_correction=0.,thickness;
+        IssmDouble* xyz_list;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        Input* mb_correction_input = element->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+        Input* mb_input            = element->GetInput(BasalforcingsMeltingRateEnum);     _assert_(mb_input);
+        Input* ms_input            = element->GetInput(SurfaceforcingsMassBalanceEnum);   _assert_(ms_input);
+        Input* thickness_input     = element->GetInput(ThicknessEnum);                    _assert_(thickness_input);
+
+        /*Initialize mb_correction to 0, do not forget!:*/
+        /* 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);
+
+                ms_input->GetInputValue(&ms,gauss);
+                mb_input->GetInputValue(&mb,gauss);
+                thickness_input->GetInputValue(&thickness,gauss);
+                if(mb_correction_input)
+                 mb_correction_input->GetInputValue(&mb_correction,gauss);
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return pe;
+}/*}}}*/
+ElementVector* MasstransportAnalysis::CreatePVectorDG(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble  Jdet,dt;
+        IssmDouble  ms,mb,mb_correction=0.,thickness;
+        IssmDouble* xyz_list;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        Input* mb_correction_input = element->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+        Input* mb_input            = element->GetInput(BasalforcingsMeltingRateEnum);     _assert_(mb_input);
+        Input* ms_input            = element->GetInput(SurfaceforcingsMassBalanceEnum);   _assert_(ms_input);
+        Input* thickness_input     = element->GetInput(ThicknessEnum);                    _assert_(thickness_input);
+
+        /*Initialize mb_correction to 0, do not forget!:*/
+        /* 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);
+
+                ms_input->GetInputValue(&ms,gauss);
+                mb_input->GetInputValue(&mb,gauss);
+                thickness_input->GetInputValue(&thickness,gauss);
+                if(mb_correction_input)
+                 mb_correction_input->GetInputValue(&mb_correction,gauss);
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return pe;
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -23 +23 @@
                 ElementMatrix* CreateKMatrix(Element* element);
                 ElementVector* CreatePVector(Element* element);
+                ElementVector* CreatePVectorCG(Element* element);
+                ElementVector* CreatePVectorDG(Element* element);
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

@@ -1975 +1975 @@
 bool Tria::IsOnBed(){
 
-        bool onbed;
         int meshtype;
         this->parameters->FindParam(&meshtype,MeshTypeEnum);
@@ -1981 +1980 @@
                 case Mesh2DverticalEnum:
                         return HasEdgeOnBed();
-                        break;
                 case Mesh2DhorizontalEnum:
-                        inputs->GetInputValue(&onbed,MeshElementonbedEnum);
-                        break;
+                        return true;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
-        return onbed;
 }
 /*}}}*/