1*c4762a1bSJed Brown static char help[] = "Simple MOAB example\n\n"; 2*c4762a1bSJed Brown 3*c4762a1bSJed Brown #include <petscdmmoab.h> 4*c4762a1bSJed Brown #include "moab/ScdInterface.hpp" 5*c4762a1bSJed Brown 6*c4762a1bSJed Brown typedef struct { 7*c4762a1bSJed Brown DM dm; /* DM implementation using the MOAB interface */ 8*c4762a1bSJed Brown PetscLogEvent createMeshEvent; 9*c4762a1bSJed Brown /* Domain and mesh definition */ 10*c4762a1bSJed Brown PetscInt dim; 11*c4762a1bSJed Brown char filename[PETSC_MAX_PATH_LEN]; 12*c4762a1bSJed Brown char tagname[PETSC_MAX_PATH_LEN]; 13*c4762a1bSJed Brown } AppCtx; 14*c4762a1bSJed Brown 15*c4762a1bSJed Brown PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 16*c4762a1bSJed Brown { 17*c4762a1bSJed Brown PetscErrorCode ierr; 18*c4762a1bSJed Brown PetscBool flg; 19*c4762a1bSJed Brown 20*c4762a1bSJed Brown PetscFunctionBegin; 21*c4762a1bSJed Brown ierr = PetscStrcpy(options->filename, "");CHKERRQ(ierr); 22*c4762a1bSJed Brown ierr = PetscStrcpy(options->tagname, "petsc_tag");CHKERRQ(ierr); 23*c4762a1bSJed Brown options->dim = -1; 24*c4762a1bSJed Brown 25*c4762a1bSJed Brown ierr = PetscOptionsBegin(comm, "", "MOAB example options", "DMMOAB");CHKERRQ(ierr); 26*c4762a1bSJed Brown ierr = PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex1.cxx", options->dim, &options->dim, NULL,PETSC_DECIDE,3);CHKERRQ(ierr); 27*c4762a1bSJed Brown ierr = PetscOptionsString("-filename", "The file containing the mesh", "ex1.cxx", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 28*c4762a1bSJed Brown ierr = PetscOptionsString("-tagname", "The tag name from which to create a vector", "ex1.cxx", options->tagname, options->tagname, sizeof(options->tagname), &flg);CHKERRQ(ierr); 29*c4762a1bSJed Brown ierr = PetscOptionsEnd(); 30*c4762a1bSJed Brown 31*c4762a1bSJed Brown ierr = PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);CHKERRQ(ierr); 32*c4762a1bSJed Brown PetscFunctionReturn(0); 33*c4762a1bSJed Brown } 34*c4762a1bSJed Brown 35*c4762a1bSJed Brown PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 36*c4762a1bSJed Brown { 37*c4762a1bSJed Brown moab::Interface *iface=NULL; 38*c4762a1bSJed Brown moab::Tag tag=NULL; 39*c4762a1bSJed Brown moab::Tag ltog_tag=NULL; 40*c4762a1bSJed Brown moab::Range range; 41*c4762a1bSJed Brown PetscInt tagsize; 42*c4762a1bSJed Brown moab::ErrorCode merr; 43*c4762a1bSJed Brown PetscErrorCode ierr; 44*c4762a1bSJed Brown 45*c4762a1bSJed Brown PetscFunctionBegin; 46*c4762a1bSJed Brown ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr); 47*c4762a1bSJed Brown ierr = DMMoabCreateMoab(comm, iface, <og_tag, &range, dm);CHKERRQ(ierr); 48*c4762a1bSJed Brown std::cout << "Created DMMoab using DMMoabCreateMoab." << std::endl; 49*c4762a1bSJed Brown ierr = DMMoabGetInterface(*dm, &iface);CHKERRQ(ierr); 50*c4762a1bSJed Brown 51*c4762a1bSJed Brown // load file and get entities of requested or max dimension 52*c4762a1bSJed Brown if (strlen(user->filename) > 0) { 53*c4762a1bSJed Brown merr = iface->load_file(user->filename);MBERRNM(merr); 54*c4762a1bSJed Brown std::cout << "Read mesh from file " << user->filename << std::endl; 55*c4762a1bSJed Brown } 56*c4762a1bSJed Brown else { 57*c4762a1bSJed Brown // make a simple structured mesh 58*c4762a1bSJed Brown moab::ScdInterface *scdi; 59*c4762a1bSJed Brown merr = iface->query_interface(scdi); 60*c4762a1bSJed Brown 61*c4762a1bSJed Brown moab::ScdBox *box; 62*c4762a1bSJed Brown merr = scdi->construct_box (moab::HomCoord(0,0,0), moab::HomCoord(5,5,5), NULL, 0, box);MBERRNM(merr); 63*c4762a1bSJed Brown user->dim = 3; 64*c4762a1bSJed Brown merr = iface->set_dimension(user->dim);MBERRNM(merr); 65*c4762a1bSJed Brown std::cout << "Created structured 5x5x5 mesh." << std::endl; 66*c4762a1bSJed Brown } 67*c4762a1bSJed Brown if (-1 == user->dim) { 68*c4762a1bSJed Brown moab::Range tmp_range; 69*c4762a1bSJed Brown merr = iface->get_entities_by_handle(0, tmp_range);MBERRNM(merr); 70*c4762a1bSJed Brown if (tmp_range.empty()) { 71*c4762a1bSJed Brown MBERRNM(moab::MB_FAILURE); 72*c4762a1bSJed Brown } 73*c4762a1bSJed Brown user->dim = iface->dimension_from_handle(*tmp_range.rbegin()); 74*c4762a1bSJed Brown } 75*c4762a1bSJed Brown merr = iface->get_entities_by_dimension(0, user->dim, range);MBERRNM(merr); 76*c4762a1bSJed Brown ierr = DMMoabSetLocalVertices(*dm, &range);CHKERRQ(ierr); 77*c4762a1bSJed Brown 78*c4762a1bSJed Brown // get the requested tag and create if necessary 79*c4762a1bSJed Brown std::cout << "Creating tag with name: " << user->tagname << ";\n"; 80*c4762a1bSJed Brown merr = iface->tag_get_handle(user->tagname, 1, moab::MB_TYPE_DOUBLE, tag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE);MBERRNM(merr); 81*c4762a1bSJed Brown { 82*c4762a1bSJed Brown // initialize new tag with gids 83*c4762a1bSJed Brown std::vector<double> tag_vals(range.size()); 84*c4762a1bSJed Brown merr = iface->tag_get_data(ltog_tag, range, tag_vals.data());MBERRNM(merr); // read them into ints 85*c4762a1bSJed Brown double *dval = tag_vals.data(); int *ival = reinterpret_cast<int*>(dval); // "stretch" them into doubles, from the end 86*c4762a1bSJed Brown for (int i = tag_vals.size()-1; i >= 0; i--) dval[i] = ival[i]; 87*c4762a1bSJed Brown merr = iface->tag_set_data(tag, range, tag_vals.data());MBERRNM(merr); // write them into doubles 88*c4762a1bSJed Brown } 89*c4762a1bSJed Brown merr = iface->tag_get_length(tag, tagsize);MBERRNM(merr); 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown ierr = DMSetUp(*dm);CHKERRQ(ierr); 92*c4762a1bSJed Brown 93*c4762a1bSJed Brown // create the dmmoab and initialize its data 94*c4762a1bSJed Brown ierr = PetscObjectSetName((PetscObject) *dm, "MOAB mesh");CHKERRQ(ierr); 95*c4762a1bSJed Brown ierr = PetscLogEventEnd(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr); 96*c4762a1bSJed Brown user->dm = *dm; 97*c4762a1bSJed Brown PetscFunctionReturn(0); 98*c4762a1bSJed Brown } 99*c4762a1bSJed Brown 100*c4762a1bSJed Brown int main(int argc, char **argv) 101*c4762a1bSJed Brown { 102*c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 103*c4762a1bSJed Brown moab::ErrorCode merr; 104*c4762a1bSJed Brown PetscErrorCode ierr; 105*c4762a1bSJed Brown Vec vec; 106*c4762a1bSJed Brown moab::Interface* mbImpl=NULL; 107*c4762a1bSJed Brown moab::Tag datatag=NULL; 108*c4762a1bSJed Brown 109*c4762a1bSJed Brown ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 110*c4762a1bSJed Brown ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 111*c4762a1bSJed Brown 112*c4762a1bSJed Brown ierr = CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);CHKERRQ(ierr); /* create the MOAB dm and the mesh */ 113*c4762a1bSJed Brown 114*c4762a1bSJed Brown ierr = DMMoabGetInterface(user.dm, &mbImpl);CHKERRQ(ierr); 115*c4762a1bSJed Brown merr = mbImpl->tag_get_handle(user.tagname, datatag);MBERRNM(merr); 116*c4762a1bSJed Brown ierr = DMMoabCreateVector(user.dm, datatag, NULL, PETSC_TRUE, PETSC_FALSE,&vec);CHKERRQ(ierr); /* create a vec from user-input tag */ 117*c4762a1bSJed Brown 118*c4762a1bSJed Brown std::cout << "Created VecMoab from existing tag." << std::endl; 119*c4762a1bSJed Brown ierr = VecDestroy(&vec);CHKERRQ(ierr); 120*c4762a1bSJed Brown std::cout << "Destroyed VecMoab." << std::endl; 121*c4762a1bSJed Brown ierr = DMDestroy(&user.dm);CHKERRQ(ierr); 122*c4762a1bSJed Brown std::cout << "Destroyed DMMoab." << std::endl; 123*c4762a1bSJed Brown ierr = PetscFinalize(); 124*c4762a1bSJed Brown return ierr; 125*c4762a1bSJed Brown } 126*c4762a1bSJed Brown 127*c4762a1bSJed Brown /*TEST 128*c4762a1bSJed Brown 129*c4762a1bSJed Brown build: 130*c4762a1bSJed Brown requires: moab 131*c4762a1bSJed Brown 132*c4762a1bSJed Brown test: 133*c4762a1bSJed Brown 134*c4762a1bSJed Brown TEST*/ 135