[MOAB-dev] commit/MOAB: iulian07: transport example in imeshp and fortran
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Dec 31 16:37:30 CST 2013
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.
More information about the moab-dev
mailing list