/[escript]/trunk/dudley/src/Mesh_optimizeDOFLabeling.c
ViewVC logotype

Contents of /trunk/dudley/src/Mesh_optimizeDOFLabeling.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 4576 byte(s)
Merging dudley and scons updates from branches

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 /************************************************************************/
15
16 /* Dudley: Mesh: optimizes the labeling of the DOFs on each processor */
17
18 /************************************************************************/
19
20 #include "Mesh.h"
21 #include "IndexList.h"
22
23 /**************************************************************/
24
25 void Dudley_Mesh_optimizeDOFLabeling(Dudley_Mesh * in, dim_t * distribution)
26 {
27
28 index_t myFirstVertex, myLastVertex, *newGlobalDOFID = NULL, firstVertex, lastVertex;
29 register index_t k;
30 dim_t mpiSize, myNumVertices, len, p, i;
31 Paso_Pattern *pattern = NULL;
32 Esys_MPI_rank myRank, dest, source, current_rank;
33 Dudley_IndexList *index_list = NULL;
34 #ifdef ESYS_MPI
35 MPI_Status status;
36 #endif
37
38 if (in == NULL)
39 return;
40 if (in->Nodes == NULL)
41 return;
42
43 myRank = in->MPIInfo->rank;
44 mpiSize = in->MPIInfo->size;
45 myFirstVertex = distribution[myRank];
46 myLastVertex = distribution[myRank + 1];
47 myNumVertices = myLastVertex - myFirstVertex;
48 len = 0;
49 for (p = 0; p < mpiSize; ++p)
50 len = MAX(len, distribution[p + 1] - distribution[p]);
51
52 index_list = TMPMEMALLOC(myNumVertices, Dudley_IndexList);
53 newGlobalDOFID = TMPMEMALLOC(len, index_t);
54 /* create the adjacency structure xadj and adjncy */
55 if (!(Dudley_checkPtr(index_list) || Dudley_checkPtr(newGlobalDOFID)))
56 {
57 #pragma omp parallel private(i)
58 {
59 #pragma omp for schedule(static)
60 for (i = 0; i < myNumVertices; ++i)
61 {
62 index_list[i].extension = NULL;
63 index_list[i].n = 0;
64 }
65 /* insert contributions from element matrices into colums index index_list: */
66 Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
67 in->Elements, in->Nodes->globalDegreesOfFreedom,
68 in->Nodes->globalDegreesOfFreedom);
69 Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
70 in->FaceElements,
71 in->Nodes->globalDegreesOfFreedom,
72 in->Nodes->globalDegreesOfFreedom);
73 Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
74 in->Points, in->Nodes->globalDegreesOfFreedom,
75 in->Nodes->globalDegreesOfFreedom);
76 }
77 /* create the local matrix pattern */
78 pattern =
79 Dudley_IndexList_createPattern(0, myNumVertices, index_list, myFirstVertex, myLastVertex, -myFirstVertex);
80
81 /* clean up index list */
82 if (index_list != NULL)
83 {
84 #pragma omp parallel for private(i)
85 for (i = 0; i < myNumVertices; ++i)
86 Dudley_IndexList_free(index_list[i].extension);
87 }
88
89 if (Dudley_noError())
90 Paso_Pattern_reduceBandwidth(pattern, newGlobalDOFID);
91
92 Paso_Pattern_free(pattern);
93 }
94 Esys_MPIInfo_noError(in->MPIInfo);
95 if (Dudley_noError())
96 {
97 /* shift new labeling to create a global id */
98 #pragma omp parallel for private(i)
99 for (i = 0; i < myNumVertices; ++i)
100 newGlobalDOFID[i] += myFirstVertex;
101
102 /* distribute new labeling to other processors */
103 dest = Esys_MPIInfo_mod(mpiSize, myRank + 1);
104 source = Esys_MPIInfo_mod(mpiSize, myRank - 1);
105 current_rank = myRank;
106 for (p = 0; p < mpiSize; ++p)
107 {
108 firstVertex = distribution[current_rank];
109 lastVertex = distribution[current_rank + 1];
110 #pragma omp parallel for private(i,k)
111 for (i = 0; i < in->Nodes->numNodes; ++i)
112 {
113 k = in->Nodes->globalDegreesOfFreedom[i];
114 if ((firstVertex <= k) && (k < lastVertex))
115 {
116 in->Nodes->globalDegreesOfFreedom[i] = newGlobalDOFID[k - firstVertex];
117 }
118 }
119
120 if (p < mpiSize - 1)
121 { /* the final send can be skipped */
122 #ifdef ESYS_MPI
123 MPI_Sendrecv_replace(newGlobalDOFID, len, MPI_INT,
124 dest, in->MPIInfo->msg_tag_counter,
125 source, in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
126 #endif
127 in->MPIInfo->msg_tag_counter++;
128 current_rank = Esys_MPIInfo_mod(mpiSize, current_rank - 1);
129 }
130 }
131 }
132 TMPMEMFREE(index_list);
133 TMPMEMFREE(newGlobalDOFID);
134 #if 0
135 for (i = 0; i < in->Nodes->numNodes; ++i)
136 printf("%d ", in->Nodes->globalDegreesOfFreedom[i]);
137 printf("\n");
138 #endif
139 return;
140 }

  ViewVC Help
Powered by ViewVC 1.1.26