GModelIO_P3D.cpp

Go to the documentation of this file.

// Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file in the Gmsh root directory for license information.
// Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
 
#include "GModel.h"
#include "OS.h"
#include "MQuadrangle.h"
#include "MHexahedron.h"
#include "discreteFace.h"
#include "discreteRegion.h"
 
int GModel::readP3D(const std::string &name)
{
  FILE *fp = Fopen(name.c_str(), "r");
  if(!fp) {
    Msg::Error("Unable to open file '%s'", name.c_str());
    return 0;
  }
 
  int numBlocks = 0;
  if(fscanf(fp, "%d", &numBlocks) != 1 || numBlocks <= 0) {
    fclose(fp);
    return 0;
  }
 
  std::vector<int> Ni(numBlocks), Nj(numBlocks), Nk(numBlocks);
  for(int n = 0; n < numBlocks; n++)
    if(fscanf(fp, "%d %d %d", &Ni[n], &Nj[n], &Nk[n]) != 3) {
      fclose(fp);
      return 0;
    }
 
  for(int n = 0; n < numBlocks; n++) {
    if(Nk[n] == 1) {
      GFace *gf = new discreteFace(this, getMaxElementaryNumber(2) + 1);
      add(gf);
      gf->transfinite_vertices.resize(Ni[n]);
      for(int i = 0; i < Ni[n]; i++) gf->transfinite_vertices[i].resize(Nj[n]);
      for(int coord = 0; coord < 3; coord++) {
        for(int j = 0; j < Nj[n]; j++) {
          for(int i = 0; i < Ni[n]; i++) {
            double d;
            if(fscanf(fp, "%lf", &d) != 1) {
              fclose(fp);
              return 0;
            }
            if(coord == 0) {
              MVertex *v = new MVertex(d, 0., 0., gf);
              gf->transfinite_vertices[i][j] = v;
              gf->mesh_vertices.push_back(v);
            }
            else if(coord == 1) {
              gf->transfinite_vertices[i][j]->y() = d;
            }
            else if(coord == 2) {
              gf->transfinite_vertices[i][j]->z() = d;
            }
          }
        }
      }
      for(std::size_t i = 0; i < gf->transfinite_vertices.size() - 1; i++)
        for(std::size_t j = 0; j < gf->transfinite_vertices[0].size() - 1; j++)
          gf->quadrangles.push_back(new MQuadrangle(
            gf->transfinite_vertices[i][j], gf->transfinite_vertices[i + 1][j],
            gf->transfinite_vertices[i + 1][j + 1],
            gf->transfinite_vertices[i][j + 1]));
    }
    else {
      GRegion *gr = new discreteRegion(this, getMaxElementaryNumber(3) + 1);
      add(gr);
      gr->transfinite_vertices.resize(Ni[n]);
      for(int i = 0; i < Ni[n]; i++) {
        gr->transfinite_vertices[i].resize(Nj[n]);
        for(int j = 0; j < Nj[n]; j++) {
          gr->transfinite_vertices[i][j].resize(Nk[n]);
        }
      }
      for(int coord = 0; coord < 3; coord++) {
        for(int k = 0; k < Nk[n]; k++) {
          for(int j = 0; j < Nj[n]; j++) {
            for(int i = 0; i < Ni[n]; i++) {
              double d;
              if(fscanf(fp, "%lf", &d) != 1) {
                fclose(fp);
                return 0;
              }
              if(coord == 0) {
                MVertex *v = new MVertex(d, 0., 0., gr);
                gr->transfinite_vertices[i][j][k] = v;
                gr->mesh_vertices.push_back(v);
              }
              else if(coord == 1) {
                gr->transfinite_vertices[i][j][k]->y() = d;
              }
              else if(coord == 2) {
                gr->transfinite_vertices[i][j][k]->z() = d;
              }
            }
          }
        }
      }
      for(std::size_t i = 0; i < gr->transfinite_vertices.size() - 1; i++)
        for(std::size_t j = 0; j < gr->transfinite_vertices[0].size() - 1; j++)
          for(std::size_t k = 0; k < gr->transfinite_vertices[0][0].size() - 1;
              k++)
            gr->hexahedra.push_back(
              new MHexahedron(gr->transfinite_vertices[i][j][k],
                              gr->transfinite_vertices[i + 1][j][k],
                              gr->transfinite_vertices[i + 1][j + 1][k],
                              gr->transfinite_vertices[i][j + 1][k],
                              gr->transfinite_vertices[i][j][k + 1],
                              gr->transfinite_vertices[i + 1][j][k + 1],
                              gr->transfinite_vertices[i + 1][j + 1][k + 1],
                              gr->transfinite_vertices[i][j + 1][k + 1]));
    }
  }
 
  fclose(fp);
  return 1;
}
 
int GModel::writeP3D(const std::string &name, bool saveAll,
                     double scalingFactor)
{
  FILE *fp = Fopen(name.c_str(), "w");
  if(!fp) {
    Msg::Error("Unable to open file '%s'", name.c_str());
    return 0;
  }
 
  if(noPhysicalGroups()) saveAll = true;
 
  std::vector<GFace *> faces;
  for(auto it = firstFace(); it != lastFace(); ++it)
    if((*it)->transfinite_vertices.size() &&
       (*it)->transfinite_vertices[0].size() &&
       ((*it)->physicals.size() || saveAll))
      faces.push_back(*it);
 
  std::vector<GRegion *> regions;
  for(auto it = firstRegion(); it != lastRegion(); ++it)
    if((*it)->transfinite_vertices.size() &&
       (*it)->transfinite_vertices[0].size() &&
       (*it)->transfinite_vertices[0][0].size() &&
       ((*it)->physicals.size() || saveAll))
      regions.push_back(*it);
 
  if(faces.empty() && regions.empty()) {
    Msg::Warning("No structured grids to save");
    fclose(fp);
    return 0;
  }
 
  fprintf(fp, "%d\n", (int)(faces.size() + regions.size()));
 
  for(std::size_t i = 0; i < faces.size(); i++)
    fprintf(fp, "%d %d 1\n", (int)faces[i]->transfinite_vertices.size(),
            (int)faces[i]->transfinite_vertices[0].size());
 
  for(std::size_t i = 0; i < regions.size(); i++)
    fprintf(fp, "%d %d %d\n", (int)regions[i]->transfinite_vertices.size(),
            (int)regions[i]->transfinite_vertices[0].size(),
            (int)regions[i]->transfinite_vertices[0][0].size());
 
  for(std::size_t i = 0; i < faces.size(); i++) {
    GFace *gf = faces[i];
    for(int coord = 0; coord < 3; coord++) {
      for(std::size_t k = 0; k < gf->transfinite_vertices[0].size(); k++) {
        for(std::size_t j = 0; j < gf->transfinite_vertices.size(); j++) {
          MVertex *v = gf->transfinite_vertices[j][k];
          double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
          fprintf(fp, "%.16g ", d * scalingFactor);
        }
        fprintf(fp, "\n");
      }
    }
  }
 
  for(std::size_t i = 0; i < regions.size(); i++) {
    GRegion *gr = regions[i];
    for(int coord = 0; coord < 3; coord++) {
      for(std::size_t l = 0; l < gr->transfinite_vertices[0][0].size(); l++) {
        for(std::size_t k = 0; k < gr->transfinite_vertices[0].size(); k++) {
          for(std::size_t j = 0; j < gr->transfinite_vertices.size(); j++) {
            MVertex *v = gr->transfinite_vertices[j][k][l];
            double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
            fprintf(fp, "%.16g ", d * scalingFactor);
          }
          fprintf(fp, "\n");
        }
      }
    }
  }
 
  fclose(fp);
  return 1;
}