MOAB: Mesh Oriented datABase  (version 5.3.0)
wrap_intx.cpp
Go to the documentation of this file.
00001 /*
00002  * wrap_intx.cpp
00003  *  Will implement the intersection method that will be callable from fortran too
00004  *  will be added to the sub-library IntxMesh
00005  *
00006  *
00007  *  Created on: Dec 14, 2013
00008  *      Author: iulian
00009  */
00010 #include "iMesh.h"
00011 #include "iMeshP.h"
00012 #include "MBiMesh.hpp"
00013 #include "moab/Core.hpp"
00014 #include "moab/Range.hpp"
00015 #include "moab/Intx2MeshOnSphere.hpp"
00016 #include "moab/ReadUtilIface.hpp"
00017 #include "moab/ParallelComm.hpp"
00018 #include "MBTagConventions.hpp"
00019 #include "moab/ParallelMergeMesh.hpp"
00020 #include "moab/IntxUtils.hpp"
00021 
00022 #include <sstream>
00023 #include <mpi.h>
00024 
00025 using namespace moab;
00026 static double radius = 1.;
00027 static double gtol   = 1.e-9;
00028 static bool debug    = false;
00029 
00030 // this mapping to coordinates will keep an index into the coords and dep_coords array
00031 // more exactly, the fine vertices are in a Range fineVerts;
00032 // then the vertex index i in fineVerts , fineVerts[i] will have 3 coordinates at
00033 // index  mapping_to_coords[i] * 3
00034 // more exactly, x = coords[ mapping_to_coords[i] * 3    ]
00035 //               y = coords[ mapping_to_coords[i] * 3 + 1]
00036 //               z = coords[ mapping_to_coords[i] * 3 + 2]
00037 
00038 // in the same way, departure position for vertex depVerts[i] will be at index
00039 //   mapping_to_coords[i] * 2 in dep_coords array (it has just latitude and longitude)
00040 //
00041 static int* mapping_to_coords = NULL;
00042 static int numVertices        = 0;
00043 // this will be instantiated at create mesh step
00044 // should be cleaned up at the end
00045 static Intx2MeshOnSphere* pworker = NULL;
00046 
00047 // should get rid of this; instead of using array[NC+1][NC+1], use  row based indexing (C-style):
00048 // parray =  new int[ (NC+1)*(NC+1) ] , and instead of array[i][j], use parray[ i*(NC+1) + j ]
00049 #define NC 3
00050 
00051 /*
00052  *  methods defined here:
00053  *  void update_tracer(iMesh_Instance instance,
00054  iBase_EntitySetHandle imesh_euler_set, int * ierr);
00055 
00056  void create_mesh(iMesh_Instance instance,
00057  iBase_EntitySetHandle * imesh_euler_set, double * coords, int * corners,
00058  int nc, int nelem, MPI_Fint comm, int * ierr) ;
00059 
00060  void intersection_at_level(iMesh_Instance instance,
00061  iBase_EntitySetHandle fine_set, iBase_EntitySetHandle * intx_set, double * dep_coords, double
00062  radius, int nc, int nelem, MPI_Fint comm, int * ierr)
00063 
00064  */
00065 #ifdef __cplusplus
00066 extern "C" {
00067 #endif
00068 
00069 void initialize_area_and_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, double* area_vals,
00070                                  int* ierr )
00071 {
00072 
00073     EntityHandle eul_set = (EntityHandle)imesh_euler_set;
00074 
00075     moab::Interface* mb = MOABI;
00076     *ierr               = 1;
00077 
00078     ErrorCode rval;
00079 
00080     Range eulQuads;
00081     rval = mb->get_entities_by_type( eul_set, MBQUAD, eulQuads );
00082     ERRORV( rval, "can't get eulerian quads" );
00083 
00084     // area of the euler element is fixed, store it; it is used to recompute the averages at each
00085     // time step
00086     Tag tagArea = 0;
00087     std::string tag_name4( "Area" );
00088     rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
00089     ERRORV( rval, "can't get area tag" );
00090 
00091     std::cout << " num quads = " << eulQuads.size() << "\n";
00092     ERRORV( rval, "can't set area for each quad" );
00093 
00094     rval = mb->tag_set_data( tagArea, eulQuads, &area_vals[0] );
00095 
00096     *ierr = 0;
00097     return;
00098 }
00099 
00100 void update_tracer_test( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set,
00101                          iBase_EntitySetHandle imesh_output_set, int numTracers, double* tracer_vals, int* ierr )
00102 {
00103 
00104     EntityHandle eul_set    = (EntityHandle)imesh_euler_set;
00105     EntityHandle output_set = (EntityHandle)imesh_output_set;
00106 
00107     moab::Interface* mb = MOABI;
00108     *ierr               = 1;
00109 
00110     ErrorCode rval;
00111 
00112     Range eulQuads;
00113     rval = mb->get_entities_by_type( eul_set, MBQUAD, eulQuads );
00114     ERRORV( rval, "can't get eulerian quads" );
00115 
00116     // tagElem is the average computed at each element, from nodal values
00117     Tag tagElem = 0;
00118     std::string tag_name2( "TracerAverage" );
00119     rval = mb->tag_get_handle( tag_name2.c_str(), numTracers, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );
00120     ERRORV( rval, "can't get tracer tag " );
00121 
00122     // area of the euler element is fixed, store it; it is used to recompute the averages at each
00123     // time step
00124     Tag tagArea = 0;
00125     std::string tag_name4( "Area" );
00126     rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
00127     ERRORV( rval, "can't get area tag" );
00128 
00129     std::cout << " num quads = " << eulQuads.size() << "\n";
00130 
00131     rval = mb->tag_set_data( tagElem, eulQuads, &tracer_vals[0] );
00132     ERRORV( rval, "can't set tracer data" );
00133 
00134     rval = pworker->update_tracer_data( output_set, tagElem, tagArea );
00135     ERRORV( rval, "can't update tracer " );
00136 
00137     rval = mb->tag_get_data( tagElem, eulQuads, &tracer_vals[0] );
00138     ERRORV( rval, "can't get tracer data" );
00139 
00140     *ierr = 0;
00141     return;
00142 }
00143 
00144 void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle imesh_euler_set, int* ierr )
00145 {
00146     using namespace moab;
00147     const double radius = 1.;
00148     const double gtol   = 1.e-9;
00149     const bool debug    = false;
00150 
00151     Range ents;
00152     moab::Interface* mb = MOABI;
00153     *ierr               = 1;
00154 
00155     EntityHandle euler_set = (EntityHandle)imesh_euler_set;
00156 
00157     Intx2MeshOnSphere worker( mb );
00158     worker.SetRadius( radius );
00159 
00160     worker.SetErrorTolerance( gtol );
00161 
00162     EntityHandle covering_lagr_set;
00163 
00164     ErrorCode rval = mb->create_meshset( MESHSET_SET, covering_lagr_set );
00165     ERRORV( rval, "can't create covering set " );
00166 
00167     // we need to update the correlation tag and remote tuples
00168     rval = worker.create_departure_mesh_2nd_alg( euler_set, covering_lagr_set );
00169     ERRORV( rval, "can't populate covering set " );
00170 
00171     if( debug )
00172     {
00173         rval = mb->write_file( "lagr.h5m", 0, 0, &covering_lagr_set, 1 );
00174         ERRORV( rval, "can't write covering set " );
00175     }
00176 
00177     //
00178     rval = enforce_convexity( mb, covering_lagr_set );
00179     ERRORV( rval, "can't write covering set " );
00180 
00181     EntityHandle outputSet;
00182     rval = mb->create_meshset( MESHSET_SET, outputSet );
00183     ERRORV( rval, "can't create output set " );
00184 
00185     rval = worker.intersect_meshes( covering_lagr_set, euler_set, outputSet );
00186     ERRORV( rval, "can't intersect " );
00187 
00188     if( debug )
00189     {
00190         rval = mb->write_file( "output.vtk", 0, 0, &outputSet, 1 );
00191         ERRORV( rval, "can't write covering set " );
00192     }
00193 
00194     // tagElem is the average computed at each element, from nodal values
00195     Tag tagElem = 0;
00196     std::string tag_name2( "TracerAverage" );
00197     rval = mb->tag_get_handle( tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT );
00198     ERRORV( rval, "can't get tracer tag " );
00199 
00200     // area of the euler element is fixed, store it; it is used to recompute the averages at each
00201     // time step
00202     Tag tagArea = 0;
00203     std::string tag_name4( "Area" );
00204     rval = mb->tag_get_handle( tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT );
00205     ERRORV( rval, "can't get area tag" );
00206 
00207     rval = worker.update_tracer_data( outputSet, tagElem, tagArea );
00208     ERRORV( rval, "can't update tracer " );
00209 
00210     // everything can be deleted now from intx data; polygons, etc.
00211 
00212     *ierr = 0;
00213     return;
00214 }
00215 
00216 ErrorCode create_coarse_mesh( Interface* mb, ParallelComm* pcomm, EntityHandle coarseSet, double* coords, int* corners,
00217                               int nc, int nelem, EntityHandle& start_vert, int& totalNumVertices,
00218                               int& numCornerVertices, std::vector< double* >& coordv )
00219 {
00220 
00221     int rank = pcomm->proc_config().proc_rank();
00222 
00223     Tag tagid;  // global id at corners
00224     ErrorCode rval;
00225     tagid = mb->globalId_tag();
00226 
00227     int dum_id = -1;
00228     Tag partitionTag;
00229     mb->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, partitionTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00230                         &dum_id );
00231 
00232     // there are nelem*4 corners, and 3*(nc2+1)*(nc2+1)*nelem coordinates
00233     // create first a coarse mesh,
00234     int size_corners = 4 * nelem;
00235     // int size_coords=3*(nc+1)*(nc+1)*nelem;
00236     // order first the corners array, and eliminate duplicates
00237     std::vector< int > corn1( size_corners );
00238     std::copy( corners, corners + size_corners, corn1.begin() );
00239     std::sort( corn1.begin(), corn1.end() );
00240     corn1.erase( std::unique( corn1.begin(), corn1.end() ), corn1.end() );
00241 
00242     numCornerVertices = (int)corn1.size();
00243 
00244     // these corners are actually dofs in spectral mesh, so they will vary from 0 to
00245     // np*np*nel_global (at most) they are shared on edges, corners of spectral elements, so number
00246     // of dofs is lower than this estimate still, the number of coarse vertices is close to number
00247     // of cells + 2
00248 
00249     std::map< int, int > index_map;
00250     for( size_t i = 0; i < corn1.size(); i++ )
00251         index_map[corn1[i]] = i;
00252 
00253     // estimate for the vertices that will get created, from euler formula
00254     // v-e+f = 2 for whole sphere
00255     // for one task (distributed) v-e+f = 1, for one connected region
00256     // if multiple connex regions (out of HFSC distribution from homme), v-e+f = k, k is number of
00257     // connectivity regs so in any case, e = v+f -k, where k is at least 1 (1, 2, 3, etc) so in any
00258     // case number of coarse edges is at most e_max = v+f-1
00259 
00260     // this will give an estimate for the number of "fine" vertices
00261     int e_max = nelem + numCornerVertices - 1;
00262     // total number of extra vertices will be
00263     int numVertsOnEdges = ( nc - 1 ) * e_max;
00264     int numExtraVerts   = numVertsOnEdges + ( nc - 1 ) * ( nc - 1 ) * nelem;  // internal fine vertices
00265 
00266     totalNumVertices = numCornerVertices + numExtraVerts;  // this could be overestimated, because we are not sure
00267     // about the number of edges
00268 
00269     // used to determine the if the nodes are matching at corners of elements
00270     std::vector< int > index_used( numCornerVertices, 0 );
00271 
00272     ReadUtilIface* read_iface;
00273     rval = mb->query_interface( read_iface );ERRORR( rval, "can't get query intf " );
00274 
00275     EntityHandle start_elem, *connect;
00276     // create verts, num is 2(nquads+1) because they're in a 1d row; will initialize coords in loop
00277     // over quads later
00278     rval = read_iface->get_node_coords( 3, totalNumVertices, 0, start_vert, coordv );ERRORR( rval, "can't get node coords " );
00279     // fill it up
00280     int stride = ( nc + 1 ) * ( nc + 1 );  // 16, for nc ==3
00281     // int order_in_coord[] = {0, 3, 15, 12}; // for nc=3
00282     /*
00283 
00284      *  first j, then i, so this is the order of the points in coords array, now:
00285      *
00286      *   nc(nc+1),   ...      (nc+1)*(nc+1)-1
00287      *
00288      *   2(nc+1),
00289      *   nc+1,   nc+2,            2*(nc+1)-1
00290      *   0,      1 ,    2, ..........., nc
00291      */
00292     int order_in_coord[] = { 0, nc, ( nc + 1 ) * ( nc + 1 ) - 1, nc * ( nc + 1 ) };
00293     for( int i = 0; i < size_corners; i++ )
00294     {
00295         int id    = corners[i];
00296         int index = index_map[id];
00297 
00298         int j         = i % 4;
00299         int ind_coord = i / 4;
00300         ind_coord     = ( ind_coord * stride + order_in_coord[j] ) * 3;
00301         if( index_used[index] )
00302         {
00303             if( fabs( coordv[0][index] - coords[ind_coord] ) > 0.000001 )
00304                 std::cout << " id:" << corners[i] << " i:" << i << " j:" << j << " " << coordv[0][index] << " "
00305                           << coords[ind_coord] << "\n";
00306             continue;
00307         }
00308         coordv[0][index]     = coords[ind_coord];
00309         coordv[1][index]     = coords[ind_coord + 1];
00310         coordv[2][index]     = coords[ind_coord + 2];
00311         index_used[index]    = 1;
00312         EntityHandle vertexh = start_vert + index;
00313         rval                 = mb->tag_set_data( tagid, &vertexh, 1, (void*)&id );ERRORR( rval, "can't set tag id on vertex " );
00314     }
00315     // create quads; one quad for each spectral element; later, we will create fine quads
00316     rval = read_iface->get_element_connect( nelem, 4, MBQUAD, 0, start_elem, connect );ERRORR( rval, "can't create elements " );
00317 
00318     for( int i = 0; i < nelem; i++ )
00319     {
00320         for( int j = 0; j < 4; j++ )
00321         {
00322             int index_v        = index_map[corners[i * 4 + j]];
00323             connect[i * 4 + j] = start_vert + index_v;
00324         }
00325     }
00326 
00327     Range quads( start_elem, start_elem + nelem - 1 );
00328 
00329     mb->add_entities( coarseSet, quads );
00330 
00331     // create all edges adjacent to quads
00332     Range coarseEdges;
00333     rval = mb->get_adjacencies( quads, 1, true, coarseEdges, Interface::UNION );ERRORR( rval, "can't create edges " );
00334 
00335     mb->add_entities( coarseSet, coarseEdges );
00336 
00337     // see how much we overestimated the number e_max
00338     std::cout << " on rank " << rank << " e_max is " << e_max << " actual number of edges: " << coarseEdges.size()
00339               << "\n";
00340     rval = pcomm->resolve_shared_ents( coarseSet, 2, 1 );  // resolve vertices and edges
00341     ERRORR( rval, "can't resolve shared vertices and edges " );
00342 
00343     mb->tag_set_data( partitionTag, &coarseSet, 1, &rank );
00344 
00345     /* // do we have any edges?
00346      Range edges;
00347      mb->get_entities_by_dimension(0, 1, edges);
00348      if (rank == 0)
00349      std::cout << " number of edges after resolve share ents:" << edges.size()
00350      << "\n";*/
00351 
00352     rval = mb->write_file( "coarse.h5m", 0, "PARALLEL=WRITE_PART", &coarseSet, 1 );ERRORR( rval, "can't write in parallel coarse mesh" );
00353     // coarse mesh, from corners
00354 
00355     return rval;
00356 }
00357 
00358 // start_v and coordv refer to all vertices, including the coarse ones
00359 ErrorCode fill_coord_on_edges( Interface* mb, std::vector< double* >& coordv, double* coords, Range& edges,
00360                                EntityHandle start_v, Range& coarseQuads, int nc, int numCornerVertices,
00361                                Tag& fineVertOnEdgeTag )
00362 {
00363     ErrorCode rval = MB_SUCCESS;
00364 
00365     double* coordv2[3];
00366     for( int k = 0; k < 3; k++ )
00367         coordv2[k] = coordv[k] + numCornerVertices;  // they will start later
00368 
00369     assert( NC == nc );
00370 
00371     int edges_index[4][NC - 1];  // indices in the coords array for vertices, oriented positively
00372     /*
00373      *  first j, then i, so this is the order of the points in coords array, now:
00374      *
00375      *   nc(nc+1),   ...      (nc+1)*(nc+1)-1
00376      *
00377      *   2(nc+1),
00378      *   nc+1,   nc+2,            2*(nc+1)-1
00379      *   0,      1 ,    2, ..........., nc
00380      */
00381     // for first edge, oriented positive, the coords indices (*3) are 1, 2, (nc-1)
00382     // second edge                                                  2*(nc+1)-1, 3*(nc+1) -1, ...,
00383     // nc*(nc+1)-1 third edge:                                                  (nc+1) * (nc+1) -2,
00384     // ..., nc*(nc+1) +1 fourth edge:                                                 (nc-1)*(nc+1),
00385     // ..., (nc+1)
00386     for( int j = 1; j <= nc - 1; j++ )
00387     {
00388         edges_index[0][j - 1] = j;                                // for nc = 3: 1, 2
00389         edges_index[1][j - 1] = ( j + 1 ) * ( nc + 1 ) - 1;       //             7, 11
00390         edges_index[2][j - 1] = ( nc + 1 ) * ( nc + 1 ) - j - 1;  //             14, 13
00391         edges_index[3][j - 1] = ( nc - j ) * ( nc + 1 );          //             8, 4
00392     }
00393     // int num_quads=(int)coarseQuads.size();
00394     int stride = ( nc + 1 ) * ( nc + 1 );
00395     int indexv = 0;
00396     for( Range::iterator eit = edges.begin(); eit != edges.end(); ++eit )
00397     {
00398         EntityHandle edge = *eit;
00399         std::vector< EntityHandle > faces;
00400         rval = mb->get_adjacencies( &edge, 1, 2, false, faces );ERRORR( rval, "can't get adjacent faces." );
00401         if( faces.size() < 1 ) return MB_FAILURE;
00402         int sense = 0, side_number = -1, offset = -1;
00403         EntityHandle quad = faces[0];  // just consider first quad
00404         rval              = mb->side_number( quad, edge, side_number, sense, offset );ERRORR( rval, "can't get side number" );
00405         int indexq = coarseQuads.index( quad );
00406 
00407         if( indexq == -1 ) return MB_FAILURE;
00408 
00409         EntityHandle firstRefinedV = start_v + numCornerVertices + indexv;
00410         rval                       = mb->tag_set_data( fineVertOnEdgeTag, &edge, 1, &firstRefinedV );ERRORR( rval, "can't set refined vertex tag" );
00411         // copy the right coordinates from the coords array to coordv2 array
00412 
00413         double* start_quad = &coords[3 * stride * indexq];
00414         if( sense > 0 )
00415         {
00416             for( int k = 1; k <= nc - 1; k++ )
00417             {
00418                 int index_in_quad  = edges_index[side_number][k - 1] * 3;
00419                 coordv2[0][indexv] = start_quad[index_in_quad];
00420                 coordv2[1][indexv] = start_quad[index_in_quad + 1];
00421                 coordv2[2][indexv] = start_quad[index_in_quad + 2];
00422                 indexv++;
00423             }
00424         }
00425         else
00426         {
00427             // sense < 0, so we will traverse the edge in inverse sense
00428             for( int k = 1; k <= nc - 1; k++ )
00429             {
00430                 int index_in_quad  = edges_index[side_number][nc - 1 - k] * 3;
00431                 coordv2[0][indexv] = start_quad[index_in_quad];
00432                 coordv2[1][indexv] = start_quad[index_in_quad + 1];
00433                 coordv2[2][indexv] = start_quad[index_in_quad + 2];
00434                 indexv++;
00435             }
00436         }
00437     }
00438     return rval;
00439 }
00440 /*
00441  ErrorCode resolve_interior_verts_on_bound_edges(Interface * mb, ParallelComm * pcomm,
00442  Range & edges)
00443  {
00444  // edges are coarse edges;
00445  ErrorCode rval;
00446  int rank=pcomm->proc_config().proc_rank();
00447  Range sharedCoarseEdges=edges;// filter the non shared ones
00448  rval = pcomm->filter_pstatus(sharedCoarseEdges, PSTATUS_SHARED, PSTATUS_AND);ERRORR(rval, "can't filter coarse edges
00449  "); ParallelMergeMesh pmerge(pcomm, 0.0001); ErrorCode rval = pmerge.merge();
00450 
00451  return rval;
00452  }*/
00453 ErrorCode create_fine_mesh( Interface* mb, ParallelComm* pcomm, EntityHandle coarseSet, EntityHandle fine_set,
00454                             double* coords, int nc, int nelem, EntityHandle start_vert, int numCornerVertices,
00455                             std::vector< double* >& coordv )
00456 {
00457     int rank   = pcomm->proc_config().proc_rank();
00458     int stride = ( nc + 1 ) * ( nc + 1 );  // 16, for nc ==3
00459     // there are stride*3 coordinates for each coarse quad, representing the fine mesh
00460 
00461     Tag fineVertTag;
00462     // will store the first refined vertex on an edge in a entity handle tag
00463     EntityHandle def = 0;
00464     ErrorCode rval =
00465         mb->tag_get_handle( "__firstvertex", 1, MB_TYPE_HANDLE, fineVertTag, MB_TAG_DENSE | MB_TAG_CREAT, &def );ERRORR( rval, "can't get tag on first vertex" );
00466 
00467     Range coarseQuads;
00468     Range edges;
00469     rval = mb->get_entities_by_dimension( coarseSet, 2, coarseQuads );ERRORR( rval, "can't get coarse quads  " );
00470 
00471     rval = mb->get_entities_by_dimension( coarseSet, 1, edges );ERRORR( rval, "can't get coarse edges  " );
00472 
00473     Range verts;
00474     rval = mb->get_connectivity( coarseQuads, verts );ERRORR( rval, "can't get coarse vertices " );
00475 
00476     /*std::cout <<" local coarse mesh on rank " << rank << "  "<< coarseQuads.size() << " quads, "
00477      << edges.size() << " edges, " << verts.size() <<  " vertices.\n";*/
00478 
00479     int dum_id = -1;
00480     Tag partitionTag;
00481     mb->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, partitionTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00482                         &dum_id );
00483     // fine mesh, with all coordinates
00484     // std::vector<double *> coordv2;
00485 
00486     ReadUtilIface* read_iface;
00487     rval = mb->query_interface( read_iface );ERRORR( rval, "can't get query intf " );
00488 
00489     ;
00490     // create verts, (nc+1)*(nc+1)*nelem - verts.size()
00491     //
00492     /* int numVertsOnEdges = (nc-1)*(int)edges.size();
00493      int numExtraVerts = numVertsOnEdges + (nc-1)*(nc-1)*nelem; // internal fine vertices
00494      rval = read_iface->get_node_coords(3, numExtraVerts, 0,
00495      start_vert, coordv2);*/
00496     // ERRORR(rval, "can't get coords fine mesh");
00497     // fill coordinates for vertices on the edges, then vertices in the interior of coarse quads
00498     // we know that all quads are in order, their index corresponds to index in coords array
00499     rval =
00500         fill_coord_on_edges( mb, coordv, coords, edges, start_vert, coarseQuads, nc, numCornerVertices, fineVertTag );ERRORR( rval, "can't fill edges vertex coords on fine mesh" );
00501 
00502     EntityHandle start_elem, *connect;
00503     rval = read_iface->get_element_connect( nelem * nc * nc, 4, MBQUAD, 0, start_elem, connect );ERRORR( rval, "can't create elements fine mesh" );
00504 
00505     int iv     = ( nc - 1 ) * edges.size() + numCornerVertices;  // iv is in the coordv array indices
00506     start_vert = start_vert + numCornerVertices + ( nc - 1 ) * edges.size();
00507 
00508     /* // add a child to the mesh set, with the ordered vertices, as they come out in the list of
00509      coordinates EntityHandle vertSet; rval = mb->create_meshset(MESHSET_ORDERED, vertSet);ERRORR(rval, "can't create
00510      vertex set ");
00511 
00512      rval = mb->add_parent_child(fine_set, vertSet);ERRORR(rval, "can't create parent child relation between fine set
00513      and vertSet ");*/
00514 
00515     std::vector< EntityHandle > vertList;
00516     vertList.reserve( nelem * ( nc + 1 ) * ( nc + 1 ) );  // will have a list of vertices, in order
00517 
00518     // now fill coordinates on interior nodes; also mark the start for each interior vertex
00519     //  in a coarse quad
00520     for( int ie = 0; ie < nelem; ie++ )
00521     {
00522         // just fill coordinates for an array of (nc-1)*(nc-1) vertices
00523         EntityHandle firstVert = start_vert + ( nc - 1 ) * ( nc - 1 ) * ie;
00524         EntityHandle eh        = coarseQuads[ie];
00525         rval                   = mb->tag_set_data( fineVertTag, &eh, 1, &firstVert );ERRORR( rval, "can't set refined vertex tag" );
00526 
00527         int index_coords = stride * ie;
00528         for( int j = 1; j <= ( nc - 1 ); j++ )
00529         {
00530             for( int i = 1; i <= ( nc - 1 ); i++ )
00531             {
00532                 int indx2     = 3 * ( index_coords + ( nc + 1 ) * j + i );
00533                 coordv[0][iv] = coords[indx2];
00534                 coordv[1][iv] = coords[indx2 + 1];
00535                 coordv[2][iv] = coords[indx2 + 2];
00536                 iv++;
00537             }
00538         }
00539     }
00540 
00541     Range quads3( start_elem, start_elem + nelem * nc * nc - 1 );
00542 
00543     int ic = 0;
00544     for( int ie = 0; ie < nelem; ie++ )
00545     {
00546         // just fill coordinates for an array of (nc-1)*(nc-1) vertices
00547         EntityHandle arr2[NC + 1][NC + 1];  //
00548         /*
00549          *     (nc,0)         (nc,nc)
00550          *
00551          *
00552          *     (1,0)           (1,nc)
00553          *     (0,0) (0,1)     (0,nc)
00554          */
00555         EntityHandle coarseQ      = coarseQuads[ie];
00556         const EntityHandle* conn4 = NULL;
00557         int nnodes                = 0;
00558         rval                      = mb->get_connectivity( coarseQ, conn4, nnodes );ERRORR( rval, "can't get conn of coarse quad" );
00559         if( nnodes != 4 ) return MB_FAILURE;
00560 
00561         arr2[0][0]   = conn4[0];
00562         arr2[nc][0]  = conn4[1];
00563         arr2[nc][nc] = conn4[2];
00564         arr2[0][nc]  = conn4[3];
00565 
00566         // get the coarse edges
00567         std::vector< EntityHandle > aedges;
00568         rval = mb->get_adjacencies( &coarseQ, 1, 1, false, aedges );ERRORR( rval, "can't get adje edges of coarse quad" );
00569         assert( (int)aedges.size() == 4 );
00570 
00571         for( int k = 0; k < 4; k++ )
00572         {
00573             EntityHandle edh = aedges[k];
00574             /*
00575              edges_index[0][j-1] = j;                        // for nc = 3: 1, 2
00576              edges_index[1][j-1] = (j+1)*(nc+1) - 1;         //             7, 11
00577              edges_index[2][j-1] = (nc+1)*(nc+1) - j - 1;    //             14, 13
00578              edges_index[3][j-1] = (nc-j) * (nc+1) ;         //             8, 4
00579              */
00580             int sense = 0, side_number = -1, offset = -1;
00581             rval = mb->side_number( coarseQ, edh, side_number, sense, offset );ERRORR( rval, "can't get side number" );
00582             EntityHandle firstV;  // first vertex on edge, if edge oriented positively
00583             rval = mb->tag_get_data( fineVertTag, &edh, 1, &firstV );ERRORR( rval, "can't get first vertex tag on edge" );
00584             if( sense > 0 )
00585             {
00586                 if( 0 == side_number )
00587                 {
00588                     for( int i = 1; i <= nc - 1; i++ )
00589                         arr2[i][0] = firstV + i - 1;
00590                 }
00591                 else if( 1 == side_number )
00592                 {
00593                     for( int j = 1; j <= nc - 1; j++ )
00594                         arr2[nc][j] = firstV + j - 1;
00595                 }
00596                 else if( 2 == side_number )
00597                 {
00598                     for( int i = nc - 1; i >= 1; i-- )
00599                         arr2[i][nc] = firstV + nc - 1 - i;
00600                 }
00601                 else if( 3 == side_number )
00602                 {
00603                     for( int j = nc - 1; j >= 1; j-- )
00604                         arr2[0][j] = firstV + nc - 1 - j;
00605                 }
00606             }
00607             else  // if (sense<0)
00608             {
00609                 if( 0 == side_number )
00610                 {
00611                     for( int i = 1; i <= nc - 1; i++ )
00612                         arr2[i][0] = firstV + nc - 1 - i;
00613                 }
00614                 else if( 1 == side_number )
00615                 {
00616                     for( int j = 1; j <= nc - 1; j++ )
00617                         arr2[nc][j] = firstV + nc - 1 - j;
00618                 }
00619                 else if( 2 == side_number )
00620                 {
00621                     for( int i = nc - 1; i >= 1; i-- )
00622                         arr2[i][nc] = firstV + i - 1;
00623                 }
00624                 else if( 3 == side_number )
00625                 {
00626                     for( int j = nc - 1; j >= 1; j-- )
00627                         arr2[0][j] = firstV + j - 1;
00628                 }
00629             }
00630         }
00631         // fill the interior points matrix
00632         EntityHandle firstV;  // first vertex on interior of coarse quad
00633         rval = mb->tag_get_data( fineVertTag, &coarseQ, 1, &firstV );ERRORR( rval, "can't get first vertex tag on coarse tag" );
00634         int inc = 0;
00635         for( int j = 1; j <= nc - 1; j++ )
00636         {
00637             for( int i = 1; i <= nc - 1; i++ )
00638             {
00639                 arr2[i][j] = firstV + inc;
00640                 inc++;
00641             }
00642         }
00643         // fill now the matrix of quads; vertices are from 0 to (nc+1)*(nc+1)-1
00644         for( int j1 = 0; j1 < nc; j1++ )
00645         {
00646             for( int i1 = 0; i1 < nc; i1++ )
00647             {
00648                 connect[ic++] = arr2[i1][j1];          // first one
00649                 connect[ic++] = arr2[i1 + 1][j1];      //
00650                 connect[ic++] = arr2[i1 + 1][j1 + 1];  // opp diagonal
00651                 connect[ic++] = arr2[i1][j1 + 1];      //
00652             }
00653         }
00654 
00655         for( int j1 = 0; j1 <= nc; j1++ )
00656         {
00657             for( int i1 = 0; i1 <= nc; i1++ )
00658             {
00659                 vertList.push_back( arr2[i1][j1] );
00660             }
00661         }
00662     }
00663     /*
00664      rval = mb->add_entities(vertSet, &vertList[0], (int)vertList.size());ERRORR(rval,"can't add to the vert set the
00665      list of ordered vertices");*/
00666 
00667     mb->add_entities( fine_set, quads3 );
00668     // notify MOAB of the new elements
00669     rval = read_iface->update_adjacencies( start_elem, nelem * nc * nc, 4, connect );ERRORR( rval, "can't update adjacencies on fine quads" );
00670 
00671     rval = mb->tag_set_data( partitionTag, &fine_set, 1, &rank );ERRORR( rval, "can't set partition tag on fine set" );
00672 
00673     /*
00674      // delete the coarse mesh, except vertices
00675      mb->delete_entities(coarseQuads);
00676      mb->delete_entities(edges);
00677      */
00678 
00679     // the vertices on the boundary edges of the partition need to be shared and resolved
00680     ParallelMergeMesh pmerge( pcomm, 0.0001 );
00681     rval = pmerge.merge();ERRORR( rval, "can't resolve vertices on interior of boundary edges " );
00682 
00683     rval = mb->get_connectivity( quads3, verts );ERRORR( rval, "can't get vertices " );
00684 
00685     Range owned_verts = verts;
00686     rval              = pcomm->filter_pstatus( owned_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT );ERRORR( rval, "can't filter for owned vertices only" );
00687 
00688     Range entities[4];
00689     entities[0] = owned_verts;
00690     entities[2] = quads3;
00691     // assign new ids only for owned entities
00692     // will eliminate gaps in global id space for vertices
00693     rval = pcomm->assign_global_ids( entities, 2, 1, true, false );ERRORR( rval, "can't assign global ids for vertices " );
00694 
00695     /*ErrorCode ParallelComm::assign_global_ids( Range entities[],
00696      const int dimension,
00697      const int start_id,
00698      const bool parallel,
00699      const bool owned_only) */
00700 
00701     std::stringstream fff;
00702     fff << "fine0" << pcomm->proc_config().proc_rank() << ".h5m";
00703     mb->write_mesh( fff.str().c_str(), &fine_set, 1 );
00704 
00705     rval = mb->write_file( "fine.h5m", 0, "PARALLEL=WRITE_PART", &fine_set, 1 );ERRORR( rval, "can't write set 3, fine " );
00706 
00707     // we need to keep a mapping index, from the coords array to the vertex handles
00708     // so, for a given vertex entity handle, at what index in the coords array the vertex
00709     // coordinates are?
00710     // in the coords array, vertices are repeated 2 times if they are interior to a coarse edge, and
00711     // repeated 3 or 4 times if they are a corner vertex in a coarse quad
00712 
00713     numVertices       = (int)verts.size();
00714     mapping_to_coords = new int[numVertices];
00715     for( int k = 0; k < numVertices; k++ )
00716         mapping_to_coords[k] = -1;  // it means it was not located yet in vertList
00717     // vertList is parallel to the coords and dep_coords array
00718 
00719     // now loop over vertsList, and see where
00720     // vertList has nelem * (nc+1)*(nc+1) vertices; loop over them, and see where are they located
00721 
00722     for( int kk = 0; kk < (int)vertList.size(); kk++ )
00723     {
00724         EntityHandle v = vertList[kk];
00725         int index      = verts.index( v );
00726         if( -1 == index )
00727         {
00728             std::cout << " can't locate vertex " << v << " in vertex Range \n";
00729             return MB_FAILURE;
00730         }
00731         if( mapping_to_coords[index] == -1 )  // it means the vertex v was not yet encountered in the vertList
00732         {
00733             mapping_to_coords[index] = kk;
00734         }
00735     }
00736     // check that every mapping has an index different from -1
00737     for( int k = 0; k < numVertices; k++ )
00738     {
00739         if( mapping_to_coords[k] == -1 )
00740         {
00741             {
00742                 std::cout << " vertex at index " << k << " in vertex Range " << verts[k] << " is not mapped \n";
00743                 return MB_FAILURE;  //
00744             }
00745         }
00746     }
00747 
00748     return rval;
00749 }
00750 
00751 // this is called from Fortran 90
00752 void create_mesh( iMesh_Instance instance, iBase_EntitySetHandle* imesh_euler_set,
00753                   iBase_EntitySetHandle* imesh_departure_set, iBase_EntitySetHandle* imesh_intx_set, double* coords,
00754                   int* corners, int nc, int nelem, MPI_Fint comm, int* ierr )
00755 {
00756     /* double * coords=(double*) icoords;
00757      int * corners = (int*) icorners;*/
00758     *ierr            = 1;
00759     Interface* mb    = MOABI;
00760     MPI_Comm mpicomm = MPI_Comm_f2c( comm );
00761     // instantiate parallel comm now or not?
00762     ParallelComm* pcomm = new ParallelComm( mb, mpicomm );
00763 
00764     EntityHandle coarseSet;
00765     ErrorCode rval = mb->create_meshset( MESHSET_SET, coarseSet );
00766     ERRORV( rval, "can't create coarse set " );
00767 
00768     EntityHandle start_vert;
00769     int totalNumVertices;
00770     int numCornerVertices;
00771     std::vector< double* > coordv;
00772     rval = create_coarse_mesh( mb, pcomm, coarseSet, coords, corners, nc, nelem, start_vert, totalNumVertices,
00773                                numCornerVertices, coordv );
00774     ERRORV( rval, "can't create coarse set " );
00775 
00776     EntityHandle fine_set;
00777     rval = mb->create_meshset( MESHSET_SET, fine_set );
00778     ERRORV( rval, "can't create fine set " );
00779 
00780     rval = create_fine_mesh( mb, pcomm, coarseSet, fine_set, coords, nc, nelem, start_vert, numCornerVertices, coordv );
00781     ERRORV( rval, "can't create fine mesh set " );
00782 
00783     *imesh_departure_set = (iBase_EntitySetHandle)fine_set;
00784 
00785     EntityHandle euler_set;
00786     rval = mb->create_meshset( MESHSET_SET, euler_set );
00787     ERRORV( rval, "can't create moab euler set " );
00788     *imesh_euler_set = (iBase_EntitySetHandle)euler_set;
00789 
00790     // call in Intx utils
00791     // it will copy the second set from the first set
00792     rval = deep_copy_set_with_quads( mb, fine_set, euler_set );
00793     ERRORV( rval, "can't populate lagrange set " );
00794 
00795     EntityHandle intx_set;
00796     rval = mb->create_meshset( MESHSET_SET, intx_set );
00797     ERRORV( rval, "can't create output set " );
00798 
00799     *imesh_intx_set = (iBase_EntitySetHandle)intx_set;
00800 
00801     pworker = new Intx2MeshOnSphere( mb );
00802 
00803     pworker->set_box_error( 100 * gtol );
00804     Range local_verts;
00805     rval = pworker->build_processor_euler_boxes( euler_set,
00806                                                  local_verts );  // output also the local_verts
00807     ERRORV( rval, "can't compute euler boxes " );
00808     pworker->SetErrorTolerance( gtol );
00809 
00810     *ierr = 0;
00811     return;
00812 }
00813 ErrorCode set_departure_points_position( Interface* mb, EntityHandle lagrSet, double* dep_coords, double radius2 )
00814 {
00815 
00816     // the departure quads are created in the same order as the fine quads
00817     // for each coarse element, there are nc*nc fine quads
00818     // their vertices are in order
00819     Range lagr_quads;
00820     ErrorCode rval = mb->get_entities_by_type( lagrSet, MBQUAD, lagr_quads );ERRORR( rval, "can't get lagrange quads" );
00821 
00822     // get all vertices from lagr_quads
00823     Range lVerts;
00824     rval = mb->get_connectivity( lagr_quads, lVerts );ERRORR( rval, "can't get lagrangian vertices (departure)" );
00825 
00826     // they are parallel to the verts Array, they must have the same number of vertices
00827     assert( numVertices == (int)lVerts.size() );
00828 
00829     for( int i = 0; i < numVertices; i++ )
00830     {
00831         EntityHandle v = lVerts[i];
00832         int index      = mapping_to_coords[i];
00833         assert( -1 != index );
00834 
00835         SphereCoords sph;
00836         sph.R   = radius2;
00837         sph.lat = dep_coords[2 * index];
00838         sph.lon = dep_coords[2 * index + 1];
00839 
00840         CartVect depPoint = spherical_to_cart( sph );
00841         rval              = mb->set_coords( &v, 1, (double*)depPoint.array() );ERRORR( rval, "can't set position of vertex" );
00842     }
00843 
00844     return MB_SUCCESS;
00845 }
00846 void intersection_at_level( iMesh_Instance instance, iBase_EntitySetHandle fine_set, iBase_EntitySetHandle lagr_set,
00847                             iBase_EntitySetHandle intx_set, double* dep_coords, double radius2, int* ierr )
00848 {
00849     *ierr         = 1;
00850     Interface* mb = MOABI;
00851     // MPI_Comm mpicomm = MPI_Comm_f2c(comm);
00852     // instantiate parallel comm now or not?
00853 
00854     EntityHandle lagrMeshSet = (EntityHandle)lagr_set;
00855 
00856     ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
00857     if( NULL == pcomm ) return;  // error is 1
00858 
00859     // set the departure tag on the fine mesh vertices
00860     ErrorCode rval = set_departure_points_position( mb, lagrMeshSet, dep_coords, radius2 );
00861     ERRORV( rval, "can't set departure tag" );
00862     if( debug )
00863     {
00864         std::stringstream fff;
00865         fff << "lagr0" << pcomm->proc_config().proc_rank() << ".vtk";
00866         rval = mb->write_mesh( fff.str().c_str(), &lagrMeshSet, 1 );
00867         ERRORV( rval, "can't write covering set " );
00868     }
00869 
00870     // it should be done earlier
00871     pworker->SetRadius( radius );
00872 
00873     EntityHandle covering_set;
00874     rval = pworker->create_departure_mesh_3rd_alg( lagrMeshSet, covering_set );
00875     ERRORV( rval, "can't compute covering set " );
00876 
00877     if( debug )
00878     {
00879         std::stringstream fff;
00880         fff << "cover" << pcomm->proc_config().proc_rank() << ".vtk";
00881         rval = mb->write_mesh( fff.str().c_str(), &covering_set, 1 );
00882 
00883         ERRORV( rval, "can't write covering set " );
00884     }
00885     EntityHandle intxSet = (EntityHandle)intx_set;
00886     rval                 = pworker->intersect_meshes( covering_set, (EntityHandle)fine_set, intxSet );
00887     ERRORV( rval, "can't intersect " );
00888 
00889     if( debug )
00890     {
00891         std::stringstream fff;
00892         fff << "intx0" << pcomm->proc_config().proc_rank() << ".vtk";
00893         rval = mb->write_mesh( fff.str().c_str(), &intxSet, 1 );
00894         ERRORV( rval, "can't write covering set " );
00895     }
00896 
00897     return;
00898 }
00899 void cleanup_after_intersection( iMesh_Instance instance, iBase_EntitySetHandle fine_set,
00900                                  iBase_EntitySetHandle lagr_set, iBase_EntitySetHandle intx_set, int* ierr )
00901 {
00902     *ierr         = 1;
00903     Interface* mb = MOABI;
00904     // delete elements
00905     // delete now the polygons and the elements of out_set
00906     // also, all verts that are not in euler set or lagr_set
00907     Range allVerts;
00908     ErrorCode rval = mb->get_entities_by_dimension( 0, 0, allVerts );
00909     ERRORV( rval, "can't get all vertices" );
00910 
00911     Range allElems;
00912     rval = mb->get_entities_by_dimension( 0, 2, allElems );
00913     ERRORV( rval, "can't get all elems" );
00914 
00915     Range polys;  // first euler polys
00916     rval = mb->get_entities_by_dimension( (EntityHandle)fine_set, 2, polys );
00917     ERRORV( rval, "can't get all polys from fine set" );
00918 
00919     // add to polys range the lagr polys
00920     rval = mb->get_entities_by_dimension( (EntityHandle)lagr_set, 2,
00921                                           polys );  // do not delete lagr set either, with its vertices
00922     ERRORV( rval, "can't get all polys from lagr set" );
00923     // add to the connecVerts range all verts, from all initial polys
00924     Range vertsToStay;
00925     rval = mb->get_connectivity( polys, vertsToStay );
00926     ERRORV( rval, "get verts that stay" );
00927 
00928     Range todeleteVerts = subtract( allVerts, vertsToStay );
00929 
00930     Range todeleteElem = subtract( allElems, polys );  // this is coarse mesh too (if still here)
00931 
00932     // empty the out mesh set
00933     EntityHandle out_set = (EntityHandle)intx_set;
00934     rval                 = mb->clear_meshset( &out_set, 1 );
00935     ERRORV( rval, "clear mesh set" );
00936 
00937     rval = mb->delete_entities( todeleteElem );
00938     ERRORV( rval, "delete intx elements" );
00939     rval = mb->delete_entities( todeleteVerts );
00940     ERRORV( rval, "failed to delete intx vertices" );
00941 
00942     *ierr = 0;
00943     return;
00944 }
00945 
00946 void cleanup_after_simulation( int* ierr )
00947 {
00948     delete[] mapping_to_coords;
00949     numVertices = 0;
00950     *ierr       = 0;
00951     return;
00952 }
00953 #ifdef __cplusplus
00954 }  // extern "C"
00955 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines