Ice Sheet System Model  4.18
Code documentation
TriaInput2.cpp
Go to the documentation of this file.
1 
5 #ifdef HAVE_CONFIG_H
6  #include <config.h>
7 #else
8 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9 #endif
10 
11 #include "../classes.h"
12 #include "../../shared/shared.h"
13 #include "./TriaInput2.h"
14 
15 /*TriaInput2 constructors and destructor*/
17 
18  this->numberofelements_local = -1;
19  this->numberofvertices_local = -1;
20  this->isserved = false;
21  this->isserved_collapsed= 0;
22  this->M = -1;
23  this->N = -1;
24  this->values = NULL;
25  this->element_values = NULL;
26 
27 }/*}}}*/
28 TriaInput2::TriaInput2(int nbe_in,int nbv_in,int interp_in){/*{{{*/
29 
30  _assert_(nbe_in>0);
31  _assert_(nbe_in<1e11);
32  _assert_(nbv_in>0);
33  _assert_(nbv_in<1e11);
34  this->numberofelements_local = nbe_in;
35  this->numberofvertices_local = nbv_in;
36  this->isserved = false;
37  this->isserved_collapsed = 0;
38 
39  /*Reset takes care of the rest*/
40  this->Reset(interp_in);
41 }/*}}}*/
43  if(this->element_values) xDelete<IssmDouble>(this->element_values);
44  if(this->values) xDelete<IssmDouble>(this->values);
45 }
46 /*}}}*/
47 void TriaInput2::Reset(int interp_in){/*{{{*/
48 
49  /*Clean up*/
50  if(this->values) xDelete<IssmDouble>(this->values);
51  if(this->element_values) xDelete<IssmDouble>(this->element_values);
52 
53  /*Set interpolation*/
54  this->interpolation = interp_in;
55 
56  /*Create Sizes*/
57  if(this->interpolation==P1Enum){
58  this->M = this->numberofvertices_local;
59  this->N = 1;
60  }
61  else{
62  this->M = this->numberofelements_local;
63  this->N = TriaRef::NumberofNodes(interp_in);
64  }
65 
66  /*Allocate Pointers*/
67  this->values = xNewZeroInit<IssmDouble>(this->M*this->N);
68  this->element_values = xNewZeroInit<IssmDouble>(TriaRef::NumberofNodes(interp_in));
69 }/*}}}*/
70 
71 /*Object virtual functions definitions:*/
73 
75 
76  xMemCpy<IssmDouble>(output->values,this->values,this->M*this->N);
77  xMemCpy<IssmDouble>(output->element_values,this->element_values,TriaRef::NumberofNodes(this->interpolation));
78 
79  return output;
80 }
81 /*}}}*/
82 void TriaInput2::DeepEcho(void){/*{{{*/
83  _printf_("TriaInput2 Echo:\n");
84  _printf_(" interpolation: "<<EnumToStringx(this->interpolation)<<"\n");
85  _printf_(" Size: "<<M<<"x"<<N<<"\n");
86  _printf_(" isserved: "<<(isserved?"true":"false") << "\n");
87  _printf_(" isserved_collapsed: "<<isserved_collapsed << "\n");
88  if(isserved){
89  _printf_(" current values: ");
90  for(int i=0;i<3;i++) _printf_(" "<<this->element_values[i]);
91  _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
92  }
93  printarray(this->values,this->M,this->N);
94  //_printf_(setw(15)<<" TriaInput2 "<<setw(25)<<left<<EnumToStringx(this->enum_type)<<" "<<(value?"true":"false") << "\n");
95 }
96 /*}}}*/
97 void TriaInput2::Echo(void){/*{{{*/
98  _printf_("TriaInput2 Echo:\n");
99  _printf_(" interpolation: "<<EnumToStringx(this->interpolation)<<"\n");
100  _printf_(" Size: "<<M<<"x"<<N<<"\n");
101  _printf_(" isserved: "<<(isserved?"true":"false") << "\n");
102  _printf_(" isserved_collapsed: "<<isserved_collapsed << "\n");
103  if(isserved){
104  _printf_(" current values: ");
105  _printf_("[ ");
106  for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) _printf_(" "<<this->element_values[i]);
107  _printf_("] ("<<EnumToStringx(this->interpolation)<<")\n");
108  }
109 }
110 /*}}}*/
111 int TriaInput2::Id(void){/*{{{*/
112  return -1;
113 }/*}}}*/
114 void TriaInput2::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
115 
119  MARSHALLING(this->interpolation);
120  MARSHALLING(this->M);
121  MARSHALLING(this->N);
122  this->isserved = false;
123  this->isserved_collapsed = 0;
124  if(this->M*this->N){
125  MARSHALLING_DYNAMIC(this->values,IssmDouble,this->M*this->N);
126  }
127  else this->values = NULL;
128 
129  if(marshall_direction == MARSHALLING_BACKWARD){
130  this->element_values = xNewZeroInit<IssmDouble>(TriaRef::NumberofNodes(this->interpolation));
131  }
132 
133 }
134 /*}}}*/
135 int TriaInput2::ObjectEnum(void){/*{{{*/
136  return TriaInput2Enum;
137 }
138 /*}}}*/
139 
140 /*TriaInput2 management*/
141 void TriaInput2::SetInput(int interp_in,int row,IssmDouble value_in){/*{{{*/
142 
143  _assert_(this);
144  _assert_(row>=0);
145  _assert_(row<this->M);
146  _assert_(this->N==1);
147 
148  this->values[row] = value_in;
149  this->isserved = false;
150 }
151 /*}}}*/
152 void TriaInput2::SetInput(int interp_in,int numindices,int* indices,IssmDouble* values_in){/*{{{*/
153 
154  _assert_(this);
155  if(interp_in==P1Enum && this->interpolation==P1Enum){
156  _assert_(this->N==1);
157  for(int i=0;i<numindices;i++){
158  int row = indices[i];
159  _assert_(row>=0);
160  _assert_(row<this->M);
161  this->values[row] = values_in[i];
162  }
163  }
164  else if(interp_in==P0Enum && this->interpolation==P0Enum){
165  _assert_(this->N==1);
166  for(int i=0;i<numindices;i++){
167  int row = indices[i];
168  _assert_(row>=0);
169  _assert_(row<this->M);
170  this->values[row] = values_in[i];
171  }
172  }
173  else if(this->interpolation!=P1Enum && interp_in==P1Enum){
174  this->Reset(interp_in);
175  for(int i=0;i<numindices;i++){
176  int row = indices[i];
177  _assert_(row>=0);
178  _assert_(row<this->M);
179  this->values[row] = values_in[i];
180  }
181  }
182  else{
183  _error_("Cannot convert "<<EnumToStringx(this->interpolation)<<" to "<<EnumToStringx(interp_in));
184  }
185  this->isserved = false;
186 }
187 /*}}}*/
188 void TriaInput2::SetInput(int interp_in,int row,int numindices,IssmDouble* values_in){/*{{{*/
189 
190  _assert_(this);
191  if(interp_in==this->interpolation){
192  _assert_(this->N==numindices);
193  }
194  else{
195  this->Reset(interp_in);
196  _assert_(this->N==numindices);
197  }
198  for(int i=0;i<numindices;i++) this->values[row*this->N+i] = values_in[i];
199  this->isserved = false;
200 }
201 /*}}}*/
202 void TriaInput2::Serve(int numindices,int* indices){/*{{{*/
203 
204  _assert_(this);
205  _assert_(this->N==1);
206 
207  for(int i=0;i<numindices;i++){
208  int row = indices[i];
209  _assert_(row>=0);
210  _assert_(row<this->M);
211  this->element_values[i] = this->values[row];
212  }
213 
214  /*Set input as served*/
215  this->isserved = true;
216  this->isserved_collapsed = 0;
217 }
218 /*}}}*/
219 void TriaInput2::Serve(int row,int numindices){/*{{{*/
220 
221  _assert_(this);
222  _assert_(this->N==numindices);
223  _assert_(row<this->M);
224  _assert_(row>=0);
225 
226  for(int i=0;i<numindices;i++){
227  this->element_values[i] = this->values[row*this->N+i];
228  }
229 
230  /*Set input as served*/
231  this->isserved = true;
232  this->isserved_collapsed = 0;
233 } /*}}}*/
234 void TriaInput2::ServeCollapsed(int row,int id1,int id2){/*{{{*/
235 
236  _assert_(this);
237  _assert_(this->N>=3);
238  _assert_(row<this->M);
239  _assert_(row>=0);
240  _assert_(id1>=0 && id1<3);
241  _assert_(id2>=0 && id2<3);
242 
243  this->element_values[0] = this->values[row*this->N+id1];
244  this->element_values[1] = this->values[row*this->N+id2];
245 
246  /*Set input as served*/
247  this->isserved = true;
248  this->isserved_collapsed = 1;
249 }/*}}}*/
250 void TriaInput2::SetServeCollapsed(bool status){/*{{{*/
251  this->isserved_collapsed = 1;
252 }/*}}}*/
254  return this->interpolation;
255 }/*}}}*/
257  _assert_(this);
258  _assert_(this->isserved);
259 
260  int numnodes = this->NumberofNodes(this->interpolation);
261  if(this->isserved_collapsed) numnodes = 2;
262  IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
263  IssmDouble value = 0.;
264 
265  for(int i=0;i<numnodes;i++) value+=this->element_values[i];
266  value = value/numnodesd;
267 
268  *pvalue=value;
269 }/*}}}*/
271  _assert_(this);
272  _assert_(this->isserved);
273 
274  int numnodes = this->NumberofNodes(this->interpolation);
275  if(this->isserved_collapsed) numnodes = 2;
276  IssmDouble min=this->element_values[0];
277 
278  for(int i=1;i<numnodes;i++){
279  if(this->element_values[i]<min) min=this->element_values[i];
280  }
281  return min;
282 }/*}}}*/
284  _assert_(this);
285  _assert_(this->isserved);
286 
287  int numnodes = this->NumberofNodes(this->interpolation);
288  if(this->isserved_collapsed) numnodes = 2;
289  IssmDouble max=this->element_values[0];
290 
291  for(int i=1;i<numnodes;i++){
292  if(this->element_values[i]>max) max=this->element_values[i];
293  }
294  return max;
295 }/*}}}*/
297  _assert_(this);
298  _assert_(this->isserved);
299 
300  int numnodes = this->NumberofNodes(this->interpolation);
301  if(this->isserved_collapsed) numnodes = 2;
302  IssmDouble maxabs=fabs(this->element_values[0]);
303 
304  for(int i=1;i<numnodes;i++){
305  if(fabs(this->element_values[i])>maxabs) maxabs=fabs(this->element_values[i]);
306  }
307  return maxabs;
308 }/*}}}*/
309 void TriaInput2::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){/*{{{*/
310  _assert_(this);
311  _assert_(this->isserved);
312 
313  if(this->isserved_collapsed){
314  _assert_(gauss->Enum()==GaussSegEnum);
315  if(this->interpolation==P0Enum){
316  derivativevalues[0] = 0.;
317  }
318  else{
319  SegRef temp;
320  temp.GetInputDerivativeValue(derivativevalues,this->element_values,xyz_list,(GaussSeg*)gauss,P1Enum);
321  }
322  }
323  else{
324  _assert_(gauss->Enum()==GaussTriaEnum);
325  TriaRef::GetInputDerivativeValue(derivativevalues,this->element_values,xyz_list,(GaussTria*)gauss,this->interpolation);
326  }
327 }/*}}}*/
328 void TriaInput2::GetInputValue(IssmDouble* pvalue,Gauss* gauss){/*{{{*/
329  _assert_(this);
330  _assert_(this->isserved);
331  if(this->isserved_collapsed){
332  _assert_(gauss->Enum()==GaussSegEnum);
333  if(this->interpolation==P0Enum){
334  *pvalue = this->element_values[0];
335  }
336  else{
337  SegRef temp;
338  temp.GetInputValue(pvalue,this->element_values,(GaussSeg*)gauss,P1Enum);
339  }
340  }
341  else{
342  _assert_(gauss->Enum()==GaussTriaEnum);
343  TriaRef::GetInputValue(pvalue,this->element_values,(GaussTria*)gauss,this->interpolation);
344  }
345 }/*}}}*/
347  return 1;
348 }
349 /*}}}*/
351  if(this->interpolation==P0Enum || this->interpolation==P0DGEnum){
352  return P0Enum;
353  }
354  return P1Enum;
355 }/*}}}*/
357  return this->N;
358 }
359 /*}}}*/
361 
362  for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*this->values[i];
363  for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*this->element_values[i];
364 }
365 /*}}}*/
367 
368  for(int i=0;i<this->M*this->N;i++) this->values[i] = pow(this->values[i],alpha);
369  for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = pow(this->element_values[i],alpha);
370 }
371 /*}}}*/
373 
374  /*xinput is of the same type, so cast it: */
375  if(xinput->ObjectEnum()!=TriaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
376  TriaInput2* xtriainput=xDynamicCast<TriaInput2*>(xinput);
377  if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
378 
379  /*Carry out the AXPY operation depending on type:*/
380  for(int i=0;i<this->M*this->N;i++) this->values[i] = alpha*xtriainput->values[i] + this->values[i];
381  for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = alpha*xtriainput->element_values[i] + this->element_values[i];
382 }
383 /*}}}*/
384 void TriaInput2::PointWiseMult(Input2* xinput){/*{{{*/
385 
386  /*xinput is of the same type, so cast it: */
387  if(xinput->ObjectEnum()!=TriaInput2Enum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
388  TriaInput2* xtriainput=xDynamicCast<TriaInput2*>(xinput);
389  if(xtriainput->GetInterpolation()!=this->interpolation) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
390 
391  /* we need to check that the vector sizes are identical*/
392  if(xtriainput->M!=this->M||xtriainput->N!=this->N) _error_("Operation not permitted because the inputs have different sizes");
393 
394  /*Carry out the AXPY operation depending on type:*/
395  for(int i=0;i<this->M*this->N;i++) this->values[i] = xtriainput->values[i] * this->values[i];
396  for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] = xtriainput->element_values[i] * this->element_values[i];
397 }
398 /*}}}*/
399 
400 /*Object functions*/
TriaInput2::Echo
void Echo()
Definition: TriaInput2.cpp:97
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
TriaInput2::SetServeCollapsed
void SetServeCollapsed(bool)
Definition: TriaInput2.cpp:250
TriaInput2::ObjectEnum
int ObjectEnum()
Definition: TriaInput2.cpp:135
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
MARSHALLING_ENUM
#define MARSHALLING_ENUM(EN)
Definition: Marshalling.h:14
SegRef
Definition: SegRef.h:12
SegRef::GetInputValue
void GetInputValue(IssmDouble *p, IssmDouble *plist, GaussSeg *gauss, int finiteelement)
Definition: SegRef.cpp:61
P0Enum
@ P0Enum
Definition: EnumDefinitions.h:661
TriaInput2::Serve
void Serve(int numindices, int *indices)
Definition: TriaInput2.cpp:202
ElementInput2::interpolation
int interpolation
Definition: ElementInput2.h:12
TriaInput2::GetInputValue
void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: TriaInput2.cpp:328
ElementInput2::numberofelements_local
int numberofelements_local
Definition: ElementInput2.h:10
P0DGEnum
@ P0DGEnum
Definition: EnumDefinitions.h:1214
MARSHALLING_DYNAMIC
#define MARSHALLING_DYNAMIC(FIELD, TYPE, SIZE)
Definition: Marshalling.h:61
TriaInput2::GetResultNumberOfNodes
int GetResultNumberOfNodes(void)
Definition: TriaInput2.cpp:356
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
TriaInput2::copy
Input2 * copy()
Definition: TriaInput2.cpp:72
GaussTria
Definition: GaussTria.h:12
TriaInput2.h
TriaInput2::Scale
void Scale(IssmDouble scalar)
Definition: TriaInput2.cpp:360
TriaInput2::ServeCollapsed
void ServeCollapsed(int row, int id0, int in1)
Definition: TriaInput2.cpp:234
TriaInput2::GetInputAverage
void GetInputAverage(IssmDouble *pvalue)
Definition: TriaInput2.cpp:256
TriaInput2Enum
@ TriaInput2Enum
Definition: EnumDefinitions.h:1124
TriaInput2::Id
int Id()
Definition: TriaInput2.cpp:111
TriaInput2::PointWiseMult
void PointWiseMult(Input2 *xinput)
Definition: TriaInput2.cpp:384
TriaInput2::Pow
void Pow(IssmDouble scalar)
Definition: TriaInput2.cpp:366
TriaInput2::GetInputMax
IssmDouble GetInputMax()
Definition: TriaInput2.cpp:283
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
TriaInput2
Definition: TriaInput2.h:8
TriaInput2::GetInputMaxAbs
IssmDouble GetInputMaxAbs()
Definition: TriaInput2.cpp:296
TriaInput2::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: TriaInput2.cpp:114
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
MARSHALLING
#define MARSHALLING(FIELD)
Definition: Marshalling.h:29
TriaInput2::DeepEcho
void DeepEcho()
Definition: TriaInput2.cpp:82
ElementInput2::values
IssmDouble * values
Definition: ElementInput2.h:15
TriaInput2::GetResultArraySize
int GetResultArraySize(void)
Definition: TriaInput2.cpp:346
Gauss::Enum
virtual int Enum(void)=0
TriaInput2::GetInputDerivativeValue
void GetInputDerivativeValue(IssmDouble *derivativevalues, IssmDouble *xyz_list, Gauss *gauss)
Definition: TriaInput2.cpp:309
TriaInput2::GetResultInterpolation
int GetResultInterpolation(void)
Definition: TriaInput2.cpp:350
TriaInput2::GetInterpolation
int GetInterpolation()
Definition: TriaInput2.cpp:253
TriaRef::GetInputDerivativeValue
void GetInputDerivativeValue(IssmDouble *pp, IssmDouble *plist, IssmDouble *xyz_list, Gauss *gauss, int finiteelement)
Definition: TriaRef.cpp:34
SegRef::GetInputDerivativeValue
void GetInputDerivativeValue(IssmDouble *p, IssmDouble *plist, IssmDouble *xyz_list, GaussSeg *gauss, int finiteelement)
Definition: SegRef.cpp:31
TriaInput2::isserved_collapsed
int isserved_collapsed
Definition: TriaInput2.h:11
Object::ObjectEnum
virtual int ObjectEnum()=0
Input2
Definition: Input2.h:18
MARSHALLING_BACKWARD
@ MARSHALLING_BACKWARD
Definition: Marshalling.h:10
TriaInput2::GetInputMin
IssmDouble GetInputMin()
Definition: TriaInput2.cpp:270
TriaInput2::~TriaInput2
~TriaInput2()
Definition: TriaInput2.cpp:42
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
printarray
void printarray(IssmPDouble *array, int lines, int cols=1)
Definition: PrintArrays.cpp:6
GaussSegEnum
@ GaussSegEnum
Definition: EnumDefinitions.h:1079
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
ElementInput2::N
int N
Definition: ElementInput2.h:13
TriaInput2::SetInput
void SetInput(int interp_in, int row, IssmDouble value_in)
Definition: TriaInput2.cpp:141
TriaRef::GetInputValue
void GetInputValue(IssmDouble *pp, IssmDouble *plist, Gauss *gauss, int finiteelement)
Definition: TriaRef.cpp:68
GaussTriaEnum
@ GaussTriaEnum
Definition: EnumDefinitions.h:1081
TriaRef::NumberofNodes
int NumberofNodes(int finiteelement)
Definition: TriaRef.cpp:463
ElementInput2::element_values
IssmDouble * element_values
Definition: ElementInput2.h:18
ElementInput2::numberofvertices_local
int numberofvertices_local
Definition: ElementInput2.h:11
TriaInput2::Reset
void Reset(int interp_in)
Definition: TriaInput2.cpp:47
TriaInput2::AXPY
void AXPY(Input2 *xinput, IssmDouble scalar)
Definition: TriaInput2.cpp:372
ElementInput2::M
int M
Definition: ElementInput2.h:13
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
GaussSeg
Definition: GaussSeg.h:12
ElementInput2::isserved
bool isserved
Definition: ElementInput2.h:14
Gauss
Definition: Gauss.h:8
TriaInput2::TriaInput2
TriaInput2()
Definition: TriaInput2.cpp:16