1 new commit in MOAB: https://bitbucket.org/fathomteam/moab/commits/1791bb594f3f/ Changeset: 1791bb594f3f Branch: master User: iulian07 Date: 2013-12-31 23:33:07 Summary: transport example in imeshp and fortran first, add a DP tag to a partitioned mesh using create_dp then, use intx_imesh or advection to interface with mbcslam intersection wrap_intx contains the callable c imesh function fortran example might need cleanup, with respect to passing the imeshp arguments Affected #: 4 files diff --git a/tools/mbcslam/advection.F90 b/tools/mbcslam/advection.F90 new file mode 100644 index 0000000..fb102ba --- /dev/null +++ b/tools/mbcslam/advection.F90 @@ -0,0 +1,122 @@ +! advection: do a one time step transport, using c binding module +! example of how to do it in parallel in fortran +! This program shows how to load in parallel in moab from Fortran90, and how +! to call the transport wrapper. It is a fortran equivalent of intx_imesh +! driver +! +! Usage: advection + +#define CHECK(a) if (ierr .ne. 0) print *, a +program advection + + use ISO_C_BINDING + implicit none + +#include "mpif.h" + +#ifdef USE_MPI +# include "iMeshP_f.h" +#else +# include "iMesh_f.h" +#endif + +! extern void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr); + INTERFACE + SUBROUTINE update_tracer ( instance , opEulerSet, ierr) bind(C) + use ISO_C_BINDING + implicit none + iMesh_Instance, INTENT(IN) , VALUE :: instance + iBase_EntityHandle, INTENT(OUT) :: opEulerSet + integer(c_int) , INTENT (OUT) :: ierr + END SUBROUTINE update_tracer + END INTERFACE + + ! declarations + ! imesh is the instance handle + iMesh_Instance imesh + + ! ents, verts will be arrays storing vertex/entity handles + iBase_EntityHandle, pointer :: ents, verts + iBase_EntitySetHandle root_set + iBase_EntitySetHandle opEulerSet + CHARACTER (LEN=200) options + CHARACTER (LEN=200) filename + CHARACTER (LEN=200) optionswrite + CHARACTER (LEN=200) outname + TYPE(C_PTR) :: vertsPtr, entsPtr + + integer rank, sz, ierr + integer lenopt, lenname + +#ifdef USE_MPI + ! local variables for parallel runs + iMeshP_PartitionHandle imeshp + IBASE_HANDLE_T mpi_comm_c +#endif + + ! init the parallel partition + call MPI_INIT(ierr) + CHECK("fail to initialize MPI") + call MPI_COMM_SIZE(MPI_COMM_WORLD, sz, ierr) + CHECK("fail to get MPI size") + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) + CHECK("fail to get MPI rank") + + ! now load the mesh; this also initializes parallel sharing + imesh = 0 + imeshp = 0 + call iMesh_newMesh("MOAB", imesh, ierr) + CHECK("fail to initialize imesh") + + call iMesh_getRootSet(%VAL(imesh), root_set, ierr) + CHECK("fail to get root set") + + call iMeshP_getCommunicator(%VAL(imesh), MPI_COMM_WORLD, mpi_comm_c, ierr) + + call iMeshP_createPartitionAll(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr) + CHECK("fail to create parallel partition ") + options = " moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION" // & + " moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE " +! " moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION " & +! " moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE ", & ! options + if (0 .eq. rank) then + print *, "load in parallel file HN16DP.h5m" + endif + lenname = 11; + lenopt = 123 + filename = "HN16DP.h5m" + call iMeshP_loadAll(%VAL(imesh), & + %VAL(imeshp), & + %VAL(root_set), & + filename, & ! filename + options, & !options + ierr, & + %VAL(lenname), & ! strlen(filename), + %VAL(lenopt) ) !119) !strlen(options)); + CHECK("fail to load mesh in parallel ") + + call update_tracer(imesh, opEulerSet, ierr) + CHECK("fail to update tracer ") + + outname = "outF.h5m"; + optionswrite = " moab:PARALLEL=WRITE_PART " ; + lenname = 8 + lenopt = 27 + call iMeshP_saveAll( %VAL(imesh), & + %VAL(imeshp), & + %VAL(opEulerSet), & + outname, & + optionswrite, & + ierr, & + %VAL(lenname), & ! strlen(filename), + %VAL(lenopt) ) !119) !strlen(options)); + CHECK(" can't save ") + + if (0==rank) then + print *, "Done" + endif + + call MPI_FINALIZE(ierr) + stop +end program advection + diff --git a/tools/mbcslam/create_dp.cpp b/tools/mbcslam/create_dp.cpp new file mode 100644 index 0000000..bd2459b --- /dev/null +++ b/tools/mbcslam/create_dp.cpp @@ -0,0 +1,273 @@ +/* + * create_dp.cpp + * + * Created on: Dec 12, 2013 + * just to add a "DP" tag to an already partitioned mesh, based on some formula + * it is launched in serial + */ +#include "moab/Core.hpp" +#include "moab/Interface.hpp" +#include <iostream> +#include <math.h> +#include <TestUtil.hpp> + +#include "CslamUtils.hpp" +#include <assert.h> +using namespace moab; + +double radius = 1.;// in m: 6371220. +int field_type = 1; + +ErrorCode add_field_value(Interface & mb) +{ + ErrorCode rval = MB_SUCCESS; + + Tag tagTracer = 0; + std::string tag_name("Tracer"); + rval = mb.tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagTracer, MB_TAG_DENSE | MB_TAG_CREAT); + CHECK_ERR(rval); + + // tagElem is the average computed at each element, from nodal values + Tag tagElem = 0; + std::string tag_name2("TracerAverage"); + rval = mb.tag_get_handle(tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT); + CHECK_ERR(rval); + + Tag tagArea = 0; + std::string tag_name4("Area"); + rval = mb.tag_get_handle(tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT); + CHECK_ERR(rval); + + /* + * get all plys first, then vertices, then move them on the surface of the sphere + * radius is 1., most of the time + * + */ + Range polygons; + rval = mb.get_entities_by_dimension(0, 2, polygons); + if (MB_SUCCESS != rval) + return rval; + + Range connecVerts; + rval = mb.get_connectivity(polygons, connecVerts); + if (MB_SUCCESS != rval) + return rval; + + + + void *data; // pointer to the LOC in memory, for each vertex + int count; + + rval = mb.tag_iterate(tagTracer, connecVerts.begin(), connecVerts.end(), count, data); + CHECK_ERR(rval); + // here we are checking contiguity + assert(count == (int) connecVerts.size()); + double * ptr_DP=(double*)data; + // lambda is for longitude, theta for latitude + // param will be: (la1, te1), (la2, te2), b, c; hmax=1, r=1/2 + // nondivergent flow, page 5, case 1, (la1, te1) = (M_PI, M_PI/3) + // (la2, te2) = (M_PI, -M_PI/3) + // la1, te1 la2 te2 b c hmax r + if (field_type==1) // quasi smooth + { + double params[] = { M_PI, M_PI/3, M_PI, -M_PI/3, 0.1, 0.9, 1., 0.5}; + for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) + { + EntityHandle oldV=*vit; + CartVect posi; + rval = mb.get_coords(&oldV, 1, &(posi[0]) ); + CHECK_ERR(rval); + + SphereCoords sphCoord = cart_to_spherical(posi); + + ptr_DP[0]=quasi_smooth_field(sphCoord.lon, sphCoord.lat, params);; + + ptr_DP++; // increment to the next node + } + } + else if (2 == field_type) // smooth + { + CartVect p1, p2; + SphereCoords spr; + spr.R = 1; + spr.lat = M_PI/3; + spr.lon= M_PI; + p1 = spherical_to_cart(spr); + spr.lat = -M_PI/3; + p2 = spherical_to_cart(spr); + // x1, y1, z1, x2, y2, z2, h_max, b0 + double params[] = { p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], 1, 5.}; + for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) + { + EntityHandle oldV=*vit; + CartVect posi; + rval = mb.get_coords(&oldV, 1, &(posi[0]) ); + CHECK_ERR(rval); + + SphereCoords sphCoord = cart_to_spherical(posi); + + ptr_DP[0]=smooth_field(sphCoord.lon, sphCoord.lat, params);; + + ptr_DP++; // increment to the next node + } + } + else if (3 == field_type) // slotted + { + // la1, te1, la2, te2, b, c, r + double params[] = { M_PI, M_PI/3, M_PI, -M_PI/3, 0.1, 0.9, 0.5};// no h_max + for (Range::iterator vit=connecVerts.begin();vit!=connecVerts.end(); vit++ ) + { + EntityHandle oldV=*vit; + CartVect posi; + rval = mb.get_coords(&oldV, 1, &(posi[0]) ); + CHECK_ERR(rval); + + SphereCoords sphCoord = cart_to_spherical(posi); + + ptr_DP[0]=slotted_cylinder_field(sphCoord.lon, sphCoord.lat, params);; + + ptr_DP++; // increment to the next node + } + } + + // add average value for quad/polygon (average corners) + // do some averages + + + Range::iterator iter = polygons.begin(); + double local_mass = 0.; // this is total mass on one proc + while (iter != polygons.end()) + { + rval = mb.tag_iterate(tagElem, iter, polygons.end(), count, data); + CHECK_ERR(rval); + double * ptr=(double*)data; + + rval = mb.tag_iterate(tagArea, iter, polygons.end(), count, data); + CHECK_ERR(rval); + double * ptrArea=(double*)data; + for (int i=0; i<count; i++, iter++, ptr++, ptrArea++) + { + const moab::EntityHandle * conn = NULL; + int num_nodes = 0; + rval = mb.get_connectivity(*iter, conn, num_nodes); + CHECK_ERR(rval); + if (num_nodes==0) + return MB_FAILURE; + std::vector<double> nodeVals(num_nodes); + double average=0.; + rval = mb.tag_get_data(tagTracer, conn, num_nodes, &nodeVals[0] ); + CHECK_ERR(rval); + for (int j=0; j<num_nodes; j++) + average+=nodeVals[j]; + average/=num_nodes; + *ptr = average; + + // now get area + std::vector<double> coords; + coords.resize(3*num_nodes); + rval = mb.get_coords(conn, num_nodes, &coords[0]); + CHECK_ERR(rval); + *ptrArea = area_spherical_polygon_lHuiller (&coords[0], num_nodes, radius); + + // we should have used some + // total mass: + local_mass += *ptrArea * average; + } + + } + + // now we can delete the tags? not yet + return MB_SUCCESS; +} + +int main(int argc, char **argv) +{ + + if (argc < 3) + { + std::cout << " usage: create_dp <input><output> -t <time> -dt <delta_t> -h \n"; + return 1; + } + + double dt=0.1; + double t=0.1; // corresponding to diffusion first step + + int index = 2; + char * input_mesh1 = argv[1]; + char * output = argv[2]; + while (index < argc) + { + if (!strcmp(argv[index], "-t")) // this is for radius to project + { + t = atof(argv[++index]); + } + if (!strcmp(argv[index], "-dt")) // delete partition sets + { + dt = atof(argv[++index]); + } + + if (!strcmp(argv[index], "-h")) + { + std::cout << " usage: create_dp <input><output> -t <time> -dt <delta_t> -h \n"; + return 1; + } + index++; + } + + Core moab; + Interface & mb = moab; + + ErrorCode rval; + + rval = mb.load_mesh(input_mesh1); + + std::cout << " -t " << t << " -dt " << dt << " input: " << input_mesh1 << + " output: " << output << "\n"; + + Range verts; + rval = mb.get_entities_by_dimension(0, 0, verts); + if (MB_SUCCESS != rval) + return 1; + + double *x_ptr, *y_ptr, *z_ptr; + int count; + rval = mb.coords_iterate(verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count); + if (MB_SUCCESS != rval) + return 1; + assert(count == (int) verts.size()); // should end up with just one contiguous chunk of vertices + + Tag tagh = 0; + std::string tag_name("DP"); + rval = mb.tag_get_handle(tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT); + CHECK_ERR(rval); + void *data; // pointer to the LOC in memory, for each vertex + int count_tag; + + rval = mb.tag_iterate(tagh, verts.begin(), verts.end(), count_tag, data); + CHECK_ERR(rval); + // here we are checking contiguity + assert(count_tag == (int) verts.size()); + double * ptr_DP=(double*)data; + + for (int v = 0; v < count; v++) { + //EntityHandle v = verts[v]; + CartVect pos( x_ptr[v], y_ptr[v] , z_ptr[v]); + CartVect newPos; + departure_point_case1(pos, t, dt, newPos); + ptr_DP[0]=newPos[0]; + ptr_DP[1]=newPos[1]; + ptr_DP[2]=newPos[2]; + ptr_DP+=3; // increment to the next vertex + } + + + rval = add_field_value(mb); + + mb.write_file(output); + + return 0; +} + + + + diff --git a/tools/mbcslam/intx_imesh.cpp b/tools/mbcslam/intx_imesh.cpp new file mode 100644 index 0000000..0c49ac4 --- /dev/null +++ b/tools/mbcslam/intx_imesh.cpp @@ -0,0 +1,92 @@ + +/* + * This program updates a manufactured tracer field from time T0 to time T1, in parallel. + * Input: arrival mesh, already distributed on processors, and a departure position for + * each vertex, saved in a tag DP + */ +#include <stdio.h> +#include <string.h> +#include "moab_mpi.h" +#include "iMeshP.h" + + +#define IMESH_ASSERT(ierr) if (ierr!=0) printf("imesh assert\n"); +#define IMESH_NULL 0 + +#ifdef __cplusplus +extern "C" { +#endif + void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr); + +#ifdef __cplusplus +} // extern "C" +#endif + + +int main(int argc, char* argv[]){ + MPI_Init(&argc, &argv); + + iMesh_Instance imesh; + iMeshP_PartitionHandle partn; + int ierr, num_sets; + + iBase_EntitySetHandle root; + imesh = IMESH_NULL; + iMesh_newMesh(0, &imesh, &ierr, 0); + IMESH_ASSERT(ierr); + iMesh_getRootSet( imesh, &root, &ierr ); + IMESH_ASSERT(ierr); + + iMeshP_createPartitionAll(imesh, MPI_COMM_WORLD, &partn, &ierr); + int rank, size; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &size); + IMESH_ASSERT(ierr); + + const char options[] = " moab:PARALLEL=READ_PART " + " moab:PARTITION=PARALLEL_PARTITION " + " moab:PARALLEL_RESOLVE_SHARED_ENTS " + " moab:PARTITION_DISTRIBUTE "; + const char * filename = "HN16DP.h5m"; // the file should have the dp tag already + + if (0==rank) + printf("load in parallel the file: %s \n", filename); + iMeshP_loadAll(imesh, + partn, + root, + filename, + options, + &ierr, + strlen(filename), + strlen(options)); + IMESH_ASSERT(ierr); + + + iMesh_getNumEntSets(imesh, + IMESH_NULL, + 1, + &num_sets, + &ierr); + IMESH_ASSERT(ierr); + printf("There's %d entity sets here on process rank %d \n", num_sets, rank); + + iBase_EntitySetHandle euler_set; + update_tracer( imesh, &euler_set, &ierr); + IMESH_ASSERT(ierr); + + // write everything + const char * out_name = "out.h5m"; + const char optionswrite[] = " moab:PARALLEL=WRITE_PART " ; + iMeshP_saveAll( imesh, partn, euler_set, + out_name, + optionswrite, + &ierr, + strlen(out_name), + strlen(optionswrite)); + IMESH_ASSERT(ierr); + + if (0==rank) + printf("Done\n"); + MPI_Finalize(); //probably the 4th time this is called.. no big deal + +} diff --git a/tools/mbcslam/wrap_intx.cpp b/tools/mbcslam/wrap_intx.cpp new file mode 100644 index 0000000..fa6bc20 --- /dev/null +++ b/tools/mbcslam/wrap_intx.cpp @@ -0,0 +1,104 @@ +/* + * wrap_intx.cpp + * Will implement the intersection method that will be callable from fortran too + * will be added to the library mbcslam.a + * + * + * Created on: Dec 14, 2013 + * Author: iulian + */ +#include "iMesh.h" +#include "iMeshP.h" +#include "MBiMesh.hpp" +#include "moab/Core.hpp" +#include "moab/Range.hpp" +#include "Intx2MeshOnSphere.hpp" + +using namespace moab; +double radius = 1.; +double gtol = 1.e-9; +bool debug = true; + +#ifdef __cplusplus +extern "C" { +#endif + +void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr) +{ + Range ents; + moab::Interface * mb =MOABI; + *ierr =1; + ErrorCode rval = mb->get_entities_by_dimension(0, 2, ents);// root 0 + + ERRORV(rval, "can't get all 2d entities from root"); + + EntityHandle euler_set; + //EntityHandle lagr_set; + + rval = mb->create_meshset(MESHSET_SET, euler_set); + ERRORV(rval , "can't create arrival mesh set"); + + *opEulerSet = (iBase_EntitySetHandle)euler_set; + + rval = mb->add_entities(euler_set, ents); + ERRORV(rval , "can't add ents to arrival set"); + + Intx2MeshOnSphere worker(mb); + worker.SetRadius(radius); + + worker.SetErrorTolerance(gtol); + + EntityHandle covering_lagr_set; + rval = mb->create_meshset(MESHSET_SET, covering_lagr_set); + ERRORV(rval , "can't create covering set "); + + // we need to update the correlation tag and remote tuples + rval = worker.create_departure_mesh_2nd_alg(euler_set, covering_lagr_set); + ERRORV(rval , "can't populate covering set "); + + if (debug) + { + rval = mb->write_file("lagr.h5m", 0, 0, &covering_lagr_set, 1 ); + ERRORV(rval , "can't write covering set "); + } + + EntityHandle outputSet; + rval = mb->create_meshset(MESHSET_SET, outputSet); + ERRORV(rval , "can't create output set "); + + rval = worker.intersect_meshes(covering_lagr_set, euler_set, outputSet); + ERRORV(rval , "can't intersect "); + + if (debug) + { + rval = mb->write_file("output.vtk", 0, 0, &outputSet, 1 ); + ERRORV(rval , "can't write covering set "); + } + + // tagElem is the average computed at each element, from nodal values + Tag tagElem = 0; + std::string tag_name2("TracerAverage"); + rval = mb->tag_get_handle(tag_name2.c_str(), 1, MB_TYPE_DOUBLE, tagElem, MB_TAG_DENSE | MB_TAG_CREAT); + ERRORV(rval , "can't get tracer tag "); + + // area of the euler element is fixed, store it; it is used to recompute the averages at each + // time step + Tag tagArea = 0; + std::string tag_name4("Area"); + rval = mb->tag_get_handle(tag_name4.c_str(), 1, MB_TYPE_DOUBLE, tagArea, MB_TAG_DENSE | MB_TAG_CREAT); + ERRORV(rval , "can't get area tag"); + + rval = worker.update_tracer_data(outputSet, tagElem, tagArea); + ERRORV(rval , "can't update tracer "); + + // everything can be deleted now from intx data; polygons, etc. + + *ierr = 0; + return; +} + +#ifdef __cplusplus +} // extern "C" +#endif + + Repository URL: https://bitbucket.org/fathomteam/moab/ -- This is a commit notification from bitbucket.org. You are receiving this because you have the service enabled, addressing the recipient of this email.