MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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