source: issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp@ 23612

Last change on this file since 23612 was 23612, checked in by Mathieu Morlighem, 6 years ago

CHG: merging Dofindexing and node
File size: 8.8 KB
RevLine 
[16126]1/*!\file AllocateSystemMatricesx
2 * \brief retrieve vector from inputs in elements
3 */
4
5#include "./AllocateSystemMatricesx.h"
6#include "../../shared/shared.h"
7#include "../../toolkits/toolkits.h"
8
9void AllocateSystemMatricesx(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf,FemModel* femmodel){
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;
[16438]25 char* toolkittype=NULL;
[16126]26
27 /*retrieve parameters: */
28 femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
29
30 /*retrieve node info*/
[23587]31 fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
32 ssize = femmodel->nodes->NumberOfDofs(SsetEnum);
33 flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
34 slocalsize = femmodel->nodes->NumberOfDofsLocal(SsetEnum);
[16126]35
[23587]36 numberofdofspernode=femmodel->nodes->MaxNumDofs(GsetEnum);
[16126]37
[16438]38 /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc
[16441]39 * does, and this routine is essentially useless. Force standard alloc in this case: */
[16438]40 toolkittype=ToolkitOptions::GetToolkitType();
41
[16126]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*/
[16441]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 }
[16126]61 }
62 if(pKfs){
63 m=flocalsize; n=slocalsize; /*local sizes*/
64 M=fsize; N=ssize; /*global sizes*/
[16441]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 }
[16126]74 }
75 if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
76 if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
77 }
[23066]78
[16438]79 /*Free ressources: */
80 xDelete<char>(toolkittype);
[16126]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}
88
89void MatrixNonzeros(int** pd_nnz,int** po_nnz,FemModel* femmodel,int set1enum,int set2enum){
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*/
[23587]109 int numnodes = femmodel->nodes->NumberOfNodes();
[16126]110 int localnumnodes = femmodel->nodes->Size();
[23587]111 int numberofdofspernode = femmodel->nodes->MaxNumDofs(GsetEnum);
112 int M = femmodel->nodes->NumberOfDofs(set1enum);
113 int N = femmodel->nodes->NumberOfDofs(set2enum);
114 int m = femmodel->nodes->NumberOfDofsLocal(set1enum);
115 int n = femmodel->nodes->NumberOfDofsLocal(set2enum);
[16126]116 int numnodesperelement = femmodel->elements->MaxNumNodes();
[23587]117 int numnodesperload = femmodel->loads->MaxNumNodes();
[16126]118
119 /*First, we are building chaining vectors so that we know what nodes are
120 * connected to what elements. These vectors are such that:
121 * for(int i=head[id];i!=-1;i=next[i])
122 * will loop over all the elements that are connected to the node number
123 * id*/
124 head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
125 next_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
126 count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
127
128 k=0;
129 for(i=0;i<femmodel->elements->Size();i++){
[18521]130 element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[16126]131 lidlist = xNew<int>(element->GetNumberOfNodes());
132 element->GetNodesLidList(lidlist);
133
134 for(j=0;j<element->GetNumberOfNodes();j++){
135 index = lidlist[j];
136 _assert_(index>=0 && index<numnodes);
137
138 count2offset_e[k]=i;
139 next_e[k]=head_e[index];
140 head_e[index]=k++;
141 }
142 for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
143
144 xDelete<int>(lidlist);
145 }
146
147 /*Chain for loads*/
148 head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
[23587]149 next_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
150 count2offset_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
[16126]151 k=0;
152 for(i=0;i<femmodel->loads->Size();i++){
[18521]153 load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
[16126]154 lidlist = xNew<int>(load->GetNumberOfNodes());
155 load->GetNodesLidList(lidlist);
156
157 for(j=0;j<load->GetNumberOfNodes();j++){
158 index = lidlist[j];
159 _assert_(index>=0 && index<numnodes);
160
161 count2offset_l[k]=i;
162 next_l[k]=head_l[index];
163 head_l[index]=k++;
164 }
165 for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
166
167 xDelete<int>(lidlist);
168 }
169
170 /*OK now count number of dofs and flag each nodes for each node i*/
171 bool *flags = xNew<bool>(localnumnodes);
172 int *flagsindices = xNew<int>(localnumnodes);
173 int *d_connectivity = xNewZeroInit<int>(numnodes);
174 int *o_connectivity = xNewZeroInit<int>(numnodes);
175 int *connectivity_clone = xNewZeroInit<int>(numnodes);
176 int *all_connectivity_clone = xNewZeroInit<int>(numnodes);
177
178 /*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
179 * to reset in flagsindices, initialized with -1*/
180 for(i = 0;i<localnumnodes;i++) flags[i] = false;
181 for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
182
183 /*Create connectivity vector*/
184 for(i=0;i<femmodel->nodes->Size();i++){
[18521]185 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[16126]186
[23499]187 /*Reinitialize flags to false*/
188 j=0;
189 while(j<localnumnodes){
190 if(flagsindices[j]>=0){
191 flags[flagsindices[j]] = false;
192 flagsindices[j] = -1;
193 j++;
[16126]194 }
[23499]195 else{
196 break;
197 }
198 }
[16126]199
[23499]200 //for(j=0;j<localnumnodes;j++) flags[j]=false;
[16126]201
[23499]202 /*Loop over elements that hold node number i*/
203 //if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
204 // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
205 //}
206 for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
207 offset=count2offset_e[j];
208 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
209 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
210 if(node->IsClone()){
211 connectivity_clone[node->Sid()]+=d_nz+o_nz;
[16126]212 }
[23499]213 else{
214 d_connectivity[node->Sid()]+=d_nz;
215 o_connectivity[node->Sid()]+=o_nz;
[16126]216 }
217 }
[23499]218 for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
219 offset=count2offset_l[j];
220 load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
221 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
222 if(node->IsClone()){
223 connectivity_clone[node->Sid()]+=d_nz+o_nz;
224 }
225 else{
226 d_connectivity[node->Sid()]+=d_nz;
227 o_connectivity[node->Sid()]+=o_nz;
228 }
229 }
[16126]230 }
231 xDelete<bool>(flags);
232 xDelete<int>(flagsindices);
233 xDelete<int>(count2offset_e);
234 xDelete<int>(head_e);
235 xDelete<int>(next_e);
236 xDelete<int>(count2offset_l);
237 xDelete<int>(head_l);
238 xDelete<int>(next_l);
239
240 /*sum over all cpus*/
241 ISSM_MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
242 xDelete<int>(connectivity_clone);
243
244 if(set1enum==FsetEnum){
245 count=0;
246 d_nnz=xNew<int>(m);
247 o_nnz=xNew<int>(m);
248 for(i=0;i<femmodel->nodes->Size();i++){
[18521]249 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[23499]250 if(!node->IsClone()){
[23612]251 for(j=0;j<node->fsize;j++){
[16126]252 _assert_(count<m);
253 d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
254 o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
255 if(d_nnz[count]>n) d_nnz[count]=n;
256 if(o_nnz[count]>N-n) o_nnz[count]=N-n;
257 count++;
258 }
259 }
260 }
261 _assert_(m==count);
262 }
263 else{
264 _error_("STOP not implemented");
265 }
266 xDelete<int>(d_connectivity);
267 xDelete<int>(o_connectivity);
268 xDelete<int>(all_connectivity_clone);
269
270 /*Allocate ouptput pointer*/
271 *pd_nnz=d_nnz;
272 *po_nnz=o_nnz;
273}
Note: See TracBrowser for help on using the repository browser.