measure.cpp

Go to the documentation of this file.

    1 /*
    2  * Copyright (c) 2005 Lawrence Livermore National Laboratory under
    3  * contract number B545069 with the University of Wisconsin - Madison.
    4  *
    5  * Redistribution and use in source and binary forms, with or without
    6  * modification, are permitted provided that the following conditions
    7  * are met:
    8  * 1. Redistributions of source code must retain the above copyright
    9  *    notice, this list of conditions and the following disclaimer.
   10  * 2. Redistributions in binary form must reproduce the above copyright
   11  *    notice, this list of conditions and the following disclaimer in the
   12  *    documentation and/or other materials provided with the distribution.
   13  *
   14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
   15  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   16  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
   17  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
   18  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
   19  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
   20  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
   21  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
   22  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
   23  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   24  */
   25  
   26 #include "moab/CartVect.hpp"
   27 #include "measure.hpp"
   28 using namespace moab;
   29  
   30 inline static double tet_volume( const CartVect& v0, const CartVect& v1, const CartVect& v2, const CartVect& v3 )
   31 {
   32     return 1. / 6. * ( ( ( v1 - v0 ) * ( v2 - v0 ) ) % ( v3 - v0 ) );
   33 }
   34  
   35 double edge_length( const double* start_vtx_coords, const double* end_vtx_coords )
   36 {
   37     const CartVect* start = reinterpret_cast< const CartVect* >( start_vtx_coords );
   38     const CartVect* end   = reinterpret_cast< const CartVect* >( end_vtx_coords );
   39     return ( *start - *end ).length();
   40 }
   41  
   42 double measure( moab::EntityType type, int num_vertices, const double* vertex_coordinates )
   43 {
   44     const CartVect* coords = reinterpret_cast< const CartVect* >( vertex_coordinates );
   45     switch( type )
   46     {
   47         case moab::MBEDGE:
   48             return ( coords[0] - coords[1] ).length();
   49         case moab::MBTRI:
   50             return 0.5 * ( ( coords[1] - coords[0] ) * ( coords[2] - coords[0] ) ).length();
   51         case moab::MBQUAD:
   52             num_vertices = 4;
   53             // fall through
   54         case moab::MBPOLYGON: {
   55             CartVect mid( 0, 0, 0 );
   56             for( int i = 0; i < num_vertices; ++i )
   57                 mid += coords[i];
   58             mid /= num_vertices;
   59  
   60             double sum = 0.0;
   61             for( int i = 0; i < num_vertices; ++i )
   62             {
   63                 int j = ( i + 1 ) % num_vertices;
   64                 sum += ( ( mid - coords[i] ) * ( mid - coords[j] ) ).length();
   65             }
   66             return 0.5 * sum;
   67         }
   68         case moab::MBTET:
   69             return tet_volume( coords[0], coords[1], coords[2], coords[3] );
   70         case moab::MBPYRAMID:
   71             return tet_volume( coords[0], coords[1], coords[2], coords[4] ) +
   72                    tet_volume( coords[0], coords[2], coords[3], coords[4] );
   73         case moab::MBPRISM:
   74             return tet_volume( coords[0], coords[1], coords[2], coords[5] ) +
   75                    tet_volume( coords[3], coords[5], coords[4], coords[0] ) +
   76                    tet_volume( coords[1], coords[4], coords[5], coords[0] );
   77         case moab::MBHEX:
   78             return tet_volume( coords[0], coords[1], coords[3], coords[4] ) +
   79                    tet_volume( coords[7], coords[3], coords[6], coords[4] ) +
   80                    tet_volume( coords[4], coords[5], coords[1], coords[6] ) +
   81                    tet_volume( coords[1], coords[6], coords[3], coords[4] ) +
   82                    tet_volume( coords[2], coords[6], coords[3], coords[1] );
   83         default:
   84             return 0.0;
   85     }
   86 }