/[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 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 16782 byte(s)
Much needed sync with 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 if (!((min_DOF <= myFirstDOF) && (myLastDOF - 1 <= max_DOF)))
89 {
90 throw DudleyException("Local elements do not span local degrees of freedom.");
91 }
92 rcv_len = new dim_t[mpiSize];
93 snd_len = new dim_t[mpiSize];
94 #ifdef ESYS_MPI
95 mpi_requests = new MPI_Request[mpiSize * 2];
96 mpi_stati = new MPI_Status[mpiSize * 2];
97 #else
98 mpi_requests = new int[mpiSize * 2];
99 mpi_stati = new int[mpiSize * 2];
100 #endif
101 wanted_DOFs = new index_t[numNodes];
102 nodeMask = new index_t[numNodes];
103 neighbor = new int[mpiSize];
104 shared = new index_t[numNodes * (p_max - p_min + 1)];
105 offsetInShared = new index_t[mpiSize + 1];
106 locDOFMask = new index_t[len_loc_dof];
107
108 memset(rcv_len, 0, sizeof(dim_t) * mpiSize);
109 #pragma omp parallel
110 {
111 #pragma omp for private(i) schedule(static)
112 for (i = 0; i < len_loc_dof; ++i)
113 locDOFMask[i] = UNUSED;
114 #pragma omp for private(i) schedule(static)
115 for (i = 0; i < numNodes; ++i)
116 nodeMask[i] = UNUSED;
117 #pragma omp for private(i,k) schedule(static)
118 for (i = 0; i < numNodes; ++i)
119 {
120 k = globalDOFIndex[i];
121 if (k > -1)
122 {
123 locDOFMask[k - min_DOF] = UNUSED - 1;
124 #ifdef BOUNDS_CHECK
125 if ((k - min_DOF) >= len_loc_dof)
126 {
127 printf("BOUNDS_CHECK %s %d i=%d k=%d min_DOF=%d\n", __FILE__, __LINE__, i, k, min_DOF);
128 exit(1);
129 }
130 #endif
131 }
132 }
133
134 #pragma omp for private(i) schedule(static)
135 for (i = myFirstDOF - min_DOF; i < myLastDOF - min_DOF; ++i)
136 {
137 locDOFMask[i] = i - myFirstDOF + min_DOF;
138 #ifdef BOUNDS_CHECK
139 if (i < 0 || i >= len_loc_dof)
140 {
141 printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
142 exit(1);
143 }
144 #endif
145 }
146 }
147
148 numNeighbors = 0;
149 n = 0;
150 lastn = n;
151 for (p = p_min; p <= p_max; ++p)
152 {
153 firstDOF = std::max(min_DOF, dof_distribution->first_component[p]);
154 lastDOF = std::min(max_DOF + 1, dof_distribution->first_component[p + 1]);
155 if (p != myRank)
156 {
157 for (i = firstDOF - min_DOF; i < lastDOF - min_DOF; ++i)
158 {
159 #ifdef BOUNDS_CHECK
160 if (i < 0 || i >= len_loc_dof)
161 {
162 printf("BOUNDS_CHECK %s %d p=%d i=%d\n", __FILE__, __LINE__, p, i);
163 exit(1);
164 }
165 #endif
166 if (locDOFMask[i] == UNUSED - 1)
167 {
168 locDOFMask[i] = myLastDOF - myFirstDOF + n;
169 wanted_DOFs[n] = i + min_DOF;
170 ++n;
171 }
172 }
173 if (n > lastn)
174 {
175 rcv_len[p] = n - lastn;
176 neighbor[numNeighbors] = p;
177 #ifdef BOUNDS_CHECK
178 if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
179 {
180 printf("BOUNDS_CHECK %s %d p=%d numNeighbors=%d n=%d\n", __FILE__, __LINE__, p, numNeighbors,
181 n);
182 exit(1);
183 }
184 #endif
185 offsetInShared[numNeighbors] = lastn;
186 numNeighbors++;
187 lastn = n;
188 }
189 }
190 }
191 #ifdef BOUNDS_CHECK
192 if (numNeighbors < 0 || numNeighbors >= mpiSize + 1)
193 {
194 printf("BOUNDS_CHECK %s %d numNeighbors=%d\n", __FILE__, __LINE__, numNeighbors);
195 exit(1);
196 }
197 #endif
198 offsetInShared[numNeighbors] = lastn;
199
200 /* assign new DOF labels to nodes */
201 #pragma omp parallel for private(i,k) schedule(static)
202 for (i = 0; i < numNodes; ++i)
203 {
204 k = globalDOFIndex[i];
205 if (k > -1)
206 nodeMask[i] = locDOFMask[k - min_DOF];
207 }
208
209 /* now we can set the mapping from nodes to local DOFs */
210 this_mapping = Dudley_NodeMapping_alloc(numNodes, nodeMask, UNUSED);
211 /* define how to get DOF values for controlled bu other processors */
212 #ifdef BOUNDS_CHECK
213 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
214 {
215 if (i < 0 || i >= numNodes * (p_max - p_min + 1))
216 {
217 printf("BOUNDS_CHECK %s %d i=%d\n", __FILE__, __LINE__, i);
218 exit(1);
219 }
220 }
221 #endif
222 #pragma omp parallel for private(i) schedule(static)
223 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
224 shared[i] = myLastDOF - myFirstDOF + i;
225
226 rcv_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,
227 numNeighbors, neighbor, shared, offsetInShared, 1, 0,
228 mpi_info));
229
230 /*
231 * now we build the sender
232 */
233 #ifdef ESYS_MPI
234 MPI_Alltoall(rcv_len, 1, MPI_INT, snd_len, 1, MPI_INT, mpi_info->comm);
235 #else
236 for (p = 0; p < mpiSize; ++p)
237 snd_len[p] = rcv_len[p];
238 #endif
239 count = 0;
240 for (p = 0; p < rcv_shcomp->numNeighbors; p++)
241 {
242 #ifdef ESYS_MPI
243 MPI_Isend(&(wanted_DOFs[rcv_shcomp->offsetInShared[p]]),
244 rcv_shcomp->offsetInShared[p + 1] - rcv_shcomp->offsetInShared[p], MPI_INT,
245 rcv_shcomp->neighbor[p], mpi_info->counter() + myRank, mpi_info->comm,
246 &mpi_requests[count]);
247 #endif
248 count++;
249 }
250 n = 0;
251 numNeighbors = 0;
252 for (p = 0; p < mpiSize; p++)
253 {
254 if (snd_len[p] > 0)
255 {
256 #ifdef ESYS_MPI
257 MPI_Irecv(&(shared[n]), snd_len[p],
258 MPI_INT, p, mpi_info->counter() + p, mpi_info->comm, &mpi_requests[count]);
259 #endif
260 count++;
261 neighbor[numNeighbors] = p;
262 offsetInShared[numNeighbors] = n;
263 numNeighbors++;
264 n += snd_len[p];
265 }
266 }
267 offsetInShared[numNeighbors] = n;
268 #ifdef ESYS_MPI
269 mpi_info->incCounter(mpi_info->size);
270 MPI_Waitall(count, mpi_requests, mpi_stati);
271 #endif
272 /* map global ids to local id's */
273 #pragma omp parallel for private(i) schedule(static)
274 for (i = 0; i < offsetInShared[numNeighbors]; ++i)
275 {
276 shared[i] = locDOFMask[shared[i] - min_DOF];
277 }
278
279 snd_shcomp.reset(new paso::SharedComponents(myLastDOF - myFirstDOF,
280 numNeighbors, neighbor, shared, offsetInShared, 1, 0,
281 dof_distribution->mpi_info));
282
283 this_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
284
285 /* assign new DOF labels to nodes */
286 delete[] rcv_len;
287 delete[] snd_len;
288 delete[] mpi_requests;
289 delete[] mpi_stati;
290 delete[] wanted_DOFs;
291 delete[] nodeMask;
292 delete[] neighbor;
293 delete[] shared;
294 delete[] offsetInShared;
295 delete[] locDOFMask;
296 if (use_reduced_elements)
297 {
298 in->Nodes->reducedDegreesOfFreedomMapping = this_mapping;
299 in->Nodes->reducedDegreesOfFreedomConnector = this_connector;
300 }
301 else
302 {
303 in->Nodes->degreesOfFreedomMapping = this_mapping;
304 in->Nodes->degreesOfFreedomConnector = this_connector;
305 }
306 }
307
308 void Dudley_Mesh_createMappings(Dudley_Mesh* mesh, index_t* dof_distribution, index_t* node_distribution)
309 {
310 int i;
311 index_t *maskReducedNodes = NULL, *indexReducedNodes = NULL;
312 dim_t numReducedNodes;
313
314 maskReducedNodes = new index_t[mesh->Nodes->numNodes];
315 indexReducedNodes = new index_t[mesh->Nodes->numNodes];
316
317 #pragma omp parallel for private(i) schedule(static)
318 for (i = 0; i < mesh->Nodes->numNodes; ++i)
319 maskReducedNodes[i] = -1;
320 Dudley_Mesh_markNodes(maskReducedNodes, 0, mesh, true);
321
322 numReducedNodes = Dudley_Util_packMask(mesh->Nodes->numNodes, maskReducedNodes, indexReducedNodes);
323 Dudley_Mesh_createNodeFileMappings(mesh, numReducedNodes,
324 indexReducedNodes, dof_distribution, node_distribution);
325 delete[] maskReducedNodes;
326 delete[] indexReducedNodes;
327 }
328
329 void Dudley_Mesh_createNodeFileMappings(Dudley_Mesh* in, dim_t numReducedNodes,
330 index_t* indexReducedNodes, index_t* dof_first_component, index_t* nodes_first_component)
331 {
332
333 index_t myFirstDOF, myLastDOF, myFirstNode, myLastNode, *reduced_dof_first_component = NULL, *nodeMask = NULL,
334 *reduced_nodes_first_component = NULL, k, *maskMyReducedDOF = NULL, *indexMyReducedDOF =
335 NULL, *maskMyReducedNodes = NULL, *indexMyReducedNodes = NULL;
336 dim_t myNumDOF, myNumNodes, myNumReducedNodes, myNumReducedDOF, globalNumReducedNodes, globalNumReducedDOF, i,
337 mpiSize;
338 int myRank;
339
340 mpiSize = in->Nodes->MPIInfo->size;
341 myRank = in->Nodes->MPIInfo->rank;
342
343 /* mark the nodes used by the reduced mesh */
344
345 reduced_dof_first_component = new index_t[mpiSize + 1];
346 reduced_nodes_first_component = new index_t[mpiSize + 1];
347
348 myFirstDOF = dof_first_component[myRank];
349 myLastDOF = dof_first_component[myRank + 1];
350 myNumDOF = myLastDOF - myFirstDOF;
351
352 myFirstNode = nodes_first_component[myRank];
353 myLastNode = nodes_first_component[myRank + 1];
354 myNumNodes = myLastNode - myFirstNode;
355
356 maskMyReducedDOF = new index_t[myNumDOF];
357 indexMyReducedDOF = new index_t[myNumDOF];
358 maskMyReducedNodes = new index_t[myNumNodes];
359 indexMyReducedNodes = new index_t[myNumNodes];
360
361 #pragma omp parallel private(i)
362 {
363 #pragma omp for schedule(static)
364 for (i = 0; i < myNumNodes; ++i)
365 maskMyReducedNodes[i] = -1;
366 #pragma omp for schedule(static)
367 for (i = 0; i < myNumDOF; ++i)
368 maskMyReducedDOF[i] = -1;
369 #pragma omp for private(k) schedule(static)
370 for (i = 0; i < numReducedNodes; ++i)
371 {
372 k = in->Nodes->globalNodesIndex[indexReducedNodes[i]];
373 if ((k >= myFirstNode) && (myLastNode > k))
374 maskMyReducedNodes[k - myFirstNode] = i;
375 k = in->Nodes->globalDegreesOfFreedom[indexReducedNodes[i]];
376 if ((k >= myFirstDOF) && (myLastDOF > k))
377 {
378 maskMyReducedDOF[k - myFirstDOF] = i;
379 }
380 }
381 }
382 myNumReducedNodes = Dudley_Util_packMask(myNumNodes, maskMyReducedNodes, indexMyReducedNodes);
383 myNumReducedDOF = Dudley_Util_packMask(myNumDOF, maskMyReducedDOF, indexMyReducedDOF);
384
385 #ifdef ESYS_MPI
386 MPI_Allgather(&myNumReducedNodes, 1, MPI_INT, reduced_nodes_first_component, 1, MPI_INT,
387 in->Nodes->MPIInfo->comm);
388 MPI_Allgather(&myNumReducedDOF, 1, MPI_INT, reduced_dof_first_component, 1, MPI_INT,
389 in->Nodes->MPIInfo->comm);
390 #else
391 reduced_nodes_first_component[0] = myNumReducedNodes;
392 reduced_dof_first_component[0] = myNumReducedDOF;
393 #endif
394 globalNumReducedNodes = 0;
395 globalNumReducedDOF = 0;
396 for (i = 0; i < mpiSize; ++i)
397 {
398 k = reduced_nodes_first_component[i];
399 reduced_nodes_first_component[i] = globalNumReducedNodes;
400 globalNumReducedNodes += k;
401
402 k = reduced_dof_first_component[i];
403 reduced_dof_first_component[i] = globalNumReducedDOF;
404 globalNumReducedDOF += k;
405 }
406 reduced_nodes_first_component[mpiSize] = globalNumReducedNodes;
407 reduced_dof_first_component[mpiSize] = globalNumReducedDOF;
408 /* ==== distribution of Nodes =============================== */
409 in->Nodes->nodesDistribution.reset(new paso::Distribution(in->Nodes->MPIInfo, nodes_first_component, 1, 0));
410
411 /* ==== distribution of DOFs =============================== */
412 in->Nodes->degreesOfFreedomDistribution.reset(
413 new paso::Distribution(in->Nodes->MPIInfo, dof_first_component, 1, 0));
414
415 /* ==== distribution of reduced Nodes =============================== */
416 in->Nodes->reducedNodesDistribution.reset(new paso::Distribution(
417 in->Nodes->MPIInfo, reduced_nodes_first_component, 1, 0));
418
419 /* ==== distribution of reduced DOF =============================== */
420 in->Nodes->reducedDegreesOfFreedomDistribution.reset(
421 new paso::Distribution(in->Nodes->MPIInfo, reduced_dof_first_component, 1, 0));
422
423 delete[] maskMyReducedDOF;
424 delete[] indexMyReducedDOF;
425 delete[] maskMyReducedNodes;
426 delete[] indexMyReducedNodes;
427 delete[] reduced_dof_first_component;
428 delete[] reduced_nodes_first_component;
429
430 nodeMask = new index_t[in->Nodes->numNodes];
431 /* ==== nodes mapping which is a dummy structure ======== */
432 #pragma omp parallel for private(i) schedule(static)
433 for (i = 0; i < in->Nodes->numNodes; ++i)
434 nodeMask[i] = i;
435 in->Nodes->nodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
436
437 /* ==== mapping between nodes and reduced nodes ========== */
438 #pragma omp parallel for private(i) schedule(static)
439 for (i = 0; i < in->Nodes->numNodes; ++i)
440 nodeMask[i] = UNUSED;
441 #pragma omp parallel for private(i) schedule(static)
442 for (i = 0; i < numReducedNodes; ++i)
443 nodeMask[indexReducedNodes[i]] = i;
444 in->Nodes->reducedNodesMapping = Dudley_NodeMapping_alloc(in->Nodes->numNodes, nodeMask, UNUSED);
445 delete[] nodeMask;
446 /* ==== mapping between nodes and DOFs + DOF connector ========== */
447 Dudley_Mesh_createDOFMappingAndCoupling(in, false);
448 /* == mapping between nodes and reduced DOFs + reduced DOF connector == */
449 Dudley_Mesh_createDOFMappingAndCoupling(in, true);
450
451 /* get the Ids for DOFs and reduced nodes */
452 #pragma omp parallel private(i)
453 {
454 #pragma omp for
455 for (i = 0; i < in->Nodes->reducedNodesMapping->numTargets; ++i)
456 in->Nodes->reducedNodesId[i] = in->Nodes->Id[in->Nodes->reducedNodesMapping->map[i]];
457 #pragma omp for
458 for (i = 0; i < in->Nodes->degreesOfFreedomMapping->numTargets; ++i)
459 in->Nodes->degreesOfFreedomId[i] = in->Nodes->Id[in->Nodes->degreesOfFreedomMapping->map[i]];
460 #pragma omp for
461 for (i = 0; i < in->Nodes->reducedDegreesOfFreedomMapping->numTargets; ++i)
462 in->Nodes->reducedDegreesOfFreedomId[i] =
463 in->Nodes->Id[in->Nodes->reducedDegreesOfFreedomMapping->map[i]];
464 }
465 }
466
467 } // namespace dudley
468

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-6007 /trunk/ripley/test/python/dudley/src/Mesh_createNodeFileMappings.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26