Changeset 23970

Timestamp:

05/31/19 11:24:02 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: adding Neumann flux to GlaDS

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (2 diffs)
• classes/Loads/Neumannflux.cpp (modified) (3 diffs)
• classes/Loads/Neumannflux.h (modified) (1 diff)

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

@@ -60 +60 @@
         }
         iomodel->DeleteData(1,"md.mesh.vertexonbase");
+
+        /*Deal with Neumann BC*/
+        int M,N;
+        int *segments = NULL;
+        iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
+
+        /*Check that the size seem right*/
+        _assert_(N==3); _assert_(M>=3);
+        for(int i=0;i<M;i++){
+                if(iomodel->my_elements[segments[i*3+2]-1]){
+                        loads->AddObject(new Neumannflux(i+1,i,iomodel,segments));
+                }
+        }
+
 }/*}}}*/
 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
@@ -108 +122 @@
         iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
         iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
+        iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum,0.);
         iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);

issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp

@@ -187 +187 @@
                         /*Nothing!*/
                         break;
+                case HydrologyGlaDSAnalysisEnum:
+                        /*Nothing!*/
+                        break;
                 default:
                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -209 +212 @@
                 case HydrologyShaktiAnalysisEnum:
                         pe=CreatePVectorHydrologyShakti();
+                        break;
+                case HydrologyGlaDSAnalysisEnum:
+                        pe=CreatePVectorHydrologyGlaDS();
                         break;
                 default:
@@ -379 +385 @@
 }
 /*}}}*/
+ElementVector* Neumannflux::CreatePVectorHydrologyGlaDS(void){/*{{{*/
+
+        /* constants*/
+        const int numdof=2;
+
+        /* Intermediaries*/
+        IssmDouble Jdet,flux;
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble basis[numdof];
+
+        /*Initialize Load Vector and return if necessary*/
+        Tria*  tria=(Tria*)element;
+        _assert_(tria->FiniteElement()==P1Enum); 
+        if(!tria->IsIceInElement() || tria->IsFloating()) return NULL;
+
+        /*Initialize Element vector and other vectors*/
+        ElementVector* pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        Input* flux_input = tria->inputs->GetInput(HydrologyNeumannfluxEnum);  _assert_(flux_input); 
+
+        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
+        int index1=tria->GetNodeIndex(nodes[0]);
+        int index2=tria->GetNodeIndex(nodes[1]);
+
+        /* Start  looping on the number of gaussian points: */
+        GaussTria* gauss=new GaussTria(index1,index2,2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
+                flux_input->GetInputValue(&flux,gauss);
+
+                for(int i=0;i<numdof;i++) pe->values[i] += gauss->weight*Jdet*flux*basis[i];
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        return pe;
+}
+/*}}}*/

issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h

@@ -72 +72 @@
                 /*Neumannflux management:{{{*/
                 ElementVector* CreatePVectorHydrologyShakti(void);
+                ElementVector* CreatePVectorHydrologyGlaDS(void);
                 /*}}}*/