/[escript]/trunk/finley/src/Assemble_CopyNodalData.cpp
ViewVC logotype

Contents of /trunk/finley/src/Assemble_CopyNodalData.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 18024 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18 #include "Assemble.h"
19 #include "Util.h"
20
21 namespace finley {
22
23 using escript::ValueError;
24 using escript::NotImplementedError;
25
26 template<typename Scalar>
27 void Assemble_CopyNodalData(const NodeFile* nodes, escript::Data& out,
28 const escript::Data& in)
29 {
30 if (!nodes)
31 return;
32
33 const int mpiSize = nodes->MPIInfo->size;
34 const int numComps = out.getDataPointSize();
35 const int in_data_type = in.getFunctionSpace().getTypeCode();
36 const int out_data_type = out.getFunctionSpace().getTypeCode();
37
38 // check out and in
39 if (numComps != in.getDataPointSize()) {
40 throw ValueError("Assemble_CopyNodalData: number of components of input and output Data do not match.");
41 } else if (!out.actsExpanded()) {
42 throw ValueError("Assemble_CopyNodalData: expanded Data object is expected for output data.");
43 } else if (in.isComplex() != out.isComplex()) {
44 throw ValueError("Assemble_CopyNodalData: complexity of input and output Data must match.");
45 }
46
47 // more sophisticated test needed for overlapping node/DOF counts
48 if (in_data_type == FINLEY_NODES) {
49 if (!in.numSamplesEqual(1, nodes->getNumNodes())) {
50 throw ValueError("Assemble_CopyNodalData: illegal number of samples of input Data object");
51 }
52 } else if (in_data_type == FINLEY_REDUCED_NODES) {
53 if (!in.numSamplesEqual(1, nodes->getNumReducedNodes())) {
54 throw ValueError("Assemble_CopyNodalData: illegal number of samples of input Data object");
55 }
56 } else if (in_data_type == FINLEY_DEGREES_OF_FREEDOM) {
57 if (!in.numSamplesEqual(1, nodes->getNumDegreesOfFreedom())) {
58 throw ValueError("Assemble_CopyNodalData: illegal number of samples of input Data object");
59 }
60 if (((out_data_type == FINLEY_NODES) || (out_data_type == FINLEY_DEGREES_OF_FREEDOM)) && !in.actsExpanded() && (mpiSize>1)) {
61 throw ValueError("Assemble_CopyNodalData: FINLEY_DEGREES_OF_FREEDOM to FINLEY_NODES or FINLEY_DEGREES_OF_FREEDOM requires expanded input data on more than one processor.");
62 }
63 } else if (in_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
64 if (!in.numSamplesEqual(1, nodes->getNumReducedDegreesOfFreedom())) {
65 throw ValueError("Assemble_CopyNodalData: illegal number of samples of input Data object");
66 }
67 if ((out_data_type == FINLEY_DEGREES_OF_FREEDOM) && !in.actsExpanded() && (mpiSize>1)) {
68 throw ValueError("Assemble_CopyNodalData: FINLEY_REDUCED_DEGREES_OF_FREEDOM to FINLEY_DEGREES_OF_FREEDOM requires expanded input data on more than one processor.");
69 }
70 } else {
71 throw ValueError( "Assemble_CopyNodalData: illegal function space type for target object");
72 }
73
74 dim_t numOut = 0;
75 switch (out_data_type) {
76 case FINLEY_NODES:
77 numOut = nodes->getNumNodes();
78 break;
79
80 case FINLEY_REDUCED_NODES:
81 numOut = nodes->getNumReducedNodes();
82 break;
83
84 case FINLEY_DEGREES_OF_FREEDOM:
85 numOut = nodes->getNumDegreesOfFreedom();
86 break;
87
88 case FINLEY_REDUCED_DEGREES_OF_FREEDOM:
89 numOut = nodes->getNumReducedDegreesOfFreedom();
90 break;
91
92 default:
93 throw ValueError("Assemble_CopyNodalData: illegal function space type for source object");
94 }
95
96 if (!out.numSamplesEqual(1, numOut)) {
97 throw ValueError("Assemble_CopyNodalData: illegal number of samples of output Data object");
98 }
99
100 const Scalar zero = static_cast<Scalar>(0);
101 const size_t numComps_size = numComps * sizeof(Scalar);
102
103 /**************************** FINLEY_NODES ******************************/
104 if (in_data_type == FINLEY_NODES) {
105 out.requireWrite();
106 if (out_data_type == FINLEY_NODES) {
107 #pragma omp parallel for
108 for (index_t n = 0; n < numOut; n++) {
109 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(n, zero), numComps_size);
110 }
111 } else if (out_data_type == FINLEY_REDUCED_NODES) {
112 const IndexVector& map = nodes->borrowReducedNodesTarget();
113 const dim_t mapSize = map.size();
114 #pragma omp parallel for
115 for (index_t n = 0; n < mapSize; n++) {
116 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(map[n], zero),
117 numComps_size);
118 }
119 } else if (out_data_type == FINLEY_DEGREES_OF_FREEDOM) {
120 const IndexVector& map = nodes->borrowDegreesOfFreedomTarget();
121 #pragma omp parallel for
122 for (index_t n = 0; n < numOut; n++) {
123 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(map[n], zero),
124 numComps_size);
125 }
126 } else if (out_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
127 const std::vector<index_t>& map = nodes->borrowReducedDegreesOfFreedomTarget();
128 #pragma omp parallel for
129 for (index_t n=0; n<numOut; n++) {
130 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(map[n], zero),
131 numComps_size);
132 }
133 }
134
135 /************************ FINLEY_REDUCED_NODES **************************/
136 } else if (in_data_type == FINLEY_REDUCED_NODES) {
137 if (out_data_type == FINLEY_NODES) {
138 throw ValueError("Assemble_CopyNodalData: cannot copy from reduced nodes to nodes.");
139 } else if (out_data_type == FINLEY_REDUCED_NODES) {
140 out.requireWrite();
141 const dim_t nNodes = nodes->getNumNodes();
142 #pragma omp parallel for
143 for (index_t n=0; n < nNodes; n++) {
144 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(n, zero), numComps_size);
145 }
146 } else if (out_data_type == FINLEY_DEGREES_OF_FREEDOM) {
147 throw ValueError("Assemble_CopyNodalData: cannot copy from reduced nodes to degrees of freedom.");
148 } else if (out_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
149 out.requireWrite();
150 const index_t* target = nodes->borrowTargetReducedNodes();
151 const IndexVector& map = nodes->borrowReducedDegreesOfFreedomTarget();
152 #pragma omp parallel for
153 for (index_t n = 0; n < numOut; n++) {
154 memcpy(out.getSampleDataRW(n, zero),
155 in.getSampleDataRO(target[map[n]], zero), numComps_size);
156 }
157 }
158 /********************** FINLEY_DEGREES_OF_FREEDOM ***********************/
159 } else if (in_data_type == FINLEY_DEGREES_OF_FREEDOM) {
160 out.requireWrite();
161 if (out_data_type == FINLEY_NODES) {
162 const_cast<escript::Data*>(&in)->resolve();
163 const index_t* target = nodes->borrowTargetDegreesOfFreedom();
164 #ifdef ESYS_HAVE_TRILINOS
165 using namespace esys_trilinos;
166
167 Teuchos::RCP<const MapType> colMap;
168 Teuchos::RCP<const MapType> rowMap;
169 MapType colPointMap;
170 MapType rowPointMap;
171 if (numComps > 1) {
172 colPointMap = BlockVectorType<Scalar>::makePointMap(
173 *nodes->trilinosColMap, numComps);
174 rowPointMap = BlockVectorType<Scalar>::makePointMap(
175 *nodes->trilinosRowMap, numComps);
176 colMap = Teuchos::rcpFromRef(colPointMap);
177 rowMap = Teuchos::rcpFromRef(rowPointMap);
178 } else {
179 colMap = nodes->trilinosColMap;
180 rowMap = nodes->trilinosRowMap;
181 }
182
183 const ImportType importer(rowMap, colMap);
184 const Teuchos::ArrayView<const Scalar> localIn(
185 in.getSampleDataRO(0, zero),
186 in.getNumDataPoints()*numComps);
187 Teuchos::RCP<VectorType<Scalar> > lclData =
188 rcp(new VectorType<Scalar>(rowMap, localIn, localIn.size(), 1));
189 Teuchos::RCP<VectorType<Scalar> > gblData = rcp(new VectorType<Scalar>(colMap, 1));
190 gblData->doImport(*lclData, importer, Tpetra::INSERT);
191 Teuchos::ArrayRCP<const Scalar> gblArray(gblData->getData(0));
192 #pragma omp parallel for
193 for (index_t i = 0; i < numOut; i++) {
194 const Scalar* src = &gblArray[target[i] * numComps];
195 std::copy(src, src+numComps, out.getSampleDataRW(i, zero));
196 }
197 #elif defined(ESYS_HAVE_PASO)
198 paso::Coupler_ptr<Scalar> coupler(new paso::Coupler<Scalar>(
199 nodes->degreesOfFreedomConnector, numComps, nodes->MPIInfo));
200 coupler->startCollect(in.getSampleDataRO(0, zero));
201 const Scalar* recv_buffer = coupler->finishCollect();
202 const index_t upperBound = nodes->getNumDegreesOfFreedom();
203 #pragma omp parallel for
204 for (index_t n = 0; n < numOut; n++) {
205 const index_t k = target[n];
206 if (k < upperBound) {
207 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(k, zero),
208 numComps_size);
209 } else {
210 memcpy(out.getSampleDataRW(n, zero),
211 &recv_buffer[(k - upperBound) * numComps],
212 numComps_size);
213 }
214 }
215 #endif // Trilinos / Paso
216 } else if (out_data_type == FINLEY_REDUCED_NODES) {
217 const_cast<escript::Data*>(&in)->resolve();
218 const index_t* target = nodes->borrowTargetDegreesOfFreedom();
219 const IndexVector& map = nodes->borrowReducedNodesTarget();
220 #ifdef ESYS_HAVE_TRILINOS
221 using namespace esys_trilinos;
222
223 Teuchos::RCP<const MapType> colMap;
224 Teuchos::RCP<const MapType> rowMap;
225 MapType colPointMap;
226 MapType rowPointMap;
227 if (numComps > 1) {
228 colPointMap = BlockVectorType<Scalar>::makePointMap(
229 *nodes->trilinosColMap, numComps);
230 rowPointMap = BlockVectorType<Scalar>::makePointMap(
231 *nodes->trilinosRowMap, numComps);
232 colMap = Teuchos::rcpFromRef(colPointMap);
233 rowMap = Teuchos::rcpFromRef(rowPointMap);
234 } else {
235 colMap = nodes->trilinosColMap;
236 rowMap = nodes->trilinosRowMap;
237 }
238
239 const ImportType importer(rowMap, colMap);
240 const Teuchos::ArrayView<const Scalar> localIn(
241 in.getSampleDataRO(0, zero),
242 in.getNumDataPoints()*numComps);
243 Teuchos::RCP<VectorType<Scalar> > lclData = rcp(new VectorType<Scalar>(
244 rowMap, localIn, localIn.size(), 1));
245 Teuchos::RCP<VectorType<Scalar> > gblData = rcp(new VectorType<Scalar>(colMap, 1));
246 gblData->doImport(*lclData, importer, Tpetra::INSERT);
247 Teuchos::ArrayRCP<const Scalar> gblArray(gblData->getData(0));
248 #pragma omp parallel for
249 for (index_t i = 0; i < numOut; i++) {
250 const Scalar* src = &gblArray[target[map[i]] * numComps];
251 std::copy(src, src+numComps, out.getSampleDataRW(i, zero));
252 }
253 #elif defined(ESYS_HAVE_PASO)
254 paso::Coupler_ptr<Scalar> coupler(new paso::Coupler<Scalar>(
255 nodes->degreesOfFreedomConnector, numComps, nodes->MPIInfo));
256 coupler->startCollect(in.getSampleDataRO(0, zero));
257 const Scalar* recv_buffer = coupler->finishCollect();
258 const index_t upperBound = nodes->getNumDegreesOfFreedom();
259 const dim_t mapSize = map.size();
260
261 #pragma omp parallel for
262 for (index_t n = 0; n < mapSize; n++) {
263 const index_t k = target[map[n]];
264 if (k < upperBound) {
265 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(k, zero),
266 numComps_size);
267 } else {
268 memcpy(out.getSampleDataRW(n, zero),
269 &recv_buffer[(k-upperBound)*numComps],
270 numComps_size);
271 }
272 }
273 #endif // Trilinos / Paso
274 } else if (out_data_type == FINLEY_DEGREES_OF_FREEDOM) {
275 #pragma omp parallel for
276 for (index_t n = 0; n < numOut; n++) {
277 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(n, zero),
278 numComps_size);
279 }
280 } else if (out_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
281 const index_t* target = nodes->borrowTargetDegreesOfFreedom();
282 const std::vector<index_t>& map = nodes->borrowReducedDegreesOfFreedomTarget();
283 #pragma omp parallel for
284 for (index_t n=0; n<numOut; n++) {
285 memcpy(out.getSampleDataRW(n, zero),
286 in.getSampleDataRO(target[map[n]], zero), numComps_size);
287 }
288 }
289
290 /****************** FINLEY_REDUCED_DEGREES_OF_FREEDOM *******************/
291 } else if (in_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
292 if (out_data_type == FINLEY_NODES) {
293 throw ValueError("Assemble_CopyNodalData: cannot copy from reduced degrees of freedom to nodes.");
294 } else if (out_data_type == FINLEY_REDUCED_NODES) {
295 const_cast<escript::Data*>(&in)->resolve();
296 const index_t* target = nodes->borrowTargetReducedDegreesOfFreedom();
297 const IndexVector& map = nodes->borrowReducedNodesTarget();
298 out.requireWrite();
299 #ifdef ESYS_HAVE_TRILINOS
300 using namespace esys_trilinos;
301
302 Teuchos::RCP<const MapType> colMap;
303 Teuchos::RCP<const MapType> rowMap;
304 MapType colPointMap;
305 MapType rowPointMap;
306 if (numComps > 1) {
307 colPointMap = BlockVectorType<Scalar>::makePointMap(
308 *nodes->trilinosReducedColMap, numComps);
309 rowPointMap = BlockVectorType<Scalar>::makePointMap(
310 *nodes->trilinosReducedRowMap, numComps);
311 colMap = Teuchos::rcpFromRef(colPointMap);
312 rowMap = Teuchos::rcpFromRef(rowPointMap);
313 } else {
314 colMap = nodes->trilinosReducedColMap;
315 rowMap = nodes->trilinosReducedRowMap;
316 }
317
318 const ImportType importer(rowMap, colMap);
319 const Teuchos::ArrayView<const Scalar> localIn(
320 in.getSampleDataRO(0, zero),
321 in.getNumDataPoints()*numComps);
322 Teuchos::RCP<VectorType<Scalar> > lclData = rcp(new VectorType<Scalar>(
323 rowMap, localIn, localIn.size(), 1));
324 Teuchos::RCP<VectorType<Scalar> > gblData = rcp(new VectorType<Scalar>(colMap, 1));
325 gblData->doImport(*lclData, importer, Tpetra::INSERT);
326 Teuchos::ArrayRCP<const Scalar> gblArray(gblData->getData(0));
327 #pragma omp parallel for
328 for (index_t i = 0; i < numOut; i++) {
329 const Scalar* src = &gblArray[target[map[i]] * numComps];
330 std::copy(src, src+numComps, out.getSampleDataRW(i, zero));
331 }
332 #elif defined(ESYS_HAVE_PASO)
333 paso::Coupler_ptr<Scalar> coupler(new paso::Coupler<Scalar>(
334 nodes->reducedDegreesOfFreedomConnector, numComps, nodes->MPIInfo));
335 coupler->startCollect(in.getSampleDataRO(0, zero));
336 const index_t upperBound = nodes->getNumReducedDegreesOfFreedom();
337 const dim_t mapSize = map.size();
338 const Scalar* recv_buffer = coupler->finishCollect();
339 #pragma omp parallel for
340 for (index_t n = 0; n < mapSize; n++) {
341 const index_t k = target[map[n]];
342 if (k < upperBound) {
343 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(k, zero),
344 numComps_size);
345 } else {
346 memcpy(out.getSampleDataRW(n, zero),
347 &recv_buffer[(k - upperBound) * numComps],
348 numComps_size);
349 }
350 }
351 #endif // Trilinos / Paso
352 } else if (out_data_type == FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
353 out.requireWrite();
354 #pragma omp parallel for
355 for (index_t n = 0; n < numOut; n++) {
356 memcpy(out.getSampleDataRW(n, zero), in.getSampleDataRO(n, zero), numComps_size);
357 }
358 } else if (out_data_type == FINLEY_DEGREES_OF_FREEDOM) {
359 throw ValueError("Assemble_CopyNodalData: cannot copy from reduced degrees of freedom to degrees of freedom.");
360 }
361 } // in_data_type
362 }
363
364 // instantiate our two supported versions
365 template void Assemble_CopyNodalData<escript::DataTypes::real_t>(
366 const NodeFile* nodes, escript::Data& out, const escript::Data& in);
367 template void Assemble_CopyNodalData<escript::DataTypes::cplx_t>(
368 const NodeFile* nodes, escript::Data& out, const escript::Data& in);
369
370 } // namespace finley
371

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Assemble_CopyNodalData.cpp:5567-5588 /branches/diaplayground/finley/src/Assemble_CopyNodalData.cpp:4940-5147 /branches/lapack2681/finley/src/Assemble_CopyNodalData.cpp:2682-2741 /branches/pasowrap/finley/src/Assemble_CopyNodalData.cpp:3661-3674 /branches/py3_attempt2/finley/src/Assemble_CopyNodalData.cpp:3871-3891 /branches/restext/finley/src/Assemble_CopyNodalData.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Assemble_CopyNodalData.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_CopyNodalData.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Assemble_CopyNodalData.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Assemble_CopyNodalData.cpp:5898-6118 /release/3.0/finley/src/Assemble_CopyNodalData.cpp:2591-2601 /release/4.0/finley/src/Assemble_CopyNodalData.cpp:5380-5406 /trunk/finley/src/Assemble_CopyNodalData.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26