1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
/*
 * proj1.cpp
 *
 *  project on a sphere of radius R, delete sets if needed, and delete edges between parts
 *  (created by resolve shared ents)
 */

#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include <iostream>
#include <cmath>

#include "moab/IntxMesh/IntxUtils.hpp"
#include <cassert>
using namespace moab;

double radius = 1.;  // in m:  6371220.

int main( int argc, char** argv )
{

    bool delete_partition_sets = false;

    if( argc < 3 ) return 1;

    int index         = 1;
    char* input_mesh1 = argv[1];
    char* output      = argv[2];
    while( index < argc )
    {
        if( !strcmp( argv[index], "-R" ) )  // this is for radius to project
        {
            radius = atof( argv[++index] );
        }
        if( !strcmp( argv[index], "-DS" ) )  // delete partition sets
        {
            delete_partition_sets = true;
        }

        if( !strcmp( argv[index], "-h" ) )
        {
            std::cout << " usage: proj1 <input> <output> -R <value>  -DS (delete partition sets)\n";
            return 1;
        }
        index++;
    }

    Core moab;
    Interface& mb = moab;

    ErrorCode rval = mb.load_mesh( input_mesh1 );<--- rval is initialized

    std::cout << " -R " << radius << " input: " << input_mesh1 << "  output: " << output << "\n";

    Range verts;
    rval = mb.get_entities_by_dimension( 0, 0, verts );<--- rval is overwritten
    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

    for( int v = 0; v < count; v++ )
    {
        // EntityHandle v = verts[v];
        CartVect pos( x_ptr[v], y_ptr[v], z_ptr[v] );
        pos      = pos / pos.length();
        pos      = radius * pos;
        x_ptr[v] = pos[0];
        y_ptr[v] = pos[1];
        z_ptr[v] = pos[2];
    }

    Range edges;
    rval = mb.get_entities_by_dimension( 0, 1, edges );
    if( MB_SUCCESS != rval ) return 1;
    // write edges to a new set, and after that, write the set, delete the edges and the set
    EntityHandle sf1;
    rval = mb.create_meshset( MESHSET_SET, sf1 );
    if( MB_SUCCESS != rval ) return 1;
    rval = mb.add_entities( sf1, edges );
    if( MB_SUCCESS != rval ) return 1;
    rval = mb.write_mesh( "edgesOnly.h5m", &sf1, 1 );
    if( MB_SUCCESS != rval ) return 1;
    rval = mb.delete_entities( &sf1, 1 );
    if( MB_SUCCESS != rval ) return 1;
    mb.delete_entities( edges );
    mb.write_file( output );

    if( delete_partition_sets )
    {
        Tag par_tag;
        rval = mb.tag_get_handle( "PARALLEL_PARTITION", par_tag );
        if( MB_SUCCESS == rval )

        {
            Range par_sets;
            rval =<--- Variable 'rval' is assigned a value that is never used.
                mb.get_entities_by_type_and_tag( 0, MBENTITYSET, &par_tag, NULL, 1, par_sets, moab::Interface::UNION );
            if( !par_sets.empty() ) mb.delete_entities( par_sets );
            mb.tag_delete( par_tag );
        }
    }

    return 0;
}