SpatialOps
FieldHelper.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2014 The University of Utah
3  *
4  * Permission is hereby granted, free of charge, to any person obtaining a copy
5  * of this software and associated documentation files (the "Software"), to
6  * deal in the Software without restriction, including without limitation the
7  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8  * sell copies of the Software, and to permit persons to whom the Software is
9  * furnished to do so, subject to the following conditions:
10  *
11  * The above copyright notice and this permission notice shall be included in
12  * all copies or substantial portions of the Software.
13  *
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20  * IN THE SOFTWARE.
21  */
22 
23 #ifndef SpatialOps_FieldHelper_h
24 #define SpatialOps_FieldHelper_h
25 
27 
28 #include <ostream>
29 #include <fstream>
30 #include <string>
31 #include <cmath>
32 #include <iomanip>
33 
36 namespace SpatialOps{
37 
71 template<typename Field>
72 inline void internal_initialize_field( typename Field::iterator fi,
73  const MemoryWindow& mw,
74  const double start,
75  const bool print,
76  const double range )
77 {
78  const int xExtent = mw.extent(0);
79  const int yExtent = mw.extent(1);
80  const int zExtent = mw.extent(2);
81 
82  for( int z = 1; z <= zExtent; ++z ){
83  for( int y = 1; y <= yExtent; ++y ){
84  for( int x = 1; x <= xExtent; ++x, ++fi ){
85  *fi = range * std::sin(start + x + y * xExtent + z * xExtent * yExtent);
86  if( print ) std::cout << *fi << " ";
87  }
88  if( print ) std::cout << std::endl;
89  }
90  if( print ) std::cout << std::endl;
91  }
92 }
93 
123 template<typename Field>
124 inline void initialize_field( Field & f,
125  const double start = 0.0,
126  const bool print = false,
127  const double range = 1.0 )
128 {
129  internal_initialize_field<Field>( f.begin(), f.window_with_ghost(), start, print, range );
130 }
131 
161 template<typename Field>
162 inline void interior_initialize_field( Field & f,
163  const double start = 0.0,
164  const bool print = false,
165  const double range = 1.0 )
166 {
167  internal_initialize_field<Field>( f.interior_begin(), f.window_without_ghost(), start, print, range );
168 }
169 
214 template<typename Field>
215 inline void internal_print_field( typename Field::const_iterator fi,
216  const MemoryWindow& mw,
217  std::ostream& os,
218  bool addFormat )
219 {
220  const int xExtent = mw.extent(0);
221  const int yExtent = mw.extent(1);
222  const int zExtent = mw.extent(2);
223 
224  for( int z = 1; z <= zExtent; ++z ){
225  for( int y = 1; y <= yExtent; ++y ){
226  for( int x = 1; x <= xExtent; ++x, ++fi ){
227  if( !addFormat ) os << *fi << " ";
228  else os << std::setprecision(2) << *fi << "\t";
229  }
230  os << std::endl;
231  }
232  os << std::endl;
233  }
234 }
235 
280 template<typename Field>
281 inline void print_field( const Field& f, std::ostream& os, const bool addFormat = false )
282 {
283  if( IS_GPU_INDEX(const_cast<Field&>(f).active_device_index()) ) const_cast<Field&>(f).add_device(CPU_INDEX);
284  internal_print_field<Field>(f.begin(), f.window_with_ghost(), os, addFormat );
285 };
286 
332 template<typename Field>
333 inline void interior_print_field( const Field& f, std::ostream& os, const bool addFormat = false )
334 {
335  internal_print_field<Field>( f.interior_begin(), f.window_without_ghost(), os, addFormat );
336 };
337 
402 template<typename Field>
403 inline bool internal_display_fields_compare( typename Field::const_iterator fi1,
404  typename Field::const_iterator fi2,
405  const MemoryWindow& mw,
406  const bool display,
407  const bool print )
408 {
409  bool result = true;
410  int xExtent = mw.extent(0);
411  int yExtent = mw.extent(1);
412  int zExtent = mw.extent(2);
413 
414  //copies of iterators for printing
415  typename Field::const_iterator cfi1 = fi1;
416  typename Field::const_iterator cif2 = fi2;
417 
418  // end condition for each test: index < axisExtent && (result || print || display)
419  // this ends the loops early if and only if the result has been found to be false in some cell
420  // AND print == false
421  // AND display == false
422  for(int z = 0; z < zExtent && (result || print || display); ++z ){
423  for(int y = 0; y < yExtent && (result || print || display); ++y ){
424  for(int x = 0; x < xExtent && (result || print || display); ++x, ++fi1, ++fi2 ){
425  const bool compare = (*fi1 == *fi2);
426  result = result && compare;
427  if( display ) std::cout << compare << " ";
428  }
429  if( print ){
430  std::cout << "\t\t";
431  for( int x = 0; x < xExtent; ++x, ++cfi1 ){
432  std::cout << *cfi1 << " ";
433  }
434  std::cout << "\t\t";
435  for( int x = 0; x < xExtent; ++x, ++cif2 ){
436  std::cout << *cif2 << " ";
437  }
438  }
439  if( print || display ) std::cout << std::endl;
440  }
441  if( print || display ) std::cout << std::endl;
442  }
443 
444  return result;
445 }
446 
511 template<typename Field>
512 inline bool display_fields_compare(const Field& field1,
513  const Field& field2,
514  const bool display = false,
515  const bool print = false)
516 {
517  return internal_display_fields_compare<Field>( field1.begin(),
518  field2.begin(),
519  field1.window_with_ghost(),
520  display,
521  print );
522 }
523 
589 template<typename Field>
590 inline bool interior_display_fields_compare( const Field& field1,
591  const Field& field2,
592  const bool display = false,
593  const bool print = false )
594 {
595  return internal_display_fields_compare<Field>( field1.interior_begin(),
596  field2.interior_begin(),
597  field1.window_without_ghost(),
598  display,
599  print );
600 }
601 
602 
612 template<typename FieldT>
613 void write_matlab( const FieldT& field,
614  const std::string& prefix,
615  const bool includeGhost=false )
616 {
617 # ifdef ENABLE_CUDA
618  // IO only works on CPU. Ensure that we have a field there.
619  const_cast<FieldT&>(field).add_device( CPU_INDEX );
620 # endif
621 
622  const std::string fname = "load_"+prefix+".m";
623 
624  std::ofstream fout( fname.c_str() );
625  fout << "function [f,n] = load_" << prefix << "()" << std::endl;
626  fout << std::scientific;
627 
628  const MemoryWindow& mw = includeGhost ? field.window_with_ghost() : field.window_without_ghost();
629 
630  typename FieldT::const_iterator fi = includeGhost ? field.begin() : field.interior_begin();
631  const typename FieldT::const_iterator fie = includeGhost ? field.end() : field.interior_end();
632 
633  const size_t nx = mw.extent(0);
634  const size_t ny = mw.extent(1);
635  const size_t nz = mw.extent(2);
636 
637  fout << "f = [\n";
638 
639 // for( ; fi!=fie; ++fi ){
640 // fout << *fi << "\n";
641 // }
642 // fout << "];\n";
643 // fout << "f = reshape(f," << mw.extent(0) << "," << mw.extent(1) << "," << mw.extent(2) << ");\n";
644 
645  if( nx>1 && ny>1 && nz>1 ){ // 3D
646  for( size_t k=1; k<=nz; ++k ){
647  if( k>1 ) fout << "f(:,:," << k << ") = [ \n";
648  for( size_t j=1; j<=ny; ++j ){
649  for( size_t i=1; i<=nx; ++i, ++fi ){
650  fout << *fi << " ";
651  }
652  fout << std::endl;
653  }
654  fout << "];\n";
655  }
656  fout << "n=[" << nx << "," << ny << "," << nz << "];\n";
657  }
658  else if( mw.extent(0) > 1 && mw.extent(1) > 1 ){ // 2D XY
659  for( size_t j=1; j<=ny; ++j ){
660  for( size_t i=1; i<=nx; ++i, ++fi ){
661  fout << *fi << " ";
662  }
663  fout << std::endl;
664  }
665  fout << "];\nn=[" << nx << "," << ny << "];\n";
666  }
667  else if( mw.extent(0) > 1 && mw.extent(2) > 1 ){ // 2D XZ
668  for( size_t j=1; j<=nz; ++j ){
669  for( size_t i=1; i<=nx; ++i, ++fi ){
670  fout << *fi << " ";
671  }
672  fout << std::endl;
673  }
674  fout << "];\nn=[" << nx << "," << nz << "];\n";
675  }
676  else if( mw.extent(1) > 1 && mw.extent(2) > 1 ){ // 2D YZ
677  for( size_t j=1; j<=nz; ++j ){
678  for( size_t i=1; i<=ny; ++i, ++fi ){
679  fout << *fi << " ";
680  }
681  fout << std::endl;
682  }
683  fout << "];\nn=[" << ny << "," << nz << "];\n";
684  }
685  else if( mw.extent(0) > 1 ){ // 1D X
686  for( ; fi!=fie; ++fi ) fout << *fi << " ";
687  fout << std::endl << "];\n";
688  fout << "n=[" << nx << "];\n";
689  }
690  else if( mw.extent(1) > 1 ){ // 1D Y
691  for( ; fi!=fie; ++fi ) fout << *fi << " ";
692  fout << std::endl << "];\n";
693  fout << "n=[" << ny << "];\n";
694  }
695  else if( mw.extent(2) > 1 ){ // 1D Z
696  for( ; fi!=fie; ++fi ) fout << *fi << " ";
697  fout << std::endl << "];\n";
698  fout << "n=[" << nz << "];\n";
699  }
700 
701  fout.close();
702 }
703 
704 
705 } // namespace SpatialOps
706 
707 #endif
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
void internal_initialize_field(typename Field::iterator fi, const MemoryWindow &mw, const double start, const bool print, const double range)
INTERNAL initialize a field with pseudorandom numbers.
Definition: FieldHelper.h:72
void print_field(const Field &f, std::ostream &os, const bool addFormat=false)
print the values of a field (and ghost cells) to standard output
Definition: FieldHelper.h:281
void initialize_field(Field &f, const double start=0.0, const bool print=false, const double range=1.0)
initialize a field (and ghost cells) with pseudorandom numbers
Definition: FieldHelper.h:124
bool internal_display_fields_compare(typename Field::const_iterator fi1, typename Field::const_iterator fi2, const MemoryWindow &mw, const bool display, const bool print)
INTERNAL compare two fields and possibly print values.
Definition: FieldHelper.h:403
void interior_initialize_field(Field &f, const double start=0.0, const bool print=false, const double range=1.0)
initialize a field (without ghost cells) with pseudorandom numbers
Definition: FieldHelper.h:162
bool display_fields_compare(const Field &field1, const Field &field2, const bool display=false, const bool print=false)
compare two fields and possibly print values (with ghost cells)
Definition: FieldHelper.h:512
void interior_print_field(const Field &f, std::ostream &os, const bool addFormat=false)
print the values of a field (without ghost cells) to standard output
Definition: FieldHelper.h:333
bool interior_display_fields_compare(const Field &field1, const Field &field2, const bool display=false, const bool print=false)
compare two fields and possibly print values (without ghost cells)
Definition: FieldHelper.h:590
void internal_print_field(typename Field::const_iterator fi, const MemoryWindow &mw, std::ostream &os, bool addFormat)
INTERNAL print the values of a field to standard output.
Definition: FieldHelper.h:215