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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4803 - (hide annotations)
Wed Mar 26 06:52:28 2014 UTC (5 years, 10 months ago) by caltinay
File size: 23795 byte(s)
Removed obsolete wrappers for malloc and friends.
Paso_Pattern -> paso::Pattern

1 jfenwick 3981 /*****************************************************************************
2 ksteube 1811 *
3 jfenwick 4657 * Copyright (c) 2003-2014 by University of Queensland
4 jfenwick 3981 * http://www.uq.edu.au
5 ksteube 1811 *
6     * Primary Business: Queensland, Australia
7     * Licensed under the Open Software License version 3.0
8     * http://www.opensource.org/licenses/osl-3.0.php
9     *
10 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
12     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 jfenwick 3981 *
14     *****************************************************************************/
15 ksteube 1312
16 ksteube 1811
17 caltinay 4496 /****************************************************************************
18 jgs 82
19 caltinay 4496 Finley: Mesh
20 jgs 82
21 caltinay 4496 *****************************************************************************/
22 jgs 82
23     #include "Mesh.h"
24 caltinay 4496 #include "IndexList.h"
25 jgs 82
26 caltinay 4496 namespace finley {
27 jgs 82
28 caltinay 4496 /// Constructor.
29     /// Allocates a Mesh with given name and dimensionality
30     Mesh::Mesh(const std::string name, int numDim, Esys_MPIInfo *mpi_info) :
31     m_name(name),
32     approximationOrder(-1),
33     reducedApproximationOrder(-1),
34     integrationOrder(-1),
35     reducedIntegrationOrder(-1),
36     Elements(NULL),
37     FaceElements(NULL),
38     ContactElements(NULL),
39     Points(NULL),
40     FullFullPattern(NULL),
41     FullReducedPattern(NULL),
42     ReducedFullPattern(NULL),
43     ReducedReducedPattern(NULL)
44     {
45     MPIInfo = Esys_MPIInfo_getReference(mpi_info);
46 jgs 82
47 caltinay 4496 // allocate node table
48     Nodes = new NodeFile(numDim, mpi_info);
49     }
50    
51     /// destructor
52     Mesh::~Mesh()
53 bcumming 730 {
54 caltinay 4496 delete Nodes;
55     delete FaceElements;
56     delete Elements;
57     delete ContactElements;
58     delete Points;
59     tagMap.clear();
60 caltinay 4800 paso::SystemMatrixPattern_free(FullFullPattern);
61     paso::SystemMatrixPattern_free(FullReducedPattern);
62     paso::SystemMatrixPattern_free(ReducedFullPattern);
63     paso::SystemMatrixPattern_free(ReducedReducedPattern);
64 caltinay 4496 Esys_MPIInfo_free(MPIInfo);
65     }
66    
67     void Mesh::setElements(ElementFile *elements)
68     {
69     delete Elements;
70     Elements=elements;
71     }
72    
73     void Mesh::setFaceElements(ElementFile *elements)
74     {
75     delete FaceElements;
76     FaceElements=elements;
77     }
78    
79     void Mesh::setContactElements(ElementFile *elements)
80     {
81     delete ContactElements;
82     ContactElements=elements;
83     }
84    
85     void Mesh::setPoints(ElementFile *elements)
86     {
87     delete Points;
88     Points=elements;
89     }
90    
91     void Mesh::setOrders()
92     {
93     const int ORDER_MAX=9999999;
94     int locals[4] = { ORDER_MAX, ORDER_MAX, ORDER_MAX, ORDER_MAX };
95    
96     if (Elements != NULL && Elements->numElements > 0) {
97     locals[0]=std::min(locals[0], Elements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
98     locals[1]=std::min(locals[1], Elements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
99     locals[2]=std::min(locals[2], Elements->referenceElementSet->referenceElement->integrationOrder);
100     locals[3]=std::min(locals[3], Elements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
101     }
102     if (FaceElements != NULL && FaceElements->numElements > 0) {
103     locals[0]=std::min(locals[0], FaceElements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
104     locals[1]=std::min(locals[1], FaceElements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
105     locals[2]=std::min(locals[2], FaceElements->referenceElementSet->referenceElement->integrationOrder);
106     locals[3]=std::min(locals[3], FaceElements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
107     }
108     if (ContactElements != NULL && ContactElements->numElements > 0) {
109     locals[0]=std::min(locals[0], ContactElements->referenceElementSet->referenceElement->BasisFunctions->Type->numOrder);
110     locals[1]=std::min(locals[1], ContactElements->referenceElementSet->referenceElement->LinearBasisFunctions->Type->numOrder);
111     locals[2]=std::min(locals[2], ContactElements->referenceElementSet->referenceElement->integrationOrder);
112     locals[3]=std::min(locals[3], ContactElements->referenceElementSet->referenceElementReducedQuadrature->integrationOrder);
113     }
114    
115     #ifdef ESYS_MPI
116     int globals[4];
117     MPI_Allreduce(locals, globals, 4, MPI_INT, MPI_MIN, MPIInfo->comm);
118     approximationOrder=(globals[0] < ORDER_MAX ? globals[0] : -1);
119     reducedApproximationOrder=(globals[1] < ORDER_MAX ? globals[1] : -1);
120     integrationOrder=(globals[2] < ORDER_MAX ? globals[2] : -1);
121     reducedIntegrationOrder=(globals[3] < ORDER_MAX ? globals[3] : -1);
122     #else
123     approximationOrder=(locals[0] < ORDER_MAX ? locals[0] : -1);
124     reducedApproximationOrder=(locals[1] < ORDER_MAX ? locals[1] : -1);
125     integrationOrder=(locals[2] < ORDER_MAX ? locals[2] : -1);
126     reducedIntegrationOrder=(locals[3] < ORDER_MAX ? locals[3] : -1);
127     #endif
128     }
129    
130     /// creates node mappings without (re-)distributing anything
131     void Mesh::createMappings(const std::vector<int>& dofDistribution,
132     const std::vector<int>& nodeDistribution)
133     {
134     std::vector<short> maskReducedNodes(Nodes->numNodes, -1);
135     markNodes(maskReducedNodes, 0, true);
136     std::vector<int> indexReducedNodes = util::packMask(maskReducedNodes);
137     if (noError())
138     Nodes->createNodeMappings(indexReducedNodes, dofDistribution,
139     nodeDistribution);
140     }
141    
142     /// redistributes the Nodes and Elements including overlap
143     /// according to the DOF distribution. It will create an element colouring
144     /// but will not create any mappings.
145     void Mesh::distributeByRankOfDOF(const std::vector<int>& dof_distribution)
146     {
147     std::vector<int> mpiRankOfDOF(Nodes->numNodes);
148     Nodes->assignMPIRankToDOFs(mpiRankOfDOF, dof_distribution);
149    
150     // first, the elements are redistributed according to mpiRankOfDOF
151     // at the input the Node tables refer to the local labeling of the nodes
152     // while at the output they refer to the global labeling which is rectified
153     // in the next step
154     if (noError())
155     Elements->distributeByRankOfDOF(mpiRankOfDOF, Nodes->Id);
156     if (noError())
157     FaceElements->distributeByRankOfDOF(mpiRankOfDOF, Nodes->Id);
158     if (noError())
159     ContactElements->distributeByRankOfDOF(mpiRankOfDOF, Nodes->Id);
160     if (noError())
161     Points->distributeByRankOfDOF(mpiRankOfDOF, Nodes->Id);
162    
163     // resolve the node ids
164     if (noError())
165     resolveNodeIds();
166    
167     // create a local labeling of the DOFs
168     const std::pair<int,int> dof_range(Nodes->getDOFRange());
169     const int len=dof_range.second-dof_range.first+1;
170     // local mask for used nodes
171     std::vector<int> localDOF_mask(len, -1);
172     std::vector<int> localDOF_map(Nodes->numNodes, -1);
173    
174     #pragma omp parallel for
175     for (int n=0; n<Nodes->numNodes; n++) {
176     #ifdef BOUNDS_CHECK
177     if ((Nodes->globalDegreesOfFreedom[n]-dof_range.first) >= len ||
178     (Nodes->globalDegreesOfFreedom[n]-dof_range.first) < 0) {
179     printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);
180     exit(1);
181     }
182     #endif
183     localDOF_mask[Nodes->globalDegreesOfFreedom[n]-dof_range.first]=n;
184     }
185    
186     int numDOFs=0;
187     for (int n=0; n<len; n++) {
188     const int k=localDOF_mask[n];
189     if (k>=0) {
190     localDOF_mask[n]=numDOFs;
191     numDOFs++;
192     }
193     }
194     #pragma omp parallel for
195     for (int n=0; n<Nodes->numNodes; n++) {
196     const int k=localDOF_mask[Nodes->globalDegreesOfFreedom[n]-dof_range.first];
197     localDOF_map[n]=k;
198     }
199     // create element coloring
200     if (noError())
201     createColoring(localDOF_map);
202     }
203    
204     /// prints the mesh details to standard output
205     void Mesh::print()
206     {
207     // write header
208     printf("Mesh name: %s\n", m_name.c_str());
209 jgs 82
210 caltinay 4496 // write nodes
211     Nodes->print();
212 jgs 82
213 caltinay 4496 // write elements
214     if (Elements) {
215     printf("=== %s:\nnumber of elements=%d\ncolor range=[%d,%d]\n",
216     Elements->referenceElementSet->referenceElement->Type->Name,
217     Elements->numElements, Elements->minColor, Elements->maxColor);
218     if (Elements->numElements > 0) {
219     const int NN=Elements->referenceElementSet->referenceElement->Type->numNodes;
220     const int NN2=Elements->numNodes;
221     printf("Id,Tag,Owner,Color,Nodes\n");
222     for (int i=0; i<Elements->numElements; i++) {
223     printf("%d,%d,%d,%d,", Elements->Id[i], Elements->Tag[i],
224     Elements->Owner[i], Elements->Color[i]);
225     for (int j=0; j<NN; j++)
226     printf(" %d", Nodes->Id[Elements->Nodes[INDEX2(j,i,NN2)]]);
227     printf("\n");
228     }
229     }
230     }
231 jgs 102
232 caltinay 4496 // write face elements
233     if (FaceElements) {
234     printf("=== %s:\nnumber of elements=%d\ncolor range=[%d,%d]\n",
235     FaceElements->referenceElementSet->referenceElement->Type->Name,
236     FaceElements->numElements, FaceElements->minColor,
237     FaceElements->maxColor);
238     if (FaceElements->numElements > 0) {
239     const int NN=FaceElements->referenceElementSet->referenceElement->Type->numNodes;
240     const int NN2=FaceElements->numNodes;
241     printf("Id,Tag,Owner,Color,Nodes\n");
242     for (int i=0; i<FaceElements->numElements; i++) {
243     printf("%d,%d,%d,%d,", FaceElements->Id[i],
244     FaceElements->Tag[i], FaceElements->Owner[i],
245     FaceElements->Color[i]);
246     for (int j=0; j<NN; j++)
247     printf(" %d", Nodes->Id[FaceElements->Nodes[INDEX2(j,i,NN2)]]);
248     printf("\n");
249     }
250     }
251     }
252    
253     // write Contact elements
254     if (ContactElements) {
255     printf("=== %s:\nnumber of elements=%d\ncolor range=[%d,%d]\n",
256     ContactElements->referenceElementSet->referenceElement->Type->Name,
257     ContactElements->numElements, ContactElements->minColor,
258     ContactElements->maxColor);
259     if (ContactElements->numElements > 0) {
260     const int NN=ContactElements->referenceElementSet->referenceElement->Type->numNodes;
261     const int NN2=ContactElements->numNodes;
262     printf("Id,Tag,Owner,Color,Nodes\n");
263     for (int i=0; i<ContactElements->numElements; i++) {
264     printf("%d,%d,%d,%d,", ContactElements->Id[i],
265     ContactElements->Tag[i], ContactElements->Owner[i],
266     ContactElements->Color[i]);
267     for (int j=0; j<NN; j++)
268     printf(" %d", Nodes->Id[ContactElements->Nodes[INDEX2(j,i,NN2)]]);
269     printf("\n");
270     }
271     }
272     }
273 jgs 82
274 caltinay 4496 // write points
275     if (Points) {
276     printf("=== %s:\nnumber of elements=%d\ncolor range=[%d,%d]\n",
277     Points->referenceElementSet->referenceElement->Type->Name,
278     Points->numElements, Points->minColor, Points->maxColor);
279     if (Points->numElements > 0) {
280     const int NN=Points->referenceElementSet->referenceElement->Type->numNodes;
281     const int NN2=Points->numNodes;
282     printf("Id,Tag,Owner,Color,Nodes\n");
283     for (int i=0; i<Points->numElements; i++) {
284     printf("%d,%d,%d,%d,", Points->Id[i], Points->Tag[i],
285     Points->Owner[i], Points->Color[i]);
286     for (int j=0; j<NN; j++)
287     printf(" %d", Nodes->Id[Points->Nodes[INDEX2(j,i,NN2)]]);
288     printf("\n");
289     }
290     }
291     }
292     }
293 gross 2856
294 caltinay 4496 void Mesh::markNodes(std::vector<short>& mask, int offset, bool useLinear)
295     {
296     Elements->markNodes(mask, offset, useLinear);
297     FaceElements->markNodes(mask, offset, useLinear);
298     ContactElements->markNodes(mask, offset, useLinear);
299     Points->markNodes(mask, offset, useLinear);
300 jgs 82 }
301    
302 caltinay 4496 void Mesh::markDOFsConnectedToRange(int* mask, int offset, int marker,
303     int firstDOF, int lastDOF, bool useLinear)
304     {
305     const int *dofIndex = (useLinear ? Nodes->globalReducedDOFIndex
306     : Nodes->globalDegreesOfFreedom);
307     Elements->markDOFsConnectedToRange(mask, offset, marker, firstDOF, lastDOF,
308     dofIndex, useLinear);
309     FaceElements->markDOFsConnectedToRange(mask, offset, marker, firstDOF,
310     lastDOF, dofIndex, useLinear);
311     ContactElements->markDOFsConnectedToRange(mask, offset, marker, firstDOF,
312     lastDOF, dofIndex, useLinear);
313     Points->markDOFsConnectedToRange(mask, offset, marker, firstDOF, lastDOF,
314     dofIndex, useLinear);
315 jgs 102 }
316    
317 caltinay 4496 /// optimizes the labeling of the DOFs on each processor
318     void Mesh::optimizeDOFLabeling(const std::vector<int>& distribution)
319     {
320     const int myRank=MPIInfo->rank;
321     const int mpiSize=MPIInfo->size;
322     const int myFirstVertex=distribution[myRank];
323     const int myLastVertex=distribution[myRank+1];
324     const int myNumVertices=myLastVertex-myFirstVertex;
325     int len=0;
326     for (int p=0; p<mpiSize; ++p)
327     len=std::max(len, distribution[p+1]-distribution[p]);
328 jgs 82
329 caltinay 4496 IndexList* index_list=new IndexList[myNumVertices];
330     std::vector<int> newGlobalDOFID(len);
331     // create the adjacency structure xadj and adjncy
332     #pragma omp parallel
333     {
334     // insert contributions from element matrices into columns index
335     IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list,
336     myFirstVertex, myLastVertex, Elements,
337     Nodes->globalDegreesOfFreedom,
338     Nodes->globalDegreesOfFreedom);
339     IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list,
340     myFirstVertex, myLastVertex, FaceElements,
341     Nodes->globalDegreesOfFreedom,
342     Nodes->globalDegreesOfFreedom);
343     IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list,
344     myFirstVertex, myLastVertex, ContactElements,
345     Nodes->globalDegreesOfFreedom,
346     Nodes->globalDegreesOfFreedom);
347     IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list,
348     myFirstVertex, myLastVertex, Points,
349     Nodes->globalDegreesOfFreedom,
350     Nodes->globalDegreesOfFreedom);
351     }
352     // create the local matrix pattern
353 caltinay 4803 paso::Pattern *pattern=IndexList_createPattern(0, myNumVertices,
354 caltinay 4496 index_list, myFirstVertex, myLastVertex, -myFirstVertex);
355    
356     if (noError())
357 caltinay 4803 paso::Pattern_reduceBandwidth(pattern, &newGlobalDOFID[0]);
358 caltinay 4496
359 caltinay 4803 paso::Pattern_free(pattern);
360 caltinay 4496 delete[] index_list;
361     Esys_MPIInfo_noError(MPIInfo);
362    
363     if (noError()) {
364     // shift new labeling to create a global id
365     #pragma omp parallel for
366     for (int i=0; i<myNumVertices; ++i)
367     newGlobalDOFID[i]+=myFirstVertex;
368    
369     // distribute new labeling to other processors
370     #ifdef ESYS_MPI
371     const int dest=Esys_MPIInfo_mod(mpiSize, myRank + 1);
372     const int source=Esys_MPIInfo_mod(mpiSize, myRank - 1);
373     #endif
374     int current_rank=myRank;
375     for (int p=0; p<mpiSize; ++p) {
376     const int firstVertex=distribution[current_rank];
377     const int lastVertex=distribution[current_rank+1];
378     #pragma omp parallel for
379     for (int i=0; i<Nodes->numNodes; ++i) {
380     const int k=Nodes->globalDegreesOfFreedom[i];
381     if (firstVertex<=k && k<lastVertex) {
382     Nodes->globalDegreesOfFreedom[i]=newGlobalDOFID[k-firstVertex];
383     }
384     }
385    
386     if (p<mpiSize-1) { // the final send can be skipped
387     #ifdef ESYS_MPI
388     MPI_Status status;
389     MPI_Sendrecv_replace(&newGlobalDOFID[0], len, MPI_INT,
390     dest, MPIInfo->msg_tag_counter,
391     source, MPIInfo->msg_tag_counter,
392     MPIInfo->comm, &status);
393     #endif
394     MPIInfo->msg_tag_counter++;
395     current_rank=Esys_MPIInfo_mod(mpiSize, current_rank-1);
396     }
397     }
398     }
399 jgs 82 }
400    
401 caltinay 4496 /// prepares the mesh for further use
402     void Mesh::prepare(bool optimize)
403     {
404     setOrders();
405 jgs 82
406 caltinay 4496 // first step is to distribute the elements according to a global
407     // distribution of DOF
408     std::vector<int> distribution(MPIInfo->size+1);
409 jgs 82
410 caltinay 4496 // first we create dense labeling for the DOFs
411     int newGlobalNumDOFs=Nodes->createDenseDOFLabeling();
412    
413     // create a distribution of the global DOFs and determine the MPI rank
414     // controlling the DOFs on this processor
415     Esys_MPIInfo_setDistribution(MPIInfo, 0, newGlobalNumDOFs-1, &distribution[0]);
416    
417     // now the mesh is re-distributed according to the distribution vector
418     // this will redistribute the Nodes and Elements including overlap and
419     // will create an element coloring but will not create any mappings
420     // (see later in this function)
421     if (noError())
422     distributeByRankOfDOF(distribution);
423    
424     // at this stage we are able to start an optimization of the DOF
425     // distribution using ParMetis. On return distribution is altered and
426     // new DOF IDs have been assigned
427     if (noError() && optimize && MPIInfo->size>1) {
428     optimizeDOFDistribution(distribution);
429     if (noError())
430     distributeByRankOfDOF(distribution);
431     }
432     // the local labelling of the degrees of freedom is optimized
433     if (noError() && optimize) {
434     optimizeDOFLabeling(distribution);
435     }
436     // rearrange elements with the aim of bringing elements closer to memory
437     // locations of the nodes (distributed shared memory!):
438     optimizeElementOrdering();
439    
440     // create the global indices
441     if (noError()) {
442     std::vector<short> maskReducedNodes(Nodes->numNodes, -1);
443     std::vector<int> nodeDistribution(MPIInfo->size+1);
444     markNodes(maskReducedNodes, 0, true);
445     std::vector<int> indexReducedNodes = util::packMask(maskReducedNodes);
446    
447     Nodes->createDenseNodeLabeling(nodeDistribution, distribution);
448     // created reduced DOF labeling
449     Nodes->createDenseReducedLabeling(maskReducedNodes, false);
450     // created reduced node labeling
451     Nodes->createDenseReducedLabeling(maskReducedNodes, true);
452    
453     // create the missing mappings
454     if (noError())
455     Nodes->createNodeMappings(indexReducedNodes, distribution, nodeDistribution);
456     }
457    
458     updateTagList();
459 jgs 82 }
460    
461 caltinay 4496 /// tries to reduce the number of colours for all element files
462     void Mesh::createColoring(const std::vector<int>& dofMap)
463     {
464     if (noError())
465     Elements->createColoring(dofMap);
466     if (noError())
467     FaceElements->createColoring(dofMap);
468     if (noError())
469     Points->createColoring(dofMap);
470     if (noError())
471     ContactElements->createColoring(dofMap);
472 jgs 82 }
473 caltinay 4496
474     /// redistributes elements to minimize communication during assemblage
475     void Mesh::optimizeElementOrdering()
476     {
477     if (noError())
478     Elements->optimizeOrdering();
479     if (noError())
480     FaceElements->optimizeOrdering();
481     if (noError())
482     Points->optimizeOrdering();
483     if (noError())
484     ContactElements->optimizeOrdering();
485 jgs 82 }
486 caltinay 4496
487     /// regenerates list of tags in use for node file and element files
488     void Mesh::updateTagList()
489     {
490     if (noError()) Nodes->updateTagList();
491     if (noError()) Elements->updateTagList();
492     if (noError()) FaceElements->updateTagList();
493     if (noError()) Points->updateTagList();
494     if (noError()) ContactElements->updateTagList();
495 jgs 82 }
496 caltinay 4496
497     /// assigns new node reference numbers to all element files
498     void Mesh::relabelElementNodes(const std::vector<int>& newNode, int offset)
499     {
500     Elements->relabelNodes(newNode, offset);
501     FaceElements->relabelNodes(newNode, offset);
502     ContactElements->relabelNodes(newNode, offset);
503     Points->relabelNodes(newNode, offset);
504 ksteube 1312 }
505 gross 2856
506 caltinay 4496 void Mesh::resolveNodeIds()
507 gross 2856 {
508 caltinay 4496 // Initially the element nodes refer to the numbering defined by the global
509     // id assigned to the nodes in the NodeFile. It is also not ensured that
510     // all nodes referred by an element are actually available on the process.
511     // At the output, a local node labeling is used and all nodes are
512     // available. In particular the numbering of the element nodes is between
513     // 0 and NodeFile->numNodes.
514     // The function does not create a distribution of the degrees of freedom.
515 gross 2856
516 caltinay 4496 // find the minimum and maximum id used by elements
517     int min_id=std::numeric_limits<int>::max();
518     int max_id=std::numeric_limits<int>::min();
519     std::pair<int,int> range(Elements->getNodeRange());
520     max_id=std::max(max_id,range.second);
521     min_id=std::min(min_id,range.first);
522     range=FaceElements->getNodeRange();
523     max_id=std::max(max_id,range.second);
524     min_id=std::min(min_id,range.first);
525     range=ContactElements->getNodeRange();
526     max_id=std::max(max_id,range.second);
527     min_id=std::min(min_id,range.first);
528     range=Points->getNodeRange();
529     max_id=std::max(max_id,range.second);
530     min_id=std::min(min_id,range.first);
531     #ifdef Finley_TRACE
532     int global_min_id, global_max_id;
533     #ifdef ESYS_MPI
534     int id_range[2], global_id_range[2];
535     id_range[0]=-min_id;
536     id_range[1]=max_id;
537     MPI_Allreduce(id_range, global_id_range, 2, MPI_INT, MPI_MAX, MPIInfo->comm);
538     global_min_id=-global_id_range[0];
539     global_max_id=global_id_range[1];
540     #else
541     global_min_id=min_id;
542     global_max_id=max_id;
543     #endif
544     printf("Node id range used by elements is %d:%d\n",global_min_id,global_max_id);
545     #endif
546     if (min_id>max_id) {
547     max_id=-1;
548     min_id=0;
549     }
550    
551     // allocate mappings for new local node labeling to global node labeling
552     // (newLocalToGlobalNodeLabels) and global node labeling to the new local
553     // node labeling (globalToNewLocalNodeLabels[i-min_id] is the new local id
554     // of global node i)
555     int len=(max_id>=min_id) ? max_id-min_id+1 : 0;
556 gross 2856
557 caltinay 4496 // mark the nodes referred by elements in usedMask
558     std::vector<short> usedMask(len, -1);
559     markNodes(usedMask, min_id, false);
560 gross 2856
561 caltinay 4496 // create a local labeling newLocalToGlobalNodeLabels of the local nodes
562     // by packing the mask usedMask
563     std::vector<int> newLocalToGlobalNodeLabels=util::packMask(usedMask);
564     const int newNumNodes=newLocalToGlobalNodeLabels.size();
565     usedMask.clear();
566 gross 2856
567 caltinay 4496 // invert the new labeling and shift the index newLocalToGlobalNodeLabels
568     // to global node ids
569     std::vector<int> globalToNewLocalNodeLabels(len, -1);
570 gross 2856
571 caltinay 4496 #pragma omp parallel for
572     for (int n=0; n<newNumNodes; n++) {
573     #ifdef BOUNDS_CHECK
574     if (newLocalToGlobalNodeLabels[n] >= len || newLocalToGlobalNodeLabels[n] < 0) {
575     printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);
576     exit(1);
577     }
578     #endif
579     globalToNewLocalNodeLabels[newLocalToGlobalNodeLabels[n]]=n;
580     newLocalToGlobalNodeLabels[n]+=min_id;
581     }
582 gross 2856
583 caltinay 4496 // create a new table
584     NodeFile *newNodeFile=new NodeFile(getDim(), MPIInfo);
585     if (noError()) {
586     newNodeFile->allocTable(newNumNodes);
587     }
588     if (noError()) {
589     newNodeFile->gather_global(newLocalToGlobalNodeLabels, Nodes);
590     }
591     if (noError()) {
592     delete Nodes;
593     Nodes=newNodeFile;
594     // relabel nodes of the elements
595     relabelElementNodes(globalToNewLocalNodeLabels, min_id);
596     }
597     }
598 gross 2856
599 caltinay 4496 /// sets new coordinates for the nodes
600     void Mesh::setCoordinates(const escript::Data& newX)
601     {
602     Nodes->setCoordinates(newX);
603 gross 2856 }
604 caltinay 4496
605     void Mesh::addTagMap(const char* name, int tag_key)
606     {
607     tagMap[std::string(name)]=tag_key;
608     }
609    
610     int Mesh::getTag(const char* name) const
611     {
612     TagMap::const_iterator it = tagMap.find(name);
613     if (it == tagMap.end()) {
614     std::stringstream ss;
615     ss << "getTag: unknown tag name " << name << ".";
616     const std::string errorMsg(ss.str());
617     setError(VALUE_ERROR, errorMsg.c_str());
618     return -1;
619     }
620     return it->second;
621     }
622    
623     bool Mesh::isValidTagName(const char* name) const
624     {
625     return (tagMap.count(std::string(name)) > 0);
626     }
627    
628     } // namespace finley
629    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26