/[escript]/trunk/finley/src/NodeFile.cpp
ViewVC logotype

Annotation of /trunk/finley/src/NodeFile.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4493 - (hide annotations)
Wed Jul 3 06:32:11 2013 UTC (6 years, 2 months ago) by caltinay
File size: 40811 byte(s)
Fixed omp typo - gcc didn't pick it up...

1 jgs 150
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
7     * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10     *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 ksteube 1312
16 ksteube 1811
17 caltinay 4428 /****************************************************************************
18 jgs 82
19 caltinay 4428 Finley: NodeFile
20    
21     *****************************************************************************/
22    
23 jgs 82 #include "NodeFile.h"
24 caltinay 4428 #include <escript/Data.h>
25 caltinay 4492 #include <paso/Coupler.h>
26 jgs 82
27 caltinay 4428 #include <limits>
28     #include <sstream>
29 jgs 82
30 caltinay 4428 namespace finley {
31 ksteube 1312
32 caltinay 4428 // helper function
33     static void scatterEntries(int n, int* index, int min_index, int max_index,
34     int* Id_out, int* Id_in, int* Tag_out, int* Tag_in,
35     int* globalDegreesOfFreedom_out,
36     int* globalDegreesOfFreedom_in,
37     int numDim, double* Coordinates_out,
38     double* Coordinates_in)
39 ksteube 1312 {
40 caltinay 4428 const int range = max_index-min_index;
41     const size_t numDim_size = numDim*sizeof(double);
42 ksteube 1312
43 caltinay 4428 #pragma omp parallel for
44     for (int i=0; i<n; i++) {
45     const int k=index[i]-min_index;
46     if ((k>=0) && (k<range)) {
47     Id_out[k]=Id_in[i];
48     Tag_out[k]=Tag_in[i];
49     globalDegreesOfFreedom_out[k]=globalDegreesOfFreedom_in[i];
50     memcpy(&(Coordinates_out[INDEX2(0,k,numDim)]),
51     &(Coordinates_in[INDEX2(0,i,numDim)]), numDim_size);
52     }
53     }
54     }
55 ksteube 1312
56 caltinay 4428 // helper function
57     static void gatherEntries(int n, int* index, int min_index, int max_index,
58 caltinay 4492 int* Id_out, int* Id_in, int* Tag_out, int* Tag_in,
59 caltinay 4428 int* globalDegreesOfFreedom_out,
60 caltinay 4492 int* globalDegreesOfFreedom_in,
61 caltinay 4428 int numDim, double* Coordinates_out,
62     double* Coordinates_in)
63     {
64     const int range = max_index-min_index;
65     const size_t numDim_size = numDim*sizeof(double);
66 ksteube 1312
67 caltinay 4428 #pragma omp parallel for
68     for (int i=0; i<n; i++) {
69     const int k=index[i]-min_index;
70     if ((k>=0) && (k<range)) {
71     Id_out[i]=Id_in[k];
72     Tag_out[i]=Tag_in[k];
73     globalDegreesOfFreedom_out[i]=globalDegreesOfFreedom_in[k];
74     memcpy(&(Coordinates_out[INDEX2(0,i,numDim)]),
75     &(Coordinates_in[INDEX2(0,k,numDim)]), numDim_size);
76     }
77     }
78 jgs 82 }
79    
80 caltinay 4428 /// constructor
81     /// use NodeFile::allocTable to allocate the node table (Id,Coordinates)
82     NodeFile::NodeFile(int nDim, Esys_MPIInfo *mpiInfo) :
83     numNodes(0),
84     numDim(nDim),
85     Id(NULL),
86     Tag(NULL),
87     globalDegreesOfFreedom(NULL),
88     Coordinates(NULL),
89     globalReducedDOFIndex(NULL),
90     globalReducedNodesIndex(NULL),
91     globalNodesIndex(NULL),
92     nodesDistribution(NULL),
93     reducedNodesDistribution(NULL),
94     degreesOfFreedomDistribution(NULL),
95     reducedDegreesOfFreedomDistribution(NULL),
96     degreesOfFreedomConnector(NULL),
97     reducedDegreesOfFreedomConnector(NULL),
98     reducedNodesId(NULL),
99     degreesOfFreedomId(NULL),
100     reducedDegreesOfFreedomId(NULL),
101     status(FINLEY_INITIAL_STATUS)
102     {
103     MPIInfo = Esys_MPIInfo_getReference(mpiInfo);
104 jgs 82 }
105 ksteube 1312
106 caltinay 4428 /// destructor
107     NodeFile::~NodeFile()
108     {
109     freeTable();
110     Esys_MPIInfo_free(MPIInfo);
111 ksteube 1312 }
112    
113 caltinay 4431 /// allocates the node table within this node file to hold NN nodes.
114 caltinay 4492 void NodeFile::allocTable(int NN)
115 caltinay 4428 {
116     if (numNodes>0)
117     freeTable();
118    
119     Id=new int[NN];
120     Coordinates=new double[NN*numDim];
121     Tag=new int[NN];
122     globalDegreesOfFreedom=new int[NN];
123     globalReducedDOFIndex=new int[NN];
124     globalReducedNodesIndex=new int[NN];
125     globalNodesIndex=new int[NN];
126     reducedNodesId=new int[NN];
127     degreesOfFreedomId=new int[NN];
128     reducedDegreesOfFreedomId=new int[NN];
129     numNodes=NN;
130    
131     // this initialization makes sure that data are located on the right
132     // processor
133     #pragma omp parallel for
134     for (int n=0; n<numNodes; n++) {
135     Id[n]=-1;
136     for (int i=0; i<numDim; i++)
137     Coordinates[INDEX2(i,n,numDim)]=0.;
138     Tag[n]=-1;
139     globalDegreesOfFreedom[n]=-1;
140     globalReducedDOFIndex[n]=-1;
141     globalReducedNodesIndex[n]=-1;
142     globalNodesIndex[n]=-1;
143     reducedNodesId[n]=-1;
144     degreesOfFreedomId[n]=-1;
145 caltinay 4492 reducedDegreesOfFreedomId[n]=-1;
146 caltinay 4428 }
147 ksteube 1312 }
148    
149 caltinay 4431 /// frees the node table within this node file
150 caltinay 4428 void NodeFile::freeTable()
151     {
152     delete[] Id;
153     delete[] Coordinates;
154     delete[] globalDegreesOfFreedom;
155     delete[] globalReducedDOFIndex;
156     delete[] globalReducedNodesIndex;
157     delete[] globalNodesIndex;
158     delete[] Tag;
159     delete[] reducedNodesId;
160     delete[] degreesOfFreedomId;
161     delete[] reducedDegreesOfFreedomId;
162     tagsInUse.clear();
163 caltinay 4492 nodesMapping.clear();
164     reducedNodesMapping.clear();
165     degreesOfFreedomMapping.clear();
166     reducedDegreesOfFreedomMapping.clear();
167 caltinay 4428 Paso_Distribution_free(nodesDistribution);
168     nodesDistribution=NULL;
169     Paso_Distribution_free(reducedNodesDistribution);
170     nodesDistribution=NULL;
171     Paso_Distribution_free(degreesOfFreedomDistribution);
172     degreesOfFreedomDistribution=NULL;
173     Paso_Distribution_free(reducedDegreesOfFreedomDistribution);
174     reducedDegreesOfFreedomDistribution=NULL;
175     Paso_Connector_free(degreesOfFreedomConnector);
176     degreesOfFreedomConnector=NULL;
177     Paso_Connector_free(reducedDegreesOfFreedomConnector);
178     reducedDegreesOfFreedomConnector=NULL;
179    
180     numNodes=0;
181 ksteube 1312 }
182 caltinay 4428
183 caltinay 4492 void NodeFile::print() const
184     {
185     std::cout << "=== " << numDim << "D-Nodes:\nnumber of nodes=" << numNodes
186     << std::endl;
187     std::cout << "Id,Tag,globalDegreesOfFreedom,degreesOfFreedom,reducedDegreesOfFeedom,node,reducedNode,Coordinates" << std::endl;
188     for (int i=0; i<numNodes; i++) {
189     std::cout << Id[i] << "," << Tag[i] << "," << globalDegreesOfFreedom[i]
190     << "," << degreesOfFreedomMapping.target[i]
191     << "," << reducedDegreesOfFreedomMapping.target[i]
192     << "," << nodesMapping.target[i] << reducedNodesMapping.target[i]
193     << " ";
194     std::cout.precision(15);
195     std::cout.setf(std::ios::scientific, std::ios::floatfield);
196     for (int j=0; j<numDim; j++)
197     std:: cout << Coordinates[INDEX2(j,i,numDim)];
198     std::cout << std::endl;
199     }
200     }
201    
202 caltinay 4428 /// copies the array newX into this->coordinates
203     void NodeFile::setCoordinates(const escript::Data& cNewX)
204     {
205     if (cNewX.getDataPointSize() != numDim) {
206     std::stringstream ss;
207     ss << "NodeFile::setCoordinates: number of dimensions of new "
208     "coordinates has to be " << numDim;
209     const std::string errorMsg(ss.str());
210     Finley_setError(VALUE_ERROR, errorMsg.c_str());
211     } else if (cNewX.getNumDataPointsPerSample() != 1 ||
212     cNewX.getNumSamples() != numNodes) {
213     std::stringstream ss;
214     ss << "NodeFile::setCoordinates: number of given nodes must be "
215     << numNodes;
216     const std::string errorMsg(ss.str());
217     Finley_setError(VALUE_ERROR, errorMsg.c_str());
218     } else {
219     const size_t numDim_size=numDim*sizeof(double);
220     Finley_increaseStatus(this);
221     escript::Data& newX = *const_cast<escript::Data*>(&cNewX);
222     #pragma omp parallel for
223     for (int n=0; n<numNodes; n++) {
224     memcpy(&(Coordinates[INDEX2(0,n,numDim)]), newX.getSampleDataRO(n), numDim_size);
225     }
226     }
227 ksteube 1312 }
228    
229 caltinay 4428 /// sets tags to newTag where mask>0
230     void NodeFile::setTags(const int newTag, const escript::Data& cMask)
231     {
232     Finley_resetError();
233    
234     if (1 != cMask.getDataPointSize()) {
235     Finley_setError(TYPE_ERROR, "NodeFile::setTags: number of components of mask must be 1.");
236     return;
237     } else if (cMask.getNumDataPointsPerSample() != 1 ||
238     cMask.getNumSamples() != numNodes) {
239     Finley_setError(TYPE_ERROR, "NodeFile::setTags: illegal number of samples of mask Data object");
240     return;
241     }
242    
243     escript::Data& mask = *const_cast<escript::Data*>(&cMask);
244     #pragma omp parallel for
245     for (int n=0; n<numNodes; n++) {
246     if (mask.getSampleDataRO(n)[0] > 0)
247     Tag[n]=newTag;
248     }
249     updateTagList();
250 ksteube 1312 }
251    
252 caltinay 4428 std::pair<int,int> NodeFile::getDOFRange() const
253     {
254 caltinay 4441 std::pair<int,int> result(util::getMinMaxInt(
255     1, numNodes, globalDegreesOfFreedom));
256 caltinay 4428 if (result.second < result.first) {
257     result.first = -1;
258     result.second = 0;
259     }
260     return result;
261 ksteube 1312 }
262 caltinay 4428
263     std::pair<int,int> NodeFile::getGlobalIdRange() const
264     {
265 caltinay 4441 std::pair<int,int> result(util::getMinMaxInt(1, numNodes, Id));
266 caltinay 4428
267     #ifdef ESYS_MPI
268     int global_id_range[2];
269     int id_range[2] = { -result.first, result.second };
270     MPI_Allreduce(id_range, global_id_range, 2, MPI_INT, MPI_MAX, MPIInfo->comm);
271     result.first = -global_id_range[0];
272     result.second = global_id_range[1];
273     #endif
274     if (result.second < result.first) {
275     result.first = -1;
276     result.second = 0;
277     }
278     return result;
279 ksteube 1312 }
280    
281 caltinay 4428 std::pair<int,int> NodeFile::getGlobalDOFRange() const
282     {
283 caltinay 4441 std::pair<int,int> result(util::getMinMaxInt(
284     1, numNodes, globalDegreesOfFreedom));
285 ksteube 1312
286 caltinay 4428 #ifdef ESYS_MPI
287     int global_id_range[2];
288     int id_range[2] = { -result.first, result.second };
289     MPI_Allreduce(id_range, global_id_range, 2, MPI_INT, MPI_MAX, MPIInfo->comm);
290     result.first = -global_id_range[0];
291     result.second = global_id_range[1];
292     #endif
293     if (result.second < result.first) {
294     result.first = -1;
295     result.second = 0;
296     }
297     return result;
298 ksteube 1312 }
299 caltinay 4428
300     std::pair<int,int> NodeFile::getGlobalNodeIDIndexRange() const
301     {
302 caltinay 4441 std::pair<int,int> result(util::getMinMaxInt(1, numNodes, globalNodesIndex));
303 caltinay 4428
304     #ifdef ESYS_MPI
305     int global_id_range[2];
306     int id_range[2] = { -result.first, result.second };
307     MPI_Allreduce(id_range, global_id_range, 2, MPI_INT, MPI_MAX, MPIInfo->comm);
308     result.first = -global_id_range[0];
309     result.second = global_id_range[1];
310     #endif
311     if (result.second < result.first) {
312     result.first = -1;
313     result.second = 0;
314     }
315     return result;
316 ksteube 1312 }
317    
318 caltinay 4428 void NodeFile::copyTable(int offset, int idOffset, int dofOffset,
319     const NodeFile* in)
320     {
321     // check number of dimensions and table size
322     if (numDim != in->numDim) {
323     Finley_setError(TYPE_ERROR, "NodeFile::copyTable: dimensions of node files don't match");
324     return;
325     }
326     if (numNodes < in->numNodes+offset) {
327     Finley_setError(MEMORY_ERROR, "NodeFile::copyTable: node table is too small.");
328     return;
329     }
330 ksteube 1312
331 caltinay 4428 #pragma omp parallel for
332     for (int n=0; n<in->numNodes; n++) {
333     Id[offset+n]=in->Id[n]+idOffset;
334     Tag[offset+n]=in->Tag[n];
335     globalDegreesOfFreedom[offset+n]=in->globalDegreesOfFreedom[n]+dofOffset;
336     for(int i=0; i<numDim; i++)
337     Coordinates[INDEX2(i, offset+n, numDim)] =
338     in->Coordinates[INDEX2(i, n, in->numDim)];
339     }
340 ksteube 1312 }
341    
342 caltinay 4428 /// scatters the NodeFile in into this NodeFile using index[0:in->numNodes-1].
343     /// index has to be between 0 and numNodes-1.
344     /// colouring is chosen for the worst case
345     void NodeFile::scatter(int* index, const NodeFile* in)
346     {
347     scatterEntries(numNodes, index, 0, in->numNodes, Id, in->Id, Tag, in->Tag,
348     globalDegreesOfFreedom, in->globalDegreesOfFreedom,
349     numDim, Coordinates, in->Coordinates);
350 ksteube 1312 }
351    
352 caltinay 4428 /// gathers this NodeFile from the NodeFile 'in' using the entries in
353     /// index[0:out->numNodes-1] which are between min_index and max_index
354     /// (exclusive)
355 caltinay 4492 void NodeFile::gather(int* index, const NodeFile* in)
356 caltinay 4428 {
357     const std::pair<int,int> id_range(in->getGlobalIdRange());
358     gatherEntries(numNodes, index, id_range.first, id_range.second, Id, in->Id,
359     Tag, in->Tag, globalDegreesOfFreedom, in->globalDegreesOfFreedom,
360     numDim, Coordinates, in->Coordinates);
361 ksteube 1312 }
362    
363 caltinay 4428 void NodeFile::gather_global(int* index, const NodeFile* in)
364     {
365     // get the global range of node ids
366     const std::pair<int,int> id_range(in->getGlobalIdRange());
367     const int undefined_node=id_range.first-1;
368     std::vector<int> distribution(in->MPIInfo->size+1);
369    
370     // distribute the range of node ids
371     int buffer_len=Esys_MPIInfo_setDistribution(in->MPIInfo,
372     id_range.first, id_range.second, &distribution[0]);
373    
374     // allocate buffers
375     int *Id_buffer=new int[buffer_len];
376     int *Tag_buffer=new int[buffer_len];
377     int *globalDegreesOfFreedom_buffer=new int[buffer_len];
378     double *Coordinates_buffer=new double[buffer_len*numDim];
379    
380     // fill Id_buffer by the undefined_node marker to check if nodes
381     // are defined
382     #pragma omp parallel for
383     for (int n=0; n<buffer_len; n++)
384     Id_buffer[n]=undefined_node;
385 caltinay 4492
386 caltinay 4428 // fill the buffer by sending portions around in a circle
387     #ifdef ESYS_MPI
388     MPI_Status status;
389     int dest=Esys_MPIInfo_mod(in->MPIInfo->size, in->MPIInfo->rank+1);
390     int source=Esys_MPIInfo_mod(in->MPIInfo->size, in->MPIInfo->rank-1);
391     #endif
392     int buffer_rank=in->MPIInfo->rank;
393     for (int p=0; p<in->MPIInfo->size; ++p) {
394     if (p>0) { // the initial send can be skipped
395     #ifdef ESYS_MPI
396     MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_INT, dest,
397     in->MPIInfo->msg_tag_counter, source,
398     in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
399     MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT, dest,
400     in->MPIInfo->msg_tag_counter+1, source,
401     in->MPIInfo->msg_tag_counter+1, in->MPIInfo->comm, &status);
402     MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len,
403     MPI_INT, dest, in->MPIInfo->msg_tag_counter+2, source,
404     in->MPIInfo->msg_tag_counter+2, in->MPIInfo->comm, &status);
405     MPI_Sendrecv_replace(Coordinates_buffer, buffer_len*numDim,
406     MPI_DOUBLE, dest, in->MPIInfo->msg_tag_counter+3, source,
407     in->MPIInfo->msg_tag_counter+3, in->MPIInfo->comm, &status);
408     #endif
409     in->MPIInfo->msg_tag_counter+=4;
410     }
411     buffer_rank=Esys_MPIInfo_mod(in->MPIInfo->size, buffer_rank-1);
412     scatterEntries(in->numNodes, in->Id, distribution[buffer_rank],
413     distribution[buffer_rank+1], Id_buffer, in->Id,
414     Tag_buffer, in->Tag, globalDegreesOfFreedom_buffer,
415     in->globalDegreesOfFreedom, numDim, Coordinates_buffer,
416     in->Coordinates);
417     }
418     // now entries are collected from the buffer again by sending the
419     // entries around in a circle
420     #ifdef ESYS_MPI
421     dest=Esys_MPIInfo_mod(in->MPIInfo->size, in->MPIInfo->rank+1);
422     source=Esys_MPIInfo_mod(in->MPIInfo->size, in->MPIInfo->rank-1);
423     #endif
424     buffer_rank=in->MPIInfo->rank;
425     for (int p=0; p<in->MPIInfo->size; ++p) {
426     gatherEntries(numNodes, index, distribution[buffer_rank],
427     distribution[buffer_rank+1], Id, Id_buffer, Tag, Tag_buffer,
428     globalDegreesOfFreedom, globalDegreesOfFreedom_buffer, numDim,
429     Coordinates, Coordinates_buffer);
430     if (p < in->MPIInfo->size-1) { // the last send can be skipped
431     #ifdef ESYS_MPI
432     MPI_Sendrecv_replace(Id_buffer, buffer_len, MPI_INT, dest,
433     in->MPIInfo->msg_tag_counter, source,
434     in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
435     MPI_Sendrecv_replace(Tag_buffer, buffer_len, MPI_INT, dest,
436     in->MPIInfo->msg_tag_counter+1, source,
437     in->MPIInfo->msg_tag_counter+1, in->MPIInfo->comm, &status);
438     MPI_Sendrecv_replace(globalDegreesOfFreedom_buffer, buffer_len,
439     MPI_INT, dest, in->MPIInfo->msg_tag_counter+2, source,
440     in->MPIInfo->msg_tag_counter+2, in->MPIInfo->comm, &status);
441     MPI_Sendrecv_replace(Coordinates_buffer, buffer_len*numDim,
442     MPI_DOUBLE, dest, in->MPIInfo->msg_tag_counter+3, source,
443     in->MPIInfo->msg_tag_counter+3, in->MPIInfo->comm, &status);
444     #endif
445     in->MPIInfo->msg_tag_counter+=4;
446     }
447     buffer_rank=Esys_MPIInfo_mod(in->MPIInfo->size, buffer_rank-1);
448     }
449     // check if all nodes are set:
450     #pragma omp parallel for
451     for (int n=0; n<numNodes; ++n) {
452     if (Id[n] == undefined_node) {
453     std::stringstream ss;
454     ss << "NodeFile::gather_global: Node id " << Id[n]
455     << " at position " << n << " is referenced but not defined.";
456     const std::string errorMsg(ss.str());
457     Finley_setError(VALUE_ERROR, errorMsg.c_str());
458     }
459     }
460     delete[] Id_buffer;
461     delete[] Tag_buffer;
462     delete[] globalDegreesOfFreedom_buffer;
463     delete[] Coordinates_buffer;
464     // make sure that the error is global
465     Esys_MPIInfo_noError(in->MPIInfo);
466 ksteube 1312 }
467    
468 caltinay 4428 void NodeFile::assignMPIRankToDOFs(Esys_MPI_rank* mpiRankOfDOF,
469     int *distribution)
470     {
471     Esys_MPI_rank p_min=MPIInfo->size, p_max=-1;
472     // first we retrieve the min and max DOF on this processor to reduce
473     // costs for searching
474     const std::pair<int,int> dof_range(getDOFRange());
475    
476     for (int p=0; p<MPIInfo->size; ++p) {
477     if (distribution[p]<=dof_range.first) p_min=p;
478     if (distribution[p]<=dof_range.second) p_max=p;
479     }
480     #pragma omp parallel for
481     for (int n=0; n<numNodes; ++n) {
482     const int k=globalDegreesOfFreedom[n];
483     for (int p=p_min; p<=p_max; ++p) {
484     if (k < distribution[p+1]) {
485     mpiRankOfDOF[n]=p;
486     break;
487     }
488     }
489     }
490 caltinay 4492 }
491 caltinay 4428
492     int NodeFile::prepareLabeling(int* mask, std::vector<int>& buffer,
493     std::vector<int>& distribution, bool useNodes)
494     {
495     const int UNSET_ID=-1,SET_ID=1;
496    
497     // get the global range of DOF/node ids
498     std::pair<int,int> idRange(useNodes ?
499     getGlobalNodeIDIndexRange() : getGlobalDOFRange());
500     const int* indexArray = (useNodes ? globalNodesIndex : globalDegreesOfFreedom);
501     // distribute the range of node ids
502     distribution.assign(MPIInfo->size+1, 0);
503     int buffer_len=Esys_MPIInfo_setDistribution(MPIInfo, idRange.first,
504     idRange.second, &distribution[0]);
505     const int myCount=distribution[MPIInfo->rank+1]-distribution[MPIInfo->rank];
506    
507     // fill buffer by the UNSET_ID marker to check if nodes are defined
508     buffer.assign(buffer_len, UNSET_ID);
509    
510     // fill the buffer by sending portions around in a circle
511     #ifdef ESYS_MPI
512     MPI_Status status;
513     int dest=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank + 1);
514     int source=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank - 1);
515     #endif
516     int buffer_rank=MPIInfo->rank;
517     for (int p=0; p<MPIInfo->size; ++p) {
518     if (p>0) { // the initial send can be skipped
519     #ifdef ESYS_MPI
520     MPI_Sendrecv_replace(&buffer[0], buffer.size(), MPI_INT, dest,
521     MPIInfo->msg_tag_counter, source, MPIInfo->msg_tag_counter,
522     MPIInfo->comm, &status);
523     #endif
524     MPIInfo->msg_tag_counter++;
525     }
526     buffer_rank=Esys_MPIInfo_mod(MPIInfo->size, buffer_rank-1);
527     const int id0=distribution[buffer_rank];
528     const int id1=distribution[buffer_rank+1];
529     #pragma omp parallel for
530     for (int n=0; n<numNodes; n++) {
531     if (!mask || mask[n]>-1) {
532     const int k=indexArray[n];
533     if (id0<=k && k<id1) {
534     buffer[k-id0] = SET_ID;
535     }
536     }
537     }
538     }
539     // count the entries in the buffer
540     // TODO: OMP parallel
541     int myNewCount=0;
542     for (int n=0; n<myCount; ++n) {
543     if (buffer[n] == SET_ID) {
544     buffer[n]=myNewCount;
545     myNewCount++;
546     }
547     }
548     return myNewCount;
549 ksteube 1312 }
550    
551 caltinay 4428 int NodeFile::createDenseDOFLabeling()
552     {
553     std::vector<int> DOF_buffer;
554     std::vector<int> distribution;
555     std::vector<int> loc_offsets(MPIInfo->size);
556     std::vector<int> offsets(MPIInfo->size);
557     int new_numGlobalDOFs=0;
558    
559     // retrieve the number of own DOFs and fill buffer
560     loc_offsets[MPIInfo->rank]=prepareLabeling(NULL, DOF_buffer, distribution,
561     false);
562     #ifdef ESYS_MPI
563     MPI_Allreduce(&loc_offsets[0], &offsets[0], MPIInfo->size, MPI_INT,
564     MPI_SUM, MPIInfo->comm);
565     for (int n=0; n<MPIInfo->size; ++n) {
566     loc_offsets[n]=new_numGlobalDOFs;
567     new_numGlobalDOFs+=offsets[n];
568     }
569     #else
570     new_numGlobalDOFs=loc_offsets[0];
571     loc_offsets[0]=0;
572     #endif
573    
574     const int myDOFs=distribution[MPIInfo->rank+1]-distribution[MPIInfo->rank];
575     #pragma omp parallel for
576     for (int n=0; n<myDOFs; ++n)
577     DOF_buffer[n]+=loc_offsets[MPIInfo->rank];
578    
579 caltinay 4429 std::vector<bool_t> set_new_DOF(numNodes, TRUE);
580 caltinay 4428
581     // now entries are collected from the buffer again by sending them around
582     // in a circle
583     #ifdef ESYS_MPI
584     int dest=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank + 1);
585     int source=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank - 1);
586     #endif
587     int buffer_rank=MPIInfo->rank;
588     for (int p=0; p<MPIInfo->size; ++p) {
589     const int dof0=distribution[buffer_rank];
590     const int dof1=distribution[buffer_rank+1];
591     #pragma omp parallel for
592     for (int n=0; n<numNodes; n++) {
593 caltinay 4429 const int k=globalDegreesOfFreedom[n];
594     if (set_new_DOF[n] && dof0<=k && k<dof1) {
595     globalDegreesOfFreedom[n]=DOF_buffer[k-dof0];
596     set_new_DOF[n]=FALSE;
597     }
598 caltinay 4428 }
599     if (p<MPIInfo->size-1) { // the last send can be skipped
600     #ifdef ESYS_MPI
601     MPI_Status status;
602     MPI_Sendrecv_replace(&DOF_buffer[0], DOF_buffer.size(), MPI_INT,
603     dest, MPIInfo->msg_tag_counter, source,
604     MPIInfo->msg_tag_counter, MPIInfo->comm, &status);
605     #endif
606     MPIInfo->msg_tag_counter+=1;
607     }
608     buffer_rank=Esys_MPIInfo_mod(MPIInfo->size, buffer_rank-1);
609     }
610    
611     return new_numGlobalDOFs;
612 ksteube 1312 }
613    
614 caltinay 4428 int NodeFile::createDenseNodeLabeling(int* node_distribution,
615 caltinay 4492 const int* dof_distribution)
616 caltinay 4428 {
617     const int UNSET_ID=-1, SET_ID=1;
618     const int myFirstDOF=dof_distribution[MPIInfo->rank];
619     const int myLastDOF=dof_distribution[MPIInfo->rank+1];
620    
621     // find the range of node ids controlled by me
622     int min_id=std::numeric_limits<int>::max();
623     int max_id=std::numeric_limits<int>::min();
624     #pragma omp parallel
625     {
626     int loc_max_id=max_id;
627     int loc_min_id=min_id;
628     #pragma omp for
629     for (int n=0; n<numNodes; n++) {
630     const int dof=globalDegreesOfFreedom[n];
631     if (myFirstDOF<=dof && dof<myLastDOF) {
632     loc_max_id=std::max(loc_max_id, Id[n]);
633     loc_min_id=std::min(loc_min_id, Id[n]);
634     }
635     }
636     #pragma omp critical
637     {
638     max_id=std::max(loc_max_id, max_id);
639     min_id=std::min(loc_min_id, min_id);
640     }
641     }
642     int my_buffer_len = (max_id>=min_id ? max_id-min_id+1 : 0);
643     int buffer_len;
644    
645     #ifdef ESYS_MPI
646     MPI_Allreduce(&my_buffer_len, &buffer_len, 1, MPI_INT, MPI_MAX,
647     MPIInfo->comm);
648     #else
649     buffer_len=my_buffer_len;
650     #endif
651    
652     const int header_len=2;
653     std::vector<int> Node_buffer(buffer_len+header_len, UNSET_ID);
654     // extra storage for these IDs
655     Node_buffer[0]=min_id;
656     Node_buffer[1]=max_id;
657    
658     // mark and count the nodes in use
659 caltinay 4429 #pragma omp parallel for
660     for (int n=0; n<numNodes; n++) {
661     globalNodesIndex[n]=-1;
662     const int dof=globalDegreesOfFreedom[n];
663     if (myFirstDOF<=dof && dof<myLastDOF)
664     Node_buffer[Id[n]-min_id+header_len]=SET_ID;
665 caltinay 4428 }
666     int myNewNumNodes=0;
667     for (int n=0; n<my_buffer_len; n++) {
668     if (Node_buffer[header_len+n]==SET_ID) {
669     Node_buffer[header_len+n]=myNewNumNodes;
670     myNewNumNodes++;
671     }
672     }
673     // make the local number of nodes globally available
674     #ifdef ESYS_MPI
675     MPI_Allgather(&myNewNumNodes, 1, MPI_INT, node_distribution, 1, MPI_INT,
676     MPIInfo->comm);
677     #else
678     node_distribution[0]=myNewNumNodes;
679     #endif
680    
681     int globalNumNodes=0;
682     for (int p=0; p<MPIInfo->size; ++p) {
683     const int itmp=node_distribution[p];
684     node_distribution[p]=globalNumNodes;
685     globalNumNodes+=itmp;
686     }
687     node_distribution[MPIInfo->size]=globalNumNodes;
688    
689     // offset node buffer
690 caltinay 4429 #pragma omp parallel for
691 caltinay 4428 for (int n=0; n<my_buffer_len; n++)
692     Node_buffer[n+header_len]+=node_distribution[MPIInfo->rank];
693    
694     // now we send this buffer around to assign global node index
695     #ifdef ESYS_MPI
696     int dest=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank + 1);
697     int source=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank - 1);
698     #endif
699     int buffer_rank=MPIInfo->rank;
700     for (int p=0; p<MPIInfo->size; ++p) {
701     const int nodeID_0=Node_buffer[0];
702     const int nodeID_1=Node_buffer[1];
703     const int dof0=dof_distribution[buffer_rank];
704     const int dof1=dof_distribution[buffer_rank+1];
705     if (nodeID_0 <= nodeID_1) {
706 caltinay 4429 #pragma omp parallel for
707 caltinay 4428 for (int n=0; n<numNodes; n++) {
708     const int dof=globalDegreesOfFreedom[n];
709     const int id=Id[n]-nodeID_0;
710     if (dof0<=dof && dof<dof1 && id>=0 && id<=nodeID_1-nodeID_0)
711     globalNodesIndex[n]=Node_buffer[id+header_len];
712     }
713     }
714     if (p<MPIInfo->size-1) { // the last send can be skipped
715     #ifdef ESYS_MPI
716     MPI_Status status;
717     MPI_Sendrecv_replace(&Node_buffer[0], Node_buffer.size(), MPI_INT,
718     dest, MPIInfo->msg_tag_counter, source,
719     MPIInfo->msg_tag_counter, MPIInfo->comm, &status);
720     #endif
721     MPIInfo->msg_tag_counter+=1;
722     }
723     buffer_rank=Esys_MPIInfo_mod(MPIInfo->size, buffer_rank-1);
724     }
725     return globalNumNodes;
726 ksteube 1312 }
727    
728 caltinay 4492 int NodeFile::createDenseReducedLabeling(int* reducedMask, bool useNodes)
729 caltinay 4428 {
730     std::vector<int> buffer;
731     std::vector<int> distribution;
732     std::vector<int> loc_offsets(MPIInfo->size);
733     std::vector<int> offsets(MPIInfo->size);
734     int new_numGlobalReduced=0;
735    
736     // retrieve the number of own DOFs/nodes and fill buffer
737     loc_offsets[MPIInfo->rank]=prepareLabeling(reducedMask, buffer,
738     distribution, useNodes);
739     #ifdef ESYS_MPI
740     MPI_Allreduce(&loc_offsets[0], &offsets[0], MPIInfo->size, MPI_INT,
741     MPI_SUM, MPIInfo->comm);
742     for (int n=0; n<MPIInfo->size; ++n) {
743     loc_offsets[n]=new_numGlobalReduced;
744     new_numGlobalReduced+=offsets[n];
745     }
746     #else
747     new_numGlobalReduced=loc_offsets[0];
748     loc_offsets[0]=0;
749     #endif
750    
751     const int myCount=distribution[MPIInfo->rank+1]-distribution[MPIInfo->rank];
752     #pragma omp parallel for
753     for (int n=0; n<myCount; ++n)
754     buffer[n]+=loc_offsets[MPIInfo->rank];
755    
756     const int* denseArray =
757     (useNodes ? globalNodesIndex : globalDegreesOfFreedom);
758     int* reducedArray =
759     (useNodes ? globalReducedNodesIndex : globalReducedDOFIndex);
760    
761     #pragma omp parallel for
762     for (int n=0; n<numNodes; ++n)
763     reducedArray[n]=loc_offsets[0]-1;
764    
765     // now entries are collected from the buffer by sending them around
766     // in a circle
767     #ifdef ESYS_MPI
768     int dest=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank + 1);
769     int source=Esys_MPIInfo_mod(MPIInfo->size, MPIInfo->rank - 1);
770     #endif
771     int buffer_rank=MPIInfo->rank;
772     for (int p=0; p<MPIInfo->size; ++p) {
773     const int id0=distribution[buffer_rank];
774     const int id1=distribution[buffer_rank+1];
775     #pragma omp parallel for
776     for (int n=0; n<numNodes; n++) {
777     if (reducedMask[n] > -1) {
778     const int k=denseArray[n];
779     if (id0<=k && k<id1)
780     reducedArray[n]=buffer[k-id0];
781     }
782     }
783     if (p<MPIInfo->size-1) { // the last send can be skipped
784     #ifdef ESYS_MPI
785     MPI_Status status;
786     MPI_Sendrecv_replace(&buffer[0], buffer.size(), MPI_INT, dest,
787     MPIInfo->msg_tag_counter, source,
788     MPIInfo->msg_tag_counter, MPIInfo->comm, &status);
789     #endif
790     MPIInfo->msg_tag_counter+=1;
791     }
792     buffer_rank=Esys_MPIInfo_mod(MPIInfo->size, buffer_rank-1);
793     }
794     return new_numGlobalReduced;
795 ksteube 1312 }
796    
797 caltinay 4492 void NodeFile::createDOFMappingAndCoupling(bool use_reduced_elements)
798     {
799     Paso_Distribution* dof_distribution;
800     const int* globalDOFIndex;
801     if (use_reduced_elements) {
802     dof_distribution=reducedDegreesOfFreedomDistribution;
803     globalDOFIndex=globalReducedDOFIndex;
804     } else {
805     dof_distribution=degreesOfFreedomDistribution;
806     globalDOFIndex=globalDegreesOfFreedom;
807     }
808     const int myFirstDOF=Paso_Distribution_getFirstComponent(dof_distribution);
809     const int myLastDOF=Paso_Distribution_getLastComponent(dof_distribution);
810     const int mpiSize=MPIInfo->size;
811     const int myRank=MPIInfo->rank;
812    
813     int min_DOF, max_DOF;
814     std::pair<int,int> DOF_range(util::getFlaggedMinMaxInt(
815     numNodes, globalDOFIndex, -1));
816    
817     if (DOF_range.second < DOF_range.first) {
818     min_DOF=myFirstDOF;
819     max_DOF=myLastDOF-1;
820     } else {
821     min_DOF=DOF_range.first;
822     max_DOF=DOF_range.second;
823     }
824    
825     int p_min=mpiSize;
826     int p_max=-1;
827     if (max_DOF >= min_DOF) {
828     for (int p=0; p<mpiSize; ++p) {
829     if (dof_distribution->first_component[p]<=min_DOF) p_min=p;
830     if (dof_distribution->first_component[p]<=max_DOF) p_max=p;
831     }
832     }
833    
834     if (!((min_DOF<=myFirstDOF) && (myLastDOF-1<=max_DOF))) {
835     Finley_setError(SYSTEM_ERROR, "Local elements do not span local degrees of freedom.");
836     return;
837     }
838     const int UNUSED = -1;
839     const int len_loc_dof=max_DOF-min_DOF+1;
840     std::vector<int> shared(numNodes*(p_max-p_min+1));
841     std::vector<int> offsetInShared(mpiSize+1);
842     std::vector<int> locDOFMask(len_loc_dof, UNUSED);
843    
844     #pragma omp parallel
845     {
846     #pragma omp for
847     for (int i=0;i<numNodes;++i) {
848     const int k=globalDOFIndex[i];
849     if (k > -1) {
850     #ifdef BOUNDS_CHECK
851     if ((k-min_DOF)>=len_loc_dof) {
852     printf("BOUNDS_CHECK %s %d i=%d k=%d min_DOF=%d\n", __FILE__, __LINE__, i, k, min_DOF);
853     exit(1);
854     }
855     #endif
856     locDOFMask[k-min_DOF]=UNUSED-1;
857     }
858     }
859     #ifdef BOUNDS_CHECK
860     if (myLastDOF-min_DOF > len_loc_dof) {
861     printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);
862     exit(1);
863     }
864     #endif
865     #pragma omp for
866     for (int i=myFirstDOF-min_DOF; i<myLastDOF-min_DOF; ++i) {
867     locDOFMask[i]=i-myFirstDOF+min_DOF;
868     }
869     }
870    
871     std::vector<int> wanted_DOFs(numNodes);
872     std::vector<int> rcv_len(mpiSize);
873     std::vector<int> snd_len(mpiSize);
874     std::vector<int> neighbor(mpiSize);
875     int numNeighbors=0;
876     int n=0;
877     int lastn=n;
878     for (int p=p_min; p<=p_max; ++p) {
879     if (p != myRank) {
880     const int firstDOF=std::max(min_DOF, dof_distribution->first_component[p]);
881     const int lastDOF=std::min(max_DOF+1, dof_distribution->first_component[p+1]);
882     #ifdef BOUNDS_CHECK
883     if (firstDOF-min_DOF<0 || lastDOF-min_DOF>len_loc_dof) {
884     printf("BOUNDS_CHECK %s %d p=%d\n", __FILE__, __LINE__, p);
885     exit(1);
886     }
887     #endif
888     for (int i=firstDOF-min_DOF; i<lastDOF-min_DOF; ++i) {
889     if (locDOFMask[i] == UNUSED-1) {
890     locDOFMask[i]=myLastDOF-myFirstDOF+n;
891     wanted_DOFs[n]=i+min_DOF;
892     ++n;
893     }
894     }
895     if (n > lastn) {
896     rcv_len[p]=n-lastn;
897     #ifdef BOUNDS_CHECK
898     if (numNeighbors >= mpiSize+1) {
899     printf("BOUNDS_CHECK %s %d p=%d numNeighbors=%d n=%d\n", __FILE__, __LINE__, p, numNeighbors, n);
900     exit(1);
901     }
902     #endif
903     neighbor[numNeighbors]=p;
904     offsetInShared[numNeighbors]=lastn;
905     numNeighbors++;
906     lastn=n;
907     }
908     } // if p!=myRank
909     } // for p
910    
911     #ifdef BOUNDS_CHECK
912     if (numNeighbors >= mpiSize+1) {
913     printf("BOUNDS_CHECK %s %d numNeighbors=%d\n", __FILE__, __LINE__, numNeighbors);
914     exit(1);
915     }
916     #endif
917     offsetInShared[numNeighbors]=lastn;
918    
919     // assign new DOF labels to nodes
920     std::vector<int> nodeMask(numNodes, UNUSED);
921     #pragma omp parallel for
922     for (int i=0; i<numNodes; ++i) {
923     const int k=globalDOFIndex[i];
924     if (k > -1)
925     nodeMask[i]=locDOFMask[k-min_DOF];
926     }
927    
928     // now we can set the mapping from nodes to local DOFs
929     if (use_reduced_elements) {
930     reducedDegreesOfFreedomMapping.assign(nodeMask, UNUSED);
931     } else {
932     degreesOfFreedomMapping.assign(nodeMask, UNUSED);
933     }
934    
935     // define how to get DOF values for controlled but other processors
936     #ifdef BOUNDS_CHECK
937     if (offsetInShared[numNeighbours] >= numNodes*(p_max-p_min+1)) {
938     printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);
939     exit(1);
940     }
941     #endif
942     #pragma omp parallel for
943     for (int i=0; i<offsetInShared[numNeighbors]; ++i)
944     shared[i]=myLastDOF-myFirstDOF+i;
945    
946     Paso_SharedComponents *rcv_shcomp=Paso_SharedComponents_alloc(
947     myLastDOF-myFirstDOF, numNeighbors, &neighbor[0], &shared[0],
948     &offsetInShared[0], 1, 0, MPIInfo);
949    
950     /////////////////////////////////
951     // now we build the sender //
952     /////////////////////////////////
953     #ifdef ESYS_MPI
954     std::vector<MPI_Request> mpi_requests(mpiSize*2);
955     std::vector<MPI_Status> mpi_stati(mpiSize*2);
956     MPI_Alltoall(&rcv_len[0], 1, MPI_INT, &snd_len[0], 1, MPI_INT, MPIInfo->comm);
957     int count=0;
958     #else
959     snd_len[0]=rcv_len[0];
960     #endif
961    
962     for (int p=0; p<rcv_shcomp->numNeighbors; p++) {
963     #ifdef ESYS_MPI
964     MPI_Isend(&(wanted_DOFs[rcv_shcomp->offsetInShared[p]]),
965     rcv_shcomp->offsetInShared[p+1]-rcv_shcomp->offsetInShared[p],
966     MPI_INT, rcv_shcomp->neighbor[p],
967     MPIInfo->msg_tag_counter+myRank, MPIInfo->comm,
968     &mpi_requests[count]);
969     count++;
970     #endif
971     }
972     n=0;
973     numNeighbors=0;
974     for (int p=0; p<mpiSize; p++) {
975     if (snd_len[p] > 0) {
976     #ifdef ESYS_MPI
977     MPI_Irecv(&shared[n], snd_len[p], MPI_INT, p,
978     MPIInfo->msg_tag_counter+p, MPIInfo->comm,
979     &mpi_requests[count]);
980     count++;
981     #endif
982     neighbor[numNeighbors]=p;
983     offsetInShared[numNeighbors]=n;
984     numNeighbors++;
985     n+=snd_len[p];
986     }
987     }
988     MPIInfo->msg_tag_counter+=MPIInfo->size;
989     offsetInShared[numNeighbors]=n;
990     #ifdef ESYS_MPI
991     MPI_Waitall(count, &mpi_requests[0], &mpi_stati[0]);
992     #endif
993     // map global ids to local id's
994     #pragma omp parallel for
995     for (int i=0; i<offsetInShared[numNeighbors]; ++i) {
996     shared[i]=locDOFMask[shared[i]-min_DOF];
997     }
998    
999     Paso_SharedComponents* snd_shcomp=Paso_SharedComponents_alloc(
1000     myLastDOF-myFirstDOF, numNeighbors, &neighbor[0], &shared[0],
1001     &offsetInShared[0], 1, 0, MPIInfo);
1002    
1003     if (Finley_noError()) {
1004     if (use_reduced_elements) {
1005     reducedDegreesOfFreedomConnector=Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1006     } else {
1007     degreesOfFreedomConnector=Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
1008     }
1009     }
1010    
1011     Paso_SharedComponents_free(rcv_shcomp);
1012     Paso_SharedComponents_free(snd_shcomp);
1013     }
1014    
1015     void NodeFile::createNodeMappings(int numReducedNodes,
1016     const std::vector<int>& indexReducedNodes,
1017     const int* dofDist, const int* nodeDist)
1018     {
1019     const int mpiSize=MPIInfo->size;
1020     const int myRank=MPIInfo->rank;
1021    
1022     const int myFirstDOF=dofDist[myRank];
1023     const int myLastDOF=dofDist[myRank+1];
1024     const int myNumDOF=myLastDOF-myFirstDOF;
1025    
1026     const int myFirstNode=nodeDist[myRank];
1027     const int myLastNode=nodeDist[myRank+1];
1028     const int myNumNodes=myLastNode-myFirstNode;
1029    
1030     std::vector<int> maskMyReducedDOF(myNumDOF, -1);
1031     std::vector<int> maskMyReducedNodes(myNumNodes, -1);
1032    
1033     // mark the nodes used by the reduced mesh
1034     #pragma omp parallel for
1035     for (int i=0; i<numReducedNodes; ++i) {
1036     int k=globalNodesIndex[indexReducedNodes[i]];
1037     if (k>=myFirstNode && myLastNode>k)
1038     maskMyReducedNodes[k-myFirstNode]=i;
1039     k=globalDegreesOfFreedom[indexReducedNodes[i]];
1040     if (k>=myFirstDOF && myLastDOF>k) {
1041     maskMyReducedDOF[k-myFirstDOF]=i;
1042     }
1043     }
1044     std::vector<int> indexMyReducedDOF(myNumDOF);
1045     std::vector<int> indexMyReducedNodes(myNumNodes);
1046     int myNumReducedDOF=util::packMask(myNumDOF, &maskMyReducedDOF[0], &indexMyReducedDOF[0]);
1047     int myNumReducedNodes=util::packMask(myNumNodes, &maskMyReducedNodes[0], &indexMyReducedNodes[0]);
1048    
1049     std::vector<int> rdofDist(mpiSize+1);
1050     std::vector<int> rnodeDist(mpiSize+1);
1051     #ifdef ESYS_MPI
1052     MPI_Allgather(&myNumReducedNodes, 1, MPI_INT, &rnodeDist[0], 1, MPI_INT, MPIInfo->comm);
1053     MPI_Allgather(&myNumReducedDOF, 1, MPI_INT, &rdofDist[0], 1, MPI_INT, MPIInfo->comm);
1054     #else
1055     rnodeDist[0]=myNumReducedNodes;
1056     rdofDist[0]=myNumReducedDOF;
1057     #endif
1058     int globalNumReducedNodes=0;
1059     int globalNumReducedDOF=0;
1060     for (int i=0; i<mpiSize;++i) {
1061     int k=rnodeDist[i];
1062     rnodeDist[i]=globalNumReducedNodes;
1063     globalNumReducedNodes+=k;
1064    
1065     k=rdofDist[i];
1066     rdofDist[i]=globalNumReducedDOF;
1067     globalNumReducedDOF+=k;
1068     }
1069     rnodeDist[mpiSize]=globalNumReducedNodes;
1070     rdofDist[mpiSize]=globalNumReducedDOF;
1071    
1072     // ==== distribution of Nodes ===============================
1073     nodesDistribution=Paso_Distribution_alloc(MPIInfo, nodeDist, 1, 0);
1074     // ==== distribution of DOFs ================================
1075     degreesOfFreedomDistribution=Paso_Distribution_alloc(MPIInfo, &dofDist[0], 1,0);
1076     // ==== distribution of reduced Nodes =======================
1077     reducedNodesDistribution=Paso_Distribution_alloc(MPIInfo, &rnodeDist[0], 1, 0);
1078     // ==== distribution of reduced DOF =========================
1079     reducedDegreesOfFreedomDistribution=Paso_Distribution_alloc(MPIInfo, &rdofDist[0], 1, 0);
1080    
1081     std::vector<int> nodeMask(numNodes);
1082    
1083     if (Finley_noError()) {
1084     const int UNUSED = -1;
1085     // ==== nodes mapping which is a dummy structure ========
1086     #pragma omp parallel for
1087     for (int i=0; i<numNodes; ++i)
1088     nodeMask[i]=i;
1089     nodesMapping.assign(nodeMask, UNUSED);
1090    
1091     // ==== mapping between nodes and reduced nodes ==========
1092     #pragma omp parallel for
1093     for (int i=0; i<numNodes; ++i)
1094     nodeMask[i]=UNUSED;
1095     #pragma omp parallel for
1096     for (int i=0; i<numReducedNodes; ++i)
1097     nodeMask[indexReducedNodes[i]]=i;
1098     reducedNodesMapping.assign(nodeMask, UNUSED);
1099     }
1100     // ==== mapping between nodes and DOFs + DOF connector
1101     if (Finley_noError())
1102     createDOFMappingAndCoupling(false);
1103     // ==== mapping between nodes and reduced DOFs + reduced DOF connector
1104     if (Finley_noError())
1105     createDOFMappingAndCoupling(true);
1106    
1107     // get the Ids for DOFs and reduced nodes
1108     if (Finley_noError()) {
1109 caltinay 4493 #pragma omp parallel
1110 caltinay 4492 {
1111     #pragma omp for
1112     for (int i=0; i<reducedNodesMapping.getNumTargets(); ++i)
1113     reducedNodesId[i]=Id[reducedNodesMapping.map[i]];
1114     #pragma omp for
1115     for (int i=0; i<degreesOfFreedomMapping.getNumTargets(); ++i)
1116     degreesOfFreedomId[i]=Id[degreesOfFreedomMapping.map[i]];
1117     #pragma omp for
1118     for (int i=0; i<reducedDegreesOfFreedomMapping.getNumTargets(); ++i)
1119     reducedDegreesOfFreedomId[i]=Id[reducedDegreesOfFreedomMapping.map[i]];
1120     }
1121     } else {
1122     Paso_Distribution_free(nodesDistribution);
1123     Paso_Distribution_free(reducedNodesDistribution);
1124     Paso_Distribution_free(degreesOfFreedomDistribution);
1125     Paso_Distribution_free(reducedDegreesOfFreedomDistribution);
1126     Paso_Connector_free(degreesOfFreedomConnector);
1127     Paso_Connector_free(reducedDegreesOfFreedomConnector);
1128     nodesDistribution=NULL;
1129     reducedNodesDistribution=NULL;
1130     degreesOfFreedomDistribution=NULL;
1131     reducedDegreesOfFreedomDistribution=NULL;
1132     degreesOfFreedomConnector=NULL;
1133     reducedDegreesOfFreedomConnector=NULL;
1134     }
1135     }
1136    
1137 caltinay 4428 } // namespace finley
1138 ksteube 1312

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/NodeFile.cpp:2682-2741 /branches/pasowrap/finley/src/NodeFile.cpp:3661-3674 /branches/py3_attempt2/finley/src/NodeFile.cpp:3871-3891 /branches/restext/finley/src/NodeFile.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/NodeFile.cpp:3669-3791 /branches/stage3.0/finley/src/NodeFile.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/NodeFile.cpp:3471-3974 /release/3.0/finley/src/NodeFile.cpp:2591-2601 /trunk/finley/src/NodeFile.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26