Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // C++ includes
      21             : #include <fstream>
      22             : #include <iomanip>
      23             : #include <iostream>
      24             : #include <deque>
      25             : #include <map>
      26             : 
      27             : // Local includes
      28             : #include "libmesh/libmesh_config.h"
      29             : #include "libmesh/fro_io.h"
      30             : #include "libmesh/mesh_base.h"
      31             : #include "libmesh/boundary_info.h"
      32             : #include "libmesh/elem.h"
      33             : 
      34             : namespace libMesh
      35             : {
      36             : 
      37             : 
      38             : 
      39             : // ------------------------------------------------------------
      40             : // FroIO  members
      41           0 : void FroIO::write (const std::string & fname)
      42             : {
      43             :   // We may need to gather a DistributedMesh to output it, making that
      44             :   // const qualifier in our constructor a dirty lie
      45           0 :   MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
      46             : 
      47           0 :   if (this->mesh().processor_id() == 0)
      48             :     {
      49             :       // Open the output file stream
      50           0 :       std::ofstream out_stream (fname.c_str());
      51           0 :       libmesh_assert (out_stream.good());
      52             : 
      53             :       // Make sure it opened correctly
      54           0 :       if (!out_stream.good())
      55           0 :         libmesh_file_error(fname.c_str());
      56             : 
      57             :       // Get a reference to the mesh
      58           0 :       const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
      59             : 
      60             :       // Write the header
      61           0 :       out_stream << the_mesh.n_elem()  << " "
      62           0 :                  << the_mesh.n_nodes() << " "
      63           0 :                  << "0 0 "
      64           0 :                  << the_mesh.get_boundary_info().n_boundary_ids()  << " 1\n";
      65             : 
      66             :       // Write the nodes -- 1-based!
      67           0 :       for (auto n : make_range(the_mesh.n_nodes()))
      68           0 :         out_stream << n+1 << " \t"
      69           0 :                    << std::scientific
      70           0 :                    << std::setprecision(this->ascii_precision())
      71           0 :                    << the_mesh.point(n)(0) << " \t"
      72           0 :                    << the_mesh.point(n)(1) << " \t"
      73           0 :                    << 0. << '\n';
      74             : 
      75             :       // Write the elements -- 1-based!
      76           0 :       unsigned int e = 0;
      77           0 :       for (const auto & elem : the_mesh.active_element_ptr_range())
      78             :         {
      79             :           // .fro likes TRI3's
      80           0 :           libmesh_error_msg_if(elem->type() != TRI3,
      81             :                                "ERROR:  .fro format only valid for triangles!\n"
      82             :                                << "  writing of " << fname << " aborted.");
      83             : 
      84           0 :           out_stream << ++e << " \t";
      85             : 
      86           0 :           for (const Node & node : elem->node_ref_range())
      87           0 :             out_stream << node.id()+1 << " \t";
      88             : 
      89             :           //   // LHS -> RHS Mapping, for inverted triangles
      90             :           //   out_stream << elem->node_id(0)+1 << " \t";
      91             :           //   out_stream << elem->node_id(2)+1 << " \t";
      92             :           //   out_stream << elem->node_id(1)+1 << " \t";
      93             : 
      94           0 :           out_stream << "1\n";
      95           0 :         }
      96             : 
      97             :       // Write BCs.
      98             :       {
      99             :         const std::set<boundary_id_type> & bc_ids =
     100           0 :           the_mesh.get_boundary_info().get_boundary_ids();
     101             : 
     102             :         // Build a list of (elem, side, bc) tuples.
     103           0 :         auto bc_triples = the_mesh.get_boundary_info().build_side_list();
     104             : 
     105             :         // Map the boundary ids into [1,n_bc_ids],
     106             :         // treat them one at a time.
     107           0 :         boundary_id_type bc_id=0;
     108           0 :         for (const auto & id : bc_ids)
     109             :           {
     110           0 :             std::deque<dof_id_type> node_list;
     111             : 
     112             :             std::map<dof_id_type, dof_id_type>
     113           0 :               forward_edges, backward_edges;
     114             : 
     115             :             // Get all sides on this element with the relevant BC id.
     116           0 :             for (const auto & t : bc_triples)
     117           0 :               if (std::get<2>(t) == id)
     118             :                 {
     119             :                   // need to build up node_list as a sorted array of edge nodes...
     120             :                   // for the following:
     121             :                   // a---b---c---d---e
     122             :                   // node_list [ a b c d e];
     123             :                   //
     124             :                   // the issue is just how to get this out of the elem/side based data structure.
     125             :                   // the approach is to build up 'chain links' like this:
     126             :                   // a---b b---c c---d d---e
     127             :                   // and piece them together.
     128             :                   //
     129             :                   // so, for an arbitrary edge n0---n1, we build the
     130             :                   // "forward_edges"  map n0-->n1
     131             :                   // "backward_edges" map n1-->n0
     132             :                   // and then start with one chain link, and add on...
     133             :                   //
     134             :                   std::unique_ptr<const Elem> side =
     135           0 :                     the_mesh.elem_ref(std::get<0>(t)).build_side_ptr(std::get<1>(t));
     136             : 
     137             :                   const dof_id_type
     138           0 :                     n0 = side->node_id(0),
     139           0 :                     n1 = side->node_id(1);
     140             : 
     141             :                   // insert into forward-edge set
     142           0 :                   forward_edges.emplace(n0, n1);
     143             : 
     144             :                   // insert into backward-edge set
     145           0 :                   backward_edges.emplace(n1, n0);
     146             : 
     147             :                   // go ahead and add one edge to the list -- this will give us the beginning of a
     148             :                   // chain to work from!
     149           0 :                   if (node_list.empty())
     150             :                     {
     151           0 :                       node_list.push_front(n0);
     152           0 :                       node_list.push_back (n1);
     153             :                     }
     154           0 :                 }
     155             : 
     156             :             // we now have the node_list with one edge, the forward_edges, and the backward_edges
     157             :             // the node_list will be filled when (node_list.size() == (n_edges+1))
     158             :             // until that is the case simply add on to the beginning and end of the node_list,
     159             :             // building up a chain of ordered nodes...
     160           0 :             const std::size_t n_edges = forward_edges.size();
     161             : 
     162           0 :             while (node_list.size() != (n_edges+1))
     163             :               {
     164             :                 const dof_id_type
     165           0 :                   front_node = node_list.front(),
     166           0 :                   back_node  = node_list.back();
     167             : 
     168             :                 // look for front_pair in the backward_edges list
     169           0 :                 if (const auto pos = backward_edges.find(front_node);
     170           0 :                     pos != backward_edges.end())
     171             :                   {
     172           0 :                     node_list.push_front(pos->second);
     173           0 :                     backward_edges.erase(pos);
     174             :                   }
     175             : 
     176             :                 // look for back_pair in the forward_edges list
     177           0 :                 if (const auto pos = forward_edges.find(back_node);
     178           0 :                     pos != forward_edges.end())
     179             :                   {
     180           0 :                     node_list.push_back(pos->second);
     181           0 :                     forward_edges.erase(pos);
     182             :                   }
     183             :               }
     184             : 
     185           0 :             out_stream << ++bc_id << " " << node_list.size() << '\n';
     186             : 
     187           0 :             for (const auto & node_id : node_list)
     188           0 :               out_stream << node_id + 1 << " \t0\n";
     189             :           }
     190             :       }
     191           0 :     }
     192           0 : }
     193             : 
     194             : } // namespace libMesh