source: issm/oecreview/Archive/16554-17801/ISSM-17515-17516.diff@ 17802

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

Added archives
File size: 73.3 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17515)
4+++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17516)
5@@ -213,8 +213,8 @@
6 IssmDouble *xyz_list = NULL;
7
8 /*Fetch number of nodes and dof for this finite element*/
9- int vnumnodes = element->GetNumberOfNodesVelocity();
10- int pnumnodes = element->GetNumberOfNodesPressure();
11+ int vnumnodes = element->NumberofNodesVelocity();
12+ int pnumnodes = element->NumberofNodesPressure();
13 int numdof = vnumnodes*3 + pnumnodes;
14
15 /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
16@@ -319,8 +319,8 @@
17 }
18
19 /*Fetch number of nodes and dof for this finite element*/
20- int vnumnodes = element->GetNumberOfNodesVelocity();
21- int pnumnodes = element->GetNumberOfNodesPressure();
22+ int vnumnodes = element->NumberofNodesVelocity();
23+ int pnumnodes = element->NumberofNodesPressure();
24
25 /*Prepare coordinate system list*/
26 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
27@@ -938,8 +938,8 @@
28 IssmDouble FSreconditioning;
29
30 /*Fetch number of nodes and dof for this finite element*/
31- int vnumnodes = element->GetNumberOfNodesVelocity();
32- int pnumnodes = element->GetNumberOfNodesPressure();
33+ int vnumnodes = element->NumberofNodesVelocity();
34+ int pnumnodes = element->NumberofNodesPressure();
35 int vnumdof = vnumnodes*3;
36 int pnumdof = pnumnodes*1;
37
38Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
39===================================================================
40--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17515)
41+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17516)
42@@ -194,7 +194,7 @@
43 IssmDouble *vertex_pairing=NULL;
44 IssmDouble *nodeonbed=NULL;
45 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
46- iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
47+ if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
48
49 for(int i=0;i<numvertex_pairing;i++){
50
51@@ -204,7 +204,9 @@
52 _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
53
54 /*Skip if one of the two is not on the bed*/
55- if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
56+ if(iomodel->meshtype!=Mesh2DhorizontalEnum){
57+ if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
58+ }
59
60 /*Get node ids*/
61 penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
62Index: ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp
63===================================================================
64--- ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17515)
65+++ ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17516)
66@@ -34,7 +34,7 @@
67 iomodel->FetchDataToInput(elements,SurfaceEnum);
68 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
69 iomodel->FetchDataToInput(elements,VxEnum);
70- iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
71+ if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
72 if(iomodel->meshtype==Mesh3DEnum){
73 iomodel->FetchDataToInput(elements,VzEnum);
74 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
75@@ -68,7 +68,7 @@
76 IssmDouble *vertex_pairing=NULL;
77 IssmDouble *nodeonsurface=NULL;
78 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
79- iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
80+ if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
81 for(int i=0;i<numvertex_pairing;i++){
82
83 if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
84@@ -77,7 +77,9 @@
85 _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
86
87 /*Skip if one of the two is not on the bed*/
88- if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
89+ if(iomodel->meshtype!=Mesh2DhorizontalEnum){
90+ if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
91+ }
92
93 /*Get node ids*/
94 penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
95Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
96===================================================================
97--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17515)
98+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17516)
99@@ -2753,8 +2753,8 @@
100 }
101
102 /*Fetch number of nodes and dof for this finite element*/
103- int vnumnodes = element->GetNumberOfNodesVelocity();
104- int pnumnodes = element->GetNumberOfNodesPressure();
105+ int vnumnodes = element->NumberofNodesVelocity();
106+ int pnumnodes = element->NumberofNodesPressure();
107
108 /*Initialize output vector*/
109 ElementVector* de = element->NewElementVector(FSvelocityEnum);
110@@ -2783,8 +2783,8 @@
111 IssmDouble *xyz_list = NULL;
112
113 /*Fetch number of nodes and dof for this finite element*/
114- int vnumnodes = element->GetNumberOfNodesVelocity();
115- int pnumnodes = element->GetNumberOfNodesPressure();
116+ int vnumnodes = element->NumberofNodesVelocity();
117+ int pnumnodes = element->NumberofNodesPressure();
118 int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1;
119
120 /*Prepare coordinate system list*/
121@@ -2885,8 +2885,8 @@
122 }
123
124 /*Fetch number of nodes and dof for this finite element*/
125- int vnumnodes = element->GetNumberOfNodesVelocity();
126- int pnumnodes = element->GetNumberOfNodesPressure();
127+ int vnumnodes = element->NumberofNodesVelocity();
128+ int pnumnodes = element->NumberofNodesPressure();
129 int numdof = vnumnodes*dim + pnumnodes;
130 int bsize = epssize + 2;
131
132@@ -2969,8 +2969,8 @@
133 }
134
135 /*Fetch number of nodes and dof for this finite element*/
136- int vnumnodes = element->GetNumberOfNodesVelocity();
137- int pnumnodes = element->GetNumberOfNodesPressure();
138+ int vnumnodes = element->NumberofNodesVelocity();
139+ int pnumnodes = element->NumberofNodesPressure();
140 int numdof = vnumnodes*dim + pnumnodes;
141
142 /*Initialize Element matrix and vectors*/
143@@ -3064,8 +3064,8 @@
144 }
145
146 /*Fetch number of nodes and dof for this finite element*/
147- int vnumnodes = element->GetNumberOfNodesVelocity();
148- int pnumnodes = element->GetNumberOfNodesPressure();
149+ int vnumnodes = element->NumberofNodesVelocity();
150+ int pnumnodes = element->NumberofNodesPressure();
151 int numdof = vnumnodes*dim + pnumnodes;
152
153 /*Initialize Element matrix and vectors*/
154@@ -3125,8 +3125,8 @@
155 }
156
157 /*Fetch number of nodes and dof for this finite element*/
158- int vnumnodes = element->GetNumberOfNodesVelocity();
159- int pnumnodes = element->GetNumberOfNodesPressure();
160+ int vnumnodes = element->NumberofNodesVelocity();
161+ int pnumnodes = element->NumberofNodesPressure();
162
163 /*Prepare coordinate system list*/
164 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
165@@ -3208,8 +3208,8 @@
166 }
167
168 /*Fetch number of nodes and dof for this finite element*/
169- int vnumnodes = element->GetNumberOfNodesVelocity();
170- int pnumnodes = element->GetNumberOfNodesPressure();
171+ int vnumnodes = element->NumberofNodesVelocity();
172+ int pnumnodes = element->NumberofNodesPressure();
173
174 /*Prepare coordinate system list*/
175 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
176@@ -3287,8 +3287,8 @@
177 }
178
179 /*Fetch number of nodes and dof for this finite element*/
180- int vnumnodes = element->GetNumberOfNodesVelocity();
181- int pnumnodes = element->GetNumberOfNodesPressure();
182+ int vnumnodes = element->NumberofNodesVelocity();
183+ int pnumnodes = element->NumberofNodesPressure();
184
185 /*Prepare coordinate system list*/
186 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
187@@ -3364,8 +3364,8 @@
188 }
189
190 /*Fetch number of nodes and dof for this finite element*/
191- int vnumnodes = element->GetNumberOfNodesVelocity();
192- int pnumnodes = element->GetNumberOfNodesPressure();
193+ int vnumnodes = element->NumberofNodesVelocity();
194+ int pnumnodes = element->NumberofNodesPressure();
195 int numvertices = element->GetNumberOfVertices();
196
197 /*Prepare coordinate system list*/
198@@ -3658,8 +3658,8 @@
199 */
200
201 /*Fetch number of nodes for this finite element*/
202- int pnumnodes = element->GetNumberOfNodesPressure();
203- int vnumnodes = element->GetNumberOfNodesVelocity();
204+ int pnumnodes = element->NumberofNodesPressure();
205+ int vnumnodes = element->NumberofNodesVelocity();
206 int pnumdof = pnumnodes;
207 int vnumdof = vnumnodes*dim;
208
209@@ -3785,8 +3785,8 @@
210 }
211
212 /*Fetch number of nodes and dof for this finite element*/
213- int vnumnodes = element->GetNumberOfNodesVelocity();
214- int pnumnodes = element->GetNumberOfNodesPressure();
215+ int vnumnodes = element->NumberofNodesVelocity();
216+ int pnumnodes = element->NumberofNodesPressure();
217 int vnumdof = vnumnodes*dim;
218 int pnumdof = pnumnodes*1;
219
220@@ -4320,8 +4320,8 @@
221 element->GetInputValue(&approximation,ApproximationEnum);
222 if(element->IsFloating() || !element->IsOnBed()) return NULL;
223
224- int vnumnodes = element->GetNumberOfNodesVelocity();
225- int pnumnodes = element->GetNumberOfNodesPressure();
226+ int vnumnodes = element->NumberofNodesVelocity();
227+ int pnumnodes = element->NumberofNodesPressure();
228 int numnodes = 2*vnumnodes-1+pnumnodes;
229
230 /*Prepare node list*/
231@@ -4438,8 +4438,8 @@
232 Element* pentabase=element->GetBasalElement();
233 Element* basaltria=pentabase->SpawnBasalElement();
234
235- int vnumnodes = element->GetNumberOfNodesVelocity();
236- int pnumnodes = element->GetNumberOfNodesPressure();
237+ int vnumnodes = element->NumberofNodesVelocity();
238+ int pnumnodes = element->NumberofNodesPressure();
239 int numnodes = 2*vnumnodes-1+pnumnodes;
240
241 /*Prepare node list*/
242@@ -4607,8 +4607,8 @@
243 element->GetInputValue(&approximation,ApproximationEnum);
244 if(approximation!=HOFSApproximationEnum) return NULL;
245
246- int vnumnodes = element->GetNumberOfNodesVelocity();
247- int pnumnodes = element->GetNumberOfNodesPressure();
248+ int vnumnodes = element->NumberofNodesVelocity();
249+ int pnumnodes = element->NumberofNodesPressure();
250 int numnodes = vnumnodes+pnumnodes;
251
252 /*Prepare coordinate system list*/
253@@ -4686,8 +4686,8 @@
254 /*Initialize Element vector and return if necessary*/
255 element->GetInputValue(&approximation,ApproximationEnum);
256 if(approximation!=HOFSApproximationEnum) return NULL;
257- int vnumnodes = element->GetNumberOfNodesVelocity();
258- int pnumnodes = element->GetNumberOfNodesPressure();
259+ int vnumnodes = element->NumberofNodesVelocity();
260+ int pnumnodes = element->NumberofNodesPressure();
261
262 /*Prepare coordinate system list*/
263 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
264@@ -4770,8 +4770,8 @@
265 if(!element->IsOnBed() || element->IsFloating()) return NULL;
266 element->GetInputValue(&approximation,ApproximationEnum);
267 if(approximation!=SSAFSApproximationEnum) return NULL;
268- int vnumnodes = element->GetNumberOfNodesVelocity();
269- int pnumnodes = element->GetNumberOfNodesPressure();
270+ int vnumnodes = element->NumberofNodesVelocity();
271+ int pnumnodes = element->NumberofNodesPressure();
272
273 /*Prepare coordinate system list*/
274 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
275@@ -4846,8 +4846,8 @@
276 /*Initialize Element vector and return if necessary*/
277 element->GetInputValue(&approximation,ApproximationEnum);
278 if(approximation!=SSAFSApproximationEnum) return NULL;
279- int vnumnodes = element->GetNumberOfNodesVelocity();
280- int pnumnodes = element->GetNumberOfNodesPressure();
281+ int vnumnodes = element->NumberofNodesVelocity();
282+ int pnumnodes = element->NumberofNodesPressure();
283
284 /*Prepare coordinate system list*/
285 int* cs_list = xNew<int>(vnumnodes+pnumnodes);
286Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
287===================================================================
288--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17515)
289+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17516)
290@@ -72,7 +72,7 @@
291 }/*}}}*/
292 void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
293
294- iomodel->FetchData(1,MeshVertexonbedEnum);
295+ if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbedEnum);
296 ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum);
297 iomodel->DeleteData(1,MeshVertexonbedEnum);
298 }/*}}}*/
299@@ -86,21 +86,8 @@
300 }/*}}}*/
301 void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
302
303- /*create penalties for nodes: no node can have a damage > 1*/
304- iomodel->FetchData(1,DamageSpcdamageEnum);
305- CreateSingleNodeToElementConnectivity(iomodel);
306+ /*Nothing for now*/
307
308- for(int i=0;i<iomodel->numberofvertices;i++){
309-
310- /*keep only this partition's nodes:*/
311- if(iomodel->my_vertices[i]){
312- if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes!
313- loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum));
314- }
315- }
316- }
317- iomodel->DeleteData(1,DamageSpcdamageEnum);
318-
319 }/*}}}*/
320
321 /*Finite Element Analysis*/
322Index: ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
323===================================================================
324--- ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17515)
325+++ ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17516)
326@@ -59,7 +59,7 @@
327 IssmDouble *vertex_pairing=NULL;
328 IssmDouble *nodeonbed=NULL;
329 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
330- iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
331+ if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
332 for(int i=0;i<numvertex_pairing;i++){
333
334 if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
335@@ -68,7 +68,9 @@
336 _assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
337
338 /*Skip if one of the two is not on the bed*/
339- if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
340+ if(iomodel->meshtype!=Mesh2DhorizontalEnum){
341+ if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
342+ }
343
344 /*Get node ids*/
345 penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
346Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp
347===================================================================
348--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp (revision 17515)
349+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp (revision 17516)
350@@ -1,80 +0,0 @@
351-/*
352- * \brief: solutionsequence_damage_nonlinear.cpp: core of the damage solution
353- */
354-
355-#include "../toolkits/toolkits.h"
356-#include "../classes/classes.h"
357-#include "../shared/shared.h"
358-#include "../modules/modules.h"
359-
360-void solutionsequence_damage_nonlinear(FemModel* femmodel){
361-
362- /*solution : */
363- Vector<IssmDouble>* Dg=NULL;
364- Vector<IssmDouble>* Df=NULL;
365- Vector<IssmDouble>* Df_old=NULL;
366- Vector<IssmDouble>* ys=NULL;
367-
368- /*intermediary: */
369- Matrix<IssmDouble>* Kff=NULL;
370- Matrix<IssmDouble>* Kfs=NULL;
371- Vector<IssmDouble>* pf=NULL;
372- Vector<IssmDouble>* df=NULL;
373-
374- bool converged;
375- int constraints_converged;
376- int num_unstable_constraints;
377- int count;
378- int damage_penalty_threshold;
379- int damage_maxiter;
380-
381- /*parameters:*/
382- int configuration_type;
383-
384- /*Recover parameters: */
385- femmodel->parameters->FindParam(&damage_penalty_threshold,DamagePenaltyThresholdEnum);
386- femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
387- femmodel->parameters->FindParam(&damage_maxiter,DamageMaxiterEnum);
388-
389- count=1;
390- converged=false;
391-
392- InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
393- InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
394- femmodel->UpdateConstraintsx();
395-
396- for(;;){
397-
398- delete Df_old; Df_old=Df;
399- SystemMatricesx(&Kff, &Kfs, &pf,&df, NULL,femmodel);
400- CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
401- Reduceloadx(pf, Kfs, ys); delete Kfs;
402- Solverx(&Df, Kff, pf,Df_old, df, femmodel->parameters);
403- delete Kff;delete pf;delete Dg; delete df;
404- Mergesolutionfromftogx(&Dg, Df,ys,femmodel->nodes,femmodel->parameters); delete ys;
405- InputUpdateFromSolutionx(femmodel,Dg);
406-
407- ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel);
408-
409- if (!converged){
410- if(VerboseConvergence()) _printf0_(" #unstable constraints = " << num_unstable_constraints << "\n");
411- if (num_unstable_constraints <= damage_penalty_threshold)converged=true;
412- if (count>=damage_maxiter){
413- converged=true;
414- _printf0_(" maximum number of iterations (" << damage_maxiter << ") exceeded\n");
415- }
416- }
417- count++;
418-
419- InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
420-
421- if(converged)break;
422- }
423-
424- InputUpdateFromSolutionx(femmodel,Dg);
425-
426- /*Free ressources: */
427- delete Dg;
428- delete Df;
429- delete Df_old;
430-}
431Index: ../trunk-jpl/src/c/solutionsequences/solutionsequences.h
432===================================================================
433--- ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17515)
434+++ ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17516)
435@@ -12,7 +12,6 @@
436 #include "../shared/Numerics/types.h"
437
438 void solutionsequence_thermal_nonlinear(FemModel* femmodel);
439-void solutionsequence_damage_nonlinear(FemModel* femmodel);
440 void solutionsequence_hydro_nonlinear(FemModel* femmodel);
441 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
442 void solutionsequence_newton(FemModel* femmodel);
443Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
444===================================================================
445--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17515)
446+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17516)
447@@ -242,9 +242,6 @@
448 case HydrologyDCInefficientAnalysisEnum:
449 Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax);
450 break;
451- case DamageEvolutionAnalysisEnum:
452- Ke=PenaltyCreateKMatrixDamageEvolution(kmax);
453- break;
454 default:
455 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
456 }
457@@ -276,9 +273,6 @@
458 case HydrologyDCInefficientAnalysisEnum:
459 pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax);
460 break;
461- case DamageEvolutionAnalysisEnum:
462- pe=PenaltyCreatePVectorDamageEvolution(kmax);
463- break;
464 default:
465 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
466 }
467@@ -416,11 +410,6 @@
468 ConstraintActivateHydrologyDCInefficient(punstable);
469 return;
470 }
471- else if (analysis_type==DamageEvolutionAnalysisEnum){
472- ConstraintActivateDamageEvolution(punstable);
473- return;
474- }
475-
476 else{
477 _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet");
478 }
479@@ -709,104 +698,6 @@
480 return pe;
481 }
482 /*}}}*/
483-/*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/
484-void Pengrid::ConstraintActivateDamageEvolution(int* punstable){
485-
486- // The penalty is stable if it doesn't change during to successive iterations.
487- IssmDouble max_damage;
488- IssmDouble damage;
489- int new_active;
490- int unstable=0;
491- int penalty_lock;
492-
493- /*check that pengrid is not a clone (penalty to be added only once)*/
494- if (node->IsClone()){
495- unstable=0;
496- *punstable=unstable;
497- return;
498- }
499-
500- //First recover damage using the element: */
501- element->GetInputValue(&damage,node,DamageDEnum);
502-
503- //Recover our data:
504- parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
505- parameters->FindParam(&max_damage,DamageMaxDamageEnum);
506-
507- //Figure out if damage>max_damage, in which case, this penalty needs to be activated.
508- //Would need to do the same for damage<0 if penalties are used. For now, ConstraintStatex
509- //is not called in solutionsequence_damage_nonlinear, so no penalties are applied.
510-
511- if (damage>max_damage){
512- new_active=1;
513- }
514- else{
515- new_active=0;
516- }
517-
518- //Figure out stability of this penalty
519- if (active==new_active){
520- unstable=0;
521- }
522- else{
523- unstable=1;
524- if(penalty_lock)zigzag_counter++;
525- }
526-
527- /*If penalty keeps zigzagging more than penalty_lock times: */
528- if(penalty_lock){
529- if(zigzag_counter>penalty_lock){
530- unstable=0;
531- active=1;
532- }
533- }
534-
535- //Set penalty flag
536- active=new_active;
537-
538- //*Assign output pointers:*/
539- *punstable=unstable;
540-}
541-/*}}}*/
542-/*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/
543-ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){
544-
545- IssmDouble penalty_factor;
546-
547- /*Initialize Element matrix and return if necessary*/
548- if(!this->active) return NULL;
549- ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);
550-
551- /*recover parameters: */
552- parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
553-
554- Ke->values[0]=kmax*pow(10.,penalty_factor);
555-
556- /*Clean up and return*/
557- return Ke;
558-}
559-/*}}}*/
560-/*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/
561-ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){
562-
563- IssmDouble penalty_factor;
564- IssmDouble max_damage;
565-
566- /*Initialize Element matrix and return if necessary*/
567- if(!this->active) return NULL;
568- ElementVector* pe=new ElementVector(&node,1,this->parameters);
569-
570- //Recover our data:
571- parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);
572- parameters->FindParam(&max_damage,DamageMaxDamageEnum);
573-
574- //right hand side penalizes to max_damage
575- pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage;
576-
577- /*Clean up and return*/
578- return pe;
579-}
580-/*}}}*/
581 /*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/
582 ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){
583
584Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.h
585===================================================================
586--- ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17515)
587+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17516)
588@@ -86,9 +86,6 @@
589 ElementVector* PenaltyCreatePVectorThermal(IssmDouble kmax);
590 ElementVector* PenaltyCreatePVectorMelting(IssmDouble kmax);
591 void ConstraintActivateThermal(int* punstable);
592- ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax);
593- ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax);
594- void ConstraintActivateDamageEvolution(int* punstable);
595 ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax);
596 ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax);
597 void ConstraintActivateHydrologyDCInefficient(int* punstable);
598Index: ../trunk-jpl/src/c/classes/Elements/Element.h
599===================================================================
600--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17515)
601+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17516)
602@@ -37,6 +37,8 @@
603 class Element: public Object,public Update{
604
605 public:
606+ int id;
607+ int sid;
608 Inputs *inputs;
609 Node **nodes;
610 Vertex **vertices;
611@@ -54,8 +56,14 @@
612 /* bool AllActive(void); */
613 /* bool AnyActive(void); */
614 void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
615+ void Echo();
616+ void DeepEcho();
617 void DeleteMaterials(void);
618 IssmDouble Divergence(void);
619+ void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
620+ void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
621+ IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
622+ IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
623 void FindParam(bool* pvalue,int paramenum);
624 void FindParam(int* pvalue,int paramenum);
625 void FindParam(IssmDouble* pvalue,int paramenum);
626@@ -80,19 +88,26 @@
627 void GetVerticesCoordinates(IssmDouble** xyz_list);
628 void GetVerticesSidList(int* sidlist);
629 void GetVerticesConnectivityList(int* connectivitylist);
630+ bool HasNodeOnBed();
631+ bool HasNodeOnSurface();
632+ int Id();
633+ int Sid();
634 void InputChangeName(int enum_type,int enum_type_old);
635 void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
636 void InputUpdateFromConstant(IssmDouble constant, int name);
637 void InputUpdateFromConstant(int constant, int name);
638 void InputUpdateFromConstant(bool constant, int name);
639+ bool IsIceInElement();
640 bool IsInput(int name);
641 bool IsFloating();
642 ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum);
643 ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum);
644 ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
645+ IssmDouble PureIceEnthalpy(IssmDouble pressure);
646 void ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
647 void ResultToVector(Vector<IssmDouble>* vector,int output_enum);
648 void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum);
649+ void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
650 void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
651 void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input);
652 void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
653@@ -134,12 +149,8 @@
654 virtual IssmDouble CharacteristicLength(void)=0;
655 virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;
656 virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0;
657- virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
658 virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
659- virtual void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0;
660- virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
661- virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
662- virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0;
663+
664 virtual int FiniteElement(void)=0;
665 virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
666 virtual void NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
667@@ -153,7 +164,6 @@
668 virtual void NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
669 virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
670 virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
671- virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
672
673 virtual Element* GetUpperElement(void)=0;
674 virtual Element* GetLowerElement(void)=0;
675@@ -168,16 +178,11 @@
676 virtual void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum)=0;
677 virtual int GetNodeIndex(Node* node)=0;
678 virtual int GetNumberOfNodes(void)=0;
679- virtual int GetNumberOfNodesVelocity(void)=0;
680- virtual int GetNumberOfNodesPressure(void)=0;
681 virtual int GetNumberOfVertices(void)=0;
682
683- virtual int Id()=0;
684- virtual int Sid()=0;
685 virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0;
686 virtual bool IsOnBed()=0;
687 virtual bool IsOnSurface()=0;
688- virtual bool IsIceInElement()=0;
689 virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
690 virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
691 virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
692Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
693===================================================================
694--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17515)
695+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17516)
696@@ -94,64 +94,6 @@
697 /*}}}*/
698
699 /*Other*/
700-/*FUNCTION Tria::SetwiseNodeConnectivity{{{*/
701-void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
702-
703- /*Intermediaries*/
704- const int numnodes = this->NumberofNodes();
705-
706- /*Output */
707- int d_nz = 0;
708- int o_nz = 0;
709-
710- /*Loop over all nodes*/
711- for(int i=0;i<numnodes;i++){
712-
713- if(!flags[this->nodes[i]->Lid()]){
714-
715- /*flag current node so that no other element processes it*/
716- flags[this->nodes[i]->Lid()]=true;
717-
718- int counter=0;
719- while(flagsindices[counter]>=0) counter++;
720- flagsindices[counter]=this->nodes[i]->Lid();
721-
722- /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
723- switch(set2_enum){
724- case FsetEnum:
725- if(nodes[i]->indexing.fsize){
726- if(this->nodes[i]->IsClone())
727- o_nz += 1;
728- else
729- d_nz += 1;
730- }
731- break;
732- case GsetEnum:
733- if(nodes[i]->indexing.gsize){
734- if(this->nodes[i]->IsClone())
735- o_nz += 1;
736- else
737- d_nz += 1;
738- }
739- break;
740- case SsetEnum:
741- if(nodes[i]->indexing.ssize){
742- if(this->nodes[i]->IsClone())
743- o_nz += 1;
744- else
745- d_nz += 1;
746- }
747- break;
748- default: _error_("not supported");
749- }
750- }
751- }
752-
753- /*Assign output pointers: */
754- *pd_nz=d_nz;
755- *po_nz=o_nz;
756-}
757-/*}}}*/
758 /*FUNCTION Tria::AddBasalInput{{{*/
759 void Tria::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){
760
761@@ -291,35 +233,6 @@
762
763 }
764 /*}}}*/
765-/*FUNCTION Tria::DeepEcho{{{*/
766-void Tria::DeepEcho(void){
767-
768- _printf_("Tria:\n");
769- _printf_(" id: " << id << "\n");
770- if(nodes){
771- nodes[0]->DeepEcho();
772- nodes[1]->DeepEcho();
773- nodes[2]->DeepEcho();
774- }
775- else _printf_("nodes = NULL\n");
776-
777- if (material) material->DeepEcho();
778- else _printf_("material = NULL\n");
779-
780- if (matpar) matpar->DeepEcho();
781- else _printf_("matpar = NULL\n");
782-
783- _printf_(" parameters\n");
784- if (parameters) parameters->DeepEcho();
785- else _printf_("parameters = NULL\n");
786-
787- _printf_(" inputs\n");
788- if (inputs) inputs->DeepEcho();
789- else _printf_("inputs=NULL\n");
790-
791- return;
792-}
793-/*}}}*/
794 /*FUNCTION Tria::Delta18oParameterization{{{*/
795 void Tria::Delta18oParameterization(void){
796
797@@ -415,32 +328,6 @@
798 *hz=0.;
799 }
800 /*}}}*/
801-/*FUNCTION Tria::Echo{{{*/
802-void Tria::Echo(void){
803- _printf_("Tria:\n");
804- _printf_(" id: " << id << "\n");
805- if(nodes){
806- nodes[0]->Echo();
807- nodes[1]->Echo();
808- nodes[2]->Echo();
809- }
810- else _printf_("nodes = NULL\n");
811-
812- if (material) material->Echo();
813- else _printf_("material = NULL\n");
814-
815- if (matpar) matpar->Echo();
816- else _printf_("matpar = NULL\n");
817-
818- _printf_(" parameters\n");
819- if (parameters) parameters->Echo();
820- else _printf_("parameters = NULL\n");
821-
822- _printf_(" inputs\n");
823- if (inputs) inputs->Echo();
824- else _printf_("inputs=NULL\n");
825-}
826-/*}}}*/
827 /*FUNCTION Tria::FiniteElement{{{*/
828 int Tria::FiniteElement(void){
829 return this->element_type;
830@@ -921,16 +808,6 @@
831 return this->NumberofNodes();
832 }
833 /*}}}*/
834-/*FUNCTION Tria::GetNumberOfNodesPressure;{{{*/
835-int Tria::GetNumberOfNodesPressure(void){
836- return this->NumberofNodesPressure();
837-}
838-/*}}}*/
839-/*FUNCTION Tria::GetNumberOfNodesVelocity;{{{*/
840-int Tria::GetNumberOfNodesVelocity(void){
841- return this->NumberofNodesVelocity();
842-}
843-/*}}}*/
844 /*FUNCTION Tria::GetNumberOfVertices;{{{*/
845 int Tria::GetNumberOfVertices(void){
846 return NUMVERTICES;
847@@ -1009,20 +886,6 @@
848 return y;
849 }
850 /*}}}*/
851-/*FUNCTION Tria::Id {{{*/
852-int Tria::Id(){
853-
854- return id;
855-
856-}
857-/*}}}*/
858-/*FUNCTION Tria::Sid {{{*/
859-int Tria::Sid(){
860-
861- return sid;
862-
863-}
864-/*}}}*/
865 /*FUNCTION Tria::InputDepthAverageAtBase {{{*/
866 void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){
867
868@@ -1350,19 +1213,6 @@
869 }
870 }
871 /*}}}*/
872-/*FUNCTION Tria::HasNodeOnBed {{{*/
873-bool Tria::HasNodeOnBed(){
874-
875- IssmDouble values[NUMVERTICES];
876- IssmDouble sum;
877-
878- /*Retrieve all inputs and parameters*/
879- GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
880- sum = values[0]+values[1]+values[2];
881-
882- return sum>0.;
883-}
884-/*}}}*/
885 /*FUNCTION Tria::HasEdgeOnSurface {{{*/
886 bool Tria::HasEdgeOnSurface(){
887
888@@ -1385,19 +1235,6 @@
889 }
890 }
891 /*}}}*/
892-/*FUNCTION Tria::HasNodeOnSurface {{{*/
893-bool Tria::HasNodeOnSurface(){
894-
895- IssmDouble values[NUMVERTICES];
896- IssmDouble sum;
897-
898- /*Retrieve all inputs and parameters*/
899- GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
900- sum = values[0]+values[1]+values[2];
901-
902- return sum>0.;
903-}
904-/*}}}*/
905 /*FUNCTION Tria::EdgeOnBedIndices{{{*/
906 void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
907
908@@ -1541,20 +1378,6 @@
909 return new GaussTria(indices[0],indices[1],order);
910 }
911 /*}}}*/
912-/*FUNCTION Tria::IsIceInElement {{{*/
913-bool Tria::IsIceInElement(){
914-
915- /*Get levelset*/
916- IssmDouble ls[NUMVERTICES];
917- GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
918-
919- /*If the level set on one of the nodes is <0, ice is present in this element*/
920- if(ls[0]<0. || ls[1]<0. || ls[2]<0.)
921- return true;
922- else
923- return false;
924-}
925-/*}}}*/
926 /*FUNCTION Tria::NodalFunctions{{{*/
927 void Tria::NodalFunctions(IssmDouble* basis, Gauss* gauss){
928
929Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
930===================================================================
931--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17515)
932+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17516)
933@@ -31,18 +31,12 @@
934
935 public:
936
937- int id;
938- int sid;
939-
940 /*Tria constructors, destructors {{{*/
941 Tria(){};
942 Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels);
943 ~Tria();
944 /*}}}*/
945 /*Object virtual functions definitions:{{{ */
946- void Echo();
947- void DeepEcho();
948- int Id();
949 int ObjectEnum();
950 Object *copy();
951 /*}}}*/
952@@ -62,11 +56,8 @@
953 void ComputeSurfaceNormalVelocity();
954 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
955 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
956- void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
957 void Delta18oParameterization(void);
958 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
959- void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
960- void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
961 int FiniteElement(void);
962 Element* GetUpperElement(void){_error_("not implemented yet");};
963 Element* GetLowerElement(void){_error_("not implemented yet");};
964@@ -76,16 +67,11 @@
965 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
966 int GetNodeIndex(Node* node);
967 int GetNumberOfNodes(void);
968- int GetNumberOfNodesPressure(void);
969- int GetNumberOfNodesVelocity(void);
970 int GetNumberOfVertices(void);
971- int Sid();
972 bool IsOnBed();
973 bool IsOnSurface();
974 bool HasEdgeOnBed();
975- bool HasNodeOnBed();
976 bool HasEdgeOnSurface();
977- bool HasNodeOnSurface();
978 void EdgeOnSurfaceIndices(int* pindex1,int* pindex);
979 void EdgeOnBedIndices(int* pindex1,int* pindex);
980 int EdgeOnBedIndex();
981@@ -93,7 +79,6 @@
982 bool IsNodeOnShelfFromFlags(IssmDouble* flags);
983 int NumberofNodesVelocity(void);
984 int NumberofNodesPressure(void);
985- bool IsIceInElement();
986 void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
987 void GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
988 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
989@@ -110,7 +95,6 @@
990 Element* SpawnBasalElement(void);
991 Element* SpawnTopElement(void);
992 int VelocityInterpolation();
993- IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
994 int PressureInterpolation();
995 IssmDouble SurfaceArea(void);
996 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
997@@ -191,8 +175,6 @@
998 /*Tria specific routines:{{{*/
999 void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
1000 void AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
1001- IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
1002- IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
1003 IssmDouble GetArea(void);
1004 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
1005 int GetElementType(void);
1006Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
1007===================================================================
1008--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17515)
1009+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17516)
1010@@ -370,26 +370,6 @@
1011 this->inputs->Configure(parameters);
1012 }
1013 /*}}}*/
1014-/*FUNCTION Penta::DeepEcho{{{*/
1015-void Penta::DeepEcho(void){
1016-
1017- _printf_("Penta:\n");
1018- _printf_(" id: " << id << "\n");
1019- nodes[0]->DeepEcho();
1020- nodes[1]->DeepEcho();
1021- nodes[2]->DeepEcho();
1022- nodes[3]->DeepEcho();
1023- nodes[4]->DeepEcho();
1024- nodes[5]->DeepEcho();
1025- material->DeepEcho();
1026- matpar->DeepEcho();
1027- _printf_(" neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id() << "\n");
1028- _printf_(" parameters\n");
1029- parameters->DeepEcho();
1030- _printf_(" inputs\n");
1031- inputs->DeepEcho();
1032-}
1033-/*}}}*/
1034 /*FUNCTION Penta::Delta18oParameterization{{{*/
1035 void Penta::Delta18oParameterization(void){
1036 /*Are we on the base? If not, return*/
1037@@ -464,32 +444,6 @@
1038 delete gauss;
1039 }
1040 /*}}}*/
1041-/*FUNCTION Penta::Echo{{{*/
1042-
1043-void Penta::Echo(void){
1044- this->DeepEcho();
1045-}
1046-/*}}}*/
1047-/*FUNCTION Penta::ThermalToEnthalpy{{{*/
1048-void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){
1049- matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
1050-}
1051-/*}}}*/
1052-/*FUNCTION Penta::EnthalpyToThermal{{{*/
1053-void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
1054- matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
1055-}
1056-/*}}}*/
1057-/*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/
1058-IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
1059- return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
1060-}
1061-/*}}}*/
1062-/*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/
1063-IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
1064- return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
1065-}
1066-/*}}}*/
1067 /*FUNCTION Penta::FiniteElement{{{*/
1068 int Penta::FiniteElement(void){
1069 return this->element_type;
1070@@ -847,16 +801,6 @@
1071 return this->NumberofNodes();
1072 }
1073 /*}}}*/
1074-/*FUNCTION Penta::GetNumberOfNodesPressure;{{{*/
1075-int Penta::GetNumberOfNodesPressure(void){
1076- return this->NumberofNodesPressure();
1077-}
1078-/*}}}*/
1079-/*FUNCTION Penta::GetNumberOfNodesVelocity;{{{*/
1080-int Penta::GetNumberOfNodesVelocity(void){
1081- return this->NumberofNodesVelocity();
1082-}
1083-/*}}}*/
1084 /*FUNCTION Penta::GetNumberOfVertices;{{{*/
1085 int Penta::GetNumberOfVertices(void){
1086 return NUMVERTICES;
1087@@ -1161,18 +1105,6 @@
1088 *pxyz_zero= xyz_zero;
1089 }
1090 /*}}}*/
1091-/*FUNCTION Penta::Sid {{{*/
1092-int Penta::Sid(){
1093-
1094- return sid;
1095-
1096-}
1097-/*}}}*/
1098-/*FUNCTION Penta::Id {{{*/
1099-int Penta::Id(void){
1100- return id;
1101-}
1102-/*}}}*/
1103 /*FUNCTION Penta::InputDepthAverageAtBase{{{*/
1104 void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
1105
1106@@ -1642,20 +1574,6 @@
1107
1108 }
1109 /*}}}*/
1110-/*FUNCTION Penta::IsIceInElement {{{*/
1111-bool Penta::IsIceInElement(){
1112-
1113- /*Get levelset*/
1114- IssmDouble ls[NUMVERTICES];
1115- GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
1116-
1117- /*If the level set one one of the nodes is <0, ice is present in this element*/
1118- if(ls[0]<0. || ls[1]<0. || ls[2]<0.)
1119- return true;
1120- else
1121- return false;
1122-}
1123-/*}}}*/
1124 /*FUNCTION Penta::MinEdgeLength{{{*/
1125 IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){
1126 /*Return the minimum lenght of the nine egdes of the penta*/
1127@@ -1812,12 +1730,11 @@
1128 /*FUNCTION Penta::NormalBase {{{*/
1129 void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
1130
1131- int i;
1132 IssmDouble v13[3],v23[3];
1133 IssmDouble normal[3];
1134 IssmDouble normal_norm;
1135
1136- for (i=0;i<3;i++){
1137+ for(int i=0;i<3;i++){
1138 v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
1139 v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
1140 }
1141@@ -1825,7 +1742,7 @@
1142 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
1143 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
1144 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
1145- normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) );
1146+ normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
1147
1148 /*Bed normal is opposite to surface normal*/
1149 bed_normal[0]=-normal[0]/normal_norm;
1150@@ -1918,12 +1835,6 @@
1151 delete gauss;
1152 }
1153 /*}}}*/
1154-/*FUNCTION Penta::PureIceEnthalpy{{{*/
1155-IssmDouble Penta::PureIceEnthalpy(IssmDouble pressure){
1156-
1157- return this->matpar->PureIceEnthalpy(pressure);
1158-}
1159-/*}}}*/
1160 /*FUNCTION Penta::ReduceMatrices{{{*/
1161 void Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){
1162
1163@@ -2071,77 +1982,6 @@
1164 this->element_type=element_type_in;
1165 }
1166 /*}}}*/
1167-/*FUNCTION Penta::SetwiseNodeConnectivity{{{*/
1168-void Penta::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
1169-
1170- /*Intermediaries*/
1171- const int numnodes = this->NumberofNodes();
1172-
1173- /*Output */
1174- int d_nz = 0;
1175- int o_nz = 0;
1176-
1177- /*Loop over all nodes*/
1178- for(int i=0;i<numnodes;i++){
1179-
1180- if(!flags[this->nodes[i]->Lid()]){
1181-
1182- /*flag current node so that no other element processes it*/
1183- flags[this->nodes[i]->Lid()]=true;
1184-
1185- int counter=0;
1186- while(flagsindices[counter]>=0) counter++;
1187- flagsindices[counter]=this->nodes[i]->Lid();
1188-
1189- /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
1190- switch(set2_enum){
1191- case FsetEnum:
1192- if(nodes[i]->indexing.fsize){
1193- if(this->nodes[i]->IsClone())
1194- o_nz += 1;
1195- else
1196- d_nz += 1;
1197- }
1198- break;
1199- case GsetEnum:
1200- if(nodes[i]->indexing.gsize){
1201- if(this->nodes[i]->IsClone())
1202- o_nz += 1;
1203- else
1204- d_nz += 1;
1205- }
1206- break;
1207- case SsetEnum:
1208- if(nodes[i]->indexing.ssize){
1209- if(this->nodes[i]->IsClone())
1210- o_nz += 1;
1211- else
1212- d_nz += 1;
1213- }
1214- break;
1215- default: _error_("not supported");
1216- }
1217- }
1218- }
1219-
1220- /*Special case: 2d/3d coupling, the node of this element might be connected
1221- *to the basal element*/
1222- int analysis_type,approximation,numlayers;
1223- parameters->FindParam(&analysis_type,AnalysisTypeEnum);
1224- if(analysis_type==StressbalanceAnalysisEnum){
1225- inputs->GetInputValue(&approximation,ApproximationEnum);
1226- if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
1227- parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
1228- o_nz += numlayers*3;
1229- d_nz += numlayers*3;
1230- }
1231- }
1232-
1233- /*Assign output pointers: */
1234- *pd_nz=d_nz;
1235- *po_nz=o_nz;
1236-}
1237-/*}}}*/
1238 /*FUNCTION Penta::SpawnTria {{{*/
1239 Tria* Penta::SpawnTria(int index1,int index2,int index3){
1240
1241Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
1242===================================================================
1243--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17515)
1244+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17516)
1245@@ -31,9 +31,6 @@
1246
1247 public:
1248
1249- int id;
1250- int sid;
1251-
1252 Penta **verticalneighbors; // 2 neighbors: first one under, second one above
1253
1254 /*Penta constructors and destructor: {{{*/
1255@@ -43,10 +40,7 @@
1256 /*}}}*/
1257 /*Object virtual functions definitions: {{{*/
1258 Object *copy();
1259- void DeepEcho();
1260- void Echo();
1261 int ObjectEnum();
1262- int Id();
1263 /*}}}*/
1264 /*Update virtual functions definitions: {{{*/
1265 void InputUpdateFromVector(IssmDouble* vector, int name, int type);
1266@@ -66,10 +60,7 @@
1267 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
1268 int FiniteElement(void);
1269 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
1270- void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
1271 void Delta18oParameterization(void);
1272- void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
1273- void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
1274 Penta* GetUpperPenta(void);
1275 Penta* GetLowerPenta(void);
1276 Penta* GetSurfacePenta(void);
1277@@ -82,8 +73,6 @@
1278 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
1279 int GetNodeIndex(Node* node);
1280 int GetNumberOfNodes(void);
1281- int GetNumberOfNodesPressure(void);
1282- int GetNumberOfNodesVelocity(void);
1283 int GetNumberOfVertices(void);
1284 void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
1285 IssmDouble GetXcoord(Gauss* gauss);
1286@@ -93,14 +82,12 @@
1287 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
1288 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
1289
1290- int Sid();
1291 void InputDepthAverageAtBase(int enum_type,int average_enum_type);
1292 void InputDuplicate(int original_enum,int new_enum);
1293 void InputScale(int enum_type,IssmDouble scale_factor);
1294 int NumberofNodesVelocity(void);
1295 int NumberofNodesPressure(void);
1296 int VelocityInterpolation();
1297- IssmDouble PureIceEnthalpy(IssmDouble pressure);
1298 int PressureInterpolation();
1299 bool IsZeroLevelset(int levelset_enum);
1300 bool IsIcefront(void);
1301@@ -188,8 +175,6 @@
1302 void NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
1303 void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list);
1304 ElementMatrix* CreateBasalMassMatrix(void);
1305- IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
1306- IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
1307 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
1308
1309 void GetVertexPidList(int* doflist);
1310@@ -207,7 +192,6 @@
1311 void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss);
1312 void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
1313 void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
1314- bool IsIceInElement(void);
1315 Gauss* NewGauss(void);
1316 Gauss* NewGauss(int order);
1317 Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
1318Index: ../trunk-jpl/src/c/classes/Elements/Seg.cpp
1319===================================================================
1320--- ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17515)
1321+++ ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17516)
1322@@ -64,64 +64,11 @@
1323 return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
1324 }
1325 /*}}}*/
1326-/*FUNCTION Seg::Echo{{{*/
1327-void Seg::Echo(void){
1328- _printf_("Seg:\n");
1329- _printf_(" id: " << id << "\n");
1330- if(nodes){
1331- nodes[0]->Echo();
1332- nodes[1]->Echo();
1333- }
1334- else _printf_("nodes = NULL\n");
1335-
1336- if (material) material->Echo();
1337- else _printf_("material = NULL\n");
1338-
1339- if (matpar) matpar->Echo();
1340- else _printf_("matpar = NULL\n");
1341-
1342- _printf_(" parameters\n");
1343- if (parameters) parameters->Echo();
1344- else _printf_("parameters = NULL\n");
1345-
1346- _printf_(" inputs\n");
1347- if (inputs) inputs->Echo();
1348- else _printf_("inputs=NULL\n");
1349-}
1350-/*}}}*/
1351 /*FUNCTION Seg::FiniteElement{{{*/
1352 int Seg::FiniteElement(void){
1353 return this->element_type;
1354 }
1355 /*}}}*/
1356-/*FUNCTION Seg::DeepEcho{{{*/
1357-void Seg::DeepEcho(void){
1358-
1359- _printf_("Seg:\n");
1360- _printf_(" id: " << id << "\n");
1361- if(nodes){
1362- nodes[0]->DeepEcho();
1363- nodes[1]->DeepEcho();
1364- }
1365- else _printf_("nodes = NULL\n");
1366-
1367- if (material) material->DeepEcho();
1368- else _printf_("material = NULL\n");
1369-
1370- if (matpar) matpar->DeepEcho();
1371- else _printf_("matpar = NULL\n");
1372-
1373- _printf_(" parameters\n");
1374- if (parameters) parameters->DeepEcho();
1375- else _printf_("parameters = NULL\n");
1376-
1377- _printf_(" inputs\n");
1378- if (inputs) inputs->DeepEcho();
1379- else _printf_("inputs=NULL\n");
1380-
1381- return;
1382-}
1383-/*}}}*/
1384 /*FUNCTION Seg::ObjectEnum{{{*/
1385 int Seg::ObjectEnum(void){
1386
1387@@ -129,14 +76,7 @@
1388
1389 }
1390 /*}}}*/
1391-/*FUNCTION Seg::Id {{{*/
1392-int Seg::Id(){
1393
1394- return id;
1395-
1396-}
1397-/*}}}*/
1398-
1399 void Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
1400
1401 /* Intermediaries */
1402@@ -185,20 +125,6 @@
1403 *pxyz_list = xyz_list;
1404
1405 }/*}}}*/
1406-/*FUNCTION Seg::IsIceInElement {{{*/
1407-bool Seg::IsIceInElement(){
1408-
1409- /*Get levelset*/
1410- IssmDouble ls[NUMVERTICES];
1411- GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
1412-
1413- /*If the level set on one of the nodes is <0, ice is present in this element*/
1414- if(ls[0]<0. || ls[1]<0.)
1415- return true;
1416- else
1417- return false;
1418-}
1419-/*}}}*/
1420 bool Seg::IsIcefront(void){/*{{{*/
1421
1422 bool isicefront;
1423Index: ../trunk-jpl/src/c/classes/Elements/Tetra.cpp
1424===================================================================
1425--- ../trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17515)
1426+++ ../trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17516)
1427@@ -51,80 +51,16 @@
1428 }
1429 /*}}}*/
1430
1431-/*FUNCTION Tetra::Echo{{{*/
1432-void Tetra::Echo(void){
1433- _printf_("Tetra:\n");
1434- _printf_(" id: " << id << "\n");
1435- if(nodes){
1436- nodes[0]->Echo();
1437- nodes[1]->Echo();
1438- nodes[2]->Echo();
1439- nodes[3]->Echo();
1440- }
1441- else _printf_("nodes = NULL\n");
1442-
1443- if (material) material->Echo();
1444- else _printf_("material = NULL\n");
1445-
1446- if (matpar) matpar->Echo();
1447- else _printf_("matpar = NULL\n");
1448-
1449- _printf_(" parameters\n");
1450- if (parameters) parameters->Echo();
1451- else _printf_("parameters = NULL\n");
1452-
1453- _printf_(" inputs\n");
1454- if (inputs) inputs->Echo();
1455- else _printf_("inputs=NULL\n");
1456-}
1457-/*}}}*/
1458 /*FUNCTION Tetra::FiniteElement{{{*/
1459 int Tetra::FiniteElement(void){
1460 return this->element_type;
1461-}
1462-/*}}}*/
1463-/*FUNCTION Tetra::DeepEcho{{{*/
1464-void Tetra::DeepEcho(void){
1465-
1466- _printf_("Tetra:\n");
1467- _printf_(" id: " << id << "\n");
1468- if(nodes){
1469- nodes[0]->DeepEcho();
1470- nodes[1]->DeepEcho();
1471- nodes[2]->DeepEcho();
1472- nodes[3]->DeepEcho();
1473- }
1474- else _printf_("nodes = NULL\n");
1475-
1476- if (material) material->DeepEcho();
1477- else _printf_("material = NULL\n");
1478-
1479- if (matpar) matpar->DeepEcho();
1480- else _printf_("matpar = NULL\n");
1481-
1482- _printf_(" parameters\n");
1483- if (parameters) parameters->DeepEcho();
1484- else _printf_("parameters = NULL\n");
1485-
1486- _printf_(" inputs\n");
1487- if (inputs) inputs->DeepEcho();
1488- else _printf_("inputs=NULL\n");
1489-
1490- return;
1491-}
1492-/*}}}*/
1493+} /*}}}*/
1494 /*FUNCTION Tetra::ObjectEnum{{{*/
1495 int Tetra::ObjectEnum(void){
1496
1497 return TetraEnum;
1498
1499-}
1500-/*}}}*/
1501-/*FUNCTION Tetra::Id {{{*/
1502-int Tetra::Id(){
1503- return id;
1504-}
1505-/*}}}*/
1506+}/*}}}*/
1507
1508 /*FUNCTION Tetra::AddInput{{{*/
1509 void Tetra::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){
1510@@ -239,16 +175,6 @@
1511 return NUMVERTICES;
1512 }
1513 /*}}}*/
1514-/*FUNCTION Tetra::GetNumberOfNodesPressure THIS ONE (and corresponding TetraRef function){{{*/
1515-int Tetra::GetNumberOfNodesPressure(void){
1516- return this->NumberofNodesPressure();
1517-}
1518-/*}}}*/
1519-/*FUNCTION Tetra::GetNumberOfNodesVelocity;{{{*/
1520-int Tetra::GetNumberOfNodesVelocity(void){
1521- return this->NumberofNodesVelocity();
1522-}
1523-/*}}}*/
1524 /*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
1525 void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
1526
1527@@ -340,19 +266,6 @@
1528 }
1529 }
1530 /*}}}*/
1531-/*FUNCTION Tetra::HasNodeOnBed THIS ONE{{{*/
1532-bool Tetra::HasNodeOnBed(){
1533-
1534- IssmDouble values[NUMVERTICES];
1535- IssmDouble sum;
1536-
1537- /*Retrieve all inputs and parameters*/
1538- GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
1539- sum = values[0]+values[1]+values[2]+values[3];
1540-
1541- return sum>0.;
1542-}
1543-/*}}}*/
1544 /*FUNCTION Tetra::InputUpdateFromIoModel {{{*/
1545 void Tetra::InputUpdateFromIoModel(int index,IoModel* iomodel){
1546
1547@@ -484,21 +397,6 @@
1548 xDelete<int>(doflist);
1549 }
1550 /*}}}*/
1551-/*FUNCTION Tetra::IsIceInElement THIS ONE{{{*/
1552-bool Tetra::IsIceInElement(){
1553-
1554- /*Get levelset*/
1555- IssmDouble ls[NUMVERTICES];
1556- GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
1557-
1558- /*If the level set on one of the nodes is <0, ice is present in this element*/
1559- for(int i=0;i<NUMVERTICES;i++){
1560- if(ls[i]<0.) return true;
1561- }
1562-
1563- return false;
1564-}
1565-/*}}}*/
1566 /*FUNCTION Tetra::IsOnBed {{{*/
1567 bool Tetra::IsOnBed(){
1568 return HasFaceOnBed();
1569@@ -650,6 +548,51 @@
1570 for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
1571 }
1572 /*}}}*/
1573+/*FUNCTION Tetra::NormalBase (THIS ONE){{{*/
1574+void Tetra::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
1575+
1576+ IssmDouble v13[3],v23[3];
1577+ IssmDouble normal[3];
1578+ IssmDouble normal_norm;
1579+
1580+ for(int i=0;i<3;i++){
1581+ v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
1582+ v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
1583+ }
1584+
1585+ normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
1586+ normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
1587+ normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
1588+ normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
1589+
1590+ /*Bed normal is opposite to surface normal*/
1591+ bed_normal[0]=-normal[0]/normal_norm;
1592+ bed_normal[1]=-normal[1]/normal_norm;
1593+ bed_normal[2]=-normal[2]/normal_norm;
1594+}
1595+/*}}}*/
1596+/*FUNCTION Tetra::NormalTop (THIS ONE){{{*/
1597+void Tetra::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){
1598+
1599+ IssmDouble v13[3],v23[3];
1600+ IssmDouble normal[3];
1601+ IssmDouble normal_norm;
1602+
1603+ for(int i=0;i<3;i++){
1604+ v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i];
1605+ v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i];
1606+ }
1607+
1608+ normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
1609+ normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
1610+ normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
1611+ normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]);
1612+
1613+ top_normal[0]=normal[0]/normal_norm;
1614+ top_normal[1]=normal[1]/normal_norm;
1615+ top_normal[2]=normal[2]/normal_norm;
1616+}
1617+/*}}}*/
1618 /*FUNCTION Tetra::NumberofNodesPressure{{{*/
1619 int Tetra::NumberofNodesPressure(void){
1620 return TetraRef::NumberofNodesPressure();
1621@@ -665,7 +608,8 @@
1622
1623 if(pe){
1624 if(this->element_type==MINIcondensedEnum){
1625- _error_("Not implemented");
1626+ int indices[3]={12,13,14};
1627+ pe->StaticCondensation(Ke,3,&indices[0]);
1628 }
1629 else if(this->element_type==P1bubblecondensedEnum){
1630 int size = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
1631@@ -680,7 +624,8 @@
1632
1633 if(Ke){
1634 if(this->element_type==MINIcondensedEnum){
1635- _error_("Not implemented");
1636+ int indices[3]={12,13,14};
1637+ Ke->StaticCondensation(3,&indices[0]);
1638 }
1639 else if(this->element_type==P1bubblecondensedEnum){
1640 int size = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
1641@@ -770,71 +715,6 @@
1642
1643 }
1644 /*}}}*/
1645-/*FUNCTION Tetra::SetwiseNodeConnectivity THIS ONE (take from Penta.cpp){{{*/
1646-void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){
1647-
1648- /*Intermediaries*/
1649- const int numnodes = this->NumberofNodes();
1650-
1651- /*Output */
1652- int d_nz = 0;
1653- int o_nz = 0;
1654-
1655- /*Loop over all nodes*/
1656- for(int i=0;i<numnodes;i++){
1657-
1658- if(!flags[this->nodes[i]->Lid()]){
1659-
1660- /*flag current node so that no other element processes it*/
1661- flags[this->nodes[i]->Lid()]=true;
1662-
1663- int counter=0;
1664- while(flagsindices[counter]>=0) counter++;
1665- flagsindices[counter]=this->nodes[i]->Lid();
1666-
1667- /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
1668- switch(set2_enum){
1669- case FsetEnum:
1670- if(nodes[i]->indexing.fsize){
1671- if(this->nodes[i]->IsClone())
1672- o_nz += 1;
1673- else
1674- d_nz += 1;
1675- }
1676- break;
1677- case GsetEnum:
1678- if(nodes[i]->indexing.gsize){
1679- if(this->nodes[i]->IsClone())
1680- o_nz += 1;
1681- else
1682- d_nz += 1;
1683- }
1684- break;
1685- case SsetEnum:
1686- if(nodes[i]->indexing.ssize){
1687- if(this->nodes[i]->IsClone())
1688- o_nz += 1;
1689- else
1690- d_nz += 1;
1691- }
1692- break;
1693- default: _error_("not supported");
1694- }
1695- }
1696- }
1697-
1698- /*Assign output pointers: */
1699- *pd_nz=d_nz;
1700- *po_nz=o_nz;
1701-}
1702-/*}}}*/
1703-/*FUNCTION Tetra::Sid THIS ONE{{{*/
1704-int Tetra::Sid(){
1705-
1706- return sid;
1707-
1708-}
1709-/*}}}*/
1710 /*FUNCTION Tetra::SpawnBasalElement{{{*/
1711 Element* Tetra::SpawnBasalElement(void){
1712
1713Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
1714===================================================================
1715--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17515)
1716+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17516)
1717@@ -29,18 +29,12 @@
1718
1719 public:
1720
1721- int id;
1722- int sid;
1723-
1724 /*Seg constructors, destructors {{{*/
1725 Seg(){};
1726 Seg(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels);
1727 ~Seg();
1728 /*}}}*/
1729 /*Object virtual functions definitions:{{{ */
1730- void Echo();
1731- void DeepEcho();
1732- int Id();
1733 int ObjectEnum();
1734 Object *copy();
1735 /*}}}*/
1736@@ -63,13 +57,8 @@
1737 void ComputeStressTensor(){_error_("not implemented yet");};
1738 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
1739 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
1740- void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};
1741 void Delta18oParameterization(void){_error_("not implemented yet");};
1742 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
1743- void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
1744- void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
1745- IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
1746- IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
1747 int FiniteElement(void);
1748 Element* GetUpperElement(void){_error_("not implemented yet");};
1749 Element* GetLowerElement(void){_error_("not implemented yet");};
1750@@ -77,13 +66,10 @@
1751 Element* GetBasalElement(void){_error_("not implemented yet");};
1752 int GetNodeIndex(Node* node){_error_("not implemented yet");};
1753 int GetNumberOfNodes(void);
1754- int GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
1755- int GetNumberOfNodesPressure(void){_error_("not implemented yet");};
1756 int GetNumberOfVertices(void);
1757 void GetVerticesCoordinates(IssmDouble** pxyz_list);
1758 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
1759 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
1760- int Sid(){_error_("not implemented yet");};
1761 bool IsOnBed(){_error_("not implemented yet");};
1762 bool IsOnSurface(){_error_("not implemented yet");};
1763 bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
1764@@ -101,7 +87,6 @@
1765 void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1766 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1767 void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1768- bool IsIceInElement();
1769 void NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
1770 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
1771 void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
1772@@ -110,7 +95,6 @@
1773 Element* SpawnBasalElement(void){_error_("not implemented yet");};
1774 Element* SpawnTopElement(void){_error_("not implemented yet");};
1775 IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
1776- IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
1777 int PressureInterpolation(void){_error_("not implemented yet");};
1778 void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
1779 void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1780Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
1781===================================================================
1782--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17515)
1783+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17516)
1784@@ -29,18 +29,12 @@
1785
1786 public:
1787
1788- int id;
1789- int sid;
1790-
1791 /*Tetra constructors, destructors {{{*/
1792 Tetra(){};
1793 Tetra(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels);
1794 ~Tetra();
1795 /*}}}*/
1796 /*Object virtual functions definitions:{{{ */
1797- void Echo();
1798- void DeepEcho();
1799- int Id();
1800 int ObjectEnum();
1801 Object *copy();
1802 /*}}}*/
1803@@ -63,13 +57,8 @@
1804 void ComputeStressTensor(){_error_("not implemented yet");};
1805 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
1806 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
1807- void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
1808 void Delta18oParameterization(void){_error_("not implemented yet");};
1809 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
1810- void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
1811- void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
1812- IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
1813- IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
1814 void FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3);
1815 void FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3);
1816 void FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3);
1817@@ -80,15 +69,11 @@
1818 Element* GetBasalElement(void){_error_("not implemented yet");};
1819 int GetNodeIndex(Node* node){_error_("not implemented yet");};
1820 int GetNumberOfNodes(void);
1821- int GetNumberOfNodesVelocity(void);
1822- int GetNumberOfNodesPressure(void);
1823 int GetNumberOfVertices(void);
1824 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
1825 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
1826 bool HasFaceOnBed();
1827 bool HasFaceOnSurface();
1828- bool HasNodeOnBed();
1829- int Sid();
1830 bool IsOnBed();
1831 bool IsOnSurface();
1832 bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
1833@@ -106,17 +91,15 @@
1834 void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1835 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1836 void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
1837- bool IsIceInElement();
1838 void NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
1839- void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
1840- void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
1841+ void NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
1842+ void NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
1843 int NumberofNodesVelocity(void);
1844 int NumberofNodesPressure(void);
1845 Element* SpawnBasalElement(void);
1846 Element* SpawnTopElement(void);
1847 Tria* SpawnTria(int index1,int index2,int index3);
1848 IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
1849- IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
1850 int PressureInterpolation(void);
1851 void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
1852 void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
1853@@ -227,4 +210,4 @@
1854 int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){_error_("not implemented yet");};
1855 /*}}}*/
1856 };
1857-#endif /* _SEG_H */
1858+#endif /* _TETRA_H_*/
1859Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
1860===================================================================
1861--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17515)
1862+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17516)
1863@@ -16,6 +16,8 @@
1864
1865 /*Constructors/destructor/copy*/
1866 Element::Element(){/*{{{*/
1867+ this->id = -1;
1868+ this->sid = -1;
1869 this->inputs = NULL;
1870 this->nodes = NULL;
1871 this->vertices = NULL;
1872@@ -129,6 +131,71 @@
1873 void Element::DeleteMaterials(void){/*{{{*/
1874 delete this->material;
1875 }/*}}}*/
1876+void Element::DeepEcho(void){/*{{{*/
1877+
1878+ _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
1879+ _printf_(" id : "<<this->id <<"\n");
1880+ _printf_(" sid: "<<this->sid<<"\n");
1881+ if(vertices){
1882+ int numvertices = this->GetNumberOfVertices();
1883+ for(int i=0;i<numvertices;i++) vertices[i]->Echo();
1884+ }
1885+ else _printf_("vertices = NULL\n");
1886+
1887+ if(nodes){
1888+ int numnodes = this->GetNumberOfNodes();
1889+ for(int i=0;i<numnodes;i++) nodes[i]->DeepEcho();
1890+ }
1891+ else _printf_("nodes = NULL\n");
1892+
1893+ if (material) material->DeepEcho();
1894+ else _printf_("material = NULL\n");
1895+
1896+ if (matpar) matpar->DeepEcho();
1897+ else _printf_("matpar = NULL\n");
1898+
1899+ _printf_(" parameters\n");
1900+ if (parameters) parameters->DeepEcho();
1901+ else _printf_("parameters = NULL\n");
1902+
1903+ _printf_(" inputs\n");
1904+ if (inputs) inputs->DeepEcho();
1905+ else _printf_("inputs=NULL\n");
1906+
1907+ return;
1908+}
1909+/*}}}*/
1910+void Element::Echo(void){/*{{{*/
1911+ _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
1912+ _printf_(" id : "<<this->id <<"\n");
1913+ _printf_(" sid: "<<this->sid<<"\n");
1914+ if(vertices){
1915+ int numvertices = this->GetNumberOfVertices();
1916+ for(int i=0;i<numvertices;i++) vertices[i]->Echo();
1917+ }
1918+ else _printf_("vertices = NULL\n");
1919+
1920+ if(nodes){
1921+ int numnodes = this->GetNumberOfNodes();
1922+ for(int i=0;i<numnodes;i++) nodes[i]->Echo();
1923+ }
1924+ else _printf_("nodes = NULL\n");
1925+
1926+ if (material) material->Echo();
1927+ else _printf_("material = NULL\n");
1928+
1929+ if (matpar) matpar->Echo();
1930+ else _printf_("matpar = NULL\n");
1931+
1932+ _printf_(" parameters\n");
1933+ if (parameters) parameters->Echo();
1934+ else _printf_("parameters = NULL\n");
1935+
1936+ _printf_(" inputs\n");
1937+ if (inputs) inputs->Echo();
1938+ else _printf_("inputs=NULL\n");
1939+}
1940+/*}}}*/
1941 IssmDouble Element::Divergence(void){/*{{{*/
1942 /*Compute element divergence*/
1943
1944@@ -161,6 +228,18 @@
1945 delete gauss;
1946 return divergence;
1947 }/*}}}*/
1948+void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
1949+ matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
1950+}/*}}}*/
1951+void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
1952+ matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
1953+}/*}}}*/
1954+IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
1955+ return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
1956+}/*}}}*/
1957+IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
1958+ return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
1959+}/*}}}*/
1960 void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
1961 this->parameters->FindParam(pvalue,paramenum);
1962 }/*}}}*/
1963@@ -199,7 +278,7 @@
1964 void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/
1965
1966 /*Fetch number of nodes and dof for this finite element*/
1967- int numnodes = this->GetNumberOfNodesVelocity();
1968+ int numnodes = this->NumberofNodesVelocity();
1969
1970 /*First, figure out size of doflist and create it: */
1971 int numberofdofs=0;
1972@@ -222,8 +301,8 @@
1973 void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
1974
1975 /*Fetch number of nodes and dof for this finite element*/
1976- int vnumnodes = this->GetNumberOfNodesVelocity();
1977- int pnumnodes = this->GetNumberOfNodesPressure();
1978+ int vnumnodes = this->NumberofNodesVelocity();
1979+ int pnumnodes = this->NumberofNodesPressure();
1980
1981 /*First, figure out size of doflist and create it: */
1982 int numberofdofs=0;
1983@@ -387,7 +466,7 @@
1984
1985 _assert_(pvalue);
1986
1987- int numnodes = this->GetNumberOfNodesVelocity();
1988+ int numnodes = this->NumberofNodesVelocity();
1989 Input *input = this->GetInput(enumtype);
1990 if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
1991
1992@@ -469,6 +548,18 @@
1993 for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
1994 }
1995 /*}}}*/
1996+bool Element::HasNodeOnBed(){/*{{{*/
1997+ return (this->inputs->Max(MeshVertexonbedEnum)>0.);
1998+}/*}}}*/
1999+bool Element::HasNodeOnSurface(){/*{{{*/
2000+ return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.);
2001+}/*}}}*/
2002+int Element::Id(){/*{{{*/
2003+
2004+ return this->id;
2005+
2006+}
2007+/*}}}*/
2008 void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/
2009 this->inputs->ChangeEnum(original_enum,new_enum);
2010 }
2011@@ -586,6 +677,10 @@
2012
2013 return shelf;
2014 }/*}}}*/
2015+bool Element::IsIceInElement(){/*{{{*/
2016+ return (this->inputs->Min(MaskIceLevelsetEnum)<0.);
2017+}
2018+/*}}}*/
2019 bool Element::IsInput(int name){/*{{{*/
2020 if (
2021 name==ThicknessEnum ||
2022@@ -674,6 +769,9 @@
2023 return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
2024 }
2025 /*}}}*/
2026+IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
2027+ return this->matpar->PureIceEnthalpy(pressure);
2028+}/*}}}*/
2029 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/
2030
2031 Input* input=this->inputs->GetInput(output_enum);
2032@@ -772,6 +870,82 @@
2033 input->ResultToPatch(values,nodesperelement,this->Sid());
2034
2035 } /*}}}*/
2036+void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
2037+
2038+ /*Intermediaries*/
2039+ const int numnodes = this->GetNumberOfNodes();
2040+
2041+ /*Output */
2042+ int d_nz = 0;
2043+ int o_nz = 0;
2044+
2045+ /*Loop over all nodes*/
2046+ for(int i=0;i<numnodes;i++){
2047+
2048+ if(!flags[this->nodes[i]->Lid()]){
2049+
2050+ /*flag current node so that no other element processes it*/
2051+ flags[this->nodes[i]->Lid()]=true;
2052+
2053+ int counter=0;
2054+ while(flagsindices[counter]>=0) counter++;
2055+ flagsindices[counter]=this->nodes[i]->Lid();
2056+
2057+ /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
2058+ switch(set2_enum){
2059+ case FsetEnum:
2060+ if(nodes[i]->indexing.fsize){
2061+ if(this->nodes[i]->IsClone())
2062+ o_nz += 1;
2063+ else
2064+ d_nz += 1;
2065+ }
2066+ break;
2067+ case GsetEnum:
2068+ if(nodes[i]->indexing.gsize){
2069+ if(this->nodes[i]->IsClone())
2070+ o_nz += 1;
2071+ else
2072+ d_nz += 1;
2073+ }
2074+ break;
2075+ case SsetEnum:
2076+ if(nodes[i]->indexing.ssize){
2077+ if(this->nodes[i]->IsClone())
2078+ o_nz += 1;
2079+ else
2080+ d_nz += 1;
2081+ }
2082+ break;
2083+ default: _error_("not supported");
2084+ }
2085+ }
2086+ }
2087+
2088+ /*Special case: 2d/3d coupling, the node of this element might be connected
2089+ *to the basal element*/
2090+ int analysis_type,approximation,numlayers;
2091+ parameters->FindParam(&analysis_type,AnalysisTypeEnum);
2092+ if(analysis_type==StressbalanceAnalysisEnum){
2093+ inputs->GetInputValue(&approximation,ApproximationEnum);
2094+ if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
2095+ parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
2096+ o_nz += numlayers*3;
2097+ d_nz += numlayers*3;
2098+ }
2099+ }
2100+
2101+ /*Assign output pointers: */
2102+ *pd_nz=d_nz;
2103+ *po_nz=o_nz;
2104+}
2105+/*}}}*/
2106+int Element::Sid(){/*{{{*/
2107+
2108+ return this->sid;
2109+
2110+}
2111+/*}}}*/
2112 IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
2113 _assert_(matpar);
2114 return this->matpar->TMeltingPoint(pressure);
Note: See TracBrowser for help on using the repository browser.