Changeset 17511

Timestamp:

03/20/14 17:16:09 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: done with depth averaging solution :)

Location:

issm/trunk-jpl/src/c

Files:

• analyses/DepthAverageAnalysis.cpp (modified) (5 diffs)
• classes/Elements/Tetra.cpp (modified) (1 diff)
• classes/Elements/Tetra.h (modified) (1 diff)
• cores/AnalysisConfiguration.cpp (modified) (2 diffs)
• cores/cores.h (modified) (1 diff)
• cores/depthaverage_core.cpp (modified) (1 diff)
• cores/surfaceslope_core.cpp (modified) (1 diff)

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

@@ -50 +50 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,D;
+        IssmDouble  Jdet,D,dt=1.e+9;
         IssmDouble *xyz_list = NULL;
 
@@ -59 +59 @@
                 case Mesh2DverticalEnum: dim = 2; break;
                 case Mesh3DEnum:         dim = 3; break;
+                case Mesh3DtetrasEnum:    dim = 3; break;
                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
         }
@@ -68 +69 @@
         ElementMatrix* Ke     = element->NewElementMatrix();
         IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -79 +79 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                element->NodalFunctions(Bprime,gauss);
                 GetB(B,element,dim,xyz_list,gauss);
-                D=gauss->weight*Jdet;
+                D=gauss->weight*Jdet*dt;
 
+                /*vertical diffusion*/
                 TripleMultiply(B,1,numnodes,1,
                                         &D,1,1,0,
-                                        Bprime,1,numnodes,0,
+                                        B,1,numnodes,0,
+                                        &Ke->values[0],1);
+
+                /*Next value*/
+                D=gauss->weight*Jdet;
+                element->NodalFunctions(B,gauss);
+                TripleMultiply(B,numnodes,1,0,
+                                        &D,1,1,0,
+                                        B,1,numnodes,0,
                                         &Ke->values[0],1);
         } 
@@ -92 +100 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
         delete gauss;
         return Ke;
 }/*}}}*/
 ElementVector* DepthAverageAnalysis::CreatePVector(Element* element){/*{{{*/
-        return NULL;
+
+        /*Intermediaries*/
+        int         input_enum;
+        IssmDouble  Jdet,scalar,value;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector*/
+        ElementVector* pe     = element->NewElementVector();
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&input_enum,InputToDepthaverageEnum);
+        Input* input = element->GetInput(input_enum); _assert_(input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                /* Build transient now */
+                input->GetInputValue(&value, gauss);
+                scalar=value*Jdet*gauss->weight;
+                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return pe;
 }/*}}}*/
 void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

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

@@ -445 +445 @@
 }
 /*}}}*/
+/*FUNCTION Tetra::InputUpdateFromSolutionOneDof{{{*/
+void  Tetra::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
+
+        /*Intermediary*/
+        int* doflist = NULL;
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Fetch dof list and allocate solution vector*/
+        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+        IssmDouble* values    = xNew<IssmDouble>(numnodes);
+
+        /*Use the dof list to index into the solution vector: */
+        for(int i=0;i<numnodes;i++){
+                values[i]=solution[doflist[i]];
+                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Add input to the element: */
+        this->inputs->AddInput(new TetraInput(enum_type,values,P1Enum));
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(values);
+        xDelete<int>(doflist);
+}
+/*}}}*/
 /*FUNCTION Tetra::IsIceInElement    THIS ONE{{{*/
 bool   Tetra::IsIceInElement(){

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

@@ -47 +47 @@
                 /*Update virtual functions resolution: {{{*/
                 void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
-                void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
+                void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum);
                 void  InputUpdateFromVector(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
 #ifdef _HAVE_DAKOTA_

issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp

@@ -26 +26 @@
 
                 case StressbalanceSolutionEnum:
-                        numanalyses=5;
+                        numanalyses=6;
                         analyses=xNew<int>(numanalyses);
                         analyses[0]=StressbalanceAnalysisEnum;
@@ -33 +33 @@
                         analyses[3]=L2ProjectionBaseAnalysisEnum;
                         analyses[4]=ExtrudeFromBaseAnalysisEnum;
+                        analyses[5]=DepthAverageAnalysisEnum;
                         break;
 

issm/trunk-jpl/src/c/cores/cores.h

@@ -30 +30 @@
 void controltao_core(FemModel* femmodel);
 void masstransport_core(FemModel* femmodel);
+void depthaverage_core(FemModel* femmodel);
 void extrudefrombase_core(FemModel* femmodel);
 void extrudefromtop_core(FemModel* femmodel);

issm/trunk-jpl/src/c/cores/depthaverage_core.cpp

@@ -12 +12 @@
 void depthaverage_core(FemModel* femmodel){
 
-        if(VerboseSolution()) _printf0_("extruding solution from base...\n");
+        if(VerboseSolution()) _printf0_("depth averaging...\n");
 
         /*Call on core computations: */

issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp

@@ -33 +33 @@
         }
         if(meshtype==Mesh2DverticalEnum){
-              femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
+                femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
                 extrudefrombase_core(femmodel);
         }