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

Last change on this file since 16126 was 16126, checked in by Mathieu Morlighem, 12 years ago

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