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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (hide annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 4803 byte(s)
Assorted spelling fixes.

1 ksteube 1315
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 1315
16 jfenwick 3981 /**********************************************************************************************/
17 ksteube 1315
18 jfenwick 3086 /* Dudley: Mesh: optimizes the labeling of the DOFs on each processor */
19 ksteube 1315
20 jfenwick 3981 /**********************************************************************************************/
21 ksteube 1315
22     #include "Mesh.h"
23     #include "IndexList.h"
24    
25 jfenwick 3981 /************************************************************************************/
26 ksteube 1315
27 jfenwick 3224 void Dudley_Mesh_optimizeDOFLabeling(Dudley_Mesh * in, dim_t * distribution)
28     {
29 ksteube 1315
30 jfenwick 3224 index_t myFirstVertex, myLastVertex, *newGlobalDOFID = NULL, firstVertex, lastVertex;
31     register index_t k;
32     dim_t mpiSize, myNumVertices, len, p, i;
33     Paso_Pattern *pattern = NULL;
34 caltinay 3603 Esys_MPI_rank myRank, current_rank;
35 jfenwick 3224 Dudley_IndexList *index_list = NULL;
36 jfenwick 3231 #ifdef ESYS_MPI
37 caltinay 3603 Esys_MPI_rank dest, source;
38 jfenwick 3224 MPI_Status status;
39     #endif
40 ksteube 1315
41 jfenwick 3224 if (in == NULL)
42     return;
43     if (in->Nodes == NULL)
44     return;
45 ksteube 1315
46 jfenwick 3224 myRank = in->MPIInfo->rank;
47     mpiSize = in->MPIInfo->size;
48     myFirstVertex = distribution[myRank];
49     myLastVertex = distribution[myRank + 1];
50     myNumVertices = myLastVertex - myFirstVertex;
51     len = 0;
52     for (p = 0; p < mpiSize; ++p)
53     len = MAX(len, distribution[p + 1] - distribution[p]);
54 ksteube 1315
55 jfenwick 3224 index_list = TMPMEMALLOC(myNumVertices, Dudley_IndexList);
56     newGlobalDOFID = TMPMEMALLOC(len, index_t);
57     /* create the adjacency structure xadj and adjncy */
58     if (!(Dudley_checkPtr(index_list) || Dudley_checkPtr(newGlobalDOFID)))
59     {
60     #pragma omp parallel private(i)
61     {
62     #pragma omp for schedule(static)
63     for (i = 0; i < myNumVertices; ++i)
64     {
65     index_list[i].extension = NULL;
66     index_list[i].n = 0;
67     }
68 caltinay 4286 /* insert contributions from element matrices into columns index index_list: */
69 jfenwick 3224 Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
70     in->Elements, in->Nodes->globalDegreesOfFreedom,
71     in->Nodes->globalDegreesOfFreedom);
72     Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
73     in->FaceElements,
74     in->Nodes->globalDegreesOfFreedom,
75     in->Nodes->globalDegreesOfFreedom);
76     Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list, myFirstVertex, myLastVertex,
77     in->Points, in->Nodes->globalDegreesOfFreedom,
78     in->Nodes->globalDegreesOfFreedom);
79     }
80     /* create the local matrix pattern */
81     pattern =
82     Dudley_IndexList_createPattern(0, myNumVertices, index_list, myFirstVertex, myLastVertex, -myFirstVertex);
83 ksteube 1315
84 jfenwick 3224 /* clean up index list */
85     if (index_list != NULL)
86     {
87     #pragma omp parallel for private(i)
88     for (i = 0; i < myNumVertices; ++i)
89     Dudley_IndexList_free(index_list[i].extension);
90     }
91 ksteube 1315
92 jfenwick 3224 if (Dudley_noError())
93     Paso_Pattern_reduceBandwidth(pattern, newGlobalDOFID);
94 ksteube 1315
95 jfenwick 3224 Paso_Pattern_free(pattern);
96     }
97 jfenwick 3227 Esys_MPIInfo_noError(in->MPIInfo);
98 jfenwick 3224 if (Dudley_noError())
99     {
100     /* shift new labeling to create a global id */
101     #pragma omp parallel for private(i)
102     for (i = 0; i < myNumVertices; ++i)
103     newGlobalDOFID[i] += myFirstVertex;
104 ksteube 1315
105 jfenwick 3224 /* distribute new labeling to other processors */
106 caltinay 3603 #ifdef ESYS_MPI
107 jfenwick 3227 dest = Esys_MPIInfo_mod(mpiSize, myRank + 1);
108     source = Esys_MPIInfo_mod(mpiSize, myRank - 1);
109 caltinay 3603 #endif
110 jfenwick 3224 current_rank = myRank;
111     for (p = 0; p < mpiSize; ++p)
112     {
113     firstVertex = distribution[current_rank];
114     lastVertex = distribution[current_rank + 1];
115     #pragma omp parallel for private(i,k)
116     for (i = 0; i < in->Nodes->numNodes; ++i)
117     {
118     k = in->Nodes->globalDegreesOfFreedom[i];
119     if ((firstVertex <= k) && (k < lastVertex))
120     {
121     in->Nodes->globalDegreesOfFreedom[i] = newGlobalDOFID[k - firstVertex];
122     }
123     }
124 ksteube 1315
125 jfenwick 3224 if (p < mpiSize - 1)
126     { /* the final send can be skipped */
127 jfenwick 3231 #ifdef ESYS_MPI
128 jfenwick 3224 MPI_Sendrecv_replace(newGlobalDOFID, len, MPI_INT,
129     dest, in->MPIInfo->msg_tag_counter,
130     source, in->MPIInfo->msg_tag_counter, in->MPIInfo->comm, &status);
131     #endif
132     in->MPIInfo->msg_tag_counter++;
133 jfenwick 3227 current_rank = Esys_MPIInfo_mod(mpiSize, current_rank - 1);
134 jfenwick 3224 }
135     }
136     }
137     TMPMEMFREE(index_list);
138     TMPMEMFREE(newGlobalDOFID);
139 ksteube 1315 #if 0
140 jfenwick 3224 for (i = 0; i < in->Nodes->numNodes; ++i)
141     printf("%d ", in->Nodes->globalDegreesOfFreedom[i]);
142     printf("\n");
143 ksteube 1315 #endif
144 jfenwick 3224 return;
145 ksteube 1315 }

  ViewVC Help
Powered by ViewVC 1.1.26