Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
linear_remap.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <cmath>
3 
4 #include "moab/Core.hpp"
5 #include "moab/Interface.hpp"
6 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
7 #include "moab/ParallelComm.hpp"
8 #include "moab/ProgOptions.hpp"
9 #include "MBParallelConventions.h"
10 #include "moab/ReadUtilIface.hpp"
11 #include "MBTagConventions.hpp"
12 
13 #include "moab/IntxMesh/IntxUtils.hpp"
14 #include "IntxUtilsCSLAM.hpp"
15 
16 #include "TestUtil.hpp"
17 
18 // std::string file_name("./uniform_15.g");
19 // std::string file_name("./eulerHomme.vtk");
20 
21 using namespace moab;
22 
23 // computes cell barycenters in Cartesian X,Y,Z coordinates
24 void get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag, moab::Tag& barycenterTag );
25 
26 // computes gnomonic plane for each cell of the Eulerian mesh
27 void get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag );
28 
29 // computes coefficients (A,B,C) for linear reconstruction in gnomonic coordinates: rho(x,y) = Ax +
30 // By + C
31 void get_linear_reconstruction( moab::Interface* mb,
32  moab::EntityHandle set,
33  moab::Tag& rhoTag,
34  moab::Tag& planeTag,
35  moab::Tag& barycenterTag,
36  moab::Tag& linearCoefTag );
37 
38 // evaluates the integral of rho(x,y) over cell, should be equal to cell average rho
39 void test_linear_reconstruction( moab::Interface* mb,
40  moab::EntityHandle set,
41  moab::Tag& rhoTag,
42  moab::Tag& planeTag,
43  moab::Tag& barycenterTag,
44  moab::Tag& linearCoefTag );
45 
46 // set density function
47 moab::ErrorCode add_field_value( moab::Interface* mb,
48  moab::EntityHandle euler_set,
49  int rank,
50  moab::Tag& tagTracer,
51  moab::Tag& tagElem,
52  moab::Tag& tagArea,
53  int field_type );
54 
55 // functions that implement gnomonic projection as in Homme (different from Moab implementation)
56 int gnomonic_projection_test( const moab::CartVect& pos, double R, int plane, double& c1, double& c2 );
57 int reverse_gnomonic_projection_test( const double& c1, const double& c2, double R, int plane, moab::CartVect& pos );
58 void decide_gnomonic_plane_test( const CartVect& pos, int& plane );
59 
60 double radius = 1.; // in m: 6371220.
61 
62 int main( int argc, char* argv[] )
63 {
64 
65  // set up MOAB interface and parallel communication
66  MPI_Init( &argc, &argv );
67  moab::Core moab;
68  moab::Interface& mb = moab;
69  moab::ParallelComm mb_pcomm( &mb, MPI_COMM_WORLD );
70 
71  // int rank = mb_pcomm->proc_config().proc_rank();
72  int rank = mb_pcomm.proc_config().proc_rank();
73 
74  // create meshset
75  moab::EntityHandle euler_set;
76  moab::ErrorCode rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
77 
78  // std::stringstream opts;
79  // opts <<
80  // "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;GATHER_SET=0;PARTITION_METHOD=TRIVIAL_PARTITION;VARIABLE=";
81  // opts << "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
82  std::string fileN = TestDir + "unittest/mbcslam/fine4.h5m";
83 
84  rval = mb.load_file( fileN.c_str(), &euler_set );CHECK_ERR( rval );
85 
86  // Create tag for cell density
87  moab::Tag rhoTag = 0;
88  rval = mb.tag_get_handle( "Density", 1, moab::MB_TYPE_DOUBLE, rhoTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
89  moab::Tag rhoNodeTag = 0;
90  rval = mb.tag_get_handle( "DensityNode", 1, moab::MB_TYPE_DOUBLE, rhoNodeTag,
91  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
92  // Create tag for cell area
93  moab::Tag areaTag = 0;
94  rval = mb.tag_get_handle( "Area", 1, moab::MB_TYPE_DOUBLE, areaTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
95  // Create tag for cell barycenters in 3D Cartesian space
96  moab::Tag barycenterTag = 0;
97  rval = mb.tag_get_handle( "CellBarycenter", 3, moab::MB_TYPE_DOUBLE, barycenterTag,
98  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
99  // Create tag for cell density reconstruction coefficients
100  moab::Tag coefRhoTag = 0;
101  rval = mb.tag_get_handle( "LinearCoefRho", 3, moab::MB_TYPE_DOUBLE, coefRhoTag,
102  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
103  // Create tag for index of gnomonic plane for each cell
104  moab::Tag planeTag = 0;
105  rval = mb.tag_get_handle( "gnomonicPlane", 1, moab::MB_TYPE_INTEGER, planeTag,
106  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );CHECK_ERR( rval );
107  // Set density distributions
108  rval = add_field_value( &mb, euler_set, rank, rhoNodeTag, rhoTag, areaTag, 2 );CHECK_ERR( rval );
109  // get cell plane
110  get_gnomonic_plane( &mb, euler_set, planeTag );
111 
112  // get cell barycenters
113  get_barycenters( &mb, euler_set, planeTag, barycenterTag );
114 
115  // get linear reconstruction
116  get_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
117 
118  // test linear reconstruction
119  test_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, coefRhoTag );
120 
121  MPI_Finalize();
122  return 0;
123 }
124 
125 void get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag )
126 {
127  // get all entities of dimension 2
128  // then get the connectivity, etc
129  moab::Range cells;
130  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
131  if( MB_SUCCESS != rval ) return;
132 
133  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
134  {
135  moab::EntityHandle icell = *it;
136  // get the nodes, then the coordinates
137  const moab::EntityHandle* verts;
138  int num_nodes;
139  rval = mb->get_connectivity( icell, verts, num_nodes );
140  if( MB_SUCCESS != rval ) return;
141 
142  std::vector< double > coords( 3 * num_nodes );
143  // get coordinates
144  rval = mb->get_coords( verts, num_nodes, &coords[0] );
145  if( MB_SUCCESS != rval ) return;
146 
147  double centerx = 0;
148  double centery = 0;
149  double centerz = 0;
150  for( int inode = 0; inode < num_nodes; inode++ )
151  {
152  centerx += coords[inode * 3] / num_nodes;
153  centery += coords[inode * 3 + 1] / num_nodes;
154  centerz += coords[inode * 3 + 2] / num_nodes;
155  }
156  double R = std::sqrt( centerx * centerx + centery * centery + centerz * centerz );
157  centerx = centerx / R;
158  centery = centery / R;
159  centerz = centerz / R;
160  moab::CartVect center( centerx, centery, centerz );
161  R = 1.0;
162 
163  int plane = 0;
164  moab::IntxUtils::decide_gnomonic_plane( center, plane );
165  // decide_gnomonic_plane_test(center,plane);
166 
167  rval = mb->tag_set_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
168  }
169  return;
170 }
171 
172 void get_barycenters( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag, moab::Tag& barycenterTag )
173 {
174  // get all entities of dimension 2
175  moab::Range cells;
176  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
177  if( MB_SUCCESS != rval ) return;
178 
179  // set sphere radius to 1
180  double R = 1.0;
181 
182  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
183  {
184  moab::EntityHandle icell = *it;
185 
186  // get the nodes
187  const moab::EntityHandle* verts;
188  int num_nodes;
189  rval = mb->get_connectivity( icell, verts, num_nodes );
190  if( MB_SUCCESS != rval ) return;
191 
192  // get coordinates
193  std::vector< double > coords( 3 * num_nodes );
194  rval = mb->get_coords( verts, num_nodes, &coords[0] );
195  if( MB_SUCCESS != rval ) return;
196 
197  // get plane for cell
198  int plane = 0;
199  rval = mb->tag_get_data( planeTag, &icell, 1, &plane );
200  if( MB_SUCCESS != rval ) return;
201  std::vector< double > x( num_nodes );
202  std::vector< double > y( num_nodes );
203  double area = 0;
204  double bary_x = 0;
205  double bary_y = 0;
206  for( int inode = 0; inode < num_nodes; inode++ )
207  {
208  // radius should be 1.0, but divide by it just in case for now
209  double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
210  coords[inode * 3 + 2] * coords[inode * 3 + 2] );
211  CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
212  moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
213  // int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
214  }
215 
216  for( int inode = 0; inode < num_nodes; inode++ )
217  {
218  int inode2 = inode + 1;
219  if( inode2 >= num_nodes ) inode2 = 0;
220  double xmid = 0.5 * ( x[inode] + x[inode2] );
221  double ymid = 0.5 * ( y[inode] + y[inode2] );
222  double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
223  double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
224  double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
225  double hx = x[inode2] - x[inode];
226 
227  area += hx *
228  ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
229  y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
230  6.0;
231 
232  bary_x += hx *
233  ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
234  4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
235  x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
236  6.0;
237 
238  bary_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
239  }
240 
241  bary_x = bary_x / area;
242  bary_y = bary_y / area;
243 
244  moab::CartVect barycent;
245  moab::IntxUtils::reverse_gnomonic_projection( bary_x, bary_y, R, plane, barycent );
246  // reverse_gnomonic_projection_test(bary_x, bary_y, R, plane, barycent);
247  std::vector< double > barycenter( 3 );
248  barycenter[0] = barycent[0];
249  barycenter[1] = barycent[1];
250  barycenter[2] = barycent[2];
251 
252  rval = mb->tag_set_data( barycenterTag, &icell, 1, &barycenter[0] );CHECK_ERR( rval );
253  }
254  return;
255 }
256 
257 void get_linear_reconstruction( moab::Interface* mb,
258  moab::EntityHandle set,
259  moab::Tag& rhoTag,
260  moab::Tag& planeTag,
261  moab::Tag& barycenterTag,
262  moab::Tag& linearCoefTag )
263 {
264  // get all entities of dimension 2
265  Range cells;
266  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
267  if( MB_SUCCESS != rval ) return;
268 
269  // set sphere radius to 1
270  double R = 1;
271 
272  // Get coefficients of reconstruction (in cubed-sphere coordinates)
273  // rho(x,y) = Ax + By + C
274  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
275  {
276  moab::EntityHandle icell = *it;
277 
278  // get the nodes, then the coordinates
279  const moab::EntityHandle* verts;
280  int num_nodes;
281  rval = mb->get_connectivity( icell, verts, num_nodes );
282  if( MB_SUCCESS != rval ) return;
283 
284  moab::Range adjacentEdges;
285  rval = mb->get_adjacencies( &icell, 1, 1, true, adjacentEdges );CHECK_ERR( rval );
286  // int num_edges = adjacentEdges.size();
287 
288  // get adjacent cells from edges
289  moab::Range adjacentCells;
290  rval = mb->get_adjacencies( adjacentEdges, 2, true, adjacentCells, Interface::UNION );
291  if( MB_SUCCESS != rval ) return;
292  // get plane for cell
293  int plane = 0;
294  rval = mb->tag_get_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
295 
296  std::vector< double > dx( adjacentCells.size() - 1 );
297  std::vector< double > dy( adjacentCells.size() - 1 );
298  std::vector< double > dr( adjacentCells.size() - 1 );
299  double bary_x;
300  double bary_y;
301 
302  // get barycenter of cell where reconstruction occurs
303  std::vector< double > barycent( 3 );
304  rval = mb->tag_get_data( barycenterTag, &icell, 1, &barycent[0] );CHECK_ERR( rval );
305  CartVect barycenter( barycent[0], barycent[1], barycent[2] );
306  double cellbaryx = 0;
307  double cellbaryy = 0;
308  moab::IntxUtils::gnomonic_projection( barycenter, R, plane, cellbaryx, cellbaryy );
309  // int rc = gnomonic_projection_test(barycenter, R, plane, cellbaryx, cellbaryy);
310 
311  // get density value
312  double rhocell = 0;
313  rval = mb->tag_get_data( rhoTag, &icell, 1, &rhocell );CHECK_ERR( rval );
314 
315  // get barycenters of surrounding cells
316  std::vector< double > cell_barys( 3 * adjacentCells.size() );
317  rval = mb->tag_get_data( barycenterTag, adjacentCells, &cell_barys[0] );CHECK_ERR( rval );
318 
319  // get density of surrounding cells
320  std::vector< double > rho_vals( adjacentCells.size() );
321  rval = mb->tag_get_data( rhoTag, adjacentCells, &rho_vals[0] );CHECK_ERR( rval );
322 
323  std::size_t jind = 0;
324  for( std::size_t i = 0; i < adjacentCells.size(); i++ )
325  {
326 
327  if( adjacentCells[i] != icell )
328  {
329  CartVect bary_xyz( cell_barys[i * 3], cell_barys[i * 3 + 1], cell_barys[i * 3 + 2] );
330  moab::IntxUtils::gnomonic_projection( bary_xyz, R, plane, bary_x, bary_y );
331  // rc = gnomonic_projection_test(bary_xyz, R, plane, bary_x, bary_y);
332 
333  dx[jind] = bary_x - cellbaryx;
334  dy[jind] = bary_y - cellbaryy;
335  dr[jind] = rho_vals[i] - rhocell;
336 
337  jind++;
338  }
339  }
340 
341  std::vector< double > linearCoef( 3 );
342  if( adjacentCells.size() == 5 )
343  {
344 
345  // compute normal equations matrix
346  double N11 = dx[1] * dx[1] + dx[2] * dx[2] + dx[3] * dx[3] + dx[0] * dx[0];
347  double N22 = dy[1] * dy[1] + dy[2] * dy[2] + dy[3] * dy[3] + dy[0] * dy[0];
348  double N12 = dx[1] * dy[1] + dx[2] * dy[2] + dx[3] * dy[3] + dx[0] * dy[0];
349 
350  // rhs
351  double Rx = dx[1] * dr[1] + dx[2] * dr[2] + dx[3] * dr[3] + dx[0] * dr[0];
352  double Ry = dy[1] * dr[1] + dy[2] * dr[2] + dy[3] * dr[3] + dy[0] * dr[0];
353 
354  // determinant
355  double Det = N11 * N22 - N12 * N12;
356 
357  // solution
358  linearCoef[0] = ( Rx * N22 - Ry * N12 ) / Det;
359  linearCoef[1] = ( Ry * N11 - Rx * N12 ) / Det;
360  linearCoef[2] = rhocell - linearCoef[0] * cellbaryx - linearCoef[1] * cellbaryy;
361  }
362  else
363  {
364  // default to first order
365  linearCoef[0] = 0.0;
366  linearCoef[1] = 0.0;
367  linearCoef[2] = rhocell;
368  std::cout << "Need 4 adjacent cells for linear reconstruction! \n";
369  }
370 
371  rval = mb->tag_set_data( linearCoefTag, &icell, 1, &linearCoef[0] );CHECK_ERR( rval );
372  }
373  return;
374 }
375 
376 void test_linear_reconstruction( moab::Interface* mb,
377  moab::EntityHandle set,
378  moab::Tag& rhoTag,
379  moab::Tag& planeTag,
380  moab::Tag& barycenterTag,
381  moab::Tag& linearCoefTag )
382 {
383  // get all entities of dimension 2
384  Range cells;
385  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );
386  if( MB_SUCCESS != rval ) return;
387 
388  // set sphere radius to 1
389  double R = 1;
390 
391  // For get coefficients for reconstruction (in cubed-sphere coordinates)
392  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
393  {
394  moab::EntityHandle icell = *it;
395 
396  // get the nodes, then the coordinates
397  const moab::EntityHandle* verts;
398  int num_nodes;
399  rval = mb->get_connectivity( icell, verts, num_nodes );
400  if( MB_SUCCESS != rval ) return;
401 
402  // get coordinates
403  std::vector< double > coords( 3 * num_nodes );
404  rval = mb->get_coords( verts, num_nodes, &coords[0] );
405  if( MB_SUCCESS != rval ) return;
406 
407  // get plane for cell
408  int plane = 0;
409  rval = mb->tag_get_data( planeTag, &icell, 1, &plane );CHECK_ERR( rval );
410 
411  // get vertex coordinates projections
412  std::vector< double > x( num_nodes );
413  std::vector< double > y( num_nodes );
414  for( int inode = 0; inode < num_nodes; inode++ )
415  {
416  double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
417  coords[inode * 3 + 2] * coords[inode * 3 + 2] );
418  CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
419  moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
420  // int dum = gnomonic_projection_test(xyzcoord, R, plane, x[inode],y[inode]);
421  }
422 
423  double area = 0;
424  double int_x = 0;
425  double int_y = 0;
426  for( int inode = 0; inode < num_nodes; inode++ )
427  {
428  int inode2 = inode + 1;
429  if( inode2 >= num_nodes ) inode2 = 0;
430  double xmid = 0.5 * ( x[inode] + x[inode2] );
431  double ymid = 0.5 * ( y[inode] + y[inode2] );
432  double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
433  double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
434  double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
435  double hx = x[inode2] - x[inode];
436 
437  area += hx *
438  ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
439  y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
440  6.0;
441 
442  int_x += hx *
443  ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
444  4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
445  x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
446  6.0;
447 
448  int_y += -hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
449  }
450 
451  // get linear coeficients
452  std::vector< double > rho_coefs( 3 );
453  rval = mb->tag_get_data( linearCoefTag, &icell, 1, &rho_coefs[0] );CHECK_ERR( rval );
454 
455  // get barycenters
456  std::vector< double > bary( 3 );
457  rval = mb->tag_get_data( barycenterTag, &icell, 1, &bary[0] );CHECK_ERR( rval );
458  double bary_x;
459  double bary_y;
460  CartVect bary_xyz( bary[0], bary[1], bary[2] );
461  moab::IntxUtils::gnomonic_projection( bary_xyz, R, plane, bary_x, bary_y );
462  // int rc = gnomonic_projection_test(bary_xyz, R, plane, bary_x, bary_y);
463 
464  // get cell average density
465  double cell_rho;
466  rval = mb->tag_get_data( rhoTag, &icell, 1, &cell_rho );CHECK_ERR( rval );
467 
468  // ave rho = \int rho^h(x,y) dV / area = (\int (Ax + By + C) dV) / area
469  double rho_test1 = ( rho_coefs[0] * int_x + rho_coefs[1] * int_y + rho_coefs[2] * area ) / area;
470 
471  // ave rho = A*bary_x + B*bary_y + C
472  double rho_test2 = rho_coefs[0] * bary_x + rho_coefs[1] * bary_y + rho_coefs[2];
473 
474  std::cout << cell_rho << " " << rho_test1 << " " << rho_test2 << " " << cell_rho - rho_test1 << "\n";
475  }
476  return;
477 }
478 
479 int reverse_gnomonic_projection_test( const double& c1, const double& c2, double R, int plane, CartVect& pos )
480 {
481 
482  double x = c1;
483  double y = c2;
484  double r = sqrt( 1.0 + x * x + y * y );
485 
486  switch( plane )
487  {
488 
489  case 1: {
490  pos[0] = R / r * R;
491  pos[1] = R / r * x;
492  pos[2] = R / r * y;
493  break;
494  }
495  case 2: {
496  pos[0] = -R / r * x;
497  pos[1] = R / r * R;
498  pos[2] = R / r * y;
499  break;
500  }
501  case 3: {
502  pos[0] = -R / r * R;
503  pos[1] = -R / r * x;
504  pos[2] = R / r * y;
505  break;
506  }
507  case 4: {
508  pos[0] = R / r * x;
509  pos[1] = -R / r * R;
510  pos[2] = R / r * y;
511  break;
512  }
513  case 5: {
514  pos[0] = R / r * y;
515  pos[1] = R / r * x;
516  pos[2] = -R / r * R;
517  break;
518  }
519  case 6: {
520  pos[0] = -R / r * y;
521  pos[1] = R / r * x;
522  pos[2] = R / r * R;
523  break;
524  }
525  }
526 
527  return 0; // no error
528 }
529 
530 void decide_gnomonic_plane_test( const CartVect& pos, int& plane )
531 {
532 
533  // This is from early version of Homme vorticity calculation in parvis
534  // Poles are reversed from Homme and Iulian version.
535 
536  // Now has been changed for consistency
537 
538  double X = pos[0];
539  double Y = pos[1];
540  double Z = pos[2];
541  double R = sqrt( X * X + Y * Y + Z * Z );
542  X = X / R;
543  Y = Y / R;
544  Z = Z / R;
545 
546  if( ( Y < X ) & ( Y > -X ) )
547  {
548  if( Z > X )
549  {
550  plane = 6;
551  }
552  else if( Z < -X )
553  {
554  plane = 5;
555  }
556  else
557  {
558  plane = 1;
559  }
560  }
561  else if( ( Y > X ) & ( Y < -X ) )
562  {
563  if( Z > -X )
564  {
565  plane = 6;
566  }
567  else if( Z < X )
568  {
569  plane = 5;
570  }
571  else
572  {
573  plane = 3;
574  }
575  }
576  else if( ( Y > X ) & ( Y > -X ) )
577  {
578  if( Z > Y )
579  {
580  plane = 6;
581  }
582  else if( Z < -Y )
583  {
584  plane = 5;
585  }
586  else
587  {
588  plane = 2;
589  }
590  }
591  else if( ( Y < X ) & ( Y < -X ) )
592  {
593  if( Z > -Y )
594  {
595  plane = 6;
596  }
597  else if( Z < Y )
598  {
599  plane = 5;
600  }
601  else
602  {
603  plane = 4;
604  }
605  }
606  else
607  {
608  if( fabs( X ) < Z )
609  {
610  plane = 6;
611  }
612  else if( Z < -fabs( X ) )
613  {
614  plane = 5;
615  }
616  else if( ( X > 0 ) & ( Y > 0 ) )
617  {
618  plane = 1;
619  }
620  else if( ( X < 0 ) & ( Y > 0 ) )
621  {
622  plane = 2;
623  }
624  else if( ( X < 0 ) & ( Y < 0 ) )
625  {
626  plane = 3;
627  }
628  else
629  {
630  plane = 4;
631  }
632  }
633 
634  return;
635 }
636 
637 moab::ErrorCode add_field_value( moab::Interface* mb,
638  moab::EntityHandle euler_set,
639  int /*rank*/,
640  moab::Tag& tagTracer,
641  moab::Tag& tagElem,
642  moab::Tag& tagArea,
643  int field_type )
644 {
645 
646  /*
647  * get all plys first, then vertices, then move them on the surface of the sphere
648  * radius is 1., most of the time
649  *
650  */
651  moab::Range polygons;
652  moab::ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polygons );
653  if( MB_SUCCESS != rval ) return rval;
654 
655  moab::Range connecVerts;
656  rval = mb->get_connectivity( polygons, connecVerts );
657  if( MB_SUCCESS != rval ) return rval;
658 
659  void* data; // pointer to the LOC in memory, for each vertex
660  int count;
661 
662  rval = mb->tag_iterate( tagTracer, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
663  // here we are checking contiguity
664  assert( count == (int)connecVerts.size() );
665  double* ptr_DP = (double*)data;
666  // lambda is for longitude, theta for latitude
667  // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2
668  // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3)
669  // (la2, te2) = (M_PI, -M_PI/3)
670  // la1, te1 la2 te2 b c hmax r
671  if( field_type == 1 ) // quasi smooth
672  {
673  double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6, 0.0, 0.1, 0.9, 1., 0.5 };
674  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
675  {
676  moab::EntityHandle oldV = *vit;
677  moab::CartVect posi;
678  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
679 
680  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
681 
682  ptr_DP[0] = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
683 
684  ptr_DP++; // increment to the next node
685  }
686  }
687  else if( 2 == field_type ) // smooth
688  {
689  moab::CartVect p1, p2;
690  moab::IntxUtils::SphereCoords spr;
691  spr.R = 1;
692  spr.lat = M_PI / 3;
693  spr.lon = M_PI;
694  p1 = moab::IntxUtils::spherical_to_cart( spr );
695  spr.lat = -M_PI / 3;
696  p2 = moab::IntxUtils::spherical_to_cart( spr );
697  // x1, y1, z1, x2, y2, z2, h_max, b0
698  double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
699  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
700  {
701  moab::EntityHandle oldV = *vit;
702  moab::CartVect posi;
703  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
704 
705  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
706 
707  ptr_DP[0] = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
708 
709  ptr_DP++; // increment to the next node
710  }
711  }
712  else if( 3 == field_type ) // slotted
713  {
714  // la1, te1, la2, te2, b, c, r
715  double params[] = { M_PI, M_PI / 3, M_PI, -M_PI / 3, 0.1, 0.9, 0.5 }; // no h_max
716  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
717  {
718  moab::EntityHandle oldV = *vit;
719  moab::CartVect posi;
720  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
721 
722  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
723 
724  ptr_DP[0] = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
725  ;
726 
727  ptr_DP++; // increment to the next node
728  }
729  }
730  else if( 4 == field_type ) // constant = 1
731  {
732  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
733  {
734  /* moab::EntityHandle oldV = *vit;
735  moab::CartVect posi;
736  rval = mb->get_coords(&oldV, 1, &(posi[0]));
737 
738  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical(posi);*/
739 
740  ptr_DP[0] = 1.0;
741 
742  ptr_DP++; // increment to the next node
743  }
744  }
745  // add average value for quad/polygon (average corners)
746  // do some averages
747 
748  Range::iterator iter = polygons.begin();
749  while( iter != polygons.end() )
750  {
751  rval = mb->tag_iterate( tagElem, iter, polygons.end(), count, data );CHECK_ERR( rval );
752  double* ptr = (double*)data;
753 
754  rval = mb->tag_iterate( tagArea, iter, polygons.end(), count, data );CHECK_ERR( rval );
755  double* ptrArea = (double*)data;
756  for( int i = 0; i < count; i++, ++iter, ptr++, ptrArea++ )
757  {
758  const moab::EntityHandle* conn = NULL;
759  int num_nodes = 0;
760  rval = mb->get_connectivity( *iter, conn, num_nodes );CHECK_ERR( rval );
761  if( num_nodes == 0 ) return MB_FAILURE;
762  std::vector< double > nodeVals( num_nodes );
763  double average = 0.;
764  rval = mb->tag_get_data( tagTracer, conn, num_nodes, &nodeVals[0] );CHECK_ERR( rval );
765  for( int j = 0; j < num_nodes; j++ )
766  average += nodeVals[j];
767  average /= num_nodes;
768  *ptr = average;
769  }
770  }
771 
772  return MB_SUCCESS;
773 }