Changeset 16899

Timestamp:

11/22/13 16:15:01 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: preparing KMatrix Hydrology DC

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyDCInefficientAnalysis.cpp (modified) (3 diffs)
• analyses/HydrologyDCInefficientAnalysis.h (modified) (1 diff)
• analyses/L2ProjectionEPLAnalysis.cpp (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)
• classes/Materials/Matpar.cpp (modified) (1 diff)

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

@@ -138 +138 @@
 /*Finite Element Analysis*/
 ElementMatrix* HydrologyDCInefficientAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
-ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
 
         /*Intermediaries*/
@@ -160 +157 @@
 
         /*Intermediaries */
+        IssmDouble  D_scalar,Jdet,dt;
+        IssmDouble *xyz_list  = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = basalelement->GetNumberOfNodes();
+
+        /*Initialize Element vector*/
+        ElementMatrix* Ke     = basalelement->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble     D[2][2]={0.};
+
+        /*Retrieve all inputs and parameters*/
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        IssmDouble sediment_storing       = SedimentStoring(basalelement);
+        IssmDouble sediment_transmitivity = basalelement->GetMaterialParameter(HydrologydcSedimentTransmitivityEnum);
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                /*Diffusivity*/
+                D_scalar=sediment_transmitivity*gauss->weight*Jdet;
+                if(dt!=0.) D_scalar=D_scalar*dt;
+                D[0][0]=D_scalar;
+                D[1][1]=D_scalar;
+                GetB(B,element,xyz_list,gauss); 
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        B,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+                /*Transient*/
+                if(dt!=0.){
+                        basalelement->NodalFunctions(&basis[0],gauss);
+                        D_scalar=sediment_storing*gauss->weight*Jdet;
+                        TripleMultiply(basis,numnodes,1,0,
+                                                &D_scalar,1,1,0,
+                                                basis,1,numnodes,0,
+                                                &Ke->values[0],1);
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        return Ke;
+}/*}}}*/
+ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        int      meshtype;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Intermediaries */
         IssmDouble dt,scalar,water_head;
         IssmDouble water_load,transfer;
@@ -210 +282 @@
         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         return pe;
+}/*}}}*/
+void HydrologyDCInefficientAnalysis::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=[ 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++){
+                B[numnodes*0+i] = dbasis[0*numnodes+i];
+                B[numnodes*1+i] = dbasis[1*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -27 +27 @@
 
                 /*Intermediaries*/
+                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 IssmDouble SedimentStoring(Element* element);
                 void ElementizeEplMask(FemModel* femmodel);

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

@@ -71 +71 @@
 /*Finite Element Analysis*/
 ElementMatrix* L2ProjectionEPLAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        /*Intermediaries*/
+        int      meshtype;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh2DverticalEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Intermediaries */
+        IssmDouble  D,Jdet;
+        IssmDouble *xyz_list  = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = basalelement->GetNumberOfNodes();
+
+        /*Initialize Element vector*/
+        ElementMatrix* Ke    = basalelement->NewElementMatrix(NoneApproximationEnum);
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        basalelement->GetVerticesCoordinates(&xyz_list);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctions(basis,gauss);
+                D=gauss->weight*Jdet;
+
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        return Ke;
 }/*}}}*/
 ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/

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

@@ -4790 +4790 @@
 
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(2,3);
+        gauss=new GaussPenta(3,3);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -260 +260 @@
                 case HydrologydcEplCompressibilityEnum:     return this->epl_compressibility;
                 case HydrologydcWaterCompressibilityEnum:   return this->water_compressibility;
+                case HydrologydcSedimentTransmitivityEnum:  return this->sediment_transmitivity;
                 case ConstantsGEnum:                        return this->g;
                 default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");