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
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
[25486]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();
[23641]110 int localmasters = femmodel->nodes->NumberOfNodesLocal();
[16126]111 int localnumnodes = femmodel->nodes->Size();
[23587]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);
[16126]117 int numnodesperelement = femmodel->elements->MaxNumNodes();
[23587]118 int numnodesperload = femmodel->loads->MaxNumNodes();
[16126]119
[25486]120 int elementssize = femmodel->elements->Size();
121 int loadssize = femmodel->loads->Size();
[16126]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;
[25486]128 next_e = xNew<int>(elementssize*numnodesperelement);
129 count2offset_e = xNew<int>(elementssize*numnodesperelement);
[16126]130
131 k=0;
[25486]132 for(i=0;i<elementssize;i++){
[18521]133 element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[25447]134 int elementnumnodes = element->GetNumberOfNodes();
135 lidlist = xNew<int>(elementnumnodes);
[16126]136 element->GetNodesLidList(lidlist);
137
[25447]138 for(j=0;j<elementnumnodes;j++){
[16126]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 }
[25447]146 k = k + (numnodesperelement-elementnumnodes);
[16126]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;
[25486]153 next_l = xNew<int>(loadssize*numnodesperload);
154 count2offset_l = xNew<int>(loadssize*numnodesperload);
[16126]155 k=0;
[25486]156 for(i=0;i<loadssize;i++){
[18521]157 load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
[25447]158 int loadnumnodes = load->GetNumberOfNodes();
159 lidlist = xNew<int>(loadnumnodes);
[16126]160 load->GetNodesLidList(lidlist);
161
[25447]162 for(j=0;j<loadnumnodes;j++){
[16126]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 }
[25447]170 k = k + (numnodesperload-loadnumnodes);
[16126]171
172 xDelete<int>(lidlist);
173 }
174
175 /*OK now count number of dofs and flag each nodes for each node i*/
[25384]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;
[25486]181 int analysis_type;
[16126]182
[23641]183 Vector<IssmDouble> *connectivity_clone= new Vector<IssmDouble>(localmasters,numnodes);
184
[25486]185 femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
[25384]186 /*Resetting flags to false at each iteration takes a lot of time, so we keep track of the flags
[16126]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++){
[18521]193 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[25447]194 int lid = node->Lid();
195 int pid = node->Pid();
[23499]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++;
[16126]203 }
[23499]204 else{
205 break;
206 }
207 }
[25384]208 flagsindices_counter = 0;
[16126]209
[25447]210 for(j=head_e[lid];j!=-1;j=next_e[j]){
[23499]211 offset=count2offset_e[j];
212 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
[25486]213 element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum,analysis_type);
[23499]214 if(node->IsClone()){
[25447]215 connectivity_clone->SetValue(pid,d_nz+o_nz,ADD_VAL);
[16126]216 }
[23499]217 else{
[25447]218 d_connectivity[lid]+=d_nz;
219 o_connectivity[lid]+=o_nz;
[16126]220 }
221 }
[25447]222 for(j=head_l[lid];j!=-1;j=next_l[j]){
[23499]223 offset=count2offset_l[j];
224 load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
[25386]225 load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum);
[23499]226 if(node->IsClone()){
[25447]227 connectivity_clone->SetValue(pid,d_nz+o_nz,ADD_VAL);
[23499]228 }
229 else{
[25447]230 d_connectivity[lid]+=d_nz;
231 o_connectivity[lid]+=o_nz;
[23499]232 }
233 }
[16126]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*/
[23641]245 connectivity_clone->Assemble();
246 IssmDouble* serial_connectivity_clone=NULL;
[23642]247 femmodel->GetLocalVectorWithClonesNodes(&serial_connectivity_clone,connectivity_clone);
[23641]248 delete connectivity_clone;
[16126]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++){
[18521]255 Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[25447]256 int lid = node->Lid();
[23499]257 if(!node->IsClone()){
[23612]258 for(j=0;j<node->fsize;j++){
[16126]259 _assert_(count<m);
[25447]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]));
[16126]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);
[23641]275 xDelete<IssmDouble>(serial_connectivity_clone);
[16126]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.