source: issm/oecreview/Archive/22819-23185/ISSM-23045-23046.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 62.9 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 23045)
4+++ ../trunk-jpl/src/c/modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (revision 23046)
5@@ -1,6 +1,6 @@
6 /*!\file: InterpFromGridToMeshx.cpp
7 * \brief "c" core code for interpolating values from a structured grid.
8- */
9+ */
10
11 #ifdef HAVE_CONFIG_H
12 #include <config.h>
13@@ -24,9 +24,6 @@
14 double* y=NULL;
15 int i;
16
17- // TEST
18- _printf_("\r N: "<<N<<" M:"<<M<<"\n");
19-
20 /*Some checks on arguments: */
21 if ((M<2) || (N<2) || (nods<=0)){
22 _error_("nothing to be done according to the dimensions of input matrices and vectors.");
23@@ -150,7 +147,7 @@
24 * | | |
25 * | | |
26 * y1 x---------+-----x Q21
27- * x1 x2
28+ * x1 x2
29 *
30 */
31 x1=x[n]; x2=x[n+1];
32@@ -224,7 +221,7 @@
33 double triangleinterp(double x1,double x2,double y1,double y2,double Q11,double Q12,double Q21,double Q22,double x,double y){
34 /*split the rectangle in 2 triangle and
35 * use Lagrange P1 interpolation
36- *
37+ *
38 * +3----+2,3' Q12----Q22
39 * | /| | /|
40 * | / | | / |
41@@ -281,7 +278,7 @@
42 _assert_(x2>x1 && y2>y1);
43 _assert_(x<=x2 && x>=x1 && y<=y2 && y>=y1);
44
45- return
46+ return
47 +Q11*(x2-x)*(y2-y)/((x2-x1)*(y2-y1))
48 +Q21*(x-x1)*(y2-y)/((x2-x1)*(y2-y1))
49 +Q12*(x2-x)*(y-y1)/((x2-x1)*(y2-y1))
50@@ -300,7 +297,7 @@
51 * | | |
52 * | | |
53 * y1 x--------x---------x Q21
54- * x1 xm x2
55+ * x1 xm x2
56 *
57 */
58 /*Checks*/
59Index: ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js
60===================================================================
61--- ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js (revision 23045)
62+++ ../trunk-jpl/src/wrappers/InterpFromGridToMesh/InterpFromGridToMesh.js (revision 23046)
63@@ -5,7 +5,7 @@
64 *
65 * Usage:
66 * var data_mesh=InterpFromGridToMesh(xIn,yIn,dataIn,xMeshIn,yMeshIn,defaultValue,interpType);\
67- *
68+ *
69 * xIn,yIn : coordinates of matrix data. (x and y must be in increasing order)
70 * dataIn : matrix holding the data to be interpolated onto the mesh
71 * xMeshIn,yMeshIn : coordinates of the points onto which we interpolate
72@@ -50,13 +50,13 @@
73 var y = {};
74 var yMesh = {};
75 //}}}
76-
77
78+
79 /*
80 Dynamic allocations
81 */
82 //{{{
83-
84+
85 /*
86 Input
87 */
88@@ -67,7 +67,7 @@
89 dxHeap = new Uint8Array(Module.HEAPU8.buffer, dxPtr, nx);
90 dxHeap.set(new Uint8Array(dx.buffer));
91 x = dxHeap.byteOffset;
92-
93+
94 dy = new Float64Array(yIn);
95 ny = dy.length * dy.BYTES_PER_ELEMENT;
96 dyPtr = Module._malloc(ny);
97@@ -74,7 +74,7 @@
98 dyHeap = new Uint8Array(Module.HEAPU8.buffer, dyPtr, ny);
99 dyHeap.set(new Uint8Array(dy.buffer));
100 y = dyHeap.byteOffset;
101-
102+
103 ddata = new Float64Array(dataIn);
104 ndata = ddata.length * ddata.BYTES_PER_ELEMENT;
105 ddataPtr = Module._malloc(ndata);
106@@ -81,7 +81,7 @@
107 ddataHeap = new Uint8Array(Module.HEAPU8.buffer, ddataPtr, ndata);
108 ddataHeap.set(new Uint8Array(ddata.buffer));
109 data = ddataHeap.byteOffset;
110-
111+
112 dxMesh = new Float64Array(xMeshIn);
113 nxMesh = dxMesh.length * dxMesh.BYTES_PER_ELEMENT;
114 dxMeshPtr = Module._malloc(nxMesh);
115@@ -88,7 +88,7 @@
116 dxMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dxMeshPtr, nxMesh);
117 dxMeshHeap.set(new Uint8Array(dxMesh.buffer));
118 xMesh = dxMeshHeap.byteOffset;
119-
120+
121 dyMesh = new Float64Array(yMeshIn);
122 nyMesh = dyMesh.length * dyMesh.BYTES_PER_ELEMENT;
123 dyMeshPtr = Module._malloc(nyMesh);
124@@ -95,12 +95,11 @@
125 dyMeshHeap = new Uint8Array(Module.HEAPU8.buffer, dyMeshPtr, nyMesh);
126 dyMeshHeap.set(new Uint8Array(dyMesh.buffer));
127 yMesh = dyMeshHeap.byteOffset;
128-
129+
130 nods = xMeshIn.length;
131 meshNumRows = xMeshIn.length;
132- console.log('InterpFromGridToMesh.js: meshNumRows => ' + meshNumRows);
133-
134-
135+
136+
137 /*
138 Retrieve interpolation type
139 */
140@@ -111,23 +110,23 @@
141 interpType = 'bilinear';
142 }
143 //}}}
144-
145+
146 /*
147 Output
148 */
149 pdataMesh = Module._malloc(4);
150 //}}}
151-
152+
153 //}}}
154-
155-
156+
157+
158 /*
159 Declare InterpFromGridToMesh module
160 */
161 //{{{
162 InterpFromGridToMeshModule = Module.cwrap(
163- 'InterpFromGridToMeshModule',
164- 'number',
165+ 'InterpFromGridToMeshModule',
166+ 'number',
167 [
168 'number', // output : pdataMesh
169 'number', // input : x
170@@ -144,20 +143,20 @@
171 ]
172 );
173 //}}}
174-
175-
176+
177+
178 /*
179 Call InterpFromGridToMesh module
180 */
181 //{{{
182 InterpFromGridToMeshModule(
183- pdataMesh,
184- x,
185- y,
186- data,
187- xMesh,
188- yMesh,
189- defaultValue,
190+ pdataMesh,
191+ x,
192+ y,
193+ data,
194+ xMesh,
195+ yMesh,
196+ defaultValue,
197 nods,
198 dataNumRowsIn,
199 dataNumColsIn,
200@@ -165,8 +164,8 @@
201 interpType
202 );
203 //}}}
204-
205-
206+
207+
208 /*
209 Dynamic copying from heap
210 */
211@@ -174,8 +173,8 @@
212 dataMeshPtr = Module.getValue(pdataMesh, 'i32');
213 dataMesh = Module.HEAPF64.slice(dataMeshPtr / 8, dataMeshPtr / 8 + nods);
214 //}}}
215-
216-
217+
218+
219 /*
220 Free resources
221 */
222@@ -182,7 +181,7 @@
223 //{{{
224 Module._free(pdataMesh);
225 //}}}
226-
227-
228+
229+
230 return dataMesh;
231 }
232Index: ../trunk-jpl/jenkins/macosx_pine-island
233===================================================================
234--- ../trunk-jpl/jenkins/macosx_pine-island (revision 23045)
235+++ ../trunk-jpl/jenkins/macosx_pine-island (revision 23046)
236@@ -4,9 +4,9 @@
237 #-------------------------------#
238
239 #MATLAB path
240-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
241+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
242
243-#ISSM CONFIGURATION
244+#ISSM CONFIGURATION
245 ISSM_CONFIG='--prefix=$ISSM_DIR \
246 --with-matlab-dir=$MATLAB_PATH \
247 --with-triangle-dir=$ISSM_DIR/externalpackages/triangle/install \
248@@ -53,6 +53,6 @@
249 #as follows: runme($MATLAB_NROPTIONS). The options must be understandable
250 #by Matlab and runme.m
251 #ex: "'id',[101 102 103]"
252-## FS
253+## FS
254 PYTHON_NROPTIONS=""
255 MATLAB_NROPTIONS="'exclude',[701,702,703,435,IdFromString('Dakota')]"
256Index: ../trunk-jpl/jenkins/macosx_pine-island_static
257===================================================================
258--- ../trunk-jpl/jenkins/macosx_pine-island_static (revision 23045)
259+++ ../trunk-jpl/jenkins/macosx_pine-island_static (revision 23046)
260@@ -4,9 +4,9 @@
261 #-------------------------------#
262
263 #MATLAB path
264-MATLAB_PATH="/Applications/MATLAB_R2017b.app/"
265+MATLAB_PATH="/Applications/MATLAB_R2015b.app/"
266
267-#ISSM CONFIGURATION
268+#ISSM CONFIGURATION
269 ISSM_CONFIG='--prefix=$ISSM_DIR \
270 --disable-static \
271 --enable-standalone-executables \
272@@ -22,7 +22,7 @@
273 --with-metis-dir=$ISSM_DIR/externalpackages/petsc/install \
274 --with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \
275 --with-math77-dir=$ISSM_DIR/externalpackages/math77/install \
276- --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin15/6.1.0/libgcc.a" \
277+ --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \
278 --with-numthreads=4'
279
280 #PYTHON and MATLAB testing
281@@ -60,6 +60,6 @@
282 #as follows: runme($MATLAB_NROPTIONS). The options must be understandable
283 #by Matlab and runme.m
284 #ex: "'id',[101 102 103]"
285-## bamg mesh FS
286+## bamg mesh FS
287 #PYTHON_NROPTIONS=""
288 #MATLAB_NROPTIONS="'exclude',[119,243,514,701,702,703,435,IdFromString('Dakota')]"
289Index: ../trunk-jpl/src/c/classes/FemModel.cpp
290===================================================================
291--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23045)
292+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23046)
293@@ -1463,7 +1463,7 @@
294
295 /*Figure out maximum across the cluster: */
296 ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
297- ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
298+ ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
299 maxvy=node_maxvy;
300
301 /*Assign output pointers:*/
302@@ -1487,7 +1487,7 @@
303
304 /*Figure out maximum across the cluster: */
305 ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
306- ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
307+ ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
308 maxvz=node_maxvz;
309
310 /*Assign output pointers:*/
311@@ -1511,7 +1511,7 @@
312
313 /*Figure out minimum across the cluster: */
314 ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
315- ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
316+ ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
317 minvel=node_minvel;
318
319 /*Assign output pointers:*/
320@@ -1535,7 +1535,7 @@
321
322 /*Figure out minimum across the cluster: */
323 ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
324- ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
325+ ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
326 minvx=node_minvx;
327
328 /*Assign output pointers:*/
329@@ -1559,7 +1559,7 @@
330
331 /*Figure out minimum across the cluster: */
332 ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
333- ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
334+ ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
335 minvy=node_minvy;
336
337 /*Assign output pointers:*/
338@@ -1583,7 +1583,7 @@
339
340 /*Figure out minimum across the cluster: */
341 ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
342- ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
343+ ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
344 minvz=node_minvz;
345
346 /*Assign output pointers:*/
347@@ -1628,7 +1628,7 @@
348 weights_input->GetInputValue(&weight,gauss,OmegaAbsGradientEnum);
349 omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
350
351- /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
352+ /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
353 //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
354 J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight;
355 }
356@@ -1964,7 +1964,7 @@
357 else{
358 IssmDouble* values = xNewZeroInit<IssmDouble>(nlines*ncols);
359 IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
360-
361+
362 /*Fill-in matrix*/
363 for(int j=0;j<elements->Size();j++){
364 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
365@@ -1973,7 +1973,7 @@
366 /*Gather from all cpus*/
367 ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
368 xDelete<IssmDouble>(values);
369-
370+
371 if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
372 xDelete<IssmDouble>(allvalues);
373 }
374@@ -2074,13 +2074,13 @@
375 case MaterialsRheologyBbarEnum: this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
376 case VelEnum: this->ElementResponsex(responses,VelEnum); break;
377 case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
378- default:
379+ default:
380 if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){
381 IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum);
382 *responses=double_result;
383 }
384- else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
385- break;
386+ else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
387+ break;
388 }
389
390 }
391@@ -2211,7 +2211,7 @@
392 weights_input->GetInputValue(&weight,gauss,ThicknessAbsGradientEnum);
393 thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
394
395- /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
396+ /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
397 J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
398 }
399
400@@ -2402,9 +2402,9 @@
401 Analysis* analysis= EnumToAnalysis(analysis_type);
402 analysis->UpdateConstraints(this);
403 delete analysis;
404-
405+
406 /*Second, constraints might be time dependent: */
407- SpcNodesx(nodes,constraints,parameters,analysis_type);
408+ SpcNodesx(nodes,constraints,parameters,analysis_type);
409
410 /*Now, update degrees of freedoms: */
411 NodesDofx(nodes,parameters,analysis_type);
412@@ -2415,7 +2415,7 @@
413
414 IssmDouble *surface = NULL;
415 IssmDouble *bed = NULL;
416-
417+
418 if(VerboseSolution()) _printf0_(" updating vertices positions\n");
419
420 /*get vertex vectors for bed and thickness: */
421@@ -2454,7 +2454,7 @@
422 /*}}}*/
423
424 /*AMR*/
425-#if !defined(_HAVE_ADOLC_)
426+#if !defined(_HAVE_ADOLC_)
427 void FemModel::ReMesh(void){/*{{{*/
428
429 /*Intermediaries*/
430@@ -2464,13 +2464,13 @@
431 int *newelementslist = NULL;
432 int newnumberofvertices = -1;
433 int newnumberofelements = -1;
434- bool* my_elements = NULL;
435+ bool* my_elements = NULL;
436 int* my_vertices = NULL;
437 int elementswidth = this->GetElementsWidth();//just tria elements in this version
438 int amrtype,basalforcing_model;
439 bool isgroundingline;
440
441- /*Branch to specific amr depending on requested method*/
442+ /*Branch to specific amr depending on requested method*/
443 this->parameters->FindParam(&amrtype,AmrTypeEnum);
444 switch(amrtype){
445 #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
446@@ -2483,7 +2483,7 @@
447
448 default: _error_("not implemented yet");
449 }
450-
451+
452 /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
453 this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
454
455@@ -2514,7 +2514,7 @@
456 int analysis_enum = this->analysis_type_list[i];
457
458 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
459- if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
460+ if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
461
462 this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
463 this->CreateConstraints(new_vertices,nodecounter,constraintcounter,analysis_enum,new_constraints);
464@@ -2544,7 +2544,7 @@
465 analysis_type=this->analysis_type_list[i];
466 //SetCurrentConfiguration(analysis_type);
467
468- this->analysis_counter=i;
469+ this->analysis_counter=i;
470 /*Now, plug analysis_counter and analysis_type inside the parameters: */
471 this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
472 this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
473@@ -2561,7 +2561,7 @@
474 }
475
476 ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
477- if(i==0){
478+ if(i==0){
479 VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
480 }
481 SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type);
482@@ -2605,7 +2605,7 @@
483
484 /*Insert bedrock from mismip+ setup*/
485 /*This was used to Misomip project/simulations*/
486-
487+
488 if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n");
489
490 IssmDouble x,y,bx,by;
491@@ -2622,7 +2622,7 @@
492 bx = -150.-728.8*pow(x/300000.,2)+343.91*pow(x/300000.,4)-50.57*pow(x/300000.,6);
493 by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.)));
494 r[i] = max(bx+by,-720.);
495- }
496+ }
497 /*insert new bedrock*/
498 element->AddInput(BedEnum,&r[0],P1Enum);
499 /*Cleanup*/
500@@ -2635,7 +2635,7 @@
501 void FemModel::AdjustBaseThicknessAndMask(void){/*{{{*/
502
503 if(VerboseSolution())_printf0_(" call adjust base and thickness module\n");
504-
505+
506 int numvertices = this->GetElementsWidth();
507 IssmDouble rho_water,rho_ice,density,base_float;
508 IssmDouble* phi = xNew<IssmDouble>(numvertices);
509@@ -2647,7 +2647,7 @@
510
511 for(int el=0;el<this->elements->Size();el++){
512 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
513-
514+
515 element->GetInputListOnVertices(&s[0],SurfaceEnum);
516 element->GetInputListOnVertices(&r[0],BedEnum);
517 element->GetInputListOnVertices(&sl[0],SealevelEnum);
518@@ -2659,16 +2659,16 @@
519 /*calculate base floatation (which supports given surface*/
520 base_float = rho_ice*s[i]/(rho_ice-rho_water);
521 if(r[i]>base_float){
522- b[i] = r[i];
523- }
524+ b[i] = r[i];
525+ }
526 else {
527- b[i] = base_float;
528- }
529+ b[i] = base_float;
530+ }
531
532 if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!");
533 /*update thickness and mask grounded ice level set*/
534 h[i] = s[i]-b[i];
535- phi[i] = h[i]+r[i]/density;
536+ phi[i] = h[i]+r[i]/density;
537 }
538
539 /*Update inputs*/
540@@ -2676,7 +2676,7 @@
541 element->AddInput(ThicknessEnum,&h[0],P1Enum);
542 element->AddInput(BaseEnum,&b[0],P1Enum);
543 }
544-
545+
546 /*Delete*/
547 xDelete<IssmDouble>(phi);
548 xDelete<IssmDouble>(h);
549@@ -2704,7 +2704,7 @@
550 int* P1input_interp = NULL;
551 Vector<IssmDouble>* input_interpolations = NULL;
552 IssmDouble* input_interpolations_serial = NULL;
553- int* pos = NULL;
554+ int* pos = NULL;
555 IssmDouble value = 0;
556
557 /*Figure out how many inputs we have and their respective interpolation*/
558@@ -2724,7 +2724,7 @@
559 P0input_enums = xNew<int>(numP0inputs);
560 P0input_interp = xNew<int>(numP0inputs);
561 P1input_enums = xNew<int>(numP1inputs);
562- P1input_interp = xNew<int>(numP1inputs);
563+ P1input_interp = xNew<int>(numP1inputs);
564 }
565 numP0inputs = 0;
566 numP1inputs = 0;
567@@ -2756,23 +2756,23 @@
568 }
569 }
570 }
571-
572- /*Get P0 and P1 inputs over the elements*/
573+
574+ /*Get P0 and P1 inputs over the elements*/
575 pos = xNew<int>(elementswidth);
576 vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
577 vP1inputs= new Vector<IssmDouble>(numberofvertices*numP1inputs);
578 for(int i=0;i<this->elements->Size();i++){
579 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
580-
581+
582 /*Get P0 inputs*/
583 for(int j=0;j<numP0inputs;j++){
584- TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));
585+ TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));
586 input->GetInputAverage(&value);
587 pos[0]=element->Sid()*numP0inputs+j;
588- /*Insert input in the vector*/
589+ /*Insert input in the vector*/
590 vP0inputs->SetValues(1,pos,&value,INS_VAL);
591 }
592-
593+
594 /*Get P1 inputs*/
595 for(int j=0;j<numP1inputs;j++){
596 TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P1input_enums[j]));
597@@ -2779,8 +2779,8 @@
598 pos[0]=element->vertices[0]->Sid()*numP1inputs+j;
599 pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
600 pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
601- /*Insert input in the vector*/
602- vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);
603+ /*Insert input in the vector*/
604+ vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);
605 }
606 }
607
608@@ -2789,16 +2789,16 @@
609 vP1inputs->Assemble();
610 P0inputs=vP0inputs->ToMPISerial();
611 P1inputs=vP1inputs->ToMPISerial();
612-
613+
614 /*Assign pointers*/
615- *pnumP0inputs = numP0inputs;
616- *pP0inputs = P0inputs;
617+ *pnumP0inputs = numP0inputs;
618+ *pP0inputs = P0inputs;
619 *pP0input_enums = P0input_enums;
620 *pP0input_interp = P0input_interp;
621- *pnumP1inputs = numP1inputs;
622- *pP1inputs = P1inputs;
623+ *pnumP1inputs = numP1inputs;
624+ *pP1inputs = P1inputs;
625 *pP1input_enums = P1input_enums;
626- *pP1input_interp = P1input_interp;
627+ *pP1input_interp = P1input_interp;
628
629 /*Cleanup*/
630 delete input_interpolations;
631@@ -2809,7 +2809,7 @@
632 }
633 /*}}}*/
634 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
635-
636+
637 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh
638 int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition
639 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh
640@@ -2825,7 +2825,7 @@
641 IssmDouble* newP1inputs = NULL; //just on the new partition
642 int* P1input_enums = NULL;
643 int* P1input_interp = NULL;
644- IssmDouble* values = NULL;
645+ IssmDouble* values = NULL;
646 IssmDouble* vector = NULL;
647 IssmDouble* x = NULL;//global, entire old mesh
648 IssmDouble* y = NULL;//global, entire old mesh
649@@ -2837,7 +2837,7 @@
650 IssmDouble* newxc = NULL;//just on the new partition
651 IssmDouble* newyc = NULL;//just on the new partition
652 int* newelementslist = NULL;//just on the new partition
653- int* sidtoindex = NULL;//global vertices sid to partition index
654+ int* sidtoindex = NULL;//global vertices sid to partition index
655
656 /*Get old P0 and P1 inputs (entire mesh)*/
657 this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp);
658@@ -2870,17 +2870,17 @@
659 newx,newy,newnumberofvertices,NULL);
660
661 /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
662- values=xNew<IssmDouble>(elementswidth);
663+ values=xNew<IssmDouble>(elementswidth);
664 for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
665 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
666 /*newP0inputs is just on the new partition*/
667 for(int j=0;j<numP0inputs;j++){
668- switch(P0input_interp[j]){
669+ switch(P0input_interp[j]){
670 case P0Enum:
671 case DoubleInputEnum:
672 element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j]));
673 break;
674- case IntInputEnum:
675+ case IntInputEnum:
676 element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j])));
677 break;
678 case BoolInputEnum:
679@@ -2898,7 +2898,7 @@
680 element->inputs->AddInput(new TriaInput(P1input_enums[j],values,P1Enum));
681 }
682 }
683-
684+
685 /*Cleanup*/
686 xDelete<IssmDouble>(P0inputs);
687 xDelete<IssmDouble>(newP0inputs);
688@@ -2937,7 +2937,7 @@
689 int* elementslist = NULL;
690
691 if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
692-
693+
694 parameters->FindParam(&step,StepEnum);
695 parameters->FindParam(&time,TimeEnum);
696 numberofelements=this->elements->NumberOfElements();
697@@ -2955,7 +2955,7 @@
698
699 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
700 y,numberofvertices,1,step,time));
701-
702+
703 /*Cleanup*/
704 xDelete<IssmDouble>(x);
705 xDelete<IssmDouble>(y);
706@@ -3016,7 +3016,7 @@
707
708 /*Assemble*/
709 vmasklevelset->Assemble();
710-
711+
712 /*Serialize and set output*/
713 (*pmasklevelset)=vmasklevelset->ToMPISerial();
714
715@@ -3030,7 +3030,7 @@
716 void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/
717
718 /*newelementslist is in Matlab indexing*/
719-
720+
721 /*Creating connectivity table*/
722 int* connectivity=NULL;
723 connectivity=xNewZeroInit<int>(newnumberofvertices);
724@@ -3041,12 +3041,12 @@
725 _assert_(vertexid>0 && vertexid-1<newnumberofvertices);//Matlab indexing
726 connectivity[vertexid-1]+=1;//Matlab to C indexing
727 }
728- }
729+ }
730
731 /*Create vertex and insert in vertices*/
732 for(int i=0;i<newnumberofvertices;i++){
733 if(my_vertices[i]){
734- Vertex *newvertex=new Vertex();
735+ Vertex *newvertex=new Vertex();
736 newvertex->id=i+1;
737 newvertex->sid=i;
738 newvertex->pid=UNDEF;
739@@ -3057,8 +3057,8 @@
740 newvertex->sigma=0.;
741 newvertex->connectivity=connectivity[i];
742 newvertex->clone=false;//itapopo check this
743- vertices->AddObject(newvertex);
744- }
745+ vertices->AddObject(newvertex);
746+ }
747 }
748
749 xDelete<int>(connectivity);
750@@ -3085,13 +3085,13 @@
751 for(int j=0;j<nummodels;j++) newtria->element_type_list[j]=0;
752 }
753 else newtria->element_type_list=NULL;
754-
755+
756 /*Element hook*/
757 int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
758 int material_id=i+1; // retrieve material_id = i+1;
759 /*retrieve vertices ids*/
760 int* vertex_ids=xNew<int>(elementswidth);
761- for(int j=0;j<elementswidth;j++) vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
762+ for(int j=0;j<elementswidth;j++) vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
763 /*Setting the hooks*/
764 newtria->numanalyses =this->nummodels;
765 newtria->hnodes =new Hook*[this->nummodels];
766@@ -3103,8 +3103,8 @@
767 for(int j=0;j<this->nummodels;j++) newtria->hnodes[j]=NULL;
768 /*Clean up*/
769 xDelete<int>(vertex_ids);
770- elements->AddObject(newtria);
771- }
772+ elements->AddObject(newtria);
773+ }
774 }
775
776 }
777@@ -3114,14 +3114,14 @@
778 /*Just Matice in this version*/
779 for(int i=0;i<newnumberofelements;i++){
780 if(my_elements[i]){
781- materials->AddObject(new Matice(i+1,i,MaticeEnum));
782- }
783+ materials->AddObject(new Matice(i+1,i,MaticeEnum));
784+ }
785 }
786-
787+
788 /*Add new constant material property to materials, at the end: */
789 Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
790 newmatpar->SetMid(newnumberofelements+1);
791- materials->AddObject(newmatpar);//put it at the end of the materials
792+ materials->AddObject(newmatpar);//put it at the end of the materials
793 }
794 /*}}}*/
795 void FemModel::CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes){/*{{{*/
796@@ -3128,23 +3128,23 @@
797
798 int lid=0;
799 for(int j=0;j<newnumberofvertices;j++){
800- if(my_vertices[j]){
801-
802- Node* newnode=new Node();
803-
804+ if(my_vertices[j]){
805+
806+ Node* newnode=new Node();
807+
808 /*id: */
809 newnode->id=nodecounter+j+1;
810 newnode->sid=j;
811 newnode->lid=lid++;
812 newnode->analysis_enum=analysis_enum;
813-
814+
815 /*Initialize coord_system: Identity matrix by default*/
816 for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
817 for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
818-
819+
820 /*indexing:*/
821 newnode->indexingupdate=true;
822-
823+
824 Analysis* analysis=EnumToAnalysis(analysis_enum);
825 int *doftypes=NULL;
826 int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum);
827@@ -3158,7 +3158,7 @@
828
829 /*Stressbalance Horiz*/
830 if(analysis_enum==StressbalanceAnalysisEnum){
831- // itapopo this code is rarely used.
832+ // itapopo this code is rarely used.
833 /*Coordinate system provided, convert to coord_system matrix*/
834 //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
835 //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);
836@@ -3173,13 +3173,13 @@
837
838 if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
839 if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
840-
841+
842 int numberofvertices = femmodel_vertices->NumberOfVertices();
843 int numberofelements = femmodel_elements->NumberOfElements();
844 int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements)
845 IssmDouble* x = NULL;
846 IssmDouble* y = NULL;
847- IssmDouble* z = NULL;
848+ IssmDouble* z = NULL;
849 int* elementslist = NULL;
850 int* elem_vertices = NULL;
851 IssmDouble *id1 = NULL;
852@@ -3188,7 +3188,7 @@
853
854 /*Get vertices coordinates*/
855 VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
856-
857+
858 /*Get element vertices*/
859 elem_vertices = xNew<int>(elementswidth);
860 Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
861@@ -3203,7 +3203,7 @@
862 vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
863 vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
864 }
865-
866+
867 /*Assemble*/
868 vid1->Assemble();
869 vid2->Assemble();
870@@ -3213,7 +3213,7 @@
871 id1 = vid1->ToMPISerial();
872 id2 = vid2->ToMPISerial();
873 id3 = vid3->ToMPISerial();
874-
875+
876 /*Construct elements list*/
877 elementslist=xNew<int>(numberofelements*elementswidth);
878 if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
879@@ -3222,7 +3222,7 @@
880 elementslist[elementswidth*i+1] = reCast<int>(id2[i])+1; //InterpMesh wants Matlab indexing
881 elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
882 }
883-
884+
885 /*Assign pointers*/
886 *px = x;
887 *py = y;
888@@ -3243,7 +3243,7 @@
889
890 if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
891 if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
892-
893+
894 int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition
895 int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh)
896 int numberofelements = femmodel_elements->Size(); //number of elements of this partition
897@@ -3250,20 +3250,20 @@
898 int elementswidth = this->GetElementsWidth(); //just 2D mesh in this version (just tria elements)
899 IssmDouble* x = NULL;
900 IssmDouble* y = NULL;
901- IssmDouble* z = NULL;
902+ IssmDouble* z = NULL;
903 int* elementslist = NULL;
904 int* sidtoindex = NULL;
905 int* elem_vertices = NULL;
906-
907+
908 /*Get vertices coordinates of this partition*/
909 sidtoindex = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
910 x = xNew<IssmDouble>(numberofvertices);//just this partition
911 y = xNew<IssmDouble>(numberofvertices);//just this partitio;
912 z = xNew<IssmDouble>(numberofvertices);//just this partitio;
913-
914+
915 /*Go through in this partition (vertices)*/
916 for(int i=0;i<numberofvertices;i++){//just this partition
917- Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
918+ Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
919 /*Attention: no spherical coordinates*/
920 x[i]=vertex->GetX();
921 y[i]=vertex->GetY();
922@@ -3276,7 +3276,7 @@
923 elem_vertices= xNew<int>(elementswidth);
924 elementslist = xNew<int>(numberofelements*elementswidth);
925 if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
926-
927+
928 for(int i=0;i<numberofelements;i++){//just this partition
929 Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
930 element->GetVerticesSidList(elem_vertices);
931@@ -3283,14 +3283,14 @@
932 elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing
933 elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
934 elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
935- }
936-
937+ }
938+
939 /*Assign pointers*/
940 *px = x;
941 *py = y;
942 *pz = z;
943 *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
944- *psidtoindex = sidtoindex; //it is ncessary to insert inputs
945+ *psidtoindex = sidtoindex; //it is ncessary to insert inputs
946
947 /*Cleanup*/
948 xDelete<int>(elem_vertices);
949@@ -3301,9 +3301,9 @@
950 /*ATTENTION: JUST SPCVX AND SPCVY*/
951 /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
952 if(analysis_enum!=StressbalanceAnalysisEnum) return;
953-
954+
955 int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum);
956- int dofpernode = 2; //vx and vy
957+ int dofpernode = 2; //vx and vy
958 int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector
959 int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh
960 int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh
961@@ -3327,7 +3327,7 @@
962 newx=xNew<IssmDouble>(newnumberofvertices);//just the new partition
963 newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
964 for(int i=0;i<newnumberofvertices;i++){//just the new partition
965- Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
966+ Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
967 /*Attention: no spherical coordinates*/
968 newx[i]=vertex->GetX();
969 newy[i]=vertex->GetY();
970@@ -3335,7 +3335,7 @@
971
972 /*Get spcvx and spcvy of old mesh*/
973 for(int i=0;i<this->constraints->Size();i++){
974-
975+
976 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
977 if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n");
978
979@@ -3342,14 +3342,14 @@
980 SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint);
981 int dof = spcstatic->GetDof();
982 int node = spcstatic->GetNodeId();
983- IssmDouble spcvalue = spcstatic->GetValue();
984+ IssmDouble spcvalue = spcstatic->GetValue();
985 int nodeindex = node-1;
986-
987+
988 /*vx and vx flag insertion*/
989 if(dof==0) {//vx
990 vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx
991 vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
992- }
993+ }
994 /*vy and vy flag insertion*/
995 if(dof==1){//vy
996 vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL); //vy
997@@ -3365,7 +3365,7 @@
998 InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements,
999 spc,numberofvertices,numberofcols,
1000 newx,newy,newnumberofvertices,NULL);
1001-
1002+
1003 /*Now, insert the interpolated constraints in the data set (constraints)*/
1004 count=0;
1005 for(int i=0;i<newnumberofvertices;i++){//just in the new partition
1006@@ -3382,7 +3382,7 @@
1007 Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
1008 /*spcvy*/
1009 if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
1010- newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
1011+ newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
1012 //add count'th spc, on node i+1, setting dof 1 to vx.
1013 count++;
1014 }
1015@@ -3439,13 +3439,13 @@
1016 int numprocs = IssmComm::GetSize();
1017 bool *my_elements = NULL;
1018 int *my_vertices = NULL;
1019-
1020- _assert_(newnumberofvertices>0);
1021- _assert_(newnumberofelements>0);
1022+
1023+ _assert_(newnumberofvertices>0);
1024+ _assert_(newnumberofelements>0);
1025 epart=xNew<int>(newnumberofelements);
1026 npart=xNew<int>(newnumberofvertices);
1027 index=xNew<int>(elementswidth*newnumberofelements);
1028-
1029+
1030 for (int i=0;i<newnumberofelements;i++){
1031 for (int j=0;j<elementswidth;j++){
1032 *(index+elementswidth*i+j)=(*(newelementslist+elementswidth*i+j))-1; //-1 for C indexing in Metis
1033@@ -3465,7 +3465,7 @@
1034 for (int i=0;i<newnumberofelements;i++) epart[i]=0;
1035 for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
1036 }
1037- else _error_("At least one processor is required");
1038+ else _error_("At least one processor is required");
1039
1040 my_vertices=xNew<int>(newnumberofvertices);
1041 my_elements=xNew<bool>(newnumberofelements);
1042@@ -3475,14 +3475,14 @@
1043 /*Start figuring out, out of the partition, which elements belong to this cpu: */
1044 for(int i=0;i<newnumberofelements;i++){
1045 /*!All elements have been partitioned above, only deal with elements for this cpu: */
1046- if(my_rank==epart[i]){
1047+ if(my_rank==epart[i]){
1048 my_elements[i]=true;
1049- /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
1050- *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
1051- into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
1052+ /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
1053+ *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
1054+ into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
1055 will hold which vertices belong to this partition*/
1056 for(int j=0;j<elementswidth;j++){
1057- _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing
1058+ _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing
1059 my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing
1060 }
1061 }
1062@@ -3494,18 +3494,18 @@
1063
1064 /*Free ressources:*/
1065 xDelete<int>(epart);
1066- xDelete<int>(npart);
1067+ xDelete<int>(npart);
1068 xDelete<int>(index);
1069 }
1070 /*}}}*/
1071 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
1072-
1073+
1074 int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements
1075 int numberofvertices = this->vertices->NumberOfVertices();
1076 IssmDouble weight = 0.;
1077- IssmDouble* tauxx = NULL;
1078- IssmDouble* tauyy = NULL;
1079- IssmDouble* tauxy = NULL;
1080+ IssmDouble* tauxx = NULL;
1081+ IssmDouble* tauyy = NULL;
1082+ IssmDouble* tauxy = NULL;
1083 IssmDouble* totalweight = NULL;
1084 IssmDouble* deviatoricstressxx = xNew<IssmDouble>(elementswidth);
1085 IssmDouble* deviatoricstressyy = xNew<IssmDouble>(elementswidth);
1086@@ -3515,10 +3515,10 @@
1087 Vector<IssmDouble>* vectauyy = new Vector<IssmDouble>(numberofvertices);
1088 Vector<IssmDouble>* vectauxy = new Vector<IssmDouble>(numberofvertices);
1089 Vector<IssmDouble>* vectotalweight = new Vector<IssmDouble>(numberofvertices);
1090-
1091+
1092 /*Update the Deviatoric Stress tensor over the elements*/
1093 this->DeviatoricStressx();
1094-
1095+
1096 /*Calculate the Smoothed Deviatoric Stress tensor*/
1097 for(int i=0;i<this->elements->Size();i++){
1098 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1099@@ -3563,7 +3563,7 @@
1100
1101 /*Divide for the total weight*/
1102 for(int i=0;i<numberofvertices;i++){
1103- _assert_(totalweight[i]>0);
1104+ _assert_(totalweight[i]>0);
1105 tauxx[i] = tauxx[i]/totalweight[i];
1106 tauyy[i] = tauyy[i]/totalweight[i];
1107 tauxy[i] = tauxy[i]/totalweight[i];
1108@@ -3588,7 +3588,7 @@
1109 /*}}}*/
1110 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
1111
1112- /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
1113+ /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
1114 * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
1115
1116 IssmDouble Jdet,error,ftxx,ftyy,ftxy;
1117@@ -3608,7 +3608,7 @@
1118
1119 /*Get smoothed deviatoric stress tensor*/
1120 this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
1121-
1122+
1123 /*Integrate the error over elements*/
1124 for(int i=0;i<this->elements->Size();i++){
1125 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1126@@ -3616,7 +3616,7 @@
1127 element->GetInputListOnVertices(tauyy,DeviatoricStressyyEnum);
1128 element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
1129 element->GetVerticesSidList(elem_vertices);
1130-
1131+
1132 /*Integrate*/
1133 element->GetVerticesCoordinates(&xyz_list);
1134 Gauss* gauss=element->NewGauss(2);
1135@@ -3631,9 +3631,9 @@
1136 ftyy+=(tauyy[n]-smoothedtauyy[elem_vertices[n]])*basis[n];
1137 ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
1138 }
1139- error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
1140+ error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
1141 }
1142- /*Set the error in the global vector*/
1143+ /*Set the error in the global vector*/
1144 sid=element->Sid();
1145 error = sqrt(error);//sqrt(e^2)
1146 velementerror->SetValue(sid,error,INS_VAL);
1147@@ -3647,7 +3647,7 @@
1148
1149 /*Serialize and set output*/
1150 (*pelementerror)=velementerror->ToMPISerial();
1151-
1152+
1153 /*Cleanup*/
1154 xDelete<IssmDouble>(smoothedtauxx);
1155 xDelete<IssmDouble>(smoothedtauyy);
1156@@ -3691,7 +3691,7 @@
1157 /*weight to calculate the smoothed grad H*/
1158 Tria* triaelement = xDynamicCast<Tria*>(element);
1159 weight = triaelement->GetArea();//the tria area is a choice for the weight
1160-
1161+
1162 /*dH/dx*/
1163 vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
1164 vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL);
1165@@ -3759,7 +3759,7 @@
1166
1167 /*Get smoothed deviatoric stress tensor*/
1168 this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
1169-
1170+
1171 /*Integrate the error over elements*/
1172 for(int i=0;i<this->elements->Size();i++){
1173 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1174@@ -3845,7 +3845,7 @@
1175 Vector<IssmDouble>* vyc = new Vector<IssmDouble>(numberofelements);
1176 IssmDouble* x = NULL;
1177 IssmDouble* y = NULL;
1178- IssmDouble* z = NULL;
1179+ IssmDouble* z = NULL;
1180 IssmDouble* xyz_list = NULL;
1181 IssmDouble x1,y1,x2,y2,x3,y3;
1182
1183@@ -3854,7 +3854,7 @@
1184 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1185 //element->GetVerticesSidList(elem_vertices);
1186 int sid = element->Sid();
1187- element->GetVerticesCoordinates(&xyz_list);
1188+ element->GetVerticesCoordinates(&xyz_list);
1189 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
1190 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
1191 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
1192@@ -3887,11 +3887,11 @@
1193 if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){
1194 _error_("level set type not implemented yet!");
1195 }
1196-
1197+
1198 /*Outputs*/
1199 IssmDouble* zerolevelset_points = NULL;
1200 int npoints = 0;
1201-
1202+
1203 /*Intermediaries*/
1204 int elementswidth = this->GetElementsWidth();
1205 int numberofelements = this->elements->NumberOfElements();
1206@@ -3904,7 +3904,7 @@
1207 IssmDouble* y_zerolevelset = NULL;
1208 int count,sid;
1209 IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
1210-
1211+
1212 /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
1213 for(int i=0;i<this->elements->Size();i++){
1214 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1215@@ -3911,15 +3911,15 @@
1216 element->GetInputListOnVertices(levelset,levelset_type);
1217 element->GetVerticesSidList(elem_vertices);
1218 sid= element->Sid();
1219- element->GetVerticesCoordinates(&xyz_list);
1220+ element->GetVerticesCoordinates(&xyz_list);
1221 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
1222 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
1223 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
1224 xc = NAN;
1225- yc = NAN;
1226+ yc = NAN;
1227 Tria* tria = xDynamicCast<Tria*>(element);
1228 if(tria->IsIceInElement()){/*verify if there is ice in the element*/
1229- if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
1230+ if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
1231 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
1232 xc=(x1+x2+x3)/3.;
1233 yc=(y1+y2+y3)/3.;
1234@@ -3949,7 +3949,7 @@
1235 count++;
1236 }
1237 }
1238-
1239+
1240 /*Assign outputs*/
1241 numberofpoints = npoints;
1242 (*pzerolevelset_points) = zerolevelset_points;
1243@@ -3989,7 +3989,7 @@
1244 /*save the d_responses pointer: */
1245 responses_pointer=d_responses;
1246
1247- //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled.
1248+ //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled.
1249 //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
1250
1251 for(i=0;i<numresponsedescriptors;i++){
1252@@ -4082,10 +4082,10 @@
1253 void FemModel::EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy){/*{{{*/
1254
1255 int ns,nsmax;
1256-
1257+
1258 /*Go through elements, and add contribution from each element to the deflection vector wg:*/
1259 ns = elements->Size();
1260-
1261+
1262 /*Figure out max of ns: */
1263 ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
1264 ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
1265@@ -4104,7 +4104,7 @@
1266 pY->Assemble();
1267 }
1268 }
1269-
1270+
1271 /*One last time: */
1272 pUp->Assemble();
1273 pNorth->Assemble();
1274@@ -4123,11 +4123,11 @@
1275 IssmDouble eartharea_cpu=0;
1276
1277 int ns,nsmax;
1278-
1279+
1280 /*Go through elements, and add contribution from each element to the deflection vector wg:*/
1281 ns = elements->Size();
1282-
1283- /*First, figure out the surface area of Earth: */
1284+
1285+ /*First, figure out the surface area of Earth: */
1286 for(int i=0;i<ns;i++){
1287 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
1288 eartharea_cpu += element->GetAreaSpherical();
1289@@ -4151,7 +4151,7 @@
1290 pEast->Assemble();
1291 }
1292 }
1293-
1294+
1295 /*One last time: */
1296 pUp->Assemble();
1297 pNorth->Assemble();
1298@@ -4213,7 +4213,7 @@
1299 if(i%loop==0)pRSLgi->Assemble();
1300 }
1301 if(VerboseConvergence())_printf0_("\n");
1302-
1303+
1304 /*One last time: */
1305 pRSLgi->Assemble();
1306
1307@@ -4231,12 +4231,12 @@
1308
1309 /*serialized vectors:*/
1310 IssmDouble* RSLg_old=NULL;
1311-
1312+
1313 IssmDouble eartharea=0;
1314 IssmDouble eartharea_cpu=0;
1315
1316 int ns,nsmax;
1317-
1318+
1319 /*Serialize vectors from previous iteration:*/
1320 RSLg_old=pRSLg_old->ToMPISerial();
1321
1322@@ -4248,7 +4248,7 @@
1323 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
1324 eartharea_cpu += element->GetAreaSpherical();
1325 }
1326-
1327+
1328 ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
1329 ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1330
1331@@ -4266,7 +4266,7 @@
1332 if(i%loop==0)pRSLgo->Assemble();
1333 }
1334 if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
1335-
1336+
1337 /*Free ressources:*/
1338 xDelete<IssmDouble>(RSLg_old);
1339 }
1340@@ -4354,8 +4354,8 @@
1341
1342 /*Free ressources:*/
1343 xDelete<IssmDouble>(RSLg_old);
1344-
1345-
1346+
1347+
1348 }
1349 /*}}}*/
1350 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
1351@@ -4362,18 +4362,18 @@
1352
1353 /*serialized vectors:*/
1354 IssmDouble* RSLg=NULL;
1355-
1356+
1357 IssmDouble eartharea=0;
1358 IssmDouble eartharea_cpu=0;
1359
1360 int ns,nsmax;
1361-
1362+
1363 /*Serialize vectors from previous iteration:*/
1364 RSLg=pRSLg->ToMPISerial();
1365
1366 /*Go through elements, and add contribution from each element to the deflection vector wg:*/
1367 ns = elements->Size();
1368-
1369+
1370 /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
1371 for(int i=0;i<ns;i++){
1372 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
1373@@ -4401,7 +4401,7 @@
1374 }
1375 }
1376 }
1377-
1378+
1379 /*One last time: */
1380 pUp->Assemble();
1381 if(horiz){
1382@@ -4435,13 +4435,13 @@
1383 }
1384 ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
1385 ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1386-
1387+
1388 ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
1389 ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1390
1391 /*Free ressources:*/
1392 xDelete<IssmDouble>(RSLg_serial);
1393-
1394+
1395 return oceanvalue/oceanarea;
1396 }
1397 /*}}}*/
1398@@ -4457,7 +4457,7 @@
1399 IssmDouble* old_active = NULL;
1400 int* eplzigzag_counter = NULL;
1401 int eplflip_lock;
1402-
1403+
1404 HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis();
1405 HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis();
1406
1407@@ -4464,10 +4464,10 @@
1408 /*Step 1: update mask, the mask might be extended by residual and/or using downstream sediment head*/
1409 mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
1410 recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
1411- this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
1412- this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
1413+ this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
1414+ this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
1415 GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
1416-
1417+
1418 for (int i=0;i<elements->Size();i++){
1419 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
1420 effanalysis->HydrologyEPLGetMask(mask,recurence,eplzigzag_counter,element);
1421@@ -4485,8 +4485,8 @@
1422 this->parameters->SetParam(eplzigzag_counter,this->nodes->Size(),EplZigZagCounterEnum);
1423 /*Assemble and serialize*/
1424 mask->Assemble();
1425- serial_mask=mask->ToMPISerial();
1426-
1427+ serial_mask=mask->ToMPISerial();
1428+
1429 xDelete<int>(eplzigzag_counter);
1430 xDelete<IssmDouble>(serial_rec);
1431 xDelete<IssmDouble>(old_active);
1432@@ -4528,7 +4528,7 @@
1433 delete inefanalysis;
1434 int sum_counter;
1435 ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
1436- ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
1437+ ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
1438 counter=sum_counter;
1439 *pEplcount = counter;
1440 if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n");
1441@@ -4649,7 +4649,7 @@
1442 xDelete<IssmDouble>(serial_active);
1443 int sum_counter;
1444 ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
1445- ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
1446+ ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
1447 counter=sum_counter;
1448 *pL2count = counter;
1449 if(VerboseSolution()) _printf0_(" Number of active nodes L2 Projection: "<< counter <<"\n");
1450@@ -4734,7 +4734,7 @@
1451 }
1452 }
1453 /*}}}*/
1454-#ifdef _HAVE_JAVASCRIPT_
1455+#ifdef _HAVE_JAVASCRIPT_
1456 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
1457 /*configuration: */
1458 int solution_type;
1459@@ -4749,12 +4749,12 @@
1460
1461 /*From command line arguments, retrieve different filenames needed to create the FemModel: */
1462 solution_type=StringToEnumx(solution);
1463-
1464+
1465 /*Create femmodel from input files: */
1466 profiler->Start(MPROCESSOR);
1467 this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
1468 profiler->Stop(MPROCESSOR);
1469-
1470+
1471 /*Save communicator in the parameters dataset: */
1472 this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
1473
1474@@ -4769,7 +4769,7 @@
1475 char** poutputbuffer;
1476 size_t* poutputbuffersize;
1477
1478-
1479+
1480 /*Before we delete the profiler, report statistics for this run: */
1481 profiler->Stop(TOTAL); //final tagging
1482 _printf0_("\n");
1483@@ -4782,7 +4782,7 @@
1484 <<profiler->TotalTimeModSec(TOTAL)<<" sec"
1485 );
1486 _printf0_("\n");
1487-
1488+
1489 /*Before we close the output file, recover the buffer and size:*/
1490 outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
1491 poutputbuffer=outputbufferparam->GetParameterValue();
1492@@ -4822,8 +4822,8 @@
1493 fclose(toolkitsoptionsfid);
1494
1495 /*Open output file once for all and add output file descriptor to parameters*/
1496- output_fid=open_memstream(&outputbuffer,&outputsize);
1497- if(output_fid==NULL||output_fid==0)_error_("could not initialize output stream");
1498+ output_fid=open_memstream(&outputbuffer,&outputsize);
1499+ if(output_fid==NULL)_error_("could not initialize output stream");
1500 this->parameters->SetParam(output_fid,OutputFilePointerEnum);
1501 this->parameters->AddObject(new GenericParam<char**>(&outputbuffer,OutputBufferPointerEnum));
1502 this->parameters->AddObject(new GenericParam<size_t*>(&outputsize,OutputBufferSizePointerEnum));
1503@@ -4864,7 +4864,7 @@
1504 this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){
1505 /*Initialize hmaxvertices with NAN*/
1506 hmaxvertices_serial=xNew<IssmDouble>(numberofvertices);
1507- for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
1508+ for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
1509 /*Fill hmaxvertices*/
1510 if(this->amrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum);
1511 if(this->amrbamg->icefront_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum);
1512@@ -4871,7 +4871,7 @@
1513 if(this->amrbamg->thicknesserror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum);
1514 if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
1515 }
1516-
1517+
1518 if(my_rank==0){
1519 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
1520 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process.");
1521@@ -4891,10 +4891,10 @@
1522 newz=xNew<IssmDouble>(newnumberofvertices);
1523 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
1524 }
1525- ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1526- ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1527- ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1528- ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
1529+ ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1530+ ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1531+ ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1532+ ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
1533
1534 /*Assign output pointers*/
1535 *pnewnumberofvertices = newnumberofvertices;
1536@@ -4927,7 +4927,7 @@
1537
1538 /*Create bamg data structures for bamg*/
1539 this->amrbamg = new AmrBamg();
1540-
1541+
1542 /*Get amr parameters*/
1543 this->parameters->FindParam(&hmin,AmrHminEnum);
1544 this->parameters->FindParam(&hmax,AmrHmaxEnum);
1545@@ -4951,7 +4951,7 @@
1546 this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation);
1547
1548 /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
1549- if(my_rank==0){
1550+ if(my_rank==0){
1551 this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
1552 }
1553
1554@@ -4965,7 +4965,7 @@
1555 void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/
1556
1557 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
1558-
1559+
1560 /*Intermediaries*/
1561 int numberofvertices = this->vertices->NumberOfVertices();
1562 IssmDouble* verticedistance = NULL;
1563@@ -4972,7 +4972,7 @@
1564 IssmDouble threshold,resolution;
1565
1566 switch(levelset_type){
1567- case MaskGroundediceLevelsetEnum:
1568+ case MaskGroundediceLevelsetEnum:
1569 threshold = this->amrbamg->groundingline_distance;
1570 resolution = this->amrbamg->groundingline_resolution;
1571 break;
1572@@ -4986,7 +4986,7 @@
1573 /*Get vertice distance to zero level set points*/
1574 this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
1575 if(!verticedistance) _error_("verticedistance is NULL!\n");
1576-
1577+
1578 /*Fill hmaxVertices*/
1579 for(int i=0;i<numberofvertices;i++){
1580 if(verticedistance[i]<threshold){
1581@@ -5000,7 +5000,7 @@
1582 }
1583 /*}}}*/
1584 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
1585-
1586+
1587 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
1588
1589 /*Intermediaries*/
1590@@ -5008,7 +5008,7 @@
1591 int numberofelements = this->elements->NumberOfElements();
1592 int numberofvertices = this->vertices->NumberOfVertices();
1593 IssmDouble* maxlength = xNew<IssmDouble>(numberofelements);
1594- IssmDouble* error_vertices = xNewZeroInit<IssmDouble>(numberofvertices);
1595+ IssmDouble* error_vertices = xNewZeroInit<IssmDouble>(numberofvertices);
1596 IssmDouble* error_elements = NULL;
1597 IssmDouble* x = NULL;
1598 IssmDouble* y = NULL;
1599@@ -5021,7 +5021,7 @@
1600
1601 /*Fill variables*/
1602 switch(errorestimator_type){
1603- case ThicknessErrorEstimatorEnum:
1604+ case ThicknessErrorEstimatorEnum:
1605 threshold = this->amrbamg->thicknesserror_threshold;
1606 groupthreshold = this->amrbamg->thicknesserror_groupthreshold;
1607 resolution = this->amrbamg->thicknesserror_resolution;
1608@@ -5046,12 +5046,12 @@
1609 switch(errorestimator_type){
1610 case ThicknessErrorEstimatorEnum: this->amrbamg->thicknesserror_maximum = maxerror;break;
1611 case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break;
1612- }
1613+ }
1614 }
1615
1616 /*Get mesh*/
1617 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
1618-
1619+
1620 /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
1621 for(int i=0;i<numberofelements;i++){
1622 v1=index[i*elementswidth+0]-1;//Matlab to C indexing
1623@@ -5065,7 +5065,7 @@
1624 error_vertices[v1]+=error_elements[i];
1625 error_vertices[v2]+=error_elements[i];
1626 error_vertices[v3]+=error_elements[i];
1627- }
1628+ }
1629
1630 /*Fill hmaxvertices with the criteria*/
1631 for(int i=0;i<numberofelements;i++){
1632@@ -5079,7 +5079,7 @@
1633 if(error_vertices[vid]>groupthreshold*maxerror) refine=true;
1634 }
1635 }
1636- /*Now, fill the hmaxvertices if requested*/
1637+ /*Now, fill the hmaxvertices if requested*/
1638 if(refine){
1639 for(int j=0;j<elementswidth;j++){
1640 vid=index[i*elementswidth+j]-1;//Matlab to C indexing
1641@@ -5109,7 +5109,7 @@
1642
1643 /*Output*/
1644 IssmDouble* verticedistance;
1645-
1646+
1647 /*Intermediaries*/
1648 int numberofvertices = this->vertices->NumberOfVertices();
1649 IssmDouble* levelset_points= NULL;
1650@@ -5121,8 +5121,8 @@
1651
1652 /*Get vertices coordinates*/
1653 VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
1654-
1655- /*Get points which level set is zero (center of elements with zero level set)*/
1656+
1657+ /*Get points which level set is zero (center of elements with zero level set)*/
1658 this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
1659
1660 /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/
1661@@ -5131,9 +5131,9 @@
1662 verticedistance[i]=INFINITY;
1663 for(int j=0;j<numberofpoints;j++){
1664 distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
1665- verticedistance[i]=min(distance,verticedistance[i]);
1666+ verticedistance[i]=min(distance,verticedistance[i]);
1667 }
1668- }
1669+ }
1670
1671 /*Assign the pointer*/
1672 (*pverticedistance)=verticedistance;
1673@@ -5167,12 +5167,12 @@
1674 /*Get fields, if requested*/
1675 if(this->amr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum);
1676 if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
1677- if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror);
1678- if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror);
1679-
1680+ if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror);
1681+ if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror);
1682+
1683 if(my_rank==0){
1684 this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
1685- &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
1686+ &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
1687 newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
1688 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
1689 }
1690@@ -5186,12 +5186,12 @@
1691 newz=xNew<IssmDouble>(newnumberofvertices);
1692 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
1693 }
1694- ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1695- ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1696- ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1697- ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
1698+ ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1699+ ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1700+ ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
1701+ ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
1702
1703- /*Assign the pointers*/
1704+ /*Assign the pointers*/
1705 (*pnewelementslist) = newelementslist; //Matlab indexing
1706 (*pnewx) = newx;
1707 (*pnewy) = newy;
1708@@ -5207,7 +5207,7 @@
1709 }
1710 /*}}}*/
1711 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
1712-
1713+
1714 /*Define variables*/
1715 int my_rank = IssmComm::GetRank();
1716 int numberofvertices = this->vertices->NumberOfVertices();
1717@@ -5216,7 +5216,7 @@
1718 IssmDouble* y = NULL;
1719 IssmDouble* z = NULL;
1720 int* elements = NULL;
1721-
1722+
1723 /*Initialize field as NULL for now*/
1724 this->amr = NULL;
1725
1726@@ -5224,12 +5224,12 @@
1727 /*elements comes in Matlab indexing*/
1728 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
1729
1730- /*Create initial mesh (coarse mesh) in neopz data structure*/
1731+ /*Create initial mesh (coarse mesh) in neopz data structure*/
1732 /*Just CPU #0 should keep AMR object*/
1733 /*Initialize refinement pattern*/
1734 this->SetRefPatterns();
1735 this->amr = new AdaptiveMeshRefinement();
1736- this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
1737+ this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
1738 /*Get amr parameters*/
1739 this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum);
1740 this->parameters->FindParam(&this->amr->gradation,AmrGradationEnum);
1741@@ -5242,7 +5242,7 @@
1742 this->parameters->FindParam(&this->amr->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum);
1743 this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
1744 this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
1745- if(my_rank==0){
1746+ if(my_rank==0){
1747 this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
1748 }
1749
1750@@ -5263,7 +5263,7 @@
1751
1752 /*Output*/
1753 IssmDouble* elementdistance;
1754-
1755+
1756 /*Intermediaries*/
1757 int numberofelements = this->elements->NumberOfElements();
1758 Vector<IssmDouble>* velementdistance = new Vector<IssmDouble>(numberofelements);
1759@@ -5272,14 +5272,14 @@
1760 IssmDouble mindistance,distance;
1761 IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
1762 int numberofpoints;
1763-
1764- /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
1765+
1766+ /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
1767 this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
1768
1769 for(int i=0;i<this->elements->Size();i++){//parallel
1770 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
1771 int sid = element->Sid();
1772- element->GetVerticesCoordinates(&xyz_list);
1773+ element->GetVerticesCoordinates(&xyz_list);
1774 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
1775 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
1776 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
1777@@ -5289,11 +5289,11 @@
1778 /*Loop over each point (where level set is zero)*/
1779 for(int j=0;j<numberofpoints;j++){
1780 distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
1781- mindistance=min(distance,mindistance);
1782+ mindistance=min(distance,mindistance);
1783 }
1784 velementdistance->SetValue(sid,mindistance,INS_VAL);
1785 xDelete<IssmDouble>(xyz_list);
1786- }
1787+ }
1788
1789 /*Assemble*/
1790 velementdistance->Assemble();
Note: See TracBrowser for help on using the repository browser.