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