/[escript]/branches/doubleplusgood/dudley/src/Mesh_createNodeFileMappings.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Mesh_createNodeFileMappings.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 17733 byte(s)
like sand though the hourglass
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 /************************************************************************************/
17
18 /* Dudley: NodeFile : creates the mappings using the indexReducedNodes */
19 /* no distribution is happening */
20
21 /************************************************************************************/
22
23 #include "Mesh.h"
24 #define UNUSED -1
25
26 /************************************************************************************/
27
28 void Dudley_Mesh_createDOFMappingAndCoupling(Dudley_Mesh * in, bool_t use_reduced_elements)
29 {
30 index_t min_DOF, max_DOF, *shared = NULL, *offsetInShared = NULL, *locDOFMask =
31 NULL, i, k, myFirstDOF, myLastDOF, *nodeMask = NULL, firstDOF, lastDOF, *globalDOFIndex, *wanted_DOFs = NULL;
32 dim_t mpiSize, len_loc_dof, numNeighbors, n, lastn, numNodes, *rcv_len = NULL, *snd_len = NULL, count;
33 Esys_MPI_rank myRank, p, p_min, p_max, *neighbor = NULL;
34 Paso_SharedComponents *rcv_shcomp = NULL, *snd_shcomp = NULL;
35 Dudley_NodeMapping *this_mapping = NULL;
36 Paso_Connector *this_connector = NULL;
37 Paso_Distribution *dof_distribution;
38 Esys_MPIInfo *mpi_info = in->MPIInfo;
39 #ifdef ESYS_MPI
40 MPI_Request *mpi_requests = NULL;
41 MPI_Status *mpi_stati = NULL;
42 #else
43 int *mpi_requests = NULL, *mpi_stati = NULL;
44 #endif
45
46 numNodes = in->Nodes->numNodes;
47 if (use_reduced_elements)
48 {
49 dof_distribution = in->Nodes->reducedDegreesOfFreedomDistribution;
50 globalDOFIndex = in->Nodes->globalReducedDOFIndex;
51 }
52 else
53 {
54 dof_distribution = in->Nodes->degreesOfFreedomDistribution;
55 globalDOFIndex = in->Nodes->globalDegreesOfFreedom;
56 }
57 myFirstDOF = Paso_Distribution_getFirstComponent(dof_distribution);
58 myLastDOF = Paso_Distribution_getLastComponent(dof_distribution);
59
60 mpiSize = mpi_info->size;
61 myRank = mpi_info->rank;
62
63 min_DOF = Dudley_Util_getFlaggedMinInt(1, numNodes, globalDOFIndex, -1);
64 max_DOF = Dudley_Util_getFlaggedMaxInt(1, numNodes, globalDOFIndex, -1);
65
66 if (max_DOF < min_DOF)
67 {
68 min_DOF = myFirstDOF;
69 max_DOF = myLastDOF - 1;
70 }
71
72 p_min = mpiSize;
73 p_max = -1;
74 if (max_DOF >= min_DOF)
75 {
76 for (p = 0; p < mpiSize; ++p)
77 {
78 if (dof_distribution->first_component[p] <= min_DOF)
79 p_min = p;
80 if (dof_distribution->first_component[p] <= max_DOF)
81 p_max = p;
82 }
83 }
84
85 len_loc_dof = max_DOF - min_DOF + 1;
86 if (!((min_DOF <= myFirstDOF) && (myLastDOF - 1 <= max_DOF)))
87 {
88 Dudley_setError(SYSTEM_ERROR, "Local elements do not span local degrees of freedom.");
89 return;
90 }
91 rcv_len = new dim_t[mpiSize];
92 snd_len = new dim_t[mpiSize];
93 #ifdef ESYS_MPI
94 mpi_requests = MEMALLOC(mpiSize * 2, MPI_Request);
95 mpi_stati = MEMALLOC(mpiSize * 2, MPI_Status);
96 #else
97 mpi_requests = MEMALLOC(mpiSize * 2, int);
98 mpi_stati = MEMALLOC(mpiSize * 2, int);
99 #endif
100 wanted_DOFs = new index_t[numNodes];
101 nodeMask = new index_t[numNodes];
102 neighbor = new Esys_MPI_rank[mpiSize];
103 shared = new index_t[numNodes * (p_max - p_min + 1)];
104 offsetInShared = new index_t[mpiSize + 1];
105 locDOFMask = new index_t[len_loc_dof];
106 if (!
107 (Dudley_checkPtr(neighbor) || Dudley_checkPtr(shared) || Dudley_checkPtr(offsetInShared)
108 || Dudley_checkPtr(locDOFMask) || Dudley_checkPtr(nodeMask) || Dudley_checkPtr(rcv_len)
109 || Dudley_checkPtr(snd_len) || Dudley_checkPtr(mpi_requests) || Dudley_checkPtr(mpi_stati)
110 || Dudley_checkPtr(mpi_stati)))
111 {
112
113 memset(rcv_len, 0, sizeof(dim_t) * mpiSize);
114 #pragma omp parallel
115 {
116 #pragma omp for private(i) schedule(static)
117 for (i = 0; i < len_loc_dof; ++i)
118 locDOFMask[i] = UNUSED;
119 #pragma omp for private(i) schedule(static)
120 for (i = 0; i < numNodes; ++i)
121 nodeMask[i] = UNUSED;
122 #pragma omp for private(i,k) schedule(static)
123 for (i = 0; i < numNodes; ++i)
124 {
125 k = globalDOFIndex[i];
126 if (k > -1)
127 {
128 locDOFMask[k - min_DOF] = UNUSED - 1;
129 #ifdef BOUNDS_CHECK
130 if ((k - min_DOF) >= len_loc_dof)
131 {
132 printf("BOUNDS_CHECK %s %d i=%d k=%d min_DOF=%d\n", __FILE__, __LINE__, i, k, min_DOF);
133 exit(1);
134 }
135 #endif
136 }
137 }
138
139 #pragma omp for private(i) schedule(static)
140 for (i = myFirstDOF - min_DOF; i < myLastDOF - min_DOF; ++i)
141 {
142 locDOFMask[i] = i - myFirstDOF + min_DOF;
143 #ifdef BOUNDS_CHECK
144 if (i < 0 || i >= len_loc_dof)
145 {
146 printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
147 exit(1);
148 }
149 #endif
150 }
151 }
152
153 numNeighbors = 0;
154 n = 0;
155 lastn = n;
156 for (p = p_min; p <= p_max; ++p)
157 {
158 firstDOF = MAX(min_DOF, dof_distribution->first_component[p]);
159 lastDOF = MIN(max_DOF + 1, dof_distribution->first_component[p + 1]);
160 if (p != myRank)
161 {
162 for (i = firstDOF - min_DOF; i < lastDOF - min_DOF; ++i)
163 {
164 #ifdef BOUNDS_CHECK
165 if (i < 0 || i >= len_loc_dof)
166 {
167 printf("BOUNDS_CHECK %s %d p=%d i=%d\n", __FILE__, __LINE__, p, i);
168 exit(1);
169 }
170 #endif
171 if (locDOFMask[i] == UNUSED - 1)
172 {
173 locDOFMask[i] = myLastDOF - myFirstDOF + n;
174 wanted_DOFs[n] = i + min_DOF;
175 ++n;
176 }
177 }
178 if (n > lastn)
179 {
180 rcv_len[p] = n - lastn;
181 neighbor[numNeighbors] = p;
182 #ifdef BOUNDS_CHECK
183 if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
184 {
185 printf("BOUNDS_CHECK %s %d p=%d numNeighbors=%d n=%d\n", __FILE__, __LINE__, p, numNeighbors,
186 n);
187 exit(1);
188 }
189 #endif
190 offsetInShared[numNeighbors] = lastn;
191 numNeighbors++;
192 lastn = n;
193 }
194 }
195 }
196 #ifdef BOUNDS_CHECK
197 if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
198 {
199 printf("BOUNDS_CHECK %s %d numNeighbors=%d\n", __FILE__, __LINE__, numNeighbors);
200 exit(1);
201 }
202 #endif
203 offsetInShared[numNeighbors] = lastn;
204
205 /* assign new DOF labels to nodes */
206 #pragma omp parallel for private(i,k) schedule(static)
207 for (i = 0; i < numNodes; ++i)
208 {
209 k = globalDOFIndex[i];
210 if (k > -1)
211 nodeMask[i] = locDOFMask[k - min_DOF];
212 }
213
214 /* now we can set the mapping from nodes to local DOFs */
215 this_mapping = Dudley_NodeMapping_alloc(numNodes, nodeMask, UNUSED);
216 /* define how to get DOF values for controlled bu other processors */
217 #ifdef BOUNDS_CHECK
218 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
219 {
220 if (i < 0 || i >= numNodes * (p_max - p_min + 1))
221 {
222 printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
223 exit(1);
224 }
225 }
226 #endif
227 #pragma omp parallel for private(i) schedule(static)
228 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
229 shared[i] = myLastDOF - myFirstDOF + i;
230
231 rcv_shcomp =
232 Paso_SharedComponents_alloc(myLastDOF - myFirstDOF, numNeighbors, neighbor, shared, offsetInShared, 1, 0,
233 mpi_info);
234
235 /*
236 * now we build the sender
237 */
238 #ifdef ESYS_MPI
239 MPI_Alltoall(rcv_len, 1, MPI_INT, snd_len, 1, MPI_INT, mpi_info->comm);
240 #else
241 for (p = 0; p < mpiSize; ++p)
242 snd_len[p] = rcv_len[p];
243 #endif
244 count = 0;
245 for (p = 0; p < rcv_shcomp->numNeighbors; p++)
246 {
247 #ifdef ESYS_MPI
248 MPI_Isend(&(wanted_DOFs[rcv_shcomp->offsetInShared[p]]),
249 rcv_shcomp->offsetInShared[p + 1] - rcv_shcomp->offsetInShared[p], MPI_INT,
250 rcv_shcomp->neighbor[p], mpi_info->msg_tag_counter + myRank, mpi_info->comm,
251 &mpi_requests[count]);
252 #endif
253 count++;
254 }
255 n = 0;
256 numNeighbors = 0;
257 for (p = 0; p < mpiSize; p++)
258 {
259 if (snd_len[p] > 0)
260 {
261 #ifdef ESYS_MPI
262 MPI_Irecv(&(shared[n]), snd_len[p],
263 MPI_INT, p, mpi_info->msg_tag_counter + p, mpi_info->comm, &mpi_requests[count]);
264 #endif
265 count++;
266 neighbor[numNeighbors] = p;
267 offsetInShared[numNeighbors] = n;
268 numNeighbors++;
269 n += snd_len[p];
270 }
271 }
272 mpi_info->msg_tag_counter += mpi_info->size;
273 offsetInShared[numNeighbors] = n;
274 #ifdef ESYS_MPI
275 MPI_Waitall(count, mpi_requests, mpi_stati);
276 #endif
277 /* map global ids to local id's */
278 #pragma omp parallel for private(i) schedule(static)
279 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
280 {
281 shared[i] = locDOFMask[shared[i] - min_DOF];
282 }
283
284 snd_shcomp =
285 Paso_SharedComponents_alloc(myLastDOF - myFirstDOF, numNeighbors, neighbor, shared, offsetInShared, 1, 0,
286 dof_distribution->mpi_info);
287
288 if (Dudley_noError())
289 this_connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
290 /* assign new DOF labels to nodes */
291 Paso_SharedComponents_free(rcv_shcomp);
292 Paso_SharedComponents_free(snd_shcomp);
293 }
294 delete[] rcv_len;
295 delete[] snd_len;
296 delete[] mpi_requests;
297 delete[] mpi_stati;
298 delete[] wanted_DOFs;
299 delete[] nodeMask;
300 delete[] neighbor;
301 delete[] shared;
302 delete[] offsetInShared;
303 delete[] locDOFMask;
304 if (Dudley_noError())
305 {
306 if (use_reduced_elements)
307 {
308 in->Nodes->reducedDegreesOfFreedomMapping = this_mapping;
309 in->Nodes->reducedDegreesOfFreedomConnector = this_connector;
310 }
311 else
312 {
313 in->Nodes->degreesOfFreedomMapping = this_mapping;
314 in->Nodes->degreesOfFreedomConnector = this_connector;
315 }
316 }
317 else
318 {
319 Dudley_NodeMapping_free(this_mapping);
320 Paso_Connector_free(this_connector);
321
322 }
323 }
324
325 void Dudley_Mesh_createMappings(Dudley_Mesh * mesh, index_t * dof_distribution, index_t * node_distribution)
326 {
327 int i;
328 index_t *maskReducedNodes = NULL, *indexReducedNodes = NULL;
329 dim_t numReducedNodes;
330
331 maskReducedNodes = new index_t[mesh->Nodes->numNodes];
332 indexReducedNodes = new index_t[mesh->Nodes->numNodes];
333
334 if (!(Dudley_checkPtr(maskReducedNodes) || Dudley_checkPtr(indexReducedNodes)))
335 {
336 #pragma omp parallel for private(i) schedule(static)
337 for (i = 0; i < mesh->Nodes->numNodes; ++i)
338 maskReducedNodes[i] = -1;
339 Dudley_Mesh_markNodes(maskReducedNodes, 0, mesh, TRUE);
340
341 numReducedNodes = Dudley_Util_packMask(mesh->Nodes->numNodes, maskReducedNodes, indexReducedNodes);
342 if (Dudley_noError())
343 Dudley_Mesh_createNodeFileMappings(mesh, numReducedNodes, indexReducedNodes, dof_distribution,
344 node_distribution);
345 }
346
347 delete[] maskReducedNodes;
348 delete[] indexReducedNodes;
349 }
350
351 void Dudley_Mesh_createNodeFileMappings(Dudley_Mesh * in, dim_t numReducedNodes, index_t * indexReducedNodes,
352 index_t * dof_first_component, index_t * nodes_first_component)
353 {
354
355 index_t myFirstDOF, myLastDOF, myFirstNode, myLastNode, *reduced_dof_first_component = NULL, *nodeMask = NULL,
356 *reduced_nodes_first_component = NULL, k, *maskMyReducedDOF = NULL, *indexMyReducedDOF =
357 NULL, *maskMyReducedNodes = NULL, *indexMyReducedNodes = NULL;
358 dim_t myNumDOF, myNumNodes, myNumReducedNodes, myNumReducedDOF, globalNumReducedNodes, globalNumReducedDOF, i,
359 mpiSize;
360 Esys_MPI_rank myRank;
361
362 mpiSize = in->Nodes->MPIInfo->size;
363 myRank = in->Nodes->MPIInfo->rank;
364
365 /* mark the nodes used by the reduced mesh */
366
367 reduced_dof_first_component = new index_t[mpiSize + 1];
368 reduced_nodes_first_component = new index_t[mpiSize + 1];
369
370 if (!(Dudley_checkPtr(reduced_dof_first_component) || Dudley_checkPtr(reduced_nodes_first_component)))
371 {
372
373 myFirstDOF = dof_first_component[myRank];
374 myLastDOF = dof_first_component[myRank + 1];
375 myNumDOF = myLastDOF - myFirstDOF;
376
377 myFirstNode = nodes_first_component[myRank];
378 myLastNode = nodes_first_component[myRank + 1];
379 myNumNodes = myLastNode - myFirstNode;
380
381 maskMyReducedDOF = new index_t[myNumDOF];
382 indexMyReducedDOF = new index_t[myNumDOF];
383 maskMyReducedNodes = new index_t[myNumNodes];
384 indexMyReducedNodes = new index_t[myNumNodes];
385
386 if (!
387 (Dudley_checkPtr(maskMyReducedDOF) || Dudley_checkPtr(indexMyReducedDOF)
388 || Dudley_checkPtr(maskMyReducedNodes) || Dudley_checkPtr(indexMyReducedNodes)))
389 {
390
391 #pragma omp parallel private(i)
392 {
393 #pragma omp for schedule(static)
394 for (i = 0; i < myNumNodes; ++i)
395 maskMyReducedNodes[i] = -1;
396 #pragma omp for schedule(static)
397 for (i = 0; i < myNumDOF; ++i)
398 maskMyReducedDOF[i] = -1;
399 #pragma omp for private(k) schedule(static)
400 for (i = 0; i < numReducedNodes; ++i)
401 {
402 k = in->Nodes->globalNodesIndex[indexReducedNodes[i]];
403 if ((k >= myFirstNode) && (myLastNode > k))
404 maskMyReducedNodes[k - myFirstNode] = i;
405 k = in->Nodes->globalDegreesOfFreedom[indexReducedNodes[i]];
406 if ((k >= myFirstDOF) && (myLastDOF > k))
407 {
408 maskMyReducedDOF[k - myFirstDOF] = i;
409 }
410 }
411 }
412 myNumReducedNodes = Dudley_Util_packMask(myNumNodes, maskMyReducedNodes, indexMyReducedNodes);
413 myNumReducedDOF = Dudley_Util_packMask(myNumDOF, maskMyReducedDOF, indexMyReducedDOF);
414
415 #ifdef ESYS_MPI
416 MPI_Allgather(&myNumReducedNodes, 1, MPI_INT, reduced_nodes_first_component, 1, MPI_INT,
417 in->Nodes->MPIInfo->comm);
418 MPI_Allgather(&myNumReducedDOF, 1, MPI_INT, reduced_dof_first_component, 1, MPI_INT,
419 in->Nodes->MPIInfo->comm);
420 #else
421 reduced_nodes_first_component[0] = myNumReducedNodes;
422 reduced_dof_first_component[0] = myNumReducedDOF;
423 #endif
424 globalNumReducedNodes = 0;
425 globalNumReducedDOF = 0;
426 for (i = 0; i < mpiSize; ++i)
427 {
428 k = reduced_nodes_first_component[i];
429 reduced_nodes_first_component[i] = globalNumReducedNodes;
430 globalNumReducedNodes += k;
431
432 k = reduced_dof_first_component[i];
433 reduced_dof_first_component[i] = globalNumReducedDOF;
434 globalNumReducedDOF += k;
435 }
436 reduced_nodes_first_component[mpiSize] = globalNumReducedNodes;
437 reduced_dof_first_component[mpiSize] = globalNumReducedDOF;
438 /* ==== distribution of Nodes =============================== */
439 in->Nodes->nodesDistribution = Paso_Distribution_alloc(in->Nodes->MPIInfo, nodes_first_component, 1, 0);
440
441 /* ==== distribution of DOFs =============================== */
442 in->Nodes->degreesOfFreedomDistribution =
443 Paso_Distribution_alloc(in->Nodes->MPIInfo, dof_first_component, 1, 0);
444
445 /* ==== distribution of reduced Nodes =============================== */
446 in->Nodes->reducedNodesDistribution =
447 Paso_Distribution_alloc(in->Nodes->MPIInfo, reduced_nodes_first_component, 1, 0);
448
449 /* ==== distribution of reduced DOF =============================== */
450 in->Nodes->reducedDegreesOfFreedomDistribution =
451 Paso_Distribution_alloc(in->Nodes->MPIInfo, reduced_dof_first_component, 1, 0);
452 }
453 delete[] maskMyReducedDOF;
454 delete[] indexMyReducedDOF;
455 delete[] maskMyReducedNodes;
456 delete[] indexMyReducedNodes;
457 }
458 delete[] reduced_dof_first_component;
459 delete[] reduced_nodes_first_component;
460
461 nodeMask = new index_t[in->Nodes->numNodes];
462 if (!Dudley_checkPtr(nodeMask) && Dudley_noError())
463 {
464
465 /* ==== nodes mapping which is a dummy structure ======== */
466 #pragma omp parallel for private(i) schedule(static)
467 for (i = 0; i < in->Nodes->numNodes; ++i)
468 nodeMask[i] = i;
469 in->Nodes->nodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
470
471 /* ==== mapping between nodes and reduced nodes ========== */
472 #pragma omp parallel for private(i) schedule(static)
473 for (i = 0; i < in->Nodes->numNodes; ++i)
474 nodeMask[i] = UNUSED;
475 #pragma omp parallel for private(i) schedule(static)
476 for (i = 0; i < numReducedNodes; ++i)
477 nodeMask[indexReducedNodes[i]] = i;
478 in->Nodes->reducedNodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
479 }
480 delete[] nodeMask;
481 /* ==== mapping between nodes and DOFs + DOF connector ========== */
482 if (Dudley_noError())
483 Dudley_Mesh_createDOFMappingAndCoupling(in, FALSE);
484 /* ==== mapping between nodes and reduced DOFs + reduced DOF connector ========== */
485 if (Dudley_noError())
486 Dudley_Mesh_createDOFMappingAndCoupling(in, TRUE);
487
488 /* get the Ids for DOFs and reduced nodes */
489 if (Dudley_noError())
490 {
491 #pragma omp parallel private(i)
492 {
493 #pragma omp for
494 for (i = 0; i < in->Nodes->reducedNodesMapping->numTargets; ++i)
495 in->Nodes->reducedNodesId[i] = in->Nodes->Id[in->Nodes->reducedNodesMapping->map[i]];
496 #pragma omp for
497 for (i = 0; i < in->Nodes->degreesOfFreedomMapping->numTargets; ++i)
498 in->Nodes->degreesOfFreedomId[i] = in->Nodes->Id[in->Nodes->degreesOfFreedomMapping->map[i]];
499 #pragma omp for
500 for (i = 0; i < in->Nodes->reducedDegreesOfFreedomMapping->numTargets; ++i)
501 in->Nodes->reducedDegreesOfFreedomId[i] =
502 in->Nodes->Id[in->Nodes->reducedDegreesOfFreedomMapping->map[i]];
503 }
504 }
505 else
506 {
507 Dudley_NodeMapping_free(in->Nodes->nodesMapping);
508 Dudley_NodeMapping_free(in->Nodes->reducedNodesMapping);
509 Dudley_NodeMapping_free(in->Nodes->degreesOfFreedomMapping);
510 Dudley_NodeMapping_free(in->Nodes->reducedDegreesOfFreedomMapping);
511 Paso_Distribution_free(in->Nodes->nodesDistribution);
512 Paso_Distribution_free(in->Nodes->reducedNodesDistribution);
513 Paso_Distribution_free(in->Nodes->degreesOfFreedomDistribution);
514 Paso_Distribution_free(in->Nodes->reducedDegreesOfFreedomDistribution);
515 Paso_Connector_free(in->Nodes->degreesOfFreedomConnector);
516 Paso_Connector_free(in->Nodes->reducedDegreesOfFreedomConnector);
517 in->Nodes->nodesMapping = NULL;
518 in->Nodes->reducedNodesMapping = NULL;
519 in->Nodes->degreesOfFreedomMapping = NULL;
520 in->Nodes->reducedDegreesOfFreedomMapping = NULL;
521 in->Nodes->nodesDistribution = NULL;
522 in->Nodes->reducedNodesDistribution = NULL;
523 in->Nodes->degreesOfFreedomDistribution = NULL;
524 in->Nodes->reducedDegreesOfFreedomDistribution = NULL;
525 in->Nodes->degreesOfFreedomConnector = NULL;
526 in->Nodes->reducedDegreesOfFreedomConnector = NULL;
527 }
528 }

  ViewVC Help
Powered by ViewVC 1.1.26