Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryMeshReader.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "ImmersedBoundaryMeshReader.hpp"
37#include "Exception.hpp"
38
39#include <sstream>
40
41template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
43 : mFilesBaseName(pathBaseName),
44 mIndexFromZero(false), // initially assume that nodes are not numbered from zero
45 mNumNodes(0),
46 mNumElements(0),
47 mNumLaminas(0),
48 mNodesRead(0),
49 mElementsRead(0),
50 mLaminasRead(0),
51 mNumElementAttributes(0)
52{
53 OpenFiles();
55}
56
57template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
59{
60 return mNumElements;
61}
62
63template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
65{
66 return mNumLaminas;
67}
68
69template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
71{
72 return mNumNodes;
73}
74
75template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
77{
78 return mNumGridPtsX;
79}
80
81template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
83{
84 return mNumGridPtsY;
85}
86
87template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
89{
90 return mNumElementAttributes;
91}
92
93template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
95{
96 return mNumLaminaAttributes;
97}
98
99template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
101{
102 return mCharacteristicNodeSpacing;
103}
104
105template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
107{
108 CloseFiles();
109 OpenFiles();
110 ReadHeaders();
111
112 mNodesRead = 0;
113 mElementsRead = 0;
114}
115
116template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
118{
119 std::vector<double> node_data;
120
121 std::string buffer = "";
122 GetNextLineFromStream(mNodesFile, buffer);
123
124 std::stringstream buffer_stream(buffer);
125
126 unsigned index;
127 buffer_stream >> index;
128
129 unsigned offset = mIndexFromZero ? 0 : 1;
130 if (index != mNodesRead + offset)
131 {
132 EXCEPTION("Data for node " << mNodesRead << " missing");
133 }
134
135 double node_value;
136 for (unsigned i=0; i<SPACE_DIM+1; i++)
137 {
138 buffer_stream >> node_value;
139 node_data.push_back(node_value);
140 }
141
142 mNodesRead++;
143 return node_data;
144}
145
146template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
148{
149 std::vector<double> grid_row;
150 grid_row.resize(mNumGridPtsX);
151
152 std::string buffer = "";
153 GetNextLineFromStream(mGridFile, buffer);
154
155 std::stringstream buffer_stream(buffer);
156
157 double value;
158 for (unsigned i=0; i<mNumGridPtsX; i++)
159 {
160 buffer_stream >> value;
161 grid_row[i] = value;
162 }
163
164 return grid_row;
165}
166
167template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
169{
170 // Create data structure for this element
171 ImmersedBoundaryElementData element_data;
172
173 std::string buffer = "";
174 GetNextLineFromStream(mElementsFile, buffer);
175
176 std::stringstream buffer_stream(buffer);
177
178 unsigned element_index;
179 buffer_stream >> element_index;
180
181 unsigned offset = mIndexFromZero ? 0 : 1;
182 if (element_index != mElementsRead + offset)
183 {
184 EXCEPTION("Data for element " << mElementsRead << " missing");
185 }
186
187 unsigned num_nodes_in_element;
188 buffer_stream >> num_nodes_in_element;
189
190 // Store node indices owned by this element
191 unsigned node_index;
192 for (unsigned i=0; i<num_nodes_in_element; i++)
193 {
194 buffer_stream >> node_index;
195 element_data.NodeIndices.push_back(node_index - offset);
196 }
197
198 if (mNumElementAttributes > 0)
199 {
200 assert(mNumElementAttributes == 1);
201
202 buffer_stream >> element_data.AttributeValue;
203 }
204 else
205 {
206 element_data.AttributeValue = 0;
207 }
208
209 // Immersed boundary specific data
210 buffer_stream >> element_data.hasFluidSource;
211
212 if (element_data.hasFluidSource) {
213 buffer_stream >> element_data.fluidSourceIndex;
214 }
215
216 unsigned num_corner_nodes;
217 buffer_stream >> num_corner_nodes;
218 for (unsigned i = 0; i < num_corner_nodes; i++) {
219 buffer_stream >> node_index;
220 element_data.cornerNodeIndices.push_back(node_index);
221 }
222
223 buffer_stream >> element_data.averageNodeSpacing;
224
225 buffer_stream >> element_data.isBoundaryElement;
226
227 mElementsRead++;
228 return element_data;
229}
230
231template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
233{
234 // Create data structure for this lamina
235 ImmersedBoundaryElementData lamina_data;
236
237 std::string buffer = "";
238 GetNextLineFromStream(mLaminasFile, buffer);
239
240 std::stringstream buffer_stream(buffer);
241
242 unsigned lamina_idx;
243 buffer_stream >> lamina_idx;
244
245 unsigned offset = mIndexFromZero ? 0 : 1;
246 if (lamina_idx != mLaminasRead + offset)
247 {
248 EXCEPTION("Data for lamina " << mLaminasRead << " missing");
249 }
250
251 unsigned num_nodes_in_lamina;
252 buffer_stream >> num_nodes_in_lamina;
253
254 // Store node indices owned by this element
255 unsigned node_index;
256 for (unsigned i=0; i<num_nodes_in_lamina; i++)
257 {
258 buffer_stream >> node_index;
259 lamina_data.NodeIndices.push_back(node_index - offset);
260 }
261
262 if (mNumLaminaAttributes > 0)
263 {
264 assert(mNumLaminaAttributes == 1);
265
266 buffer_stream >> lamina_data.AttributeValue;
267 }
268 else
269 {
270 lamina_data.AttributeValue = 0;
271 }
272
273 // Immersed boundary specific data
274 buffer_stream >> lamina_data.hasFluidSource;
275
276 if (lamina_data.hasFluidSource) {
277 buffer_stream >> lamina_data.fluidSourceIndex; //LCOV_EXCL_LINE
278 }
279
280 unsigned num_corner_nodes;
281 buffer_stream >> num_corner_nodes;
282 for (unsigned i = 0; i < num_corner_nodes; i++) {
283 buffer_stream >> node_index;
284 lamina_data.cornerNodeIndices.push_back(node_index);
285 }
286
287 buffer_stream >> lamina_data.averageNodeSpacing;
288
289 buffer_stream >> lamina_data.isBoundaryElement;
290
291 mLaminasRead++;
292 return lamina_data;
293}
294
295
296template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
298{
299 OpenNodeFile();
300 OpenElementsFile();
301 OpenGridFile();
302 OpenLaminasFile();
303}
304
305template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
307{
308 // Nodes definition
309 std::string file_name = mFilesBaseName + ".node";
310 mNodesFile.open(file_name.c_str());
311 if (!mNodesFile.is_open())
312 {
313 EXCEPTION("Could not open data file: " + file_name);
314 }
315}
316
317template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
319{
320 // Elements definition
321 std::string file_name;
322 file_name = mFilesBaseName + ".elem";
323
324 mElementsFile.open(file_name.c_str());
325 if (!mElementsFile.is_open())
326 {
327 EXCEPTION("Could not open data file: " + file_name);
328 }
329}
330
331template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
333{
334 // Elements definition
335 std::string file_name;
336 file_name = mFilesBaseName + ".lam";
337
338 mLaminasFile.open(file_name.c_str());
339 if (!mLaminasFile.is_open())
340 {
341 EXCEPTION("Could not open data file: " + file_name);
342 }
343}
344
345template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
347{
348 // Grid definition
349 std::string file_name;
350 file_name = mFilesBaseName + ".grid";
351
352 mGridFile.open(file_name.c_str());
353 if (!mGridFile.is_open())
354 {
355 EXCEPTION("Could not open data file: " + file_name);
356 }
357}
358
359template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
361{
362 std::string buffer = "";
363 GetNextLineFromStream(mNodesFile, buffer);
364
365 unsigned local_space_dim;
366
367 std::stringstream buffer_stream(buffer);
368 buffer_stream >> mNumNodes >> local_space_dim >> mNumNodeAttributes >> mCharacteristicNodeSpacing;
369
370 // Get the next line to see if nodes are indexed from zero or not
371 GetNextLineFromStream(mNodesFile, buffer);
372 std::stringstream node_buffer_stream(buffer);
373
374 unsigned first_index;
375 node_buffer_stream >> first_index;
376 assert(first_index == 0 || first_index == 1);
377 mIndexFromZero = (first_index == 0);
378
379 // Close, reopen, skip header
380 mNodesFile.close();
381 OpenNodeFile();
382 GetNextLineFromStream(mNodesFile, buffer);
383
384 // Element header
385 GetNextLineFromStream(mElementsFile, buffer);
386 std::stringstream element_buffer_stream(buffer);
387 element_buffer_stream >> mNumElements >> mNumElementAttributes;
388
389 // Lamina header
390 GetNextLineFromStream(mLaminasFile, buffer);
391 std::stringstream lamina_buffer_stream(buffer);
392 lamina_buffer_stream >> mNumLaminas >> mNumLaminaAttributes;
393
394 // Grid header
395 GetNextLineFromStream(mGridFile, buffer);
396 std::stringstream grid_buffer_stream(buffer);
397 grid_buffer_stream >> mNumGridPtsX >> mNumGridPtsY;
398}
399
400template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
402{
403 mNodesFile.close();
404 mElementsFile.close();
405 mLaminasFile.close();
406 mGridFile.close();
407}
408
409template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
410void ImmersedBoundaryMeshReader<ELEMENT_DIM, SPACE_DIM>::GetNextLineFromStream(std::ifstream& fileStream, std::string& rawLine)
411{
412 bool line_is_blank;
413
414 do
415 {
416 getline(fileStream, rawLine);
417
418 if (fileStream.eof())
419 {
420 EXCEPTION("Cannot get the next line from node, element, lamina or grid file due to incomplete data");
421 }
422
423 // Get rid of any comment
424 rawLine = rawLine.substr(0,rawLine.find('#', 0));
425
426 line_is_blank = (rawLine.find_first_not_of(" \t", 0) == std::string::npos);
427 }
428 while (line_is_blank);
429}
430
431template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
436
437template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
439{
440 ElementData ret;
441 ret.NodeIndices = std::vector<unsigned>();
442 ret.AttributeValue = 0;
443 return ret;
444}
445
446template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
448{
449 ElementData ret;
450 ret.NodeIndices = std::vector<unsigned>();
451 ret.AttributeValue = 0;
452 return ret;
453}
454
455// Explicit instantiation
#define EXCEPTION(message)
ImmersedBoundaryElementData GetNextImmersedBoundaryElementData()
ImmersedBoundaryMeshReader(std::string pathBaseName)
void GetNextLineFromStream(std::ifstream &fileStream, std::string &rawLine)
ImmersedBoundaryElementData GetNextImmersedBoundaryLaminaData()
std::vector< unsigned > NodeIndices