Ice Sheet System Model  4.18
Code documentation
Functions
ModelProcessorx.h File Reference
#include "../../classes/classes.h"
#include "../../analyses/analyses.h"

Go to the source code of this file.

Functions

void ModelProcessorx (Elements **pelements, Nodes ***pnodes, Vertices **pvertices, Materials **pmaterials, Constraints ***pconstraints, Loads ***ploads, Parameters **pparameters, Inputs2 **pinputs2, IoModel *iomodel, FILE *toolkitfile, char *rootpath, const int solution_type, const int nummodels, const int *analysis_type_listh)
 
void CreateElements (Elements *elements, IoModel *iomodel, int nummodels)
 
void CreateMaterials (Elements *elements, Inputs2 *inputs2, Materials *materials, IoModel *iomodel, int nummodels)
 
void CreateVertices (Elements *elements, Vertices *vertices, IoModel *iomodel, int solution_type, bool isamr=false)
 
void CreateParameters (Parameters *parameters, IoModel *iomodel, char *rootpath, FILE *toolkitfile, const int solution_type)
 
void CreateParametersAutodiff (Parameters *parameters, IoModel *iomodel)
 
void CreateParametersControl (Parameters *parameters, IoModel *iomodel, int solution_type)
 
void CreateParametersDakota (Parameters *parameters, IoModel *iomodel, char *rootpath)
 
void CreateOutputDefinitions (Elements *elements, Parameters *parameters, Inputs2 *inputs2, IoModel *iomodel)
 
void UpdateElementsAndMaterialsControl (Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
 
void UpdateElementsAndMaterialsControlAD (Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
 
void UpdateElementsAndMaterialsDakota (Elements *elements, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
 
void UpdateElementsTransient (Elements *elements, Parameters *parameters, Inputs2 *inputs2, IoModel *iomodel)
 
void CreateNodes (Nodes *nodes, IoModel *iomodel, int analysis, int finite_element, bool isamr=false, int approximation=NoneApproximationEnum, int *approximations=NULL)
 
void ElementsAndVerticesPartitioning (IoModel *iomodel)
 
void DiscontinuousGalerkinNodesPartitioning (bool **pmy_nodes, bool *my_elements, bool *my_vertices, IoModel *iomodel)
 
void FacesPartitioning (IoModel *iomodel)
 
void EdgesPartitioning (IoModel *iomodel)
 
void CreateEdges (IoModel *iomodel)
 
void CreateFaces (IoModel *iomodel)
 
void CreateFaces3d (IoModel *iomodel)
 
void EdgeOnBoundaryFlags (bool **pflags, IoModel *iomodel)
 
void CreateSingleNodeToElementConnectivity (IoModel *iomodel)
 
void CreateNumberNodeToElementConnectivity (IoModel *iomodel)
 

Function Documentation

◆ ModelProcessorx()

void ModelProcessorx ( Elements **  pelements,
Nodes ***  pnodes,
Vertices **  pvertices,
Materials **  pmaterials,
Constraints ***  pconstraints,
Loads ***  ploads,
Parameters **  pparameters,
Inputs2 **  pinputs2,
IoModel iomodel,
FILE *  toolkitfile,
char *  rootpath,
const int  solution_type,
const int  nummodels,
const int *  analysis_type_listh 
)

Definition at line 15 of file ModelProcessorx.cpp.

15  {
16  _assert_(nummodels>0);
17 
18  /*Set Verbosity once for all*/
19  int verbose;
20  iomodel->FindConstant(&verbose,"md.verbose");
21  SetVerbosityLevel(verbose);
22 
23  if(VerboseMProcessor()) _printf0_(" starting model processor \n");
24 
25  /*Initialize datasets*/
26  Elements *elements = new Elements();
27  Vertices *vertices = new Vertices();
28  Materials *materials = new Materials();
29  Parameters *parameters = new Parameters();
30  Constraints **constraints = xNew<Constraints*>(nummodels);
31  Loads **loads = xNew<Loads*>(nummodels);
32  Nodes **nodes = xNew<Nodes*>(nummodels);
33  for(int i = 0;i<nummodels;i++) constraints[i] = new Constraints();
34  for(int i = 0;i<nummodels;i++) loads[i] = new Loads();
35  for(int i = 0;i<nummodels;i++) nodes[i] = new Nodes();
36 
37  /*Partition Elements and Nodes*/
39 
40  /*Create elements, vertices and materials, independent of analysis_enum: */
41  CreateElements(elements,iomodel,nummodels);
42  CreateVertices(elements,vertices,iomodel,solution_enum);
43  CreateParameters(parameters,iomodel,rootpath,toolkitfile,solution_enum);
44 
45  /*Should move to CreateInputs2*/
46  Inputs2 *inputs2 = new Inputs2(elements->Size(),vertices->Size());
47  if (iomodel->domaintype != Domain3DsurfaceEnum) iomodel->FetchDataToInput(inputs2,elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
48 
49  /*Can now do Materials since we have created Inputs*/
50  CreateMaterials(elements,inputs2,materials,iomodel,nummodels);
51 
52  /*Update datasets based on each analysis (and add nodes, constrains and loads)*/
53  for(int i=0;i<nummodels;i++){
54 
55  int analysis_enum=analysis_enum_list[i];
56  parameters->AddObject(new IntParam(AnalysisCounterEnum,i));
57 
58  if(VerboseMProcessor()) _printf0_(" creating datasets for analysis " << EnumToStringx(analysis_enum) << "\n");
59  Analysis* analysis = EnumToAnalysis(analysis_enum);
60  analysis->UpdateParameters(parameters,iomodel,solution_enum,analysis_enum);
61  analysis->CreateNodes(nodes[i],iomodel);
62  analysis->UpdateElements(elements,inputs2,iomodel,i,analysis_enum);
63  analysis->CreateConstraints(constraints[i],iomodel);
64  analysis->CreateLoads(loads[i],iomodel);
65  delete analysis;
66 
67  /*Tell datasets that Ids are already sorted*/
68  constraints[i]->Presort();
69  loads[i]->Presort();
70  nodes[i]->Presort();
71 
72  /*Finalize loads (count pengrids,penpairs,rifts,etc)*/
73  loads[i]->Finalize();
74  }
75 
76  /*Solution specific updates*/
77  if(VerboseMProcessor()) _printf0_(" updating elements and materials for control parameters" << "\n");
78  UpdateElementsAndMaterialsControl(elements,parameters,inputs2,materials,iomodel);
79  #ifdef _HAVE_DAKOTA_
80  if(VerboseMProcessor()) _printf0_(" updating elements and materials for uncertainty quantification" << "\n");
81  UpdateElementsAndMaterialsDakota(elements,inputs2,materials,iomodel);
82  #endif
83  if(solution_enum==TransientSolutionEnum) UpdateElementsTransient(elements,parameters,inputs2,iomodel);
84 
85  /*Output definitions dataset: */
86  if(VerboseMProcessor()) _printf0_(" creating output definitions" << "\n");
87  CreateOutputDefinitions(elements,parameters,inputs2,iomodel);
88 
89  /* Sort datasets:
90  * All our datasets are already ordered by ids. Set presort flag so that
91  * later on, when sorting is requested on these datasets, it will not be
92  * redone: */
93  elements->Presort();
94  vertices->Presort();
95  materials->Presort();
96 
97  /*Assign output pointers:*/
98  *pelements = elements;
99  *pnodes = nodes;
100  *pvertices = vertices;
101  *pmaterials = materials;
102  *pconstraints = constraints;
103  *ploads = loads;
104  *pparameters = parameters;
105  *pinputs2 = inputs2;
106 
107  if(VerboseMProcessor()) _printf0_(" done with model processor \n");
108 }

◆ CreateElements()

void CreateElements ( Elements elements,
IoModel iomodel,
int  nummodels 
)

Definition at line 35 of file CreateElementsVerticesAndMaterials.cpp.

35  {/*{{{*/
36 
37  /*Intermediary*/
38  bool control_analysis;
39  bool adolc_analysis;
40 
41  /*Fetch parameters: */
42  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
43  iomodel->FindConstant(&adolc_analysis,"md.autodiff.isautodiff");
44 
45  /*Did we already create the elements? : */
46  _assert_(elements->Size()==0);
47 
48  /*Create elements*/
49  if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
50  if(iomodel->domaintype==Domain2DverticalEnum || iomodel->domaindim==3) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
51 
52  int count = 0;
53  switch(iomodel->meshelementtype){
54  case TriaEnum:
55  for(int i=0;i<iomodel->numberofelements;i++){
56  if(iomodel->my_elements[i]){
57  elements->AddObject(new Tria(i+1,i,count,iomodel,nummodels));
58  count++;
59  }
60  }
61  break;
62  case TetraEnum:
63  for(int i=0;i<iomodel->numberofelements;i++){
64  if(iomodel->my_elements[i]){
65  elements->AddObject(new Tetra(i+1,i,count,iomodel,nummodels));
66  count++;
67  }
68  }
69  break;
70  case PentaEnum:
71  iomodel->FetchData(2,"md.mesh.upperelements","md.mesh.lowerelements");
72  for(int i=0;i<iomodel->numberofelements;i++){
73  if(iomodel->my_elements[i]){
74  elements->AddObject(new Penta(i+1,i,count,iomodel,nummodels));
75  count++;
76  }
77  }
78  break;
79  default:
80  _error_("Mesh not supported yet");
81  }
82 
83  /*Free data: */
84  iomodel->DeleteData(6,"md.mesh.upperelements","md.mesh.lowerelements","md.inversion.min_parameters","md.inversion.max_parameters","md.mesh.vertexonbase","md.mesh.vertexonsurface");
85 
86 }/*}}}*/

◆ CreateMaterials()

void CreateMaterials ( Elements elements,
Inputs2 inputs2,
Materials materials,
IoModel iomodel,
int  nummodels 
)

Definition at line 87 of file CreateElementsVerticesAndMaterials.cpp.

87  {/*{{{*/
88 
89  /*Intermediary*/
90  int i;
91  int nnat,dummy;
92  int* nature=NULL;
93 
94  /*Fetch parameters: */
95  int materials_type;
96  iomodel->FindConstant(&materials_type,"md.materials.type");
97 
98  /*Did we already create the materials? : */
99  _assert_(materials->Size()==0);
100 
101  /*Create materials*/
102  switch(materials_type){
103  case MaticeEnum:
104  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
105  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
106  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
107  switch(iomodel->domaindim){
108  case 2:
110  break;
111  case 3:
112  break;
113  default:
114  _error_("Mesh not supported yet");
115  }
116  break;
117  case MatenhancediceEnum:
118  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
119  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
120  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
121  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
122  switch(iomodel->domaindim){
123  case 2:
126  break;
127  case 3:
128  break;
129  default:
130  _error_("Mesh not supported yet");
131  }
132  break;
133  case MatdamageiceEnum:
134  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
135  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
136  iomodel->FetchDataToInput(inputs2,elements,"md.damage.D",DamageDEnum);
137  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
138  switch(iomodel->domaindim){
139  case 2:
142  break;
143  case 3:
144  break;
145  default:
146  _error_("Mesh not supported yet");
147  }
148  break;
149  case MatestarEnum:
150  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
151  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
152  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
153  for(i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matestar(i+1,i,iomodel));
154  switch(iomodel->domaindim){
155  case 2:
159  break;
160  case 3:
161  break;
162  default:
163  _error_("Mesh not supported yet");
164  }
165  break;
166  case MaterialsEnum:
167 
168  //we have several types of materials. Retrieve this info first:
169  iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
170 
171  //make sure materials that are not tied to elements come last: for now, only Matlitho qualifies.
172  for(int i=0;i<nnat;i++){
173  if (IoCodeToEnumNature(nature[i])==MatlithoEnum){
174  int temp=nature[nnat-1];
175  nature[nnat-1]=nature[i];
176  nature[i]=temp;
177  }
178  }
179 
180  //go through list of materials, and create them:
181  for(int i=0;i<nnat;i++){
182  switch(IoCodeToEnumNature(nature[i])){ //{{{
183  case MaticeEnum:
184  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
185  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
186  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
187  switch(iomodel->domaindim){
188  case 2:
190  break;
191  case 3:
192  break;
193  default:
194  _error_("Mesh not supported yet");
195  }
196  break;
197  case MatlithoEnum:
198  iomodel->FetchData(9,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.isburgers","md.materials.issolid","md.materials.density");
199  materials->AddObject(new Matlitho(materials->Size()+1,iomodel));
200  iomodel->DeleteData(9,"md.materials.radius","md.materials.viscosity","md.materials.lame_lambda","md.materials.lame_mu","md.materials.burgers_viscosity","md.materials.burgers_mu","md.materials.isburgers","md.materials.issolid","md.materials.density");
201  break;
202 
203  case MathydroEnum:
204  {
205  /*If we don't have any materials pointed to by elements (meaning, if we are running only litho or hydro),
206  * then we need to zero out the hmaterial pointers inside the elements dataset so that it won't error out
207  * during configuration: */
208  bool isice=false;
209  for (int j=0;j<nnat;j++){
210  if((IoCodeToEnumNature(nature[i])==MaticeEnum)||
211  (IoCodeToEnumNature(nature[i])==MatenhancediceEnum)||
212  (IoCodeToEnumNature(nature[i])==MatestarEnum)||
213  (IoCodeToEnumNature(nature[i])==MatdamageiceEnum)){
214  isice=true; break; }
215  }
216  if (!isice){
217  /*go through elements, and zero the hmaterials pointers: */
218  for(int j=0;j<elements->Size();j++){
219  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
220  switch(element->ObjectEnum()){
221  case TriaEnum:
222  {
223  Tria* tria= xDynamicCast<Tria*>(element);
224  tria->hmaterial=NULL;
225  }
226  break;
227  default:
228  _error_("Not implemented yet");
229  }
230  }
231  }
232  }
233  break;
234 
235  case MatenhancediceEnum:
236  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
237  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
238  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_E",MaterialsRheologyEEnum);
239  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
240  switch(iomodel->domaindim){
241  case 2:
244  break;
245  case 3:
246  break;
247  default:
248  _error_("Mesh not supported yet");
249  }
250  break;
251  case MatdamageiceEnum:
252  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
253  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_n",MaterialsRheologyNEnum);
254  iomodel->FetchDataToInput(inputs2,elements,"md.damage.D",DamageDEnum);
255  for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel));
256  switch(iomodel->domaindim){
257  case 2:
260  break;
261  case 3:
262  break;
263  default:
264  _error_("Mesh not supported yet");
265  }
266  break;
267  case MatestarEnum:
268  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",MaterialsRheologyBEnum);
269  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Ec",MaterialsRheologyEcEnum);
270  iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_Es",MaterialsRheologyEsEnum);
271  for(i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matestar(i+1,i,iomodel));
272  switch(iomodel->domaindim){
273  case 2:
277  break;
278  case 3:
279  break;
280  default:
281  _error_("Mesh not supported yet");
282  }
283  break;
284 
285  default:
286  _error_("Materials nature type "<<EnumToStringx(IoCodeToEnumNature(nature[i]))<<" not supported");
287 
288  } //}}}
289  }
290  //Free ressources:
291  xDelete<int>(nature);
292  break;
293 
294  default:
295  _error_("Materials "<<EnumToStringx(materials_type)<<" not supported");
296  }
297 
298  /*Free data: */
299  iomodel->DeleteData(3,"md.material.rheology_B","md.material.rheology_n","md.damage.D");
300 }/*}}}*/

◆ CreateVertices()

void CreateVertices ( Elements elements,
Vertices vertices,
IoModel iomodel,
int  solution_type,
bool  isamr = false 
)

Definition at line 301 of file CreateElementsVerticesAndMaterials.cpp.

301  {/*{{{*/
302 
303  /*Get element partitionning*/
304  int* epart = iomodel->epart;
305 
306  /*Determine element width*/
307  int elements_width;
308  switch(iomodel->meshelementtype){
309  case TriaEnum: elements_width=3; break;
310  case TetraEnum: elements_width=4; break;
311  case PentaEnum: elements_width=6; break;
312  default: _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
313  }
314 
315  /*Get my_rank:*/
316  int my_rank = IssmComm::GetRank();
317  int num_procs = IssmComm::GetSize();
318 
319  /*create matrix that keeps track of all ranks that have vertex i, and initialize as -1 (Common to all CPUs)*/
320  int* vertices_ranks = xNew<int>(MAXCONNECTIVITY*iomodel->numberofvertices);
321  for(int i=0;i<MAXCONNECTIVITY*iomodel->numberofvertices;i++) vertices_ranks[i] = -1;
322 
323  /*For all vertices, count how many cpus hold vertex i (initialize with 0)*/
324  int* vertices_proc_count = xNewZeroInit<int>(iomodel->numberofvertices);
325 
326  /*Go through all elements and mark all vertices for all partitions*/
327  for(int i=0;i<iomodel->numberofelements;i++){
328  for(int j=0;j<elements_width;j++){
329  /*Get current vertex sid*/
330  int vid = iomodel->elements[elements_width*i+j]-1;
331  AddVertexToRank(vertices_ranks,vertices_proc_count,vid,epart[i]);
332  }
333  }
334 
335  /*Take care of penalties (only in non-AMR for now)*/
336  if(!isamr){
337  int numvertex_pairing;
338  int *vertex_pairing = NULL;
339  iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.stressbalance.vertex_pairing");
340  for(int i=0;i<numvertex_pairing;i++){
341  int id1 = vertex_pairing[2*i+0]-1;
342  int id2 = vertex_pairing[2*i+1]-1;
343  for(int e=0;e<num_procs;e++){
344  if(IsVertexInRank(vertices_ranks,vertices_proc_count,id1,e)){
345  AddVertexToRank(vertices_ranks,vertices_proc_count,id2,e);
346  }
347  }
348  }
349  xDelete<int>(vertex_pairing);
350  iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.masstransport.vertex_pairing");
351  for(int i=0;i<numvertex_pairing;i++){
352  int id1 = vertex_pairing[2*i+0]-1;
353  int id2 = vertex_pairing[2*i+1]-1;
354  for(int e=0;e<num_procs;e++){
355  if(IsVertexInRank(vertices_ranks,vertices_proc_count,id1,e)){
356  AddVertexToRank(vertices_ranks,vertices_proc_count,id2,e);
357  }
358  }
359  }
360  xDelete<int>(vertex_pairing);
361  }
362 
363  /*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/
364  int offset = 0;
365  int* vertices_offsets = xNew<int>(iomodel->numberofvertices);
366  for(int i=0;i<iomodel->numberofvertices;i++){
367  if(IsVertexInRank(vertices_ranks,vertices_proc_count,i,my_rank)){
368  vertices_offsets[i] = offset++;
369  }
370  else{
371  vertices_offsets[i] = -1;
372  }
373  }
374 
375  /*Now, Count how many clones we have with other partitions*/
376  int* common_send = xNew<int>(num_procs);
377  int* common_recv = xNew<int>(num_procs);
378  int** common_send_ids = xNew<int*>(num_procs);
379  int** common_recv_ids = xNew<int*>(num_procs);
380 
381  /*First step: allocate, Step 2: populate*/
382  for(int step=0;step<2;step++){
383 
384  if(step==1){
385  /*Allocate send and receive arrays of ids now*/
386  for(int i=0;i<num_procs;i++){
387  _assert_(common_send[i]>=0 && common_recv[i]>=0);
388  common_send_ids[i] = xNew<int>(common_send[i]);
389  common_recv_ids[i] = xNew<int>(common_recv[i]);
390  }
391  }
392 
393  /*Re/Initialize counters to 0*/
394  for(int i=0;i<num_procs;i++){
395  common_recv[i]=0;
396  common_send[i]=0;
397  }
398 
399  /*Go through table and find clones/masters etc*/
400  for(int i=0;i<iomodel->numberofvertices;i++){
401 
402  /*If we did not find this vertex in our current partition, go to next vertex*/
403  if(vertices_offsets[i] == -1) continue;
404 
405  /*Find in what column this rank belongs*/
406  int col = -1;
407  for(int j=0;j<MAXCONNECTIVITY;j++){
408  if(vertices_ranks[MAXCONNECTIVITY*i+j] == my_rank){
409  col = j;
410  break;
411  }
412  }
413  _assert_(col!=-1);
414 
415  /*If col==0, it is either not on boundary, or a master*/
416  if(col==0){
417  /*1. is this vertex on the boundary? Skip if not*/
418  if(vertices_ranks[MAXCONNECTIVITY*i+col+1]==-1){
419  continue;
420  }
421  else{
422  for(int j=1;j<vertices_proc_count[i];j++){
423  _assert_(vertices_ranks[MAXCONNECTIVITY*i+j]>=0);
424  int rank = vertices_ranks[MAXCONNECTIVITY*i+j];
425  if(step==1){
426  common_send_ids[rank][common_send[rank]] = vertices_offsets[i];
427  }
428  common_send[rank]++;
429  }
430  }
431  }
432  else{
433  /*3. It is a slave, record that we need to receive for this cpu*/
434  int rank = vertices_ranks[MAXCONNECTIVITY*i+0];
435  if(step==1){
436  common_recv_ids[rank][common_recv[rank]] = vertices_offsets[i];
437  }
438  common_recv[rank]++;
439  }
440  }
441  }
442 
443  /*Create Vertices, depending on the constructor type: */
444  if(solution_type!=LoveSolutionEnum) CreateNumberNodeToElementConnectivity(iomodel);
445  if(!isamr){
446  bool isoceancoupling;
447  iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling");
448 
449  iomodel->FetchData(6,"md.mesh.x","md.mesh.y","md.mesh.z","md.geometry.base","md.geometry.thickness","md.mask.ice_levelset");
450  if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
451  if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
452 
453  for(int i=0;i<iomodel->numberofvertices;i++){
454  if(vertices_offsets[i]!=-1){
455  bool isclone = (vertices_ranks[MAXCONNECTIVITY*i+0]!=my_rank);
456  vertices->AddObject(new Vertex(i+1,i,isclone,iomodel,isamr));
457  }
458  }
459 
460  /*Free data: */
461  iomodel->DeleteData(6,"md.mesh.x","md.mesh.y","md.mesh.z","md.geometry.base","md.geometry.thickness","md.mask.ice_levelset");
462  if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
463  if (isoceancoupling) iomodel->DeleteData(2,"md.mesh.lat","md.mesh.long");
464  }
465  else{
466  for(int i=0;i<iomodel->numberofvertices;i++){
467  if(vertices_offsets[i]!=-1){
468  bool isclone = (vertices_ranks[MAXCONNECTIVITY*i+0]!=my_rank);
469  vertices->AddObject(new Vertex(i+1,i,isclone,iomodel,isamr));
470  }
471  }
472  }
473  xDelete<int>(vertices_offsets);
474 
475  /*Final step, create my_vertices*/
476  _assert_(!iomodel->my_vertices);
477  iomodel->my_vertices = xNew<bool>(iomodel->numberofvertices);
478  for(int i=0;i<iomodel->numberofvertices;i++){
479  if(IsVertexInRank(vertices_ranks,vertices_proc_count,i,my_rank)){
480  iomodel->my_vertices[i] = true;
481  }
482  else{
483  iomodel->my_vertices[i] = false;
484  }
485  }
486 
487  /*Free data: */
488  xDelete<int>(vertices_ranks);
489  xDelete<int>(vertices_proc_count);
490 
491  /*Assign communicators*/
492  vertices->common_send=common_send;
493  vertices->common_recv=common_recv;
494  vertices->common_send_ids=common_send_ids;
495  vertices->common_recv_ids=common_recv_ids;
496 
497  /*Finalize Initialization*/
498  vertices->Finalize(iomodel);
499 }/*}}}*/

◆ CreateParameters()

void CreateParameters ( Parameters parameters,
IoModel iomodel,
char *  rootpath,
FILE *  toolkitfile,
const int  solution_type 
)

Definition at line 18 of file CreateParameters.cpp.

18  {
19 
20  int i,j,m,k;
21  int numoutputs,basalforcing_model,timestepping_type;
22  char** requestedoutputs = NULL;
23  char** outputonnodes = NULL;
24  char* fieldname = NULL;
25  IssmDouble time;
26 
27  /*parameters for mass flux:*/
28  int mass_flux_num_profiles = 0;
29  bool qmu_mass_flux_present = false;
30  bool autodiff_mass_flux_present = false;
31  bool mass_flux_present = false;
32  bool interp;
33  IssmDouble **array = NULL;
34  int *mdims_array = NULL;
35  int *ndims_array = NULL;
36  IssmDouble *temp_matrix = NULL;
37  int temp_m,temp_n;
38  IssmDouble *matrix = NULL;
39  int count;
40 
41  IssmDouble *temp = NULL;
42  IssmDouble *transparam = NULL;
43  IssmDouble yts;
44  int N,M;
45 
46  /*Copy some constants from iomodel */
47  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.domain_type",DomainTypeEnum));
48  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.domain_dimension",DomainDimensionEnum));
49  parameters->AddObject(iomodel->CopyConstantObject("md.settings.output_frequency",SettingsOutputFrequencyEnum));
50  parameters->AddObject(iomodel->CopyConstantObject("md.settings.sb_coupling_frequency",SettingsSbCouplingFrequencyEnum));
51  parameters->AddObject(iomodel->CopyConstantObject("md.settings.recording_frequency",SettingsRecordingFrequencyEnum));
52  parameters->AddObject(iomodel->CopyConstantObject("md.constants.yts",ConstantsYtsEnum));
53  parameters->AddObject(iomodel->CopyConstantObject("md.debug.profiling",DebugProfilingEnum));
54  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.average_vertex_connectivity",MeshAverageVertexConnectivityEnum));
55  parameters->AddObject(iomodel->CopyConstantObject("md.settings.waitonlock",SettingsWaitonlockEnum));
56  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberofvertices",MeshNumberofverticesEnum));
57  parameters->AddObject(iomodel->CopyConstantObject("md.settings.io_gather",SettingsIoGatherEnum));
58  parameters->AddObject(iomodel->CopyConstantObject("md.settings.solver_residue_threshold",SettingsSolverResidueThresholdEnum));
59  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.isautodiff",AutodiffIsautodiffEnum));
60  parameters->AddObject(iomodel->CopyConstantObject("md.qmu.isdakota",QmuIsdakotaEnum));
61  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.iscontrol",InversionIscontrolEnum));
62  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
63  parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
64  parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.parameterization",FrontalForcingsParamEnum));
65  parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1));
66 
67  {/*This is specific to ice...*/
68  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));
69  parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.reltol",SteadystateReltolEnum));
70  parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.maxiter",SteadystateMaxiterEnum));
71  parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.migration",GroundinglineMigrationEnum));
72  parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.friction_interpolation",GroundinglineFrictionInterpolationEnum));
73  parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.melt_interpolation",GroundinglineMeltInterpolationEnum));
74  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isstressbalance",TransientIsstressbalanceEnum));
75  parameters->AddObject(iomodel->CopyConstantObject("md.transient.ismasstransport",TransientIsmasstransportEnum));
76  parameters->AddObject(iomodel->CopyConstantObject("md.transient.issmb",TransientIssmbEnum));
77  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isthermal",TransientIsthermalEnum));
78  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgroundingline",TransientIsgroundinglineEnum));
79  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgia",TransientIsgiaEnum));
80  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isesa",TransientIsesaEnum));
81  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isdamageevolution",TransientIsdamageevolutionEnum));
82  parameters->AddObject(iomodel->CopyConstantObject("md.transient.ishydrology",TransientIshydrologyEnum));
83  parameters->AddObject(iomodel->CopyConstantObject("md.transient.ismovingfront",TransientIsmovingfrontEnum));
84  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isslr",TransientIsslrEnum));
85  parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum));
86  parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum));
87  parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));
88 
89  /*For stress balance only*/
90  parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));
91  parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.rift_penalty_threshold",StressbalanceRiftPenaltyThresholdEnum));
92  parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.maxiter",StressbalanceMaxiterEnum));
93  parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.restol",StressbalanceRestolEnum));
94  parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.reltol",StressbalanceReltolEnum));
95  parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.abstol",StressbalanceAbstolEnum));
96  if(iomodel->domaintype==Domain3DEnum)
97  parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberoflayers",MeshNumberoflayersEnum));
98  }
99 
100  /*amr properties*/
101  int amrtype,amr_frequency;
102  iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency");
103  if(solution_type==TransientSolutionEnum && amr_frequency){
104  /*Load common amr parameters*/
105  parameters->AddObject(iomodel->CopyConstantObject("md.amr.type",AmrTypeEnum));
106  parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum));
107  parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum));
108  parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum));
109  parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_groupthreshold",AmrThicknessErrorGroupThresholdEnum));
110  parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum));
111  parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum));
112  parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_groupthreshold",AmrDeviatoricErrorGroupThresholdEnum));
113  parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum));
114  parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum));
115  /*Load specific amr parameters*/
116  iomodel->FindConstant(&amrtype,"md.amr.type");
117  switch(amrtype){
118  #ifdef _HAVE_NEOPZ_
119  case AmrNeopzEnum:
120  parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum));
121  parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum));
122  parameters->AddObject(iomodel->CopyConstantObject("md.amr.lag",AmrLagEnum));
123  break;
124  #endif
125 
126  #ifdef _HAVE_BAMG_
127  case AmrBamgEnum:
128  parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmin",AmrHminEnum));
129  parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum));
130  parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum));
131  parameters->AddObject(iomodel->CopyConstantObject("md.amr.keepmetric",AmrKeepMetricEnum));
132  parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum));
133  parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_resolution",AmrGroundingLineResolutionEnum));
134  parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_resolution",AmrIceFrontResolutionEnum));
135  parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum));
136  parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum));
137  /*Convert fieldname to enum and put it in params*/
138  iomodel->FindConstant(&fieldname,"md.amr.fieldname");
139  parameters->AddObject(new IntParam(AmrFieldEnum,StringToEnumx(fieldname)));
140  xDelete<char>(fieldname);
141  break;
142  #endif
143 
144  default:
145  _error_("Adaptive mesh refinement "<<EnumToStringx(amrtype)<<" not implemented yet");
146  }
147  }
148 
149  /*Basal forcing parameters*/
150  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.model",BasalforcingsEnum));
151  iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
152  switch(basalforcing_model){
154  /*Nothing to add to parameters*/
155  break;
157  iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
158  iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate");
159  if(N==1){
160  _assert_(M==1);
161  parameters->AddObject(new DoubleParam(BasalforcingsDeepwaterMeltingRateEnum,transparam[0]));
162  }
163  else{
164  _assert_(N==2);
165  parameters->AddObject(new TransientParam(BasalforcingsDeepwaterMeltingRateEnum,&transparam[0],&transparam[M],interp,M));
166  }
167  xDelete<IssmDouble>(transparam);
168  iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate");
169  if(N==1){
170  _assert_(M==1);
171  parameters->AddObject(new DoubleParam(BasalforcingsUpperwaterMeltingRateEnum,transparam[0]));
172  }
173  else{
174  _assert_(N==2);
175  parameters->AddObject(new TransientParam(BasalforcingsUpperwaterMeltingRateEnum,&transparam[0],&transparam[M],interp,M));
176  }
177  xDelete<IssmDouble>(transparam);
178  iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation");
179  if(N==1){
180  _assert_(M==1);
181  parameters->AddObject(new DoubleParam(BasalforcingsDeepwaterElevationEnum,transparam[0]));
182  }
183  else{
184  _assert_(N==2);
185  parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M));
186  }
187  xDelete<IssmDouble>(transparam);
188  iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation");
189  if(N==1){
190  _assert_(M==1);
191  parameters->AddObject(new DoubleParam(BasalforcingsUpperwaterElevationEnum,transparam[0]));
192  }
193  else{
194  _assert_(N==2);
195  parameters->AddObject(new TransientParam(BasalforcingsUpperwaterElevationEnum,&transparam[0],&transparam[M],interp,M));
196  }
197  xDelete<IssmDouble>(transparam);
198  break;
200  /*Nothing to add to parameters:*/
201  break;
203  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));
204  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.threshold_thickness",BasalforcingsThresholdThicknessEnum));
205  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.upperdepth_melt",BasalforcingsUpperdepthMeltEnum));
206  break;
208  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.mantleconductivity",BasalforcingsMantleconductivityEnum));
209  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.nusselt",BasalforcingsNusseltEnum));
210  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.dtbg",BasalforcingsDtbgEnum));
211  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumeradius",BasalforcingsPlumeradiusEnum));
212  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.topplumedepth",BasalforcingsTopplumedepthEnum));
213  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.bottomplumedepth",BasalforcingsBottomplumedepthEnum));
214  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumex",BasalforcingsPlumexEnum));
215  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumey",BasalforcingsPlumeyEnum));
216  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.crustthickness",BasalforcingsCrustthicknessEnum));
217  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.uppercrustthickness",BasalforcingsUppercrustthicknessEnum));
218  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.uppercrustheat",BasalforcingsUppercrustheatEnum));
219  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.lowercrustheat",BasalforcingsLowercrustheatEnum));
220  break;
222  iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
223  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsPicoNumBasinsEnum));
224  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.maxboxcount",BasalforcingsPicoMaxboxcountEnum));
225  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_T",BasalforcingsPicoGammaTEnum));
226  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isplume",BasalforcingsPicoIsplumeEnum));
227  iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature");
228  _assert_(M>=1 && N>=1);
229  parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M));
230  xDelete<IssmDouble>(transparam);
231  iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity");
232  _assert_(M>=1 && N>=1);
233  parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M));
234  xDelete<IssmDouble>(transparam);
235  break;
237  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum));
238  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum));
239  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum));
240  iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t");
241  parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N));
242  xDelete<IssmDouble>(transparam);
243  iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.tf_depths");
244  parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6TfDepthsEnum,transparam,N));
245  xDelete<IssmDouble>(transparam);
246  break;
248  parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));
249  break;
250  default:
251  _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
252  }
253 
254  /*some parameters that did not come with the iomodel: */
255  parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type));
256 
257  /*Time stepping*/
258  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.type",TimesteppingTypeEnum));
259  iomodel->FindConstant(&timestepping_type,"md.timestepping.type");
260  switch(timestepping_type){
262  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.start_time",TimesteppingStartTimeEnum));
263  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.final_time",TimesteppingFinalTimeEnum));
264  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step",TimesteppingTimeStepEnum));
265  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.interp_forcings",TimesteppingInterpForcingsEnum));
266  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.coupling_time",TimesteppingCouplingTimeEnum));
267  break;
269  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.start_time",TimesteppingStartTimeEnum));
270  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.final_time",TimesteppingFinalTimeEnum));
271  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step_min",TimesteppingTimeStepMinEnum));
272  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step_max",TimesteppingTimeStepMaxEnum));
273  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.cfl_coefficient",TimesteppingCflCoefficientEnum));
274  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.interp_forcings",TimesteppingInterpForcingsEnum));
275  parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.coupling_time",TimesteppingCouplingTimeEnum));
276  break;
277  default:
278  _error_("Time stepping \""<<EnumToStringx(timestepping_type)<<"\" not supported yet");
279  }
280  iomodel->FindConstant(&time,"md.timestepping.start_time");
281  parameters->AddObject(new DoubleParam(TimeEnum,time));
282  parameters->AddObject(new IntParam(StepEnum,0));
283 
284  /*By default, save all results*/
285  parameters->AddObject(new BoolParam(SaveResultsEnum,true));
286 
287  /*Should we output results on nodes?*/
288  iomodel->FindConstant(&outputonnodes,&numoutputs,"md.settings.results_on_nodes");
289  parameters->AddObject(new IntParam(SettingsNumResultsOnNodesEnum,numoutputs));
290  if(numoutputs)parameters->AddObject(new StringArrayParam(SettingsResultsOnNodesEnum,outputonnodes,numoutputs));
291  iomodel->DeleteData(&outputonnodes,numoutputs,"md.settings.results_on_nodes");
292 
293  /*Requested outputs */
294  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.transient.requested_outputs");
295  parameters->AddObject(new IntParam(TransientNumRequestedOutputsEnum,numoutputs));
296  if(numoutputs)parameters->AddObject(new StringArrayParam(TransientRequestedOutputsEnum,requestedoutputs,numoutputs));
297  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.transient.requested_outputs");
298 
299  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.steadystate.requested_outputs");
300  parameters->AddObject(new IntParam(SteadystateNumRequestedOutputsEnum,numoutputs));
301  if(numoutputs)parameters->AddObject(new StringArrayParam(SteadystateRequestedOutputsEnum,requestedoutputs,numoutputs));
302  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.steadystate.requested_outputs");
303 
304  int materialstype;
305  iomodel->FindConstant(&materialstype,"md.materials.type");
306 
307  switch(materialstype){
308  case MaticeEnum:
309  case MatdamageiceEnum:
310  case MatenhancediceEnum:
311  case MatestarEnum:
312  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
313  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
314  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
315  parameters->AddObject(iomodel->CopyConstantObject("md.materials.mu_water",MaterialsMuWaterEnum));
316  parameters->AddObject(iomodel->CopyConstantObject("md.materials.heatcapacity",MaterialsHeatcapacityEnum));
317  parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
318  parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
319  parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
320  parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
321  parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));
322  parameters->AddObject(iomodel->CopyConstantObject("md.materials.meltingpoint",MaterialsMeltingpointEnum));
323  parameters->AddObject(iomodel->CopyConstantObject("md.constants.referencetemperature",ConstantsReferencetemperatureEnum));
324  parameters->AddObject(iomodel->CopyConstantObject("md.materials.mixed_layer_capacity",MaterialsMixedLayerCapacityEnum));
325  parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum));
326  parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum));
327  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));
328  parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
329 
330  break;
331  case MaterialsEnum:{
332  int nnat,dummy;
333  int* nature=NULL;
334  iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
335  for(int i=0;i<nnat;i++){
336  switch(IoCodeToEnumNature(nature[i])){
337  case MatlithoEnum:
338  break;
339  case MaticeEnum:
340  case MatdamageiceEnum:
341  case MatenhancediceEnum:
342  case MatestarEnum:
343  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
344  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
345  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
346  parameters->AddObject(iomodel->CopyConstantObject("md.materials.mu_water",MaterialsMuWaterEnum));
347  parameters->AddObject(iomodel->CopyConstantObject("md.materials.heatcapacity",MaterialsHeatcapacityEnum));
348  parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
349  parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
350  parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
351  parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
352  parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));
353  parameters->AddObject(iomodel->CopyConstantObject("md.materials.meltingpoint",MaterialsMeltingpointEnum));
354  parameters->AddObject(iomodel->CopyConstantObject("md.constants.referencetemperature",ConstantsReferencetemperatureEnum));
355  parameters->AddObject(iomodel->CopyConstantObject("md.materials.mixed_layer_capacity",MaterialsMixedLayerCapacityEnum));
356  parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum));
357  parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum));
358  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));
359 
360  /*slr:*/
361  parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
362  break;
363  case MathydroEnum:
364  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
365  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
366  parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
367  parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
368  break;
369  }
370  }
371  /*Free rssources:*/
372  xDelete<int>(nature);
373  break;
374 
375  }
376  default:
377  _error_("Material "<< EnumToStringx(materialstype) <<" not supported yet");
378  }
379 
380  int smb_model;
381  iomodel->FindConstant(&smb_model,"md.smb.model");
382  switch(smb_model){
383  case SMBforcingEnum:
384  case SMBgradientsEnum:
385  case SMBgradientselaEnum:
386  case SMBhenningEnum:
387  case SMBcomponentsEnum:
390  /*Nothing to add*/
391  break;
392  case SMBgembEnum:
393  parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
394  parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
395  break;
396  case SMBpddEnum:
397  parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
398  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
399  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
400  break;
401  case SMBpddSicopolisEnum:
402  parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
403  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
404  break;
405  case SMBd18opddEnum:
406  parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
407  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
408  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
409  parameters->AddObject(iomodel->CopyConstantObject("md.smb.dpermil",SmbDpermilEnum));
410  break;
411  case SMBsemicEnum:
412  parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
413  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
414  parameters->AddObject(iomodel->CopyConstantObject("md.smb.rdl",SmbRdlEnum));
415  break;
416  default:
417  _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
418  }
419 
420  int hydrology_model;
421  iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
422  if(hydrology_model==HydrologydcEnum){
423  /*FIXME: this cshould go to Analysis!!!*/
424  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum));
425  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_porosity",HydrologydcSedimentPorosityEnum));
426  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_thickness",HydrologydcSedimentThicknessEnum));
427  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum));
428  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isefficientlayer",HydrologydcIsefficientlayerEnum));
429 
430  bool isefficientlayer;
431  iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
432  if(isefficientlayer){
433  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_compressibility",HydrologydcEplCompressibilityEnum));
434  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum));
435  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_initial_thickness",HydrologydcEplInitialThicknessEnum));
436  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_colapse_thickness",HydrologydcEplColapseThicknessEnum));
437  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_max_thickness",HydrologydcEplMaxThicknessEnum));
438  parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_conductivity",HydrologydcEplConductivityEnum));
439  }
440  }
441  else if(hydrology_model==HydrologyshreveEnum){
442  /*Nothing to add*/
443  }
444  else if(hydrology_model==HydrologyshaktiEnum){
445  /*Nothing to add*/
446  }
447  else if(hydrology_model==HydrologypismEnum){
448  /*Nothing to add*/
449  }
450  else if(hydrology_model==HydrologyGlaDSEnum){
451  /*Nothing to add*/
452  }
453  else{
454  _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
455  }
456 
457  if(materialstype==MatdamageiceEnum){
458  iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.damage.requested_outputs");
459  parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
460  if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
461  iomodel->DeleteData(&requestedoutputs,numoutputs,"md.damage.requested_outputs");
462  }
463 
464  /*Deal with mass flux segments: {{{*/
465  iomodel->FetchData(&qmu_mass_flux_present,"md.qmu.mass_flux_segments_present");
466  iomodel->FetchData(&autodiff_mass_flux_present,"md.autodiff.mass_flux_segments_present");
467 
468  if(qmu_mass_flux_present || autodiff_mass_flux_present)mass_flux_present=true;
469  else mass_flux_present=false;
470  parameters->AddObject(new BoolParam(MassFluxSegmentsPresentEnum,mass_flux_present));
471 
472  if(mass_flux_present){
473 
474  /*Fetch the mass flux segments necessary to compute the mass fluxes. Build a DoubleMatArrayParam object out of them: */
475  iomodel->FetchData(&array,&mdims_array,&ndims_array,&mass_flux_num_profiles,"md.qmu.mass_flux_segments");
476  if(mass_flux_num_profiles==0)_error_("mass_flux_num_profiles is 0, when MassFlux computations were requested!");
477 
478  /*Go through segments, and extract those that belong to this cpu: */
479  for(i=0;i<mass_flux_num_profiles;i++){
480  temp_matrix=array[i];
481  temp_m=mdims_array[i];
482  temp_n=ndims_array[i];
483  _assert_(temp_n==5);
484 
485  m=0;
486  for(j=0;j<temp_m;j++){
487  if ( iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1] )m++;
488  }
489  if(m){
490  matrix=xNewZeroInit<IssmDouble>(5*m);
491  count=0;
492  for(j=0;j<temp_m;j++){
493  if (iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1]){
494  for(k=0;k<5;k++)*(matrix+5*count+k)=*(temp_matrix+5*j+k);
495  count++;
496  }
497  }
498  }
499  else{
500  matrix=NULL;
501  }
502 
503  /*Assign: */
504  array[i]=matrix;
505  mdims_array[i]=m;
506  ndims_array[i]=5;
507 
508  /*Free temporary matrix: */
509  xDelete<IssmDouble>(temp_matrix);
510  }
511 
512  /*Ok, we have an array of segments, different on every cpu. Create a DoubleMatArrayParam object with it: */
513  parameters->AddObject(new DoubleMatArrayParam(MassFluxSegmentsEnum,array,mass_flux_num_profiles,mdims_array,ndims_array));
514 
515  /*Free data: */
516  for(i=0;i<mass_flux_num_profiles;i++){
517  IssmDouble* matrix=array[i];
518  xDelete<IssmDouble>(matrix);
519  }
520  xDelete<int>(mdims_array);
521  xDelete<int>(ndims_array);
522  xDelete<IssmDouble*>(array);
523  }
524  /*}}}*/
525 
526  /*Before returning, create parameters in case we are running Qmu or control types runs: */
527  CreateParametersControl(parameters,iomodel,solution_type);
528 
529  #ifdef _HAVE_DAKOTA_
530  CreateParametersDakota(parameters,iomodel,rootpath);
531  #endif
532 
533  /*Now, deal with toolkits options, which need to be put into the parameters dataset: */
534  ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);
535 
536  #ifdef _HAVE_AD_
537  if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n");
538  CreateParametersAutodiff(parameters,iomodel);
539  if(VerboseMProcessor()) _printf0_(" ending autodiff parameters \n");
540  #endif
541 }

◆ CreateParametersAutodiff()

void CreateParametersAutodiff ( Parameters parameters,
IoModel iomodel 
)

Definition at line 10 of file CreateParametersAutodiff.cpp.

10  {
11 
12  int i;
13  bool isautodiff;
14  int num_dependent_objects;
15  int num_dep=0;
16  char** names=NULL;
17  int* types=NULL;
18  int dummy;
19  char* autodiff_driver=NULL;
20  int* indices=NULL;
21  int num_indices;
22  char* options=NULL;
23  char* toolkit=NULL;
24 
25  IssmDouble* xp=NULL;
26  IssmDouble* xp_backup=NULL;
27  int num_ind,local_num_ind;
28  DataSet* dependent_objects=NULL;
29 
30  /*retrieve some parameters: */
31  iomodel->FindConstant(&isautodiff,"md.autodiff.isautodiff");
32 
33  #ifdef _HAVE_ADOLC_
34  /*initialize a placeholder to store solver pointers: {{{*/
36 
37  /*Solver pointers depend on what type of solver we are implementing: */
38  options=OptionsFromAnalysis(&toolkit,parameters,DefaultAnalysisEnum);
39  ToolkitOptions::Init(toolkit,options);
40  xDelete<char>(toolkit);
41 
43  case MumpsEnum:{
44  #ifdef _HAVE_MUMPS_
45  theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(mumpsSolveEDF);
46  #else
47  _error_("requesting mumps solver without MUMPS being compiled in!");
48  #endif
49  break;
50  }
51  case GslEnum: {
52  #ifdef _HAVE_GSL_
53  theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx);
54  #else
55  _error_("requesting GSL solver without GSL being compiled in!");
56  #endif
57  break;
58  }
59  default:
60  _error_("solver type not supported yet!");
61  }
62 
63  // to save some space:
64  // we know we won't use adolc inside of the solver:
65  theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->nestedAdolc=false;
66  // the solution vector is just allocated and doesn't have a meaningful prior value
67  theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->dp_y_priorRequired=false;
68  // the solver wrapper makes sure the matrix and the right hand side don't change
69  theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->dp_x_changes=false;
70  parameters->AddObject(theAdolcEDF_p);
71 
72  /*Free ressources: */
73  xDelete<char>(options);
74  /*}}}*/
75  #elif _HAVE_CODIPACK_
76  //fprintf(stderr, "*** Codipack CreateParametersAutodiff()\n");
77  /*initialize a placeholder to store solver pointers: {{{*/
78  /*Solver pointers depend on what type of solver we are implementing: */
79  options=OptionsFromAnalysis(&toolkit,parameters,DefaultAnalysisEnum);
80  ToolkitOptions::Init(toolkit,options);
81  xDelete<char>(toolkit);
82 
84  case MumpsEnum:{
85  #ifndef _HAVE_MUMPS_
86  _error_("CoDiPack: requesting mumps solver without MUMPS being compiled in!");
87  #endif
88  break;
89  }
90  case GslEnum: {
91  #ifndef _HAVE_GSL_
92  _error_("CoDiPack: requesting GSL solver without GSL being compiled in!");
93  #endif
94  break;
95  }
96  default:
97  _error_("solver type not supported yet!");
98  }
99  /*Free ressources: */
100  xDelete<char>(options);
101  #endif
102  #if defined(_HAVE_AD_)
103 
104  if(isautodiff){
105  #if defined(_HAVE_ADOLC_)
106  /*Copy some parameters from IoModel to parameters dataset*/
107  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum));
108  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum));
109  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.lbufsize",AutodiffLbufsizeEnum));
110  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tbufsize",AutodiffTbufsizeEnum));
111  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum));
112  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum));
113 
114  #elif defined(_HAVE_CODIPACK_)
115  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.tapeAlloc",AutodiffTapeAllocEnum));
116 
117  #else
118  _error_("not supported yet");
119  #endif
120 
121  /*retrieve driver: {{{*/
122  iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver");
123  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum));
124 
125  if(strcmp(autodiff_driver,"fos_forward")==0){
126 #if _HAVE_CODIPACK_
127  // FIXME codi support Foward Mode (scalar)
128  _error_("Foward Mode (scalar) not supported yet!");
129 #endif
130  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_forward_index",AutodiffFosForwardIndexEnum));
131  }
132  else if(strcmp(autodiff_driver,"fos_reverse")==0){
133  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_reverse_index",AutodiffFosReverseIndexEnum));
134  }
135  else if(strcmp(autodiff_driver,"fov_forward")==0){
136 #if _HAVE_CODIPACK_
137  // FIXME codi support Foward Mode (vector)
138  _error_("Foward Mode (vector) not supported yet!");
139 #endif
140  /*Retrieve list of indices: */
141  iomodel->FetchData(&indices,&num_indices,&dummy,"md.autodiff.fov_forward_indices");
142  parameters->AddObject(new IntMatParam(AutodiffFovForwardIndicesEnum,indices,num_indices,1));
143  xDelete<int>(indices);
144  }
145  xDelete<char>(autodiff_driver);
146  /*}}}*/
147  /*Deal with dependents first: {{{*/
148  iomodel->FindConstant(&num_dependent_objects,"md.autodiff.num_dependent_objects");
149  dependent_objects=new DataSet();
150  num_dep=0;
151 
152  if(num_dependent_objects){
153  iomodel->FindConstant(&names,&dummy,"md.autodiff.dependent_object_names");
154  iomodel->FetchData(&types,&dummy,&dummy,"md.autodiff.dependent_object_types");
155  iomodel->FetchData(&indices,&dummy,&dummy,"md.autodiff.dependent_object_indices");
156 
157  for(i=0;i<num_dependent_objects;i++){
158  DependentObject* dep=new DependentObject(names[i],types[i],indices[i]);
159  dependent_objects->AddObject(dep);
160  num_dep+=dep->NumDependents();
161  }
162 
163  /*Free ressources:*/
164  for(i=0;i<num_dependent_objects;i++){
165  char* string=names[i]; xDelete<char>(string);
166  }
167  xDelete<char*>(names);
168  xDelete<int>(types);
169  xDelete<int>(indices);
170  }
171  parameters->AddObject(new DataSetParam(AutodiffDependentObjectsEnum,dependent_objects));
172  parameters->AddObject(new IntParam(AutodiffNumDependentsEnum,num_dep));
173 
174  delete dependent_objects;
175  /*}}}*/
176  /*Deal with independents: {{{*/
177 
178  /*Independents have already been recovered in iomodel->DeclareIndependents. Just do some more processing.
179  *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/
180  num_ind=iomodel->NumIndependents();
181  parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind));
182 
183  if(num_ind){
184  xp=xNew<IssmDouble>(num_ind);
185  iomodel->FillIndependents(xp);
186  parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind));
187  xDelete<IssmDouble>(xp);
188  }
189  /*}}}*/
190  }
191  #endif
192 }

◆ CreateParametersControl()

void CreateParametersControl ( Parameters parameters,
IoModel iomodel,
int  solution_type 
)

Definition at line 10 of file CreateParametersControl.cpp.

10  {
11 
12  bool control_analysis;
13  int inversiontype;
14  int nsteps;
15  int num_controls;
16  int num_costfunc;
17  char** controls = NULL;
18  int *maxiter = NULL;
19  char** cm_responses = NULL;
20  IssmDouble *cm_jump = NULL;
21  IssmDouble *optscal = NULL;
22  IssmDouble *control_scaling_factors = NULL;
23 
24  /*retrieve some parameters: */
25  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
26  iomodel->FindConstant(&inversiontype,"md.inversion.type");
27 
28  if(control_analysis){
29 
30  switch(inversiontype){
31  {
32  case 0:/*Brent Search*/
33  case 1:/*TAO*/
34  case 2:/*M1QN3*/
35  case 3:/*Validation*/
36  /*How many controls and how many responses?*/
37  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_control_parameters",InversionNumControlParametersEnum));
38  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.num_cost_functions",InversionNumCostFunctionsEnum));
39 
40  /*recover controls and convert to Enums*/
41  iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
42  if(num_controls<1) _error_("no controls found");
43  int* control_enums=xNew<int>(num_controls);
44  for(int i=0;i<num_controls;i++){
45  control_enums[i]=StringToEnumx(controls[i]);
46  xDelete<char>(controls[i]);
47  }
48  xDelete<char*>(controls);
49  parameters->AddObject(new IntVecParam(InversionControlParametersEnum,control_enums,num_controls));
50 
51  iomodel->FindConstant(&cm_responses,&num_costfunc,"md.inversion.cost_functions");
52  if(num_costfunc<1) _error_ ("no cost functions found");
53  int* costfunc_enums=xNew<int>(num_costfunc);
54  for(int i=0;i<num_costfunc;i++){
55  costfunc_enums[i]=StringToEnumx(cm_responses[i]);
56  xDelete<char>(cm_responses[i]);
57  }
58  xDelete<char*>(cm_responses);
59  parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
60 
61  xDelete<int>(control_enums);
62  xDelete<int>(costfunc_enums);
63 
64  break;
65  }
66  case 4:/*AD M1QN3*/
67  {
68  /*Intermediaries*/
69  int num_independent_objects,M;
70  char** names = NULL;
71 
72  /*this is done somewhere else*/
73  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum));
74  parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum));
75 
76  /*Step1: create controls (independents)*/
77  iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
78  _assert_(num_independent_objects>0);
79  iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names");
80  _assert_(M==num_independent_objects);
81  int* ind_enums=xNew<int>(num_independent_objects);
82  for(int i=0;i<num_independent_objects;i++){
83  ind_enums[i]=StringToEnumx(names[i]);
84  xDelete<char>(names[i]);
85  }
86 
87  parameters->AddObject(new IntVecParam(InversionControlParametersEnum,ind_enums,num_independent_objects));
88  iomodel->FindConstant(&cm_responses,&num_costfunc,"md.autodiff.dependent_object_names");
89  _assert_(num_costfunc>0);
90  if(num_costfunc<1) _error_ ("no cost functions found");
91  int* costfunc_enums=xNew<int>(num_costfunc);
92  for(int i=0;i<num_costfunc;i++){
93  costfunc_enums[i]=StringToEnumx(cm_responses[i]);
94  xDelete<char>(cm_responses[i]);
95  }
96  xDelete<char*>(cm_responses);
97  parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
98 
99  iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors");
100  parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects));
101 
102  /*cleanup*/
103  for(int i=0;i<num_independent_objects;i++){
104  xDelete<char>(names[i]);
105  }
106  xDelete<char*>(names);
107  xDelete<int>(ind_enums);
108  xDelete<int>(costfunc_enums);
109  break;
110  }
111  default:
112  _error_("not supported");
113  }
114 
115  /*Inversion type specifics*/
116  switch(inversiontype){
117  case 0:/*Brent Search*/
118  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.incomplete_adjoint",InversionIncompleteAdjointEnum));
119  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.nsteps",InversionNstepsEnum));
120  iomodel->FetchData(&cm_jump,&nsteps,NULL,"md.inversion.step_threshold");
121  iomodel->FetchData(&optscal,NULL,NULL,"md.inversion.gradient_scaling");
122  iomodel->FetchData(&maxiter,NULL,NULL,"md.inversion.maxiter_per_step");
123  parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_controls));
124  parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
125  parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
126  break;
127  case 1:/*TAO*/
128  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.incomplete_adjoint",InversionIncompleteAdjointEnum));
129  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.gatol",InversionGatolEnum));
130  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.grtol",InversionGrtolEnum));
131  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.gttol",InversionGttolEnum));
132  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum));
133  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
134  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.algorithm",InversionAlgorithmEnum));
135  break;
136  case 2:/*M1QN3*/
137  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.incomplete_adjoint",InversionIncompleteAdjointEnum));
138  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.dxmin",InversionDxminEnum));
139  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.gttol",InversionGttolEnum));
140  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum));
141  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
142  iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.inversion.control_scaling_factors");
143  parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
144  break;
145  case 3:/*Validation*/
146  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.incomplete_adjoint",InversionIncompleteAdjointEnum));
147  iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.inversion.control_scaling_factors");
148  parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_controls));
149  break;
150  case 4:/*AD M1QN3*/
151  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.dxmin",InversionDxminEnum));
152  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.gttol",InversionGttolEnum));
153  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxsteps",InversionMaxstepsEnum));
154  parameters->AddObject(iomodel->CopyConstantObject("md.inversion.maxiter",InversionMaxiterEnum));
155  break;
156  default:
157  _error_("not supported");
158  }
159 
160  xDelete<int>(maxiter);
161  xDelete<IssmDouble>(control_scaling_factors);
162  iomodel->DeleteData(cm_jump,"md.inversion.step_threshold");
163  iomodel->DeleteData(optscal,"md.inversion.gradient_scaling");
164  }
165 }

◆ CreateParametersDakota()

void CreateParametersDakota ( Parameters parameters,
IoModel iomodel,
char *  rootpath 
)

Definition at line 11 of file CreateParametersDakota.cpp.

11  {
12 
13  /*variable declarations*/
14  int i;
15  char **responsedescriptors = NULL;
16  int numresponsedescriptors;
17  char **variabledescriptors = NULL;
18  int numvariabledescriptors;
19  char *descriptor = NULL;
20  double *dakota_parameter = NULL;
21 
22  //qmu files
23  char *qmuinname = NULL;
24  char *qmuerrname = NULL;
25  char *qmuoutname = NULL;
26 
27  //descriptors:
28  char tag[50];
29 
30  bool dakota_analysis = false;
31  char *name = NULL;
32  int numberofresponses;
33  int nrows,ncols;
34 
35  //variable partitions:
36  IssmDouble **array = NULL;
37  int *mdims_array = NULL;
38  int *ndims_array = NULL;
39  int num_partitions;
40  int* intarray = NULL;
41  int M,N;
42 
43  /*recover parameters: */
44  iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
45 
46  if(dakota_analysis){
47 
48  parameters->AddObject(iomodel->CopyConstantObject("md.qmu.output",QmuOutputEnum));
49 
50  iomodel->FindConstant(&name,"md.miscellaneous.name");
51  iomodel->FindConstant(&numberofresponses,"md.qmu.numberofresponses");
52 
53  /*name of qmu input, error and output files*/
54  qmuinname=xNew<char>((strlen(rootpath)+strlen(name)+strlen(".qmu.in")+1));
55  sprintf(qmuinname,"%s%s%s",rootpath,name,".qmu.in");
56  parameters->AddObject(new StringParam(QmuInNameEnum,qmuinname));
57 
58  qmuoutname=xNew<char>((strlen(rootpath)+strlen(name)+strlen(".qmu.out")+1));
59  sprintf(qmuoutname,"%s%s%s",rootpath,name,".qmu.out");
60  parameters->AddObject(new StringParam(QmuOutNameEnum,qmuoutname));
61 
62  qmuerrname=xNew<char>((strlen(rootpath)+strlen(name)+strlen(".qmu.err")+1));
63  sprintf(qmuerrname,"%s%s%s",rootpath,name,".qmu.err");
64  parameters->AddObject(new StringParam(QmuErrNameEnum,qmuerrname));
65 
66  /*Fetch variable descriptors*/
67  iomodel->FindConstant(&variabledescriptors,&numvariabledescriptors,"md.qmu.variabledescriptors");
68 
69  /*Fetch response descriptors*/
70  iomodel->FindConstant(&responsedescriptors,&numresponsedescriptors,"md.qmu.responsedescriptors");
71 
72  /*Ok, we have all the response descriptors. Build a parameter with it: */
73  parameters->AddObject(new StringArrayParam(QmuResponsedescriptorsEnum,responsedescriptors,numresponsedescriptors));
74 
75  /*Load partitioning vectors specific to variables:*/
76  iomodel->FetchData(&array,&mdims_array,&ndims_array,&num_partitions,"md.qmu.variablepartitions");
77  parameters->AddObject(new DoubleMatArrayParam(QmuVariablePartitionsEnum,array,num_partitions,mdims_array,ndims_array));
78  iomodel->FetchData(&intarray,&M,&N,"md.qmu.variablepartitions_npart");
79  parameters->AddObject(new IntMatParam(QmuVariablePartitionsNpartEnum,intarray,M,N));
80  xDelete<int>(intarray); iomodel->FetchData(&intarray,&M,&N,"md.qmu.variablepartitions_nt");
81  parameters->AddObject(new IntMatParam(QmuVariablePartitionsNtEnum,intarray,M,N));
82 
83  /*free arrays: {{{*/
84  for(i=0;i<num_partitions;i++){
85  IssmDouble* matrix=array[i];
86  xDelete<IssmDouble>(matrix);
87  }
88  xDelete<int>(mdims_array);
89  xDelete<int>(ndims_array);
90  xDelete<IssmDouble*>(array);
91  xDelete<int>(intarray);
92  /*}}}*/
93 
94  /*Load partitioning vectors specific to responses:*/
95  iomodel->FetchData(&array,&mdims_array,&ndims_array,&num_partitions,"md.qmu.responsepartitions");
96  parameters->AddObject(new DoubleMatArrayParam(QmuResponsePartitionsEnum,array,num_partitions,mdims_array,ndims_array));
97  iomodel->FetchData(&intarray,&M,&N,"md.qmu.responsepartitions_npart");
98  parameters->AddObject(new IntMatParam(QmuResponsePartitionsNpartEnum,intarray,M,N));
99 
100  /*free arrays: {{{*/
101  for(i=0;i<num_partitions;i++){
102  IssmDouble* matrix=array[i];
103  xDelete<IssmDouble>(matrix);
104  }
105  xDelete<int>(mdims_array);
106  xDelete<int>(ndims_array);
107  xDelete<IssmDouble*>(array);
108  xDelete<int>(intarray);
109  /*}}}*/
110 
111 
112  /*Deal with data needed because of qmu variables*/
113  DataSet* dataset_variable_descriptors = new DataSet(QmuVariableDescriptorsEnum);
114  for(i=0;i<numvariabledescriptors;i++){
115  if (strncmp(variabledescriptors[i],"scaled_",7)==0){
116  /*Ok, we are dealing with a variable that is distributed over nodes. Recover the name of the variable (ex: scaled_Thickness): */
117  sscanf(variabledescriptors[i],"scaled_%s",tag);
118 
119  /*Get field name and input enum from tag*/
120  char* fieldname = NULL;
121  int param_enum = -1;
122  FieldAndEnumFromCode(&param_enum,&fieldname,tag);
123 
124  /*Recover data: */
125  iomodel->FetchData(&dakota_parameter,&nrows,&ncols,fieldname);
126  if(nrows==iomodel->numberofvertices){
127  dataset_variable_descriptors->AddObject(new DoubleMatParam(param_enum,dakota_parameter,nrows,ncols));
128  }
129  else{
130  dataset_variable_descriptors->AddObject(new DoubleTransientMatParam(param_enum,dakota_parameter,nrows,ncols));
131  }
132  xDelete<double>(dakota_parameter);
133  xDelete<char>(fieldname);
134  }
135  }
136  parameters->AddObject(new DataSetParam(QmuVariableDescriptorsEnum,dataset_variable_descriptors));
137  delete dataset_variable_descriptors;
138 
139  /*clean-up*/
140  for(i=0;i<numresponsedescriptors;i++){
141  descriptor=responsedescriptors[i];
142  xDelete<char>(descriptor);
143  }
144  xDelete<char*>(responsedescriptors);
145  for(i=0;i<numvariabledescriptors;i++){
146  descriptor=variabledescriptors[i];
147  xDelete<char>(descriptor);
148  }
149  xDelete<char*>(variabledescriptors);
150  xDelete<char>(qmuinname);
151  xDelete<char>(qmuerrname);
152  xDelete<char>(qmuoutname);
153 
154 
155 
156  }
157 
158  /*Free data*/
159  xDelete<char>(name);
160 }

◆ CreateOutputDefinitions()

void CreateOutputDefinitions ( Elements elements,
Parameters parameters,
Inputs2 inputs2,
IoModel iomodel 
)

Definition at line 10 of file CreateOutputDefinitions.cpp.

10  {
11 
12  int i,j;
13 
14  DataSet* output_definitions = NULL;
15  int* output_definition_enums = NULL;
16  int num_output_definitions;
17 
18  /*Create output_definitions dataset: */
19  output_definitions=new DataSet();
20  char** out_strings = NULL;
21  iomodel->FetchData(&out_strings,&num_output_definitions,"md.outputdefinition.list");
22  if(num_output_definitions>0){
23  output_definition_enums=xNew<int>(num_output_definitions);
24  for(int i=0;i<num_output_definitions;i++){
25  output_definition_enums[i]=StringToEnumx(out_strings[i]);
26  }
27  }
28  // free data:
29  for(int i=0;i<num_output_definitions;i++) xDelete<char>(out_strings[i]);
30  xDelete<char*>(out_strings);
31 
32  if(num_output_definitions){
33  for (i=0;i<num_output_definitions;i++){
34  if (output_definition_enums[i]==MassfluxatgateEnum){
35  /*Deal with mass flux gates:{{{ */
36 
37  /*massfluxatgate variables: */
38  int temp,numgates;
39  char **gatenames = NULL;
40  char **gatedefinitionstrings = NULL;
41  IssmDouble **gatesegments = NULL;
42  int *gatesegments_M = NULL;
43 
44  /*Fetch segments and names: */
45  iomodel->FetchMultipleData(&gatenames,&numgates, "md.massfluxatgate.name");
46  iomodel->FetchMultipleData(&gatedefinitionstrings,&temp, "md.massfluxatgate.definitionstring"); _assert_(temp==numgates);
47  iomodel->FetchMultipleData(&gatesegments,&gatesegments_M,NULL,&temp, "md.massfluxatgate.segments"); _assert_(temp==numgates);
48 
49  for(j=0;j<numgates;j++){
50  output_definitions->AddObject(new Massfluxatgate<IssmDouble>(gatenames[j],StringToEnumx(gatedefinitionstrings[j]),gatesegments_M[j],gatesegments[j]));
51  }
52  /*Free ressources:*/
53  for(j=0;j<numgates;j++){
54  char* string = gatenames[j]; xDelete<char>(string);
55  char* string2 = gatedefinitionstrings[j]; xDelete<char>(string2);
56  IssmDouble* gate = gatesegments[j]; xDelete<IssmDouble>(gate);
57  }
58  xDelete<char*>(gatenames);
59  xDelete<IssmDouble*>(gatesegments);
60  xDelete<int>(gatesegments_M);
61  xDelete<char*>(gatedefinitionstrings);
62  /*}}}*/
63  }
64  else if (output_definition_enums[i]==MisfitEnum){
65  /*Deal with misfits: {{{*/
66 
67  /*misfit variables: */
68  int nummisfits;
69  char** misfit_name_s = NULL;
70  char** misfit_definitionstring_s = NULL;
71  char** misfit_model_string_s = NULL;
72  IssmDouble** misfit_observation_s = NULL;
73  char** misfit_observation_string_s = NULL;
74  int* misfit_observation_M_s = NULL;
75  int* misfit_observation_N_s = NULL;
76  int* misfit_local_s = NULL;
77  char** misfit_timeinterpolation_s = NULL;
78  IssmDouble** misfit_weights_s = NULL;
79  int* misfit_weights_M_s = NULL;
80  int* misfit_weights_N_s = NULL;
81  char** misfit_weights_string_s = NULL;
82 
83  /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/misfit.m): */
84  iomodel->FetchMultipleData(&misfit_name_s,&nummisfits, "md.misfit.name");
85  iomodel->FetchMultipleData(&misfit_definitionstring_s,&nummisfits, "md.misfit.definitionstring");
86  iomodel->FetchMultipleData(&misfit_model_string_s,&nummisfits, "md.misfit.model_string");
87  iomodel->FetchMultipleData(&misfit_observation_s,&misfit_observation_M_s,&misfit_observation_N_s,&nummisfits, "md.misfit.observation");
88  iomodel->FetchMultipleData(&misfit_observation_string_s,&nummisfits, "md.misfit.observation_string");
89  iomodel->FetchMultipleData(&misfit_timeinterpolation_s,&nummisfits, "md.misfit.timeinterpolation");
90  iomodel->FetchMultipleData(&misfit_local_s,&nummisfits, "md.misfit.local");
91  iomodel->FetchMultipleData(&misfit_weights_s,&misfit_weights_M_s,&misfit_weights_N_s,&nummisfits, "md.misfit.weights");
92  iomodel->FetchMultipleData(&misfit_weights_string_s,&nummisfits, "md.misfit.weights_string");
93 
94  for(j=0;j<nummisfits;j++){
95 
96  int obs_vector_type=0;
97  if ((misfit_observation_M_s[j]==iomodel->numberofvertices) || (misfit_observation_M_s[j]==iomodel->numberofvertices+1)){
98  obs_vector_type=1;
99  }
100  else if ((misfit_observation_M_s[j]==iomodel->numberofelements) || (misfit_observation_M_s[j]==iomodel->numberofelements+1)){
101  obs_vector_type=2;
102  }
103  else
104  _error_("misfit observation size not supported yet");
105 
106  int weight_vector_type=0;
107  if ((misfit_weights_M_s[j]==iomodel->numberofvertices) || (misfit_weights_M_s[j]==iomodel->numberofvertices+1)){
108  weight_vector_type=1;
109  }
110  else if ((misfit_weights_M_s[j]==iomodel->numberofelements) || (misfit_weights_M_s[j]==iomodel->numberofelements+1)){
111  weight_vector_type=2;
112  }
113  else
114  _error_("misfit weight size not supported yet");
115 
116  /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
117  output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
118 
119  /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/
120  for(int k=0;k<elements->Size();k++){
121  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
122  element->InputCreate(misfit_observation_s[j],inputs2,iomodel,misfit_observation_M_s[j],misfit_observation_N_s[j],obs_vector_type,StringToEnumx(misfit_observation_string_s[j]),7);
123  element->InputCreate(misfit_weights_s[j],inputs2,iomodel,misfit_weights_M_s[j],misfit_weights_N_s[j],weight_vector_type,StringToEnumx(misfit_weights_string_s[j]),7);
124  }
125 
126  }
127 
128  /*Free ressources:*/
129  for(j=0;j<nummisfits;j++){
130  char* string=NULL;
131  IssmDouble* matrix = NULL;
132  string = misfit_definitionstring_s[j]; xDelete<char>(string);
133  string = misfit_observation_string_s[j]; xDelete<char>(string);
134  string = misfit_model_string_s[j]; xDelete<char>(string);
135  string = misfit_weights_string_s[j]; xDelete<char>(string);
136  string = misfit_name_s[j]; xDelete<char>(string);
137  string = misfit_timeinterpolation_s[j]; xDelete<char>(string);
138  matrix = misfit_observation_s[j]; xDelete<IssmDouble>(matrix);
139  matrix = misfit_weights_s[j]; xDelete<IssmDouble>(matrix);
140  }
141  xDelete<char*>(misfit_name_s);
142  xDelete<char*>(misfit_model_string_s);
143  xDelete<char*>(misfit_definitionstring_s);
144  xDelete<IssmDouble*>(misfit_observation_s);
145  xDelete<char*>(misfit_observation_string_s);
146  xDelete<int>(misfit_observation_M_s);
147  xDelete<int>(misfit_observation_N_s);
148  xDelete<int>(misfit_local_s);
149  xDelete<char*>(misfit_timeinterpolation_s);
150  xDelete<IssmDouble*>(misfit_weights_s);
151  xDelete<int>(misfit_weights_M_s);
152  xDelete<int>(misfit_weights_N_s);
153  xDelete<char*>(misfit_weights_string_s);
154  /*}}}*/
155  }
156  else if (output_definition_enums[i]==CfsurfacesquareEnum){
157  /*Deal with cfsurfacesquare: {{{*/
158 
159  /*cfsurfacesquare variables: */
160  int num_cfsurfacesquares;
161  char** cfsurfacesquare_name_s = NULL;
162  char** cfsurfacesquare_definitionstring_s = NULL;
163  char** cfsurfacesquare_model_string_s = NULL;
164  IssmDouble** cfsurfacesquare_observation_s = NULL;
165  char** cfsurfacesquare_observation_string_s = NULL;
166  int* cfsurfacesquare_observation_M_s = NULL;
167  int* cfsurfacesquare_observation_N_s = NULL;
168  IssmDouble** cfsurfacesquare_weights_s = NULL;
169  int* cfsurfacesquare_weights_M_s = NULL;
170  int* cfsurfacesquare_weights_N_s = NULL;
171  char** cfsurfacesquare_weights_string_s = NULL;
172  int* cfsurfacesquare_datatime_s = NULL;
173 
174  /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfsurfacesquare.m): */
175  iomodel->FetchMultipleData(&cfsurfacesquare_name_s,&num_cfsurfacesquares, "md.cfsurfacesquare.name");
176  iomodel->FetchMultipleData(&cfsurfacesquare_definitionstring_s,&num_cfsurfacesquares, "md.cfsurfacesquare.definitionstring");
177  iomodel->FetchMultipleData(&cfsurfacesquare_model_string_s,&num_cfsurfacesquares, "md.cfsurfacesquare.model_string");
178  iomodel->FetchMultipleData(&cfsurfacesquare_observation_s,&cfsurfacesquare_observation_M_s,&cfsurfacesquare_observation_N_s,&num_cfsurfacesquares, "md.cfsurfacesquare.observation");
179  iomodel->FetchMultipleData(&cfsurfacesquare_observation_string_s,&num_cfsurfacesquares, "md.cfsurfacesquare.observation_string");
180  iomodel->FetchMultipleData(&cfsurfacesquare_weights_s,&cfsurfacesquare_weights_M_s,&cfsurfacesquare_weights_N_s,&num_cfsurfacesquares, "md.cfsurfacesquare.weights");
181  iomodel->FetchMultipleData(&cfsurfacesquare_weights_string_s,&num_cfsurfacesquares, "md.cfsurfacesquare.weights_string");
182  iomodel->FetchMultipleData(&cfsurfacesquare_datatime_s,&num_cfsurfacesquares, "md.cfsurfacesquare.datatime");
183 
184  for(j=0;j<num_cfsurfacesquares;j++){
185 
186  int obs_vector_type=0;
187  if ((cfsurfacesquare_observation_M_s[j]==iomodel->numberofvertices) || (cfsurfacesquare_observation_M_s[j]==iomodel->numberofvertices+1)){
188  obs_vector_type=1;
189  }
190  else if ((cfsurfacesquare_observation_M_s[j]==iomodel->numberofelements) || (cfsurfacesquare_observation_M_s[j]==iomodel->numberofelements+1)){
191  obs_vector_type=2;
192  }
193  else
194  _error_("cfsurfacesquare observation size not supported yet");
195 
196  int weight_vector_type=0;
197  if ((cfsurfacesquare_weights_M_s[j]==iomodel->numberofvertices) || (cfsurfacesquare_weights_M_s[j]==iomodel->numberofvertices+1)){
198  weight_vector_type=1;
199  }
200  else if ((cfsurfacesquare_weights_M_s[j]==iomodel->numberofelements) || (cfsurfacesquare_weights_M_s[j]==iomodel->numberofelements+1)){
201  weight_vector_type=2;
202  }
203  else
204  _error_("cfsurfacesquare weight size not supported yet");
205 
206  /*First create a cfsurfacesquare object for that specific string (cfsurfacesquare_model_string_s[j]):*/
207  output_definitions->AddObject(new Cfsurfacesquare(cfsurfacesquare_name_s[j],StringToEnumx(cfsurfacesquare_definitionstring_s[j]),StringToEnumx(cfsurfacesquare_model_string_s[j]),StringToEnumx(cfsurfacesquare_observation_string_s[j]),StringToEnumx(cfsurfacesquare_weights_string_s[j]),cfsurfacesquare_datatime_s[j],false));
208 
209  /*Now, for this particular cfsurfacesquare object, make sure we plug into the elements: the observation, and the weights.*/
210  for(int k=0;k<elements->Size();k++){
211  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
212  element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j],inputs2,iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum);
213  element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j],inputs2,iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
214 
215  }
216 
217  }
218 
219  /*Free ressources:*/
220  for(j=0;j<num_cfsurfacesquares;j++){
221  char* string=NULL;
222  IssmDouble* matrix = NULL;
223 
224  string = cfsurfacesquare_definitionstring_s[j]; xDelete<char>(string);
225  string = cfsurfacesquare_observation_string_s[j]; xDelete<char>(string);
226  string = cfsurfacesquare_model_string_s[j]; xDelete<char>(string);
227  string = cfsurfacesquare_weights_string_s[j]; xDelete<char>(string);
228  string = cfsurfacesquare_name_s[j]; xDelete<char>(string);
229  matrix = cfsurfacesquare_observation_s[j]; xDelete<IssmDouble>(matrix);
230  matrix = cfsurfacesquare_weights_s[j]; xDelete<IssmDouble>(matrix);
231  }
232  xDelete<char*>(cfsurfacesquare_name_s);
233  xDelete<char*>(cfsurfacesquare_model_string_s);
234  xDelete<char*>(cfsurfacesquare_definitionstring_s);
235  xDelete<IssmDouble*>(cfsurfacesquare_observation_s);
236  xDelete<char*>(cfsurfacesquare_observation_string_s);
237  xDelete<int>(cfsurfacesquare_observation_M_s);
238  xDelete<int>(cfsurfacesquare_observation_N_s);
239  xDelete<IssmDouble*>(cfsurfacesquare_weights_s);
240  xDelete<int>(cfsurfacesquare_weights_M_s);
241  xDelete<int>(cfsurfacesquare_weights_N_s);
242  xDelete<char*>(cfsurfacesquare_weights_string_s);
243  xDelete<int>(cfsurfacesquare_datatime_s);
244  /*}}}*/
245  }
246  else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){
247  /*Deal with cfdragcoeffabsgrad: {{{*/
248 
249  /*cfdragcoeffabsgrad variables: */
250  int num_cfdragcoeffabsgrads;
251  char** cfdragcoeffabsgrad_name_s = NULL;
252  char** cfdragcoeffabsgrad_definitionstring_s = NULL;
253  IssmDouble** cfdragcoeffabsgrad_weights_s = NULL;
254  int* cfdragcoeffabsgrad_weights_M_s = NULL;
255  int* cfdragcoeffabsgrad_weights_N_s = NULL;
256  char** cfdragcoeffabsgrad_weights_string_s = NULL;
257 
258  /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfdragcoeffabsgrad.m): */
259  iomodel->FetchMultipleData(&cfdragcoeffabsgrad_name_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.name");
260  iomodel->FetchMultipleData(&cfdragcoeffabsgrad_definitionstring_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.definitionstring");
261  iomodel->FetchMultipleData(&cfdragcoeffabsgrad_weights_s,&cfdragcoeffabsgrad_weights_M_s,&cfdragcoeffabsgrad_weights_N_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.weights");
262  iomodel->FetchMultipleData(&cfdragcoeffabsgrad_weights_string_s,&num_cfdragcoeffabsgrads, "md.cfdragcoeffabsgrad.weights_string");
263 
264  for(j=0;j<num_cfdragcoeffabsgrads;j++){
265 
266  int weight_vector_type=0;
267  if ((cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofvertices) || (cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofvertices+1)){
268  weight_vector_type=1;
269  }
270  else if ((cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofelements) || (cfdragcoeffabsgrad_weights_M_s[j]==iomodel->numberofelements+1)){
271  weight_vector_type=2;
272  }
273  else
274  _error_("cfdragcoeffabsgrad weight size not supported yet");
275 
276  /*First create a cfdragcoeffabsgrad object for that specific string (cfdragcoeffabsgrad_model_string_s[j]):*/
277  output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),false));
278 
279  /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/
280  for(int k=0;k<elements->Size();k++){
281 
282  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
283 
284  element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j],inputs2,iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
285 
286  }
287 
288  }
289 
290  /*Free ressources:*/
291  for(j=0;j<num_cfdragcoeffabsgrads;j++){
292  char* string=NULL;
293  IssmDouble* matrix = NULL;
294 
295  string = cfdragcoeffabsgrad_definitionstring_s[j]; xDelete<char>(string);
296  string = cfdragcoeffabsgrad_weights_string_s[j]; xDelete<char>(string);
297  string = cfdragcoeffabsgrad_name_s[j]; xDelete<char>(string);
298  matrix = cfdragcoeffabsgrad_weights_s[j]; xDelete<IssmDouble>(matrix);
299  }
300  xDelete<char*>(cfdragcoeffabsgrad_name_s);
301  xDelete<char*>(cfdragcoeffabsgrad_definitionstring_s);
302  xDelete<IssmDouble*>(cfdragcoeffabsgrad_weights_s);
303  xDelete<int>(cfdragcoeffabsgrad_weights_M_s);
304  xDelete<int>(cfdragcoeffabsgrad_weights_N_s);
305  xDelete<char*>(cfdragcoeffabsgrad_weights_string_s);
306  /*}}}*/
307  }
308  else if (output_definition_enums[i]==CfsurfacelogvelEnum){
309  /*Deal with cfsurfacelogvel: {{{*/
310 
311  /*cfsurfacelogvel variables: */
312  int num_cfsurfacelogvels;
313  char** cfsurfacelogvel_name = NULL;
314  char** cfsurfacelogvel_definitionstring = NULL;
315  IssmDouble** cfsurfacelogvel_vxobs = NULL;
316  IssmDouble** cfsurfacelogvel_vyobs = NULL;
317  char** cfsurfacelogvel_vxobs_string = NULL;
318  char** cfsurfacelogvel_vyobs_string = NULL;
319  int* cfsurfacelogvel_observation_M = NULL;
320  int* cfsurfacelogvel_observation_N = NULL;
321  IssmDouble** cfsurfacelogvel_weights = NULL;
322  int* cfsurfacelogvel_weights_M = NULL;
323  int* cfsurfacelogvel_weights_N = NULL;
324  char** cfsurfacelogvel_weightstring = NULL;
325  int* cfsurfacelogvel_datatime = NULL;
326 
327  /*Fetch name, modeltring, observation, observationtring, etc ... (see src/m/classes/cfsurfacelogvel.m): */
328  iomodel->FetchMultipleData(&cfsurfacelogvel_name,&num_cfsurfacelogvels, "md.cfsurfacelogvel.name");
329  iomodel->FetchMultipleData(&cfsurfacelogvel_definitionstring,&num_cfsurfacelogvels, "md.cfsurfacelogvel.definitionstring");
330  iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vxobs");
331  iomodel->FetchMultipleData(&cfsurfacelogvel_vxobs_string,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vxobs_string");
332  iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs,&cfsurfacelogvel_observation_M,&cfsurfacelogvel_observation_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vyobs");
333  iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vyobs_string"); iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.weights");
334  iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels, "md.cfsurfacelogvel.weights_string");
335  iomodel->FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels, "md.cfsurfacelogvel.datatime");
336 
337  for(j=0;j<num_cfsurfacelogvels;j++){
338 
339  int obs_vector_type=0;
340  if ((cfsurfacelogvel_observation_M[j]==iomodel->numberofvertices) || (cfsurfacelogvel_observation_M[j]==iomodel->numberofvertices+1)){
341  obs_vector_type=1;
342  }
343  else if ((cfsurfacelogvel_observation_M[j]==iomodel->numberofelements) || (cfsurfacelogvel_observation_M[j]==iomodel->numberofelements+1)){
344  obs_vector_type=2;
345  }
346  else
347  _error_("cfsurfacelogvel observation size not supported yet");
348 
349  int weight_vector_type=0;
350  if ((cfsurfacelogvel_weights_M[j]==iomodel->numberofvertices) || (cfsurfacelogvel_weights_M[j]==iomodel->numberofvertices+1)){
351  weight_vector_type=1;
352  }
353  else if ((cfsurfacelogvel_weights_M[j]==iomodel->numberofelements) || (cfsurfacelogvel_weights_M[j]==iomodel->numberofelements+1)){
354  weight_vector_type=2;
355  }
356  else
357  _error_("cfsurfacelogvel weight size not supported yet");
358 
359  /*First create a cfsurfacelogvel object for that specific string (cfsurfacelogvel_modeltring[j]):*/
360  output_definitions->AddObject(new Cfsurfacelogvel(cfsurfacelogvel_name[j],StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_datatime[j],false));
361 
362  /*Now, for this particular cfsurfacelogvel object, make sure we plug into the elements: the observation, and the weights.*/
363  for(int k=0;k<elements->Size();k++){
364 
365  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
366 
367  element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j],inputs2,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum);
368  element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j],inputs2,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum);
369  element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j],inputs2,iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum);
370 
371  }
372 
373  }
374 
375  /*Free ressources:*/
376  for(j=0;j<num_cfsurfacelogvels;j++){
377  char* string=NULL;
378  IssmDouble* matrix = NULL;
379 
380  string = cfsurfacelogvel_definitionstring[j]; xDelete<char>(string);
381  string = cfsurfacelogvel_vxobs_string[j]; xDelete<char>(string);
382  string = cfsurfacelogvel_vyobs_string[j]; xDelete<char>(string);
383  string = cfsurfacelogvel_weightstring[j]; xDelete<char>(string);
384  string = cfsurfacelogvel_name[j]; xDelete<char>(string);
385  matrix = cfsurfacelogvel_weights[j]; xDelete<IssmDouble>(matrix);
386  matrix = cfsurfacelogvel_vxobs[j]; xDelete<IssmDouble>(matrix);
387  matrix = cfsurfacelogvel_vyobs[j]; xDelete<IssmDouble>(matrix);
388  }
389  xDelete<char*>(cfsurfacelogvel_name);
390  xDelete<char*>(cfsurfacelogvel_definitionstring);
391  xDelete<int>(cfsurfacelogvel_observation_M);
392  xDelete<IssmDouble*>(cfsurfacelogvel_vxobs);
393  xDelete<IssmDouble*>(cfsurfacelogvel_vyobs);
394  xDelete<char*>(cfsurfacelogvel_vxobs_string);
395  xDelete<char*>(cfsurfacelogvel_vyobs_string);
396  xDelete<int>(cfsurfacelogvel_observation_N);
397  xDelete<IssmDouble*>(cfsurfacelogvel_weights);
398  xDelete<int>(cfsurfacelogvel_weights_M);
399  xDelete<int>(cfsurfacelogvel_weights_N);
400  xDelete<char*>(cfsurfacelogvel_weightstring);
401  xDelete<int>(cfsurfacelogvel_datatime);
402  /*}}}*/
403  }
404  else if (output_definition_enums[i]==NodalvalueEnum){
405  /*Deal with nodal values: {{{*/
406 
407  /*nodal value variables: */
408  int numnodalvalues;
409  char** nodalvalue_name_s = NULL;
410  char** nodalvalue_definitionstrings = NULL;
411  char** nodalvalue_modelstrings = NULL;
412  int* nodalvalue_node_s = NULL;
413 
414  /*Fetch name, model_enum, etc ... (see src/m/classes/nodalvalue.m): */
415  iomodel->FetchMultipleData(&nodalvalue_name_s,&numnodalvalues, "md.nodalvalue.name");
416  iomodel->FetchMultipleData(&nodalvalue_definitionstrings,&numnodalvalues, "md.nodalvalue.definitionenum");
417  iomodel->FetchMultipleData(&nodalvalue_modelstrings,&numnodalvalues, "md.nodalvalue.model_enum");
418  iomodel->FetchMultipleData(&nodalvalue_node_s,&numnodalvalues, "md.nodalvalue.node");
419 
420  for(j=0;j<numnodalvalues;j++){
421 
422  /*First create a nodalvalue object for that specific enum (nodalvalue_model_enum_s[j]):*/
423  output_definitions->AddObject(new Nodalvalue(nodalvalue_name_s[j],StringToEnumx(nodalvalue_definitionstrings[j]),StringToEnumx(nodalvalue_modelstrings[j]),nodalvalue_node_s[j]-1)); //-1 because matlab to c indexing.
424  }
425 
426  /*Free ressources:*/
427  for(j=0;j<numnodalvalues;j++){
428  char* string=NULL;
429  string = nodalvalue_name_s[j]; xDelete<char>(string);
430  }
431  xDelete<char*>(nodalvalue_name_s);
432  xDelete<char*>(nodalvalue_modelstrings);
433  xDelete<char*>(nodalvalue_definitionstrings);
434  xDelete<int>(nodalvalue_node_s);
435  /*}}}*/
436  }
437  else if (output_definition_enums[i]==MassconEnum){
438  /*Deal with masscons: {{{*/
439 
440  /*masscon variables: */
441  int nummasscons;
442  char** masscon_name_s = NULL;
443  char** masscon_definitionstring_s = NULL;
444  IssmDouble** masscon_levelset_s = NULL;
445  int* masscon_levelset_M_s = NULL;
446  int* masscon_levelset_N_s = NULL;
447 
448  /*Fetch name and levelset, etc ... (see src/m/classes/masscon.m): */
449  iomodel->FetchMultipleData(&masscon_name_s,&nummasscons, "md.masscon.name");
450  iomodel->FetchMultipleData(&masscon_definitionstring_s,&nummasscons, "md.masscon.definitionstring");
451  iomodel->FetchMultipleData(&masscon_levelset_s,&masscon_levelset_M_s,&masscon_levelset_N_s,&nummasscons,"md.masscon.levelset");
452 
453  for(j=0;j<nummasscons;j++){
454 
455  /*Create a masscon object: */
456  output_definitions->AddObject(new Masscon(masscon_name_s[j],StringToEnumx(masscon_definitionstring_s[j]),masscon_levelset_s[j],masscon_levelset_M_s[j]));
457 
458  }
459 
460  /*Free ressources:*/
461  for(j=0;j<nummasscons;j++){
462  char* string=NULL;
463  IssmDouble* matrix = NULL;
464 
465  string = masscon_name_s[j]; xDelete<char>(string);
466  string = masscon_definitionstring_s[j]; xDelete<char>(string);
467  matrix = masscon_levelset_s[j]; xDelete<IssmDouble>(matrix);
468  }
469  xDelete<char*>(masscon_name_s);
470  xDelete<IssmDouble*>(masscon_levelset_s);
471  xDelete<int>(masscon_levelset_M_s);
472  xDelete<int>(masscon_levelset_N_s);
473  xDelete<char*>(masscon_definitionstring_s);
474 
475  /*}}}*/
476  }
477  else if (output_definition_enums[i]==MassconaxpbyEnum){
478  /*Deal with masscon combinations: {{{*/
479 
480  /*masscon variables: */
481  char** masscon_name_s = NULL;
482  char** masscon_definitionstring_s = NULL;
483  char** masscon_namex_s = NULL;
484  char** masscon_namey_s = NULL;
485  IssmDouble* masscon_alpha_s = NULL;
486  IssmDouble* masscon_beta_s = NULL;
487  int num;
488 
489  /*Fetch names and multiplicators, etc ... (see src/m/classes/masscon_axpby.m): */
490  iomodel->FetchMultipleData(&masscon_name_s,&num, "md.massconaxpby.name");
491  iomodel->FetchMultipleData(&masscon_definitionstring_s,&num,"md.massconaxpby.definitionstring");
492  iomodel->FetchMultipleData(&masscon_namex_s,&num, "md.massconaxpby.namex");
493  iomodel->FetchMultipleData(&masscon_namey_s,&num, "md.massconaxpby.namey");
494  iomodel->FetchMultipleData(&masscon_alpha_s,&num, "md.massconaxpby.alpha");
495  iomodel->FetchMultipleData(&masscon_beta_s,&num, "md.massconaxpby.beta");
496  for(j=0;j<num;j++){
497 
498  /*Create a masscon axpyb object: */
499  output_definitions->AddObject(new Massconaxpby(masscon_name_s[j],StringToEnumx(masscon_definitionstring_s[j]),masscon_namex_s[j],masscon_namey_s[j],masscon_alpha_s[j],masscon_beta_s[j]));
500 
501  }
502 
503  /*Free ressources:*/
504  for(j=0;j<num;j++){
505  char* string=NULL;
506  string = masscon_definitionstring_s[j]; xDelete<char>(string);
507  string = masscon_name_s[j]; xDelete<char>(string);
508  string = masscon_namex_s[j]; xDelete<char>(string);
509  string = masscon_namey_s[j]; xDelete<char>(string);
510  }
511  xDelete<char*>(masscon_definitionstring_s);
512  xDelete<char*>(masscon_name_s);
513  xDelete<char*>(masscon_namex_s);
514  xDelete<char*>(masscon_namey_s);
515  xDelete<IssmDouble>(masscon_alpha_s);
516  xDelete<IssmDouble>(masscon_beta_s);
517  /*}}}*/
518  }
519  else if (output_definition_enums[i]==RegionaloutputEnum){
520  /*Deal with regional output: {{{*/
521 
522  /*regional output variables: */
523  int numout;
524  char** reg_name_s = NULL;
525  char** reg_definitionstring_s = NULL;
526  char** reg_outputnamestring_s = NULL;
527  IssmDouble** reg_mask_s = NULL;
528  int* reg_mask_M_s = NULL;
529  int* reg_mask_N_s = NULL;
530 
531  /*Fetch name and mask, etc ... (see src/m/classes/regionaloutput.m): */
532  iomodel->FetchMultipleData(&reg_name_s,&numout, "md.regionaloutput.name");
533  iomodel->FetchMultipleData(&reg_definitionstring_s,&numout, "md.regionaloutput.definitionstring");
534  iomodel->FetchMultipleData(&reg_outputnamestring_s,&numout, "md.regionaloutput.outputnamestring");
535  iomodel->FetchMultipleData(&reg_mask_s,&reg_mask_M_s,&reg_mask_N_s,&numout, "md.regionaloutput.mask");
536  for(j=0;j<numout;j++){
537 
538  /*Create a regional output object: */
539  output_definitions->AddObject(new Regionaloutput(reg_name_s[j],StringToEnumx(reg_definitionstring_s[j]),reg_outputnamestring_s[j],reg_mask_s[j],reg_mask_M_s[j]));
540 
541  }
542 
543  /*Free ressources:*/
544  for(j=0;j<numout;j++){
545  char* string=NULL;
546  IssmDouble* matrix = NULL;
547 
548  string = reg_name_s[j]; xDelete<char>(string);
549  string = reg_definitionstring_s[j]; xDelete<char>(string);
550  string = reg_outputnamestring_s[j]; xDelete<char>(string);
551  matrix = reg_mask_s[j]; xDelete<IssmDouble>(matrix);
552  }
553  xDelete<char*>(reg_name_s);
554  xDelete<IssmDouble*>(reg_mask_s);
555  xDelete<int>(reg_mask_M_s);
556  xDelete<int>(reg_mask_N_s);
557  xDelete<char*>(reg_outputnamestring_s);
558  xDelete<char*>(reg_definitionstring_s);
559  /*}}}*/
560  }
561  else if (output_definition_enums[i]==NumberedcostfunctionEnum){
562  /*Deal with numbered cost function: {{{*/
563 
564  /*Intermediary*/
565  int numout,numout2;
566  char **ncf_name_s = NULL;
567  char **ncf_definitionstring_s = NULL;
568  char **cost_functions = NULL;
569  IssmDouble **cost_functions_weights = NULL;
570  int* cost_functions_weights_M = NULL;
571  int* cost_functions_weights_N = NULL;
572  int cost_function,domaintype;
573  int num_cost_functions;
574 
575  /*Process cost functions and convert from string to enums*/
576  iomodel->FindConstant(&num_cost_functions,"md.numberedcostfunction.num_cost_functions");
577  iomodel->FindConstant(&cost_functions,&num_cost_functions,"md.numberedcostfunction.cost_functions");
578  if(num_cost_functions<1) _error_("No cost functions found");
579  int* cost_function_enums=xNew<int>(num_cost_functions);
580  for(int i=0;i<num_cost_functions;++i){
581  cost_function_enums[i]=StringToEnumx(cost_functions[i]);
582  }
583 
584  iomodel->FetchMultipleData(&ncf_name_s,&numout,"md.numberedcostfunction.name");
585  iomodel->FetchMultipleData(&ncf_definitionstring_s,&numout2,"md.numberedcostfunction.definitionstring"); _assert_(numout2 == numout);
586  iomodel->FetchMultipleData(&cost_functions_weights,&cost_functions_weights_M,&cost_functions_weights_N,&numout2,"md.numberedcostfunction.cost_functions_coefficients"); _assert_(numout2 == numout);
587  if(numout!=1) _error_("not implemented yet, check code here");
588 
589  /*Fetch Observations */
590  iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
591  for(int i=0;i<num_cost_functions;i++){
592  cost_function=cost_function_enums[i];
593  if( cost_function==ThicknessAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.numberedcostfunction.thickness_obs",InversionThicknessObsEnum);
594  else if(cost_function==SurfaceAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.numberedcostfunction.surface_obs",InversionSurfaceObsEnum);
595  else if(cost_function==SurfaceAbsVelMisfitEnum
596  || cost_function==SurfaceRelVelMisfitEnum
597  || cost_function==SurfaceLogVelMisfitEnum
598  || cost_function==SurfaceLogVxVyMisfitEnum
599  || cost_function==SurfaceAverageVelMisfitEnum){
600  iomodel->FetchDataToInput(inputs2,elements,"md.numberedcostfunction.vx_obs",InversionVxObsEnum);
601  if(domaintype!=Domain2DverticalEnum) iomodel->FetchDataToInput(inputs2,elements,"md.numberedcostfunction.vy_obs",InversionVyObsEnum);
602  }
603  }
604 
605  for(j=0;j<numout;j++){
606 
607  /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/
608  for(int k=0;k<elements->Size();k++){
609  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
610  element->DatasetInputCreate(cost_functions_weights[j],cost_functions_weights_M[j],cost_functions_weights_N[j],cost_function_enums,num_cost_functions,inputs2,iomodel,InversionCostFunctionsCoefficientsEnum);
611  }
612  output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums));
613  }
614 
615  /*Free data: */
616  iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring");
617  xDelete<int>(cost_function_enums);
618  for(int i=0;i<num_cost_functions;i++) xDelete<char>(cost_functions[i]);
619  xDelete<char*>(cost_functions);
620 
621  /*Free ressources:*/
622  for(j=0;j<numout;j++){
623  xDelete<char>(ncf_name_s[j]);
624  xDelete<char>(ncf_definitionstring_s[j]);
625  xDelete<IssmDouble>(cost_functions_weights[j]);
626  }
627  xDelete<char*>(ncf_name_s);
628  xDelete<char*>(ncf_definitionstring_s);
629  xDelete<int>(cost_functions_weights_M);
630  xDelete<int>(cost_functions_weights_N);
631  xDelete<IssmDouble*>(cost_functions_weights);
632 
633  /*}}}*/
634  }
635  else if (output_definition_enums[i]==RadarEnum){
636  /*Deal with radar: {{{*/
637  int numout;
638  char **radar_name_s = NULL;
639  char **radar_definitionstring_s = NULL;
640  int **radar_ice_period_s = NULL;
641 
642  /*Fetch name and definition, etc ... (see src/m/classes/radar.m): */
643  iomodel->FetchMultipleData(&radar_definitionstring_s,&numout,"md.radar.definitionstring");
644  iomodel->FetchMultipleData(&radar_name_s,&numout,"md.radar.name");
645  if(numout>1) _error_("not suppored yet");
646  /*Fetch necessary inputs for calculation*/
647  //iomodel->FetchDataToInput(elements,"md.ice_period",RadarIcePeriodEnum);
648 
649  /*Add to output definitions*/
650  output_definitions->AddObject(new Radar(radar_name_s[0],StringToEnumx(radar_definitionstring_s[0])));
651  /*}}}*/
652  }
653  else _error_("output definition enum " << EnumToStringx(output_definition_enums[i]) << " not supported yet!");
654  }
655  }
656  parameters->AddObject(new DataSetParam(OutputdefinitionEnum,output_definitions));
657 
658  /*Free ressources:*/
659  delete output_definitions;
660  xDelete<int>(output_definition_enums);
661 
662 }

◆ UpdateElementsAndMaterialsControl()

void UpdateElementsAndMaterialsControl ( Elements elements,
Parameters parameters,
Inputs2 inputs2,
Materials materials,
IoModel iomodel 
)

Definition at line 11 of file UpdateElementsAndMaterialsControl.cpp.

11  {
12  /*Intermediary*/
13  bool control_analysis;
14  int M,N;
15  int control,cost_function,domaintype;
16  int num_controls,num_cost_functions;
17  IssmDouble yts,scale;
18  Element *element = NULL;
19  Material *material = NULL;
20  int *control_enums = NULL;
21  char **controls = NULL;
22  char **cost_functions = NULL;
23  IssmDouble *independent = NULL;
24  IssmDouble *independents_min = NULL;
25  IssmDouble *independents_max = NULL;
26  IssmDouble *weights = NULL;
27 
28  /*Fetch parameters: */
29  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
30  if(!control_analysis) return;
31 
32  /*Fetch parameters: */
33  bool isautodiff;
34  iomodel->FindConstant(&isautodiff,"md.autodiff.isautodiff");
35  if(isautodiff){
36  UpdateElementsAndMaterialsControlAD(elements,parameters,inputs2,materials,iomodel);
37  return;
38  }
39 
40  /*Process controls and convert from string to enums*/
41  iomodel->FindConstant(&num_controls,"md.inversion.num_control_parameters");
42  iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
43  if(num_controls<1) _error_("no controls found");
44  control_enums=xNew<int>(num_controls);
45  for(int i=0;i<num_controls;i++){
46  control_enums[i]=StringToEnumx(controls[i]);
47  }
48 
49  /*Process cost functions and convert from string to enums*/
50  iomodel->FindConstant(&num_cost_functions,"md.inversion.num_cost_functions");
51  iomodel->FindConstant(&cost_functions,&num_cost_functions,"md.inversion.cost_functions");
52  if(num_cost_functions<1) _error_("No cost functions found");
53  int* cost_function_enums=xNew<int>(num_cost_functions);
54  for(int i=0;i<num_cost_functions;++i){
55  cost_function_enums[i]=StringToEnumx(cost_functions[i]);
56  }
57 
58  /*Fetch Observations and add to inputs*/
59  iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
60  iomodel->FindConstant(&yts,"md.constants.yts");
61  iomodel->FetchData(&weights,&M,&N,"md.inversion.cost_functions_coefficients");
62 
63  /*Transpose weights for simplicity!*/
64  if(M*N && N>1){
65  IssmDouble* weights_transp = xNew<IssmDouble>(M*N);
66  for(int i=0;i<M;i++) for(int j=0;j<N;j++) weights_transp[j*M+i] = weights[i*N+j];
67  xDelete<IssmDouble>(weights);
68  weights = weights_transp;
69  }
70 
71  if(M!=iomodel->numberofvertices && N!=num_cost_functions) _error_("not supported");
72  for(int i=0;i<num_cost_functions;i++){
73  cost_function=cost_function_enums[i];
74  if( cost_function==ThicknessAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.thickness_obs",InversionThicknessObsEnum);
75  else if(cost_function==SurfaceAbsMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.surface_obs",InversionSurfaceObsEnum);
76  else if(cost_function==RheologyBInitialguessMisfitEnum) iomodel->FetchDataToInput(inputs2,elements,"md.materials.rheology_B",RheologyBInitialguessEnum);
77  else if(cost_function==SurfaceAbsVelMisfitEnum
78  || cost_function==SurfaceRelVelMisfitEnum
79  || cost_function==SurfaceLogVelMisfitEnum
80  || cost_function==SurfaceLogVxVyMisfitEnum
81  || cost_function==SurfaceAverageVelMisfitEnum){
82  iomodel->FetchDataToInput(inputs2,elements,"md.inversion.vx_obs",InversionVxObsEnum);
83  if(domaintype!=Domain2DverticalEnum) iomodel->FetchDataToInput(inputs2,elements,"md.inversion.vy_obs",InversionVyObsEnum);
84  }
85  for(int j=0;j<elements->Size();j++){
86  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
87  element->DatasetInputAdd(InversionCostFunctionsCoefficientsEnum,&weights[i*iomodel->numberofvertices],inputs2,iomodel,M,1,1,cost_function,7,cost_function);
88  }
89  }
90  parameters->AddObject(new IntParam(ControlInputSizeMEnum,iomodel->numberofvertices));
91  xDelete<IssmDouble>(weights);
92 
93  /*Get controls*/
94  iomodel->FetchData(&independents_min,&M,&N,"md.inversion.min_parameters");
95  if(M!=iomodel->numberofvertices && N!=num_controls) _error_("not supported");
96  iomodel->FetchData(&independents_max,&M,&N,"md.inversion.max_parameters");
97  if(M!=iomodel->numberofvertices && N!=num_controls) _error_("not supported");
98 
99  /*Transpose weights for simplicity!*/
100  if(M*N && N>1){
101  IssmDouble* independents_min_transp = xNew<IssmDouble>(M*N);
102  for(int i=0;i<M;i++) for(int j=0;j<N;j++) independents_min_transp[j*M+i] = independents_min[i*N+j];
103  xDelete<IssmDouble>(independents_min);
104  independents_min = independents_min_transp;
105 
106  IssmDouble* independents_max_transp = xNew<IssmDouble>(M*N);
107  for(int i=0;i<M;i++) for(int j=0;j<N;j++) independents_max_transp[j*M+i] = independents_max[i*N+j];
108  xDelete<IssmDouble>(independents_max);
109  independents_max = independents_max_transp;
110  }
111 
112  for(int i=0;i<num_controls;i++){
113  control = control_enums[i];
114  scale = 1.;
115 
116  switch(control){
117  /*List of supported controls*/
118  case BalancethicknessThickeningRateEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.thickening_rate");scale = 1./yts; break;
119  case BalancethicknessSpcthicknessEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.spcthickness"); break;
120  case VxEnum: iomodel->FetchData(&independent,&M,&N,"md.initialization.vx");scale = 1./yts; break;
121  case VyEnum: iomodel->FetchData(&independent,&M,&N,"md.initialization.vy");scale = 1./yts; break;
122  case ThicknessEnum: iomodel->FetchData(&independent,&M,&N,"md.geometry.thickness"); break;
123  case FrictionCoefficientEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.coefficient"); break;
124  case FrictionAsEnum: iomodel->FetchData(&independent,&M,&N,"md.friction.As"); break;
125  case BalancethicknessApparentMassbalanceEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.apparent_massbalance"); break;
126  case BalancethicknessOmegaEnum: iomodel->FetchData(&independent,&M,&N,"md.balancethickness.omega"); break;
127  case MaterialsRheologyBEnum: iomodel->FetchData(&independent,&M,&N,"md.materials.rheology_B"); break;
128  /*Special cases*/
129  case MaterialsRheologyBbarEnum: iomodel->FetchData(&independent,&M,&N,"md.materials.rheology_B"); break;
130  case DamageDbarEnum: iomodel->FetchData(&independent,&M,&N,"md.damage.D"); break;
131  default:
132  _error_("Control " << EnumToStringx(control) << " not implemented yet");
133  }
134  if(M!=iomodel->numberofvertices && N!=1) _error_("not supported yet");
135 
136  /*Special case if 3d*/
137  if(iomodel->domaintype==Domain3DEnum){
139  if(control==DamageDbarEnum) control=DamageDEnum;
140  }
141 
142  for(int j=0;j<elements->Size();j++){
143  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
144  element->ControlInputCreate(independent,&independents_min[i*iomodel->numberofvertices],&independents_max[i*iomodel->numberofvertices],inputs2,iomodel,M,N,scale,control,i+1);
145  }
146  xDelete<IssmDouble>(independent);
147  }
148  xDelete<IssmDouble>(independents_min);
149  xDelete<IssmDouble>(independents_max);
150 
151 
152  /*Free data: */
153  for(int i=0;i<num_controls;i++){
154  switch(control_enums[i]){
155  /*List of supported controls*/
156  case BalancethicknessThickeningRateEnum: iomodel->DeleteData(1,"md.balancethickness.thickening_rate"); break;
157  case BalancethicknessSpcthicknessEnum: iomodel->DeleteData(1,"md.balancethickness.spcthickness"); break;
158  case VxEnum: iomodel->DeleteData(1,"md.initialization.vx"); break;
159  case VyEnum: iomodel->DeleteData(1,"md.initialization.vy"); break;
160  case ThicknessEnum: iomodel->DeleteData(1,"md.geometry.thickness"); break;
161  case FrictionCoefficientEnum: iomodel->DeleteData(1,"md.friction.coefficient"); break;
162  case FrictionAsEnum: iomodel->DeleteData(1,"md.friction.As"); break;
163  case BalancethicknessApparentMassbalanceEnum: iomodel->DeleteData(1,"md.balancethickness.apparent_massbalance"); break;
164  case BalancethicknessOmegaEnum: iomodel->DeleteData(1,"md.balancethickness.omega"); break;
165  case MaterialsRheologyBEnum: iomodel->DeleteData(1,"md.materials.rheology_B"); break;
166  /*Special cases*/
167  case MaterialsRheologyBbarEnum: iomodel->DeleteData(1,"md.materials.rheology_B"); break;
168  case DamageDbarEnum: iomodel->DeleteData(1,"md.damage.D"); break;
169  default:
170  _error_("Control " << EnumToStringx(control_enums[i]) << " not implemented yet");
171  }
172  }
173 
174  xDelete<int>(control_enums);
175  xDelete<int>(cost_function_enums);
176  for(int i=0;i<num_cost_functions;i++) xDelete<char>(cost_functions[i]);
177  xDelete<char*>(cost_functions);
178  for(int i=0;i<num_controls;i++) xDelete<char>(controls[i]);
179  xDelete<char*>(controls);
180 }

◆ UpdateElementsAndMaterialsControlAD()

void UpdateElementsAndMaterialsControlAD ( Elements elements,
Parameters parameters,
Inputs2 inputs2,
Materials materials,
IoModel iomodel 
)

Definition at line 181 of file UpdateElementsAndMaterialsControl.cpp.

181  {
182 
183  #if defined(_HAVE_AD_)
184  /*Intermediaries*/
185  int num_independent_objects,M,N,M_par,N_par;
186  char** names = NULL;
187  int* types = NULL;
188  int* control_sizes = NULL;
189  int* M_all = NULL;
190  IssmDouble* independent = NULL;
191  IssmDouble* independents_fullmin = NULL;
192  IssmDouble* independents_fullmax = NULL;
193  bool control_analysis =false;
194 
195  iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
196 
197  /*Now, return if no control*/
198  if(!control_analysis) return;
199 
200  /*Step1: create controls (independents)*/
201  iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
202  _assert_(num_independent_objects>0);
203  iomodel->FetchData(&names,&M,"md.autodiff.independent_object_names");
204  _assert_(M==num_independent_objects);
205  iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
206 
207  M_all = xNew<int>(num_independent_objects);
208 
209  /*create independent objects, and at the same time, fetch the corresponding independent variables,
210  *and declare them as such in ADOLC: */
211  iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters");
212  iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
213  iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
214 
215  int* start_point = NULL;
216  start_point = xNew<int>(num_independent_objects);
217  int counter = 0;
218  for(int i=0;i<num_independent_objects;i++){
219  start_point[i]=counter;
220  counter+=control_sizes[i];
221  }
222 
223  for(int i=0;i<num_independent_objects;i++){
224 
225  if(types[i]==1){ /* vector:*/
226 
227  /*Get field name and input Enum from independent name*/
228  char* iofieldname = NULL;
229  int input_enum;
230  IssmDouble* independents_min = NULL;
231  IssmDouble* independents_max = NULL;
232 
233  FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
234 
235  /*Fetch required data*/
236  iomodel->FetchData(&independent,&M,&N,iofieldname);
237  _assert_(independent);
238  _assert_(N==control_sizes[i]);
239 
240  independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
241  independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
242  for(int m=0;m<M;m++){
243  for(int n=0;n<N;n++){
244  independents_min[N*m+n]=independents_fullmin[N_par*m+start_point[i]+n];
245  independents_max[N*m+n]=independents_fullmax[N_par*m+start_point[i]+n];
246  }
247  }
248  if(N!=1) M_all[i]=M-1;
249  else M_all[i]=M;
250 
251  for(int j=0;j<elements->Size();j++){
252  Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
253  element->ControlInputCreate(independent,independents_min,independents_max,inputs2,iomodel,M,N,1.,input_enum,i+1);
254  }
255  xDelete<IssmDouble>(independent);
256  xDelete<IssmDouble>(independents_min);
257  xDelete<IssmDouble>(independents_max);
258 
259  }
260  else{
261  _error_("Independent cannot be of size " << types[i]);
262  }
263  }
264  parameters->AddObject(new IntVecParam(ControlInputSizeNEnum,control_sizes,num_independent_objects));
265  parameters->AddObject(new IntVecParam(ControlInputSizeMEnum,M_all,num_independent_objects));
266 
267  /*cleanup*/
268  for(int i=0;i<num_independent_objects;i++){
269  xDelete<char>(names[i]);
270  }
271  xDelete<char*>(names);
272  xDelete<int>(types);
273  xDelete<int>(M_all);
274  xDelete<IssmDouble>(independents_fullmin);
275  xDelete<IssmDouble>(independents_fullmax);
276  xDelete<int>(start_point);
277  xDelete<int>(control_sizes);
278  /*Step2: create cost functions (dependents)*/
279 
280  return;
281 #else
282  _error_("AD not compiled");
283 #endif
284 }

◆ UpdateElementsAndMaterialsDakota()

void UpdateElementsAndMaterialsDakota ( Elements elements,
Inputs2 inputs2,
Materials materials,
IoModel iomodel 
)

Definition at line 11 of file UpdateElementsAndMaterialsDakota.cpp.

11  {
12 
13  /*recover parameters: */
14  bool dakota_analysis;
15  iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
16 
17  if(dakota_analysis) iomodel->FetchDataToInput(inputs2,elements,"md.geometry.hydrostatic_ratio",GeometryHydrostaticRatioEnum,0.);
18 }

◆ UpdateElementsTransient()

void UpdateElementsTransient ( Elements elements,
Parameters parameters,
Inputs2 inputs2,
IoModel iomodel 
)

Definition at line 11 of file UpdateElementsTransient.cpp.

11  {
12 
13  /*FIXME: this should go into parameterization update*/
14 
15  bool isgroundingline;
16  parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
17 
18  if(isgroundingline){
19  iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
20  }
21 }

◆ CreateNodes()

void CreateNodes ( Nodes nodes,
IoModel iomodel,
int  analysis,
int  finite_element,
bool  isamr = false,
int  approximation = NoneApproximationEnum,
int *  approximations = NULL 
)

Definition at line 36 of file CreateNodes.cpp.

36  {
37 
38  /*Intermediaries*/
39  int numnodes;
40  int element_numnodes;
41  int element_node_ids[40] = {0};
42 
43  /*Get partitioning variables*/
44  int my_rank = IssmComm::GetRank();
45  int num_procs = IssmComm::GetSize();
46  int* epart = iomodel->epart;
47 
48  /*Determine how many nodes we have in total (which depends on the type of finite element)*/
49  /*{{{*/
50  if(iomodel->meshelementtype==TriaEnum){
51  switch(finite_element){
52  case P0DGEnum:
53  numnodes = iomodel->numberofelements;
54  break;
55  case P1Enum:
56  numnodes = iomodel->numberofvertices;
57  break;
58  case P1DGEnum:
59  numnodes = 3*iomodel->numberofelements;
60  break;
62  numnodes = iomodel->numberofvertices+iomodel->numberofelements;
63  break;
64  case P2Enum:
65  EdgesPartitioning(iomodel);
66  numnodes = iomodel->numberofvertices+iomodel->numberofedges;
67  break;
68  case MINIEnum: case MINIcondensedEnum:
69  /*P1+ velocity (bubble statically condensed), P1 pressure*/
70  numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices);
71  break;
73  /*P2 velocity, P1 pressure*/
74  EdgesPartitioning(iomodel);
75  numnodes = (iomodel->numberofvertices+iomodel->numberofedges) + (iomodel->numberofvertices);
76  break;
77  case LATaylorHoodEnum:
78  /*P2 velocity*/
79  EdgesPartitioning(iomodel);
80  numnodes = (iomodel->numberofvertices+iomodel->numberofedges);
81  break;
83  /*P2b velocity, P1 DG pressure*/
84  EdgesPartitioning(iomodel);
85  numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements) + (3*iomodel->numberofelements);
86  break;
88  /*P2b velocity*/
89  EdgesPartitioning(iomodel);
90  numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements);
91  break;
92  default:
93  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
94  }
95  }
96  else if(iomodel->meshelementtype==PentaEnum){
97  switch(finite_element){
98  case P1Enum:
99  numnodes = iomodel->numberofvertices;
100  break;
102  numnodes = iomodel->numberofvertices+iomodel->numberofelements;
103  break;
104  case P1xP2Enum:
105  EdgesPartitioning(iomodel);
106  numnodes = iomodel->numberofvertices+iomodel->numberofverticaledges;
107  break;
108  case P1xP3Enum:
109  EdgesPartitioning(iomodel);
110  numnodes = iomodel->numberofvertices+2*iomodel->numberofverticaledges;
111  break;
112  case P1xP4Enum:
113  EdgesPartitioning(iomodel);
114  numnodes = iomodel->numberofvertices+3*iomodel->numberofverticaledges;
115  break;
116  case P2xP1Enum:
117  EdgesPartitioning(iomodel);
118  numnodes = iomodel->numberofvertices+iomodel->numberofhorizontaledges;
119  break;
120  case P2Enum:
121  EdgesPartitioning(iomodel);
122  FacesPartitioning(iomodel);
123  numnodes = iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces;
124  break;
125  case P2xP4Enum:
126  EdgesPartitioning(iomodel);
127  FacesPartitioning(iomodel);
128  numnodes = iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces;
129  break;
130  case MINIEnum: case MINIcondensedEnum:
131  /*P1+ velocity (bubble statically condensed), P1 pressure*/
132  numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices);
133  break;
134  case TaylorHoodEnum: case XTaylorHoodEnum:
135  /*P2 velocity, P1 pressure*/
136  EdgesPartitioning(iomodel);
137  FacesPartitioning(iomodel);
138  numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces) + (iomodel->numberofvertices);
139  break;
140  case OneLayerP4zEnum:
141  /*P2xP4 velocity, P1 pressure*/
142  EdgesPartitioning(iomodel);
143  FacesPartitioning(iomodel);
144  numnodes = (iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces) + (iomodel->numberofvertices);
145  break;
146  default:
147  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
148  }
149  }
150  else{
151  _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
152  }
153  /*}}}*/
154 
155  /*create matrix that keeps track of all ranks that have node i, and initialize as -1 (Common to all CPUs)*/
156  int* nodes_ranks = xNew<int>(MAXCONNECTIVITY*numnodes);
157  for(int i=0;i<MAXCONNECTIVITY*numnodes;i++) nodes_ranks[i] = -1;
158 
159  /*For all nodes, count how many cpus have node i (initialize with 0)*/
160  int* nodes_proc_count = xNewZeroInit<int>(numnodes);
161 
162  /*Create vector of approximation per node (used for FS: vel or pressure)*/
163  int* nodes_approx = xNew<int>(numnodes);
164  if(approximations){
165  for(int i=0;i<numnodes;i++) nodes_approx[i] = approximations[i];
166  }
167  else{
168  for(int i=0;i<numnodes;i++) nodes_approx[i] = approximation;
169  }
170 
171  /*Go through all elements and mark all vertices for all partitions*/
172  /*{{{*/
173  for(int i=0;i<iomodel->numberofelements;i++){
174 
175  /*Define nodes sids for each element*/
176  if(iomodel->meshelementtype==TriaEnum){
177  switch(finite_element){
178  case P0DGEnum:
179  element_numnodes=1;
180  element_node_ids[0]=i;
181  break;
182  case P1Enum:
183  element_numnodes=3;
184  element_node_ids[0]=iomodel->elements[3*i+0]-1;
185  element_node_ids[1]=iomodel->elements[3*i+1]-1;
186  element_node_ids[2]=iomodel->elements[3*i+2]-1;
187  break;
189  element_numnodes=4;
190  element_node_ids[0]=iomodel->elements[3*i+0]-1;
191  element_node_ids[1]=iomodel->elements[3*i+1]-1;
192  element_node_ids[2]=iomodel->elements[3*i+2]-1;
193  element_node_ids[3]=iomodel->numberofvertices+i;
194  break;
195  case P1DGEnum:
196  element_numnodes=3;
197  element_node_ids[0]=3*i+0;
198  element_node_ids[1]=3*i+1;
199  element_node_ids[2]=3*i+2;
200  break;
201  case P2Enum:
202  element_numnodes = 6;
203  element_node_ids[0]=iomodel->elements[3*i+0]-1;
204  element_node_ids[1]=iomodel->elements[3*i+1]-1;
205  element_node_ids[2]=iomodel->elements[3*i+2]-1;
206  element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
207  element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
208  element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
209  break;
210  case MINIEnum: case MINIcondensedEnum:
211  element_numnodes = 4+3;
212  element_node_ids[0]=iomodel->elements[3*i+0]-1;
213  element_node_ids[1]=iomodel->elements[3*i+1]-1;
214  element_node_ids[2]=iomodel->elements[3*i+2]-1;
215  element_node_ids[3]=iomodel->numberofvertices+i;
216  for(int n=0;n<4;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
217  element_node_ids[4]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+0]-1;
218  element_node_ids[5]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+1]-1;
219  element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+2]-1;
220  for(int n=4;n<7;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
221  break;
222  case TaylorHoodEnum: case XTaylorHoodEnum:
223  element_numnodes = 6+3;
224  element_node_ids[0]=iomodel->elements[3*i+0]-1;
225  element_node_ids[1]=iomodel->elements[3*i+1]-1;
226  element_node_ids[2]=iomodel->elements[3*i+2]-1;
227  element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
228  element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
229  element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
230  for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
231  element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+0]-1;
232  element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+1]-1;
233  element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+2]-1;
234  for(int n=6;n<9;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
235  break;
236  case LATaylorHoodEnum:
237  element_numnodes = 6;
238  element_node_ids[0]=iomodel->elements[3*i+0]-1;
239  element_node_ids[1]=iomodel->elements[3*i+1]-1;
240  element_node_ids[2]=iomodel->elements[3*i+2]-1;
241  element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
242  element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
243  element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
244  for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
245  break;
246  case CrouzeixRaviartEnum:
247  element_numnodes = 7+3;
248  element_node_ids[0]=iomodel->elements[3*i+0]-1;
249  element_node_ids[1]=iomodel->elements[3*i+1]-1;
250  element_node_ids[2]=iomodel->elements[3*i+2]-1;
251  element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
252  element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
253  element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
254  element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i;
255  for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
256  element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+0;
257  element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+1;
258  element_node_ids[9]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+2;
259  for(int n=7;n<10;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
260  break;
262  element_numnodes = 7;
263  element_node_ids[0]=iomodel->elements[3*i+0]-1;
264  element_node_ids[1]=iomodel->elements[3*i+1]-1;
265  element_node_ids[2]=iomodel->elements[3*i+2]-1;
266  element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
267  element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
268  element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
269  element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i;
270  for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
271  break;
272  default:
273  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
274  }
275  }
276  else if(iomodel->meshelementtype==PentaEnum){
277  switch(finite_element){
278  case P1Enum:
279  element_numnodes=6;
280  element_node_ids[0]=iomodel->elements[6*i+0]-1;
281  element_node_ids[1]=iomodel->elements[6*i+1]-1;
282  element_node_ids[2]=iomodel->elements[6*i+2]-1;
283  element_node_ids[3]=iomodel->elements[6*i+3]-1;
284  element_node_ids[4]=iomodel->elements[6*i+4]-1;
285  element_node_ids[5]=iomodel->elements[6*i+5]-1;
286  break;
288  element_numnodes=7;
289  element_node_ids[0]=iomodel->elements[6*i+0]-1;
290  element_node_ids[1]=iomodel->elements[6*i+1]-1;
291  element_node_ids[2]=iomodel->elements[6*i+2]-1;
292  element_node_ids[3]=iomodel->elements[6*i+3]-1;
293  element_node_ids[4]=iomodel->elements[6*i+4]-1;
294  element_node_ids[5]=iomodel->elements[6*i+5]-1;
295  element_node_ids[6]=iomodel->numberofvertices+i;
296  break;
297  case P1xP2Enum:
298  element_numnodes=9;
299  element_node_ids[0]=iomodel->elements[6*i+0]-1;
300  element_node_ids[1]=iomodel->elements[6*i+1]-1;
301  element_node_ids[2]=iomodel->elements[6*i+2]-1;
302  element_node_ids[3]=iomodel->elements[6*i+3]-1;
303  element_node_ids[4]=iomodel->elements[6*i+4]-1;
304  element_node_ids[5]=iomodel->elements[6*i+5]-1;
305  element_node_ids[6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+0];
306  element_node_ids[7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+1];
307  element_node_ids[8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+2];
308  break;
309  case P1xP3Enum:
310  element_numnodes=12;
311  element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
312  element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
313  element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
314  element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
315  element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
316  element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
317  element_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0;
318  element_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0;
319  element_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0;
320  element_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1;
321  element_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1;
322  element_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1;
323  break;
324  case P1xP4Enum:
325  element_numnodes=15;
326  element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
327  element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
328  element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
329  element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
330  element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
331  element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
332  element_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0; /*mid vertical edge 1*/
333  element_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0; /*mid vertical edge 2*/
334  element_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0; /*mid vertical edge 3*/
335  element_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 1/4 vertical edge 1*/
336  element_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 1/4 vertical edge 2*/
337  element_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 1/4 vertical edge 3*/
338  element_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+2; /* 3/4 vertical edge 1*/
339  element_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+2; /* 3/4 vertical edge 2*/
340  element_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+2; /* 3/4 vertical edge 3*/
341  break;
342  case P2xP1Enum:
343  element_numnodes=12;
344  element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
345  element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
346  element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
347  element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
348  element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
349  element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
350  element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+0];
351  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+1];
352  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+2];
353  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+3];
354  element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+4];
355  element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+5];
356  break;
357  case P2Enum:
358  element_numnodes = 18;
359  element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
360  element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
361  element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
362  element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
363  element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
364  element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
365  element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0];
366  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1];
367  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2];
368  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3];
369  element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4];
370  element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5];
371  element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6];
372  element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7];
373  element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8];
374  element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0];
375  element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1];
376  element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2];
377  break;
378  case P2xP4Enum:
379  element_numnodes = 30;
380  element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
381  element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
382  element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
383  element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
384  element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
385  element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
386  element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/
387  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/
388  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/
389  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/
390  element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/
391  element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/
392  element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/
393  element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/
394  element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/
395  element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/
396  element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/
397  element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/
398  element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/
399  element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/
400  element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/
401  element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/
402  element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/
403  element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/
404  element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/
405  element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/
406  element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/
407  element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/
408  element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/
409  element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/
410  break;
411  case MINIEnum: case MINIcondensedEnum:
412  element_numnodes = 7+6;
413  element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
414  element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
415  element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
416  element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
417  element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
418  element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
419  element_node_ids[ 6]=iomodel->numberofvertices+i;
420  if(!approximations) for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
421  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+0]-1;
422  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+1]-1;
423  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+2]-1;
424  element_node_ids[10]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+3]-1;
425  element_node_ids[11]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+4]-1;
426  element_node_ids[12]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+5]-1;
427  if(!approximations) for(int n=7;n<13;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
428  break;
429  case TaylorHoodEnum: case XTaylorHoodEnum:
430  element_numnodes = 18+6;
431  element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
432  element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
433  element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
434  element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
435  element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
436  element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
437  element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0];
438  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1];
439  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2];
440  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3];
441  element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4];
442  element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5];
443  element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6];
444  element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7];
445  element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8];
446  element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0];
447  element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1];
448  element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2];
449  for(int n=0;n<18;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
450  element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1;
451  element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1;
452  element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1;
453  element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1;
454  element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1;
455  element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1;
456  for(int n=18;n<24;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
457  break;
458  case OneLayerP4zEnum:
459  element_numnodes = 30+6;
460  element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
461  element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
462  element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
463  element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
464  element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
465  element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
466  element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/
467  element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/
468  element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/
469  element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/
470  element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/
471  element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/
472  element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/
473  element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/
474  element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/
475  element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/
476  element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/
477  element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/
478  element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/
479  element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/
480  element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/
481  element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/
482  element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/
483  element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/
484  element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/
485  element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/
486  element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/
487  element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/
488  element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/
489  element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/
490  for(int n=0;n<30;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
491  element_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1;
492  element_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1;
493  element_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1;
494  element_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1;
495  element_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1;
496  element_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1;
497  for(int n=30;n<36;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
498  break;
499  default:
500  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
501  }
502  }
503  else{
504  _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
505  }
506 
507  /*Add rank epart[i] for all nodes belonging to this element*/
508  for(int j=0;j<element_numnodes;j++){
509  int nid = element_node_ids[j]; _assert_(nid<numnodes);
510  AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[i]);
511  }
512  }
513  /*}}}*/
514  if((finite_element==P0DGEnum || finite_element==P1DGEnum) && analysis!=UzawaPressureAnalysisEnum){/*Special case for DG...{{{*/
515  if(finite_element==P1DGEnum){
516  int node_list[4];
517  if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet");
518  CreateEdges(iomodel);
519  CreateFaces(iomodel);
520  for(int i=0;i<iomodel->numberoffaces;i++){
521  int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
522  int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
523  if(e2!=-2){
524  if(epart[e1]!=epart[e2]){
525  int i1=iomodel->faces[4*i+0];
526  int i2=iomodel->faces[4*i+1];
527  int pos=-1;
528  for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j;
529  if( pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;}
530  else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;}
531  else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;}
532  else _error_("not supposed to happen");
533  pos=-1;
534  for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j;
535  if( pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;}
536  else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;}
537  else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;}
538  else _error_("not supposed to happen");
539  for(int j=0;j<4;j++){
540  int nid = node_list[j];
541  AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]);
542  AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]);
543  }
544  }
545  }
546  }
547  }
548  else if(finite_element==P0DGEnum){
549  int node_list[2];
550  if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet");
551  CreateEdges(iomodel);
552  CreateFaces(iomodel);
553  for(int i=0;i<iomodel->numberoffaces;i++){
554  int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
555  int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
556  if(e2!=-2){
557  if(epart[e1]!=epart[e2]){
558  AddNodeToRank(nodes_ranks,nodes_proc_count,e2,epart[e1]);
559  //AddNodeToRank(nodes_ranks,nodes_proc_count,e1,epart[e2]);
560  }
561  }
562  }
563  }
564  else{
565  _error_("not supported");
566  }
567  }/*}}}*/
568  /*Vertex pairing for stressbalance{{{*/
569  if(!isamr && (analysis==StressbalanceAnalysisEnum || analysis==StressbalanceVerticalAnalysisEnum)){
570  int *vertex_pairing = NULL;
571  int numvertex_pairing;
572  iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.stressbalance.vertex_pairing");
573  _assert_(numvertex_pairing==0 || finite_element==P1Enum);
574  for(int i=0;i<numvertex_pairing;i++){
575  int nid1 = vertex_pairing[2*i+0]-1;
576  int nid2 = vertex_pairing[2*i+1]-1;
577  for(int j=0;j<nodes_proc_count[nid1];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid2,nodes_ranks[MAXCONNECTIVITY*nid1+j]);
578  for(int j=0;j<nodes_proc_count[nid2];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid1,nodes_ranks[MAXCONNECTIVITY*nid2+j]);
579  }
580  xDelete<int>(vertex_pairing);
581  }
582  if(!isamr && (analysis==MasstransportAnalysisEnum || analysis==FreeSurfaceBaseAnalysisEnum || analysis==FreeSurfaceTopAnalysisEnum)){
583  int *vertex_pairing = NULL;
584  int numvertex_pairing;
585  iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.masstransport.vertex_pairing");
586  _assert_(numvertex_pairing==0 || finite_element==P1Enum);
587  for(int i=0;i<numvertex_pairing;i++){
588  int nid1 = vertex_pairing[2*i+0]-1;
589  int nid2 = vertex_pairing[2*i+1]-1;
590  for(int j=0;j<nodes_proc_count[nid1];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid2,nodes_ranks[MAXCONNECTIVITY*nid1+j]);
591  for(int j=0;j<nodes_proc_count[nid2];j++) AddNodeToRank(nodes_ranks,nodes_proc_count,nid1,nodes_ranks[MAXCONNECTIVITY*nid2+j]);
592  }
593  xDelete<int>(vertex_pairing);
594  }
595  /*}}}*/
596 
597  /*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/
598  int offset = 0;
599  int* nodes_offsets = xNew<int>(numnodes);
600  for(int i=0;i<numnodes;i++){
601  if(IsNodeInRank(nodes_ranks,nodes_proc_count,i,my_rank)){
602  nodes_offsets[i] = offset++;
603  }
604  else{
605  nodes_offsets[i] = -1;
606  }
607  }
608 
609  /*Now, Count how many clones we have with other partitions*/
610  int* common_send = xNew<int>(num_procs);
611  int* common_recv = xNew<int>(num_procs);
612  int** common_send_ids = xNew<int*>(num_procs);
613  int** common_recv_ids = xNew<int*>(num_procs);
614 
615  /*First step: allocate, Step 2: populate*/
616  for(int step=0;step<2;step++){
617  if(step==1){
618  /*Allocate send and receive arrays of ids now*/
619  for(int i=0;i<num_procs;i++){
620  _assert_(common_send[i]>=0 && common_recv[i]>=0);
621  common_send_ids[i] = xNew<int>(common_send[i]);
622  common_recv_ids[i] = xNew<int>(common_recv[i]);
623  }
624  }
625  /*Re/Initialize counters to 0*/
626  for(int i=0;i<num_procs;i++){
627  common_recv[i]=0;
628  common_send[i]=0;
629  }
630  /*Go through table and find clones/masters etc*/
631  for(int i=0;i<numnodes;i++){
632  /*If we did not find this vertex in our current partition, go to next vertex*/
633  if(nodes_offsets[i] == -1) continue;
634  /*Find in what column this rank belongs*/
635  int col = -1;
636  for(int j=0;j<MAXCONNECTIVITY;j++){
637  if(nodes_ranks[MAXCONNECTIVITY*i+j] == my_rank){
638  col = j;
639  break;
640  }
641  }
642  _assert_(col!=-1);
643 
644  /*If col==0, it is either not on boundary, or a master*/
645  if(col==0){
646  /*1. is this vertex on the boundary? Skip if not*/
647  if(nodes_ranks[MAXCONNECTIVITY*i+col+1]==-1){
648  continue;
649  }
650  else{
651  for(int j=1;j<nodes_proc_count[i];j++){
652  _assert_(nodes_ranks[MAXCONNECTIVITY*i+j]>=0);
653  int rank = nodes_ranks[MAXCONNECTIVITY*i+j];
654  if(step==1){
655  common_send_ids[rank][common_send[rank]] = nodes_offsets[i];
656  }
657  common_send[rank]++;
658  }
659  }
660  }
661  else{
662  /*3. It is a slave, record that we need to receive for this cpu*/
663  int rank = nodes_ranks[MAXCONNECTIVITY*i+0];
664  if(step==1){
665  common_recv_ids[rank][common_recv[rank]] = nodes_offsets[i];
666  }
667  common_recv[rank]++;
668  }
669  }
670  }
671  xDelete<int>(nodes_proc_count);
672 
673  /*Go ahead and create vertices now that we have all we need*/
674  for(int i=0;i<numnodes;i++){
675  if(nodes_offsets[i]!=-1){
676  bool isclone = (nodes_ranks[MAXCONNECTIVITY*i+0]!=my_rank);
677  int io_index = 0;
678  if(i<iomodel->numberofvertices) io_index = i;
679  Node* node=new Node(i+1,i,io_index,isclone,iomodel,analysis,nodes_approx[i],isamr);
680  if(finite_element==MINIcondensedEnum || finite_element==P1bubblecondensedEnum){
681  /*Bubble function is collapsed, needs to constrain it, maybe this is not the best place to do this, but that's life!*/
682  if(i>=iomodel->numberofvertices && i<iomodel->numberofvertices+iomodel->numberofelements){
683  node->HardDeactivate();
684  }
685  }
686  nodes->AddObject(node);
687  }
688  }
689 
690  /*Free data: */
691  xDelete<int>(nodes_approx);
692  xDelete<int>(nodes_ranks);
693  xDelete<int>(nodes_offsets);
694 
695  /*Assign communicators*/
696  nodes->common_send=common_send;
697  nodes->common_recv=common_recv;
698  nodes->common_send_ids=common_send_ids;
699  nodes->common_recv_ids=common_recv_ids;
700 
701  /*Finalize Initialization*/
702  nodes->Finalize();
703 }

◆ ElementsAndVerticesPartitioning()

void ElementsAndVerticesPartitioning ( IoModel iomodel)

All elements have been partitioned above, only deal with elements for this cpu:

Definition at line 17 of file ElementsAndVerticesPartitioning.cpp.

17  {
18 
19  int numberofelements2d;
20  int numberofvertices2d;
21  int numlayers;
22 
23  /*intermediary: */
24  int *epart = NULL; //element partitioning.
25  int *npart = NULL; //node partitioning.
26  int elements_width; //number of columns in elements (2d->3, 3d->6)
27  int *elements2d = NULL;
28 
29  /*Get my_rank:*/
30  int my_rank = IssmComm::GetRank();
31  int num_procs = IssmComm::GetSize();
32 
33  /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
34  if(iomodel->my_elements) return;
35 
36  /*Number of vertices per elements, needed to correctly retrieve data: */
37  /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
38  switch(iomodel->meshelementtype){
39  case TriaEnum:
40  elements_width=3;
41  numberofelements2d = 0;
42  numberofvertices2d = 0;
43  numlayers = 0;
44  break;
45  case TetraEnum:
46  elements_width=4;
47  numberofelements2d = 0;
48  numberofvertices2d = 0;
49  numlayers = 0;
50  break;
51  case PentaEnum:
52  elements_width=6;
53  iomodel->FetchData(&elements2d,NULL,NULL,"md.mesh.elements2d");
54  iomodel->FindConstant(&numberofelements2d,"md.mesh.numberofelements2d");
55  iomodel->FindConstant(&numberofvertices2d,"md.mesh.numberofvertices2d");
56  iomodel->FindConstant(&numlayers,"md.mesh.numberoflayers");
57  break;
58  default:
59  _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
60  }
61 
62  /*Use ice levelset for weights*/
63  int fordan = 0;
64  int* weights = NULL;
65  if(fordan){
66  IssmDouble* icelevelset = NULL;
67  iomodel->FetchData(&icelevelset,NULL,NULL,"md.mask.ice_levelset");
68 
69  weights = xNew<int>(iomodel->numberofvertices);
70  for(int i=0;i<iomodel->numberofvertices;i++){
71  if(icelevelset[i]>=0) weights[i] = 1;
72  if(icelevelset[i]<0) weights[i] = 100;
73  }
74  xDelete<IssmDouble>(icelevelset);
75  }
76 
77  /*Partition and free resouces*/
78  MeshPartitionx(&epart,&npart,iomodel->numberofelements,iomodel->numberofvertices,iomodel->elements,numberofelements2d,numberofvertices2d,elements2d,weights,numlayers,elements_width,iomodel->meshelementtype,num_procs);
79  xDelete<int>(elements2d);
80  xDelete<int>(npart);
81  xDelete<int>(weights);
82 
83  if(fordan){
84  for(int i=0;i<IssmComm::GetSize();i++){
85  if(i==IssmComm::GetRank()){
86  int temp =0;
87  for(int j=0;j<iomodel->numberofelements;j++) if(epart[j]==i) temp++;
88  _printf_("Partition #"<<i<<" number of elements: "<<temp<<"\n");
89  }
91  }
92  }
93 
94  /*Deal with rifts, they have to be included into one partition only, not several: */
95  int numrifts;
96  iomodel->FindConstant(&numrifts,"md.rifts.numrifts");
97  if(numrifts){
98  IssmDouble *riftinfo = NULL;
99  iomodel->FetchData(&riftinfo,&numrifts,NULL,"md.rifts.riftstruct");
100  for(int i=0;i<numrifts;i++){
101  const int RIFTINFOSIZE = 12;
102  int el1=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+2))-1; //matlab indexing to c indexing
103  int el2=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+3))-1; //matlab indexing to c indexing
104  epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding vertices;
105  }
106  iomodel->DeleteData(riftinfo,"md.rifts.riftstruct");
107  }
108 
109  /*Create my_elements, used by each partition */
110  bool *my_elements = xNewZeroInit<bool>(iomodel->numberofelements);
111 
112  /*Start figuring out, out of the partition, which elements belong to this cpu: */
113  for(int i=0;i<iomodel->numberofelements;i++){
114 
116  if(my_rank==epart[i]) my_elements[i]=true;
117  }
118 
119  /*Assign pointers to iomodel*/
120  iomodel->epart =epart;
121  iomodel->my_elements=my_elements;
122 }

◆ DiscontinuousGalerkinNodesPartitioning()

void DiscontinuousGalerkinNodesPartitioning ( bool **  pmy_nodes,
bool *  my_elements,
bool *  my_vertices,
IoModel iomodel 
)

All elements have been partitioned above, only create elements for this CPU:

Definition at line 17 of file NodesPartitioning.cpp.

17  {
18 
19  /* Each element has it own nodes (as many as vertices) + additional nodes
20  * from neighboring elements for each face. This yields to a very different
21  * partition for the nodes and the vertices. The vertices are similar to
22  * continuous galerkin, but the nodes partitioning involves faces, which
23  * messes up sorting of ids. */
24 
25  /*Intermediaries*/
26  int i,i1,i2;
27  int e1,e2;
28  int pos;
29 
30  /*Get faces and elements*/
31  CreateEdges(iomodel);
32 
33  /*Build discontinuous node partitioning
34  * - there are three nodes per element (discontinous)
35  * - for each element present of each partition, its three nodes will be in this partition
36  * - the faces require the dofs of the 2 nodes of each elements sharing the face.
37  * if the 2 elements sharing the face are on 2 different cpus, we must duplicate
38  * the two nodes that are not on the cpus so that the face can access the dofs of
39  * all its 4 nodes
40  */
41 
42  /*Allocate*/
43  bool* my_nodes=xNewZeroInit<bool>(3*iomodel->numberofelements);
44 
45  /*First: add all the nodes of all the elements belonging to this cpu*/
47  for (i=0;i<iomodel->numberofelements;i++){
48  if (my_elements[i]){
49  my_nodes[3*i+0]=true;
50  my_nodes[3*i+1]=true;
51  my_nodes[3*i+2]=true;
52  }
53  }
54  }
55  else{
56  _error_("not implemented yet");
57  }
58 
59  /*Second: add all missing nodes*/
60 
61  /*Get faces and elements*/
62  CreateFaces(iomodel);
63 
64  if(iomodel->domaintype==Domain2DhorizontalEnum){
66  for(int i=0;i<iomodel->numberoffaces;i++){
67 
68  /*Get left and right elements*/
69  e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
70  e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
71 
72  /* 1) If the element e1 is in the current partition
73  * 2) and if the face of the element is shared by another element (internal face)
74  * 3) and if this element is not in the same partition:
75  * we must clone the nodes on this partition so that the loads (Numericalflux)
76  * will have access to their properties (dofs,...)*/
77  if(my_elements[e1] && e2!=-2 && !my_elements[e2]){
78 
79  /*1: Get vertices ids*/
80  i1=iomodel->faces[4*i+0];
81  i2=iomodel->faces[4*i+1];
82 
83  /*2: Get the column where these ids are located in the index*/
84  pos=UNDEF;
85  for(int j=0;j<3;j++){
86  if(iomodel->elements[3*e2+j]==i1) pos=j;
87  }
88 
89  /*3: We have the id of the elements and the position of the vertices in the index
90  * we can now create the corresponding nodes:*/
91  if(pos==0){
92  my_nodes[e2*3+0]=true;
93  my_nodes[e2*3+2]=true;
94  }
95  else if(pos==1){
96  my_nodes[e2*3+1]=true;
97  my_nodes[e2*3+0]=true;
98  }
99  else if(pos==2){
100  my_nodes[e2*3+2]=true;
101  my_nodes[e2*3+1]=true;
102  }
103  else{
104  _error_("Problem in faces creation");
105  }
106  }
107  }
108  }
109 
110  /*Free data and assign output pointers */
111  *pmy_nodes=my_nodes;
112 }

◆ FacesPartitioning()

void FacesPartitioning ( IoModel iomodel)

Definition at line 10 of file FacesPartitioning.cpp.

10  {
11 
12  /*If faces are already present, exit*/
13  if(iomodel->my_faces) return;
14 
15  /*Get faces and elements*/
16  CreateFaces(iomodel);
18 
19  /*Mesh dependent variables*/
20  int elementnbf;
21  if(iomodel->domaintype==Domain2DhorizontalEnum){
22  elementnbf = 3;
23  }
24  else if(iomodel->domaintype==Domain2DverticalEnum){
25  elementnbf = 3;
26  }
27  else if(iomodel->domaintype==Domain3DEnum){
28  elementnbf = 5;
29  }
30  else{
31  _error_("mesh dimension not supported yet");
32  }
33  /*output: */
34  iomodel->my_faces=xNewZeroInit<bool>(iomodel->numberoffaces);
35 
36  for(int i=0;i<iomodel->numberofelements;i++){
37  if(iomodel->my_elements[i]){
38  for(int j=0;j<elementnbf;j++){
39  _assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] >= 0);
40  _assert_(iomodel->elementtofaceconnectivity[i*elementnbf+j] < iomodel->numberoffaces);
41  iomodel->my_faces[iomodel->elementtofaceconnectivity[i*elementnbf+j]] = true;
42  }
43  }
44  }
45 
46  if(iomodel->meshelementtype==PentaEnum){
47  iomodel->my_vfaces = xNewZeroInit<bool>(iomodel->numberofverticalfaces);
48  for(int i=0;i<iomodel->numberofelements;i++){
49  if(iomodel->my_elements[i]){
50  for(int j=0;j<3;j++){
52  iomodel->my_vfaces[iomodel->elementtoverticalfaceconnectivity[i*3+j]] = true;
53  }
54  }
55  }
56  }
57 }

◆ EdgesPartitioning()

void EdgesPartitioning ( IoModel iomodel)

Definition at line 10 of file EdgesPartitioning.cpp.

10  {
11 
12  /*If faces are already present, exit*/
13  if(iomodel->my_edges) return;
14 
15  /*Get edges and elements*/
16  CreateEdges(iomodel);
18 
19  /*Mesh dependent variables*/
20  int elementnbe;
21  switch(iomodel->meshelementtype){
22  case TriaEnum: elementnbe = 3; break;
23  case TetraEnum: elementnbe = 6; break;
24  case PentaEnum: elementnbe = 9; break;
25  default: _error_("mesh dimension not supported yet");
26  }
27 
28  /*output: */
29  iomodel->my_edges = xNewZeroInit<bool>(iomodel->numberofedges);
30 
31  for(int i=0;i<iomodel->numberofelements;i++){
32  if(iomodel->my_elements[i]){
33  for(int j=0;j<elementnbe;j++){
34  iomodel->my_edges[iomodel->elementtoedgeconnectivity[i*elementnbe+j]] = true;
35  }
36  }
37  }
38 
39  if(iomodel->meshelementtype==PentaEnum){
40  iomodel->my_vedges = xNewZeroInit<bool>(iomodel->numberofverticaledges);
41  iomodel->my_hedges = xNewZeroInit<bool>(iomodel->numberofhorizontaledges);
42  for(int i=0;i<iomodel->numberofelements;i++){
43  if(iomodel->my_elements[i]){
44  for(int j=0;j<3;j++) iomodel->my_vedges[iomodel->elementtoverticaledgeconnectivity[i*3+j]] = true;
45  for(int j=0;j<6;j++) iomodel->my_hedges[iomodel->elementtohorizontaledgeconnectivity[i*6+j]] = true;
46  }
47  }
48  }
49 }

◆ CreateEdges()

void CreateEdges ( IoModel iomodel)

Definition at line 9 of file CreateEdges.cpp.

9  {/*{{{*/
10 
11  /*If edges are already present, exit*/
12  if(iomodel->edges) return;
13 
14  /*Check Iomodel properties*/
15  if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
16  _assert_(iomodel->elements);
17 
18  /*Intermediaries*/
19  int i,j,v1,v2,v3;
20  int elementnbe,elementnbv;
21  int *elementedges = NULL;
22  int *elementedges_markers = NULL;
23 
24  /*Mesh dependent variables*/
25  switch(iomodel->meshelementtype){
26  case TriaEnum:
27  elementnbv = 3;
28  elementnbe = 3;
29  elementedges = xNew<int>(elementnbe*2);
30  elementedges_markers = xNew<int>(elementnbe);
31  elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = -1;
32  elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = -1;
33  elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = -1;
34  break;
35  case TetraEnum:
36  elementnbv = 4;
37  elementnbe = 6;
38  elementedges = xNew<int>(elementnbe*2);
39  elementedges_markers = xNew<int>(elementnbe);
40  elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = -1;
41  elementedges[2*1+0] = 0; elementedges[2*1+1] = 2; elementedges_markers[1] = -1;
42  elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = -1;
43  elementedges[2*3+0] = 1; elementedges[2*3+1] = 3; elementedges_markers[3] = -1;
44  elementedges[2*4+0] = 2; elementedges[2*4+1] = 3; elementedges_markers[4] = -1;
45  elementedges[2*5+0] = 0; elementedges[2*5+1] = 3; elementedges_markers[5] = -1;
46  break;
47  case PentaEnum:
48  elementnbv = 6;
49  elementnbe = 9;
50  elementedges = xNew<int>(elementnbe*2);
51  elementedges_markers = xNew<int>(elementnbe);
52  elementedges[2*0+0] = 0; elementedges[2*0+1] = 3; elementedges_markers[0] = 2;
53  elementedges[2*1+0] = 1; elementedges[2*1+1] = 4; elementedges_markers[1] = 2;
54  elementedges[2*2+0] = 2; elementedges[2*2+1] = 5; elementedges_markers[2] = 2;
55  elementedges[2*3+0] = 1; elementedges[2*3+1] = 2; elementedges_markers[3] = 1;
56  elementedges[2*4+0] = 2; elementedges[2*4+1] = 0; elementedges_markers[4] = 1;
57  elementedges[2*5+0] = 0; elementedges[2*5+1] = 1; elementedges_markers[5] = 1;
58  elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 1;
59  elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 1;
60  elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 1;
61  break;
62  default:
63  _error_("mesh dimension not supported yet");
64  }
65 
66  /*Maximum number of edges*/
67  int maxnbe = elementnbe*iomodel->numberofelements;
68 
69  /*Initialize intermediaries*/
70  int *edgestemp = xNew<int>(maxnbe*3); /*format: [vertex1 vertex2 marker]*/
71  int *vedgestemp = xNew<int>(maxnbe*2); /*format: [vertex1 vertex2] */
72  int *hedgestemp = xNew<int>(maxnbe*2); /*format: [vertex1 vertex2] */
73  int *element_edge_connectivity = xNew<int>(iomodel->numberofelements*elementnbe); /*format: [edge1 edge2 ... edgen] */
74  int *element_vedge_connectivity = NULL;
75  int *element_hedge_connectivity = NULL;
76  if(iomodel->meshelementtype==PentaEnum){
77  element_vedge_connectivity = xNew<int>(iomodel->numberofelements*3); /*format: [edge1 edge2 ... edgen] */
78  element_hedge_connectivity = xNew<int>(iomodel->numberofelements*6); /*format: [edge1 edge2 ... edgen] */
79  }
80 
81  /*Initialize chain*/
82  int* head_minv = xNew<int>(iomodel->numberofvertices);
83  int* next_edge = xNew<int>(maxnbe);
84  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
85 
86  /*Initialize number of edges*/
87  int nbe = 0;
88 
89  for(i=0;i<iomodel->numberofelements;i++){
90  for(j=0;j<elementnbe;j++){
91 
92  /*Get the two indices of the edge number j of the ith element*/
93  v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
94  v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
95 
96  /*v1 and v2 must be sorted*/
97  if(v2<v1){
98  v3=v2; v2=v1; v1=v3;
99  }
100 
101  /*This edge a priori has not been processed yet*/
102  bool exist = false;
103 
104  /*Go through all processed edges connected to v1 and check whether we have seen this edge yet*/
105  for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
106  if(edgestemp[e*3+1]==v2+1){
107  exist = true;
108  element_edge_connectivity[i*elementnbe+j]=e;
109  break;
110  }
111  }
112 
113  /*If this edge is new, add it to the lists*/
114  if(!exist){
115  _assert_(nbe<maxnbe);
116 
117  /*Update edges*/
118  edgestemp[nbe*3+0] = v1+1;
119  edgestemp[nbe*3+1] = v2+1;
120  edgestemp[nbe*3+2] = elementedges_markers[j];
121 
122  /*Update Connectivity*/
123  element_edge_connectivity[i*elementnbe+j]=nbe;
124 
125  /*Update chain*/
126  next_edge[nbe] = head_minv[v1];
127  head_minv[v1] = nbe;
128 
129  /*Increase number of edges*/
130  nbe++;
131  }
132  }
133  }
134 
135  /*vertical/horizontal edges*/
136  int nbve = 0;
137  int nbhe = 0;
138  if(iomodel->meshelementtype==PentaEnum){
139  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
140  for(i=0;i<iomodel->numberofelements;i++){
141  for(j=0;j<3;j++){
142  v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
143  v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
144  if(v2<v1){ v3=v2; v2=v1; v1=v3;}
145  bool exist = false;
146  for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
147  if(vedgestemp[e*2+1]==v2+1){
148  exist = true;
149  element_vedge_connectivity[i*3+j]=e;
150  break;
151  }
152  }
153  if(!exist){
154  vedgestemp[nbve*2+0] = v1+1;
155  vedgestemp[nbve*2+1] = v2+1;
156  element_vedge_connectivity[i*3+j]=nbve;
157  next_edge[nbve] = head_minv[v1];
158  head_minv[v1] = nbve;
159  nbve++;
160  }
161  }
162  }
163  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
164  for(i=0;i<iomodel->numberofelements;i++){
165  for(j=3;j<9;j++){
166  v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
167  v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
168  if(v2<v1){ v3=v2; v2=v1; v1=v3;}
169  bool exist = false;
170  for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
171  if(hedgestemp[e*2+1]==v2+1){
172  exist = true;
173  element_hedge_connectivity[i*6+(j-3)]=e;
174  break;
175  }
176  }
177  if(!exist){
178  hedgestemp[nbhe*2+0] = v1+1;
179  hedgestemp[nbhe*2+1] = v2+1;
180  element_hedge_connectivity[i*6+(j-3)]=nbhe;
181  next_edge[nbhe] = head_minv[v1];
182  head_minv[v1] = nbhe;
183  nbhe++;
184  }
185  }
186  }
187  }
188 
189  /*Clean up*/
190  xDelete<int>(head_minv);
191  xDelete<int>(next_edge);
192  xDelete<int>(elementedges);
193  xDelete<int>(elementedges_markers);
194 
195  /*Create final edges*/
196  int* edges = xNew<int>(nbe*3);
197  for(int i=0;i<3*nbe;i++) edges[i] = edgestemp[i];
198  xDelete<int>(edgestemp);
199  int* vedges = xNew<int>(nbve*2);
200  for(int i=0;i<2*nbve;i++) vedges[i] = vedgestemp[i];
201  xDelete<int>(vedgestemp);
202  int* hedges = xNew<int>(nbhe*2);
203  for(int i=0;i<2*nbhe;i++) hedges[i] = hedgestemp[i];
204  xDelete<int>(hedgestemp);
205 
206  /*Assign output pointers*/
207  iomodel->edges = edges;
208  iomodel->verticaledges = vedges;
209  iomodel->horizontaledges = hedges;
210  iomodel->elementtoedgeconnectivity = element_edge_connectivity;
211  iomodel->elementtoverticaledgeconnectivity = element_vedge_connectivity;
212  iomodel->elementtohorizontaledgeconnectivity = element_hedge_connectivity;
213  iomodel->numberofedges = nbe;
214  iomodel->numberofverticaledges = nbve;
215  iomodel->numberofhorizontaledges = nbhe;
216 }/*}}}*/

◆ CreateFaces()

void CreateFaces ( IoModel iomodel)

Definition at line 9 of file CreateFaces.cpp.

9  {/*{{{*/
10 
11  /*If faces are already present, exit*/
12  if(iomodel->faces) return;
13 
14  /*Some checks*/
15  if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
16  _assert_(iomodel->elements);
17 
18  /*Check Iomodel properties*/
20  /*Keep going*/
21  }
22  else if(iomodel->domaintype==Domain3DEnum){
23  CreateFaces3d(iomodel);
24  return;
25  }
26  else{
27  _error_("mesh dimension not supported yet");
28  }
29 
30  /*Intermediaries*/
31  bool exist;
32  int i,j,v1,v2,v3;
33  int maxnbf,nbf;
34 
35  /*Maximum number of faces*/
36  maxnbf = 3*iomodel->numberofelements;
37 
38  /*Initialize intermediaries*/
39  int* facestemp = xNew<int>(maxnbf*4); /*format: [vertex1 vertex2 element1 element2] */
40  bool* exchange = xNewZeroInit<bool>(maxnbf); /*Faces are ordered, we need to keep track of vertex swapping*/
41  for(i=0;i<maxnbf;i++) facestemp[i*4+3]=-1; /*Initialize last column of faces as -1 (boundary edge) */
42 
43  /*Initialize chain*/
44  int* head_minv = xNew<int>(iomodel->numberofvertices);
45  int* next_face = xNew<int>(maxnbf);
46  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
47 
48  /*Initialize number of faces*/
49  nbf = 0;
50 
51  for(i=0;i<iomodel->numberofelements;i++){
52  for(j=0;j<3;j++){
53 
54  /*Get the two indices of the edge number j of the ith triangle*/
55  v1 = iomodel->elements[i*3+j];
56  if(j==2)
57  v2 = iomodel->elements[i*3+0];
58  else
59  v2 = iomodel->elements[i*3+j+1];
60 
61  /*v1 and v2 must be sorted*/
62  if(v2<v1){
63  v3=v2; v2=v1; v1=v3;
64  }
65 
66  /*This edge a priori has not been processed yet*/
67  exist = false;
68 
69  /*Go through all processed faces connected to v1 and check whether we have seen this edge yet*/
70  _assert_(v1>=0 && v1<iomodel->numberofvertices);
71  for(int e=head_minv[v1]; e!=-1; e=next_face[e]){
72  if(facestemp[e*4+1]==v2){
73  exist = true;
74  facestemp[e*4+3]=i+1;
75  break;
76  }
77  }
78 
79  /*If this edge is new, add it to the lists*/
80  if(!exist){
81  _assert_(nbf<maxnbf);
82 
83  /*Update faces*/
84  facestemp[nbf*4+0] = v1;
85  facestemp[nbf*4+1] = v2;
86  facestemp[nbf*4+2] = i+1;
87  if(v1!=iomodel->elements[i*3+j]) exchange[nbf]=true;
88 
89  /*Update chain*/
90  next_face[nbf] = head_minv[v1];
91  head_minv[v1] = nbf;
92 
93  /*Increase number of faces*/
94  nbf++;
95  }
96  }
97  }
98 
99  /*Clean up*/
100  xDelete<int>(head_minv);
101  xDelete<int>(next_face);
102 
103  /*Create final faces*/
104  int* faces = xNew<int>(nbf*4); /*vertex1 vertex2 element1 element2*/
105  for(int i=0;i<nbf;i++){
106  if(exchange[i]){
107  faces[i*4+0]=facestemp[i*4+1];
108  faces[i*4+1]=facestemp[i*4+0];
109  }
110  else{
111  faces[i*4+0]=facestemp[i*4+0];
112  faces[i*4+1]=facestemp[i*4+1];
113  }
114  faces[i*4+2]=facestemp[i*4+2];
115  faces[i*4+3]=facestemp[i*4+3];
116  }
117  xDelete<int>(facestemp);
118  xDelete<bool>(exchange);
119 
120  /*Assign output pointers*/
121  iomodel->faces = faces;
122  iomodel->numberoffaces = nbf;
123 }/*}}}*/

◆ CreateFaces3d()

void CreateFaces3d ( IoModel iomodel)

Definition at line 124 of file CreateFaces.cpp.

124  {/*{{{*/
125 
126  /*Intermediaries*/
127  int i,j,k,v0,cols,facemaxnbv;
128  int elementnbf,elementnbv,facenbv;
129  int *elementfaces = NULL;
130  int *elementfaces_markers = NULL;
131 
132  /*Mesh specific face indexing per element*/
133  switch(iomodel->meshelementtype){
134  case PentaEnum:
135  elementnbv = 6; /*Number of vertices per element*/
136  elementnbf = 5; /*Number of faces per element*/
137  facemaxnbv = 4; /*Maximum number of vertices per face*/
138  cols = facemaxnbv + 1;
139  elementfaces = xNew<int>(elementnbf*cols);
140  elementfaces_markers = xNew<int>(elementnbf);
141  /*2 triangles*/
142  elementfaces_markers[0] = 1;
143  elementfaces_markers[1] = 1;
144  elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0; elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
145  elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 3; elementfaces[cols*1+2] = 4; elementfaces[cols*1+3] = 5;
146  /*3 quads*/
147  elementfaces_markers[2] = 2;
148  elementfaces_markers[3] = 2;
149  elementfaces_markers[4] = 2;
150  elementfaces[cols*2+0] = 4; elementfaces[cols*2+1] = 1; elementfaces[cols*2+2] = 2; elementfaces[cols*2+3] = 5; elementfaces[cols*2+4] = 4;
151  elementfaces[cols*3+0] = 4; elementfaces[cols*3+1] = 2; elementfaces[cols*3+2] = 0; elementfaces[cols*3+3] = 3; elementfaces[cols*3+4] = 5;
152  elementfaces[cols*4+0] = 4; elementfaces[cols*4+1] = 0; elementfaces[cols*4+2] = 1; elementfaces[cols*4+3] = 4; elementfaces[cols*4+4] = 3;
153  break;
154  case TetraEnum:
155  elementnbv = 4; /*Number of vertices per element*/
156  elementnbf = 4; /*Number of faces per element*/
157  facemaxnbv = 3; /*Maximum number of vertices per face*/
158  cols = facemaxnbv + 1;
159  elementfaces = xNew<int>(elementnbf*cols);
160  elementfaces_markers = xNew<int>(elementnbf);
161  /*4 triangles*/
162  elementfaces_markers[0] = 1;
163  elementfaces_markers[1] = 1;
164  elementfaces_markers[2] = 1;
165  elementfaces_markers[3] = 1;
166  elementfaces[cols*0+0] = 3; elementfaces[cols*0+1] = 0; elementfaces[cols*0+2] = 1; elementfaces[cols*0+3] = 2;
167  elementfaces[cols*1+0] = 3; elementfaces[cols*1+1] = 0; elementfaces[cols*1+2] = 3; elementfaces[cols*1+3] = 1;
168  elementfaces[cols*2+0] = 3; elementfaces[cols*2+1] = 1; elementfaces[cols*2+2] = 3; elementfaces[cols*2+3] = 2;
169  elementfaces[cols*3+0] = 3; elementfaces[cols*3+1] = 0; elementfaces[cols*3+2] = 2; elementfaces[cols*3+3] = 3;
170  break;
171  default:
172  _error_("mesh "<< EnumToStringx(iomodel->meshelementtype) <<" not supported");
173  }
174 
175  /*Allocate connectivity*/
176  int *element_face_connectivity = xNew<int>(iomodel->numberofelements*elementnbf); /*format: [face1 face2 ...] */
177  int *element_vface_connectivity = NULL;
178  if(iomodel->meshelementtype==PentaEnum){
179  element_vface_connectivity = xNew<int>(iomodel->numberofelements*3); /*format: [face1 face2 face3] */
180  }
181 
182  /*Maximum number of faces for initial allocation*/
183  int maxnbf = elementnbf*iomodel->numberofelements;
184  int facescols = 4+facemaxnbv; _assert_(facescols>6);
185 
186  /*Initialize intermediaries*/
187  int* facestemp = xNew<int>(maxnbf*facescols); /*format: [element1 element2 marker nbv vertex1 vertex2 vertex3 ...] */
188  int* vfacestemp = xNew<int>(maxnbf*4);
189  for(i=0;i<maxnbf;i++) facestemp[i*facescols+1]=-1; /*Initialize second column of faces as -1 (boundary face) */
190 
191  /*Initialize chain*/
192  int* head_minv = xNew<int>(iomodel->numberofvertices);
193  int* next_face = xNew<int>(maxnbf);
194  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
195 
196  /*Initialize number of faces and list of vertex indices*/
197  int nbf = 0;
198  int* v = xNew<int>(facemaxnbv);
199  for(i=0;i<iomodel->numberofelements;i++){
200  for(j=0;j<elementnbf;j++){
201 
202  /*Get indices of current face*/
203  facenbv = elementfaces[cols*j+0];
204  for(k=0;k<facenbv;k++){
205  v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
206  }
207 
208  /*Sort list of vertices*/
209  HeapSort(v,elementfaces[cols*j+0]);
210  v0 = v[0]; _assert_(v0>=0 && v0<iomodel->numberofvertices);
211 
212  /*This face a priori has not been processed yet*/
213  bool exist = false;
214 
215  /*Go through all processed faces connected to v0 and check whether we have seen this face yet*/
216  for(int f=head_minv[v0]; f!=-1; f=next_face[f]){
217  if(facestemp[f*facescols+5]==v[1]+1 && facestemp[f*facescols+6]==v[2]+1){
218  exist = true;
219  facestemp[f*facescols+1]=i+1;
220  element_face_connectivity[i*elementnbf+j]=f;
221  break;
222  }
223  }
224 
225  /*If this face is new, add it to the lists*/
226  if(!exist){
227  _assert_(nbf<maxnbf);
228 
229  /*Update faces*/
230  facestemp[nbf*facescols+0] = i+1;
231  facestemp[nbf*facescols+2] = elementfaces_markers[j];
232  facestemp[nbf*facescols+3] = facenbv;
233  for(k=0;k<facenbv;k++) facestemp[nbf*facescols+4+k] = v[k]+1;
234 
235  /*Update Connectivity*/
236  element_face_connectivity[i*elementnbf+j]=nbf;
237 
238  /*Update chain*/
239  next_face[nbf] = head_minv[v0];
240  head_minv[v0] = nbf;
241 
242  /*Increase number of faces*/
243  nbf++;
244  }
245  }
246  }
247 
248  /*Vertical faces*/
249  int nbvf = 0;
250  if(iomodel->meshelementtype==PentaEnum){
251  for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;
252  for(i=0;i<iomodel->numberofelements;i++){
253  for(j=2;j<5;j++){
254  for(k=0;k<4;k++) v[k] = iomodel->elements[i*elementnbv + elementfaces[cols*j+k+1]] - 1;
255  HeapSort(v,4);
256  v0 = v[0]; _assert_(v0>=0 && v0<iomodel->numberofvertices);
257  bool exist = false;
258  for(int f=head_minv[v0]; f!=-1; f=next_face[f]){
259  if(vfacestemp[f*4+1]==v[1]+1 && vfacestemp[f*4+2]==v[2]+1){
260  exist = true;
261  element_vface_connectivity[i*3+(j-2)]=f;
262  break;
263  }
264  }
265  if(!exist){ _assert_(nbvf<maxnbf);
266  for(k=0;k<4;k++) vfacestemp[nbvf*4+k] = v[k]+1;
267  element_vface_connectivity[i*3+(j-2)]=nbvf;
268  next_face[nbvf] = head_minv[v0];
269  head_minv[v0] = nbvf;
270  nbvf++;
271  }
272  }
273  }
274  }
275 
276  /*Clean up*/
277  xDelete<int>(head_minv);
278  xDelete<int>(next_face);
279  xDelete<int>(v);
280  xDelete<int>(elementfaces);
281  xDelete<int>(elementfaces_markers);
282 
283  /*Create final faces (now that we have the correct size)*/
284  int* faces = xNew<int>(nbf*facescols);
285  xMemCpy<int>(faces,facestemp,nbf*facescols);
286  xDelete<int>(facestemp);
287  int* vfaces = xNew<int>(nbvf*4);
288  xMemCpy<int>(vfaces,vfacestemp,nbvf*4);
289  xDelete<int>(vfacestemp);
290 
291  /*Assign output pointers*/
292  iomodel->faces = faces;
293  iomodel->verticalfaces = vfaces;
294  iomodel->numberoffaces = nbf;
295  iomodel->numberofverticalfaces = nbvf;
296  iomodel->facescols = facescols;
297  iomodel->elementtofaceconnectivity = element_face_connectivity;
298  iomodel->elementtoverticalfaceconnectivity = element_vface_connectivity;
299 }/*}}}*/

◆ EdgeOnBoundaryFlags()

void EdgeOnBoundaryFlags ( bool **  pflags,
IoModel iomodel 
)

Definition at line 217 of file CreateEdges.cpp.

217  {/*{{{*/
218 
219  /*Intermediaries*/
220  bool isv1,isv2;
221  int facenbv,v1,v2;
222  int id_edge,id_element;
223  int elementnbe;
224 
225  /*Mesh dependent variables*/
226  switch(iomodel->meshelementtype){
227  case TriaEnum: elementnbe = 3; break;
228  case TetraEnum: elementnbe = 6; break;
229  case PentaEnum: elementnbe = 9; break;
230  default: _error_("mesh dimension not supported yet");
231  }
232 
233  /*Get edges and allocate output*/
234  if(!iomodel->edges) CreateEdges(iomodel);
235  bool* flags = xNewZeroInit<bool>(iomodel->numberofedges);
236 
237  if(iomodel->domaindim==2){
238 
239  /*Count how many times an edge is found in elementtoedgeconnectivity*/
240  int* counter = xNewZeroInit<int>(iomodel->numberofedges);
241  for(int i=0;i<iomodel->numberofelements;i++){
242  for(int j=0;j<elementnbe;j++){
243  counter[iomodel->elementtoedgeconnectivity[elementnbe*i+j]] += 1;
244  }
245  }
246 
247  /*Now, loop over the egdes, whenever it is not connected to a second element, the edge is on boundary*/
248  for(int i=0;i<iomodel->numberofedges;i++){
249  if(counter[i]==1) flags[i]=true;
250  }
251 
252  /*Clean up*/
253  xDelete<int>(counter);
254  }
255  else if(iomodel->domaindim==3){
256 
257  /*Get faces*/
258  if(!iomodel->faces) CreateFaces(iomodel);
259 
260  /*Now, loop over the faces, whenever it is not connected to a second element, all edges are on boundary*/
261  for(int id_face=0;id_face<iomodel->numberoffaces;id_face++){
262 
263  if(iomodel->faces[id_face*iomodel->facescols+1]==-1){
264 
265  /*The face is connected to the element e only*/
266  id_element = iomodel->faces[id_face*iomodel->facescols+0]-1;
267  facenbv = iomodel->faces[id_face*iomodel->facescols+3];
268 
269  /*Get all edges for this element*/
270  for(int edge = 0; edge<elementnbe; edge++){
271 
272  id_edge = iomodel->elementtoedgeconnectivity[elementnbe*id_element+edge];
273  v1 = iomodel->edges[id_edge*3+0];
274  v2 = iomodel->edges[id_edge*3+1];
275 
276  /*Test if v1 is in the face*/
277  isv1=false;
278  for(int i=0;i<facenbv;i++){
279  if(iomodel->faces[id_face*iomodel->facescols+4+i] == v1){
280  isv1 = true; break;
281  }
282  }
283  if(!isv1) continue;
284 
285  /*test if v2 is in the face*/
286  isv2=false;
287  for(int i=0;i<facenbv;i++){
288  if(iomodel->faces[id_face*iomodel->facescols+4+i] == v2){
289  isv2 = true; break;
290  }
291  }
292 
293  /*If v1 and v2 are found, this edge is on boundary*/
294  if(isv2) flags[id_edge] = true;
295  }
296  }
297  }
298  }
299  else{
300  _error_("dimension not supported");
301  }
302 
303  /*Clean up and return*/
304  *pflags = flags;
305 }/*}}}*/

◆ CreateSingleNodeToElementConnectivity()

void CreateSingleNodeToElementConnectivity ( IoModel iomodel)

! in parallel we do not want the vertex to be connected to an element that is not in its partition!!

Definition at line 10 of file CreateSingleNodeToElementConnectivity.cpp.

10  {
11 
12  /*Intermediary*/
13  int vertexid;
14  int elementswidth;
15 
16  /*output*/
17  int* connectivity=NULL;
18 
19  /*Return if connectivity already present*/
20  if(iomodel->singlenodetoelementconnectivity) return;
21 
22  /*Some checks if debugging*/
23  _assert_(iomodel->numberofvertices);
24  _assert_(iomodel->numberofelements);
25  _assert_(iomodel->my_elements);
26  _assert_(iomodel->elements);
27 
28  /*Allocate ouput*/
29  connectivity=xNewZeroInit<int>(iomodel->numberofvertices);
30 
31  /*Get element width*/
32  switch(iomodel->meshelementtype){
33  case TriaEnum: elementswidth=3; break;
34  case PentaEnum: elementswidth=6; break;
35  case TetraEnum: elementswidth=4; break;
36  default: _error_("mesh type "<<EnumToStringx(iomodel->domaintype)<<" not supported yet");
37  }
38 
39  /*Create connectivity table*/
40  for(int i=0;i<iomodel->numberofelements;i++){
42  if(iomodel->my_elements[i]){
43  for(int j=0;j<elementswidth;j++){
44  vertexid=iomodel->elements[elementswidth*i+j];
45  _assert_(vertexid>0 && vertexid-1<iomodel->numberofvertices);
46  connectivity[vertexid-1]=i+1;
47  }
48  }
49  }
50 
51  /*Assign to iomodel*/
52  iomodel->singlenodetoelementconnectivity=connectivity;
53 }

◆ CreateNumberNodeToElementConnectivity()

void CreateNumberNodeToElementConnectivity ( IoModel iomodel)

Definition at line 16 of file CreateNumberNodeToElementConnectivity.cpp.

16  {
17 
18  /*Intermediary*/
19  int i,j;
20  int vertexid;
21  int elementswidth;
22 
23  /*output*/
24  int* connectivity=NULL;
25 
26  /*Check that this has not been done yet*/
27  if(iomodel->numbernodetoelementconnectivity) return;
28 
29  /*Some checks if debugging*/
30  _assert_(iomodel->numberofvertices);
31  _assert_(iomodel->numberofelements);
32  _assert_(iomodel->elements);
33 
34  /*Allocate ouput*/
35  connectivity=xNewZeroInit<int>(iomodel->numberofvertices);
36 
37  /*Get element width*/
38  switch(iomodel->meshelementtype){
39  case TriaEnum: elementswidth=3; break;
40  case TetraEnum: elementswidth=4; break;
41  case PentaEnum: elementswidth=6; break;
42  default: _error_("mesh not supported yet");
43  }
44 
45  /*Create connectivity table*/
46  for (i=0;i<iomodel->numberofelements;i++){
47  for (j=0;j<elementswidth;j++){
48  vertexid=iomodel->elements[elementswidth*i+j];
49  _assert_(vertexid>0 && vertexid-1<iomodel->numberofvertices);
50  connectivity[vertexid-1]+=1;
51  }
52  }
53 
54  /*Assign to iomodel*/
55  iomodel->numbernodetoelementconnectivity=connectivity;
56 }
AmrNeopzEnum
@ AmrNeopzEnum
Definition: EnumDefinitions.h:975
DataSet::Size
int Size()
Definition: DataSet.cpp:399
FreeSurfaceTopAnalysisEnum
@ FreeSurfaceTopAnalysisEnum
Definition: EnumDefinitions.h:1072
MaterialsThermalconductivityEnum
@ MaterialsThermalconductivityEnum
Definition: EnumDefinitions.h:268
InversionThicknessObsEnum
@ InversionThicknessObsEnum
Definition: EnumDefinitions.h:631
SmbAIceEnum
@ SmbAIceEnum
Definition: EnumDefinitions.h:342
Massfluxatgate
Definition: Massfluxatgate.h:18
FrictionCoefficientEnum
@ FrictionCoefficientEnum
Definition: EnumDefinitions.h:574
SurfaceObservationEnum
@ SurfaceObservationEnum
Definition: EnumDefinitions.h:827
TimesteppingFinalTimeEnum
@ TimesteppingFinalTimeEnum
Definition: EnumDefinitions.h:430
QmuResponsePartitionsNpartEnum
@ QmuResponsePartitionsNpartEnum
Definition: EnumDefinitions.h:298
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
SMBgradientselaEnum
@ SMBgradientselaEnum
Definition: EnumDefinitions.h:1249
CrouzeixRaviartEnum
@ CrouzeixRaviartEnum
Definition: EnumDefinitions.h:1023
Vertices::common_recv
int * common_recv
Definition: Vertices.h:22
Vertices
Declaration of Vertices class.
Definition: Vertices.h:15
SMBd18opddEnum
@ SMBd18opddEnum
Definition: EnumDefinitions.h:1243
TransientNumRequestedOutputsEnum
@ TransientNumRequestedOutputsEnum
Definition: EnumDefinitions.h:456
AmrGroundingLineDistanceEnum
@ AmrGroundingLineDistanceEnum
Definition: EnumDefinitions.h:20
QmuErrNameEnum
@ QmuErrNameEnum
Definition: EnumDefinitions.h:286
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IoModel::verticalfaces
int * verticalfaces
Definition: IoModel.h:89
TransientIsmasstransportEnum
@ TransientIsmasstransportEnum
Definition: EnumDefinitions.h:449
IssmDouble
double IssmDouble
Definition: types.h:37
CreateParametersDakota
void CreateParametersDakota(Parameters *parameters, IoModel *iomodel, char *rootpath)
Definition: CreateParametersDakota.cpp:11
HydrologydcEplColapseThicknessEnum
@ HydrologydcEplColapseThicknessEnum
Definition: EnumDefinitions.h:177
Nodes
Declaration of Nodes class.
Definition: Nodes.h:19
MisfitEnum
@ MisfitEnum
Definition: EnumDefinitions.h:656
AutodiffDependentObjectsEnum
@ AutodiffDependentObjectsEnum
Definition: EnumDefinitions.h:43
MaterialsRheologyEsbarEnum
@ MaterialsRheologyEsbarEnum
Definition: EnumDefinitions.h:650
TetraEnum
@ TetraEnum
Definition: EnumDefinitions.h:1300
InversionNumControlParametersEnum
@ InversionNumControlParametersEnum
Definition: EnumDefinitions.h:223
Cfdragcoeffabsgrad
Definition: Cfdragcoeffabsgrad.h:15
GroundinglineMigrationEnum
@ GroundinglineMigrationEnum
Definition: EnumDefinitions.h:161
QmuVariablePartitionsNtEnum
@ QmuVariablePartitionsNtEnum
Definition: EnumDefinitions.h:296
CreateEdges
void CreateEdges(IoModel *iomodel)
Definition: CreateEdges.cpp:9
BasalforcingsIsmip6IsLocalEnum
@ BasalforcingsIsmip6IsLocalEnum
Definition: EnumDefinitions.h:69
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
QmuOutNameEnum
@ QmuOutNameEnum
Definition: EnumDefinitions.h:289
Vertices::common_send
int * common_send
Definition: Vertices.h:24
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
TransientRequestedOutputsEnum
@ TransientRequestedOutputsEnum
Definition: EnumDefinitions.h:457
AmrThicknessErrorResolutionEnum
@ AmrThicknessErrorResolutionEnum
Definition: EnumDefinitions.h:32
InversionVyObsEnum
@ InversionVyObsEnum
Definition: EnumDefinitions.h:634
IoModel::numbernodetoelementconnectivity
int * numbernodetoelementconnectivity
Definition: IoModel.h:92
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
SMBcomponentsEnum
@ SMBcomponentsEnum
Definition: EnumDefinitions.h:1242
BasalforcingsEnum
@ BasalforcingsEnum
Definition: EnumDefinitions.h:64
BalancethicknessApparentMassbalanceEnum
@ BalancethicknessApparentMassbalanceEnum
Definition: EnumDefinitions.h:982
RegionaloutputEnum
@ RegionaloutputEnum
Definition: EnumDefinitions.h:1237
BalancethicknessThickeningRateEnum
@ BalancethicknessThickeningRateEnum
Definition: EnumDefinitions.h:474
StepEnum
@ StepEnum
Definition: EnumDefinitions.h:403
IoModel::elementtohorizontaledgeconnectivity
int * elementtohorizontaledgeconnectivity
Definition: IoModel.h:85
BasalforcingsThresholdThicknessEnum
@ BasalforcingsThresholdThicknessEnum
Definition: EnumDefinitions.h:89
HydrologydcEnum
@ HydrologydcEnum
Definition: EnumDefinitions.h:1106
DebugProfilingEnum
@ DebugProfilingEnum
Definition: EnumDefinitions.h:122
InversionNstepsEnum
@ InversionNstepsEnum
Definition: EnumDefinitions.h:222
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
Nodes::common_send
int * common_send
Definition: Nodes.h:28
MaterialsRheologyEcEnum
@ MaterialsRheologyEcEnum
Definition: EnumDefinitions.h:647
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
IoModel::elementtofaceconnectivity
int * elementtofaceconnectivity
Definition: IoModel.h:86
MathydroEnum
@ MathydroEnum
Definition: EnumDefinitions.h:1171
HydrologydcEplInitialThicknessEnum
@ HydrologydcEplInitialThicknessEnum
Definition: EnumDefinitions.h:180
MaterialsEarthDensityEnum
@ MaterialsEarthDensityEnum
Definition: EnumDefinitions.h:251
InversionControlParametersEnum
@ InversionControlParametersEnum
Definition: EnumDefinitions.h:209
Analysis::CreateConstraints
virtual void CreateConstraints(Constraints *constraints, IoModel *iomodel)=0
AutodiffXpEnum
@ AutodiffXpEnum
Definition: EnumDefinitions.h:57
IoModel::my_vfaces
bool * my_vfaces
Definition: IoModel.h:68
MeshNumberoflayersEnum
@ MeshNumberoflayersEnum
Definition: EnumDefinitions.h:272
Nodalvalue
Definition: Nodalvalue.h:17
NodalvalueEnum
@ NodalvalueEnum
Definition: EnumDefinitions.h:1199
AutodiffObufsizeEnum
@ AutodiffObufsizeEnum
Definition: EnumDefinitions.h:54
InversionCostFunctionsCoefficientsEnum
@ InversionCostFunctionsCoefficientsEnum
Definition: EnumDefinitions.h:629
Parameters
Declaration of Parameters class.
Definition: Parameters.h:18
IoModel::my_hedges
bool * my_hedges
Definition: IoModel.h:71
TransientSolutionEnum
@ TransientSolutionEnum
Definition: EnumDefinitions.h:1317
VerboseMProcessor
bool VerboseMProcessor(void)
Definition: Verbosity.cpp:22
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
CreateVertices
void CreateVertices(Elements *elements, Vertices *vertices, IoModel *iomodel, int solution_type, bool isamr)
Definition: CreateElementsVerticesAndMaterials.cpp:301
MaterialsRhoFreshwaterEnum
@ MaterialsRhoFreshwaterEnum
Definition: EnumDefinitions.h:263
MaterialsMeltingpointEnum
@ MaterialsMeltingpointEnum
Definition: EnumDefinitions.h:259
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
CfdragcoeffabsgradEnum
@ CfdragcoeffabsgradEnum
Definition: EnumDefinitions.h:1005
AutodiffFosForwardIndexEnum
@ AutodiffFosForwardIndexEnum
Definition: EnumDefinitions.h:45
BasalforcingsIsmip6TfDepthsEnum
@ BasalforcingsIsmip6TfDepthsEnum
Definition: EnumDefinitions.h:71
AdaptiveTimesteppingEnum
@ AdaptiveTimesteppingEnum
Definition: EnumDefinitions.h:969
UpdateElementsTransient
void UpdateElementsTransient(Elements *elements, Parameters *parameters, Inputs2 *inputs2, IoModel *iomodel)
Definition: UpdateElementsTransient.cpp:11
MINIEnum
@ MINIEnum
Definition: EnumDefinitions.h:1156
Constraints
Declaration of Constraints class.
Definition: Constraints.h:13
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
AmrThicknessErrorThresholdEnum
@ AmrThicknessErrorThresholdEnum
Definition: EnumDefinitions.h:33
OptionsFromAnalysis
char * OptionsFromAnalysis(char **pouttoolkit, Parameters *parameters, int analysis_type)
Definition: Parameters.cpp:595
InversionIscontrolEnum
@ InversionIscontrolEnum
Definition: EnumDefinitions.h:218
MeshAverageVertexConnectivityEnum
@ MeshAverageVertexConnectivityEnum
Definition: EnumDefinitions.h:270
MaterialsRheologyEbarEnum
@ MaterialsRheologyEbarEnum
Definition: EnumDefinitions.h:646
AmrLagEnum
@ AmrLagEnum
Definition: EnumDefinitions.h:27
BasalforcingsPicoGammaTEnum
@ BasalforcingsPicoGammaTEnum
Definition: EnumDefinitions.h:82
BasalforcingsCrustthicknessEnum
@ BasalforcingsCrustthicknessEnum
Definition: EnumDefinitions.h:60
HydrologydcIsefficientlayerEnum
@ HydrologydcIsefficientlayerEnum
Definition: EnumDefinitions.h:185
Elements
Declaration of Elements class.
Definition: Elements.h:17
AddNodeToRank
void AddNodeToRank(int *nodes_ranks, int *nodes_proc_count, int nid, int rank)
Definition: CreateNodes.cpp:25
AutodiffDriverEnum
@ AutodiffDriverEnum
Definition: EnumDefinitions.h:44
TransientIsslrEnum
@ TransientIsslrEnum
Definition: EnumDefinitions.h:452
Analysis::UpdateElements
virtual void UpdateElements(Elements *elements, Inputs2 *inputs2, IoModel *iomodel, int analysis_counter, int analysis_type)=0
InversionMaxstepsEnum
@ InversionMaxstepsEnum
Definition: EnumDefinitions.h:221
CreateParameters
void CreateParameters(Parameters *parameters, IoModel *iomodel, char *rootpath, FILE *toolkitsoptionsfid, const int solution_type)
Definition: CreateParameters.cpp:18
MaterialsLatentheatEnum
@ MaterialsLatentheatEnum
Definition: EnumDefinitions.h:254
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
BasalforcingsUppercrustheatEnum
@ BasalforcingsUppercrustheatEnum
Definition: EnumDefinitions.h:91
MismipFloatingMeltRateEnum
@ MismipFloatingMeltRateEnum
Definition: EnumDefinitions.h:1190
MeshNumberofverticesEnum
@ MeshNumberofverticesEnum
Definition: EnumDefinitions.h:273
IoModel::numberoffaces
int numberoffaces
Definition: IoModel.h:97
P1DGEnum
@ P1DGEnum
Definition: EnumDefinitions.h:1215
SmbDpermilEnum
@ SmbDpermilEnum
Definition: EnumDefinitions.h:351
SurfaceRelVelMisfitEnum
@ SurfaceRelVelMisfitEnum
Definition: EnumDefinitions.h:828
EdgesPartitioning
void EdgesPartitioning(IoModel *iomodel)
Definition: EdgesPartitioning.cpp:10
FreeSurfaceBaseAnalysisEnum
@ FreeSurfaceBaseAnalysisEnum
Definition: EnumDefinitions.h:1071
Matestar
Definition: Matestar.h:24
Material
Definition: Material.h:21
SMBhenningEnum
@ SMBhenningEnum
Definition: EnumDefinitions.h:1250
AnalysisCounterEnum
@ AnalysisCounterEnum
Definition: EnumDefinitions.h:35
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
AmrGroundingLineResolutionEnum
@ AmrGroundingLineResolutionEnum
Definition: EnumDefinitions.h:21
BasalforcingsPicoEnum
@ BasalforcingsPicoEnum
Definition: EnumDefinitions.h:990
HydrologydcSedimentPorosityEnum
@ HydrologydcSedimentPorosityEnum
Definition: EnumDefinitions.h:194
P2xP1Enum
@ P2xP1Enum
Definition: EnumDefinitions.h:1226
TransientArrayParam
Definition: TransientArrayParam.h:20
Nodes::common_recv_ids
int ** common_recv_ids
Definition: Nodes.h:27
SettingsSbCouplingFrequencyEnum
@ SettingsSbCouplingFrequencyEnum
Definition: EnumDefinitions.h:339
QmuVariableDescriptorsEnum
@ QmuVariableDescriptorsEnum
Definition: EnumDefinitions.h:293
Radar
Definition: Radar.h:15
ControlInputSizeMEnum
@ ControlInputSizeMEnum
Definition: EnumDefinitions.h:105
AmrRestartEnum
@ AmrRestartEnum
Definition: EnumDefinitions.h:29
P1bubblecondensedEnum
@ P1bubblecondensedEnum
Definition: EnumDefinitions.h:1219
IoModel::FetchMultipleData
void FetchMultipleData(char ***pstringarray, int *pnumstrings, const char *data_name)
Definition: IoModel.cpp:1960
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
SteadystateNumRequestedOutputsEnum
@ SteadystateNumRequestedOutputsEnum
Definition: EnumDefinitions.h:400
MaterialsRheologyNEnum
@ MaterialsRheologyNEnum
Definition: EnumDefinitions.h:651
NumberedcostfunctionEnum
@ NumberedcostfunctionEnum
Definition: EnumDefinitions.h:1203
AmrLevelMaxEnum
@ AmrLevelMaxEnum
Definition: EnumDefinitions.h:28
AmrHmaxEnum
@ AmrHmaxEnum
Definition: EnumDefinitions.h:22
MassfluxatgateEnum
@ MassfluxatgateEnum
Definition: EnumDefinitions.h:1162
CfsurfacelogvelEnum
@ CfsurfacelogvelEnum
Definition: EnumDefinitions.h:1006
FSpressureEnum
@ FSpressureEnum
Definition: EnumDefinitions.h:1062
UpdateElementsAndMaterialsControl
void UpdateElementsAndMaterialsControl(Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
Definition: UpdateElementsAndMaterialsControl.cpp:11
AutodiffFovForwardIndicesEnum
@ AutodiffFovForwardIndicesEnum
Definition: EnumDefinitions.h:47
BasalforcingsUpperwaterElevationEnum
@ BasalforcingsUpperwaterElevationEnum
Definition: EnumDefinitions.h:94
IntVecParam
Definition: IntVecParam.h:20
StressbalanceRiftPenaltyThresholdEnum
@ StressbalanceRiftPenaltyThresholdEnum
Definition: EnumDefinitions.h:413
SettingsResultsOnNodesEnum
@ SettingsResultsOnNodesEnum
Definition: EnumDefinitions.h:338
TransientIsgroundinglineEnum
@ TransientIsgroundinglineEnum
Definition: EnumDefinitions.h:447
BasalforcingsIsmip6NumBasinsEnum
@ BasalforcingsIsmip6NumBasinsEnum
Definition: EnumDefinitions.h:70
BasalforcingsTopplumedepthEnum
@ BasalforcingsTopplumedepthEnum
Definition: EnumDefinitions.h:90
SealevelriseRunCountEnum
@ SealevelriseRunCountEnum
Definition: EnumDefinitions.h:331
AmrDeviatoricErrorThresholdEnum
@ AmrDeviatoricErrorThresholdEnum
Definition: EnumDefinitions.h:16
SteadystateRequestedOutputsEnum
@ SteadystateRequestedOutputsEnum
Definition: EnumDefinitions.h:402
IsVertexInRank
bool IsVertexInRank(int *vertices_ranks, int *vertices_proc_count, int vid, int rank)
Definition: CreateElementsVerticesAndMaterials.cpp:12
DamageEvolutionNumRequestedOutputsEnum
@ DamageEvolutionNumRequestedOutputsEnum
Definition: EnumDefinitions.h:113
OutputdefinitionEnum
@ OutputdefinitionEnum
Definition: EnumDefinitions.h:285
SettingsSolverResidueThresholdEnum
@ SettingsSolverResidueThresholdEnum
Definition: EnumDefinitions.h:340
Vertices::common_recv_ids
int ** common_recv_ids
Definition: Vertices.h:23
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
ConstantsReferencetemperatureEnum
@ ConstantsReferencetemperatureEnum
Definition: EnumDefinitions.h:103
Penta
Definition: Penta.h:29
BasalforcingsDeepwaterMeltingRateEnum
@ BasalforcingsDeepwaterMeltingRateEnum
Definition: EnumDefinitions.h:62
SMBgradientsEnum
@ SMBgradientsEnum
Definition: EnumDefinitions.h:1247
InversionMaxiterPerStepEnum
@ InversionMaxiterPerStepEnum
Definition: EnumDefinitions.h:220
QmuOutputEnum
@ QmuOutputEnum
Definition: EnumDefinitions.h:290
BasalforcingsLowercrustheatEnum
@ BasalforcingsLowercrustheatEnum
Definition: EnumDefinitions.h:72
QmuResponsePartitionsEnum
@ QmuResponsePartitionsEnum
Definition: EnumDefinitions.h:297
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
UpdateElementsAndMaterialsControlAD
void UpdateElementsAndMaterialsControlAD(Elements *elements, Parameters *parameters, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
Definition: UpdateElementsAndMaterialsControl.cpp:181
CalvingLawEnum
@ CalvingLawEnum
Definition: EnumDefinitions.h:98
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
StressbalanceMaxiterEnum
@ StressbalanceMaxiterEnum
Definition: EnumDefinitions.h:407
FrictionAsEnum
@ FrictionAsEnum
Definition: EnumDefinitions.h:571
SurfaceLogVxVyMisfitEnum
@ SurfaceLogVxVyMisfitEnum
Definition: EnumDefinitions.h:826
ElementHook::hmaterial
Hook * hmaterial
Definition: ElementHook.h:17
FieldAndEnumFromCode
void FieldAndEnumFromCode(int *out_enum, char **pfield, const char *string_in)
Definition: IoCodeConversions.cpp:9
MaterialsRheologyBbarEnum
@ MaterialsRheologyBbarEnum
Definition: EnumDefinitions.h:644
QmuInNameEnum
@ QmuInNameEnum
Definition: EnumDefinitions.h:287
SurfaceAbsVelMisfitEnum
@ SurfaceAbsVelMisfitEnum
Definition: EnumDefinitions.h:818
DefaultAnalysisEnum
@ DefaultAnalysisEnum
Definition: EnumDefinitions.h:1032
AmrHminEnum
@ AmrHminEnum
Definition: EnumDefinitions.h:23
P0DGEnum
@ P0DGEnum
Definition: EnumDefinitions.h:1214
Element
Definition: Element.h:41
IssmSolverTypeFromToolkitOptions
int IssmSolverTypeFromToolkitOptions(void)
Definition: IssmToolkitUtils.cpp:86
Numberedcostfunction
Definition: Numberedcostfunction.h:15
BasalforcingsPicoNumBasinsEnum
@ BasalforcingsPicoNumBasinsEnum
Definition: EnumDefinitions.h:85
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
RheologyBInitialguessEnum
@ RheologyBInitialguessEnum
Definition: EnumDefinitions.h:672
MassconEnum
@ MassconEnum
Definition: EnumDefinitions.h:1160
TaylorHoodEnum
@ TaylorHoodEnum
Definition: EnumDefinitions.h:1299
InversionIncompleteAdjointEnum
@ InversionIncompleteAdjointEnum
Definition: EnumDefinitions.h:217
Element::InputCreate
void InputCreate(IssmDouble *vector, Inputs2 *inputs2, IoModel *iomodel, int M, int N, int vector_type, int vector_enum, int code)
Definition: Element.cpp:1626
SMBsemicEnum
@ SMBsemicEnum
Definition: EnumDefinitions.h:1254
HydrologydcSedimentCompressibilityEnum
@ HydrologydcSedimentCompressibilityEnum
Definition: EnumDefinitions.h:191
UzawaPressureAnalysisEnum
@ UzawaPressureAnalysisEnum
Definition: EnumDefinitions.h:1320
LATaylorHoodEnum
@ LATaylorHoodEnum
Definition: EnumDefinitions.h:1139
GenericParam::GetParameterValue
P & GetParameterValue()
Definition: GenericParam.h:62
IoModel::DeleteData
void DeleteData(int num,...)
Definition: IoModel.cpp:500
BalancethicknessOmegaEnum
@ BalancethicknessOmegaEnum
Definition: EnumDefinitions.h:473
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
IoModel::facescols
int facescols
Definition: IoModel.h:90
DependentObject
Definition: DependentObject.h:15
SurfaceAbsMisfitEnum
@ SurfaceAbsMisfitEnum
Definition: EnumDefinitions.h:817
DoubleParam
Definition: DoubleParam.h:20
Tetra
Definition: Tetra.h:27
Nodes::common_send_ids
int ** common_send_ids
Definition: Nodes.h:29
BasalforcingsNusseltEnum
@ BasalforcingsNusseltEnum
Definition: EnumDefinitions.h:75
IntMatParam
Definition: IntMatParam.h:20
Analysis::CreateLoads
virtual void CreateLoads(Loads *loads, IoModel *iomodel)=0
IoModel::CopyConstantObject
Param * CopyConstantObject(const char *constant_name, int param_enum)
Definition: IoModel.cpp:418
IoCodeToEnumNature
int IoCodeToEnumNature(int enum_in)
Definition: IoCodeConversions.cpp:270
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
IoModel::elementtoverticaledgeconnectivity
int * elementtoverticaledgeconnectivity
Definition: IoModel.h:84
AutodiffLbufsizeEnum
@ AutodiffLbufsizeEnum
Definition: EnumDefinitions.h:51
Materials
Declaration of Materials class.
Definition: Materials.h:16
MumpsEnum
@ MumpsEnum
Definition: EnumDefinitions.h:1195
BasalforcingsPlumeyEnum
@ BasalforcingsPlumeyEnum
Definition: EnumDefinitions.h:88
VyObsEnum
@ VyObsEnum
Definition: EnumDefinitions.h:852
TimesteppingTypeEnum
@ TimesteppingTypeEnum
Definition: EnumDefinitions.h:436
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
SetVerbosityLevel
void SetVerbosityLevel(int level)
Definition: Verbosity.cpp:34
MaterialsMuWaterEnum
@ MaterialsMuWaterEnum
Definition: EnumDefinitions.h:261
MAXCONNECTIVITY
#define MAXCONNECTIVITY
Definition: CreateNodes.cpp:11
MatestarEnum
@ MatestarEnum
Definition: EnumDefinitions.h:1168
Node::HardDeactivate
void HardDeactivate(void)
Definition: Node.cpp:788
MAXCONNECTIVITY
#define MAXCONNECTIVITY
Definition: CreateElementsVerticesAndMaterials.cpp:10
MaterialsEffectiveconductivityAveragingEnum
@ MaterialsEffectiveconductivityAveragingEnum
Definition: EnumDefinitions.h:252
InversionCostFunctionsEnum
@ InversionCostFunctionsEnum
Definition: EnumDefinitions.h:211
TimesteppingCflCoefficientEnum
@ TimesteppingCflCoefficientEnum
Definition: EnumDefinitions.h:428
BasalforcingsIsmip6DeltaTEnum
@ BasalforcingsIsmip6DeltaTEnum
Definition: EnumDefinitions.h:67
FloatingMeltRateEnum
@ FloatingMeltRateEnum
Definition: EnumDefinitions.h:1069
CreateElements
void CreateElements(Elements *elements, IoModel *iomodel, const int nummodels)
Definition: CreateElementsVerticesAndMaterials.cpp:35
MaterialsThermalExchangeVelocityEnum
@ MaterialsThermalExchangeVelocityEnum
Definition: EnumDefinitions.h:267
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
LinearFloatingMeltRateEnum
@ LinearFloatingMeltRateEnum
Definition: EnumDefinitions.h:1143
MassconaxpbyEnum
@ MassconaxpbyEnum
Definition: EnumDefinitions.h:1161
OneLayerP4zEnum
@ OneLayerP4zEnum
Definition: EnumDefinitions.h:1208
TransientIssmbEnum
@ TransientIssmbEnum
Definition: EnumDefinitions.h:453
MaterialsBetaEnum
@ MaterialsBetaEnum
Definition: EnumDefinitions.h:250
StressbalanceAbstolEnum
@ StressbalanceAbstolEnum
Definition: EnumDefinitions.h:404
BeckmannGoosseFloatingMeltRateEnum
@ BeckmannGoosseFloatingMeltRateEnum
Definition: EnumDefinitions.h:991
QmuResponsedescriptorsEnum
@ QmuResponsedescriptorsEnum
Definition: EnumDefinitions.h:292
MantlePlumeGeothermalFluxEnum
@ MantlePlumeGeothermalFluxEnum
Definition: EnumDefinitions.h:1158
Cfsurfacesquare
Definition: Cfsurfacesquare.h:15
Regionaloutput
Definition: Regionaloutput.h:15
DomainDimensionEnum
@ DomainDimensionEnum
Definition: EnumDefinitions.h:123
CreateNumberNodeToElementConnectivity
void CreateNumberNodeToElementConnectivity(IoModel *iomodel)
Definition: CreateNumberNodeToElementConnectivity.cpp:16
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
SMBgradientscomponentsEnum
@ SMBgradientscomponentsEnum
Definition: EnumDefinitions.h:1248
AmrIceFrontDistanceEnum
@ AmrIceFrontDistanceEnum
Definition: EnumDefinitions.h:24
IoModel::verticaledges
int * verticaledges
Definition: IoModel.h:81
IoModel::elementtoedgeconnectivity
int * elementtoedgeconnectivity
Definition: IoModel.h:83
MatdamageiceEnum
@ MatdamageiceEnum
Definition: EnumDefinitions.h:1165
CreateFaces3d
void CreateFaces3d(IoModel *iomodel)
Definition: CreateFaces.cpp:124
Domain3DsurfaceEnum
@ Domain3DsurfaceEnum
Definition: EnumDefinitions.h:1039
TransientIsgiaEnum
@ TransientIsgiaEnum
Definition: EnumDefinitions.h:446
StressbalanceReltolEnum
@ StressbalanceReltolEnum
Definition: EnumDefinitions.h:410
IoModel::FindConstant
void FindConstant(bool *pvalue, const char *constant_name)
Definition: IoModel.cpp:2362
UNDEF
#define UNDEF
Definition: constants.h:8
TransientIsoceancouplingEnum
@ TransientIsoceancouplingEnum
Definition: EnumDefinitions.h:451
DamageDEnum
@ DamageDEnum
Definition: EnumDefinitions.h:516
StringArrayParam
Definition: StringArrayParam.h:20
P2xP4Enum
@ P2xP4Enum
Definition: EnumDefinitions.h:1227
GroundinglineMeltInterpolationEnum
@ GroundinglineMeltInterpolationEnum
Definition: EnumDefinitions.h:160
AmrErrEnum
@ AmrErrEnum
Definition: EnumDefinitions.h:17
AutodiffFosReverseIndexEnum
@ AutodiffFosReverseIndexEnum
Definition: EnumDefinitions.h:46
IoModel::horizontaledges
int * horizontaledges
Definition: IoModel.h:82
MaterialsHeatcapacityEnum
@ MaterialsHeatcapacityEnum
Definition: EnumDefinitions.h:253
MeshElementtypeEnum
@ MeshElementtypeEnum
Definition: EnumDefinitions.h:271
CreateParametersAutodiff
void CreateParametersAutodiff(Parameters *parameters, IoModel *iomodel)
Definition: CreateParametersAutodiff.cpp:10
Nodes::common_recv
int * common_recv
Definition: Nodes.h:26
IoModel::my_edges
bool * my_edges
Definition: IoModel.h:69
SettingsNumResultsOnNodesEnum
@ SettingsNumResultsOnNodesEnum
Definition: EnumDefinitions.h:335
QmuIsdakotaEnum
@ QmuIsdakotaEnum
Definition: EnumDefinitions.h:288
MINIcondensedEnum
@ MINIcondensedEnum
Definition: EnumDefinitions.h:1157
GeometryHydrostaticRatioEnum
@ GeometryHydrostaticRatioEnum
Definition: EnumDefinitions.h:588
IntParam
Definition: IntParam.h:20
AutodiffIsautodiffEnum
@ AutodiffIsautodiffEnum
Definition: EnumDefinitions.h:50
SpatialLinearFloatingMeltRateEnum
@ SpatialLinearFloatingMeltRateEnum
Definition: EnumDefinitions.h:1278
FixedTimesteppingEnum
@ FixedTimesteppingEnum
Definition: EnumDefinitions.h:1066
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
ThicknessAbsMisfitEnum
@ ThicknessAbsMisfitEnum
Definition: EnumDefinitions.h:837
HydrologydcEplMaxThicknessEnum
@ HydrologydcEplMaxThicknessEnum
Definition: EnumDefinitions.h:181
HydrologydcEplConductivityEnum
@ HydrologydcEplConductivityEnum
Definition: EnumDefinitions.h:179
AutodiffCbufsizeEnum
@ AutodiffCbufsizeEnum
Definition: EnumDefinitions.h:42
TransientIshydrologyEnum
@ TransientIshydrologyEnum
Definition: EnumDefinitions.h:448
MaterialsRheologyEEnum
@ MaterialsRheologyEEnum
Definition: EnumDefinitions.h:645
MaterialsRheologyLawEnum
@ MaterialsRheologyLawEnum
Definition: EnumDefinitions.h:262
InversionSurfaceObsEnum
@ InversionSurfaceObsEnum
Definition: EnumDefinitions.h:630
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
IoModel::numberofhorizontaledges
int numberofhorizontaledges
Definition: IoModel.h:95
TransientAmrFrequencyEnum
@ TransientAmrFrequencyEnum
Definition: EnumDefinitions.h:442
Inputs2
Declaration of Inputs class.
Definition: Inputs2.h:23
AmrThicknessErrorMaximumEnum
@ AmrThicknessErrorMaximumEnum
Definition: EnumDefinitions.h:31
QmuVariablePartitionsNpartEnum
@ QmuVariablePartitionsNpartEnum
Definition: EnumDefinitions.h:295
TransientIscouplerEnum
@ TransientIscouplerEnum
Definition: EnumDefinitions.h:443
MaterialsRheologyEsEnum
@ MaterialsRheologyEsEnum
Definition: EnumDefinitions.h:649
IoModel::domaindim
int domaindim
Definition: IoModel.h:77
BasalforcingsPlumexEnum
@ BasalforcingsPlumexEnum
Definition: EnumDefinitions.h:87
SettingsRecordingFrequencyEnum
@ SettingsRecordingFrequencyEnum
Definition: EnumDefinitions.h:337
BasalforcingsUpperwaterMeltingRateEnum
@ BasalforcingsUpperwaterMeltingRateEnum
Definition: EnumDefinitions.h:95
BasalforcingsPicoFarOceantemperatureEnum
@ BasalforcingsPicoFarOceantemperatureEnum
Definition: EnumDefinitions.h:81
TransientIsthermalEnum
@ TransientIsthermalEnum
Definition: EnumDefinitions.h:455
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
MaterialsRhoSeawaterEnum
@ MaterialsRhoSeawaterEnum
Definition: EnumDefinitions.h:265
BasalforcingsIsmip6Enum
@ BasalforcingsIsmip6Enum
Definition: EnumDefinitions.h:989
SMBgembEnum
@ SMBgembEnum
Definition: EnumDefinitions.h:1246
AmrGradationEnum
@ AmrGradationEnum
Definition: EnumDefinitions.h:19
FSvelocityEnum
@ FSvelocityEnum
Definition: EnumDefinitions.h:1063
P1bubbleEnum
@ P1bubbleEnum
Definition: EnumDefinitions.h:1218
BasalforcingsIsmip6Gamma0Enum
@ BasalforcingsIsmip6Gamma0Enum
Definition: EnumDefinitions.h:68
TransientIsstressbalanceEnum
@ TransientIsstressbalanceEnum
Definition: EnumDefinitions.h:454
HydrologydcWaterCompressibilityEnum
@ HydrologydcWaterCompressibilityEnum
Definition: EnumDefinitions.h:198
DamageDbarEnum
@ DamageDbarEnum
Definition: EnumDefinitions.h:518
Object::ObjectEnum
virtual int ObjectEnum()=0
IsNodeInRank
bool IsNodeInRank(int *nodes_ranks, int *nodes_proc_count, int nid, int rank)
Definition: CreateNodes.cpp:13
AmrDeviatoricErrorMaximumEnum
@ AmrDeviatoricErrorMaximumEnum
Definition: EnumDefinitions.h:14
TimesteppingStartTimeEnum
@ TimesteppingStartTimeEnum
Definition: EnumDefinitions.h:432
EnumToAnalysis
Analysis * EnumToAnalysis(int analysis_enum)
Definition: EnumToAnalysis.cpp:13
StringToEnumx
int StringToEnumx(const char *string_in, bool notfounderror=true)
Definition: StringToEnumx.cpp:14
TimesteppingCouplingTimeEnum
@ TimesteppingCouplingTimeEnum
Definition: EnumDefinitions.h:429
DataSetParam
Definition: DataSetParam.h:20
AmrDeviatoricErrorResolutionEnum
@ AmrDeviatoricErrorResolutionEnum
Definition: EnumDefinitions.h:15
Cfsurfacelogvel
Definition: Cfsurfacelogvel.h:15
MaterialsEnum
@ MaterialsEnum
Definition: EnumDefinitions.h:1167
AmrIceFrontResolutionEnum
@ AmrIceFrontResolutionEnum
Definition: EnumDefinitions.h:25
IoModel::singlenodetoelementconnectivity
int * singlenodetoelementconnectivity
Definition: IoModel.h:100
LoveSolutionEnum
@ LoveSolutionEnum
Definition: EnumDefinitions.h:1155
InversionStepThresholdEnum
@ InversionStepThresholdEnum
Definition: EnumDefinitions.h:225
TimesteppingInterpForcingsEnum
@ TimesteppingInterpForcingsEnum
Definition: EnumDefinitions.h:431
InversionGradientScalingEnum
@ InversionGradientScalingEnum
Definition: EnumDefinitions.h:214
MatlithoEnum
@ MatlithoEnum
Definition: EnumDefinitions.h:1170
Matice
Definition: Matice.h:24
MeshPartitionx
int MeshPartitionx(int **pepart, int **pnpart, int numberofelements, int numberofnodes, int *elements, int numberofelements2d, int numberofnodes2d, doubletype *elements2d, int *vweights, int numlayers, int elements_width, int meshelementtype, int num_procs)
Definition: MeshPartitionx.h:12
DataSet::Presort
void Presort()
Definition: DataSet.cpp:368
RadarEnum
@ RadarEnum
Definition: EnumDefinitions.h:665
Nodes::Finalize
void Finalize(void)
Definition: Nodes.cpp:217
FrontalForcingsParamEnum
@ FrontalForcingsParamEnum
Definition: EnumDefinitions.h:154
HydrologyshaktiEnum
@ HydrologyshaktiEnum
Definition: EnumDefinitions.h:1108
BasalforcingsPicoMaxboxcountEnum
@ BasalforcingsPicoMaxboxcountEnum
Definition: EnumDefinitions.h:84
Element::DatasetInputCreate
virtual void DatasetInputCreate(IssmDouble *array, int M, int N, int *individual_enums, int num_inputs, Inputs2 *inputs2, IoModel *iomodel, int input_enum)
Definition: Element.h:218
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
LACrouzeixRaviartEnum
@ LACrouzeixRaviartEnum
Definition: EnumDefinitions.h:1138
Masscon
Definition: Masscon.h:18
Loads
Declaration of Loads class.
Definition: Loads.h:16
TransientIsesaEnum
@ TransientIsesaEnum
Definition: EnumDefinitions.h:445
Node
Definition: Node.h:23
InversionControlScalingFactorsEnum
@ InversionControlScalingFactorsEnum
Definition: EnumDefinitions.h:210
TimesteppingTimeStepMinEnum
@ TimesteppingTimeStepMinEnum
Definition: EnumDefinitions.h:435
Vertices::Finalize
void Finalize(IoModel *iomodel)
Definition: Vertices.cpp:173
StressbalanceVerticalAnalysisEnum
@ StressbalanceVerticalAnalysisEnum
Definition: EnumDefinitions.h:1289
BasalforcingsDeepwaterElevationEnum
@ BasalforcingsDeepwaterElevationEnum
Definition: EnumDefinitions.h:61
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
CfsurfacesquareEnum
@ CfsurfacesquareEnum
Definition: EnumDefinitions.h:1007
Tria
Definition: Tria.h:29
AmrTypeEnum
@ AmrTypeEnum
Definition: EnumDefinitions.h:34
MeshScaleFactorEnum
@ MeshScaleFactorEnum
Definition: EnumDefinitions.h:652
FlowequationIsFSEnum
@ FlowequationIsFSEnum
Definition: EnumDefinitions.h:138
MaterialsRheologyBEnum
@ MaterialsRheologyBEnum
Definition: EnumDefinitions.h:643
IoModel::faces
int * faces
Definition: IoModel.h:88
IoModel::numberofverticalfaces
int numberofverticalfaces
Definition: IoModel.h:98
StressbalanceRestolEnum
@ StressbalanceRestolEnum
Definition: EnumDefinitions.h:412
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
HydrologypismEnum
@ HydrologypismEnum
Definition: EnumDefinitions.h:1107
HydrologydcSedimentThicknessEnum
@ HydrologydcSedimentThicknessEnum
Definition: EnumDefinitions.h:195
IoModel::numberofverticaledges
int numberofverticaledges
Definition: IoModel.h:94
BasalforcingsDtbgEnum
@ BasalforcingsDtbgEnum
Definition: EnumDefinitions.h:63
TransientIsmovingfrontEnum
@ TransientIsmovingfrontEnum
Definition: EnumDefinitions.h:450
DoubleVecParam
Definition: DoubleVecParam.h:20
UpdateElementsAndMaterialsDakota
void UpdateElementsAndMaterialsDakota(Elements *elements, Inputs2 *inputs2, Materials *materials, IoModel *iomodel)
Definition: UpdateElementsAndMaterialsDakota.cpp:11
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
WeightsSurfaceObservationEnum
@ WeightsSurfaceObservationEnum
Definition: EnumDefinitions.h:864
SettingsOutputFrequencyEnum
@ SettingsOutputFrequencyEnum
Definition: EnumDefinitions.h:336
PentaEnum
@ PentaEnum
Definition: EnumDefinitions.h:1231
IoModel::NumIndependents
int NumIndependents()
Definition: IoModel.cpp:2448
SurfaceLogVelMisfitEnum
@ SurfaceLogVelMisfitEnum
Definition: EnumDefinitions.h:825
AdolcParamEnum
@ AdolcParamEnum
Definition: EnumDefinitions.h:12
InversionMaxiterEnum
@ InversionMaxiterEnum
Definition: EnumDefinitions.h:219
CreateMaterials
void CreateMaterials(Elements *elements, Inputs2 *inputs2, Materials *materials, IoModel *iomodel, const int nummodels)
Definition: CreateElementsVerticesAndMaterials.cpp:87
AmrBamgEnum
@ AmrBamgEnum
Definition: EnumDefinitions.h:974
DamageEvolutionRequestedOutputsEnum
@ DamageEvolutionRequestedOutputsEnum
Definition: EnumDefinitions.h:114
Massconaxpby
Definition: Massconaxpby.h:18
AutodiffNumDependentsEnum
@ AutodiffNumDependentsEnum
Definition: EnumDefinitions.h:52
StringParam
Definition: StringParam.h:20
BasalforcingsPicoFarOceansalinityEnum
@ BasalforcingsPicoFarOceansalinityEnum
Definition: EnumDefinitions.h:80
CreateOutputDefinitions
void CreateOutputDefinitions(Elements *elements, Parameters *parameters, Inputs2 *inputs2, IoModel *iomodel)
Definition: CreateOutputDefinitions.cpp:10
Matlitho
Definition: Matlitho.h:14
Element::ControlInputCreate
void ControlInputCreate(IssmDouble *doublearray, IssmDouble *independents_min, IssmDouble *independents_max, Inputs2 *inputs2, IoModel *iomodel, int M, int N, IssmDouble scale, int input_enum, int id)
Definition: Element.cpp:1753
ISSM_MPI_Barrier
int ISSM_MPI_Barrier(ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:148
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
InversionTypeEnum
@ InversionTypeEnum
Definition: EnumDefinitions.h:226
InversionGatolEnum
@ InversionGatolEnum
Definition: EnumDefinitions.h:213
AmrKeepMetricEnum
@ AmrKeepMetricEnum
Definition: EnumDefinitions.h:26
IoModel::FetchDataToInput
void FetchDataToInput(Inputs2 *inputs2, Elements *elements, const char *vector_name, int input_enum)
Definition: IoModel.cpp:1651
InversionDxminEnum
@ InversionDxminEnum
Definition: EnumDefinitions.h:212
FacesPartitioning
void FacesPartitioning(IoModel *iomodel)
Definition: FacesPartitioning.cpp:10
SteadystateMaxiterEnum
@ SteadystateMaxiterEnum
Definition: EnumDefinitions.h:399
BasalforcingsUppercrustthicknessEnum
@ BasalforcingsUppercrustthicknessEnum
Definition: EnumDefinitions.h:92
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
HydrologyGlaDSEnum
@ HydrologyGlaDSEnum
Definition: EnumDefinitions.h:1101
MatenhancediceEnum
@ MatenhancediceEnum
Definition: EnumDefinitions.h:1166
DoubleMatParam
Definition: DoubleMatParam.h:20
AmrThicknessErrorGroupThresholdEnum
@ AmrThicknessErrorGroupThresholdEnum
Definition: EnumDefinitions.h:30
IoModel::epart
int * epart
Definition: IoModel.h:74
CreateParametersControl
void CreateParametersControl(Parameters *parameters, IoModel *iomodel, int solution_type)
Definition: CreateParametersControl.cpp:10
IoModel::elements
int * elements
Definition: IoModel.h:79
GroundinglineFrictionInterpolationEnum
@ GroundinglineFrictionInterpolationEnum
Definition: EnumDefinitions.h:159
MaticeEnum
@ MaticeEnum
Definition: EnumDefinitions.h:1169
IoModel::edges
int * edges
Definition: IoModel.h:80
SMBpddSicopolisEnum
@ SMBpddSicopolisEnum
Definition: EnumDefinitions.h:1253
InversionVxObsEnum
@ InversionVxObsEnum
Definition: EnumDefinitions.h:633
Vertex
Definition: Vertex.h:19
SurfaceAverageVelMisfitEnum
@ SurfaceAverageVelMisfitEnum
Definition: EnumDefinitions.h:821
DependentObject::NumDependents
int NumDependents(void)
Definition: DependentObject.cpp:81
AutodiffGcTriggerRatioEnum
@ AutodiffGcTriggerRatioEnum
Definition: EnumDefinitions.h:49
IoModel::meshelementtype
int meshelementtype
Definition: IoModel.h:91
MaterialsRheologyEcbarEnum
@ MaterialsRheologyEcbarEnum
Definition: EnumDefinitions.h:648
CreateFaces
void CreateFaces(IoModel *iomodel)
Definition: CreateFaces.cpp:9
ParseToolkitsOptionsx
void ParseToolkitsOptionsx(Parameters *parameters, FILE *fid)
Definition: ParseToolkitsOptionsx.cpp:19
SettingsWaitonlockEnum
@ SettingsWaitonlockEnum
Definition: EnumDefinitions.h:341
Element::DatasetInputAdd
void DatasetInputAdd(int enum_type, IssmDouble *vector, Inputs2 *inputs2, IoModel *iomodel, int M, int N, int vector_type, int vector_enum, int code, int input_enum)
Definition: Element.cpp:1831
HydrologydcEplCompressibilityEnum
@ HydrologydcEplCompressibilityEnum
Definition: EnumDefinitions.h:178
P1xP4Enum
@ P1xP4Enum
Definition: EnumDefinitions.h:1222
BoolParam
Definition: BoolParam.h:20
InversionGttolEnum
@ InversionGttolEnum
Definition: EnumDefinitions.h:216
MaterialsTemperateiceconductivityEnum
@ MaterialsTemperateiceconductivityEnum
Definition: EnumDefinitions.h:266
ControlInputSizeNEnum
@ ControlInputSizeNEnum
Definition: EnumDefinitions.h:106
IoModel::elementtoverticalfaceconnectivity
int * elementtoverticalfaceconnectivity
Definition: IoModel.h:87
HydrologyshreveEnum
@ HydrologyshreveEnum
Definition: EnumDefinitions.h:1109
AutodiffTapeAllocEnum
@ AutodiffTapeAllocEnum
Definition: EnumDefinitions.h:55
TimesteppingTimeStepMaxEnum
@ TimesteppingTimeStepMaxEnum
Definition: EnumDefinitions.h:434
Loads::Finalize
void Finalize()
Definition: Loads.cpp:97
BasalforcingsUpperdepthMeltEnum
@ BasalforcingsUpperdepthMeltEnum
Definition: EnumDefinitions.h:93
MassFluxSegmentsPresentEnum
@ MassFluxSegmentsPresentEnum
Definition: EnumDefinitions.h:242
P1xP3Enum
@ P1xP3Enum
Definition: EnumDefinitions.h:1221
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
SmbASnowEnum
@ SmbASnowEnum
Definition: EnumDefinitions.h:344
SmbRlapslgmEnum
@ SmbRlapslgmEnum
Definition: EnumDefinitions.h:384
IoModel::my_faces
bool * my_faces
Definition: IoModel.h:67
MassFluxSegmentsEnum
@ MassFluxSegmentsEnum
Definition: EnumDefinitions.h:241
SteadystateReltolEnum
@ SteadystateReltolEnum
Definition: EnumDefinitions.h:401
InversionAlgorithmEnum
@ InversionAlgorithmEnum
Definition: EnumDefinitions.h:208
IoModel::domaintype
int domaintype
Definition: IoModel.h:78
SmbRlapsEnum
@ SmbRlapsEnum
Definition: EnumDefinitions.h:383
TransientParam
Definition: TransientParam.h:20
InversionGrtolEnum
@ InversionGrtolEnum
Definition: EnumDefinitions.h:215
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
AutodiffNumIndependentsEnum
@ AutodiffNumIndependentsEnum
Definition: EnumDefinitions.h:53
AddVertexToRank
void AddVertexToRank(int *vertices_ranks, int *vertices_proc_count, int vid, int rank)
Definition: CreateElementsVerticesAndMaterials.cpp:24
AutodiffGcTriggerMaxSizeEnum
@ AutodiffGcTriggerMaxSizeEnum
Definition: EnumDefinitions.h:48
ToolkitOptions::Init
static void Init(void)
Definition: ToolkitOptions.cpp:25
DoubleTransientMatParam
Definition: DoubleTransientMatParam.h:20
SMBmeltcomponentsEnum
@ SMBmeltcomponentsEnum
Definition: EnumDefinitions.h:1251
GslEnum
@ GslEnum
Definition: EnumDefinitions.h:1094
SmbDesfacEnum
@ SmbDesfacEnum
Definition: EnumDefinitions.h:350
HydrologydcEplPorosityEnum
@ HydrologydcEplPorosityEnum
Definition: EnumDefinitions.h:182
SmbRdlEnum
@ SmbRdlEnum
Definition: EnumDefinitions.h:381
Analysis::UpdateParameters
virtual void UpdateParameters(Parameters *parameters, IoModel *iomodel, int solution_enum, int analysis_enum)=0
AutodiffTbufsizeEnum
@ AutodiffTbufsizeEnum
Definition: EnumDefinitions.h:56
AmrDeviatoricErrorGroupThresholdEnum
@ AmrDeviatoricErrorGroupThresholdEnum
Definition: EnumDefinitions.h:13
IoModel::FillIndependents
void FillIndependents(IssmDouble *xp)
Definition: IoModel.cpp:2334
RheologyBInitialguessMisfitEnum
@ RheologyBInitialguessMisfitEnum
Definition: EnumDefinitions.h:673
BalancethicknessSpcthicknessEnum
@ BalancethicknessSpcthicknessEnum
Definition: EnumDefinitions.h:986
DoubleMatArrayParam
Definition: DoubleMatArrayParam.h:20
TriaEnum
@ TriaEnum
Definition: EnumDefinitions.h:1318
ElementsAndVerticesPartitioning
void ElementsAndVerticesPartitioning(IoModel *iomodel)
Definition: ElementsAndVerticesPartitioning.cpp:17
BasalforcingsMeltrateFactorEnum
@ BasalforcingsMeltrateFactorEnum
Definition: EnumDefinitions.h:74
SettingsIoGatherEnum
@ SettingsIoGatherEnum
Definition: EnumDefinitions.h:334
Analysis
Definition: Analysis.h:30
Vertices::common_send_ids
int ** common_send_ids
Definition: Vertices.h:25
XTaylorHoodEnum
@ XTaylorHoodEnum
Definition: EnumDefinitions.h:1329
P1xP2Enum
@ P1xP2Enum
Definition: EnumDefinitions.h:1220
Inputs2::DuplicateInput
void DuplicateInput(int original_enum, int new_enum)
Definition: Inputs2.cpp:206
AmrFieldEnum
@ AmrFieldEnum
Definition: EnumDefinitions.h:18
SMBforcingEnum
@ SMBforcingEnum
Definition: EnumDefinitions.h:1244
SMBpddEnum
@ SMBpddEnum
Definition: EnumDefinitions.h:1252
MaterialsMixedLayerCapacityEnum
@ MaterialsMixedLayerCapacityEnum
Definition: EnumDefinitions.h:260
IoModel::numberofedges
int numberofedges
Definition: IoModel.h:93
HeapSort
void HeapSort(T *c, long n)
Definition: HeapSort.h:5
QmuVariablePartitionsEnum
@ QmuVariablePartitionsEnum
Definition: EnumDefinitions.h:294
GenericParam
Definition: GenericParam.h:26
Analysis::CreateNodes
virtual void CreateNodes(Nodes *nodes, IoModel *iomodel, bool isamr=false)=0
BasalforcingsPicoIsplumeEnum
@ BasalforcingsPicoIsplumeEnum
Definition: EnumDefinitions.h:83
BasalforcingsBottomplumedepthEnum
@ BasalforcingsBottomplumedepthEnum
Definition: EnumDefinitions.h:59
BasalforcingsPlumeradiusEnum
@ BasalforcingsPlumeradiusEnum
Definition: EnumDefinitions.h:86
VxObsEnum
@ VxObsEnum
Definition: EnumDefinitions.h:848
IoModel::my_vedges
bool * my_vedges
Definition: IoModel.h:70
TransientIsdamageevolutionEnum
@ TransientIsdamageevolutionEnum
Definition: EnumDefinitions.h:444
BasalforcingsMantleconductivityEnum
@ BasalforcingsMantleconductivityEnum
Definition: EnumDefinitions.h:73
Misfit
Definition: Misfit.h:15