gmv_io.c

一个用来实现偏微分方程中网格的计算库

// $Id: gmv_io.C 2917 2008-07-08 19:41:53Z jwpeterson $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// Changes: 
// o no more subelements, all elements are written down to GMV directly
// o Some nodes have to be left out, eg node 8 in QUAD9 
// o


// C++ includes
#include <iomanip>
#include <fstream>
#include <cstring> // for std::strcpy, std::memcpy
#include <cstdio>  // for std::sprintf
#include <vector>

// Local includes
#include "libmesh_config.h"
#include "libmesh_logging.h"
#include "gmv_io.h"
#include "mesh_base.h"
#include "elem.h"
#include "equation_systems.h"
#include "numeric_vector.h"

// Wrap everything in a GMV namespace and
// use extern "C" to avoid name mangling.
#ifdef HAVE_GMV
namespace GMV
{
  extern "C"
  {
#include "gmvread.h"
  }
}
#endif

// anonymous namespace to hold local data
namespace
{
  /**
   * Defines mapping from libMesh element types to GMV element types.
   */
  struct elementDefinition {
    std::string label;
    std::vector<unsigned int> nodes;
  };


  // maps from a libMesh element type to the proper
  // GMV elementDefinition.  Placing the data structure
  // here in this anonymous namespace gives us the
  // benefits of a global variable without the nasty
  // side-effects
  std::map<ElemType, elementDefinition> eletypes;



  // ------------------------------------------------------------
  // helper function to initialize the eletypes map
  void init_eletypes ()
  {
    if (eletypes.empty())
      {
	// This should happen only once.  The first time this method
	// is called the eletypes data struture will be empty, and
	// we will fill it.  Any subsequent calls will find an initialized
	// eletypes map and will do nothing.

	//==============================
	// setup the element definitions
	elementDefinition eledef;

	// use "swap trick" from Scott Meyer's "Effective STL" to initialize
	// eledef.nodes vector
	
	// EDGE2
	{
	  eledef.label = "line 2";
	  const unsigned int nodes[] = {0,1};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	
	  eletypes[EDGE2] = eledef;
	}
  
	// LINE3
	{
	  eledef.label = "3line 3";
	  const unsigned int nodes[] = {0,1,2};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	  
	  eletypes[EDGE3] = eledef;
	}
      
	// TRI3
	{
	  eledef.label = "tri3 3";
	  const unsigned int nodes[] = {0,1,2};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	  
	  eletypes[TRI3] = eledef;
	}
      
	// TRI6
	{
	  eledef.label = "6tri 6";
	  const unsigned int nodes[] = {0,1,2,3,4,5};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[TRI6] = eledef;
	}
      
	// QUAD4
	{
	  eledef.label = "quad 4";
	  const unsigned int nodes[] = {0,1,2,3};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[QUAD4] = eledef;
	}
      
	// QUAD8, QUAD9
	{
	  eledef.label = "8quad 8";
	  const unsigned int nodes[] = {0,1,2,3,4,5,6,7};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[QUAD8] = eledef;
	  eletypes[QUAD9] = eledef;
	}
      
	// HEX8
	{
	  eledef.label = "phex8 8";
	  const unsigned int nodes[] = {0,1,2,3,4,5,6,7};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[HEX8] = eledef;
	}
      
	// HEX20, HEX27
	{
	  eledef.label = "phex20 20";
	  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[HEX20] = eledef;
	  eletypes[HEX27] = eledef;
	}
      
	// TET4
	{
	  eledef.label = "tet 4";
	  const unsigned int nodes[] = {0,1,2,3};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);

	  eletypes[TET4] = eledef;
	}
      
	// TET10
	{
	  eledef.label = "ptet10 10";
	  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,9};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	  
	  eletypes[TET10] = eledef;
	}
      
	// PRISM6
	{
	  eledef.label = "pprism6 6";
	  const unsigned int nodes[] = {0,1,2,3,4,5};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	  
	  eletypes[PRISM6] = eledef;
	}
      
	// PRISM15, PRISM18
	{
	  eledef.label = "pprism15 15";
	  const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,12,13,14,9,10,11};
	  const unsigned int nnodes = sizeof(nodes)/sizeof(nodes[0]);
	  std::vector<unsigned int>(nodes, nodes+nnodes).swap(eledef.nodes);
	  
	  eletypes[PRISM15] = eledef;
	  eletypes[PRISM18] = eledef;
	}
	//==============================      
      }
  }
  
} // end anonymous namespace


// ------------------------------------------------------------
// GMVIO  members
void GMVIO::write (const std::string& fname)
{
  if (libMesh::processor_id() == 0)
    {
      if (this->binary())
	this->write_binary (fname);
      else
	this->write_ascii_old_impl  (fname);
    }
}



void GMVIO::write_nodal_data (const std::string& fname,
                              const std::vector<Number>& soln,
                              const std::vector<std::string>& names)
{
  START_LOG("write_nodal_data()", "GMVIO");

  if (libMesh::processor_id() == 0)
    {
      if (this->binary())
	this->write_binary (fname, &soln, &names);
      else
	this->write_ascii_old_impl  (fname, &soln, &names);
    }
  
  STOP_LOG("write_nodal_data()", "GMVIO");
}



void GMVIO::write_ascii_new_impl (const std::string& fname,
				  const std::vector<Number>* v,
				  const std::vector<std::string>* solution_names)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  std::cerr << "WARNING:  GMVIO::write_ascii_new_impl() not infinite-element aware!"
	    << std::endl;
  here();

  this->write_ascii_old_impl (fname, v, solution_names);

#else
  
  // Open the output file stream
  std::ofstream out (fname.c_str());
  
  libmesh_assert (out.good());

  // Get a reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  unsigned int mesh_max_p_level = 0;

  // Begin interfacing with the GMV data file
  {
    out << "gmvinput ascii\n\n";

    // write the nodes
    out << "nodes " << mesh.n_nodes() << "\n";
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(0) << " ";
    out << "\n";
    
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(1) << " ";
    out << "\n";
    
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(2) << " ";
    out << "\n\n";
  }

  {
    // write the connectivity
    out << "cells " << mesh.n_active_elem() << "\n";
    
    // initialize the eletypes map
    init_eletypes();

    MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
    
    for ( ; it != end; ++it)
      {
        const Elem* elem = *it;

        mesh_max_p_level = std::max(mesh_max_p_level,
                                    elem->p_level());

	// Make sure we have a valid entry for
	// the current element type.
	libmesh_assert (eletypes.count(elem->type()));

        const elementDefinition& ele = eletypes[elem->type()];

	// The element mapper better not require any more nodes
	// than are present in the current element!
	libmesh_assert (ele.nodes.size() <= elem->n_nodes());
	
        out << ele.label << "\n";
        for (unsigned int i=0; i < ele.nodes.size(); i++)
          out << elem->node(ele.nodes[i])+1 << " ";
        out << "\n";
      }
    out << "\n";
  }
  
  // optionally write the partition information
  if (this->partitioning())
    {
      out << "material "
	  << mesh.n_partitions()
	// Note: GMV may give you errors like
	// Error, material for cell 1 is greater than 1
        // Error, material for cell 2 is greater than 1
        // Error, material for cell 3 is greater than 1
	// ... because you put the wrong number of partitions here.
	// To ensure you write the correct number of materials, call
	// mesh.recalculate_n_partitions() before writing out the
	// mesh.
	// Note: we can't call it now because the Mesh is const here and
	// it is a non-const function.
          << " 0\n";

      for (unsigned int proc=0; proc<mesh.n_partitions(); proc++)
        out << "proc_" << proc << "\n";
      
      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      // FIXME - don't we need to use an elementDefinition here? - RHS
      for ( ; it != end; ++it)
        out << (*it)->processor_id()+1 << "\n";
      out << "\n";
    }

  // If there are *any* variables at all in the system (including
  // p level, or arbitrary cell-based data)
  // to write, the gmv file needs to contain the word "variable"
  // on a line by itself.
  bool write_variable = false;

  // 1.) p-levels
  if (this->p_levels() && mesh_max_p_level)
    write_variable = true;

  // 2.) solution data
  if ((solution_names != NULL) && (v != NULL))
    write_variable = true;

  // 3.) cell-centered data
  if ( !(this->_cell_centered_data.empty()) )
    write_variable = true;
  
  if (write_variable)
    out << "variable\n";
  
//   if ((this->p_levels() && mesh_max_p_level) || 
//     ((solution_names != NULL) && (v != NULL)))
//     out << "variable\n";

  // optionally write the polynomial degree information
  if (this->p_levels() && mesh_max_p_level)
    {
      out << "p_level 0\n";

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      for ( ; it != end; ++it)
        {
          const Elem* elem = *it;

          const elementDefinition& ele = eletypes[elem->type()];

	  // The element mapper better not require any more nodes
	  // than are present in the current element!
	  libmesh_assert (ele.nodes.size() <= elem->n_nodes());

          for (unsigned int i=0; i < ele.nodes.size(); i++)
            out << elem->p_level() << " ";
        }
      out << "\n\n";
    }


  // optionally write cell-centered data
  if ( !(this->_cell_centered_data.empty()) )
    {
      std::map<std::string, const std::vector<Real>* >::iterator       it  = this->_cell_centered_data.begin();
      const std::map<std::string, const std::vector<Real>* >::iterator end = this->_cell_centered_data.end();      

      for (; it != end; ++it)
	{
	  // write out the variable name, followed by a zero.
	  out << (*it).first << " 0\n";

	  const std::vector<Real>* the_array = (*it).second;
	  
	  // Loop over active elements, write out cell data.  If second-order cells
	  // are split into sub-elements, the sub-elements inherit their parent's
	  // cell-centered data.
	  MeshBase::const_element_iterator       elem_it  = mesh.active_elements_begin();
	  const MeshBase::const_element_iterator elem_end = mesh.active_elements_end(); 
	  
	  for (; elem_it != elem_end; ++elem_it)
	    {
	      const Elem* e = *elem_it;
	      
	      // Use the element's ID to find the value.
	      libmesh_assert (e->id() < the_array->size());
	      const Real the_value = the_array->operator[](e->id());
	      
	      if (this->subdivide_second_order())
		for (unsigned int se=0; se < e->n_sub_elem(); se++)
		  out << the_value << " ";
	      else
		out << the_value << " ";
	    }
	  
	  out << "\n\n";
	}
    }

  
  // optionally write the data
  if ((solution_names != NULL) && (v != NULL))
    {      
      const unsigned int n_vars = solution_names->size();

      if (!(v->size() == mesh.n_nodes()*n_vars))
        std::cerr << "ERROR: v->size()=" << v->size()
                  << ", mesh.n_nodes()=" << mesh.n_nodes()
                  << ", n_vars=" << n_vars
                  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
                  << "\n";
      
      libmesh_assert (v->size() == mesh.n_nodes()*n_vars);

      for (unsigned int c=0; c<n_vars; c++)
        {