source: issm/oecreview/Archive/16554-17801/ISSM-17443-17444.diff

Last change on this file was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 5.0 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17443)
4+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 17444)
5@@ -123,10 +123,10 @@
6
7 iomodel->FetchData(1,MeshVertexonbedEnum);
8
9- //create penalties for nodes: no node can have a temperature over the melting point
10+ //create penalties for nodes: no node can have water above the max
11 CreateSingleNodeToElementConnectivity(iomodel);
12 for(int i=0;i<iomodel->numberofvertices;i++){
13- if (iomodel->meshtype==Mesh3DEnum){
14+ if (iomodel->meshtype!=Mesh3DEnum){
15 /*keep only this partition's nodes:*/
16 if(iomodel->my_vertices[i]){
17 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,HydrologyDCInefficientAnalysisEnum));
18@@ -279,7 +279,7 @@
19 /*Intermediaries */
20 bool active_element,isefficientlayer;
21 IssmDouble dt,scalar;
22- IssmDouble moulin_load,water_head;
23+ IssmDouble water_head;
24 IssmDouble water_load,transfer;
25 IssmDouble Jdet;
26
27@@ -306,7 +306,6 @@
28 Input* thickness_input = basalelement->GetInput(ThicknessEnum);
29 Input* bed_input = basalelement->GetInput(BedEnum);
30 Input* water_input = basalelement->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input);
31- Input* moulin_input = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
32 if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);}
33
34 IssmDouble sediment_storing = SedimentStoring(basalelement);
35@@ -326,10 +325,8 @@
36
37 /*Loading term*/
38 water_input->GetInputValue(&water_load,gauss);
39- moulin_input->GetInputValue(&moulin_load,gauss);
40-
41+
42 scalar = Jdet*gauss->weight*(water_load);
43- scalar = scalar + Jdet* moulin_load;
44 if(dt!=0.) scalar = scalar*dt;
45 for(int i=0;i<numnodes;i++){
46 pe->values[i]+=scalar*basis[i];
47Index: ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
48===================================================================
49--- ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17443)
50+++ ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17444)
51@@ -75,7 +75,11 @@
52 ElementVector* pe = analysis->CreatePVector(element);
53 element->ReduceMatrices(Ke,pe);
54 if(Ke) Ke->AddToGlobal(Kff,Kfs);
55- if(pe) pe->AddToGlobal(pf);
56+ if(pe){
57+ pe->AddToGlobal(pf);
58+ /* printf("-------------------USUAL \n"); */
59+ /* pf->Echo(); */
60+ }
61 delete Ke;
62 delete pe;
63 }
64@@ -86,6 +90,8 @@
65 if(load->InAnalysis(configuration_type)){
66 load->CreateKMatrix(Kff,Kfs);
67 load->CreatePVector(pf);
68+ /* printf("-------------------LOADING \n"); */
69+ /* pf->Echo(); */
70 }
71 }
72
73@@ -96,6 +102,8 @@
74 if(load->InAnalysis(configuration_type)){
75 load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
76 load->PenaltyCreatePVector(pf,kmax);
77+ /* printf("-------------------PENALTY \n"); */
78+ /* pf->Echo(); */
79 }
80 }
81 }
82Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
83===================================================================
84--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17443)
85+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17444)
86@@ -180,9 +180,24 @@
87 /*FUNCTION Pengrid::CreatePVector {{{*/
88 void Pengrid::CreatePVector(Vector<IssmDouble>* pf){
89
90- /*No loads applied, do nothing: */
91- return;
92+ ElementVector* pe=NULL;
93+ int analysis_type;
94+ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
95
96+ switch(analysis_type){
97+
98+ case HydrologyDCInefficientAnalysisEnum:
99+ pe = CreatePVectorHydrologyDCInefficient();
100+ break;
101+ default:
102+ /*No loads applied, do nothing: */
103+ return;
104+ }
105+ if(pe){
106+ pe->AddToGlobal(pf);
107+ delete pe;
108+ }
109+
110 }
111 /*}}}*/
112 /*FUNCTION Pengrid::GetNodesSidList{{{*/
113@@ -792,6 +807,23 @@
114 return pe;
115 }
116 /*}}}*/
117+/*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
118+ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){
119+
120+ IssmDouble moulin_load,dt;
121+
122+ /*Initialize Element matrix*/
123+ ElementVector* pe=new ElementVector(&node,1,this->parameters);
124+
125+ this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
126+ parameters->FindParam(&dt,TimesteppingTimeStepEnum);
127+
128+ if(dt!=0.0) pe->values[0]=moulin_load*dt;
129+
130+ /*Clean up and return*/
131+ return pe;
132+ }
133+/*}}}*/
134 /*FUNCTION Pengrid::ResetConstraint {{{*/
135 void Pengrid::ResetConstraint(void){
136 active = 0;
137Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.h
138===================================================================
139--- ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17443)
140+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17444)
141@@ -93,6 +93,7 @@
142 ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
143 void ConstraintActivateHydrologyDCInefficient(int* punstable);
144 void ConstraintActivate(int* punstable);
145+ ElementVector* CreatePVectorHydrologyDCInefficient(void);
146 void ResetConstraint(void);
147 /*}}}*/
148
Note: See TracBrowser for help on using the repository browser.