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

Last change on this file since 25486 was 25486, checked in by bdef, 5 years ago

CHG:passing analysis_type to SetwiseNodeConnectivity to avoid overcalling FindParam
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 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}
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*/
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 int elementssize = femmodel->elements->Size();
121 int loadssize = femmodel->loads->Size();
122 /*First, we are building chaining vectors so that we know what nodes are
123 * connected to what elements. These vectors are such that:
124 * for(int i=head[id];i!=-1;i=next[i])
125 * will loop over all the elements that are connected to the node number
126 * id*/
127 head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
128 next_e = xNew<int>(elementssize*numnodesperelement);
129 count2offset_e = xNew<int>(elementssize*numnodesperelement);
130
131 k=0;
132 for(i=0;i<elementssize;i++){
133 element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
134 int elementnumnodes = element->GetNumberOfNodes();
135 lidlist = xNew<int>(elementnumnodes);
136 element->GetNodesLidList(lidlist);
137
138 for(j=0;j<elementnumnodes;j++){
139 index = lidlist[j];
140 _assert_(index>=0 && index<numnodes);
141
142 count2offset_e[k]=i;
143 next_e[k]=head_e[index];
144 head_e[index]=k++;
145 }
146 k = k + (numnodesperelement-elementnumnodes);
147
148 xDelete<int>(lidlist);
149 }
150
151 /*Chain for loads*/
152 head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
153 next_l = xNew<int>(loadssize*numnodesperload);
154 count2offset_l = xNew<int>(loadssize*numnodesperload);
155 k=0;
156 for(i=0;i<loadssize;i++){
157 load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
158 int loadnumnodes = load->GetNumberOfNodes();
159 lidlist = xNew<int>(loadnumnodes);
160 load->GetNodesLidList(lidlist);
161
162 for(j=0;j<loadnumnodes;j++){
163 index = lidlist[j];
164 _assert_(index>=0 && index<numnodes);
165
166 count2offset_l[k]=i;
167 next_l[k]=head_l[index];
168 head_l[index]=k++;
169 }
170 k = k + (numnodesperload-loadnumnodes);
171
172 xDelete<int>(lidlist);
173 }
174
175 /*OK now count number of dofs and flag each nodes for each node i*/
176 bool *flags = xNew<bool>(localnumnodes);
177 int *flagsindices = xNew<int>(localnumnodes);
178 int *d_connectivity = xNewZeroInit<int>(localnumnodes);
179 int *o_connectivity = xNewZeroInit<int>(localnumnodes);
180 int flagsindices_counter;
181 int analysis_type;
182
183 Vector<IssmDouble> *connectivity_clone= new Vector<IssmDouble>(localmasters,numnodes);
184
185 femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
186 /*Resetting flags to false at each iteration takes a lot of time, so we keep track of the flags
187 * to reset in flagsindices, initialized with -1*/
188 for(i = 0;i<localnumnodes;i++) flags[i] = false;
189 for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
190
191 /*Create connectivity vector*/
192 for(i=0;i<femmodel->nodes->Size();i++){
193 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
194 int lid = node->Lid();
195 int pid = node->Pid();
196 /*Reinitialize flags to false*/
197 j=0;
198 while(j<localnumnodes){
199 if(flagsindices[j]>=0){
200 flags[flagsindices[j]] = false;
201 flagsindices[j] = -1;
202 j++;
203 }
204 else{
205 break;
206 }
207 }
208 flagsindices_counter = 0;
209
210 for(j=head_e[lid];j!=-1;j=next_e[j]){
211 offset=count2offset_e[j];
212 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
213 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum,analysis_type);
214 if(node->IsClone()){
215 connectivity_clone->SetValue(pid,d_nz+o_nz,ADD_VAL);
216 }
217 else{
218 d_connectivity[lid]+=d_nz;
219 o_connectivity[lid]+=o_nz;
220 }
221 }
222 for(j=head_l[lid];j!=-1;j=next_l[j]){
223 offset=count2offset_l[j];
224 load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
225 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum);
226 if(node->IsClone()){
227 connectivity_clone->SetValue(pid,d_nz+o_nz,ADD_VAL);
228 }
229 else{
230 d_connectivity[lid]+=d_nz;
231 o_connectivity[lid]+=o_nz;
232 }
233 }
234 }
235 xDelete<bool>(flags);
236 xDelete<int>(flagsindices);
237 xDelete<int>(count2offset_e);
238 xDelete<int>(head_e);
239 xDelete<int>(next_e);
240 xDelete<int>(count2offset_l);
241 xDelete<int>(head_l);
242 xDelete<int>(next_l);
243
244 /*sum over all cpus*/
245 connectivity_clone->Assemble();
246 IssmDouble* serial_connectivity_clone=NULL;
247 femmodel->GetLocalVectorWithClonesNodes(&serial_connectivity_clone,connectivity_clone);
248 delete connectivity_clone;
249
250 if(set1enum==FsetEnum){
251 count=0;
252 d_nnz=xNew<int>(m);
253 o_nnz=xNew<int>(m);
254 for(i=0;i<femmodel->nodes->Size();i++){
255 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
256 int lid = node->Lid();
257 if(!node->IsClone()){
258 for(j=0;j<node->fsize;j++){
259 _assert_(count<m);
260 d_nnz[count]=numberofdofspernode*(d_connectivity[lid] + reCast<int>(serial_connectivity_clone[lid]));
261 o_nnz[count]=numberofdofspernode*(o_connectivity[lid] + reCast<int>(serial_connectivity_clone[lid]));
262 if(d_nnz[count]>n) d_nnz[count]=n;
263 if(o_nnz[count]>N-n) o_nnz[count]=N-n;
264 count++;
265 }
266 }
267 }
268 _assert_(m==count);
269 }
270 else{
271 _error_("STOP not implemented");
272 }
273 xDelete<int>(d_connectivity);
274 xDelete<int>(o_connectivity);
275 xDelete<IssmDouble>(serial_connectivity_clone);
276
277 /*Allocate ouptput pointer*/
278 *pd_nnz=d_nnz;
279 *po_nnz=o_nnz;
280}
Note: See TracBrowser for help on using the repository browser.