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

Go to the source code of this file.

Functions

void AllocateSystemMatricesx (Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, Vector< IssmDouble > **pdf, Vector< IssmDouble > **ppf, FemModel *femmodel)
 
void MatrixNonzeros (int **pd_nnz, int **po_nnz, FemModel *femmodel, int set1enum, int set2enum)
 

Function Documentation

◆ AllocateSystemMatricesx()

void AllocateSystemMatricesx ( Matrix< IssmDouble > **  pKff,
Matrix< IssmDouble > **  pKfs,
Vector< IssmDouble > **  pdf,
Vector< IssmDouble > **  ppf,
FemModel *  femmodel 
)

Definition at line 9 of file AllocateSystemMatricesx.cpp.

9  {
10 
11  /*Intermediary*/
12  int fsize,ssize,flocalsize,slocalsize;
13  int connectivity, numberofdofspernode;
14  int m,n,M,N;
15  int *d_nnz = NULL;
16  int *o_nnz = NULL;
17 
18  /*output*/
19  Matrix<IssmDouble> *Kff = NULL;
20  Matrix<IssmDouble> *Kfs = NULL;
21  Vector<IssmDouble> *pf = NULL;
22  Vector<IssmDouble> *df = NULL;
23 
24  bool oldalloc=false;
25  char* toolkittype=NULL;
26 
27  /*retrieve parameters: */
28  femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
29 
30  /*retrieve node info*/
31  fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
32  ssize = femmodel->nodes->NumberOfDofs(SsetEnum);
33  flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
34  slocalsize = femmodel->nodes->NumberOfDofsLocal(SsetEnum);
35 
36  numberofdofspernode=femmodel->nodes->MaxNumDofs(GsetEnum);
37 
38  /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc
39  * does, and this routine is essentially useless. Force standard alloc in this case: */
40  toolkittype=ToolkitOptions::GetToolkitType();
41 
42  if(oldalloc){
43  if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
44  if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
45  if(pdf) df =new Vector<IssmDouble>(fsize);
46  if(ppf) pf =new Vector<IssmDouble>(fsize);
47  }
48  else{
49  if(pKff){
50  m=flocalsize; n=flocalsize; /*local sizes*/
51  M=fsize; N=fsize; /*global sizes*/
52  if(strcmp(toolkittype,"issm")==0){
53  Kff=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
54  }
55  else{
56  MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,FsetEnum);
57  Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
58  xDelete<int>(d_nnz);
59  xDelete<int>(o_nnz);
60  }
61  }
62  if(pKfs){
63  m=flocalsize; n=slocalsize; /*local sizes*/
64  M=fsize; N=ssize; /*global sizes*/
65  if(strcmp(toolkittype,"issm")==0){
66  Kfs=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
67  }
68  else{
69  MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,SsetEnum);
70  Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
71  xDelete<int>(d_nnz);
72  xDelete<int>(o_nnz);
73  }
74  }
75  if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
76  if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
77  }
78 
79  /*Free ressources: */
80  xDelete<char>(toolkittype);
81 
82  /*Allocate output pointers*/
83  if(pKff) *pKff = Kff;
84  if(pKfs) *pKfs = Kfs;
85  if(pdf) *pdf = df;
86  if(ppf) *ppf = pf;
87 }

◆ MatrixNonzeros()

void MatrixNonzeros ( int **  pd_nnz,
int **  po_nnz,
FemModel *  femmodel,
int  set1enum,
int  set2enum 
)

Definition at line 89 of file AllocateSystemMatricesx.cpp.

89  {
90 
91  /*Intermediary*/
92  int i,j,k,index,offset,count;
93  int d_nz,o_nz;
94  Element *element = NULL;
95  Load *load = NULL;
96  int *head_e = NULL;
97  int *next_e = NULL;
98  int *count2offset_e = NULL;
99  int *head_l = NULL;
100  int *next_l = NULL;
101  int *count2offset_l = NULL;
102  int *lidlist = NULL;
103 
104  /*output*/
105  int *d_nnz = NULL;
106  int *o_nnz = NULL;
107 
108  /*Get vector size and number of nodes*/
109  int numnodes = femmodel->nodes->NumberOfNodes();
110  int localmasters = femmodel->nodes->NumberOfNodesLocal();
111  int localnumnodes = femmodel->nodes->Size();
112  int numberofdofspernode = femmodel->nodes->MaxNumDofs(GsetEnum);
113  int M = femmodel->nodes->NumberOfDofs(set1enum);
114  int N = femmodel->nodes->NumberOfDofs(set2enum);
115  int m = femmodel->nodes->NumberOfDofsLocal(set1enum);
116  int n = femmodel->nodes->NumberOfDofsLocal(set2enum);
117  int numnodesperelement = femmodel->elements->MaxNumNodes();
118  int numnodesperload = femmodel->loads->MaxNumNodes();
119 
120  /*First, we are building chaining vectors so that we know what nodes are
121  * connected to what elements. These vectors are such that:
122  * for(int i=head[id];i!=-1;i=next[i])
123  * will loop over all the elements that are connected to the node number
124  * id*/
125  head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
126  next_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
127  count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
128 
129  k=0;
130  for(i=0;i<femmodel->elements->Size();i++){
131  element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
132  lidlist = xNew<int>(element->GetNumberOfNodes());
133  element->GetNodesLidList(lidlist);
134 
135  for(j=0;j<element->GetNumberOfNodes();j++){
136  index = lidlist[j];
137  _assert_(index>=0 && index<numnodes);
138 
139  count2offset_e[k]=i;
140  next_e[k]=head_e[index];
141  head_e[index]=k++;
142  }
143  for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
144 
145  xDelete<int>(lidlist);
146  }
147 
148  /*Chain for loads*/
149  head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
150  next_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
151  count2offset_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
152  k=0;
153  for(i=0;i<femmodel->loads->Size();i++){
154  load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
155  lidlist = xNew<int>(load->GetNumberOfNodes());
156  load->GetNodesLidList(lidlist);
157 
158  for(j=0;j<load->GetNumberOfNodes();j++){
159  index = lidlist[j];
160  _assert_(index>=0 && index<numnodes);
161 
162  count2offset_l[k]=i;
163  next_l[k]=head_l[index];
164  head_l[index]=k++;
165  }
166  for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
167 
168  xDelete<int>(lidlist);
169  }
170 
171  /*OK now count number of dofs and flag each nodes for each node i*/
172  bool *flags = xNew<bool>(localnumnodes);
173  int *flagsindices = xNew<int>(localnumnodes);
174  int *d_connectivity = xNewZeroInit<int>(localnumnodes);
175  int *o_connectivity = xNewZeroInit<int>(localnumnodes);
176 
177  Vector<IssmDouble> *connectivity_clone= new Vector<IssmDouble>(localmasters,numnodes);
178 
179  /*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
180  * to reset in flagsindices, initialized with -1*/
181  for(i = 0;i<localnumnodes;i++) flags[i] = false;
182  for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
183 
184  /*Create connectivity vector*/
185  for(i=0;i<femmodel->nodes->Size();i++){
186  Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
187 
188  /*Reinitialize flags to false*/
189  j=0;
190  while(j<localnumnodes){
191  if(flagsindices[j]>=0){
192  flags[flagsindices[j]] = false;
193  flagsindices[j] = -1;
194  j++;
195  }
196  else{
197  break;
198  }
199  }
200 
201  //for(j=0;j<localnumnodes;j++) flags[j]=false;
202 
203  /*Loop over elements that hold node number i*/
204  //if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
205  // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
206  //}
207  for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
208  offset=count2offset_e[j];
209  element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
210  element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
211  if(node->IsClone()){
212  connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
213  }
214  else{
215  d_connectivity[node->Lid()]+=d_nz;
216  o_connectivity[node->Lid()]+=o_nz;
217  }
218  }
219  for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
220  offset=count2offset_l[j];
221  load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
222  load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
223  if(node->IsClone()){
224  connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
225  }
226  else{
227  d_connectivity[node->Lid()]+=d_nz;
228  o_connectivity[node->Lid()]+=o_nz;
229  }
230  }
231  }
232  xDelete<bool>(flags);
233  xDelete<int>(flagsindices);
234  xDelete<int>(count2offset_e);
235  xDelete<int>(head_e);
236  xDelete<int>(next_e);
237  xDelete<int>(count2offset_l);
238  xDelete<int>(head_l);
239  xDelete<int>(next_l);
240 
241  /*sum over all cpus*/
242  connectivity_clone->Assemble();
243  IssmDouble* serial_connectivity_clone=NULL;
244  femmodel->GetLocalVectorWithClonesNodes(&serial_connectivity_clone,connectivity_clone);
245  delete connectivity_clone;
246 
247  if(set1enum==FsetEnum){
248  count=0;
249  d_nnz=xNew<int>(m);
250  o_nnz=xNew<int>(m);
251  for(i=0;i<femmodel->nodes->Size();i++){
252  Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
253  if(!node->IsClone()){
254  for(j=0;j<node->fsize;j++){
255  _assert_(count<m);
256  d_nnz[count]=numberofdofspernode*(d_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
257  o_nnz[count]=numberofdofspernode*(o_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
258  if(d_nnz[count]>n) d_nnz[count]=n;
259  if(o_nnz[count]>N-n) o_nnz[count]=N-n;
260  count++;
261  }
262  }
263  }
264  _assert_(m==count);
265  }
266  else{
267  _error_("STOP not implemented");
268  }
269  xDelete<int>(d_connectivity);
270  xDelete<int>(o_connectivity);
271  xDelete<IssmDouble>(serial_connectivity_clone);
272 
273  /*Allocate ouptput pointer*/
274  *pd_nnz=d_nnz;
275  *po_nnz=o_nnz;
276 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
Element::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
Nodes::MaxNumDofs
int MaxNumDofs(int setenum)
Definition: Nodes.cpp:294
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
ToolkitOptions::GetToolkitType
static char * GetToolkitType(void)
Definition: ToolkitOptions.cpp:29
MeshAverageVertexConnectivityEnum
@ MeshAverageVertexConnectivityEnum
Definition: EnumDefinitions.h:270
SsetEnum
@ SsetEnum
Definition: EnumDefinitions.h:1282
Load::SetwiseNodeConnectivity
virtual void SetwiseNodeConnectivity(int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)=0
MatrixNonzeros
void MatrixNonzeros(int **pd_nnz, int **po_nnz, FemModel *femmodel, int set1enum, int set2enum)
Definition: AllocateSystemMatricesx.cpp:89
Load
Definition: Load.h:22
Loads::MaxNumNodes
int MaxNumNodes()
Definition: Loads.cpp:132
Elements::MaxNumNodes
int MaxNumNodes(void)
Definition: Elements.cpp:46
Element
Definition: Element.h:41
Element::SetwiseNodeConnectivity
void SetwiseNodeConnectivity(int *d_nz, int *o_nz, Node *node, bool *flags, int *flagsindices, int set1_enum, int set2_enum)
Definition: Element.cpp:3369
Nodes::NumberOfNodes
int NumberOfNodes(void)
Definition: Nodes.cpp:354
Nodes::NumberOfDofsLocal
int NumberOfDofsLocal(int setenum)
Definition: Nodes.cpp:326
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
Load::GetNumberOfNodes
virtual int GetNumberOfNodes(void)=0
FemModel::GetLocalVectorWithClonesNodes
void GetLocalVectorWithClonesNodes(IssmDouble **plocal_vector, Vector< IssmDouble > *vector)
Definition: FemModel.cpp:1472
Node::IsClone
int IsClone()
Definition: Node.cpp:801
FemModel::loads
Loads * loads
Definition: FemModel.h:54
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
Node
Definition: Node.h:23
Node::Lid
int Lid(void)
Definition: Node.cpp:618
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
FsetEnum
@ FsetEnum
Definition: EnumDefinitions.h:1075
Nodes::NumberOfDofs
int NumberOfDofs(int setenum)
Definition: Nodes.cpp:314
Load::GetNodesLidList
virtual void GetNodesLidList(int *lidlist)=0
Vector< IssmDouble >
Node::Pid
int Pid(void)
Definition: Node.cpp:626
Nodes::NumberOfNodesLocal
int NumberOfNodesLocal(void)
Definition: Nodes.cpp:359
Vector::SetValue
void SetValue(int dof, doubletype value, InsMode mode)
Definition: Vector.h:163
Node::fsize
int fsize
Definition: Node.h:45
Element::GetNodesLidList
void GetNodesLidList(int *lidlist)
Definition: Element.cpp:1224
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19