MOAB: Mesh Oriented datABase  (version 5.5.0)
par_intx_sph.cpp
Go to the documentation of this file.
1 /*
2  * par_intx_sph.cpp
3  * test to trigger intersection on a sphere in parallel
4  * it will start from an eulerian mesh
5  * the mesh is read in parallel; lagrangian mesh is manufactured on the fly (part of the test), in
6  * a different set.
7  *
8  * intersections will be performed on the euler mesh
9  * Created on: Nov 14, 2012
10  */
11 
12 #include <iostream>
13 #include <sstream>
14 #include <ctime>
15 #include <cstdlib>
16 #include <cstdio>
17 #include <cstring>
18 #include "moab/Core.hpp"
19 #include "moab/Interface.hpp"
20 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
21 #include <cmath>
22 #include "TestUtil.hpp"
23 #include "moab/ParallelComm.hpp"
24 #include "moab/ProgOptions.hpp"
25 #include "MBParallelConventions.h"
26 #include "moab/ReadUtilIface.hpp"
27 #include "MBTagConventions.hpp"
28 
29 #include "moab/IntxMesh/IntxUtils.hpp"
30 
31 // for M_PI
32 #include <cmath>
33 
34 using namespace moab;
35 // some input data
36 double EPS1 = 0.2; // this is for box error
37 std::string input_mesh_file( "Homme_2pt.h5m" ); // input file, partitioned correctly
38 std::string mpas_file( "mpas_p8.h5m" );
39 double CubeSide = 6.; // the above file starts with cube side 6; radius depends on cube side
40 double radius;
41 void test_intx_in_parallel_elem_based();
42 void test_intx_mpas();
43 
44 int main( int argc, char** argv )
45 {
46  MPI_Init( &argc, &argv );
47  EPS1 = 0.2;
48  int result = 0;
49 
50  if( argc > 1 )
51  {
52  int index = 1;
53  while( index < argc )
54  {
55  if( !strcmp( argv[index], "-eps" ) ) // this is for box error
56  {
57  EPS1 = atof( argv[++index] );
58  }
59  if( !strcmp( argv[index], "-input" ) )
60  {
61  input_mesh_file = argv[++index];
62  }
63  if( !strcmp( argv[index], "-cube" ) )
64  {
65  CubeSide = atof( argv[++index] );
66  }
67  index++;
68  }
69  }
70 
71  radius = CubeSide / 2 * sqrt( 3. );
72  result += RUN_TEST( test_intx_in_parallel_elem_based );
73  radius = 1.;
74  result += RUN_TEST( test_intx_mpas );
75 
76  MPI_Finalize();
77  return result;
78 }
79 // will save the LOC tag on the euler nodes
80 ErrorCode manufacture_lagrange_mesh_on_sphere( Interface* mb, EntityHandle euler_set )
81 {
82  /*
83  * get all quads first, then vertices, then move them on the surface of the sphere
84  * radius is in, it comes from MeshKit/python/examples/manufHomme.py :
85  * length = 6.
86  * each edge of the cube will be divided using this meshcount
87  * meshcount = 11
88  * circumscribed sphere radius
89  * radius = length * math.sqrt(3) /2
90  */
91  // radius = CubeSide/2*sqrt(3.);// our value depends on cube side
92  Range quads;
93  ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, quads );CHECK_ERR( rval );
94 
95  Range connecVerts;
96  rval = mb->get_connectivity( quads, connecVerts );CHECK_ERR( rval );
97 
98  // the LOC tag, should be provided by the user?
99  Tag tagh = 0;
100  std::string tag_name( "DP" );
101  rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );CHECK_ERR( rval );
102  void* data; // pointer to the LOC in memory, for each vertex
103  int count;
104 
105  rval = mb->tag_iterate( tagh, connecVerts.begin(), connecVerts.end(), count, data );CHECK_ERR( rval );
106  // here we are checking contiguity
107  assert( count == (int)connecVerts.size() );
108  double* ptr_DP = (double*)data;
109  // get the coordinates of the old mesh, and move it around the sphere in the same way as in the
110  // python script
111 
112  // now put the vertices in the right place....
113  // int vix=0; // vertex index in new array
114  double t = 0.1, T = 5; // check the script
115  double time = 0.05;
116  double rot = M_PI / 10;
117  for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit )
118  {
119  EntityHandle oldV = *vit;
120  CartVect posi;
121  rval = mb->get_coords( &oldV, 1, &( posi[0] ) );CHECK_ERR( rval );
122  // do some mumbo jumbo, as in python script
123  moab::IntxUtils::SphereCoords sphCoord = moab::IntxUtils::cart_to_spherical( posi );
124  double lat1 = sphCoord.lat - 2 * M_PI * t / T; // 0.1/5
125  double uu = 3 * radius / T * pow( sin( lat1 ), 2 ) * sin( 2 * sphCoord.lon ) * cos( M_PI * t / T );
126  uu += 2 * M_PI * cos( sphCoord.lon ) / T;
127  double vv = 3 * radius / T * ( sin( 2 * lat1 ) ) * cos( sphCoord.lon ) * cos( M_PI * t / T );
128  double vx = -uu * sin( sphCoord.lon ) - vv * sin( sphCoord.lat ) * cos( sphCoord.lon );
129  double vy = -uu * cos( sphCoord.lon ) - vv * sin( sphCoord.lat ) * sin( sphCoord.lon );
130  double vz = vv * cos( sphCoord.lat );
131  posi = posi + time * CartVect( vx, vy, vz );
132  double x2 = posi[0] * cos( rot ) - posi[1] * sin( rot );
133  double y2 = posi[0] * sin( rot ) + posi[1] * cos( rot );
134  CartVect newPos( x2, y2, posi[2] );
135  double len1 = newPos.length();
136  newPos = radius * newPos / len1;
137 
138  ptr_DP[0] = newPos[0];
139  ptr_DP[1] = newPos[1];
140  ptr_DP[2] = newPos[2];
141  ptr_DP += 3; // increment to the next node
142  }
143 
144  return rval;
145 }
146 
147 void test_intx_in_parallel_elem_based()
148 {
149  std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
150  std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
151  Core moab;
152  Interface& mb = moab;
153  EntityHandle euler_set;
154  ErrorCode rval;
155  rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
156  std::string example( TestDir + "unittest/" + input_mesh_file );
157 
158  rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );CHECK_ERR( rval );
159 
160  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
161  CHECK( NULL != pcomm );
162 
163  rval = pcomm->check_all_shared_handles();CHECK_ERR( rval );
164 
165  // everybody will get a DP tag, including the non owned entities; so exchange tags is not
166  // required for LOC (here)
167  rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set );CHECK_ERR( rval );
168 
169  int rank = pcomm->proc_config().proc_rank();
170 
171  std::stringstream ste;
172  ste << "initial" << rank << ".vtk";
173  mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
174 
175  Intx2MeshOnSphere worker( &mb );
176 
177  // double radius= CubeSide/2 * sqrt(3.) ; // input
178  worker.set_radius_source_mesh( radius );
179  worker.set_radius_destination_mesh( radius );
180  worker.set_box_error( EPS1 ); //
181  // worker.SetEntityType(MBQUAD);
182 
183  worker.set_error_tolerance( radius * 1.e-8 );
184  std::cout << "error tolerance epsilon_1=" << radius * 1.e-8 << "\n";
185  // worker.locate_departure_points(euler_set);
186  // set pcomm
187  worker.set_parallel_comm( pcomm );
188 
189  rval = worker.FindMaxEdges( euler_set, euler_set ); // use euler set for both, it is just finding max_edges_1 and 2
190  // we need to make sure the covering set is bigger than the euler mesh
191  EntityHandle covering_lagr_set;
192  rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );CHECK_ERR( rval );
193 
194  rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );CHECK_ERR( rval );
195 
196  std::stringstream ss;
197  ss << "partial" << rank << ".vtk";
198  mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
199  EntityHandle outputSet;
200  rval = mb.create_meshset( MESHSET_SET, outputSet );CHECK_ERR( rval );
201  rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );CHECK_ERR( rval );
202 
203  IntxAreaUtils areaAdaptor;
204  // std::string opts_write("PARALLEL=WRITE_PART");
205  // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
206  // std::string opts_write("");
207  std::stringstream outf;
208  outf << "intersect" << rank << ".h5m";
209  rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
210  double intx_area = areaAdaptor.area_on_sphere( &mb, outputSet, radius );
211  double arrival_area = areaAdaptor.area_on_sphere( &mb, euler_set, radius );
212  std::cout << "On rank : " << rank << " arrival area: " << arrival_area << " intersection area:" << intx_area
213  << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";CHECK_ERR( rval );
214 }
215 
216 void test_intx_mpas()
217 {
218  std::string opts = std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
219  std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" );
220  Core moab;
221  Interface& mb = moab;
222  EntityHandle euler_set;
223  ErrorCode rval;
224  rval = mb.create_meshset( MESHSET_SET, euler_set );CHECK_ERR( rval );
225  std::string example( TestDir + "unittest/" + mpas_file );
226 
227  rval = mb.load_file( example.c_str(), &euler_set, opts.c_str() );
228 
229  ParallelComm* pcomm = ParallelComm::get_pcomm( &mb, 0 );
230 
231  rval = pcomm->check_all_shared_handles();CHECK_ERR( rval );
232 
233  // everybody will get a DP tag, including the non owned entities; so exchange tags is not
234  // required for LOC (here)
235  rval = manufacture_lagrange_mesh_on_sphere( &mb, euler_set );CHECK_ERR( rval );
236 
237  int rank = pcomm->proc_config().proc_rank();
238 
239  std::stringstream ste;
240  ste << "initial" << rank << ".vtk";
241  mb.write_file( ste.str().c_str(), 0, 0, &euler_set, 1 );
242 
243  Intx2MeshOnSphere worker( &mb );
244 
245  // double radius= CubeSide/2 * sqrt(3.) ; // input
246  worker.set_radius_source_mesh( radius );
247  worker.set_radius_destination_mesh( radius );
248  worker.set_box_error( EPS1 ); //
249  // worker.SetEntityType(MBQUAD);
250 
251  worker.set_error_tolerance( radius * 1.e-8 );
252  std::cout << "error tolerance epsilon_1=" << radius * 1.e-8 << "\n";
253  // worker.locate_departure_points(euler_set);
254 
255  rval = worker.FindMaxEdges( euler_set, euler_set );
256 
257  // we need to make sure the covering set is bigger than the euler mesh
258  EntityHandle covering_lagr_set;
259  rval = mb.create_meshset( MESHSET_SET, covering_lagr_set );CHECK_ERR( rval );
260 
261  // set pcomm
262  worker.set_parallel_comm( pcomm );
263 
264  rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );CHECK_ERR( rval );
265 
266  std::stringstream ss;
267  ss << "partial" << rank << ".vtk";
268  mb.write_file( ss.str().c_str(), 0, 0, &covering_lagr_set, 1 );
269  rval = moab::IntxUtils::enforce_convexity( &mb, covering_lagr_set, rank );CHECK_ERR( rval );
270  std::stringstream ss2;
271  ss2 << "partialConvex" << rank << ".vtk";
272  mb.write_file( ss2.str().c_str(), 0, 0, &covering_lagr_set, 1 );
273  EntityHandle outputSet;
274  rval = mb.create_meshset( MESHSET_SET, outputSet );CHECK_ERR( rval );
275  rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );CHECK_ERR( rval );
276 
277  // std::string opts_write("PARALLEL=WRITE_PART");
278  // rval = mb.write_file("manuf.h5m", 0, opts_write.c_str(), &outputSet, 1);
279  // std::string opts_write("");
280  std::stringstream outf;
281  outf << "intersect" << rank << ".h5m";
282  rval = mb.write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
283 
284  IntxAreaUtils areaAdaptor;
285  double intx_area = areaAdaptor.area_on_sphere( &mb, outputSet, radius );
286  double arrival_area = areaAdaptor.area_on_sphere( &mb, euler_set, radius );
287  std::cout << "On rank : " << rank << " arrival area: " << arrival_area << " intersection area:" << intx_area
288  << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";CHECK_ERR( rval );
289 }