Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
linear_advection.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <sstream>
3 #include <ctime>
4 #include <cstdlib>
5 #include <string>
6 #include <cmath>
7 #include "moab/Core.hpp"
8 #include "moab/Interface.hpp"
9 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
10 #include "moab/ParallelComm.hpp"
11 #include "moab/ProgOptions.hpp"
12 #include "MBParallelConventions.h"
13 #include "moab/ReadUtilIface.hpp"
14 #include "MBTagConventions.hpp"
15 #include "IntxUtilsCSLAM.hpp"
16 
17 #include "TestUtil.hpp"
18 
19 // std::string file_name("./uniform_30.g");
20 // std::string file_name("./uniform_120.g");
21 // std::string file_name("./eulerHomme.vtk");
22 
23 using namespace moab;
24 
25 moab::ErrorCode update_density( moab::Interface* mb,
26  moab::EntityHandle euler_set,
27  moab::EntityHandle lagr_set,
28  moab::EntityHandle out_set,
29  moab::Tag& rhoTag,
30  moab::Tag& areaTag,
31  moab::Tag& rhoCoefsTag,
32  moab::Tag& weightsTag,
33  moab::Tag& planeTag );
34 
35 moab::ErrorCode get_departure_grid( moab::Interface* mb,
36  moab::EntityHandle euler_set,
37  moab::EntityHandle lagr_set,
38  moab::EntityHandle covering_set,
39  int tStep,
40  moab::Range& connecVerts );
41 
42 moab::ErrorCode get_barycenters( moab::Interface* mb,
43  moab::EntityHandle set,
44  moab::Tag& planeTag,
45  moab::Tag& barycenterTag,
46  moab::Tag& areaTag );
47 moab::ErrorCode get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag );
48 moab::ErrorCode get_linear_reconstruction( moab::Interface* mb,
49  moab::EntityHandle set,
50  moab::Tag& rhoTag,
51  moab::Tag& planeTag,
52  moab::Tag& barycenterTag,
53  moab::Tag& linearCoefTag );
54 moab::ErrorCode get_intersection_weights( moab::Interface* mb,
55  moab::EntityHandle euler_set,
56  moab::EntityHandle lagr_set,
57  moab::EntityHandle intx_set,
58  moab::Tag& planeTag,
59  moab::Tag& weightsTag );
60 
61 moab::ErrorCode set_density( moab::Interface* mb,
62  moab::EntityHandle euler_set,
63  moab::Tag& barycenterTag,
64  moab::Tag& rhoTag,
65  int field_type );
66 
67 moab::ErrorCode create_lagr_mesh( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set );
68 
69 // functions to compute departure point locations
70 void departure_point_swirl( moab::CartVect& arrival_point, double t, double delta_t, moab::CartVect& departure_point );
71 void departure_point_swirl_rot( moab::CartVect& arrival_point,
72  double t,
73  double delta_t,
74  moab::CartVect& departure_point );
75 
76 double gtol = 1.e-9; // this is for geometry tolerance
77 
78 double radius = 1.;
79 
80 bool writeFiles = true;
81 bool parallelWrite = false;
82 bool velocity = false;
83 
84 int numSteps = 200; // number of times with velocity displayed at points
85 double T = 5;
86 
87 Intx2MeshOnSphere* pworker = NULL;
88 
89 int main( int argc, char* argv[] )
90 {
91 
92  MPI_Init( &argc, &argv );
93 
94  // set up MOAB interface and parallel communication
95  moab::Core moab;
96  moab::Interface& mb = moab;
97  moab::ParallelComm mb_pcomm( &mb, MPI_COMM_WORLD );
98 
99  // int rank = mb_pcomm->proc_config().proc_rank();
100  int rank = mb_pcomm.proc_config().proc_rank();
101 
102  // create meshset
103  moab::EntityHandle euler_set;
104  moab::ErrorCode rval = mb.create_meshset( MESHSET_SET, euler_set );MB_CHK_ERR( rval );
105 
106  // std::stringstream opts;
107  // opts <<
108  // "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;GATHER_SET=0;PARTITION_METHOD=TRIVIAL_PARTITION;VARIABLE=";
109  // opts << "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS";
110  // rval = mb.load_file(file_name.c_str(), &euler_set, opts.str().c_str());
111  std::string fileN = TestDir + "unittest/mbcslam/fine4.h5m";
112 
113  rval = mb.load_file( fileN.c_str(), &euler_set );MB_CHK_ERR( rval );
114 
115  // Create tag for cell density
116  moab::Tag rhoTag = 0;
117  rval = mb.tag_get_handle( "Density", 1, moab::MB_TYPE_DOUBLE, rhoTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
118 
119  // Create tag for cell area
120  moab::Tag areaTag = 0;
121  rval = mb.tag_get_handle( "Area", 1, moab::MB_TYPE_DOUBLE, areaTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
122  // Create tag for cell barycenters in 3D Cartesian space
123  moab::Tag barycenterTag = 0;
124  rval = mb.tag_get_handle( "CellBarycenter", 3, moab::MB_TYPE_DOUBLE, barycenterTag,
125  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
126  // Create tag for cell density reconstruction coefficients
127  moab::Tag rhoCoefTag = 0;
128  rval = mb.tag_get_handle( "LinearCoefRho", 3, moab::MB_TYPE_DOUBLE, rhoCoefTag,
129  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
130  // Create tag for index of gnomonic plane for each cell
131  moab::Tag planeTag = 0;
132  rval = mb.tag_get_handle( "gnomonicPlane", 1, moab::MB_TYPE_INTEGER, planeTag,
133  moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
134  // Create tag for intersection weights
135  moab::Tag weightsTag = 0;
136  rval = mb.tag_get_handle( "Weights", 3, moab::MB_TYPE_DOUBLE, weightsTag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE );MB_CHK_ERR( rval );
137 
138  // get cell plane
139  rval = get_gnomonic_plane( &mb, euler_set, planeTag );MB_CHK_ERR( rval );
140 
141  // get cell barycenters (and cell area)
142  rval = get_barycenters( &mb, euler_set, planeTag, areaTag, barycenterTag );MB_CHK_ERR( rval );
143 
144  // Set density distributions
145  rval = set_density( &mb, euler_set, barycenterTag, rhoTag, 1 );MB_CHK_ERR( rval );
146 
147  // Get initial values for use in error computation
148  moab::Range redEls;
149  rval = mb.get_entities_by_dimension( euler_set, 2, redEls );MB_CHK_ERR( rval );
150  std::vector< double > iniValsRho( redEls.size() );
151  rval = mb.tag_get_data( rhoTag, redEls, &iniValsRho[0] );MB_CHK_ERR( rval );
152 
153  // Get Lagrangian set
154  moab::EntityHandle out_set, lagrange_set, covering_set;
155  rval = mb.create_meshset( MESHSET_SET, out_set );MB_CHK_ERR( rval );
156  rval = mb.create_meshset( MESHSET_SET, lagrange_set );MB_CHK_ERR( rval );
157  rval = mb.create_meshset( MESHSET_SET, covering_set );MB_CHK_ERR( rval );
158  // rval = create_lagr_mesh(&mb, euler_set, lagrange_set);
159  rval = IntxUtilsCSLAM::deep_copy_set( &mb, euler_set, lagrange_set );MB_CHK_ERR( rval );
160  moab::EntityHandle dum = 0;
161  moab::Tag corrTag;
162  rval = mb.tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
163  // Set up intersection of two meshes
164 
165  /*
166  moab::Intx2MeshOnSphere worker(&mb);
167  worker.SetRadius(radius);
168  worker.SetErrorTolerance(gtol);
169  */
170 
171  pworker = new Intx2MeshOnSphere( &mb );
172  pworker->set_error_tolerance( gtol );
173  pworker->set_radius_source_mesh( radius );
174  pworker->set_radius_destination_mesh( radius );
175  pworker->set_box_error( 100 * gtol );
176 
177  rval = pworker->FindMaxEdges( euler_set, euler_set );MB_CHK_ERR( rval );
178 
179  // these stay fixed for one run
180  moab::Range local_verts;
181  rval = pworker->build_processor_euler_boxes( euler_set,
182  local_verts ); // output also the local_verts
183  MB_CHK_ERR( rval );
184  // rval = worker.build_processor_euler_boxes(euler_set, local_verts);// output also the
185  // local_verts
186 
187  // loop over time to update density
188  for( int ts = 1; ts < numSteps + 1; ts++ )
189  {
190 
191  if( ts == 1 ) // output initial condition
192  {
193  std::stringstream newDensity;
194  newDensity << "Density" << rank << "_" << ts - 1 << ".vtk";
195  rval = mb.write_file( newDensity.str().c_str(), 0, 0, &euler_set, 1 );
196  }
197 
198  // get linear reconstruction coefficients
199  get_linear_reconstruction( &mb, euler_set, rhoTag, planeTag, barycenterTag, rhoCoefTag );
200 
201  // get depature grid
202  rval = get_departure_grid( &mb, euler_set, lagrange_set, covering_set, ts, local_verts );MB_CHK_ERR( rval );
203 
204  // intersect the meshes
205  rval = pworker->intersect_meshes( lagrange_set, euler_set, out_set );MB_CHK_ERR( rval );
206 
207  // intersection weights (i.e. area, x integral, and y integral over cell intersections)
208  get_intersection_weights( &mb, euler_set, lagrange_set, out_set, planeTag, weightsTag );
209 
210  // update the density
211  rval =
212  update_density( &mb, euler_set, lagrange_set, out_set, rhoTag, areaTag, rhoCoefTag, weightsTag, planeTag );MB_CHK_ERR( rval );
213 
214  if( writeFiles && ( ts % 5 == 0 ) ) // so if write
215  {
216  std::stringstream newDensity;
217  newDensity << "Density" << rank << "_" << ts << ".vtk";
218  rval = mb.write_file( newDensity.str().c_str(), 0, 0, &euler_set, 1 );MB_CHK_ERR( rval );
219  }
220 
221  // delete the polygons and elements of out_set
222  moab::Range allVerts;
223  rval = mb.get_entities_by_dimension( 0, 0, allVerts );MB_CHK_ERR( rval );
224 
225  moab::Range allElems;
226  rval = mb.get_entities_by_dimension( 0, 2, allElems );MB_CHK_ERR( rval );
227 
228  // get Eulerian and lagrangian cells
229  moab::Range polys;
230  rval = mb.get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
231  rval =
232  mb.get_entities_by_dimension( lagrange_set, 2, polys ); // do not delete lagr set either, with its vertices
233  MB_CHK_ERR( rval );
234 
235  // add to the connecVerts range all verts, from all initial polys
236  moab::Range vertsToStay;
237  rval = mb.get_connectivity( polys, vertsToStay );MB_CHK_ERR( rval );
238 
239  moab::Range todeleteVerts = subtract( allVerts, vertsToStay );
240 
241  moab::Range todeleteElem = subtract( allElems, polys );
242  // empty the out mesh set
243  rval = mb.clear_meshset( &out_set, 1 );MB_CHK_ERR( rval );
244 
245  rval = mb.delete_entities( todeleteElem );MB_CHK_ERR( rval );
246  rval = mb.delete_entities( todeleteVerts );MB_CHK_ERR( rval );
247  if( rank == 0 ) std::cout << " step: " << ts << "\n";
248  }
249 
250  // final vals and errors
251  moab::Range::iterator iter = redEls.begin();
252  double norm1 = 0.;
253  double norm2 = 0.;
254  double exact2 = 0.;
255  double exact1 = 0.;
256  int count = 0;
257  void* data;
258  int j = 0; // index in iniVals
259  while( iter != redEls.end() )
260  {
261  rval = mb.tag_iterate( rhoTag, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
262  double* ptrTracer = (double*)data;
263 
264  rval = mb.tag_iterate( areaTag, iter, redEls.end(), count, data );MB_CHK_ERR( rval );
265  double* ptrArea = (double*)data;
266  for( int i = 0; i < count; i++, ++iter, ptrTracer++, ptrArea++, j++ )
267  {
268  // double area = *ptrArea;
269  norm1 += fabs( *ptrTracer - iniValsRho[j] ) * ( *ptrArea );
270  norm2 += ( *ptrTracer - iniValsRho[j] ) * ( *ptrTracer - iniValsRho[j] ) * ( *ptrArea );
271  exact1 += ( iniValsRho[j] ) * ( *ptrArea );
272  exact2 += ( iniValsRho[j] ) * ( iniValsRho[j] ) * ( *ptrArea );
273  }
274  }
275 
276  double total_norm1 = 0;
277  double total_norm2 = 0;
278  double total_exact1 = 0;
279  double total_exact2 = 0;
280  int mpi_err = MPI_Reduce( &norm1, &total_norm1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
281  if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
282  mpi_err = MPI_Reduce( &norm2, &total_norm2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
283  if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
284  mpi_err = MPI_Reduce( &exact1, &total_exact1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
285  if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
286  mpi_err = MPI_Reduce( &exact2, &total_exact2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
287  if( mpi_err ) std::cout << " error in MPI_reduce:" << mpi_err << "\n";
288  if( 0 == rank )
289  std::cout << " numSteps:" << numSteps << " 1-norm:" << total_norm1 / total_exact1
290  << " 2-norm:" << total_norm2 / total_exact2 << "\n";
291 
292  MPI_Finalize();
293  return 0;
294 }
295 
296 moab::ErrorCode get_gnomonic_plane( moab::Interface* mb, moab::EntityHandle set, moab::Tag& planeTag )
297 {
298  // get all entities of dimension 2
299  moab::Range cells;
300  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
301 
302  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
303  {
304  moab::EntityHandle icell = *it;
305 
306  // get the nodes
307  const moab::EntityHandle* verts;
308  int num_nodes;
309  rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
310 
311  // get coordinates
312  std::vector< double > coords( 3 * num_nodes );
313  rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
314 
315  // get cell center
316  double centerx = 0;
317  double centery = 0;
318  double centerz = 0;
319  for( int inode = 0; inode < num_nodes; inode++ )
320  {
321  centerx += coords[inode * 3] / num_nodes;
322  centery += coords[inode * 3 + 1] / num_nodes;
323  centerz += coords[inode * 3 + 2] / num_nodes;
324  }
325  double rad = std::sqrt( centerx * centerx + centery * centery + centerz * centerz );
326  centerx = centerx / rad;
327  centery = centery / rad;
328  centerz = centerz / rad;
329  moab::CartVect center( centerx, centery, centerz );
330 
331  // define gnomonic plane based on cell center coordinates
332  int plane = 0;
333  moab::IntxUtils::decide_gnomonic_plane( center, plane );
334 
335  rval = mb->tag_set_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
336  }
337 
338  return moab::MB_SUCCESS;
339 }
340 
341 moab::ErrorCode get_barycenters( moab::Interface* mb,
342  moab::EntityHandle set,
343  moab::Tag& planeTag,
344  moab::Tag& areaTag,
345  moab::Tag& barycenterTag )
346 {
347 
348  // get all entities of dimension 2
349  moab::Range cells;
350  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
351 
352  // set sphere radius to 1
353  double R = 1.0;
354 
355  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
356  {
357  moab::EntityHandle icell = *it;
358 
359  // get the nodes
360  const moab::EntityHandle* verts;
361  int num_nodes;
362  rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
363 
364  // get coordinates
365  std::vector< double > coords( 3 * num_nodes );
366  rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
367 
368  // get gnomonic plane
369  int plane = 0;
370  rval = mb->tag_get_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
371 
372  // get vertex coordinates and project onto gnomonic plane
373  std::vector< double > x( num_nodes );
374  std::vector< double > y( num_nodes );
375  double area = 0;
376  double bary_x = 0;
377  double bary_y = 0;
378  for( int inode = 0; inode < num_nodes; inode++ )
379  {
380  double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
381  coords[inode * 3 + 2] * coords[inode * 3 + 2] );
382  CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad, coords[inode * 3 + 2] / rad );
383  moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane, x[inode], y[inode] );
384  }
385 
386  // integrate over cell to get barycenter in gnomonic coordinates
387  for( int inode = 0; inode < num_nodes; inode++ )
388  {
389  int inode2 = inode + 1;
390  if( inode2 >= num_nodes ) inode2 = 0;
391  double xmid = 0.5 * ( x[inode] + x[inode2] );
392  double ymid = 0.5 * ( y[inode] + y[inode2] );
393  double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
394  double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
395  double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
396  double hx = x[inode2] - x[inode];
397 
398  area += -hx *
399  ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
400  y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
401  6.0;
402 
403  bary_x += -hx *
404  ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
405  4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
406  x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
407  6.0;
408 
409  bary_y += hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
410  }
411  bary_x = bary_x / area;
412  bary_y = bary_y / area;
413 
414  // barycenter in Cartesian X,Y,Z coordinates
415  moab::CartVect barycent;
416  moab::IntxUtils::reverse_gnomonic_projection( bary_x, bary_y, R, plane, barycent );
417 
418  // set barycenter
419  std::vector< double > barycenter( 3 );
420  barycenter[0] = barycent[0];
421  barycenter[1] = barycent[1];
422  barycenter[2] = barycent[2];
423  rval = mb->tag_set_data( barycenterTag, &icell, 1, &barycenter[0] );MB_CHK_ERR( rval );
424 
425  // set area
426  rval = mb->tag_set_data( areaTag, &icell, 1, &area );MB_CHK_ERR( rval );
427  }
428 
429  return moab::MB_SUCCESS;
430 }
431 
432 moab::ErrorCode set_density( moab::Interface* mb,
433  moab::EntityHandle euler_set,
434  moab::Tag& barycenterTag,
435  moab::Tag& rhoTag,
436  int field_type )
437 {
438  // get cells
439  moab::Range cells;
440  moab::ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, cells );MB_CHK_ERR( rval );
441 
442  // get barycenters
443  std::vector< double > cell_barys( 3 * cells.size() );
444  rval = mb->tag_get_data( barycenterTag, cells, &cell_barys[0] );MB_CHK_ERR( rval );
445 
446  // loop over cells
447  int cell_ind = 0;
448  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
449  {
450  moab::EntityHandle icell = *it;
451 
452  // convert barycenter from 3-D Cartesian to lat/lon
453  moab::CartVect bary_xyz( cell_barys[cell_ind * 3], cell_barys[cell_ind * 3 + 1], cell_barys[cell_ind * 3 + 2] );
454  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( bary_xyz );
455 
456  if( field_type == 1 ) // cosine bells
457  {
458  // lon1, lat1 lon2 lat2 b c hmax r
459  double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6, 0.0, 0.1, 0.9, 1., 0.5 };
460 
461  double rho_barycent = IntxUtilsCSLAM::quasi_smooth_field( sphCoord.lon, sphCoord.lat, params );
462  rval = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
463  }
464 
465  if( field_type == 2 ) // Gaussian hills
466  {
467  moab::CartVect p1, p2;
468  moab::IntxUtils::SphereCoords spr;
469  spr.R = 1;
470  spr.lat = M_PI / 3;
471  spr.lon = M_PI;
472  p1 = moab::IntxUtils::spherical_to_cart( spr );
473  spr.lat = -M_PI / 3;
474  p2 = moab::IntxUtils::spherical_to_cart( spr );
475  // X1, Y1, Z1, X2, Y2, Z2, ?, ?
476  double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5. };
477 
478  double rho_barycent = IntxUtilsCSLAM::smooth_field( sphCoord.lon, sphCoord.lat, params );
479  rval = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
480  }
481 
482  if( field_type == 3 ) // Zalesak cylinders
483  {
484  // lon1, lat1, lon2, lat2, b, c, r
485  double params[] = { 5 * M_PI / 6.0, 0.0, 7 * M_PI / 6.0, 0.0, 0.1, 0.9, 0.5 };
486 
487  double rho_barycent = IntxUtilsCSLAM::slotted_cylinder_field( sphCoord.lon, sphCoord.lat, params );
488  rval = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
489  }
490 
491  if( field_type == 4 ) // constant
492  {
493  double rho_barycent = 1.0;
494  rval = mb->tag_set_data( rhoTag, &icell, 1, &rho_barycent );MB_CHK_ERR( rval );
495  }
496 
497  cell_ind++;
498  }
499 
500  return moab::MB_SUCCESS;
501 }
502 
503 moab::ErrorCode get_linear_reconstruction( moab::Interface* mb,
504  moab::EntityHandle set,
505  moab::Tag& rhoTag,
506  moab::Tag& planeTag,
507  moab::Tag& barycenterTag,
508  moab::Tag& linearCoefTag )
509 {
510  // get all entities of dimension 2
511  Range cells;
512  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells );MB_CHK_ERR( rval );
513 
514  // Get coefficients for reconstruction (in cubed-sphere coordinates)
515  for( Range::iterator it = cells.begin(); it != cells.end(); ++it )
516  {
517  moab::EntityHandle icell = *it;
518 
519  // get the nodes, then the coordinates
520  const moab::EntityHandle* verts;
521  int num_nodes;
522  rval = mb->get_connectivity( icell, verts, num_nodes );MB_CHK_ERR( rval );
523 
524  moab::Range adjacentEdges;
525  rval = mb->get_adjacencies( &icell, 1, 1, true, adjacentEdges );MB_CHK_ERR( rval );
526 
527  // get adjacent cells from edges
528  moab::Range adjacentCells;
529  rval = mb->get_adjacencies( adjacentEdges, 2, true, adjacentCells, Interface::UNION );MB_CHK_ERR( rval );
530 
531  // get gnomonic plane
532  int plane = 0;
533  rval = mb->tag_get_data( planeTag, &icell, 1, &plane );MB_CHK_ERR( rval );
534 
535  std::vector< double > dx( adjacentCells.size() - 1 );
536  std::vector< double > dy( adjacentCells.size() - 1 );
537  std::vector< double > dr( adjacentCells.size() - 1 );
538  double bary_x;
539  double bary_y;
540 
541  // get barycenter of cell where reconstruction occurs
542  double rad = 1;
543  std::vector< double > barycent( 3 );
544  rval = mb->tag_get_data( barycenterTag, &icell, 1, &barycent[0] );MB_CHK_ERR( rval );
545  CartVect barycenter( barycent[0], barycent[1], barycent[2] );
546  double cellbaryx = 0;
547  double cellbaryy = 0;
548  moab::IntxUtils::gnomonic_projection( barycenter, rad, plane, cellbaryx, cellbaryy );
549 
550  // get density value
551  double rhocell = 0;
552  rval = mb->tag_get_data( rhoTag, &icell, 1, &rhocell );MB_CHK_ERR( rval );
553 
554  // get barycenters of surrounding cells
555  std::vector< double > cell_barys( 3 * adjacentCells.size() );
556  rval = mb->tag_get_data( barycenterTag, adjacentCells, &cell_barys[0] );MB_CHK_ERR( rval );
557 
558  // get density of surrounding cells
559  std::vector< double > rho_vals( adjacentCells.size() );
560  rval = mb->tag_get_data( rhoTag, adjacentCells, &rho_vals[0] );MB_CHK_ERR( rval );
561 
562  std::size_t jind = 0;
563  for( std::size_t i = 0; i < adjacentCells.size(); i++ )
564  {
565 
566  if( adjacentCells[i] != icell )
567  {
568 
569  CartVect bary_xyz( cell_barys[i * 3], cell_barys[i * 3 + 1], cell_barys[i * 3 + 2] );
570  moab::IntxUtils::gnomonic_projection( bary_xyz, rad, plane, bary_x, bary_y );
571 
572  dx[jind] = bary_x - cellbaryx;
573  dy[jind] = bary_y - cellbaryy;
574  dr[jind] = rho_vals[i] - rhocell;
575 
576  jind++;
577  }
578  }
579 
580  std::vector< double > linearCoef( 3 );
581  if( adjacentCells.size() == 5 )
582  {
583 
584  // compute normal equations matrix
585  double N11 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2] + dx[3] * dx[3];
586  double N22 = dy[0] * dy[0] + dy[1] * dy[1] + dy[2] * dy[2] + dy[3] * dy[3];
587  double N12 = dx[0] * dy[0] + dx[1] * dy[1] + dx[2] * dy[2] + dx[3] * dy[3];
588 
589  // rhs
590  double Rx = dx[0] * dr[0] + dx[1] * dr[1] + dx[2] * dr[2] + dx[3] * dr[3];
591  double Ry = dy[0] * dr[0] + dy[1] * dr[1] + dy[2] * dr[2] + dy[3] * dr[3];
592 
593  // determinant
594  double Det = N11 * N22 - N12 * N12;
595 
596  // solution
597  linearCoef[0] = ( Rx * N22 - Ry * N12 ) / Det;
598  linearCoef[1] = ( Ry * N11 - Rx * N12 ) / Det;
599  linearCoef[2] = rhocell - linearCoef[0] * cellbaryx - linearCoef[1] * cellbaryy;
600  }
601  else
602  {
603  // default to first order
604  linearCoef[0] = 0.0;
605  linearCoef[1] = 0.0;
606  linearCoef[2] = rhocell;
607  std::cout << "Need 4 adjacent cells for linear reconstruction! \n";
608  }
609 
610  rval = mb->tag_set_data( linearCoefTag, &icell, 1, &linearCoef[0] );MB_CHK_ERR( rval );
611  }
612 
613  return moab::MB_SUCCESS;
614 }
615 
616 moab::ErrorCode get_departure_grid( moab::Interface* mb,
617  moab::EntityHandle euler_set,
618  moab::EntityHandle lagr_set,
619  moab::EntityHandle covering_set,
620  int tStep,
621  Range& connecVerts )
622 {
623  EntityHandle dum = 0;
624  Tag corrTag;
625  mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );
626 
627  double t = tStep * T / numSteps; // numSteps is global; so is T
628  double delta_t = T / numSteps; // this is global too, actually
629  // double delta_t = 0.0001;
630  // double t = delta_t;
631 
632  Range polys;
633  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
634 
635  // change coordinates of lagr mesh vertices
636  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
637  {
638  moab::EntityHandle oldV = *vit;
639  CartVect posi;
640  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
641 
642  // Intx utils, case 1
643  CartVect newPos;
644  departure_point_swirl_rot( posi, t, delta_t, newPos );
645  newPos = radius * newPos; // do we need this? the radius should be 1
646  moab::EntityHandle new_vert;
647  rval = mb->tag_get_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
648 
649  // set the new position for the new vertex
650  rval = mb->set_coords( &new_vert, 1, &( newPos[0] ) );MB_CHK_ERR( rval );
651  }
652 
653  // if in parallel, we have to move some elements to another proc, and receive other cells
654  // from other procs
655  rval = pworker->create_departure_mesh_3rd_alg( lagr_set, covering_set );MB_CHK_ERR( rval );
656 
657  return rval;
658 }
659 
660 // !!! For now serial !!!
661 moab::ErrorCode update_density( moab::Interface* mb,
662  moab::EntityHandle euler_set,
663  moab::EntityHandle lagr_set,
664  moab::EntityHandle out_set,
665  moab::Tag& rhoTag,
666  moab::Tag& areaTag,
667  moab::Tag& rhoCoefsTag,
668  moab::Tag& weightsTag,
669  moab::Tag& planeTag )
670 {
671  // moab::ParallelComm * parcomm = ParallelComm::get_pcomm(mb, 0);
672  ErrorCode rval;
673  double R = 1.0;
674  moab::EntityHandle dum = 0;
675  Tag corrTag;
676 
677  rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dum );MB_CHK_ERR( rval );
678 
679  // get all polygons out of out_set; then see where are they coming from
680  moab::Range polys;
681  rval = mb->get_entities_by_dimension( out_set, 2, polys );MB_CHK_ERR( rval );
682 
683  // get all Lagrangian cells
684  moab::Range rs1;
685  rval = mb->get_entities_by_dimension( lagr_set, 2, rs1 );MB_CHK_ERR( rval );
686 
687  // get all Eulerian cells
688  moab::Range rs2;
689  rval = mb->get_entities_by_dimension( euler_set, 2, rs2 );MB_CHK_ERR( rval );
690 
691  // get gnomonic plane for Eulerian cells
692  std::vector< int > plane( rs2.size() );
693  rval = mb->tag_get_data( planeTag, rs2, &plane[0] );MB_CHK_ERR( rval );
694 
695  // get Eulerian cell reconstruction coefficients
696  std::vector< double > rhoCoefs( 3 * rs2.size() );
697  rval = mb->tag_get_data( rhoCoefsTag, rs2, &rhoCoefs[0] );MB_CHK_ERR( rval );
698 
699  // get intersection weights
700  std::vector< double > weights( 3 * polys.size() );
701  rval = mb->tag_get_data( weightsTag, polys, &weights[0] );MB_CHK_ERR( rval );
702 
703  // Initialize the new values
704  std::vector< double > newValues( rs2.size(), 0. ); // initialize with 0 all of them
705 
706  // For each polygon get red/blue parent
707  moab::Tag tgtParentTag;
708  moab::Tag srcParentTag;
709  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
710  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
711 
712  // mass_lagr = (\sum_intx \int rho^h(x,y) dV)
713  // rho_eul^n+1 = mass_lagr/area_eul
714  double check_intx_area = 0.;
715  int polyIndex = 0;
716  for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
717  {
718  moab::EntityHandle poly = *it;
719  int blueIndex, redIndex;
720  rval = mb->tag_get_data( srcParentTag, &poly, 1, &blueIndex );MB_CHK_ERR( rval );
721 
722  moab::EntityHandle blue = rs1[blueIndex];
723  rval = mb->tag_get_data( tgtParentTag, &poly, 1, &redIndex );MB_CHK_ERR( rval );
724 
725  moab::EntityHandle redArr;
726  rval = mb->tag_get_data( corrTag, &blue, 1, &redArr );MB_CHK_ERR( rval );
727  int arrRedIndex = rs2.index( redArr );
728 
729  // sum into new density values
730  newValues[arrRedIndex] += rhoCoefs[redIndex * 3] * weights[polyIndex * 3] +
731  rhoCoefs[redIndex * 3 + 1] * weights[polyIndex * 3 + 1] +
732  rhoCoefs[redIndex * 3 + 2] * weights[polyIndex * 3 + 2];
733 
734  check_intx_area += weights[polyIndex * 3 + 2];
735 
736  polyIndex++;
737  }
738 
739  // now divide by red area (current)
740  int j = 0;
741  Range::iterator iter = rs2.begin();
742  void* data = NULL; // used for stored area
743  int count = 0;
744  double total_mass_local = 0.;
745  while( iter != rs2.end() )
746  {
747  rval = mb->tag_iterate( areaTag, iter, rs2.end(), count, data );MB_CHK_ERR( rval );
748 
749  double* ptrArea = (double*)data;
750  for( int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
751  {
752  total_mass_local += newValues[j];
753  newValues[j] /= ( *ptrArea );
754  }
755  }
756 
757  rval = mb->tag_set_data( rhoTag, rs2, &newValues[0] );MB_CHK_ERR( rval );
758 
759  std::cout << "total mass now:" << total_mass_local << "\n";
760  std::cout << "check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI * R * R << " " << check_intx_area
761  << "\n";
762 
763  return MB_SUCCESS;
764 }
765 
766 /*
767  * Deformational flow
768  */
769 void departure_point_swirl( moab::CartVect& arrival_point, double t, double delta_t, moab::CartVect& departure_point )
770 {
771 
772  // always assume radius is 1 here?
773  moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
774  double k = 2.4; // flow parameter
775  /* radius needs to be within some range */
776  double sl2 = sin( sph.lon / 2 );
777  double pit = M_PI * t / T;
778  double omega = M_PI / T;
779  double costheta = cos( sph.lat );
780  // double u = k * sl2*sl2 * sin(2*sph.lat) * cos(pit);
781  double v = k * sin( sph.lon ) * costheta * cos( pit );
782  // double psi = k * sl2 * sl2 *costheta * costheta * cos(pit);
783  double u_tilda = 2 * k * sl2 * sl2 * sin( sph.lat ) * cos( pit );
784 
785  // formula 35, page 8
786  // this will approximate dep point using a Taylor series with up to second derivative
787  // this will be O(delta_t^3) exact.
788  double lon_dep =
789  sph.lon - delta_t * u_tilda -
790  delta_t * delta_t * k * sl2 *
791  ( sl2 * sin( sph.lat ) * sin( pit ) * omega - u_tilda * sin( sph.lat ) * cos( pit ) * cos( sph.lon / 2 ) -
792  v * sl2 * costheta * cos( pit ) );
793  // formula 36, page 8 again
794  double lat_dep = sph.lat - delta_t * v -
795  delta_t * delta_t / 4 * k *
796  ( sin( sph.lon ) * cos( sph.lat ) * sin( pit ) * omega -
797  u_tilda * cos( sph.lon ) * cos( sph.lat ) * cos( pit ) +
798  v * sin( sph.lon ) * sin( sph.lat ) * cos( pit ) );
799  moab::IntxUtils::SphereCoords sph_dep;
800  sph_dep.R = 1.; // radius
801  sph_dep.lat = lat_dep;
802  sph_dep.lon = lon_dep;
803 
804  departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
805  return;
806 }
807 
808 /*
809  * Deformational flow with rotation
810  */
811 void departure_point_swirl_rot( moab::CartVect& arrival_point,
812  double t,
813  double delta_t,
814  moab::CartVect& departure_point )
815 {
816 
817  moab::IntxUtils::SphereCoords sph = moab::IntxUtils::cart_to_spherical( arrival_point );
818  double omega = M_PI / T;
819  double gt = cos( M_PI * t / T );
820 
821  double lambda = sph.lon - 2.0 * omega * t;
822  double u_tilda = 4.0 * sin( lambda ) * sin( lambda ) * sin( sph.lat ) * gt + 2.0 * omega;
823  double v = 2.0 * sin( 2.0 * lambda ) * cos( sph.lat ) * gt;
824 
825  double lon_dep = sph.lon - delta_t * u_tilda -
826  delta_t * delta_t * 2.0 * sin( lambda ) *
827  ( sin( lambda ) * sin( sph.lat ) * sin( omega * t ) * omega -
828  sin( lambda ) * cos( sph.lat ) * cos( omega * t ) * v -
829  2.0 * cos( lambda ) * sin( sph.lat ) * cos( omega * t ) * u_tilda );
830 
831  double lat_dep = sph.lat - delta_t * v -
832  delta_t * delta_t * 2.0 *
833  ( cos( sph.lat ) * sin( omega * t ) * omega * sin( lambda ) * cos( lambda ) -
834  2.0 * u_tilda * cos( sph.lat ) * cos( omega * t ) * cos( lambda ) * cos( lambda ) +
835  u_tilda * cos( sph.lat ) * cos( omega * t ) +
836  v * sin( sph.lat ) * cos( omega * t ) * sin( lambda ) * cos( lambda ) );
837 
838  moab::IntxUtils::SphereCoords sph_dep;
839  sph_dep.R = 1.; // radius
840  sph_dep.lat = lat_dep;
841  sph_dep.lon = lon_dep;
842 
843  departure_point = moab::IntxUtils::spherical_to_cart( sph_dep );
844  return;
845 }
846 
847 /*
848  * Zonal flow
849  */
850 moab::ErrorCode get_intersection_weights( moab::Interface* mb,
851  moab::EntityHandle euler_set,
852  moab::EntityHandle lagr_set,
853  moab::EntityHandle intx_set,
854  moab::Tag& planeTag,
855  moab::Tag& weightsTag )
856 {
857  // get all intersection polygons
858  moab::Range polys;
859  ErrorCode rval = mb->get_entities_by_dimension( intx_set, 2, polys );MB_CHK_ERR( rval );
860 
861  // get all Eulerian cells
862  moab::Range eul_cells;
863  rval = mb->get_entities_by_dimension( euler_set, 2, eul_cells );MB_CHK_ERR( rval );
864 
865  // get all Lagrangian cells
866  moab::Range lagr_cells;
867  rval = mb->get_entities_by_dimension( lagr_set, 2, lagr_cells );MB_CHK_ERR( rval );
868 
869  // get tag for Eulerian parent cell of intersection polygon
870  moab::Tag tgtParentTag;
871  rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE );
872 
873  // get tag for Lagrangian parent cell of intersection polygon
874  moab::Tag srcParentTag;
875  rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
876 
877  // get gnomonic plane for Eulerian cells
878  std::vector< int > plane( eul_cells.size() );
879  rval = mb->tag_get_data( planeTag, eul_cells, &plane[0] );MB_CHK_ERR( rval );
880 
881  double total_area = 0.;
882  for( moab::Range::iterator it = polys.begin(); it != polys.end(); ++it )
883  {
884  moab::EntityHandle poly = *it;
885 
886  // get the nodes
887  const moab::EntityHandle* verts;
888  int num_nodes;
889  rval = mb->get_connectivity( poly, verts, num_nodes );MB_CHK_ERR( rval );
890 
891  // get coordinates
892  std::vector< double > coords( 3 * num_nodes );
893  rval = mb->get_coords( verts, num_nodes, &coords[0] );MB_CHK_ERR( rval );
894 
895  // get index of Eulerian parent cell for polygon
896  int redIndex;
897  rval = mb->tag_get_data( tgtParentTag, &poly, 1, &redIndex );MB_CHK_ERR( rval );
898 
899  // get index of Lagrangian parent cell for polygon
900  int blueIndex;
901  rval = mb->tag_get_data( srcParentTag, &poly, 1, &blueIndex );MB_CHK_ERR( rval );
902 
903  std::vector< double > x( num_nodes );
904  std::vector< double > y( num_nodes );
905  double poly_area = 0;
906  double poly_intx = 0;
907  double poly_inty = 0;
908  double R = 1.0;
909  for( int inode = 0; inode < num_nodes; inode++ )
910  {
911  double rad = sqrt( coords[inode * 3] * coords[inode * 3] + coords[inode * 3 + 1] * coords[inode * 3 + 1] +
912  coords[inode * 3 + 2] * coords[inode * 3 + 2] );
913  moab::CartVect xyzcoord( coords[inode * 3] / rad, coords[inode * 3 + 1] / rad,
914  coords[inode * 3 + 2] / rad );
915  moab::IntxUtils::gnomonic_projection( xyzcoord, R, plane[redIndex], x[inode], y[inode] );
916  }
917 
918  std::vector< double > weights( 3 );
919  for( int inode = 0; inode < num_nodes; inode++ )
920  {
921  int inode2 = inode + 1;
922  if( inode2 >= num_nodes ) inode2 = 0;
923  double xmid = 0.5 * ( x[inode] + x[inode2] );
924  double ymid = 0.5 * ( y[inode] + y[inode2] );
925  double r1 = sqrt( 1 + x[inode] * x[inode] + y[inode] * y[inode] );
926  double rm = sqrt( 1 + xmid * xmid + ymid * ymid );
927  double r2 = sqrt( 1 + x[inode2] * x[inode2] + y[inode2] * y[inode2] );
928  double hx = x[inode2] - x[inode];
929 
930  poly_area += -hx *
931  ( y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) + 4.0 * ymid / ( rm * ( 1 + xmid * xmid ) ) +
932  y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
933  6.0;
934 
935  poly_intx += -hx *
936  ( x[inode] * y[inode] / ( r1 * ( 1 + x[inode] * x[inode] ) ) +
937  4.0 * xmid * ymid / ( rm * ( 1 + xmid * xmid ) ) +
938  x[inode2] * y[inode2] / ( r2 * ( 1 + x[inode2] * x[inode2] ) ) ) /
939  6.0;
940 
941  poly_inty += hx * ( 1.0 / r1 + 4.0 / rm + 1.0 / r2 ) / 6.0;
942  }
943  weights[0] = poly_intx;
944  weights[1] = poly_inty;
945  weights[2] = poly_area;
946 
947  total_area += poly_area;
948 
949  rval = mb->tag_set_data( weightsTag, &poly, 1, &weights[0] );MB_CHK_ERR( rval );
950  }
951 
952  std::cout << "polygon area = " << total_area << "\n";
953 
954  return moab::MB_SUCCESS;
955 }
956 
957 ErrorCode create_lagr_mesh( moab::Interface* mb, moab::EntityHandle euler_set, moab::EntityHandle lagr_set )
958 {
959  // create the handle tag for the corresponding element / vertex
960 
961  moab::EntityHandle dum = 0;
962 
963  moab::Tag corrTag;
964  ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
965 
966  MB_CHK_ERR( rval );
967  // give the same global id to new verts and cells created in the lagr(departure) mesh
968  moab::Tag gid = mb->globalId_tag();
969 
970  moab::Range polys;
971  rval = mb->get_entities_by_dimension( euler_set, 2, polys );MB_CHK_ERR( rval );
972 
973  moab::Range connecVerts;
974  rval = mb->get_connectivity( polys, connecVerts );MB_CHK_ERR( rval );
975 
976  std::map< moab::EntityHandle, moab::EntityHandle > newNodes;
977  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
978  {
979  moab::EntityHandle oldV = *vit;
980  moab::CartVect posi;
981  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
982  int global_id;
983  rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
984  moab::EntityHandle new_vert;
985  // duplicate the position
986  rval = mb->create_vertex( &( posi[0] ), new_vert );MB_CHK_ERR( rval );
987  newNodes[oldV] = new_vert;
988  // set also the correspondent tag :)
989  rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
990 
991  // also the other side
992  // need to check if we really need this; the new vertex will never need the old vertex
993  // we have the global id which is the same
994  /*rval = mb->tag_set_data(corrTag, &new_vert, 1, &oldV);MB_CHK_ERR(rval);*/
995  // set the global id on the corresponding vertex the same as the initial vertex
996  rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
997  }
998 
999  for( Range::iterator it = polys.begin(); it != polys.end(); ++it )
1000  {
1001  moab::EntityHandle q = *it;
1002  int global_id;
1003  int nnodes;
1004  const moab::EntityHandle* conn;
1005 
1006  rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
1007  rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
1008 
1009  EntityType typeElem = mb->type_from_handle( q );
1010  std::vector< moab::EntityHandle > new_conn( nnodes );
1011  for( int i = 0; i < nnodes; i++ )
1012  {
1013  moab::EntityHandle v1 = conn[i];
1014  new_conn[i] = newNodes[v1];
1015  }
1016 
1017  moab::EntityHandle newElement;
1018  rval = mb->create_element( typeElem, &new_conn[0], nnodes, newElement );MB_CHK_ERR( rval );
1019  // set the corresponding tag; not sure we need this one, from old to new
1020  /*rval = mb->tag_set_data(corrTag, &q, 1, &newElement);MB_CHK_ERR(rval);*/
1021  rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
1022  // set the global id
1023  rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
1024 
1025  rval = mb->add_entities( lagr_set, &newElement, 1 );MB_CHK_ERR( rval );
1026  }
1027 
1028  return MB_SUCCESS;
1029 }