Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryMeshWriter.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 "ImmersedBoundaryMeshWriter.hpp"
37
38#include "ImmersedBoundaryArray.hpp"
41#include "Version.hpp"
42
46template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
56
58// Implementation
60
61template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
63 const std::string& rBaseName,
64 const bool clearOutputDir)
65 : AbstractMeshWriter<ELEMENT_DIM, SPACE_DIM>(rDirectory, rBaseName, clearOutputDir),
66 mpMesh(nullptr),
67 mpIters(new MeshWriterIterators<ELEMENT_DIM, SPACE_DIM>)
68{
69 mpIters->pNodeIter = nullptr;
70 mpIters->pElemIter = nullptr;
71 mpIters->pLamIter = nullptr;
72
73 switch (SPACE_DIM)
74 {
75 case 2:
76 {
77 geom_point corner_0(0.0, 0.0);
78 geom_point corner_1(1.0, 0.0);
79 geom_point corner_2(0.0, 1.0);
80 geom_point corner_3(1.0, 1.0);
81
82 mBoundaryEdges[0] = geom_segment(corner_0, corner_2); // left edge
83 mBoundaryEdges[1] = geom_segment(corner_1, corner_3); // right edge
84 mBoundaryEdges[2] = geom_segment(corner_0, corner_1); // bottom edge
85 mBoundaryEdges[3] = geom_segment(corner_2, corner_3); // top edge
86
87 break;
88 }
89
90 default:
92 }
93
94#ifdef CHASTE_VTK
95 // Dubious, since we shouldn't yet know what any details of the mesh are.
96 mpVtkUnstructedMesh = vtkUnstructuredGrid::New();
97#endif //CHASTE_VTK
98}
99
100template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
102{
103 if (mpIters->pNodeIter)
104 {
105 delete mpIters->pNodeIter;
106 delete mpIters->pElemIter;
107 delete mpIters->pLamIter;
108 }
109
110 delete mpIters;
111
112#ifdef CHASTE_VTK
113// Dubious, since we shouldn't yet know what any details of the mesh are.
114 mpVtkUnstructedMesh->Delete(); // Reference counted
115#endif //CHASTE_VTK
116}
117
118template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
120{
121 if (mpMesh)
122 {
123 // Sanity check
124 assert(this->mNumNodes == mpMesh->GetNumNodes());
125
126 std::vector<double> coordinates(SPACE_DIM+1);
127
128 // Get the node coordinates using the node iterator (thus skipping deleted nodes)
129 for (unsigned j=0; j<SPACE_DIM; j++)
130 {
131 coordinates[j] = (*(mpIters->pNodeIter))->GetPoint()[j];
132 }
133 coordinates[SPACE_DIM] = (*(mpIters->pNodeIter))->IsBoundaryNode();
134
135 ++(*(mpIters->pNodeIter));
136
137 return coordinates;
138 }
139 else
140 {
141 return AbstractMeshWriter<ELEMENT_DIM,SPACE_DIM>::GetNextNode(); //LCOV_EXCL_LINE fairly sure this is unreachable
142 }
143}
144
145template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147{
148 // Can only be constructed in 2D currently
149 //LCOV_EXCL_START
150 if constexpr (SPACE_DIM != 2) {
151 EXCEPTION("ImmersedBoundaryMeshWriter::GetNextImmersedBoundaryElement is not yet implemetned in 3D!");
152 }
153 //LCOV_EXCL_STOP
154
155 assert(this->mNumElements == mpMesh->GetNumElements());
156
158 elem_data.NodeIndices.resize((*(mpIters->pElemIter))->GetNumNodes());
159 for (unsigned j=0; j<elem_data.NodeIndices.size(); j++)
160 {
161 elem_data.NodeIndices[j] = (*(mpIters->pElemIter))->GetNodeGlobalIndex(j);
162 }
163
164 // Set attribute
165 elem_data.AttributeValue = (*(mpIters->pElemIter))->GetAttribute();
166
167 // Immersed boundary specific data
168 auto& element = *(*mpIters->pElemIter);
169
170 if (element.GetFluidSource() != nullptr) {
171 elem_data.hasFluidSource = true;
172 elem_data.fluidSourceIndex = element.GetFluidSource()->GetIndex();
173 } else {
174 elem_data.hasFluidSource = false;
175 }
176
177 for (auto cornerNode : element.rGetCornerNodes()) {
178 elem_data.cornerNodeIndices.push_back(cornerNode->GetIndex());
179 }
180
181 elem_data.averageNodeSpacing = element.GetAverageNodeSpacing();
182
183 elem_data.isBoundaryElement = element.IsElementOnBoundary();
184
185 ++(*(mpIters->pElemIter));
186
187 return elem_data;
188} // LCOV_EXCL_LINE
189
190template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
192{
193 //LCOV_EXCL_START
194 if constexpr (SPACE_DIM != 2) {
195 EXCEPTION("ImmersedBoundaryMeshWriter::GetNextImmersedBoundaryLamina is not yet implemetned in 3D!");
196 }
197 //LCOV_EXCL_STOP
198
199 assert(mNumLaminas == mpMesh->GetNumLaminas());
200
201 ImmersedBoundaryElementData lamina_data;
202 lamina_data.NodeIndices.resize((*(mpIters->pLamIter))->GetNumNodes());
203 for (unsigned j=0; j<lamina_data.NodeIndices.size(); j++)
204 {
205 lamina_data.NodeIndices[j] = (*(mpIters->pLamIter))->GetNodeGlobalIndex(j);
206 }
207
208 // Set attribute
209 lamina_data.AttributeValue = (*(mpIters->pLamIter))->GetAttribute();
210
211 // Immersed boundary specific data
212 auto& lamina = *(*mpIters->pLamIter);
213
214 try {
215 lamina.GetFluidSource();
216 // Can't really support 2D laminas at the moment
217 //LCOV_EXCL_START
218 if (lamina.GetFluidSource() != nullptr) {
219 lamina_data.hasFluidSource = true;
220 lamina_data.fluidSourceIndex = lamina.GetFluidSource()->GetIndex();
221 } else {
222 lamina_data.hasFluidSource = false;
223 }
224 //LCOV_EXCL_STOP
225 } catch(Exception& e) {
226 lamina_data.hasFluidSource = false;
227 }
228
229 for (auto cornerNode : lamina.rGetCornerNodes()) {
230 lamina_data.cornerNodeIndices.push_back(cornerNode->GetIndex()); // LCOV_EXCL_LINE
231 }
232
233 lamina_data.averageNodeSpacing = lamina.GetAverageNodeSpacing();
234
235 lamina_data.isBoundaryElement = lamina.IsElementOnBoundary();
236
237 ++(*(mpIters->pLamIter));
238
239 return lamina_data;
240} // LCOV_EXCL_LINE
241
242template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
244{
245#ifdef CHASTE_VTK
246 if constexpr (SPACE_DIM == 2)
247 {
248 // Create VTK mesh
249 MakeVtkMesh(rMesh);
250 // Now write VTK mesh to file
251 //assert(mpVtkUnstructedMesh->CheckAttributes() == 0);
252 vtkXMLUnstructuredGridWriter* p_writer = vtkXMLUnstructuredGridWriter::New();
253#if VTK_MAJOR_VERSION >= 6
254 p_writer->SetInputData(mpVtkUnstructedMesh);
255#else
256 p_writer->SetInput(mpVtkUnstructedMesh);
257#endif
258 // Uninitialised stuff arises (see #1079), but you can remove valgrind problems by removing compression:
259 // **** REMOVE WITH CAUTION *****
260 p_writer->SetCompressor(NULL);
261 // **** REMOVE WITH CAUTION *****
262
263 std::string vtk_file_name = this->mpOutputFileHandler->GetOutputDirectoryFullPath() + this->mBaseName;
264 if (stamp != "")
265 {
266 vtk_file_name += "_" + stamp;
267 }
268 vtk_file_name += ".vtu";
269
270 p_writer->SetFileName(vtk_file_name.c_str());
271 p_writer->Write();
272 p_writer->Delete(); // Reference counted
273 }
274 else
275 {
277 }
278#endif //CHASTE_VTK
279}
280
281template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
283{
284#ifdef CHASTE_VTK
298 // Assert mesh is 2D to simplify vtk creation
299 if constexpr (SPACE_DIM == 2)
300 {
301 FindElementOverlaps(rMesh);
302
303 // Make the Vtk mesh
304 vtkPoints* p_pts = vtkPoints::New(VTK_DOUBLE);
305 p_pts->GetData()->SetName("Vertex positions");
306
307 // Keep a track of the number of extra points that have been added to the purpose of visualisation
308 unsigned num_pts_added = 0;
309
310 // Next, we decide how to output the VTK data for elements depending on the type of overlap
311 for (auto iter = rMesh.GetElementIteratorBegin(); iter != rMesh.GetElementIteratorEnd(); ++iter)
312 {
313 unsigned elem_idx = iter->GetIndex();
314 unsigned num_nodes = iter->GetNumNodes();
315
316 // Case 1: no overlap
317 if (mElementParts.size() == 0 || mElementParts[elem_idx].empty())
318 {
319 vtkCell* p_cell = vtkPolygon::New();
320 vtkIdList* p_cell_id_list = p_cell->GetPointIds();
321 p_cell_id_list->SetNumberOfIds(num_nodes);
322
323 for (unsigned node_local_idx = 0; node_local_idx < num_nodes; ++node_local_idx)
324 {
325 // Get node, index and location
326 unsigned global_idx = iter->GetNode(node_local_idx)->GetIndex();
327 const auto& r_location = iter->GetNode(node_local_idx)->rGetLocation();
328
329 p_pts->InsertPoint(global_idx, r_location[0], r_location[1], 0.0);
330 p_cell_id_list->SetId(node_local_idx, global_idx);
331 }
332
333 mpVtkUnstructedMesh->InsertNextCell(p_cell->GetCellType(), p_cell_id_list);
334 p_cell->Delete(); // Reference counted
335 }
336 else // Case 2: at least one overlap
337 {
338 // Get the node index at the start of each part
339 const std::vector<unsigned>& start_idx_each_part = mElementParts[iter->GetIndex()];
340
341 // Get the number of parts, and put the first index onto the back of the vector for convenience
342 const auto num_parts = start_idx_each_part.size();
343
344 for (unsigned part = 0; part < num_parts; ++part)
345 {
346 const long this_start = start_idx_each_part[part];
347 const long next_start = start_idx_each_part[AdvanceMod(part, 1, num_parts)];
348
349 const long num_nodes_this_part = next_start > this_start ? next_start - this_start : num_nodes + next_start - this_start;
350
351 // Identify the extra points that need to be added
352 std::vector<c_vector<double, SPACE_DIM>> extra_locations;
353
354 // Start of the part
355 {
356 const auto& r_start_pos = iter->GetNode(this_start)->rGetLocation();
357 const auto& r_end_pos = iter->GetNode(AdvanceMod(this_start, -1, num_nodes))->rGetLocation();
358
359 const c_vector<double, SPACE_DIM> vec_a2b = rMesh.GetVectorFromAtoB(r_start_pos, r_end_pos);
360 extra_locations.emplace_back(GetIntersectionOfEdgeWithBoundary(r_start_pos, r_start_pos + vec_a2b));
361 }
362
363 // End of the part
364 {
365 const auto& r_start_pos = iter->GetNode(AdvanceMod(next_start, -1, num_nodes))->rGetLocation();
366 const auto& r_end_pos = iter->GetNode(next_start)->rGetLocation();
367
368 const c_vector<double, SPACE_DIM> vec_a2b = rMesh.GetVectorFromAtoB(r_start_pos, r_end_pos);
369 extra_locations.emplace_back(GetIntersectionOfEdgeWithBoundary(r_start_pos, r_start_pos + vec_a2b));
370 }
371
372 // If the two additional points are not on the same edge of the boundary, we also need to add a corner
373 if (std::fabs(extra_locations.front()[0] - extra_locations.back()[0]) > DBL_EPSILON &&
374 std::fabs(extra_locations.front()[1] - extra_locations.back()[1]) > DBL_EPSILON)
375 {
376 extra_locations.emplace_back(GetNearestCorner(extra_locations.front(), extra_locations.back()));
377 }
378
379 vtkCell* p_cell = vtkPolygon::New();
380 vtkIdList* p_cell_id_list = p_cell->GetPointIds();
381 p_cell_id_list->SetNumberOfIds(num_nodes_this_part + extra_locations.size());
382
383 for (long idx = 0; idx < num_nodes_this_part; ++idx)
384 {
385 unsigned node_idx = AdvanceMod(idx, this_start, num_nodes);
386
387 // Get index and location
388 unsigned global_idx = iter->GetNode(node_idx)->GetIndex();
389 const auto& r_location = iter->GetNode(node_idx)->rGetLocation();
390
391 p_pts->InsertPoint(global_idx, r_location[0], r_location[1], 0.0);
392 p_cell_id_list->SetId(idx, global_idx);
393 }
394
395 // Now add the extra locations
396 if (extra_locations.size() == 2)
397 {
398 const unsigned global_index = rMesh.GetNumNodes() + num_pts_added;
399
400 p_pts->InsertPoint(global_index, extra_locations[1][0], extra_locations[1][1], 0.0); // end
401 p_pts->InsertPoint(global_index + 1, extra_locations[0][0], extra_locations[0][1], 0.0); // start
402
403 p_cell_id_list->SetId(num_nodes_this_part, global_index);
404 p_cell_id_list->SetId(num_nodes_this_part + 1, global_index + 1);
405
406 num_pts_added += 2;
407 }
408 else if (extra_locations.size() == 3)
409 {
410 const unsigned global_index = rMesh.GetNumNodes() + num_pts_added;
411
412 p_pts->InsertPoint(global_index, extra_locations[1][0], extra_locations[1][1], 0.0); // end
413 p_pts->InsertPoint(global_index + 1, extra_locations[2][0], extra_locations[2][1], 0.0); // corner
414 p_pts->InsertPoint(global_index + 2, extra_locations[0][0], extra_locations[0][1], 0.0); // start
415
416 p_cell_id_list->SetId(num_nodes_this_part, global_index);
417 p_cell_id_list->SetId(num_nodes_this_part + 1, global_index + 1);
418 p_cell_id_list->SetId(num_nodes_this_part + 2, global_index + 2);
419
420 num_pts_added += 3;
421 }
422 else
423 {
425 }
426
427 mpVtkUnstructedMesh->InsertNextCell(p_cell->GetCellType(), p_cell_id_list);
428 p_cell->Delete(); // Reference counted
429 }
430 }
431 }
432
433 // Finally, output the VTK data for laminas
435 iter != rMesh.GetLaminaIteratorEnd();
436 ++iter)
437 {
438 unsigned num_nodes = iter->GetNumNodes();
439
440 vtkCell* p_cell = vtkPolygon::New();
441 vtkIdList* p_cell_id_list = p_cell->GetPointIds();
442 p_cell_id_list->SetNumberOfIds(num_nodes);
443
444 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
445 {
446 // Get node, index and location
447 Node<SPACE_DIM>* p_node = iter->GetNode(node_idx);
448 unsigned global_idx = p_node->GetIndex();
449
450 c_vector<double, SPACE_DIM> position = p_node->rGetLocation();
451 p_pts->InsertPoint(global_idx, position[0], position[1], 0.0);
452
453 p_cell_id_list->SetId(node_idx, global_idx);
454 }
455
456 mpVtkUnstructedMesh->InsertNextCell(3, p_cell_id_list);
457 p_cell->Delete(); // Reference counted
458 }
459
460 mpVtkUnstructedMesh->SetPoints(p_pts);
461 p_pts->Delete(); // Reference counted
462 }
463 else
464 {
466 }
467#endif //CHASTE_VTK
468}
469
470//LCOV_EXCL_START
475template <>
479//LCOV_EXCL_STOP
480
481template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
482void ImmersedBoundaryMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddCellData(std::string dataName, std::vector<double> dataPayload)
483{
484#ifdef CHASTE_VTK
485 vtkDoubleArray* p_scalars = vtkDoubleArray::New();
486 p_scalars->SetName(dataName.c_str());
487 for (unsigned i=0; i<dataPayload.size(); i++)
488 {
489 p_scalars->InsertNextValue(dataPayload[i]);
490 }
491
492 vtkCellData* p_cell_data = mpVtkUnstructedMesh->GetCellData();
493 p_cell_data->AddArray(p_scalars);
494 p_scalars->Delete(); // Reference counted
495#endif //CHASTE_VTK
496}
497
498template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
499void ImmersedBoundaryMeshWriter<ELEMENT_DIM, SPACE_DIM>::AddPointData(std::string dataName, std::vector<double> dataPayload)
500{
501#ifdef CHASTE_VTK
502 vtkDoubleArray* p_scalars = vtkDoubleArray::New();
503 p_scalars->SetName(dataName.c_str());
504 for (double scalar : dataPayload)
505 {
506 p_scalars->InsertNextValue(scalar);
507 }
508
509 vtkPointData* p_point_data = mpVtkUnstructedMesh->GetPointData();
510 p_point_data->AddArray(p_scalars);
511 p_scalars->Delete(); // Reference counted
512#endif //CHASTE_VTK
513}
514
515template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
517{
518 this->mpMeshReader = nullptr;
519 mpMesh = &rMesh;
520
521 this->mNumNodes = mpMesh->GetNumNodes();
522 this->mNumElements = mpMesh->GetNumElements();
523 mNumLaminas = mpMesh->GetNumLaminas();
524
525 typedef typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator NodeIterType;
526 mpIters->pNodeIter = new NodeIterType(mpMesh->GetNodeIteratorBegin());
527
529 mpIters->pElemIter = new ElemIterType(mpMesh->GetElementIteratorBegin());
530
532 mpIters->pLamIter = new LamIterType(mpMesh->GetLaminaIteratorBegin());
533
534 WriteFiles();
535}
536
537template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
539{
540 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
541
542 // Write node file
543 std::string node_file_name = this->mBaseName + ".node";
544 out_stream p_node_file = this->mpOutputFileHandler->OpenOutputFile(node_file_name);
545
546 // Write the node header
547 unsigned num_attr = 0;
548 unsigned max_bdy_marker = 1; // as we include boundary node information in the node file
549 unsigned num_nodes = this->GetNumNodes();
550
551 *p_node_file << num_nodes << "\t";
552 *p_node_file << SPACE_DIM << "\t";
553 *p_node_file << num_attr << "\t";
554 *p_node_file << mpMesh->GetCharacteristicNodeSpacing() << "\t";
555 *p_node_file << max_bdy_marker << "\n";
556 *p_node_file << std::setprecision(6);
557
558 // Write each node's data
559 for (unsigned item_num=0; item_num<num_nodes; item_num++)
560 {
561 std::vector<double> current_item = this->GetNextNode();
562 *p_node_file << item_num;
563 for (unsigned i=0; i<SPACE_DIM+1; i++)
564 {
565 *p_node_file << "\t" << current_item[i];
566 }
567 *p_node_file << "\n";
568 }
569 *p_node_file << comment << "\n";
570 p_node_file->close();
571
572 // Write element file
573 std::string element_file_name = this->mBaseName + ".elem";
574 out_stream p_element_file = this->mpOutputFileHandler->OpenOutputFile(element_file_name);
575
576 // Write the element header
577 num_attr = 1; //Always write element attributes
578 unsigned num_elements = this->GetNumElements();
579 *p_element_file << num_elements << "\t" << num_attr << "\n";
580
581 // Write each element's data
582 for (unsigned item_num=0; item_num<num_elements; item_num++)
583 {
584 if (SPACE_DIM == 2) // In 2D, write the node indices owned by this element
585 {
586 // Get data for this element
587 ImmersedBoundaryElementData elem_data = this->GetNextImmersedBoundaryElement();
588
589 // Get the node indices owned by this element
590 std::vector<unsigned> node_indices = elem_data.NodeIndices;
591
592 // Write this element's index and the number of nodes owned by it to file
593 *p_element_file << item_num << "\t" << node_indices.size();
594
595 // Write the node indices owned by this element to file
596 for (unsigned node_index : node_indices)
597 {
598 *p_element_file << "\t" << node_index;
599 }
600
601 *p_element_file << "\t" << elem_data.AttributeValue;
602
603 // Has fluid source
604 *p_element_file << "\t" << elem_data.hasFluidSource;
605
606 // Fluid source index
607 if (elem_data.hasFluidSource) {
608 *p_element_file << "\t" << elem_data.fluidSourceIndex;
609 }
610
611 // Corner node indices
612 *p_element_file << "\t" << elem_data.cornerNodeIndices.size();
613
614 for (auto nodeIndex : elem_data.cornerNodeIndices) {
615 *p_element_file << "\t" << nodeIndex;
616 }
617
618 // Average node spacing
619 *p_element_file << "\t" << elem_data.averageNodeSpacing;
620
621 // Is boundary element
622 *p_element_file << "\t" << elem_data.isBoundaryElement;
623
624 // New line
625 *p_element_file << "\n";
626 }
627 else // 3D
628 {
629 }
630 }
631 *p_element_file << comment << "\n";
632 p_element_file->close();
633
634 // Write lamina file
635 std::string lamina_file_name = this->mBaseName + ".lam";
636 out_stream p_lamina_file = this->mpOutputFileHandler->OpenOutputFile(lamina_file_name);
637
638 // Write the lamina header
639 num_attr = 1; //Always write element attributes
640 unsigned num_laminas = mpMesh->GetNumLaminas();
641 *p_lamina_file << num_laminas << "\t" << num_attr << "\n";
642
643 // Write each lamina's data
644 for (unsigned item_num=0; item_num<num_laminas; item_num++)
645 {
646 if (SPACE_DIM == 2) // In 2D, write the node indices owned by this element
647 {
648 // Get data for this element
649 ImmersedBoundaryElementData lamina_data = this->GetNextImmersedBoundaryLamina();
650
651 // Get the node indices owned by this element
652 std::vector<unsigned> node_indices = lamina_data.NodeIndices;
653
654 // Write this element's index and the number of nodes owned by it to file
655 *p_lamina_file << item_num << "\t" << node_indices.size();
656
657 // Write the node indices owned by this element to file
658 for (unsigned i=0; i<node_indices.size(); i++)
659 {
660 *p_lamina_file << "\t" << node_indices[i];
661 }
662
663 *p_lamina_file << "\t" << lamina_data.AttributeValue;
664
665 // Has fluid source
666 *p_lamina_file << "\t" << lamina_data.hasFluidSource;
667
668 // Fluid source index
669 if (lamina_data.hasFluidSource) {
670 *p_lamina_file << "\t" << lamina_data.fluidSourceIndex; //LCOV_EXCL_LINE
671 }
672
673 // Corner node indices
674 *p_lamina_file << "\t" << lamina_data.cornerNodeIndices.size();
675
676 for (auto nodeIndex : lamina_data.cornerNodeIndices) {
677 *p_lamina_file << "\t" << nodeIndex; //LCOV_EXCL_LINE
678 }
679
680 // Average node spacing
681 *p_lamina_file << "\t" << lamina_data.averageNodeSpacing;
682
683 // Is boundary lamina
684 *p_lamina_file << "\t" << lamina_data.isBoundaryElement;
685
686
687 // New line
688 *p_lamina_file << "\n";
689 }
690 else // 3D
691 {
692 }
693 }
694 *p_lamina_file << comment << "\n";
695 p_lamina_file->close();
696
697 // Write grid file
698 std::string grid_file_name = this->mBaseName + ".grid";
699 out_stream p_grid_file = this->mpOutputFileHandler->OpenOutputFile(grid_file_name);
700
701 // Write the element header
702 unsigned num_gridpts_x = mpMesh->GetNumGridPtsX();
703 unsigned num_gridpts_y = mpMesh->GetNumGridPtsY();
704
705 *p_grid_file << num_gridpts_x << "\t" << num_gridpts_y << "\n";
706
707 // Write grid data
708 const multi_array<double, 3>& vel_grids = mpMesh->rGet2dVelocityGrids();
709
710 for (unsigned y_idx = 0; y_idx < num_gridpts_y; y_idx ++)
711 {
712 for (unsigned x_idx = 0; x_idx < num_gridpts_x; x_idx ++)
713 {
714 *p_grid_file << vel_grids[0][x_idx][y_idx] << "\t";
715 }
716 *p_grid_file << "\n";
717 }
718
719 for (unsigned y_idx = 0; y_idx < num_gridpts_y; y_idx ++)
720 {
721 for (unsigned x_idx = 0; x_idx < num_gridpts_x; x_idx ++)
722 {
723 *p_grid_file << vel_grids[1][x_idx][y_idx] << "\t";
724 }
725 *p_grid_file << "\n";
726 }
727
728 *p_grid_file << comment << "\n";
729 p_grid_file->close();
730}
731
732template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
734{
735 if constexpr (SPACE_DIM == 2)
736 {
737 // Resize and initialise the vector of overlaps (bools; whether each element has an overlap)
738 mElementParts.resize(rMesh.GetNumAllElements());
739 for (auto& parts : mElementParts)
740 {
741 parts.clear();
742 }
743
744 // We loop over each element and the node index at each discontinuity due to periodic boundaries
745 for (auto iter = rMesh.GetElementIteratorBegin(); iter != rMesh.GetElementIteratorEnd(); ++iter)
746 {
747 for (unsigned node_idx = 0; node_idx < iter->GetNumNodes(); ++node_idx)
748 {
749 const unsigned prev_idx = AdvanceMod(node_idx, -1, iter->GetNumNodes());
750
751 const auto& this_location = iter->GetNode(node_idx)->rGetLocation();
752 const auto& prev_location = iter->GetNode(prev_idx)->rGetLocation();
753
754 if (norm_inf(this_location - prev_location) > 0.5)
755 {
756 mElementParts[iter->GetIndex()].emplace_back(node_idx);
757 }
758 }
759 }
760 }
761 else
762 {
764 }
765}
766
767template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
769 const c_vector<double, SPACE_DIM>& rStart,
770 const c_vector<double, SPACE_DIM>& rEnd)
771{
772 // Note that this function relies on SPACE_DIM == 2
773
774 // Turn the c_vector start and end into a boost::geometry segment object
775 geom_point start(rStart[0], rStart[1]);
776 geom_point end(rEnd[0], rEnd[1]);
777 geom_segment edge(start, end);
778
779 // Identify which boundary edge the intersection is with
780 std::vector<geom_point> intersections;
781 for (const auto& boundary_edge : mBoundaryEdges)
782 {
783 if (boost::geometry::intersects(edge, boundary_edge))
784 {
785 boost::geometry::intersection(edge, boundary_edge, intersections);
786 break;
787 }
788 }
789
790 // There should be exactly one intersection
791 if (intersections.size() != 1)
792 {
794 }
795
796 return Create_c_vector(intersections.front().get<0>(), intersections.front().get<1>());
797}
798
799template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
801 const c_vector<double, SPACE_DIM>& rA, const c_vector<double, SPACE_DIM>& rB) const
802{
803 // Identify the nearest corner to the average location of the two vectors
804 const double x = 0.5 * (rA[0] + rB[0]) < 0.5 ? 0.0 : 1.0;
805 const double y = 0.5 * (rA[1] + rB[1]) < 0.5 ? 0.0 : 1.0;
806
807 return Create_c_vector(x, y);
808}
809
810template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
811const std::vector<std::vector<unsigned>>& ImmersedBoundaryMeshWriter<ELEMENT_DIM, SPACE_DIM>::rGetElementParts() const
812{
813 return mElementParts;
814}
815
816// Explicit instantiation
#define EXCEPTION(message)
#define NEVER_REACHED
virtual std::vector< double > GetNextNode()
static std::string GetProvenanceString()
ImmersedBoundaryMeshWriter(const std::string &rDirectory, const std::string &rBaseName, bool clearOutputDir=true)
boost::geometry::model::point< double, 2, boost::geometry::cs::cartesian > geom_point
void FindElementOverlaps(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
ImmersedBoundaryElementData GetNextImmersedBoundaryElement()
c_vector< double, SPACE_DIM > GetNearestCorner(const c_vector< double, SPACE_DIM > &rA, const c_vector< double, SPACE_DIM > &rB) const
ImmersedBoundaryElementData GetNextImmersedBoundaryLamina()
const std::vector< std::vector< unsigned > > & rGetElementParts() const
void WriteVtkUsingMesh(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, std::string stamp="")
void AddPointData(std::string dataName, std::vector< double > dataPayload)
void MakeVtkMesh(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
std::array< geom_segment, 4 > mBoundaryEdges
void WriteFilesUsingMesh(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
boost::geometry::model::segment< geom_point > geom_segment
c_vector< double, SPACE_DIM > GetIntersectionOfEdgeWithBoundary(const c_vector< double, SPACE_DIM > &rStart, const c_vector< double, SPACE_DIM > &rEnd)
void AddCellData(std::string dataName, std::vector< double > dataPayload)
MeshWriterIterators< ELEMENT_DIM, SPACE_DIM > * mpIters
ImmersedBoundaryElementIterator GetElementIteratorEnd()
ImmersedBoundaryLaminaIterator GetLaminaIteratorBegin(bool skipDeletedLaminas=true)
ImmersedBoundaryElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
ImmersedBoundaryLaminaIterator GetLaminaIteratorEnd()
virtual unsigned GetNumNodes() const
c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocation1, const c_vector< double, SPACE_DIM > &rLocation2)
unsigned GetNumAllElements() const
Definition Node.hpp:59
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
unsigned GetIndex() const
Definition Node.cpp:158
ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM >::ImmersedBoundaryElementIterator * pElemIter
ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM >::ImmersedBoundaryLaminaIterator * pLamIter
AbstractMesh< ELEMENT_DIM, SPACE_DIM >::NodeIterator * pNodeIter