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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6030 - (show annotations)
Tue Mar 8 02:07:16 2016 UTC (2 years, 11 months ago) by caltinay
File size: 17323 byte(s)
sync with the now stable trunk.

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

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_createNodeFileMappings.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_createNodeFileMappings.cpp:2682-2741 /branches/pasowrap/dudley/src/Mesh_createNodeFileMappings.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_createNodeFileMappings.cpp:3871-3891 /branches/restext/finley/src/Mesh_createNodeFileMappings.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Mesh_createNodeFileMappings.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_createNodeFileMappings.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Mesh_createNodeFileMappings.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_createNodeFileMappings.cpp:3517-3974 /release/3.0/finley/src/Mesh_createNodeFileMappings.cpp:2591-2601 /release/4.0/dudley/src/Mesh_createNodeFileMappings.cpp:5380-5406 /trunk/dudley/src/Mesh_createNodeFileMappings.cpp:4257-4344,5898-6029 /trunk/ripley/test/python/dudley/src/Mesh_createNodeFileMappings.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26