MOAB: Mesh Oriented datABase  (version 5.2.1)
viscous_cfd.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2
00028 // -*-
00029 //
00030 //   SUMMARY:
00031 //     USAGE:
00032 //
00033 // ORIG-DATE: 19-Feb-02 at 10:57:52
00034 //  LAST-MOD: 23-Jul-03 at 18:04:37 by Thomas Leurent
00035 //
00036 //
00037 // DESCRIPTION:
00038 // ============
00039 /*! \file main.cpp
00040 
00041 describe main.cpp here
00042 
00043  */
00044 // DESCRIP-END.
00045 //
00046 #include "TestUtil.hpp"
00047 
00048 #include <iostream>
00049 using std::cout;
00050 using std::endl;
00051 #include <cstdlib>
00052 #include <vector>
00053 #include <algorithm>
00054 
00055 #include "Mesquite.hpp"
00056 #include "MsqError.hpp"
00057 #include "MeshImpl.hpp"
00058 #include "ViscousCFDTetShapeWrapper.hpp"
00059 #include "MsqVertex.hpp"
00060 
00061 using namespace MBMesquite;
00062 
00063 std::string DEFAULT_INPUT = TestDir + "/3D/vtk/tets/untangled/flat-tet-sphere.vtk";
00064 
00065 void help( const char* argv0 )
00066 {
00067     std::cerr << "Usage: " << argv0 << " [<input_file>] [<output_file>]" << std::endl
00068               << "  default input file is: " << DEFAULT_INPUT << std::endl
00069               << "  defualt is no output file" << std::endl
00070               << "  Warning: input mesh is assumed to lie in Z=5 plane" << std::endl;
00071     exit( 1 );
00072 }
00073 
00074 /* For each tet, calculate ratio of min/max dihedral angle.
00075  * Return stats for all tets */
00076 void tet_dihedral_angle_ratios( Mesh& mesh, double& ratio_min, double& ratio_avg, double& ratio_max, MsqError& err );
00077 
00078 int main( int argc, char* argv[] )
00079 {
00080     const char* input_file  = DEFAULT_INPUT.c_str();
00081     const char* output_file = NULL;
00082     switch( argc )
00083     {
00084         default:
00085             help( argv[0] );
00086         case 3:
00087             if( !strcmp( argv[2], "-h" ) ) help( argv[0] );
00088             output_file = argv[2];
00089         case 2:
00090             if( !strcmp( argv[1], "-h" ) ) help( argv[0] );
00091             input_file = argv[1];
00092         case 1:;
00093     }
00094 
00095     /* Read a VTK Mesh file */
00096     MsqPrintError err( cout );
00097     MBMesquite::MeshImpl mesh;
00098     mesh.read_vtk( input_file, err );
00099     if( err )
00100     {
00101         std::cerr << input_file << ": file read failed" << endl;
00102         return 1;
00103     }
00104 
00105     double min_ini, avg_ini, max_ini;
00106     tet_dihedral_angle_ratios( mesh, min_ini, avg_ini, max_ini, err );
00107 
00108     /* Run optimizer */
00109     ViscousCFDTetShapeWrapper smoother( 1e-2 );
00110     smoother.run_instructions( &mesh, err );
00111     if( err ) return 1;
00112 
00113     double min_fin, avg_fin, max_fin;
00114     tet_dihedral_angle_ratios( mesh, min_fin, avg_fin, max_fin, err );
00115 
00116     cout << "\tMinimum \tAverage \tMaximum " << endl
00117          << "Init\t" << min_ini << '\t' << avg_ini << '\t' << max_ini << endl
00118          << "Fini\t" << min_fin << '\t' << avg_fin << '\t' << max_fin << endl;
00119 
00120     if( output_file )
00121     {
00122         mesh.write_vtk( output_file, err );
00123         if( err )
00124         {
00125             std::cerr << output_file << ": file write failed" << endl;
00126             return 1;
00127         }
00128         else
00129         {
00130             std::cout << "Wrote file: " << output_file << endl;
00131         }
00132     }
00133 
00134     if( min_ini > min_fin || avg_ini > avg_fin || max_ini > max_fin )
00135     {
00136         std::cerr << "Dihedral handle ratio decreased" << endl;
00137         return 1;
00138     }
00139 
00140     return 0;
00141 }
00142 
00143 static inline double da( double dot )
00144 {
00145     return 180 - ( 180 / M_PI ) * acos( dot );
00146 }
00147 
00148 void tet_dihedral_angle_ratios( Mesh& mesh, double& ratio_min, double& ratio_avg, double& ratio_max, MsqError& err )
00149 {
00150     std::vector< Mesh::VertexHandle > verts;
00151     std::vector< Mesh::ElementHandle > tets;
00152     std::vector< MsqVertex > coords( 4 );
00153     std::vector< size_t > junk;
00154     mesh.get_all_elements( tets, err );
00155 
00156     ratio_min    = HUGE_VAL;
00157     ratio_max    = -HUGE_VAL;
00158     ratio_avg    = 0;
00159     size_t count = 0;
00160 
00161     for( std::vector< Mesh::ElementHandle >::iterator i = tets.begin(); i != tets.end(); ++i )
00162     {
00163 
00164         Mesh::ElementHandle e = *i;
00165         EntityTopology type;
00166         mesh.elements_get_topologies( &e, &type, 1, err );
00167         assert( !err );
00168         if( type != TETRAHEDRON ) continue;
00169 
00170         verts.clear();
00171         mesh.elements_get_attached_vertices( &e, 1, verts, junk, err );
00172         assert( !err );
00173         assert( verts.size() == 4 );
00174         mesh.vertices_get_coordinates( arrptr( verts ), arrptr( coords ), 4, err );
00175         assert( !err );
00176 
00177         Vector3D v01 = coords[1] - coords[0];
00178         Vector3D v02 = coords[2] - coords[0];
00179         Vector3D v31 = coords[1] - coords[3];
00180         Vector3D v32 = coords[2] - coords[3];
00181 
00182         Vector3D n012 = ~( v02 * v01 );
00183         Vector3D n013 = ~( v31 * v01 );
00184         Vector3D n023 = ~( v02 * v32 );
00185         Vector3D n123 = ~( v31 * v32 );
00186 
00187         double ds[] = { da( n012 % n013 ), da( n012 % n123 ), da( n012 % n023 ),
00188                         da( n013 % n023 ), da( n013 % n123 ), da( n023 % n123 ) };
00189 
00190         double ratio = *std::min_element( ds, ds + 6 ) / *std::max_element( ds, ds + 6 );
00191         if( ratio < ratio_min ) ratio_min = ratio;
00192         if( ratio > ratio_max ) ratio_max = ratio;
00193         ratio_avg += ratio;
00194         count++;
00195     }
00196 
00197     if( count ) ratio_avg /= count;
00198 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines