intx_imesh.cpp File Reference

#include <stdio.h>
#include <string.h>
#include "moab_mpi.h"
#include "iMeshP.h"

Go to the source code of this file.

Defines
#define	IMESH_ASSERT(ierr)   if (ierr!=0) printf("imesh assert\n");
#define	IMESH_NULL   0
Functions
void	update_tracer (iMesh_Instance instance, iBase_EntitySetHandle eulerSet, int *ierr)
int	main (int argc, char *argv[])

---

Define Documentation

#define IMESH_ASSERT

(

ierr

)

if (ierr!=0) printf("imesh assert\n");

Definition at line 13 of file intx_imesh.cpp.

#define IMESH_NULL   0

Definition at line 14 of file intx_imesh.cpp.

---

Function Documentation

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 26 of file intx_imesh.cpp.

                                {
  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;

  iMesh_createEntSet(imesh, 0, &euler_set, &ierr);
  IMESH_ASSERT(ierr);

  iBase_EntityHandle *cells = NULL;
  int ents_alloc = 0;
  int ents_size = 0;

  iMesh_getEntities(imesh, root, iBase_FACE, iMesh_ALL_TOPOLOGIES, &cells,
      &ents_alloc, &ents_size, &ierr);
  IMESH_ASSERT(ierr);

  iMesh_addEntArrToSet(imesh, cells, ents_size, euler_set, &ierr);

  IMESH_ASSERT(ierr);

  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

}

void update_tracer	(	iMesh_Instance	instance,
		iBase_EntitySetHandle	eulerSet,
		int *	ierr
	)

Definition at line 26 of file wrap_intx.cpp.

{
  Range ents;
  moab::Interface * mb =MOABI;
  *ierr =1;


  EntityHandle euler_set = (EntityHandle) imesh_euler_set;

  Intx2MeshOnSphere worker(mb);
  worker.SetRadius(radius);

  worker.SetErrorTolerance(gtol);

  EntityHandle covering_lagr_set;

  ErrorCode 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 ");
  }

  //
  rval = enforce_convexity(mb, covering_lagr_set);
  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;
}