MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/MOABConfig.h" 00002 00003 #ifdef MOAB_HAVE_MPI 00004 #include "moab_mpi.h" 00005 #endif 00006 #include "moab/iMOAB.h" 00007 #include "TestUtil.hpp" 00008 #include "moab/CpuTimer.hpp" 00009 #include "moab/ProgOptions.hpp" 00010 00011 // for malloc, free: 00012 #include <iostream> 00013 00014 #define CHECKIERR( ierr, message ) \ 00015 if( 0 != ierr ) \ 00016 { \ 00017 printf( "%s", message ); \ 00018 return 1; \ 00019 } 00020 00021 // this test will be run in serial only 00022 int main( int argc, char* argv[] ) 00023 { 00024 ErrCode ierr; 00025 std::string atmFilename, ocnFilename, lndFilename, mapFilename; 00026 00027 atmFilename = TestDir + "/wholeATM_T.h5m"; 00028 ocnFilename = TestDir + "/recMeshOcn.h5m"; 00029 lndFilename = TestDir + "/wholeLnd.h5m"; 00030 // mapFilename = TestDir + "/outCS5ICOD5_map.nc"; 00031 00032 ProgOptions opts; 00033 opts.addOpt< std::string >( "atmosphere,t", "atm mesh filename (source)", &atmFilename ); 00034 opts.addOpt< std::string >( "ocean,m", "ocean mesh filename (target)", &ocnFilename ); 00035 opts.addOpt< std::string >( "land,l", "land mesh filename (target)", &lndFilename ); 00036 00037 opts.parseCommandLine( argc, argv ); 00038 00039 { 00040 std::cout << " atm file: " << atmFilename << "\n ocn file: " << ocnFilename << "\n lnd file: " << lndFilename 00041 << "\n"; 00042 } 00043 00044 #ifdef MOAB_HAVE_MPI 00045 MPI_Init( &argc, &argv ); 00046 MPI_Comm comm = MPI_COMM_WORLD; 00047 #endif 00048 /* 00049 * MOAB needs to be initialized; A MOAB instance will be created, and will be used by each 00050 * application in this framework. There is no parallel context yet. 00051 */ 00052 ierr = iMOAB_Initialize( argc, argv ); 00053 CHECKIERR( ierr, "failed to initialize MOAB" ); 00054 00055 int atmAppID, ocnAppID, lndAppID, atmocnAppID, atmlndAppID, lndatmAppID; 00056 iMOAB_AppID atmPID = &atmAppID; 00057 iMOAB_AppID ocnPID = &ocnAppID; 00058 iMOAB_AppID lndPID = &lndAppID; 00059 iMOAB_AppID atmocnPID = &atmocnAppID; 00060 iMOAB_AppID atmlndPID = &atmlndAppID; 00061 iMOAB_AppID lndatmPID = &lndatmAppID; 00062 00063 /* 00064 * Each application has to be registered once. A mesh set and a parallel communicator will be 00065 * associated with each application. A unique application id will be returned, and will be used 00066 * for all future mesh operations/queries. 00067 */ 00068 int atmCompID = 10, ocnCompID = 20, lndCompID = 30, atmocnCompID = 100, atmlndCompID = 102, lndatmCompID = 103; 00069 ierr = iMOAB_RegisterApplication( "ATM-APP", 00070 #ifdef MOAB_HAVE_MPI 00071 &comm, 00072 #endif 00073 &atmCompID, atmPID ); 00074 CHECKIERR( ierr, "failed to register application1" ); 00075 00076 ierr = iMOAB_RegisterApplication( "OCN-APP", 00077 #ifdef MOAB_HAVE_MPI 00078 &comm, 00079 #endif 00080 &ocnCompID, ocnPID ); 00081 CHECKIERR( ierr, "failed to register application2" ); 00082 00083 ierr = iMOAB_RegisterApplication( "LND-APP", 00084 #ifdef MOAB_HAVE_MPI 00085 &comm, 00086 #endif 00087 &lndCompID, lndPID ); 00088 CHECKIERR( ierr, "failed to register application3" ); 00089 00090 ierr = iMOAB_RegisterApplication( "ATM-OCN-CPL", 00091 #ifdef MOAB_HAVE_MPI 00092 &comm, 00093 #endif 00094 &atmocnCompID, atmocnPID ); 00095 CHECKIERR( ierr, "failed to register application4" ); 00096 00097 ierr = iMOAB_RegisterApplication( "ATM-LND-CPL", 00098 #ifdef MOAB_HAVE_MPI 00099 &comm, 00100 #endif 00101 &atmlndCompID, atmlndPID ); 00102 CHECKIERR( ierr, "failed to register application5" ); 00103 00104 ierr = iMOAB_RegisterApplication( "LND-ATM-CPL", 00105 #ifdef MOAB_HAVE_MPI 00106 &comm, 00107 #endif 00108 &lndatmCompID, lndatmPID ); 00109 CHECKIERR( ierr, "failed to register application5" ); 00110 00111 const char* read_opts = ""; 00112 int num_ghost_layers = 0; 00113 /* 00114 * Loading the mesh is a parallel IO operation. Ghost layers can be exchanged too, and default 00115 * MOAB sets are augmented with ghost elements. By convention, blocks correspond to MATERIAL_SET 00116 * sets, side sets with NEUMANN_SET sets, node sets with DIRICHLET_SET sets. Each element and 00117 * vertex entity should have a GLOBAL ID tag in the file, which will be available for visible 00118 * entities 00119 */ 00120 ierr = iMOAB_LoadMesh( atmPID, atmFilename.c_str(), read_opts, &num_ghost_layers, atmFilename.size(), 00121 strlen( read_opts ) ); 00122 CHECKIERR( ierr, "failed to load mesh" ); 00123 ierr = iMOAB_LoadMesh( ocnPID, ocnFilename.c_str(), read_opts, &num_ghost_layers, ocnFilename.size(), 00124 strlen( read_opts ) ); 00125 CHECKIERR( ierr, "failed to load mesh" ); 00126 ierr = iMOAB_LoadMesh( lndPID, lndFilename.c_str(), read_opts, &num_ghost_layers, lndFilename.size(), 00127 strlen( read_opts ) ); 00128 CHECKIERR( ierr, "failed to load mesh" ); 00129 00130 int nverts[3], nelem[3]; 00131 /* 00132 * Each process in the communicator will have access to a local mesh instance, which will 00133 * contain the original cells in the local partition and ghost entities. Number of vertices, 00134 * primary cells, visible blocks, number of sidesets and nodesets boundary conditions will be 00135 * returned in size 3 arrays, for local, ghost and total numbers. 00136 */ 00137 ierr = iMOAB_GetMeshInfo( atmPID, nverts, nelem, 0, 0, 0 ); 00138 CHECKIERR( ierr, "failed to get mesh info" ); 00139 printf( "Atmosphere Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] ); 00140 00141 ierr = iMOAB_GetMeshInfo( ocnPID, nverts, nelem, 0, 0, 0 ); 00142 CHECKIERR( ierr, "failed to get mesh info" ); 00143 printf( "Ocean Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] ); 00144 00145 ierr = iMOAB_GetMeshInfo( lndPID, nverts, nelem, 0, 0, 0 ); 00146 CHECKIERR( ierr, "failed to get mesh info" ); 00147 printf( "Land Component Mesh: %d vertices and %d elements\n", nverts[0], nelem[0] ); 00148 00149 /* 00150 * The 2 tags used in this example exist in the file, already. 00151 * If a user needs a new tag, it can be defined with the same call as this one 00152 * this method, iMOAB_DefineTagStorage, will return a local index for the tag. 00153 * The name of the tag is case sensitive. 00154 * This method is collective. 00155 */ 00156 // int disc_orders[2] = {1, 1}; 00157 // const char* disc_methods[2] = {"fv", "fv"}; 00158 // const char* dof_tag_names[2] = {"GLOBAL_ID", "GLOBAL_ID"}; 00159 int disc_orders[3] = { 4, 1, 1 }; 00160 int fMonotoneTypeID = 0, fVolumetric = 0, fValidate = 1, fNoConserve = 0; 00161 00162 const std::string disc_methods[3] = { "cgll", "fv", "pcloud" }; 00163 const std::string dof_tag_names[3] = { "GLOBAL_DOFS", "GLOBAL_ID", "GLOBAL_ID" }; 00164 const std::string weights_identifiers[3] = { "scalar", "scalar_pointcloud", "scalar_conservative" }; 00165 00166 const std::string bottomTempField = "a2oTbot"; 00167 const std::string bottomTempFieldATM = "a2oTbotATM"; 00168 const std::string bottomTempProjectedField = "a2oTbot_proj"; 00169 const std::string bottomTempProjectedNCField = "a2oTbot_projnocons"; 00170 00171 // Attributes for tag data storage 00172 int tagIndex[4]; 00173 // int entTypes[2] = {1, 1}; /* both on elements; */ 00174 int tagTypes[2] = { DENSE_DOUBLE, DENSE_DOUBLE }; 00175 int atmCompNDoFs = disc_orders[0] * disc_orders[0], ocnCompNDoFs = 1 /*FV*/; 00176 00177 ierr = iMOAB_DefineTagStorage( atmPID, bottomTempField.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0], 00178 bottomTempField.size() ); 00179 CHECKIERR( ierr, "failed to define the field tag on ATM" ); 00180 00181 ierr = iMOAB_DefineTagStorage( atmPID, bottomTempFieldATM.c_str(), &tagTypes[0], &atmCompNDoFs, &tagIndex[0], 00182 bottomTempFieldATM.size() ); 00183 CHECKIERR( ierr, "failed to define the field tag on ATM" ); 00184 00185 ierr = iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[1], 00186 bottomTempProjectedField.size() ); 00187 CHECKIERR( ierr, "failed to define the field tag on OCN" ); 00188 00189 ierr = iMOAB_DefineTagStorage( ocnPID, bottomTempProjectedNCField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[2], 00190 bottomTempProjectedNCField.size() ); 00191 CHECKIERR( ierr, "failed to define the field tag on OCN" ); 00192 00193 ierr = iMOAB_DefineTagStorage( lndPID, bottomTempProjectedField.c_str(), &tagTypes[1], &ocnCompNDoFs, &tagIndex[3], 00194 bottomTempProjectedField.size() ); 00195 CHECKIERR( ierr, "failed to define the field tag on LND" ); 00196 00197 /* Next compute the mesh intersection on the sphere between the source (ATM) and target (OCN) 00198 * meshes */ 00199 ierr = iMOAB_ComputeMeshIntersectionOnSphere( atmPID, ocnPID, atmocnPID ); 00200 CHECKIERR( ierr, "failed to compute mesh intersection between ATM and OCN" ); 00201 00202 /* Next compute the mesh intersection on the sphere between the source (ATM) and target (LND) 00203 * meshes */ 00204 ierr = iMOAB_ComputePointDoFIntersection( atmPID, lndPID, atmlndPID ); 00205 CHECKIERR( ierr, "failed to compute point-cloud mapping ATM-LND" ); 00206 00207 /* Next compute the mesh intersection on the sphere between the source (LND) and target (ATM) 00208 * meshes */ 00209 ierr = iMOAB_ComputePointDoFIntersection( lndPID, atmPID, lndatmPID ); 00210 CHECKIERR( ierr, "failed to compute point-cloud mapping LND-ATM" ); 00211 00212 /* We have the mesh intersection now. Let us compute the remapping weights */ 00213 fNoConserve = 0; 00214 /* Compute the weights to preoject the solution from ATM component to OCN compoenent */ 00215 ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), disc_methods[0].c_str(), &disc_orders[0], 00216 disc_methods[1].c_str(), &disc_orders[1], &fMonotoneTypeID, &fVolumetric, 00217 &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str(), 00218 weights_identifiers[0].size(), disc_methods[0].size(), 00219 disc_methods[1].size(), dof_tag_names[0].size(), 00220 dof_tag_names[1].size() ); 00221 CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-OCN scalar " 00222 "non-conservative field" ); 00223 00224 // Let us now write the map file to disk and then read it back to test the I/O API in iMOAB 00225 #ifdef MOAB_HAVE_NETCDF 00226 { 00227 const std::string atmocn_map_file_name = "atm_ocn_map.nc"; 00228 ierr = iMOAB_WriteMappingWeightsToFile( atmocnPID, weights_identifiers[0].c_str(), atmocn_map_file_name.c_str(), 00229 weights_identifiers[0].size(), atmocn_map_file_name.size() ); 00230 CHECKIERR( ierr, "failed to write map file to disk" ); 00231 00232 const std::string intx_from_file_identifier = "map-from-file"; 00233 ierr = iMOAB_LoadMappingWeightsFromFile( atmocnPID, intx_from_file_identifier.c_str(), 00234 atmocn_map_file_name.c_str(), NULL, NULL, NULL, 00235 intx_from_file_identifier.size(), atmocn_map_file_name.size() ); 00236 CHECKIERR( ierr, "failed to load map file from disk" ); 00237 } 00238 #endif 00239 00240 /* Compute the weights to preoject the solution from ATM component to LND compoenent */ 00241 ierr = iMOAB_ComputeScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), disc_methods[0].c_str(), &disc_orders[0], 00242 disc_methods[2].c_str(), &disc_orders[2], &fMonotoneTypeID, &fVolumetric, 00243 &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[2].c_str(), 00244 weights_identifiers[1].size(), disc_methods[0].size(), 00245 disc_methods[2].size(), dof_tag_names[0].size(), 00246 dof_tag_names[2].size() ); 00247 CHECKIERR( ierr, "failed to compute remapping projection weights for ATM-LND scalar " 00248 "non-conservative field" ); 00249 00250 /* Compute the weights to preoject the solution from ATM component to LND compoenent */ 00251 ierr = iMOAB_ComputeScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(), disc_methods[2].c_str(), &disc_orders[2], 00252 disc_methods[0].c_str(), &disc_orders[0], &fMonotoneTypeID, &fVolumetric, 00253 &fNoConserve, &fValidate, dof_tag_names[2].c_str(), dof_tag_names[0].c_str(), 00254 weights_identifiers[1].size(), disc_methods[2].size(), 00255 disc_methods[0].size(), dof_tag_names[2].size(), 00256 dof_tag_names[0].size() ); 00257 CHECKIERR( ierr, "failed to compute remapping projection weights for LND-ATM scalar " 00258 "non-conservative field" ); 00259 00260 /* We have the mesh intersection now. Let us compute the remapping weights */ 00261 fNoConserve = 0; 00262 ierr = iMOAB_ComputeScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), disc_methods[0].c_str(), &disc_orders[0], 00263 disc_methods[1].c_str(), &disc_orders[1], &fMonotoneTypeID, &fVolumetric, 00264 &fNoConserve, &fValidate, dof_tag_names[0].c_str(), dof_tag_names[1].c_str(), 00265 weights_identifiers[1].size(), disc_methods[0].size(), 00266 disc_methods[1].size(), dof_tag_names[0].size(), 00267 dof_tag_names[1].size() ); 00268 CHECKIERR( ierr, "failed to compute remapping projection weights for scalar conservative field" ); 00269 00270 /* We have the remapping weights now. Let us apply the weights onto the tag we defined 00271 on the srouce mesh and get the projection on the target mesh */ 00272 ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[0].c_str(), bottomTempField.c_str(), 00273 bottomTempProjectedNCField.c_str(), weights_identifiers[0].size(), 00274 bottomTempField.size(), bottomTempProjectedNCField.size() ); 00275 CHECKIERR( ierr, "failed to apply projection weights for scalar non-conservative field" ); 00276 00277 /* We have the remapping weights now. Let us apply the weights onto the tag we defined 00278 on the srouce mesh and get the projection on the target mesh */ 00279 ierr = iMOAB_ApplyScalarProjectionWeights( atmocnPID, weights_identifiers[2].c_str(), bottomTempField.c_str(), 00280 bottomTempProjectedField.c_str(), weights_identifiers[1].size(), 00281 bottomTempField.size(), bottomTempProjectedField.size() ); 00282 CHECKIERR( ierr, "failed to apply projection weights for scalar conservative field" ); 00283 00284 /* We have the remapping weights now. Let us apply the weights onto the tag we defined 00285 on the srouce mesh and get the projection on the target mesh */ 00286 ierr = iMOAB_ApplyScalarProjectionWeights( atmlndPID, weights_identifiers[1].c_str(), bottomTempField.c_str(), 00287 bottomTempProjectedField.c_str(), weights_identifiers[1].size(), 00288 bottomTempField.size(), bottomTempProjectedField.size() ); 00289 CHECKIERR( ierr, "failed to apply projection weights for ATM-LND scalar field" ); 00290 00291 /* We have the remapping weights now. Let us apply the weights onto the tag we defined 00292 on the srouce mesh and get the projection on the target mesh */ 00293 ierr = iMOAB_ApplyScalarProjectionWeights( lndatmPID, weights_identifiers[1].c_str(), bottomTempProjectedField.c_str(), 00294 bottomTempFieldATM.c_str(), weights_identifiers[1].size(), 00295 bottomTempField.size(), bottomTempProjectedField.size() ); 00296 CHECKIERR( ierr, "failed to apply projection weights for LND-ATM scalar field" ); 00297 00298 /* 00299 * the file can be written in parallel, and it will contain additional tags defined by the user 00300 * we may extend the method to write only desired tags to the file 00301 */ 00302 { 00303 // free allocated data 00304 char outputFileAtmTgt[] = "fIntxAtmTarget.h5m"; 00305 char outputFileOcnTgt[] = "fIntxOcnTarget.h5m"; 00306 char outputFileLndTgt[] = "fIntxLndTarget.h5m"; 00307 char writeOptions[] = ""; 00308 00309 ierr = iMOAB_WriteMesh( ocnPID, outputFileOcnTgt, writeOptions, strlen( outputFileOcnTgt ), 00310 strlen( writeOptions ) ); 00311 00312 ierr = iMOAB_WriteMesh( lndPID, outputFileLndTgt, writeOptions, strlen( outputFileLndTgt ), 00313 strlen( writeOptions ) ); 00314 00315 ierr = iMOAB_WriteMesh( atmPID, outputFileAtmTgt, writeOptions, strlen( outputFileAtmTgt ), 00316 strlen( writeOptions ) ); 00317 } 00318 00319 /* 00320 * deregistering application will delete all mesh entities associated with the application and 00321 * will free allocated tag storage. 00322 */ 00323 ierr = iMOAB_DeregisterApplication( lndPID ); 00324 CHECKIERR( ierr, "failed to de-register application3" ); 00325 ierr = iMOAB_DeregisterApplication( ocnPID ); 00326 CHECKIERR( ierr, "failed to de-register application2" ); 00327 ierr = iMOAB_DeregisterApplication( atmPID ); 00328 CHECKIERR( ierr, "failed to de-register application1" ); 00329 00330 /* 00331 * this method will delete MOAB instance 00332 */ 00333 ierr = iMOAB_Finalize(); 00334 CHECKIERR( ierr, "failed to finalize MOAB" ); 00335 00336 #ifdef MOAB_HAVE_MPI 00337 MPI_Finalize(); 00338 #endif 00339 00340 return 0; 00341 }