Ice Sheet System Model  4.18
Code documentation
Functions
IoModelToConstraintsx.cpp File Reference
#include "./IoModelToConstraintsx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../ModelProcessorx/ModelProcessorx.h"

Go to the source code of this file.

Functions

void IoModelToConstraintsx (Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
 
void IoModelToDynamicConstraintsx (Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
 
void IoModelToConstraintsx (Constraints *constraints, IoModel *iomodel, IssmDouble *spcdata, int M, int N, int analysis_type, int finite_element, int dof)
 
void IoModelToDynamicConstraintsx (Constraints *constraints, IoModel *iomodel, IssmDouble *spcdata, int M, int N, int analysis_type, int finite_element, int dof)
 

Function Documentation

◆ IoModelToConstraintsx() [1/2]

void IoModelToConstraintsx ( Constraints *  constraints,
IoModel *  iomodel,
const char *  spc_name,
int  analysis_type,
int  finite_element,
int  dof 
)

Definition at line 10 of file IoModelToConstraintsx.cpp.

10  {
11 
12  /*intermediary: */
13  int code,vector_layout;
14  IssmDouble *spcdata = NULL;
15  int M,N;
16 
17  /*First of, find the record for the enum, and get code of data type: */
18  iomodel->SetFilePointerToData(&code, &vector_layout,spc_name);
19  if(code!=7)_error_("expecting a IssmDouble vector for constraints " << spc_name);
20  if(vector_layout!=1)_error_("expecting a nodal vector for constraints " << spc_name);
21 
22  /*Fetch vector:*/
23  iomodel->FetchData(&spcdata,&M,&N,spc_name);
24 
25  /*Call IoModelToConstraintsx*/
26  IoModelToConstraintsx(constraints,iomodel,spcdata,M,N,analysis_type,finite_element,dof);
27 
28  /*Clean up*/
29  xDelete<IssmDouble>(spcdata);
30 }

◆ IoModelToDynamicConstraintsx() [1/2]

void IoModelToDynamicConstraintsx ( Constraints *  constraints,
IoModel *  iomodel,
const char *  spc_name,
int  analysis_type,
int  finite_element,
int  dof 
)

Definition at line 31 of file IoModelToConstraintsx.cpp.

31  {
32 
33  /*intermediary: */
34  int code,vector_layout;
35  IssmDouble *spcdata = NULL;
36  int M,N;
37 
38  /*First of, find the record for the enum, and get code of data type: */
39  iomodel->SetFilePointerToData(&code, &vector_layout,spc_name);
40  if(code!=7)_error_("expecting a IssmDouble vector for constraints " << spc_name);
41  if(vector_layout!=1)_error_("expecting a nodal vector for constraints " << spc_name);
42 
43  /*Fetch vector:*/
44  iomodel->FetchData(&spcdata,&M,&N,spc_name);
45 
46  /*Call IoModelToConstraintsx*/
47  IoModelToDynamicConstraintsx(constraints,iomodel,spcdata,M,N,analysis_type,finite_element,dof);
48 
49  /*Clean up*/
50  xDelete<IssmDouble>(spcdata);
51 }

◆ IoModelToConstraintsx() [2/2]

void IoModelToConstraintsx ( Constraints *  constraints,
IoModel *  iomodel,
IssmDouble *  spcdata,
int  M,
int  N,
int  analysis_type,
int  finite_element,
int  dof 
)

Definition at line 53 of file IoModelToConstraintsx.cpp.

53  {/*{{{*/
54 
55  /*intermediary: */
56  int i,j,elementnbv,numfacevertices;
57  IssmDouble value;
58  IssmDouble *times = NULL;
59  IssmDouble *values = NULL;
60  bool spcpresent = false;
61 
62  /*Higher-order finite elements*/
63  int v1,v2;
64  bool *boundaryedge = NULL;
65 
66  switch(finite_element){
67  case P1Enum:
68  /*Nothing else to do*/
69  break;
70  case P1bubbleEnum:
71  switch(iomodel->meshelementtype){
72  case TriaEnum: elementnbv = 3; break;
73  case TetraEnum: elementnbv = 4; break;
74  case PentaEnum: elementnbv = 6; break;
75  default: _error_("mesh type not supported yet");
76  }
77  break;
78  case P1bubblecondensedEnum:
79  /*Nothing else to do*/
80  break;
81  case P1xP2Enum:
82  case P1xP3Enum:
83  case P1xP4Enum:
84  EdgesPartitioning(iomodel);
85  break;
86  case P2xP1Enum:
87  EdgesPartitioning(iomodel);
88  break;
89  case P2Enum:
90  EdgesPartitioning(iomodel);
91  if(iomodel->meshelementtype==PentaEnum) FacesPartitioning(iomodel);
92  EdgeOnBoundaryFlags(&boundaryedge,iomodel);
93  break;
94  case P2bubbleEnum:
95  EdgesPartitioning(iomodel);
96  if(iomodel->meshelementtype==PentaEnum){
97  FacesPartitioning(iomodel);
98  }
99  EdgeOnBoundaryFlags(&boundaryedge,iomodel);
100  switch(iomodel->meshelementtype){
101  case TriaEnum: elementnbv = 3; break;
102  case TetraEnum: elementnbv = 4; break;
103  case PentaEnum: elementnbv = 6; break;
104  default: _error_("mesh type not supported yet");
105  }
106  break;
107  case P2xP4Enum:
108  EdgesPartitioning(iomodel);
109  FacesPartitioning(iomodel);
110  break;
111  default:
112  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
113  }
114 
115  int count =0;
116  if(M==iomodel->numberofvertices){
117  switch(finite_element){
118  case P1Enum:
119  for(i=0;i<iomodel->numberofvertices;i++){
120  if((iomodel->my_vertices[i])){
121  if (!xIsNan<IssmDouble>(spcdata[i])){
122  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
123  count++;
124  }
125  }
126  }
127  break;
128  case P2Enum:
129  for(i=0;i<iomodel->numberofvertices;i++){
130  if((iomodel->my_vertices[i])){
131  if (!xIsNan<IssmDouble>(spcdata[i])){
132  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
133  count++;
134  }
135  }
136  }
137  for(i=0;i<iomodel->numberofedges;i++){
138  if(iomodel->my_edges[i] && boundaryedge[i]){
139  v1 = iomodel->edges[3*i+0]-1;
140  v2 = iomodel->edges[3*i+1]-1;
141  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
142  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
143  dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
144  count++;
145  }
146  }
147  }
148  if(iomodel->meshelementtype==PentaEnum){
149  for(i=0;i<iomodel->numberofverticalfaces;i++){
150  if(iomodel->my_vfaces[i]){
151  value=0.;
152  for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j] -1]/4.;
153  if(!xIsNan<IssmDouble>(value)){
154  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+i+1,dof,value,analysis_type));
155  count++;
156  }
157  }
158  }
159  }
160  break;
161  case P2bubbleEnum:
162  for(i=0;i<iomodel->numberofvertices;i++){
163  if((iomodel->my_vertices[i])){
164  if (!xIsNan<IssmDouble>(spcdata[i])){
165  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
166  count++;
167  }
168  }
169  }
170  for(i=0;i<iomodel->numberofedges;i++){
171  if(iomodel->my_edges[i] && boundaryedge[i]){
172  v1 = iomodel->edges[3*i+0]-1;
173  v2 = iomodel->edges[3*i+1]-1;
174  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
175  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
176  dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
177  count++;
178  }
179  }
180  }
181  if(iomodel->meshelementtype==PentaEnum){
182  for(i=0;i<iomodel->numberoffaces;i++){
183  if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
184  if(iomodel->my_faces[i]){
185  numfacevertices = iomodel->faces[i*iomodel->facescols+3];
186  value=0.;
187  for(j=0;j<numfacevertices;j++){
188  value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1];
189  }
190  value = value/reCast<IssmDouble>(numfacevertices);
191  if(!xIsNan<IssmDouble>(value)){
192  constraints->AddObject(new SpcStatic(count+1,
193  iomodel->numberofvertices+iomodel->numberofedges+i+1,
194  dof,value,analysis_type));
195  count++;
196  }
197  }
198  }
199  }
200  }
201  for(i=0;i<iomodel->numberofelements;i++){
202  if(iomodel->my_elements[i]){
203  value = spcdata[iomodel->elements[i*elementnbv+0]-1];
204  for(j=1;j<elementnbv;j++) value += spcdata[iomodel->elements[i*elementnbv+j]-1];
205  value = value/reCast<IssmDouble,int>(elementnbv+0);
206  if(!xIsNan<IssmDouble>(value)){
207  int nodeid = iomodel->numberofvertices+iomodel->numberofedges+i+1;
208  if(iomodel->meshelementtype==PentaEnum){
209  nodeid += iomodel->numberoffaces;
210  }
211  constraints->AddObject(new SpcStatic(count+1,nodeid,dof,value,analysis_type));
212  count++;
213  }
214  }
215  }
216  break;
217  case P2xP4Enum:
218  for(i=0;i<iomodel->numberofvertices;i++){
219  if((iomodel->my_vertices[i])){
220  if (!xIsNan<IssmDouble>(spcdata[i])){
221  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
222  count++;
223  }
224  }
225  }
226  for(i=0;i<iomodel->numberofedges;i++){
227  if(iomodel->my_edges[i]){
228  v1 = iomodel->edges[3*i+0]-1;
229  v2 = iomodel->edges[3*i+1]-1;
230  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
231  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
232  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
233  count++;
234  }
235  }
236  }
237  for(i=0;i<iomodel->numberofverticaledges;i++){
238  if(iomodel->my_vedges[i]){
239  v1 = iomodel->verticaledges[2*i+0]-1;
240  v2 = iomodel->verticaledges[2*i+1]-1;
241  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
242  /*FIXME: should not be 1/2 but 1/4 and 3/4!*/
243  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+2*i+1,
244  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
245  constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+iomodel->numberofedges+2*i+2,
246  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
247  count=count+2;
248  }
249  }
250  }
251  for(i=0;i<iomodel->numberofverticalfaces;i++){
252  if(iomodel->my_vfaces[i]){
253  value=0.;
254  for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j]-1]/4.;
255  if(!xIsNan<IssmDouble>(value)){
256  constraints->AddObject(new SpcStatic(count+1,
257  iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+1,
258  dof,value,analysis_type));
259  constraints->AddObject(new SpcStatic(count+2,
260  iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+2,
261  dof,value,analysis_type));
262  constraints->AddObject(new SpcStatic(count+3,
263  iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+3,
264  dof,value,analysis_type));
265  count=count+3;
266  }
267  }
268  }
269  break;
270  case P1bubbleEnum:
271  for(i=0;i<iomodel->numberofvertices;i++){
272  if((iomodel->my_vertices[i])){
273  if (!xIsNan<IssmDouble>(spcdata[i])){
274  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
275  count++;
276  }
277  }
278  }
279  for(i=0;i<iomodel->numberofelements;i++){
280  if(iomodel->my_elements[i]){
281  value = spcdata[iomodel->elements[i*elementnbv+0]-1];
282  for(j=1;j<elementnbv;j++) value += spcdata[iomodel->elements[i*elementnbv+j]-1];
283  value = value/reCast<IssmDouble,int>(elementnbv+0);
284  if(!xIsNan<IssmDouble>(value)){
285  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
286  dof,value,analysis_type));
287  count++;
288  }
289  }
290  }
291  break;
292  case P1bubblecondensedEnum:
293  for(i=0;i<iomodel->numberofvertices;i++){
294  if((iomodel->my_vertices[i])){
295  if (!xIsNan<IssmDouble>(spcdata[i])){
296  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
297  count++;
298  }
299  }
300  }
301  break;
302  case P1xP2Enum:
303  for(i=0;i<iomodel->numberofvertices;i++){
304  if((iomodel->my_vertices[i])){
305  if (!xIsNan<IssmDouble>(spcdata[i])){
306  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
307  count++;
308  }
309  }
310  }
311  for(i=0;i<iomodel->numberofverticaledges;i++){
312  if(iomodel->my_vedges[i]){
313  v1 = iomodel->verticaledges[2*i+0]-1;
314  v2 = iomodel->verticaledges[2*i+1]-1;
315  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
316  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
317  count++;
318  }
319  }
320  }
321  break;
322  case P1xP3Enum:
323  for(i=0;i<iomodel->numberofvertices;i++){
324  if((iomodel->my_vertices[i])){
325  if (!xIsNan<IssmDouble>(spcdata[i])){
326  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
327  count++;
328  }
329  }
330  }
331  for(i=0;i<iomodel->numberofverticaledges;i++){
332  if(iomodel->my_vedges[i]){
333  v1 = iomodel->verticaledges[2*i+0]-1;
334  v2 = iomodel->verticaledges[2*i+1]-1;
335  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
336  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+2*i+1,
337  dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type));
338  constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+2*i+2,
339  dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type));
340  count=count+2;
341  }
342  }
343  }
344  break;
345  case P1xP4Enum:
346  for(i=0;i<iomodel->numberofvertices;i++){
347  if((iomodel->my_vertices[i])){
348  if (!xIsNan<IssmDouble>(spcdata[i])){
349  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
350  count++;
351  }
352  }
353  }
354  for(i=0;i<iomodel->numberofverticaledges;i++){
355  if(iomodel->my_vedges[i]){
356  v1 = iomodel->verticaledges[2*i+0]-1;
357  v2 = iomodel->verticaledges[2*i+1]-1;
358  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
359  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1,
360  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
361  constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2,
362  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
363  constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3,
364  dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
365  count=count+3;
366  }
367  }
368  }
369  break;
370  case P2xP1Enum:
371  for(i=0;i<iomodel->numberofvertices;i++){
372  if((iomodel->my_vertices[i])){
373  if (!xIsNan<IssmDouble>(spcdata[i])){
374  constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
375  count++;
376  }
377  }
378  }
379  for(i=0;i<iomodel->numberofhorizontaledges;i++){
380  if(iomodel->my_hedges[i]){
381  v1 = iomodel->horizontaledges[2*i+0]-1;
382  v2 = iomodel->horizontaledges[2*i+1]-1;
383  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
384  constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
385  count++;
386  }
387  }
388  }
389  break;
390  default:
391  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
392  }
393  }
394  else if (M==(iomodel->numberofvertices+1)){
395  /*transient: create transient SpcTransient objects. Same logic, except we need to retrieve
396  * various times and values to initialize an SpcTransient object: */
397 
398  /*figure out times: */
399  times=xNew<IssmDouble>(N);
400  for(j=0;j<N;j++) times[j]=spcdata[(M-1)*N+j];
401 
402  switch(finite_element){
403  case P1Enum:
404  for(i=0;i<iomodel->numberofvertices;i++){
405  if((iomodel->my_vertices[i])){
406 
407  /*figure out times and values: */
408  values=xNew<IssmDouble>(N);
409  spcpresent=false;
410  for(j=0;j<N;j++){
411  values[j]=spcdata[i*N+j];
412  if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
413  }
414 
415  if(spcpresent){
416  constraints->AddObject(new SpcTransient(count+1,i+1,dof,N,times,values,analysis_type));
417  count++;
418  }
419  xDelete<IssmDouble>(values);
420  }
421  }
422  break;
423  case P2Enum:
424  for(i=0;i<iomodel->numberofvertices;i++){
425  if((iomodel->my_vertices[i])){
426 
427  /*figure out times and values: */
428  values=xNew<IssmDouble>(N);
429  spcpresent=false;
430  for(j=0;j<N;j++){
431  values[j]=spcdata[i*N+j];
432  if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
433  }
434  if(iomodel->meshelementtype==PentaEnum){ _error_("need to add nodes on faces! (See Spcstatic above)");}
435 
436  if(spcpresent){
437  constraints->AddObject(new SpcTransient(count+1,i+1,dof,N,times,values,analysis_type));
438  count++;
439  }
440  xDelete<IssmDouble>(values);
441  }
442  }
443  for(i=0;i<iomodel->numberofedges;i++){
444  if(iomodel->my_edges[i]){
445  v1 = iomodel->edges[3*i+0]-1;
446  v2 = iomodel->edges[3*i+1]-1;
447  values=xNew<IssmDouble>(N);
448  spcpresent=false;
449  for(j=0;j<N;j++){
450  values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
451  if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
452  }
453  if(spcpresent){
454  constraints->AddObject(new SpcTransient(count+1,iomodel->numberofvertices+i+1,dof,N,times,values,analysis_type));
455  count++;
456  }
457  xDelete<IssmDouble>(values);
458  }
459  }
460  break;
461  case P1xP2Enum:
462  for(i=0;i<iomodel->numberofvertices;i++){
463  if((iomodel->my_vertices[i])){
464 
465  /*figure out times and values: */
466  values=xNew<IssmDouble>(N);
467  spcpresent=false;
468  for(j=0;j<N;j++){
469  values[j]=spcdata[i*N+j];
470  if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
471  }
472 
473  if(spcpresent){
474  constraints->AddObject(new SpcTransient(count+1,i+1,dof,N,times,values,analysis_type));
475  count++;
476  }
477  xDelete<IssmDouble>(values);
478  }
479  }
480  for(i=0;i<iomodel->numberofverticaledges;i++){
481  if(iomodel->my_vedges[i]){
482  v1 = iomodel->verticaledges[2*i+0]-1;
483  v2 = iomodel->verticaledges[2*i+1]-1;
484  values=xNew<IssmDouble>(N);
485  spcpresent=false;
486  for(j=0;j<N;j++){
487  values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
488  if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
489  }
490  if(spcpresent){
491  constraints->AddObject(new SpcTransient(count+1,iomodel->numberofvertices+i+1,dof,N,times,values,analysis_type));
492  count++;
493  }
494  xDelete<IssmDouble>(values);
495  }
496  }
497  break;
498  case P1xP3Enum:
499  for(i=0;i<iomodel->numberofvertices;i++){
500  if((iomodel->my_vertices[i])){
501 
502  /*figure out times and values: */
503  values=xNew<IssmDouble>(N);
504  spcpresent=false;
505  for(j=0;j<N;j++){
506  values[j]=spcdata[i*N+j];
507  if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
508  }
509 
510  if(spcpresent){
511  constraints->AddObject(new SpcTransient(count+1,i+1,dof,N,times,values,analysis_type));
512  count++;
513  }
514  xDelete<IssmDouble>(values);
515  }
516  }
517  for(i=0;i<iomodel->numberofverticaledges;i++){
518  if(iomodel->my_vedges[i]){
519  v1 = iomodel->verticaledges[2*i+0]-1;
520  v2 = iomodel->verticaledges[2*i+1]-1;
521  values=xNew<IssmDouble>(N);
522  spcpresent=false;
523  for(j=0;j<N;j++){
524  values[j]=2./3.*spcdata[v1*N+j]+1./3.*spcdata[v2*N+j];
525  if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
526  }
527  if(spcpresent){
528  constraints->AddObject(new SpcTransient(count+1,iomodel->numberofvertices+2*i+1,dof,N,times,values,analysis_type));
529  count++;
530  }
531  spcpresent=false;
532  for(j=0;j<N;j++){
533  values[j]=1./3.*spcdata[v1*N+j]+2./3.*spcdata[v2*N+j];
534  if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
535  }
536  if(spcpresent){
537  constraints->AddObject(new SpcTransient(count+1,iomodel->numberofvertices+2*i+2,dof,N,times,values,analysis_type));
538  count++;
539  }
540  xDelete<IssmDouble>(values);
541  }
542  }
543  break;
544  case P2xP1Enum:
545  for(i=0;i<iomodel->numberofvertices;i++){
546  if((iomodel->my_vertices[i])){
547 
548  /*figure out times and values: */
549  values=xNew<IssmDouble>(N);
550  spcpresent=false;
551  for(j=0;j<N;j++){
552  values[j]=spcdata[i*N+j];
553  if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
554  }
555 
556  if(spcpresent){
557  constraints->AddObject(new SpcTransient(count+1,i+1,dof,N,times,values,analysis_type));
558  count++;
559  }
560  xDelete<IssmDouble>(values);
561  }
562  }
563  for(i=0;i<iomodel->numberofedges;i++){
564  if(iomodel->edges[i*3+2]!=2){
565  if(iomodel->my_edges[i]){
566  v1 = iomodel->edges[3*i+0]-1;
567  v2 = iomodel->edges[3*i+1]-1;
568  values=xNew<IssmDouble>(N);
569  spcpresent=false;
570  for(j=0;j<N;j++){
571  values[j]=(spcdata[v1*N+j]+spcdata[v2*N+j])/2.;
572  if(!xIsNan<IssmDouble>(values[j])) spcpresent=true; //NaN means no spc by default
573  }
574  if(spcpresent){
575  constraints->AddObject(new SpcTransient(count+1,iomodel->numberofvertices+i+1,dof,N,times,values,analysis_type));
576  count++;
577  }
578  xDelete<IssmDouble>(values);
579  }
580  }
581  }
582  break;
583  default:
584  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
585  }
586  }
587  else{
588  _error_("Size of spc field not supported");
589  }
590 
591  /*Free ressources:*/
592  xDelete<IssmDouble>(times);
593  xDelete<IssmDouble>(values);
594  xDelete<bool>(boundaryedge);
595 }/*}}}*/

◆ IoModelToDynamicConstraintsx() [2/2]

void IoModelToDynamicConstraintsx ( Constraints *  constraints,
IoModel *  iomodel,
IssmDouble *  spcdata,
int  M,
int  N,
int  analysis_type,
int  finite_element,
int  dof 
)

Definition at line 596 of file IoModelToConstraintsx.cpp.

596  {/*{{{*/
597 
598  /*intermediary: */
599  int i,j,elementnbv,numfacevertices;
600  IssmDouble value;
601  IssmDouble *times = NULL;
602  IssmDouble *values = NULL;
603  bool spcpresent = false;
604 
605  /*Higher-order finite elements*/
606  int v1,v2;
607  bool *boundaryedge = NULL;
608 
609  switch(finite_element){
610  case P1Enum:
611  /*Nothing else to do*/
612  break;
613  case P1bubbleEnum:
614  switch(iomodel->meshelementtype){
615  case TriaEnum: elementnbv = 3; break;
616  case TetraEnum: elementnbv = 4; break;
617  case PentaEnum: elementnbv = 6; break;
618  default: _error_("mesh type not supported yet");
619  }
620  break;
621  case P1bubblecondensedEnum:
622  /*Nothing else to do*/
623  break;
624  case P1xP2Enum:
625  EdgesPartitioning(iomodel);
626  break;
627  case P1xP3Enum:
628  EdgesPartitioning(iomodel);
629  break;
630  case P2xP1Enum:
631  EdgesPartitioning(iomodel);
632  break;
633  case P2Enum:
634  EdgesPartitioning(iomodel);
635  if(iomodel->meshelementtype==PentaEnum){
636  FacesPartitioning(iomodel);
637  }
638  EdgeOnBoundaryFlags(&boundaryedge,iomodel);
639  break;
640  case P2bubbleEnum:
641  EdgesPartitioning(iomodel);
642  if(iomodel->meshelementtype==PentaEnum){
643  FacesPartitioning(iomodel);
644  }
645  EdgeOnBoundaryFlags(&boundaryedge,iomodel);
646  switch(iomodel->meshelementtype){
647  case TriaEnum: elementnbv = 3; break;
648  case TetraEnum: elementnbv = 4; break;
649  case PentaEnum: elementnbv = 6; break;
650  default: _error_("mesh type not supported yet");
651  }
652  break;
653  case P2xP4Enum:
654  EdgesPartitioning(iomodel);
655  FacesPartitioning(iomodel);
656  break;
657  default:
658  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
659  }
660 
661  int count=0;
662  if(M==iomodel->numberofvertices){
663  switch(finite_element){
664  case P1Enum:
665  for(i=0;i<iomodel->numberofvertices;i++){
666  if((iomodel->my_vertices[i])){
667  if (!xIsNan<IssmDouble>(spcdata[i])){
668  constraints->AddObject(new SpcDynamic(count+1,i+1,dof,analysis_type));
669  count++;
670  }
671  }
672  }
673  break;
674  case P1xP2Enum:
675  for(i=0;i<iomodel->numberofvertices;i++){
676  if((iomodel->my_vertices[i])){
677  if (!xIsNan<IssmDouble>(spcdata[i])){
678  constraints->AddObject(new SpcDynamic(count+1,i+1,dof,analysis_type));
679  count++;
680  }
681  }
682  }
683  for(i=0;i<iomodel->numberofedges;i++){
684  if(iomodel->edges[i*3+2]==2){
685  if(iomodel->my_edges[i]){
686  v1 = iomodel->edges[3*i+0]-1;
687  v2 = iomodel->edges[3*i+1]-1;
688  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
689  constraints->AddObject(new SpcDynamic(count+1,iomodel->numberofvertices+i+1,
690  dof,analysis_type));
691  count++;
692  }
693  }
694  }
695  }
696  break;
697  case P1xP3Enum:
698  for(i=0;i<iomodel->numberofvertices;i++){
699  if((iomodel->my_vertices[i])){
700  if (!xIsNan<IssmDouble>(spcdata[i])){
701  constraints->AddObject(new SpcDynamic(count+1,i+1,dof,analysis_type));
702  count++;
703  }
704  }
705  }
706  for(i=0;i<iomodel->numberofedges;i++){
707  if(iomodel->edges[i*3+2]==2){
708  if(iomodel->my_edges[i]){
709  v1 = iomodel->edges[3*i+0]-1;
710  v2 = iomodel->edges[3*i+1]-1;
711  if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
712  constraints->AddObject(new SpcDynamic(count+1,iomodel->numberofvertices+2*i+1,dof,analysis_type));
713  constraints->AddObject(new SpcDynamic(count+2,iomodel->numberofvertices+2*i+2,dof,analysis_type));
714  count=count+2;
715  }
716  }
717  }
718  }
719  break;
720  default:
721  _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
722  }
723  }
724  else{
725  _error_("Size of spc field not supported");
726  }
727 
728  /*Free ressources:*/
729  xDelete<IssmDouble>(times);
730  xDelete<IssmDouble>(values);
731  xDelete<bool>(boundaryedge);
732 }/*}}}*/
EdgeOnBoundaryFlags
void EdgeOnBoundaryFlags(bool **pflags, IoModel *iomodel)
Definition: CreateEdges.cpp:217
IoModel::verticalfaces
int * verticalfaces
Definition: IoModel.h:89
IssmDouble
double IssmDouble
Definition: types.h:37
TetraEnum
@ TetraEnum
Definition: EnumDefinitions.h:1300
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
IoModel::my_vfaces
bool * my_vfaces
Definition: IoModel.h:68
IoModel::my_hedges
bool * my_hedges
Definition: IoModel.h:71
IoModel::numberoffaces
int numberoffaces
Definition: IoModel.h:97
EdgesPartitioning
void EdgesPartitioning(IoModel *iomodel)
Definition: EdgesPartitioning.cpp:10
IoModel::my_elements
bool * my_elements
Definition: IoModel.h:66
P2xP1Enum
@ P2xP1Enum
Definition: EnumDefinitions.h:1226
P1bubblecondensedEnum
@ P1bubblecondensedEnum
Definition: EnumDefinitions.h:1219
SpcStatic
Definition: SpcStatic.h:13
IoModel::my_vertices
bool * my_vertices
Definition: IoModel.h:72
IoModelToDynamicConstraintsx
void IoModelToDynamicConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:31
IoModel::numberofvertices
int numberofvertices
Definition: IoModel.h:99
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
IoModel::numberofelements
int numberofelements
Definition: IoModel.h:96
IoModel::facescols
int facescols
Definition: IoModel.h:90
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
IoModel::verticaledges
int * verticaledges
Definition: IoModel.h:81
P2bubbleEnum
@ P2bubbleEnum
Definition: EnumDefinitions.h:1224
P2xP4Enum
@ P2xP4Enum
Definition: EnumDefinitions.h:1227
IoModel::horizontaledges
int * horizontaledges
Definition: IoModel.h:82
IoModel::my_edges
bool * my_edges
Definition: IoModel.h:69
IoModel::FetchData
void FetchData(bool *pboolean, const char *data_name)
Definition: IoModel.cpp:933
IoModel::numberofhorizontaledges
int numberofhorizontaledges
Definition: IoModel.h:95
P1bubbleEnum
@ P1bubbleEnum
Definition: EnumDefinitions.h:1218
SpcTransient
Definition: SpcTransient.h:13
SpcDynamic
Definition: SpcDynamic.h:13
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
IoModel::faces
int * faces
Definition: IoModel.h:88
IoModel::numberofverticalfaces
int numberofverticalfaces
Definition: IoModel.h:98
IoModel::numberofverticaledges
int numberofverticaledges
Definition: IoModel.h:94
PentaEnum
@ PentaEnum
Definition: EnumDefinitions.h:1231
IoModel::SetFilePointerToData
FILE * SetFilePointerToData(int *pcode, int *pvector_type, const char *data_name)
Definition: IoModel.cpp:2612
FacesPartitioning
void FacesPartitioning(IoModel *iomodel)
Definition: FacesPartitioning.cpp:10
P2Enum
@ P2Enum
Definition: EnumDefinitions.h:1223
IoModelToConstraintsx
void IoModelToConstraintsx(Constraints *constraints, IoModel *iomodel, const char *spc_name, int analysis_type, int finite_element, int dof)
Definition: IoModelToConstraintsx.cpp:10
IoModel::elements
int * elements
Definition: IoModel.h:79
IoModel::edges
int * edges
Definition: IoModel.h:80
IoModel::meshelementtype
int meshelementtype
Definition: IoModel.h:91
P1xP4Enum
@ P1xP4Enum
Definition: EnumDefinitions.h:1222
P1xP3Enum
@ P1xP3Enum
Definition: EnumDefinitions.h:1221
IoModel::my_faces
bool * my_faces
Definition: IoModel.h:67
TriaEnum
@ TriaEnum
Definition: EnumDefinitions.h:1318
P1xP2Enum
@ P1xP2Enum
Definition: EnumDefinitions.h:1220
IoModel::numberofedges
int numberofedges
Definition: IoModel.h:93
IoModel::my_vedges
bool * my_vedges
Definition: IoModel.h:70
Generated on Thu Jul 2 2020 08:09:24 for Ice Sheet System Model by   doxygen 1.8.19