Ice Sheet System Model  4.18
Code documentation
Public Member Functions | Data Fields
ElementVector Class Reference

#include <ElementVector.h>

Public Member Functions

 ElementVector ()
 
 ElementVector (ElementVector *pe1, ElementVector *pe2)
 
 ElementVector (ElementVector *pe1, ElementVector *pe2, ElementVector *pe3)
 
 ElementVector (Node **nodes, int numnodes, Parameters *parameters, int approximation=NoneApproximationEnum)
 
 ~ElementVector ()
 
void AddToGlobal (Vector< IssmDouble > *pf)
 
void CheckConsistency (void)
 
void Echo (void)
 
void Init (ElementVector *pe)
 
void InsertIntoGlobal (Vector< IssmDouble > *pf)
 
void SetValue (IssmDouble scalar)
 
void StaticCondensation (ElementMatrix *Ke, int numindices, int *indices)
 

Data Fields

int nrows
 
IssmDouble * values
 
int * gglobaldoflist
 
int fsize
 
int * flocaldoflist
 
int * fglobaldoflist
 

Detailed Description

Definition at line 20 of file ElementVector.h.

Constructor & Destructor Documentation

◆ ElementVector() [1/4]

ElementVector::ElementVector ( )

Definition at line 20 of file ElementVector.cpp.

20  {/*{{{*/
21 
22  this->nrows=0;
23  this->values=NULL;
24  this->fsize=0;
25  this->flocaldoflist=NULL;
26  this->fglobaldoflist=NULL;
27 
28 }

◆ ElementVector() [2/4]

ElementVector::ElementVector ( ElementVector *  pe1,
ElementVector *  pe2 
)

Definition at line 30 of file ElementVector.cpp.

30  {/*{{{*/
31 
32  /*intermediaries*/
33  int i,j,counter;
34  int gsize,fsize;
35  int* P=NULL;
36  bool found;
37 
38  /*If one of the two matrix is NULL, we copy the other one*/
39  if(!pe1 && !pe2){
40  _error_("Two input element matrices are NULL");
41  }
42  else if(!pe1){
43  this->Init(pe2);
44  return;
45  }
46  else if(!pe2){
47  this->Init(pe1);
48  return;
49  }
50 
51  /*Initialize itransformation matrix pe[P[i]] = pe2[i]*/
52  P=xNew<int>(pe2->nrows);
53 
54  /*1: Get the new numbering of pe2 and get size of the new matrix*/
55  gsize=pe1->nrows;
56  for(i=0;i<pe2->nrows;i++){
57  found=false;
58  for(j=0;j<pe1->nrows;j++){
59  if(pe2->gglobaldoflist[i]==pe1->gglobaldoflist[j]){
60  found=true; P[i]=j; break;
61  }
62  }
63  if(!found){
64  P[i]=gsize; gsize++;
65  }
66  }
67 
68  /*2: Initialize static fields*/
69  this->nrows=gsize;
70 
71  /*Gset and values*/
72  this->gglobaldoflist=xNew<int>(this->nrows);
73  this->values=xNewZeroInit<IssmDouble>(this->nrows);
74  for(i=0;i<pe1->nrows;i++){
75  this->values[i] += pe1->values[i];
76  this->gglobaldoflist[i]=pe1->gglobaldoflist[i];
77  }
78  for(i=0;i<pe2->nrows;i++){
79  this->values[P[i]] += pe2->values[i];
80  this->gglobaldoflist[P[i]]=pe2->gglobaldoflist[i];
81  }
82 
83  /*Fset*/
84  fsize=pe1->fsize;
85  for(i=0;i<pe2->fsize;i++){
86  if(P[pe2->flocaldoflist[i]] >= pe1->nrows) fsize++;
87  }
88  this->fsize=fsize;
89  if(fsize){
90  this->flocaldoflist =xNew<int>(fsize);
91  this->fglobaldoflist=xNew<int>(fsize);
92  for(i=0;i<pe1->fsize;i++){
93  this->flocaldoflist[i] =pe1->flocaldoflist[i];
94  this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
95  }
96  counter=pe1->fsize;
97  for(i=0;i<pe2->fsize;i++){
98  if(P[pe2->flocaldoflist[i]] >= pe1->nrows){
99  this->flocaldoflist[counter] =P[pe2->flocaldoflist[i]];
100  this->fglobaldoflist[counter]=pe2->fglobaldoflist[i];
101  counter++;
102  }
103  }
104  }
105  else{
106  this->flocaldoflist=NULL;
107  this->fglobaldoflist=NULL;
108  }
109 
110  /*clean-up*/
111  xDelete<int>(P);
112 }

◆ ElementVector() [3/4]

ElementVector::ElementVector ( ElementVector *  pe1,
ElementVector *  pe2,
ElementVector *  pe3 
)

Definition at line 114 of file ElementVector.cpp.

114  {/*{{{*/
115 
116  /*Concatenate all matrices*/
117  ElementVector* pe12 =new ElementVector(pe1,pe2);
118  ElementVector* pe123=new ElementVector(pe12,pe3);
119 
120  /*Initialize current object with this matrix*/
121  this->Init(pe123);
122 
123  /*clean-up*/
124  delete pe12;
125  delete pe123;
126 }

◆ ElementVector() [4/4]

ElementVector::ElementVector ( Node **  nodes,
int  numnodes,
Parameters *  parameters,
int  approximation = NoneApproximationEnum 
)

Definition at line 128 of file ElementVector.cpp.

128  {/*{{{*/
129 
130  /*get Vector size and properties*/
131  this->nrows=GetNumberOfDofs(nodes,numnodes,GsetEnum,approximation);
132 
133  /*fill values with 0: */
134  this->values=xNewZeroInit<IssmDouble>(this->nrows);
135 
136  /*g list*/
137  this->gglobaldoflist=GetGlobalDofList(nodes,numnodes,GsetEnum,approximation);
138 
139  /*Get fsize*/
140  this->fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
141  this->flocaldoflist =GetLocalDofList( nodes,numnodes,FsetEnum,approximation);
142  this->fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
143 }

◆ ~ElementVector()

ElementVector::~ElementVector ( )

Definition at line 145 of file ElementVector.cpp.

145  {/*{{{*/
146 
147  xDelete<IssmDouble>(this->values);
148  xDelete<int>(this->gglobaldoflist);
149  xDelete<int>(this->flocaldoflist);
150  xDelete<int>(this->fglobaldoflist);
151 }

Member Function Documentation

◆ AddToGlobal()

void ElementVector::AddToGlobal ( Vector< IssmDouble > *  pf)

Definition at line 155 of file ElementVector.cpp.

155  {/*{{{*/
156 
157  int i;
158  IssmDouble* localvalues=NULL;
159 
160  /*In debugging mode, check consistency (no NaN, and values not too big)*/
161  this->CheckConsistency();
162 
163  if(this->fsize){
164  /*first, retrieve values that are in the f-set from the g-set values vector: */
165  localvalues=xNew<IssmDouble>(this->fsize);
166  for(i=0;i<this->fsize;i++){
167  localvalues[i]=this->values[this->flocaldoflist[i]];
168  }
169  /*add local values into global vector, using the fglobaldoflist: */
170  pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
171 
172  /*Free ressources:*/
173  xDelete<IssmDouble>(localvalues);
174  }
175 
176 }

◆ CheckConsistency()

void ElementVector::CheckConsistency ( void  )

Definition at line 178 of file ElementVector.cpp.

178  {/*{{{*/
179  /*Check element matrix values, only in debugging mode*/
180 #ifdef _ISSM_DEBUG_
181  for (int i=0;i<this->nrows;i++){
182  if (xIsNan<IssmDouble>(this->values[i])) _error_("NaN found in Element Vector");
183  if (xIsInf<IssmDouble>(this->values[i])) _error_("Inf found in Element Vector");
184  if (fabs( this->values[i])>1.e+50) _error_("Element Vector values exceeds 1.e+50");
185  }
186 #endif
187 }

◆ Echo()

void ElementVector::Echo ( void  )

Definition at line 189 of file ElementVector.cpp.

189  {/*{{{*/
190 
191  int i;
192 
193  _printf_("Element Vector echo:\n");
194  _printf_(" nrows: " << nrows << "\n");
195  _printf_(" values:\n");
196  for(i=0;i<nrows;i++){
197  _printf_(setw(4) << right << i << ": " << setw(10) << values[i] << "\n");
198  }
199 
200  _printf_(" gglobaldoflist (" << gglobaldoflist << "): ");
201  if(gglobaldoflist) for(i=0;i<nrows;i++) _printf_(" " << gglobaldoflist[i] );
202  _printf_(" \n");
203 
204  _printf_(" fsize: " << fsize << "\n");
205  _printf_(" flocaldoflist (" << flocaldoflist << "): ");
206  if(flocaldoflist) for(i=0;i<fsize;i++) _printf_(" " << flocaldoflist[i] );
207  _printf_(" \n");
208  _printf_(" fglobaldoflist (" << fglobaldoflist << "): ");
209  if(fglobaldoflist) for(i=0;i<fsize;i++) _printf_(" " << fglobaldoflist[i] );
210  _printf_(" \n");
211 }

◆ Init()

void ElementVector::Init ( ElementVector *  pe)

Definition at line 213 of file ElementVector.cpp.

213  {/*{{{*/
214 
215  _assert_(pe);
216 
217  this->nrows =pe->nrows;
218 
219  this->values=xNew<IssmDouble>(this->nrows);
220  xMemCpy<IssmDouble>(this->values,pe->values,this->nrows);
221 
222  this->gglobaldoflist=xNew<int>(this->nrows);
223  xMemCpy<int>(this->gglobaldoflist,pe->gglobaldoflist,this->nrows);
224 
225  this->fsize=pe->fsize;
226  if(this->fsize){
227  this->flocaldoflist=xNew<int>(this->fsize);
228  xMemCpy<int>(this->flocaldoflist,pe->flocaldoflist,this->fsize);
229  this->fglobaldoflist=xNew<int>(this->fsize);
230  xMemCpy<int>(this->fglobaldoflist,pe->fglobaldoflist,this->fsize);
231  }
232  else{
233  this->flocaldoflist=NULL;
234  this->fglobaldoflist=NULL;
235  }
236 }

◆ InsertIntoGlobal()

void ElementVector::InsertIntoGlobal ( Vector< IssmDouble > *  pf)

Definition at line 238 of file ElementVector.cpp.

238  {/*{{{*/
239 
240  int i;
241  IssmDouble* localvalues=NULL;
242 
243  if(this->fsize){
244  /*first, retrieve values that are in the f-set from the g-set values vector: */
245  localvalues=xNew<IssmDouble>(this->fsize);
246  for(i=0;i<this->fsize;i++){
247  localvalues[i]=this->values[this->flocaldoflist[i]];
248  }
249  /*add local values into global vector, using the fglobaldoflist: */
250  pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,INS_VAL);
251 
252  /*Free ressources:*/
253  xDelete<IssmDouble>(localvalues);
254  }
255 
256 }

◆ SetValue()

void ElementVector::SetValue ( IssmDouble  scalar)

Definition at line 258 of file ElementVector.cpp.

258  {/*{{{*/
259 
260  int i;
261 
262  for(i=0;i<this->nrows;i++)this->values[i]=scalar;
263 
264 }

◆ StaticCondensation()

void ElementVector::StaticCondensation ( ElementMatrix *  Ke,
int  numindices,
int *  indices 
)

Definition at line 266 of file ElementVector.cpp.

266  {/*{{{*/
267  /*
268  * | Kii Kib | | Ui | |Fi|
269  * | Kbi Kbb | | Ub | = |Fb|
270  *
271  * Kii Ui + Kib Ub = Fi
272  * Kbi Ui + Kbb Ub = Fb
273  *
274  * We want to remove Ub from the equation:
275  *
276  * Kii Ui + Kib inv(Kbb) (Fb - Kbi Ui) = Fi
277  *
278  * which gives:
279  *
280  * (Kii - Kib inv(Kbb) Kbi) Ui = Fi - Kib inv(Kbb) Fb
281  */
282 
283  /*Checks in debugging mode*/
284  _assert_(bsize>0 && bsize<this->nrows && this->values && Ke);
285  _assert_(Ke->nrows==Ke->ncols && this->nrows==Ke->nrows);
286 
287  /*Intermediaries*/
288  int counter,i,j,isize;
289  IssmDouble *Fb = NULL;
290  IssmDouble *Fi = NULL;
291  IssmDouble *Kib = NULL;
292  IssmDouble *Kbb = NULL;
293  IssmDouble *Kbbinv = NULL;
294  IssmDouble *Ftemp = NULL;
295  int *iindices = NULL;
296  bool flag;
297 
298  /*Get new sizes and indices*/
299  isize = this->nrows - bsize;
300  iindices = xNew<int>(isize);
301  counter = 0;
302  for(i=0;i<this->nrows;i++){
303  flag = true;
304  for(j=0;j<bsize;j++){
305  if(i==bindices[j]){
306  flag = false;
307  break;
308  }
309  }
310  if(flag){
311  _assert_(counter<isize);
312  iindices[counter++] = i;
313  }
314  }
315  _assert_(counter == isize);
316 
317  /*Get submatrices*/
318  Kib = xNew<IssmDouble>(isize*bsize);
319  Kbb = xNew<IssmDouble>(bsize*bsize);
320  Fb = xNew<IssmDouble>(bsize);
321  Fi = xNew<IssmDouble>(isize);
322  for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->ncols + bindices[j]];
323  for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->ncols + bindices[j]];
324  for(i=0;i<bsize;i++) Fb[i] = this->values[bindices[i]];
325  for(i=0;i<isize;i++) Fi[i] = this->values[iindices[i]];
326 
327  /*Invert Kbb*/
328  Kbbinv = xNew<IssmDouble>(bsize*bsize);
329  switch(bsize){
330  case 1:
331  Kbbinv[0] = 1./Kbb[0];
332  break;
333  case 2:
334  Matrix2x2Invert(Kbbinv,Kbb);
335  break;
336  case 3:
337  Matrix3x3Invert(Kbbinv,Kbb);
338  break;
339  default:
340  MatrixInverse(Kbbinv,bsize,bsize,NULL,0,NULL);
341  break;
342  }
343 
344  /*Calculate Kib inv(Kbb) Fb*/
345  Ftemp = xNew<IssmDouble>(isize);
346  TripleMultiply(Kib,isize,bsize,0, Kbbinv,bsize,bsize,0, Fb,bsize,1,0, Ftemp,0);
347 
348  /*New Pe*/
349  for(i=0;i<isize;i++) Ftemp[i] = Fi[i] - Ftemp[i];
350 
351  /*Update matrix values*/
352  for(i=0;i<this->nrows;i++) this->values[i]=0.;
353  for(i=0;i<isize;i++){
354  this->values[iindices[i]] = Ftemp[i];
355  }
356 
357  /*Clean up and return*/
358  xDelete<IssmDouble>(Kib);
359  xDelete<IssmDouble>(Kbb);
360  xDelete<IssmDouble>(Kbbinv);
361  xDelete<IssmDouble>(Fb);
362  xDelete<IssmDouble>(Fi);
363  xDelete<IssmDouble>(Ftemp);
364  xDelete<int>(iindices);
365  return;
366 }

Field Documentation

◆ nrows

int ElementVector::nrows

Definition at line 23 of file ElementVector.h.

◆ values

IssmDouble* ElementVector::values

Definition at line 24 of file ElementVector.h.

◆ gglobaldoflist

int* ElementVector::gglobaldoflist

Definition at line 27 of file ElementVector.h.

◆ fsize

int ElementVector::fsize

Definition at line 30 of file ElementVector.h.

◆ flocaldoflist

int* ElementVector::flocaldoflist

Definition at line 31 of file ElementVector.h.

◆ fglobaldoflist

int* ElementVector::fglobaldoflist

Definition at line 32 of file ElementVector.h.

The documentation for this class was generated from the following files:
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
MatrixInverse
int MatrixInverse(IssmDouble *a, int ndim, int nrow, IssmDouble *b, int nvec, IssmDouble *pdet)
Definition: MatrixUtils.cpp:149
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
ElementVector::nrows
int nrows
Definition: ElementVector.h:23
Matrix2x2Invert
void Matrix2x2Invert(IssmDouble *Ainv, IssmDouble *A)
Definition: MatrixUtils.cpp:329
TripleMultiply
int TripleMultiply(IssmDouble *a, int nrowa, int ncola, int itrna, IssmDouble *b, int nrowb, int ncolb, int itrnb, IssmDouble *c, int nrowc, int ncolc, int itrnc, IssmDouble *d, int iaddd)
Definition: MatrixUtils.cpp:20
ElementVector::values
IssmDouble * values
Definition: ElementVector.h:24
ElementMatrix::ncols
int ncols
Definition: ElementMatrix.h:24
GetLocalDofList
int * GetLocalDofList(Node **nodes, int numnodes, int setenum, int approximation)
Definition: Node.cpp:1072
Matrix3x3Invert
void Matrix3x3Invert(IssmDouble *Ainv, IssmDouble *A)
Definition: MatrixUtils.cpp:448
ElementVector::fglobaldoflist
int * fglobaldoflist
Definition: ElementVector.h:32
GetGlobalDofList
int * GetGlobalDofList(Node **nodes, int numnodes, int setenum, int approximation)
Definition: Node.cpp:1035
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
ElementVector::gglobaldoflist
int * gglobaldoflist
Definition: ElementVector.h:27
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
ElementVector::flocaldoflist
int * flocaldoflist
Definition: ElementVector.h:31
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
ElementVector::Init
void Init(ElementVector *pe)
Definition: ElementVector.cpp:213
ElementVector::fsize
int fsize
Definition: ElementVector.h:30
ElementVector::CheckConsistency
void CheckConsistency(void)
Definition: ElementVector.cpp:178
ElementMatrix::nrows
int nrows
Definition: ElementMatrix.h:23
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
ElementVector
Definition: ElementVector.h:20
ElementVector::ElementVector
ElementVector()
Definition: ElementVector.cpp:20
FsetEnum
@ FsetEnum
Definition: EnumDefinitions.h:1075
GetNumberOfDofs
int GetNumberOfDofs(Node **nodes, int numnodes, int setenum, int approximation)
Definition: Node.cpp:1129
ElementMatrix::values
IssmDouble * values
Definition: ElementMatrix.h:26
Generated on Thu Jul 2 2020 08:09:27 for Ice Sheet System Model by   doxygen 1.8.19