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

Last change on this file since 16441 was 16441, checked in by Eric.Larour, 11 years ago

CHG: fixed bug in module when allocating for ISSM matrices instead of PETSC
File size: 9.4 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 configuration_type;
15 int m,n,M,N;
16 int *d_nnz = NULL;
17 int *o_nnz = NULL;
18
19 /*output*/
20 Matrix<IssmDouble> *Kff = NULL;
21 Matrix<IssmDouble> *Kfs = NULL;
22 Vector<IssmDouble> *pf = NULL;
23 Vector<IssmDouble> *df = NULL;
24
25 bool oldalloc=false;
[16438]26 char* toolkittype=NULL;
[16126]27
28 /*retrieve parameters: */
29 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
30 femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
31
32 /*retrieve node info*/
33 fsize = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
34 ssize = femmodel->nodes->NumberOfDofs(configuration_type,SsetEnum);
35 flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
36 slocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,SsetEnum);
37
38 numberofdofspernode=femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);
39
[16438]40 /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc
[16441]41 * does, and this routine is essentially useless. Force standard alloc in this case: */
[16438]42 toolkittype=ToolkitOptions::GetToolkitType();
43
[16126]44 if(oldalloc){
45 if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
46 if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
47 if(pdf) df =new Vector<IssmDouble>(fsize);
48 if(ppf) pf =new Vector<IssmDouble>(fsize);
49 }
50 else{
51 if(pKff){
52 m=flocalsize; n=flocalsize; /*local sizes*/
53 M=fsize; N=fsize; /*global sizes*/
[16441]54 if(strcmp(toolkittype,"issm")==0){
55 Kff=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
56 }
57 else{
58 MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,FsetEnum);
59 Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
60 xDelete<int>(d_nnz);
61 xDelete<int>(o_nnz);
62 }
[16126]63 }
64 if(pKfs){
65 m=flocalsize; n=slocalsize; /*local sizes*/
66 M=fsize; N=ssize; /*global sizes*/
[16441]67 if(strcmp(toolkittype,"issm")==0){
68 Kfs=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
69 }
70 else{
71 MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,SsetEnum);
72 Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
73 xDelete<int>(d_nnz);
74 xDelete<int>(o_nnz);
75 }
[16126]76 }
77 if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
78 if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
79 }
[16438]80
81 /*Free ressources: */
82 xDelete<char>(toolkittype);
[16126]83
84 /*Allocate output pointers*/
85 if(pKff) *pKff = Kff;
86 if(pKfs) *pKfs = Kfs;
87 if(pdf) *pdf = df;
88 if(ppf) *ppf = pf;
89}
90
91void MatrixNonzeros(int** pd_nnz,int** po_nnz,FemModel* femmodel,int set1enum,int set2enum){
92
93 /*Intermediary*/
94 int i,j,k,index,offset,count;
95 int configuration_type;
96 int d_nz,o_nz;
97 Element *element = NULL;
98 Load *load = NULL;
99 int *head_e = NULL;
100 int *next_e = NULL;
101 int *count2offset_e = NULL;
102 int *head_l = NULL;
103 int *next_l = NULL;
104 int *count2offset_l = NULL;
105 int *lidlist = NULL;
106
107 /*output*/
108 int *d_nnz = NULL;
109 int *o_nnz = NULL;
110
111 /*retrive parameters: */
112 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
113
114 /*Get vector size and number of nodes*/
115 int numnodes = femmodel->nodes->NumberOfNodes(configuration_type);
116 int localnumnodes = femmodel->nodes->Size();
117 int numberofdofspernode = femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);
118 int M = femmodel->nodes->NumberOfDofs(configuration_type,set1enum);
119 int N = femmodel->nodes->NumberOfDofs(configuration_type,set2enum);
120 int m = femmodel->nodes->NumberOfDofsLocal(configuration_type,set1enum);
121 int n = femmodel->nodes->NumberOfDofsLocal(configuration_type,set2enum);
122 int numnodesperelement = femmodel->elements->MaxNumNodes();
123 int numnodesperload = femmodel->loads->MaxNumNodes(configuration_type);
124
125 /*First, we are building chaining vectors so that we know what nodes are
126 * connected to what elements. These vectors are such that:
127 * for(int i=head[id];i!=-1;i=next[i])
128 * will loop over all the elements that are connected to the node number
129 * id*/
130 head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
131 next_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
132 count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
133
134 k=0;
135 for(i=0;i<femmodel->elements->Size();i++){
136 element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
137 lidlist = xNew<int>(element->GetNumberOfNodes());
138 element->GetNodesLidList(lidlist);
139
140 for(j=0;j<element->GetNumberOfNodes();j++){
141 index = lidlist[j];
142 _assert_(index>=0 && index<numnodes);
143
144 count2offset_e[k]=i;
145 next_e[k]=head_e[index];
146 head_e[index]=k++;
147 }
148 for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
149
150 xDelete<int>(lidlist);
151 }
152
153 /*Chain for loads*/
154 head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
155 next_l = xNew<int>(femmodel->loads->Size(configuration_type)*numnodesperload);
156 count2offset_l = xNew<int>(femmodel->loads->Size(configuration_type)*numnodesperload);
157 k=0;
158 for(i=0;i<femmodel->loads->Size();i++){
159 load = dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
160 if(!load->InAnalysis(configuration_type)) continue;
161 lidlist = xNew<int>(load->GetNumberOfNodes());
162 load->GetNodesLidList(lidlist);
163
164 for(j=0;j<load->GetNumberOfNodes();j++){
165 index = lidlist[j];
166 _assert_(index>=0 && index<numnodes);
167
168 count2offset_l[k]=i;
169 next_l[k]=head_l[index];
170 head_l[index]=k++;
171 }
172 for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
173
174 xDelete<int>(lidlist);
175 }
176
177 /*OK now count number of dofs and flag each nodes for each node i*/
178 bool *flags = xNew<bool>(localnumnodes);
179 int *flagsindices = xNew<int>(localnumnodes);
180 int *d_connectivity = xNewZeroInit<int>(numnodes);
181 int *o_connectivity = xNewZeroInit<int>(numnodes);
182 int *connectivity_clone = xNewZeroInit<int>(numnodes);
183 int *all_connectivity_clone = xNewZeroInit<int>(numnodes);
184
185 /*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
186 * to reset in flagsindices, initialized with -1*/
187 for(i = 0;i<localnumnodes;i++) flags[i] = false;
188 for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
189
190 /*Create connectivity vector*/
191 for(i=0;i<femmodel->nodes->Size();i++){
192 Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
193 if(node->InAnalysis(configuration_type)){
194
195 /*Reinitialize flags to false*/
196 j=0;
197 while(true){
198 if(flagsindices[j]>=0){
199 flags[flagsindices[j]] = false;
200 flagsindices[j] = -1;
201 j++;
202 }
203 else{
204 break;
205 }
206 }
207
208 //for(j=0;j<localnumnodes;j++) flags[j]=false;
209
210 /*Loop over elements that hold node number i*/
211 //if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
212 // printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
213 //}
214 for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
215 offset=count2offset_e[j];
216 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(offset));
217 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
218 if(node->IsClone()){
219 connectivity_clone[node->Sid()]+=d_nz+o_nz;
220 }
221 else{
222 d_connectivity[node->Sid()]+=d_nz;
223 o_connectivity[node->Sid()]+=o_nz;
224 }
225 }
226 for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
227 offset=count2offset_l[j];
228 load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(offset));
229 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
230 if(node->IsClone()){
231 connectivity_clone[node->Sid()]+=d_nz+o_nz;
232 }
233 else{
234 d_connectivity[node->Sid()]+=d_nz;
235 o_connectivity[node->Sid()]+=o_nz;
236 }
237 }
238 }
239 }
240 xDelete<bool>(flags);
241 xDelete<int>(flagsindices);
242 xDelete<int>(count2offset_e);
243 xDelete<int>(head_e);
244 xDelete<int>(next_e);
245 xDelete<int>(count2offset_l);
246 xDelete<int>(head_l);
247 xDelete<int>(next_l);
248
249 /*sum over all cpus*/
250 ISSM_MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
251 xDelete<int>(connectivity_clone);
252
253 if(set1enum==FsetEnum){
254 count=0;
255 d_nnz=xNew<int>(m);
256 o_nnz=xNew<int>(m);
257 for(i=0;i<femmodel->nodes->Size();i++){
258 Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
259 if(node->InAnalysis(configuration_type) && !node->IsClone()){
260 for(j=0;j<node->indexing.fsize;j++){
261 _assert_(count<m);
262 d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
263 o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
264 if(d_nnz[count]>n) d_nnz[count]=n;
265 if(o_nnz[count]>N-n) o_nnz[count]=N-n;
266 count++;
267 }
268 }
269 }
270 _assert_(m==count);
271 }
272 else{
273 _error_("STOP not implemented");
274 }
275 xDelete<int>(d_connectivity);
276 xDelete<int>(o_connectivity);
277 xDelete<int>(all_connectivity_clone);
278
279 /*Allocate ouptput pointer*/
280 *pd_nnz=d_nnz;
281 *po_nnz=o_nnz;
282}
Note: See TracBrowser for help on using the repository browser.