DUNE PDELab (git)

starcdreader.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_STARCD_READER_HH
6#define DUNE_STARCD_READER_HH
7
9
10#include <dune/geometry/type.hh>
12#include <iostream>
13#include <fstream>
14#include <memory>
15
16namespace Dune {
17
51 template <class GridType>
53
54 public:
55
65 static std::unique_ptr<GridType> read(const std::string& fileName, bool verbose = true)
66 {
67 // extract the grid dimension
68 const int dim = GridType::dimension;
69
70 // currently only dim = 3 is implemented
71 if (dim != 3)
73 "Reading Star-CD format is not implemented for dimension " << dim);
74
75 // set up the grid factory
77
78 // set the name of the vertex file
79 std::string vertexFileName = fileName + ".vrt";
80
81 // set the vertex input stream
82 std::ifstream vertexFile(vertexFileName.c_str());
83 if (!vertexFile)
84 DUNE_THROW(Dune::IOError, "Could not open " << vertexFileName);
85
86 // read the vertices
87 int dummyIdx;
88 int numberOfVertices = 0;
89 while (vertexFile >> dummyIdx) {
90 numberOfVertices++;
91
93
94 for (int k = 0; k < dim; k++)
95 vertexFile >> position[k];
96
97 factory.insertVertex(position);
98 }
99 if (verbose)
100 std::cout << numberOfVertices << " vertices read." << std::endl;
101
102 // set the name of the element file
103 std::string elementFileName = fileName + ".cel";
104
105 // set the element input stream
106 std::ifstream elementFile(elementFileName.c_str());
107 if (!elementFile)
108 DUNE_THROW(Dune::IOError, "Could not open " << elementFileName);
109
110 // read the elements
111 int numberOfElements = 0;
112 int numberOfSimplices = 0;
113 int numberOfPyramids = 0;
114 int numberOfPrisms = 0;
115 int numberOfCubes = 0;;
116 int maxNumberOfVertices = (int)pow(2, dim);
117 int isVolume = 1;
118 while (elementFile >> dummyIdx) {
119 std::vector<unsigned int> vertices(maxNumberOfVertices);
120 for (int k = 0; k < maxNumberOfVertices; k++)
121 elementFile >> vertices[k];
122
123 int boundaryId;
124 elementFile >> boundaryId;
125
126 int volumeOrSurface[2];
127 elementFile >> volumeOrSurface[0] >> volumeOrSurface[1];
128
129 if (volumeOrSurface[0] == isVolume) {
130 numberOfElements++;
131
132 if (vertices[2] == vertices[3]) { // simplex or prism
133 if (vertices[4] == vertices[5]) { // simplex
134 numberOfSimplices++;
135 std::vector<unsigned int> simplexVertices(4);
136 for (int k = 0; k < 3; k++)
137 simplexVertices[k] = vertices[k] - 1;
138 simplexVertices[3] = vertices[4] - 1;
139 factory.insertElement(Dune::GeometryTypes::tetrahedron, simplexVertices);
140 }
141 else { // prism
142 numberOfPrisms++;
143 std::vector<unsigned int> prismVertices(6);
144 for (int k = 0; k < 3; k++)
145 prismVertices[k] = vertices[k] - 1;
146 for (int k = 3; k < 6; k++)
147 prismVertices[k] = vertices[k+1] - 1;
148 factory.insertElement(Dune::GeometryTypes::prism, prismVertices);
149 }
150 }
151 else { // cube or pyramid
152 if (vertices[4] == vertices[5]) { // pyramid
153 numberOfPyramids++;
154 std::vector<unsigned int> pyramidVertices(5);
155 for (int k = 0; k < 5; k++)
156 pyramidVertices[k] = vertices[k] - 1;
157 factory.insertElement(Dune::GeometryTypes::pyramid, pyramidVertices);
158 }
159 else { // cube
160 numberOfCubes++;
161 std::vector<unsigned int> cubeVertices(8);
162 for (int k = 0; k < 8; k++)
163 cubeVertices[k] = vertices[k] - 1;
164 std::swap(cubeVertices[2], cubeVertices[3]);
165 std::swap(cubeVertices[6], cubeVertices[7]);
166 factory.insertElement(Dune::GeometryTypes::hexahedron, cubeVertices);
167 }
168 }
169 }
170 }
171 if (verbose)
172 std::cout << numberOfElements << " elements read: "
173 << numberOfSimplices << " simplices, " << numberOfPyramids << " pyramids, "
174 << numberOfPrisms << " prisms, " << numberOfCubes << " cubes." << std::endl;
175
176 // finish off the construction of the grid object
177 if (verbose)
178 std::cout << "Starting createGrid() ... " << std::flush;
179
180 return factory.createGrid();
181
182 }
183
184 };
185
186}
187
188#endif
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Provide a generic factory class for unstructured grids.
Definition: gridfactory.hh:275
virtual void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices)
Insert an element into the coarse grid.
Definition: gridfactory.hh:307
virtual void insertVertex(const FieldVector< ctype, dimworld > &pos)
Insert a vertex into the coarse grid.
Definition: gridfactory.hh:296
virtual std::unique_ptr< GridType > createGrid()
Finalize grid creation and hand over the grid.
Definition: gridfactory.hh:333
Default exception class for I/O errors.
Definition: exceptions.hh:231
Default exception for dummy implementations.
Definition: exceptions.hh:263
File reader for the Star-CD format.
Definition: starcdreader.hh:52
static std::unique_ptr< GridType > read(const std::string &fileName, bool verbose=true)
Read grid from a Star-CD file.
Definition: starcdreader.hh:65
Provide a generic factory class for unstructured grids.
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType prism
GeometryType representing a 3D prism.
Definition: type.hh:528
constexpr GeometryType hexahedron
GeometryType representing a hexahedron.
Definition: type.hh:534
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:522
constexpr GeometryType tetrahedron
GeometryType representing a tetrahedron.
Definition: type.hh:516
Dune namespace.
Definition: alignedallocator.hh:13
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)