Changeset 14650

Timestamp:

04/18/13 13:42:43 (12 years ago)

Author:

Eric.Larour

Message:

Finished implementation of the solution core and Tria routines.
We now have a module GiaDeflectionCorex which will implement the physics
of PGR for a single disk, loaded using a history of thicknesses, for a
certain radial distance ri.

Location:

issm/trunk-jpl

Files:

• src/c/Makefile.am (modified) (1 diff)
• src/c/classes/FemModel.cpp (modified) (1 diff)
• src/c/classes/objects/Elements/Element.h (modified) (1 diff)
• src/c/classes/objects/Elements/Penta.cpp (modified) (1 diff)
• src/c/classes/objects/Elements/Penta.h (modified) (1 diff)
• src/c/classes/objects/Elements/Tria.cpp (modified) (4 diffs)
• src/c/classes/objects/Elements/Tria.h (modified) (1 diff)
• src/c/classes/objects/Inputs/TransientInput.cpp (modified) (1 diff)
• src/c/classes/objects/Inputs/TriaP1Input.cpp (modified) (1 diff)
• src/c/classes/objects/Inputs/TriaP1Input.h (modified) (1 diff)
• src/c/modules/GiaDeflectionCorex (added)
• src/c/modules/GiaDeflectionCorex/CMakeLists.txt (added)
• src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp (added)
• src/c/modules/GiaDeflectionCorex/GiaDeflectionCorex.h (added)
• src/c/modules/ModelProcessorx/Gia/UpdateElementsGia.cpp (modified) (1 diff)
• src/c/modules/modules.h (modified) (1 diff)
• src/c/solutions/gia_core.cpp (modified) (1 diff)
• test/Archives/Archive330.nc (modified) ( previous)
• test/NightlyRun/test330.m (modified) (1 diff)

issm/trunk-jpl/src/c/Makefile.am

@@ -500 +500 @@
                            ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \
                            ./modules/ModelProcessorx/Gia/CreateConstraintsGia.cpp\
-                           ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp
+                           ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\
+                           ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\
+                           ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h
 
 #}}}

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1552 +1552 @@
         for (i=0;i<elements->Size();i++){
                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                element->Deflection(wg,x,y);
-        }
-}
-/*}}}*/
-#endif
+                element->GiaDeflection(wg,x,y);
+        }
+}
+/*}}}*/
+#endif

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

@@ -101 +101 @@
 
                 #ifdef _HAVE_GIA_
-                virtual void   Deflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y)=0;
+                virtual void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y)=0;
                 #endif
 

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

@@ -3486 +3486 @@
 
 #ifdef _HAVE_GIA_
-/*FUNCTION Penta::Deflection {{{*/
-void Penta::Deflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
+/*FUNCTION Penta::GiaDeflection {{{*/
+void Penta::GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
         _error_("GIA deflection not implemented yet!");
 }

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

@@ -143 +143 @@
 
                 #ifdef _HAVE_GIA_
-                void   Deflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
+                void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
                 #endif
 

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

@@ -17 +17 @@
 #include "../../../Container/Container.h"
 #include "../../../include/include.h"
+#ifdef _HAVE_GIA_
+#include "../../../modules/GiaDeflectionCorex/GiaDeflectionCorex.h"
+#endif
 /*}}}*/
 
@@ -3040 +3043 @@
 
 #ifdef _HAVE_GIA_
-/*FUNCTION Tria::Deflection {{{*/
-void Tria::Deflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
+/*FUNCTION Tria::GiaDeflection {{{*/
+void Tria::GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x, IssmDouble* y){
 
         int i;
@@ -3048 +3051 @@
         IssmDouble x0,y0;
         IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble he;
+
+        /*thickness averages: */
+        IssmDouble* hes=NULL;
+        IssmDouble* times=NULL;
+        IssmDouble  currenttime;
+        int         numtimes;
+
+        /*output: */
+        IssmDouble  wi;
 
         /*how many dofs are we working with here? */
-        gsize=this->nodes[0]->indexing.gsize;
-
-        /*pull thickness average and area: */
+        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
+        
+        /*what time is it? :*/
+        this->parameters->FindParam(&currenttime,TimeEnum);
+
+        /*pull thickness averages: */
         Input* thickness_input=inputs->GetInput(ThicknessEnum); 
         if (!thickness_input)_error_("thickness input needed to compute gia deflection!");
-        thickness_input->GetInputAverage(&he);
+        thickness_input->GetInputAllTimeAverages(&hes,&times,&numtimes);
         
+        /*pull area of this Tria: */
         area=this->GetArea();
 
@@ -3069 +3084 @@
 
         for(i=0;i<gsize;i++){
+                /*compute distance from the center of the tria to the vertex i: */
                 xi=x[i]; yi=y[i];
                 ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2));
-        }
+
+                /*for this Tria, compute contribution to rebound at vertex i: */
+                GiaDeflectionCorex(&wi,ri,re,hes,times,currenttime,numtimes);
+
+                /*plug value into solution vector: */
+                wg->SetValue(i,wi,ADD_VAL);
+        }
+
+        /*Free ressources: */
+        xDelete<IssmDouble>(hes);
+        xDelete<IssmDouble>(times);
 
         return;

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

@@ -141 +141 @@
 
                 #ifdef _HAVE_GIA_
-                void   Deflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
+                void   GiaDeflection(Vector<IssmDouble>* wg,IssmDouble* x,IssmDouble* y);
                 #endif
 

issm/trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp

@@ -250 +250 @@
 void TransientInput::GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){
 
-        _error_("not implemented yet!");
-
+        int i;
+        IssmDouble* times=NULL;
+        IssmDouble* values=NULL;
+
+        /*allocate: */
+        times=xNew<IssmDouble>(this->numtimesteps);
+        values=xNew<IssmDouble>(this->numtimesteps);
+
+        for(i=0;i<numtimesteps;i++){
+                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+                input->GetInputAverage(values+i);
+                times[i]=this->timesteps[i];
+        }
+
+        *pvalues=values;
+        *ptimes=times;
+        *pnumtimes=numtimesteps;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/objects/Inputs/TriaP1Input.cpp

@@ -181 +181 @@
 }
 /*}}}*/
+/*FUNCTION TriaP1Input::GetInputAllTimeAverages{{{*/
+void TriaP1Input::GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){
+
+        IssmDouble* outvalues=NULL;
+        IssmDouble* times=NULL;
+        IssmDouble  numtimes;
+
+        /*this is not a transient forcing, so we only have 1 value, steady state: */
+        numtimes=1;
+        outvalues=xNew<IssmDouble>(1);
+        times=xNew<IssmDouble>(1);
+        
+        outvalues[0]=1./3.*(values[0]+values[1]+values[2]);
+        times[0]=0; /*we don't have a time*/
+
+        *pvalues=outvalues;
+        *ptimes=times;
+        *pnumtimes=numtimes;
+}
+/*}}}*/
+
 
 /*Intermediary*/

issm/trunk-jpl/src/c/classes/objects/Inputs/TriaP1Input.h

@@ -55 +55 @@
                 void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
                 void GetInputAverage(IssmDouble* pvalue);
-                void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
+                void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes);
                 void GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss);
                 void GetVyStrainRate2d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussTria* gauss);

issm/trunk-jpl/src/c/modules/ModelProcessorx/Gia/UpdateElementsGia.cpp

@@ -34 +34 @@
         }
 
-        iomodel->FetchDataToInput(elements,SurfaceEnum);
-        iomodel->FetchDataToInput(elements,BedEnum);
+        iomodel->FetchDataToInput(elements,ThicknessEnum);
 
         /*Free data: */

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

@@ -30 +30 @@
 #include "./GetVectorFromInputsx/GetVectorFromInputsx.h"
 #include "./GetVectorFromControlInputsx/GetVectorFromControlInputsx.h"
+#include "./GiaDeflectionCorex/GiaDeflectionCorex.h"
 #include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h"
 #include "./Gradjx/Gradjx.h"

issm/trunk-jpl/src/c/solutions/gia_core.cpp

@@ -41 +41 @@
         /*call the main module: */
         femmodel->Deflection(wg,x,y);
-        
+
+        /*assemble vector: */
+        wg->Assemble();
+
         InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,wg); 
 

issm/trunk-jpl/test/NightlyRun/test330.m

@@ -5 +5 @@
 md=setflowequation(md,'macayeal','all');
 md.cluster=generic('name',oshostname(),'np',3);
+%md.verbose=verbose('convergence',true,'solution',true,'module',true,'solver',true,'mprocessor',true);
 md=solve(md,GiaSolutionEnum());