1 /*
2 * =====================================================================================
3 *
4 * Filename: TempestOfflineMap.hpp
5 *
6 * Description: Interface to the TempestRemap library to compute the consistent,
7 * and accurate high-order conservative remapping weights for overlap
8 * grids on the sphere in climate simulations.
9 *
10 * Author: Vijay S. Mahadevan (vijaysm), mahadevan@anl.gov
11 *
12 * =====================================================================================
13 */
14
15 #include "Announce.h"
16 #include "DataMatrix3D.h"
17 #include "FiniteElementTools.h"
18 #include "SparseMatrix.h"
19 #include "STLStringHelper.h"
20
21 #include "moab/Remapping/TempestOfflineMap.hpp"
22 #include "DebugOutput.hpp"
23
24 #include <fstream>
25 #include <cmath>
26 #include <cstdlib>
27
28 ///////////////////////////////////////////////////////////////////////////////
29
30 // #define VERBOSE
31 // #define VVERBOSE
32
33 ///////////////////////////////////////////////////////////////////////////////
34
TempestOfflineMap(moab::TempestRemapper * remapper)35 moab::TempestOfflineMap::TempestOfflineMap ( moab::TempestRemapper* remapper ) : OfflineMap(), m_remapper ( remapper )
36 {
37 // Get the references for the MOAB core objects
38 mbCore = m_remapper->get_interface();
39 #ifdef MOAB_HAVE_MPI
40 pcomm = m_remapper->get_parallel_communicator();
41 #endif
42
43 // Update the references to the meshes
44 m_meshInput = remapper->GetMesh ( moab::Remapper::SourceMesh );
45 m_meshInputCov = remapper->GetCoveringMesh();
46 m_meshOutput = remapper->GetMesh ( moab::Remapper::TargetMesh );
47 m_meshOverlap = remapper->GetMesh ( moab::Remapper::IntersectedMesh );
48 m_globalMapAvailable = false;
49
50 is_parallel = false;
51 is_root = true;
52 rank = 0;
53 size = 1;
54 #ifdef MOAB_HAVE_MPI
55 int flagInit;
56 MPI_Initialized( &flagInit );
57 if (flagInit) {
58 is_parallel = true;
59 assert(pcomm != NULL);
60 rank = pcomm->rank();
61 size = pcomm->size();
62 is_root = (rank == 0);
63 }
64 #endif
65
66 moab::DebugOutput dbgprint ( std::cout, rank );
67
68 // Compute and store the total number of source and target DoFs corresponding
69 // to number of rows and columns in the mapping.
70
71 // Initialize dimension information from file
72 std::vector<std::string> dimNames(1);
73 std::vector<int> dimSizes(1);
74 dimNames[0] = "num_elem";
75
76 dbgprint.printf ( 0, "Initializing dimensions of map\n" );
77 dbgprint.printf ( 0, "Input mesh\n" );
78 dimSizes[0] = m_meshInputCov->faces.size();
79 this->InitializeSourceDimensions(dimNames, dimSizes);
80 dbgprint.printf ( 0, "Output mesh\n" );
81 dimSizes[0] = m_meshOutput->faces.size();
82 this->InitializeTargetDimensions(dimNames, dimSizes);
83
84 m_weightMapGlobal = NULL;
85
86 // Build a matrix of source and target discretization so that we know how to assign
87 // the global DoFs in parallel for the mapping weights
88 // For example, FV->FV: rows X cols = faces_source X faces_target
89 }
90
91 ///////////////////////////////////////////////////////////////////////////////
92
~TempestOfflineMap()93 moab::TempestOfflineMap::~TempestOfflineMap()
94 {
95 delete m_weightMapGlobal;
96 mbCore = NULL;
97 #ifdef MOAB_HAVE_MPI
98 pcomm = NULL;
99 #endif
100 m_meshInput = NULL;
101 m_meshOutput = NULL;
102 m_meshOverlap = NULL;
103 }
104
105 ///////////////////////////////////////////////////////////////////////////////
106
ParseVariableList(const std::string & strVariables,std::vector<std::string> & vecVariableStrings)107 static void ParseVariableList (
108 const std::string & strVariables,
109 std::vector< std::string > & vecVariableStrings
110 )
111 {
112 unsigned iVarBegin = 0;
113 unsigned iVarCurrent = 0;
114
115 // Parse variable name
116 for ( ;; )
117 {
118 if ( ( iVarCurrent >= strVariables.length() ) ||
119 ( strVariables[iVarCurrent] == ',' ) ||
120 ( strVariables[iVarCurrent] == ' ' )
121 )
122 {
123 if ( iVarCurrent == iVarBegin )
124 {
125 if ( iVarCurrent >= strVariables.length() )
126 {
127 break;
128 }
129 continue;
130 }
131
132 vecVariableStrings.push_back (
133 strVariables.substr ( iVarBegin, iVarCurrent - iVarBegin ) );
134
135 iVarBegin = iVarCurrent + 1;
136 }
137
138 iVarCurrent++;
139 }
140 }
141
142 ///////////////////////////////////////////////////////////////////////////////
143
SetDofMapTags(const std::string srcDofTagName,const std::string tgtDofTagName)144 moab::ErrorCode moab::TempestOfflineMap::SetDofMapTags(const std::string srcDofTagName, const std::string tgtDofTagName)
145 {
146 moab::ErrorCode rval;
147
148 rval = mbCore->tag_get_handle ( srcDofTagName.c_str(), m_nDofsPEl_Src*m_nDofsPEl_Src, MB_TYPE_INTEGER,
149 this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR(rval);
150 rval = mbCore->tag_get_handle ( tgtDofTagName.c_str(), m_nDofsPEl_Dest*m_nDofsPEl_Dest, MB_TYPE_INTEGER,
151 this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR(rval);
152
153 return moab::MB_SUCCESS;
154 }
155
156 ///////////////////////////////////////////////////////////////////////////////
157
SetDofMapAssociation(DiscretizationType srcType,bool isSrcContinuous,DataMatrix3D<int> * srcdataGLLNodes,DataMatrix3D<int> * srcdataGLLNodesSrc,DiscretizationType destType,bool isTgtContinuous,DataMatrix3D<int> * tgtdataGLLNodes)158 moab::ErrorCode moab::TempestOfflineMap::SetDofMapAssociation(DiscretizationType srcType, bool isSrcContinuous, DataMatrix3D<int>* srcdataGLLNodes, DataMatrix3D<int>* srcdataGLLNodesSrc,
159 DiscretizationType destType, bool isTgtContinuous, DataMatrix3D<int>* tgtdataGLLNodes)
160 {
161 moab::ErrorCode rval;
162 std::vector<bool> dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
163 std::vector<int> src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
164
165 // We are assuming that these are element based tags that are sized: np * np
166 m_srcDiscType = srcType;
167 m_destDiscType = destType;
168
169 bool vprint = is_root && false;
170
171 #ifdef VVERBOSE
172 {
173 src_soln_gdofs.resize(m_remapper->m_covering_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, -1);
174 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR(rval);
175 locsrc_soln_gdofs.resize(m_remapper->m_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src);
176 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR(rval);
177 tgt_soln_gdofs.resize(m_remapper->m_target_entities.size()*m_nDofsPEl_Dest*m_nDofsPEl_Dest);
178 rval = mbCore->tag_get_data ( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR(rval);
179
180 if (is_root)
181 {
182 {
183 std::ofstream output_file ( "sourcecov-gids-0.txt" );
184 output_file << "I, GDOF\n";
185 for (unsigned i=0; i < src_soln_gdofs.size(); ++i)
186 output_file << i << ", " << src_soln_gdofs[i] << "\n";
187
188 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
189 m_nTotDofs_SrcCov=0;
190 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
191 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
192 {
193 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
194 {
195 for ( int q = 0; q < m_nDofsPEl_Src; q++)
196 {
197 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1;
198 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
199 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
200 m_nTotDofs_SrcCov++;
201 dgll_cgll_covcol_ldofmap[ldof] = true;
202 }
203 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << m_nTotDofs_SrcCov << "\n";
204 }
205 }
206 }
207 output_file.flush(); // required here
208 output_file.close();
209 dgll_cgll_covcol_ldofmap.clear();
210 }
211
212 {
213 std::ofstream output_file ( "source-gids-0.txt" );
214 output_file << "I, GDOF\n";
215 for (unsigned i=0; i < locsrc_soln_gdofs.size(); ++i)
216 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
217
218 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
219 m_nTotDofs_Src=0;
220 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
221 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
222 {
223 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
224 {
225 for ( int q = 0; q < m_nDofsPEl_Src; q++)
226 {
227 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
228 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
229 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
230 m_nTotDofs_Src++;
231 dgll_cgll_col_ldofmap[ldof] = true;
232 }
233 output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << m_nTotDofs_Src << "\n";
234 }
235 }
236 }
237 output_file.flush(); // required here
238 output_file.close();
239 dgll_cgll_col_ldofmap.clear();
240 }
241
242 {
243 std::ofstream output_file ( "target-gids-0.txt" );
244 output_file << "I, GDOF\n";
245 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i)
246 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
247
248 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
249 m_nTotDofs_Dest=0;
250
251 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) {
252 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", " << m_nTotDofs_Dest << "\n";
253 m_nTotDofs_Dest++;
254 }
255
256 output_file.flush(); // required here
257 output_file.close();
258 }
259 }
260 else
261 {
262 {
263 std::ofstream output_file ( "sourcecov-gids-1.txt" );
264 output_file << "I, GDOF\n";
265 for (unsigned i=0; i < src_soln_gdofs.size(); ++i)
266 output_file << i << ", " << src_soln_gdofs[i] << "\n";
267
268 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
269 m_nTotDofs_SrcCov=0;
270 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
271 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
272 {
273 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
274 {
275 for ( int q = 0; q < m_nDofsPEl_Src; q++)
276 {
277 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1;
278 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
279 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
280 m_nTotDofs_SrcCov++;
281 dgll_cgll_covcol_ldofmap[ldof] = true;
282 }
283 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << m_nTotDofs_SrcCov << "\n";
284 }
285 }
286 }
287 output_file.flush(); // required here
288 output_file.close();
289 dgll_cgll_covcol_ldofmap.clear();
290 }
291
292 {
293 std::ofstream output_file ( "source-gids-1.txt" );
294 output_file << "I, GDOF\n";
295 for (unsigned i=0; i < locsrc_soln_gdofs.size(); ++i)
296 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
297
298 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
299 m_nTotDofs_Src=0;
300 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
301 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
302 {
303 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
304 {
305 for ( int q = 0; q < m_nDofsPEl_Src; q++)
306 {
307 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
308 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
309 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
310 m_nTotDofs_Src++;
311 dgll_cgll_col_ldofmap[ldof] = true;
312 }
313 output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << m_nTotDofs_Src << "\n";
314 }
315 }
316 }
317 output_file.flush(); // required here
318 output_file.close();
319 dgll_cgll_col_ldofmap.clear();
320 }
321
322 {
323 std::ofstream output_file ( "target-gids-1.txt" );
324 output_file << "I, GDOF\n";
325 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i)
326 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
327
328 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
329 m_nTotDofs_Dest=0;
330
331 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) {
332 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", " << m_nTotDofs_Dest << "\n";
333 m_nTotDofs_Dest++;
334 }
335
336 output_file.flush(); // required here
337 output_file.close();
338 }
339 }
340 }
341 #endif
342
343 // Now compute the mapping and store it for the covering mesh
344 col_dofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
345 col_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
346 col_gdofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
347 src_soln_gdofs.resize(m_remapper->m_covering_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX);
348 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR(rval);
349 m_nTotDofs_SrcCov = 0;
350 if (srcdataGLLNodes == NULL || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
351 for (unsigned i=0; i < col_dofmap.size(); ++i) {
352 col_dofmap[i] = i;
353 col_ldofmap[i] = i;
354 col_gdofmap[i] = src_soln_gdofs[i];
355 if (vprint) std::cout << "Col: " << i << ", " << src_soln_gdofs[i] << "\n";
356 m_nTotDofs_SrcCov++;
357 }
358 }
359 else {
360 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
361 // Put these remap coefficients into the SparseMatrix map
362 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
363 {
364 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
365 {
366 for ( int q = 0; q < m_nDofsPEl_Src; q++)
367 {
368 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1;
369 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
370 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
371 m_nTotDofs_SrcCov++;
372 dgll_cgll_covcol_ldofmap[ldof] = true;
373 }
374 if ( !isSrcContinuous ) m_nTotDofs_SrcCov++;
375 col_dofmap[ idof ] = ldof;
376 col_ldofmap[ ldof ] = idof;
377 assert(src_soln_gdofs[idof] > 0);
378 col_gdofmap[ idof ] = src_soln_gdofs[idof] - 1;
379 if (vprint) std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << col_gdofmap[idof] << ", " << m_nTotDofs_SrcCov << "\n";
380 }
381 }
382 }
383 }
384
385 srccol_dofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
386 srccol_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
387 srccol_gdofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
388 locsrc_soln_gdofs.resize(m_remapper->m_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX);
389 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR(rval);
390
391 // Now compute the mapping and store it for the original source mesh
392 m_nTotDofs_Src = 0;
393 if (srcdataGLLNodesSrc == NULL || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
394 for (unsigned i=0; i < srccol_dofmap.size(); ++i) {
395 srccol_dofmap[i] = i;
396 srccol_ldofmap[i] = i;
397 srccol_gdofmap[i] = locsrc_soln_gdofs[i];
398 m_nTotDofs_Src++;
399 }
400 }
401 else {
402 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize(m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false);
403 // Put these remap coefficients into the SparseMatrix map
404 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
405 {
406 for ( int p = 0; p < m_nDofsPEl_Src; p++ )
407 {
408 for ( int q = 0; q < m_nDofsPEl_Src; q++ )
409 {
410 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
411 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
412 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
413 m_nTotDofs_Src++;
414 dgll_cgll_col_ldofmap[ldof] = true;
415 }
416 if ( !isSrcContinuous ) m_nTotDofs_Src++;
417 srccol_dofmap[ idof ] = ldof;
418 srccol_ldofmap[ ldof ] = idof;
419 srccol_gdofmap[ idof ] = locsrc_soln_gdofs[idof] - 1;
420 }
421 }
422 }
423 }
424
425 row_dofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
426 row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
427 row_gdofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
428 tgt_soln_gdofs.resize(m_remapper->m_target_entities.size()*m_nDofsPEl_Dest*m_nDofsPEl_Dest, INT_MAX);
429 rval = mbCore->tag_get_data ( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR(rval);
430
431 // Now compute the mapping and store it for the target mesh
432 m_nTotDofs_Dest = 0;
433 if (tgtdataGLLNodes == NULL || destType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
434 for (unsigned i=0; i < row_dofmap.size(); ++i) {
435 row_dofmap[i] = i;
436 row_ldofmap[i] = i;
437 row_gdofmap[i] = tgt_soln_gdofs[i];
438 // if (vprint) std::cout << "Row: " << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << row_dofmap[i] << ", " << row_ldofmap[i] << ", " << tgt_soln_gdofs[i] << ", " << row_gdofmap[i] << "\n";
439 m_nTotDofs_Dest++;
440 }
441 }
442 else {
443 if (isTgtContinuous) dgll_cgll_row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, false);
444 // Put these remap coefficients into the SparseMatrix map
445 for ( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
446 {
447 for ( int p = 0; p < m_nDofsPEl_Dest; p++ )
448 {
449 for ( int q = 0; q < m_nDofsPEl_Dest; q++ )
450 {
451 const int ldof = (*tgtdataGLLNodes)[p][q][j] - 1;
452 const int idof = j * m_nDofsPEl_Dest * m_nDofsPEl_Dest + p * m_nDofsPEl_Dest + q;
453 if ( isTgtContinuous && !dgll_cgll_row_ldofmap[ldof] ) {
454 m_nTotDofs_Dest++;
455 dgll_cgll_row_ldofmap[ldof] = true;
456 }
457 if ( !isTgtContinuous ) m_nTotDofs_Dest++;
458 row_dofmap[ idof ] = ldof;
459 row_ldofmap[ ldof ] = idof;
460 row_gdofmap[ idof ] = tgt_soln_gdofs[idof] - 1;
461 if (vprint) std::cout << "Row: " << idof << ", " << ldof << ", " << tgt_soln_gdofs[idof] - 1 << "\n";
462 }
463 }
464 }
465 }
466
467 // Let us also allocate the local representation of the sparse matrix
468 #ifdef MOAB_HAVE_EIGEN
469 // if (vprint)
470 {
471 std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_dofmap.size() << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_dofmap.size() << "\n";
472 // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
473 }
474 #endif
475
476 return moab::MB_SUCCESS;
477 }
478
479 ///////////////////////////////////////////////////////////////////////////////
480
GenerateOfflineMap(std::string strInputType,std::string strOutputType,const int nPin,const int nPout,bool fBubble,int fMonotoneTypeID,bool fVolumetric,bool fNoConservation,bool fNoCheck,const std::string srcDofTagName,const std::string tgtDofTagName,const std::string strVariables,const std::string strInputData,const std::string strOutputData,const std::string strNColName,const bool fOutputDouble,const std::string strPreserveVariables,const bool fPreserveAll,const double dFillValueOverride,const bool fInputConcave,const bool fOutputConcave)481 moab::ErrorCode moab::TempestOfflineMap::GenerateOfflineMap ( std::string strInputType, std::string strOutputType,
482 const int nPin, const int nPout,
483 bool fBubble, int fMonotoneTypeID,
484 bool fVolumetric, bool fNoConservation, bool fNoCheck,
485 const std::string srcDofTagName, const std::string tgtDofTagName,
486 const std::string strVariables,
487 const std::string strInputData, const std::string strOutputData,
488 const std::string strNColName, const bool fOutputDouble,
489 const std::string strPreserveVariables, const bool fPreserveAll, const double dFillValueOverride,
490 const bool fInputConcave, const bool fOutputConcave )
491 {
492 NcError error ( NcError::silent_nonfatal );
493
494 moab::DebugOutput dbgprint ( std::cout, ( rank ) );
495 moab::ErrorCode rval;
496
497 try
498 {
499 // Check command line parameters (data arguments)
500 if ( ( strInputData != "" ) && ( strOutputData == "" ) )
501 {
502 _EXCEPTIONT ( "--in_data specified without --out_data" );
503 }
504 if ( ( strInputData == "" ) && ( strOutputData != "" ) )
505 {
506 _EXCEPTIONT ( "--out_data specified without --in_data" );
507 }
508
509 // Check command line parameters (data type arguments)
510 STLStringHelper::ToLower ( strInputType );
511 STLStringHelper::ToLower ( strOutputType );
512
513 DiscretizationType eInputType;
514 DiscretizationType eOutputType;
515
516 if ( strInputType == "fv" )
517 {
518 eInputType = DiscretizationType_FV;
519 }
520 else if ( strInputType == "cgll" )
521 {
522 eInputType = DiscretizationType_CGLL;
523 }
524 else if ( strInputType == "dgll" )
525 {
526 eInputType = DiscretizationType_DGLL;
527 }
528 else
529 {
530 _EXCEPTION1 ( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]",
531 strInputType.c_str() );
532 }
533
534 if ( strOutputType == "fv" )
535 {
536 eOutputType = DiscretizationType_FV;
537 }
538 else if ( strOutputType == "cgll" )
539 {
540 eOutputType = DiscretizationType_CGLL;
541 }
542 else if ( strOutputType == "dgll" )
543 {
544 eOutputType = DiscretizationType_DGLL;
545 }
546 else
547 {
548 _EXCEPTION1 ( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]",
549 strOutputType.c_str() );
550 }
551
552 // Monotonicity flags
553 int nMonotoneType = fMonotoneTypeID;
554
555 // Parse variable list
556 std::vector< std::string > vecVariableStrings;
557 ParseVariableList ( strVariables, vecVariableStrings );
558
559 // Parse preserve variable list
560 std::vector< std::string > vecPreserveVariableStrings;
561 ParseVariableList ( strPreserveVariables, vecPreserveVariableStrings );
562
563 if ( fPreserveAll && ( vecPreserveVariableStrings.size() != 0 ) )
564 {
565 _EXCEPTIONT ( "--preserveall and --preserve cannot both be specified" );
566 }
567
568 m_nDofsPEl_Src = nPin;
569 m_nDofsPEl_Dest = nPout;
570
571 rval = SetDofMapTags(srcDofTagName, tgtDofTagName);
572
573 // Calculate Face areas
574 if ( is_root ) dbgprint.printf ( 0, "Calculating input mesh Face areas\n" );
575 double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas(fInputConcave);
576 double dTotalAreaInput = dTotalAreaInput_loc;
577 #ifdef MOAB_HAVE_MPI
578 if (pcomm) MPI_Allreduce ( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
579 #endif
580 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
581
582 // Input mesh areas
583 m_meshInputCov->CalculateFaceAreas(fInputConcave);
584 if ( eInputType == DiscretizationType_FV )
585 {
586 this->SetSourceAreas ( m_meshInputCov->vecFaceArea );
587 }
588
589 // Calculate Face areas
590 if ( is_root ) dbgprint.printf ( 0, "Calculating output mesh Face areas\n" );
591 double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas(fOutputConcave);
592 double dTotalAreaOutput = dTotalAreaOutput_loc;
593 #ifdef MOAB_HAVE_MPI
594 if (pcomm) MPI_Allreduce ( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
595 #endif
596 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
597
598 // Output mesh areas
599 if ( eOutputType == DiscretizationType_FV )
600 {
601 this->SetTargetAreas ( m_meshOutput->vecFaceArea );
602 }
603
604 // Verify that overlap mesh is in the correct order
605 int ixSourceFaceMax = ( -1 );
606 int ixTargetFaceMax = ( -1 );
607
608 if ( m_meshOverlap->vecSourceFaceIx.size() !=
609 m_meshOverlap->vecTargetFaceIx.size()
610 )
611 {
612 _EXCEPTIONT ( "Invalid overlap mesh:\n"
613 " Possible mesh file corruption?" );
614 }
615
616 for ( unsigned i = 0; i < m_meshOverlap->vecSourceFaceIx.size(); i++ )
617 {
618 if ( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
619 {
620 ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;
621 }
622 if ( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
623 {
624 ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
625 }
626 }
627
628 /*
629 // Check for forward correspondence in overlap mesh
630 if ( // m_meshInputCov->faces.size() - ixSourceFaceMax == 0 //&&
631 ( m_meshOutput->faces.size() - ixTargetFaceMax == 0 )
632 )
633 {
634 if ( is_root ) dbgprint.printf ( 0, "Overlap mesh forward correspondence found\n" );
635 }
636 else if (
637 // m_meshOutput->faces.size() - ixSourceFaceMax == 0 //&&
638 ( m_meshInputCov->faces.size() - ixTargetFaceMax == 0 )
639 )
640 { // Check for reverse correspondence in overlap mesh
641 if ( is_root ) dbgprint.printf ( 0, "Overlap mesh reverse correspondence found (reversing)\n" );
642
643 // Reorder overlap mesh
644 m_meshOverlap->ExchangeFirstAndSecondMesh();
645 }
646 else
647 { // No correspondence found
648 _EXCEPTION4 ( "Invalid overlap mesh:\n"
649 " No correspondence found with input and output meshes (%i,%i) vs (%i,%i)",
650 m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax, ixTargetFaceMax );
651 }
652 */
653
654 // Calculate Face areas
655 if ( is_root ) dbgprint.printf ( 0, "Calculating overlap mesh Face areas\n" );
656 double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas(false);
657 double dTotalAreaOverlap = dTotalAreaOverlap_loc;
658 #ifdef MOAB_HAVE_MPI
659 if (pcomm) MPI_Allreduce ( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
660 #endif
661 if ( is_root ) dbgprint.printf ( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
662
663 // Partial cover
664 if ( fabs ( dTotalAreaOverlap - dTotalAreaInput ) > 1.0e-10 )
665 {
666 if ( !fNoCheck )
667 {
668 if ( is_root ) dbgprint.printf ( 0, "WARNING: Significant mismatch between overlap mesh area "
669 "and input mesh area.\n Automatically enabling --nocheck\n" );
670 fNoCheck = true;
671 }
672 }
673
674 /*
675 // Recalculate input mesh area from overlap mesh
676 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
677 dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
678 dbgprint.printf(0, "Recalculating source mesh areas\n");
679 dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
680 dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
681 }
682 */
683 // Finite volume input / Finite volume output
684 if ( ( eInputType == DiscretizationType_FV ) &&
685 ( eOutputType == DiscretizationType_FV )
686 )
687 {
688 // Generate reverse node array and edge map
689 m_meshInputCov->ConstructReverseNodeArray();
690 m_meshInputCov->ConstructEdgeMap();
691
692 // Initialize coordinates for map
693 this->InitializeSourceCoordinatesFromMeshFV ( *m_meshInputCov );
694 this->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput );
695
696 // Finite volume input / Finite element output
697 rval = this->SetDofMapAssociation(eInputType, false, NULL, NULL, eOutputType, false, NULL);MB_CHK_ERR(rval);
698
699 // Construct OfflineMap
700 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
701 LinearRemapFVtoFV_Tempest_MOAB ( nPin );
702 }
703 else if ( eInputType == DiscretizationType_FV )
704 {
705 DataMatrix3D<double> dataGLLJacobian;
706
707 if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data\n" );
708 double dNumericalArea_loc =
709 GenerateMetaData (
710 *m_meshOutput,
711 nPout,
712 fBubble,
713 dataGLLNodesDest,
714 dataGLLJacobian );
715
716 double dNumericalArea = dNumericalArea_loc;
717 #ifdef MOAB_HAVE_MPI
718 if (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
719 #endif
720 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
721
722 // Initialize coordinates for map
723 this->InitializeSourceCoordinatesFromMeshFV ( *m_meshInputCov );
724 this->InitializeTargetCoordinatesFromMeshFE (
725 *m_meshOutput, nPout, dataGLLNodesDest );
726
727 // Generate the continuous Jacobian
728 bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
729
730 if ( eOutputType == DiscretizationType_CGLL )
731 {
732 GenerateUniqueJacobian (
733 dataGLLNodesDest,
734 dataGLLJacobian,
735 this->GetTargetAreas() );
736 }
737 else
738 {
739 GenerateDiscontinuousJacobian (
740 dataGLLJacobian,
741 this->GetTargetAreas() );
742 }
743
744 // Generate reverse node array and edge map
745 m_meshInputCov->ConstructReverseNodeArray();
746 m_meshInputCov->ConstructEdgeMap();
747
748 // Finite volume input / Finite element output
749 rval = this->SetDofMapAssociation(eInputType, false, NULL, NULL,
750 eOutputType, (eOutputType == DiscretizationType_CGLL), &dataGLLNodesDest);MB_CHK_ERR(rval);
751
752 // Generate remap weights
753 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
754
755 if ( fVolumetric )
756 {
757 LinearRemapFVtoGLL_Volumetric_MOAB (
758 dataGLLNodesDest,
759 dataGLLJacobian,
760 this->GetTargetAreas(),
761 nPin,
762 nMonotoneType,
763 fContinuous,
764 fNoConservation );
765 }
766 else
767 {
768 LinearRemapFVtoGLL_MOAB (
769 dataGLLNodesDest,
770 dataGLLJacobian,
771 this->GetTargetAreas(),
772 nPin,
773 nMonotoneType,
774 fContinuous,
775 fNoConservation );
776 }
777 }
778 else if (
779 ( eInputType != DiscretizationType_FV ) &&
780 ( eOutputType == DiscretizationType_FV )
781 )
782 {
783 DataMatrix3D<double> dataGLLJacobianSrc, dataGLLJacobian;
784
785 if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data\n" );
786 // double dNumericalAreaCov_loc =
787 GenerateMetaData (
788 *m_meshInputCov,
789 nPin,
790 fBubble,
791 dataGLLNodesSrcCov,
792 dataGLLJacobian );
793
794 double dNumericalArea_loc =
795 GenerateMetaData (
796 *m_meshInput,
797 nPin,
798 fBubble,
799 dataGLLNodesSrc,
800 dataGLLJacobianSrc );
801
802 // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area: %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
803 // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);
804
805 double dNumericalArea = dNumericalArea_loc;
806 #ifdef MOAB_HAVE_MPI
807 if (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
808 #endif
809 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );
810
811 if ( fabs ( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
812 {
813 dbgprint.printf ( 0, "WARNING: Significant mismatch between input mesh "
814 "numerical area and geometric area\n" );
815 }
816
817 if ( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
818 {
819 _EXCEPTIONT ( "Number of element does not match between metadata and "
820 "input mesh" );
821 }
822
823 // Initialize coordinates for map
824 this->InitializeSourceCoordinatesFromMeshFE (
825 *m_meshInputCov, nPin, dataGLLNodesSrcCov );
826 this->InitializeSourceCoordinatesFromMeshFE (
827 *m_meshInput, nPin, dataGLLNodesSrc );
828 this->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput );
829
830 // Generate the continuous Jacobian for input mesh
831 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
832
833 if ( eInputType == DiscretizationType_CGLL )
834 {
835 GenerateUniqueJacobian (
836 dataGLLNodesSrcCov,
837 dataGLLJacobian,
838 this->GetSourceAreas() );
839 }
840 else
841 {
842 GenerateDiscontinuousJacobian (
843 dataGLLJacobian,
844 this->GetSourceAreas() );
845 }
846
847 // Finite element input / Finite volume output
848 rval = this->SetDofMapAssociation(eInputType, (eInputType == DiscretizationType_CGLL), &dataGLLNodesSrcCov, &dataGLLNodesSrc,
849 eOutputType, false, NULL);MB_CHK_ERR(rval);
850
851 // Generate offline map
852 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
853
854 if ( fVolumetric )
855 {
856 _EXCEPTIONT ( "Unimplemented: Volumetric currently unavailable for"
857 "GLL input mesh" );
858 }
859
860 LinearRemapSE4_Tempest_MOAB (
861 dataGLLNodesSrcCov,
862 dataGLLJacobian,
863 nMonotoneType,
864 fContinuousIn,
865 fNoConservation
866 );
867
868 }
869 else if (
870 ( eInputType != DiscretizationType_FV ) &&
871 ( eOutputType != DiscretizationType_FV )
872 )
873 {
874 DataMatrix3D<double> dataGLLJacobianIn, dataGLLJacobianSrc;
875 DataMatrix3D<double> dataGLLJacobianOut;
876
877 // Input metadata
878 if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data" );
879 double dNumericalAreaIn_loc =
880 GenerateMetaData (
881 *m_meshInputCov,
882 nPin,
883 fBubble,
884 dataGLLNodesSrcCov,
885 dataGLLJacobianIn );
886
887 double dNumericalAreaSrc_loc =
888 GenerateMetaData (
889 *m_meshInput,
890 nPin,
891 fBubble,
892 dataGLLNodesSrc,
893 dataGLLJacobianSrc );
894
895 assert(dNumericalAreaIn_loc >= dNumericalAreaSrc_loc);
896
897 double dNumericalAreaIn = dNumericalAreaSrc_loc;
898 #ifdef MOAB_HAVE_MPI
899 if (pcomm) MPI_Allreduce ( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
900 #endif
901 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Numerical Area: %1.15e", dNumericalAreaIn );
902
903 if ( fabs ( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
904 {
905 dbgprint.printf ( 0, "WARNING: Significant mismatch between input mesh "
906 "numerical area and geometric area" );
907 }
908
909 // Output metadata
910 if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data" );
911 double dNumericalAreaOut_loc =
912 GenerateMetaData (
913 *m_meshOutput,
914 nPout,
915 fBubble,
916 dataGLLNodesDest,
917 dataGLLJacobianOut );
918
919 double dNumericalAreaOut = dNumericalAreaOut_loc;
920 #ifdef MOAB_HAVE_MPI
921 if (pcomm) MPI_Allreduce ( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
922 #endif
923 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Numerical Area: %1.15e", dNumericalAreaOut );
924
925 if ( fabs ( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
926 {
927 if ( is_root ) dbgprint.printf ( 0, "WARNING: Significant mismatch between output mesh "
928 "numerical area and geometric area" );
929 }
930
931 // Initialize coordinates for map
932 this->InitializeSourceCoordinatesFromMeshFE (
933 *m_meshInputCov, nPin, dataGLLNodesSrcCov );
934 this->InitializeSourceCoordinatesFromMeshFE (
935 *m_meshInput, nPin, dataGLLNodesSrc );
936 this->InitializeTargetCoordinatesFromMeshFE (
937 *m_meshOutput, nPout, dataGLLNodesDest );
938
939 // Generate the continuous Jacobian for input mesh
940 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
941
942 if ( eInputType == DiscretizationType_CGLL )
943 {
944 GenerateUniqueJacobian (
945 dataGLLNodesSrcCov,
946 dataGLLJacobianIn,
947 this->GetSourceAreas() );
948 }
949 else
950 {
951 GenerateDiscontinuousJacobian (
952 dataGLLJacobianIn,
953 this->GetSourceAreas() );
954 }
955
956 // Generate the continuous Jacobian for output mesh
957 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
958
959 if ( eOutputType == DiscretizationType_CGLL )
960 {
961 GenerateUniqueJacobian (
962 dataGLLNodesDest,
963 dataGLLJacobianOut,
964 this->GetTargetAreas() );
965 }
966 else
967 {
968 GenerateDiscontinuousJacobian (
969 dataGLLJacobianOut,
970 this->GetTargetAreas() );
971 }
972
973 // Input Finite Element to Output Finite Element
974 rval = this->SetDofMapAssociation(eInputType, (eInputType == DiscretizationType_CGLL), &dataGLLNodesSrcCov, &dataGLLNodesSrc,
975 eOutputType, (eOutputType == DiscretizationType_CGLL), &dataGLLNodesDest);MB_CHK_ERR(rval);
976
977 // Generate offline map
978 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map" );
979
980 LinearRemapGLLtoGLL2_MOAB (
981 dataGLLNodesSrcCov,
982 dataGLLJacobianIn,
983 dataGLLNodesDest,
984 dataGLLJacobianOut,
985 this->GetTargetAreas(),
986 nPin,
987 nPout,
988 nMonotoneType,
989 fContinuousIn,
990 fContinuousOut,
991 fNoConservation
992 );
993
994 }
995 else
996 {
997 _EXCEPTIONT ( "Not implemented" );
998 }
999
1000 #ifdef MOAB_HAVE_EIGEN
1001 CopyTempestSparseMat_Eigen();
1002 #endif
1003
1004 // Verify consistency, conservation and monotonicity
1005 // gather weights to root process to perform consistency/conservation checks
1006 if ( !fNoCheck && false)
1007 {
1008 if ( !m_globalMapAvailable && size > 1 ) {
1009 // gather weights to root process to perform consistency/conservation checks
1010 rval = this->GatherAllToRoot();MB_CHK_ERR(rval);
1011 }
1012
1013 if ( is_root ) dbgprint.printf ( 0, "Verifying map" );
1014 this->IsConsistent ( 1.0e-8 );
1015 if ( !fNoConservation ) this->IsConservative ( 1.0e-8 );
1016
1017 if ( nMonotoneType != 0 )
1018 {
1019 this->IsMonotone ( 1.0e-12 );
1020 }
1021 }
1022
1023 // Apply Offline Map to data
1024 if ( strInputData != "" )
1025 {
1026 if ( is_root ) dbgprint.printf ( 0, "Applying offline map to data\n" );
1027
1028 this->SetFillValueOverride ( static_cast<float> ( dFillValueOverride ) );
1029 this->Apply (
1030 strInputData,
1031 strOutputData,
1032 vecVariableStrings,
1033 strNColName,
1034 fOutputDouble,
1035 false );
1036 }
1037
1038 // Copy variables from input file to output file
1039 if ( ( strInputData != "" ) && ( strOutputData != "" ) )
1040 {
1041 if ( fPreserveAll )
1042 {
1043 if ( is_root ) dbgprint.printf ( 0, "Preserving variables" );
1044 this->PreserveAllVariables ( strInputData, strOutputData );
1045
1046 }
1047 else if ( vecPreserveVariableStrings.size() != 0 )
1048 {
1049 if ( is_root ) dbgprint.printf ( 0, "Preserving variables" );
1050 this->PreserveVariables (
1051 strInputData,
1052 strOutputData,
1053 vecPreserveVariableStrings );
1054 }
1055 }
1056
1057 }
1058 catch ( Exception & e )
1059 {
1060 dbgprint.printf ( 0, "%s", e.ToString().c_str() );
1061 return ( moab::MB_FAILURE );
1062
1063 }
1064 catch ( ... )
1065 {
1066 return ( moab::MB_FAILURE );
1067 }
1068 return moab::MB_SUCCESS;
1069 }
1070
1071 ///////////////////////////////////////////////////////////////////////////////
1072 #ifdef MOAB_HAVE_EIGEN
1073
GetSourceGlobalNDofs()1074 int moab::TempestOfflineMap::GetSourceGlobalNDofs()
1075 {
1076 return m_weightMatrix.cols(); // return the global number of rows from the weight matrix
1077 }
1078
1079 // ///////////////////////////////////////////////////////////////////////////////
1080
GetDestinationGlobalNDofs()1081 int moab::TempestOfflineMap::GetDestinationGlobalNDofs()
1082 {
1083 return m_weightMatrix.rows(); // return the global number of columns from the weight matrix
1084 }
1085
1086 ///////////////////////////////////////////////////////////////////////////////
1087
GetSourceLocalNDofs()1088 int moab::TempestOfflineMap::GetSourceLocalNDofs()
1089 {
1090 return m_weightMatrix.cols(); // return the local number of rows from the weight matrix
1091 }
1092
1093 ///////////////////////////////////////////////////////////////////////////////
1094
GetDestinationLocalNDofs()1095 int moab::TempestOfflineMap::GetDestinationLocalNDofs()
1096 {
1097 return m_weightMatrix.rows(); // return the local number of columns from the weight matrix
1098 }
1099
1100 ///////////////////////////////////////////////////////////////////////////////
1101
GetSourceNDofsPerElement()1102 int moab::TempestOfflineMap::GetSourceNDofsPerElement()
1103 {
1104 return m_nDofsPEl_Src;
1105 }
1106
1107 ///////////////////////////////////////////////////////////////////////////////
1108
GetDestinationNDofsPerElement()1109 int moab::TempestOfflineMap::GetDestinationNDofsPerElement()
1110 {
1111 return m_nDofsPEl_Dest;
1112 }
1113
1114 ///////////////////////////////////////////////////////////////////////////////
1115
InitVectors()1116 void moab::TempestOfflineMap::InitVectors()
1117 {
1118 assert(m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0);
1119 m_rowVector.resize( m_weightMatrix.rows() );
1120 m_colVector.resize( m_weightMatrix.cols() );
1121 }
1122
1123
GetWeightMatrix()1124 moab::TempestOfflineMap::WeightMatrix& moab::TempestOfflineMap::GetWeightMatrix()
1125 {
1126 assert(m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0);
1127 return m_weightMatrix;
1128 }
1129
GetRowVector()1130 moab::TempestOfflineMap::WeightRowVector& moab::TempestOfflineMap::GetRowVector()
1131 {
1132 assert(m_rowVector.size() != 0);
1133 return m_rowVector;
1134 }
1135
GetColVector()1136 moab::TempestOfflineMap::WeightColVector& moab::TempestOfflineMap::GetColVector()
1137 {
1138 assert(m_colVector.size() != 0);
1139 return m_colVector;
1140 }
1141
1142 #endif
1143
1144 ///////////////////////////////////////////////////////////////////////////////
1145
IsConsistent(double dTolerance)1146 bool moab::TempestOfflineMap::IsConsistent (
1147 double dTolerance
1148 )
1149 {
1150 // Get map entries
1151 DataVector<int> dataRows;
1152 DataVector<int> dataCols;
1153 DataVector<double> dataEntries;
1154
1155 // Calculate row sums
1156 DataVector<double> dRowSums;
1157 if ( size > 1 )
1158 {
1159 if ( rank ) return true;
1160 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
1161 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
1162 dRowSums.Initialize ( m_mapRemapGlobal.GetRows() );
1163 }
1164 else
1165 {
1166 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
1167 dRowSums.Initialize ( m_mapRemap.GetRows() );
1168 }
1169
1170 for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
1171 {
1172 dRowSums[dataRows[i]] += dataEntries[i];
1173 }
1174
1175 // Verify all row sums are equal to 1
1176 bool fConsistent = true;
1177 for ( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1178 {
1179 if ( fabs ( dRowSums[i] - 1.0 ) > dTolerance )
1180 {
1181 fConsistent = false;
1182 Announce ( "TempestOfflineMap is not consistent in row %i (%1.15e)",
1183 i, dRowSums[i] );
1184 }
1185 }
1186
1187 return fConsistent;
1188 }
1189
1190 ///////////////////////////////////////////////////////////////////////////////
1191
IsConservative(double dTolerance)1192 bool moab::TempestOfflineMap::IsConservative (
1193 double dTolerance
1194 )
1195 {
1196 // Get map entries
1197 DataVector<int> dataRows;
1198 DataVector<int> dataCols;
1199 DataVector<double> dataEntries;
1200 const DataVector<double>& dTargetAreas = this->GetGlobalTargetAreas();
1201 const DataVector<double>& dSourceAreas = this->GetGlobalSourceAreas();
1202
1203 // Calculate column sums
1204 DataVector<double> dColumnSums;
1205 if ( size > 1 )
1206 {
1207 if ( rank ) return true;
1208 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
1209 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
1210 dColumnSums.Initialize ( m_mapRemapGlobal.GetColumns() );
1211 }
1212 else
1213 {
1214 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
1215 dColumnSums.Initialize ( m_mapRemap.GetColumns() );
1216 }
1217
1218 for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
1219 {
1220 dColumnSums[dataCols[i]] +=
1221 dataEntries[i] * dTargetAreas[dataRows[i]];
1222 }
1223
1224 // Verify all column sums equal the input Jacobian
1225 bool fConservative = true;
1226 for ( unsigned i = 0; i < dColumnSums.GetRows(); i++ )
1227 {
1228 if ( fabs ( dColumnSums[i] - dSourceAreas[i] ) > dTolerance )
1229 {
1230 fConservative = false;
1231 Announce ( "TempestOfflineMap is not conservative in column "
1232 "%i (%1.15e / %1.15e)",
1233 i, dColumnSums[i], dSourceAreas[i] );
1234 }
1235 }
1236
1237 return fConservative;
1238 }
1239
1240 ///////////////////////////////////////////////////////////////////////////////
1241
IsMonotone(double dTolerance)1242 bool moab::TempestOfflineMap::IsMonotone (
1243 double dTolerance
1244 )
1245 {
1246 // Get map entries
1247 DataVector<int> dataRows;
1248 DataVector<int> dataCols;
1249 DataVector<double> dataEntries;
1250
1251 if ( size > 1 )
1252 {
1253 if ( rank ) return true;
1254 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
1255 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
1256 }
1257 else
1258 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
1259
1260 // Verify all entries are in the range [0,1]
1261 bool fMonotone = true;
1262 for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
1263 {
1264 if ( ( dataEntries[i] < -dTolerance ) ||
1265 ( dataEntries[i] > 1.0 + dTolerance )
1266 )
1267 {
1268 fMonotone = false;
1269
1270 Announce ( "TempestOfflineMap is not monotone in entry (%i): %1.15e",
1271 i, dataEntries[i] );
1272 }
1273 }
1274
1275 return fMonotone;
1276 }
1277
1278
1279 ///////////////////////////////////////////////////////////////////////////////
1280 #ifdef MOAB_HAVE_MPI
GatherAllToRoot()1281 moab::ErrorCode moab::TempestOfflineMap::GatherAllToRoot() // Collective
1282 {
1283 Mesh globalMesh;
1284 int ierr, rootProc = 0, nprocs = size;
1285 moab::ErrorCode rval;
1286
1287 // Write SparseMatrix entries
1288 DataVector<int> vecRow;
1289 DataVector<int> vecCol;
1290 DataVector<double> vecS;
1291
1292 moab::DebugOutput dbgprint ( std::cout, ( rank ) );
1293
1294 m_mapRemap.GetEntries ( vecRow, vecCol, vecS );
1295 const DataVector<double>& dOrigSourceAreas = m_meshInput->vecFaceArea;
1296 const DataVector<double>& dSourceAreas = m_meshInputCov->vecFaceArea;
1297 const DataVector<double>& dTargetAreas = m_meshOutput->vecFaceArea;
1298
1299 // Translate the index in Row and Col to global_id and dump it out
1300
1301 // Communicate the necessary data
1302 const int NDATA = 7;
1303 int gnnz = 0, gsrc = 0, gsrccov = 0, gtar = 0, gsrcdofs = 0, gsrccovdofs = 0, gtgtdofs = 0;
1304 std::vector<int> rootSizesData, rowcolsv;
1305 DataVector<int> rows, cols, srcelmindx, tgtelmindx;
1306 {
1307 // First, accumulate the sizes of rows and columns of the matrix
1308 if ( !rank ) rootSizesData.resize ( size * NDATA );
1309
1310 int sendarray[NDATA];
1311 sendarray[0] = vecS.GetRows();
1312 sendarray[1] = dOrigSourceAreas.GetRows();
1313 sendarray[2] = dTargetAreas.GetRows();
1314 sendarray[3] = dSourceAreas.GetRows();
1315 sendarray[4] = m_nTotDofs_Src;
1316 sendarray[5] = m_nTotDofs_Dest;
1317 sendarray[6] = m_nTotDofs_SrcCov;
1318
1319 ierr = MPI_Gather ( sendarray, NDATA, MPI_INTEGER, rootSizesData.data(), NDATA, MPI_INTEGER, rootProc, pcomm->comm() );
1320 if ( ierr != MPI_SUCCESS ) return moab::MB_FAILURE;
1321
1322 if ( !rank )
1323 {
1324 for ( int i = 0; i < nprocs; ++i )
1325 {
1326 unsigned offset = NDATA * i;
1327 gnnz += rootSizesData[offset];
1328 gsrc += rootSizesData[offset + 1];
1329 gtar += rootSizesData[offset + 2];
1330 gsrccov += rootSizesData[offset + 3];
1331 gsrcdofs += rootSizesData[offset + 4];
1332 gtgtdofs += rootSizesData[offset + 5];
1333 gsrccovdofs += rootSizesData[offset + 6];
1334 }
1335 rowcolsv.resize ( 2 * gnnz + gsrc + gtar );
1336 rows.Initialize ( gnnz ); // we are assuming rows = cols
1337 cols.Initialize ( gnnz ); // we are assuming rows = cols
1338
1339 // Let us allocate our global offline map object
1340 m_weightMapGlobal = new OfflineMap();
1341 std::vector<std::string> dimNames(1);
1342 std::vector<int> dimSizes(1);
1343 dimNames[0] = "num_elem";
1344 dimSizes[0] = gsrc;
1345 m_weightMapGlobal->InitializeSourceDimensions(dimNames, dimSizes);
1346 dimSizes[0] = gtar;
1347 m_weightMapGlobal->InitializeTargetDimensions(dimNames, dimSizes);
1348
1349 /*
1350 VSM: do we need to gather the entire mesh ?!?!
1351 The following initialization can't work correctly without the mesh on the root process
1352 */
1353 if (m_srcDiscType == DiscretizationType_FV) /* unsure if we actually care about the type for this */
1354 m_weightMapGlobal->InitializeSourceCoordinatesFromMeshFV ( *m_meshInput );
1355 else
1356 m_weightMapGlobal->InitializeSourceCoordinatesFromMeshFE ( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc );
1357
1358 if (m_destDiscType == DiscretizationType_FV) /* unsure if we actually care about the type for this */
1359 m_weightMapGlobal->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput );
1360 else
1361 m_weightMapGlobal->InitializeTargetCoordinatesFromMeshFE ( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest );
1362
1363 DataVector<double>& m_areasSrcGlobal = m_weightMapGlobal->GetSourceAreas();
1364 m_areasSrcGlobal.Initialize ( gsrcdofs ); srcelmindx.Initialize ( gsrc );
1365 DataVector<double>& m_areasTgtGlobal = m_weightMapGlobal->GetTargetAreas();
1366 m_areasTgtGlobal.Initialize ( gtgtdofs ); tgtelmindx.Initialize ( gtar );
1367
1368 #ifdef VERBOSE
1369 dbgprint.printf ( 0, "Received global dimensions: %d, %d\n", vecRow.GetRows(), rows.GetRows() );
1370 dbgprint.printf ( 0, "Global: n(source) = %d, n(srccov) = %d, and n(target) = %d\n", gsrc, gsrccov, gtar );
1371 dbgprint.printf ( 0, "Operator size = %d X %d and NNZ = %d\n", gsrcdofs, gtgtdofs, gnnz );
1372 #endif
1373 }
1374 }
1375
1376 #ifdef VERBOSE
1377 {
1378 std::stringstream sstr;
1379 sstr << "rowscols_" << rank << ".txt";
1380 std::ofstream output_file ( sstr.str().c_str() );
1381 output_file << "VALUES\n";
1382 for ( unsigned ip = 0; ip < vecS.GetRows(); ++ip )
1383 {
1384 output_file << ip << " (" << row_gdofmap[row_ldofmap[vecRow[ip]]] << ", " << col_gdofmap[col_ldofmap[vecCol[ip]]] << ") = " << vecS[ip] << "\n";
1385
1386 }
1387 output_file.flush(); // required here
1388 output_file.close();
1389 }
1390 #endif
1391
1392 {
1393 // Next, accumulate the row and column values for the matrices (in local indexing)
1394 // Let's gather the data in the following order:
1395 // 1) row indices,
1396 // 2) column indices,
1397 // 3) GID for source elements,
1398 // 4) GID for target elements
1399 //
1400 const int nR = 2 * vecS.GetRows() + dOrigSourceAreas.GetRows() + dTargetAreas.GetRows();
1401 std::vector<int> sendarray ( nR );
1402 for ( unsigned ix = 0; ix < vecS.GetRows(); ++ix )
1403 {
1404 sendarray[ix] = row_gdofmap[row_ldofmap[vecRow[ix]]];
1405 }
1406 for ( unsigned ix = 0, offset = vecS.GetRows(); ix < vecS.GetRows(); ++ix )
1407 {
1408 sendarray[offset + ix] = col_gdofmap[col_ldofmap[vecCol[ix]]];
1409 }
1410
1411 {
1412 moab::Tag gidtag;
1413 rval = mbCore->tag_get_handle ( "GLOBAL_ID", gidtag ); MB_CHK_ERR ( rval );
1414 rval = mbCore->tag_get_data ( gidtag, m_remapper->m_source_entities, &sendarray[2 * vecS.GetRows()] ); MB_CHK_ERR ( rval );
1415 rval = mbCore->tag_get_data ( gidtag, m_remapper->m_target_entities, &sendarray[2 * vecS.GetRows() + dOrigSourceAreas.GetRows()] ); MB_CHK_ERR ( rval );
1416 }
1417
1418 std::vector<int> displs, rcount;
1419 if ( !rank )
1420 {
1421 displs.resize ( size, 0 );
1422 rcount.resize ( size, 0 );
1423 int gsum = 0;
1424 for ( int i = 0; i < nprocs; ++i )
1425 {
1426 displs[i] = gsum;
1427 rcount[i] = 2 * rootSizesData[NDATA * i] + rootSizesData[NDATA * i + 1] + rootSizesData[NDATA * i + 2];
1428 /* + rootSizesData[NDATA * i + 3] + rootSizesData[NDATA * i + 4] */
1429 gsum += rcount[i];
1430 }
1431 assert ( rowcolsv.size() - gsum == 0 );
1432 }
1433
1434 // Both rows and columns have a size of "rowsize"
1435 ierr = MPI_Gatherv ( &sendarray[0], nR, MPI_INTEGER, &rowcolsv[0], &rcount[0], &displs[0], MPI_INTEGER, rootProc, pcomm->comm() );
1436
1437 if ( !rank )
1438 {
1439 #ifdef VERBOSE
1440 std::ofstream output_file ( "rows-cols.txt", std::ios::out );
1441 output_file << "ROWS\n";
1442 #endif
1443 for ( int ip = 0, offset = 0; ip < nprocs; ++ip )
1444 {
1445 int istart = displs[ip], iend = istart + rootSizesData[NDATA * ip];
1446 for ( int i = istart; i < iend; ++i, ++offset )
1447 {
1448 rows[offset] = rowcolsv[i];
1449 #ifdef VERBOSE
1450 output_file << offset << " " << rows[offset] << "\n";
1451 #endif
1452 }
1453 }
1454 #ifdef VERBOSE
1455 output_file << "COLS\n";
1456 #endif
1457 for ( int ip = 0, offset = 0; ip < nprocs; ++ip )
1458 {
1459 int istart = displs[ip] + rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip];
1460 for ( int i = istart; i < iend; ++i, ++offset )
1461 {
1462 cols[offset] = rowcolsv[i];
1463 #ifdef VERBOSE
1464 output_file << offset << " " << cols[offset] << "\n";
1465 #endif
1466 }
1467 }
1468 #ifdef VERBOSE
1469 output_file.flush(); // required here
1470 output_file.close();
1471 #endif
1472 for ( int ip = 0, offset = 0; ip < nprocs; ++ip )
1473 {
1474 int istart = displs[ip] + 2 * rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip + 1];
1475 for ( int i = istart; i < iend; ++i, ++offset )
1476 {
1477 srcelmindx[offset] = rowcolsv[i];
1478 }
1479 }
1480 for ( int ip = 0, offset = 0; ip < nprocs; ++ip )
1481 {
1482 int istart = displs[ip] + 2 * rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1], iend = istart + rootSizesData[NDATA * ip + 2];
1483 for ( int i = istart; i < iend; ++i, ++offset )
1484 {
1485 tgtelmindx[offset] = rowcolsv[i];
1486 }
1487 }
1488 } /* (!rank) */
1489 rowcolsv.clear();
1490 }
1491
1492 int nSc = m_dSourceVertexLon.GetColumns();
1493 int nTc = m_dTargetVertexLon.GetColumns();
1494 {
1495 // Next, accumulate the row and column values for the matrices (in local indexing)
1496 // Let's gather the data in the following order:
1497 // 1) matrix weights,
1498 // 2) source element areas,
1499 // 3) target element areas,
1500 // 4) source element: m_dSourceCenterLon
1501 // 5) source element: m_dSourceCenterLat
1502 // 6) target element: m_dTargetCenterLon
1503 // 7) target element: m_dTargetCenterLat
1504 //
1505 const int nR = vecS.GetRows() + dOrigSourceAreas.GetRows() + dTargetAreas.GetRows() + 2*m_nTotDofs_Src*(1+nSc) + 2*m_nTotDofs_Dest*(1+nTc);
1506
1507 std::vector<double> sendarray ( nR );
1508 int locoffset = 0;
1509 std::copy ( ( double* ) vecS, ( double* ) vecS + vecS.GetRows(), sendarray.begin() + locoffset); locoffset += vecS.GetRows();
1510 std::copy ( ( const double* ) dOrigSourceAreas, ( const double* ) dOrigSourceAreas + dOrigSourceAreas.GetRows(), sendarray.begin() + locoffset ); locoffset += dOrigSourceAreas.GetRows();
1511 std::copy ( ( const double* ) dTargetAreas, ( const double* ) dTargetAreas + dTargetAreas.GetRows(), sendarray.begin() + locoffset ); locoffset += dTargetAreas.GetRows();
1512
1513 std::cout << "[" << rank << "] " << m_nTotDofs_Src << ", m_dSourceCenterLon.size() = " << m_dSourceCenterLon.GetRows() << " and " << m_nTotDofs_Dest << ", m_dTargetCenterLon.size() = " << m_dTargetCenterLon.GetRows() << "\n";
1514
1515 std::copy ( ( const double* ) m_dSourceCenterLon, ( const double* ) m_dSourceCenterLon + m_dSourceCenterLon.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dSourceCenterLon.GetRows();
1516 std::copy ( ( const double* ) m_dSourceCenterLat, ( const double* ) m_dSourceCenterLat + m_dSourceCenterLat.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dSourceCenterLat.GetRows();
1517
1518 double** dSCLon = m_dSourceVertexLon;
1519 std::copy ( &dSCLon[0][0], &dSCLon[0][0] + m_dSourceVertexLon.GetRows()*nSc, sendarray.begin() + locoffset ); locoffset += m_dSourceVertexLon.GetRows()*nSc;
1520 double** dSCLat = m_dSourceVertexLat;
1521 std::copy ( &dSCLat[0][0], &dSCLat[0][0] + m_dSourceVertexLat.GetRows()*nSc, sendarray.begin() + locoffset ); locoffset += m_dSourceVertexLat.GetRows()*nSc;
1522
1523 std::copy ( ( const double* ) m_dTargetCenterLon, ( const double* ) m_dTargetCenterLon + m_dTargetCenterLon.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dTargetCenterLon.GetRows();
1524 std::copy ( ( const double* ) m_dTargetCenterLat, ( const double* ) m_dTargetCenterLat + m_dTargetCenterLat.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dTargetCenterLat.GetRows();
1525
1526 double** dTCLon = m_dTargetVertexLon;
1527 std::copy ( &dTCLon[0][0], &dTCLon[0][0] + m_dTargetVertexLon.GetRows()*nTc, sendarray.begin() + locoffset ); locoffset += m_dTargetVertexLon.GetRows()*nTc;
1528 double** dTCLat = m_dTargetVertexLat;
1529 std::copy ( &dTCLat[0][0], &dTCLat[0][0] + m_dTargetVertexLat.GetRows()*nTc, sendarray.begin() + locoffset ); locoffset += m_dTargetVertexLat.GetRows()*nTc;
1530
1531 std::vector<int> displs, rcount;
1532 int gsum = 0;
1533 if ( !rank )
1534 {
1535 displs.resize ( size, 0 );
1536 rcount.resize ( size, 0 );
1537 for ( int i = 0; i < nprocs; ++i )
1538 {
1539 displs[i] = gsum;
1540 rcount[i] = rootSizesData[NDATA * i] + rootSizesData[NDATA * i + 1] + rootSizesData[NDATA * i + 2] +
1541 2 * rootSizesData[NDATA * i + 4] * (1 + nSc) +
1542 2 * rootSizesData[NDATA * i + 5] * (1 + nTc);
1543 gsum += rcount[i];
1544 }
1545 }
1546
1547 std::vector<double> rowcolsvals ( (!rank ? gnnz + gsrc + gtar + 2*gsrcdofs*(1+nSc) + 2*gtgtdofs*(1+nTc) : 0 ) );
1548 // Both rows and columns have a size of "rowsize"
1549 ierr = MPI_Gatherv ( sendarray.data(), sendarray.size(), MPI_DOUBLE, rowcolsvals.data(), &rcount[0], &displs[0], MPI_DOUBLE, rootProc, pcomm->comm() );
1550
1551 if ( !rank )
1552 {
1553 SparseMatrix<double>& spmat = m_weightMapGlobal->GetSparseMatrix();
1554 {
1555 #ifdef VERBOSE
1556 std::ofstream output_file ( "rows-cols.txt", std::ios::app );
1557 output_file << "VALUES\n";
1558 #endif
1559 for ( int ip = 0, offset = 0; ip < nprocs; ++ip )
1560 {
1561 for ( int i = displs[ip]; i < displs[ip] + rootSizesData[NDATA * ip]; ++i, ++offset )
1562 {
1563 // globvalues[offset] = rowcolsvals[i];
1564 spmat(rows[offset], cols[offset]) += rowcolsvals[i];
1565 #ifdef VERBOSE
1566 output_file << offset << " (" << rows[offset] << ", " << cols[offset] << ") = " << rowcolsvals[i] << "\n";
1567 #endif
1568 }
1569 }
1570 #ifdef VERBOSE
1571 output_file.flush(); // required here
1572 output_file.close();
1573 #endif
1574
1575 }
1576 // m_mapRemapGlobal.SetEntries(rows, cols, globvalues);
1577 // m_weightMapGlobal->GetSparseMatrix().AddEntries ( rows, cols, globvalues );
1578
1579 {
1580 DataVector<double>& m_areasSrcGlobal = m_weightMapGlobal->GetSourceAreas();
1581 DataVector<double>& m_areasTgtGlobal = m_weightMapGlobal->GetTargetAreas();
1582 // Store the global source and target elements areas
1583 #ifdef VERBOSE
1584 std::ofstream output_file ( "source-target-areas.txt" );
1585 output_file << "Source areas (nelems = " << gsrc << ")\n";
1586 #endif
1587 int offset = 0;
1588 for ( int ip = 0; ip < nprocs; ++ip )
1589 {
1590 int istart = displs[ip] + rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip + 1];
1591 for ( int i = istart; i < iend; ++i, ++offset )
1592 {
1593 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs);
1594 m_areasSrcGlobal[srcelmindx[offset]] = rowcolsvals[i];
1595 }
1596 }
1597
1598 #ifdef VERBOSE
1599 for ( unsigned i = 0; i < m_areasSrcGlobal.GetRows(); ++i )
1600 output_file << srcelmindx[i] << " " << m_areasSrcGlobal[i] << "\n";
1601 #endif
1602
1603 offset = 0;
1604 for ( int ip = 0; ip < nprocs; ++ip )
1605 {
1606 int istart = displs[ip] + rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1], iend = istart + rootSizesData[NDATA * ip + 2];
1607 for ( int i = istart; i < iend; ++i, ++offset )
1608 {
1609 // assert(offset < gtgtdofs && tgtelmindx[offset] < gtgtdofs);
1610 m_areasTgtGlobal[tgtelmindx[offset]] = rowcolsvals[i];
1611 }
1612 }
1613 #ifdef VERBOSE
1614 output_file << "Target areas (nelems = " << gtar << ")\n";
1615 for ( unsigned i = 0; i < m_areasTgtGlobal.GetRows(); ++i )
1616 output_file << tgtelmindx[i] << " " << m_areasTgtGlobal[i] << "\n";
1617 #endif
1618 #ifdef VERBOSE
1619 output_file.flush(); // required here
1620 output_file.close();
1621 #endif
1622 }
1623
1624 DataVector<double>& dSourceCenterLon = m_weightMapGlobal->GetSourceCenterLon();
1625 DataVector<double>& dSourceCenterLat = m_weightMapGlobal->GetSourceCenterLat();
1626 dSourceCenterLon.Initialize(gsrcdofs);
1627 dSourceCenterLat.Initialize(gsrcdofs);
1628 {
1629 int offset = 0;
1630 for ( int ip = 0; ip < nprocs; ++ip )
1631 {
1632 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2];
1633 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4];
1634 for ( int i = istart; i < iend; ++i, ++offset )
1635 {
1636 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs);
1637 // dSourceCenterLon[srcelmindx[offset]] = rowcolsvals[i];
1638 dSourceCenterLon[offset] = rowcolsvals[i];
1639 }
1640 }
1641
1642 offset = 0;
1643 for ( int ip = 0; ip < nprocs; ++ip )
1644 {
1645 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + rootSizesData[NDATA * ip + 4];
1646 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4];
1647 for ( int i = istart; i < iend; ++i, ++offset )
1648 {
1649 // dSourceCenterLat[srcelmindx[offset]] = rowcolsvals[i];
1650 dSourceCenterLat[offset] = rowcolsvals[i];
1651 }
1652 }
1653 }
1654
1655
1656 DataMatrix<double>& dSourceVertexLon = m_weightMapGlobal->GetSourceVertexLon();
1657 DataMatrix<double>& dSourceVertexLat = m_weightMapGlobal->GetSourceVertexLat();
1658 dSourceVertexLon.Initialize(gsrcdofs, nSc);
1659 dSourceVertexLat.Initialize(gsrcdofs, nSc);
1660 {
1661 int ioffset = 0;
1662 for ( int ip = 0; ip < nprocs; ++ip )
1663 {
1664 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4];
1665 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4] * nSc;
1666 for ( int i = istart; i < iend; ++i, ++ioffset )
1667 {
1668 for ( int joffset = 0; joffset < nSc; ++joffset)
1669 dSourceVertexLon[ioffset][joffset] = rowcolsvals[i+joffset];
1670 i += nSc;
1671 }
1672 }
1673
1674 ioffset = 0;
1675 for ( int ip = 0; ip < nprocs; ++ip )
1676 {
1677 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] + rootSizesData[NDATA * ip + 4] * nSc;
1678 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4] * nSc;
1679 for ( int i = istart; i < iend; ++i, ++ioffset )
1680 {
1681 for ( int joffset = 0; joffset < nSc; ++joffset)
1682 dSourceVertexLat[ioffset][joffset] = rowcolsvals[i+joffset];
1683 i += nSc;
1684 }
1685 }
1686 }
1687
1688 DataVector<double>& dTargetCenterLon = m_weightMapGlobal->GetTargetCenterLon();
1689 DataVector<double>& dTargetCenterLat = m_weightMapGlobal->GetTargetCenterLat();
1690 dTargetCenterLon.Initialize(gtgtdofs);
1691 dTargetCenterLat.Initialize(gtgtdofs);
1692 {
1693 int offset = 0;
1694 for ( int ip = 0; ip < nprocs; ++ip )
1695 {
1696 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc);
1697 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5];
1698 for ( int i = istart; i < iend; ++i, ++offset )
1699 {
1700 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs);
1701 // dTargetCenterLon[tgtelmindx[offset]] = rowcolsvals[i];
1702 dTargetCenterLon[offset] = rowcolsvals[i];
1703 }
1704 }
1705
1706 offset = 0;
1707 for ( int ip = 0; ip < nprocs; ++ip )
1708 {
1709 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + rootSizesData[NDATA * ip + 5];
1710 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5];
1711 for ( int i = istart; i < iend; ++i, ++offset )
1712 {
1713 // dTargetCenterLat[tgtelmindx[offset]] = rowcolsvals[i];
1714 dTargetCenterLat[offset] = rowcolsvals[i];
1715 }
1716 }
1717 }
1718
1719 DataMatrix<double>& dTargetVertexLon = m_weightMapGlobal->GetTargetVertexLon();
1720 DataMatrix<double>& dTargetVertexLat = m_weightMapGlobal->GetTargetVertexLat();
1721 dTargetVertexLon.Initialize(gtgtdofs, nTc);
1722 dTargetVertexLat.Initialize(gtgtdofs, nTc);
1723 {
1724 int ioffset = 0;
1725 int nc = dTargetVertexLon.GetColumns();
1726 for ( int ip = 0; ip < nprocs; ++ip )
1727 {
1728 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + 2 * rootSizesData[NDATA * ip + 5];
1729 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]*nTc;
1730 for ( int i = istart; i < iend; ++i, ++ioffset )
1731 {
1732 for ( int joffset = 0; joffset < nc; ++joffset)
1733 dTargetVertexLon[ioffset][joffset] = rowcolsvals[i+joffset];
1734 i += nc;
1735 }
1736 }
1737
1738 ioffset = 0;
1739 for ( int ip = 0; ip < nprocs; ++ip )
1740 {
1741 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + 2 * rootSizesData[NDATA * ip + 5] + rootSizesData[NDATA * ip + 5]*nTc;
1742 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]*nTc;
1743 for ( int i = istart; i < iend; ++i, ++ioffset )
1744 {
1745 for ( int joffset = 0; joffset < nc; ++joffset)
1746 dTargetVertexLat[ioffset][joffset] = rowcolsvals[i+joffset];
1747 i += nc;
1748 }
1749 }
1750 }
1751 }
1752
1753 }
1754 rootSizesData.clear();
1755
1756 // Update all processes that we have a global view of the map available
1757 // on the root process
1758 m_globalMapAvailable = true;
1759
1760 #ifdef VERBOSE
1761 if ( !rank )
1762 {
1763 dbgprint.printf ( 0, "Writing out file outGlobalView.nc\n" );
1764 m_weightMapGlobal->Write ( "outGlobalView.nc" );
1765 }
1766 #endif
1767 return moab::MB_SUCCESS;
1768 }
1769
1770 ///////////////////////////////////////////////////////////////////////////////
1771 #endif
1772
1773