Skip to content
Snippets Groups Projects
Commit efc2d907 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add VTK output for polynomial meshes of degree > 2

parent 3f96e465
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/output/VTKWriter.cpp
+ 98
− 15
View file @ efc2d907
......@@ -477,7 +477,8 @@ VTKWriter::_write(const MeshType& mesh,
if constexpr (is_polygonal_mesh_v<MeshType>) {
fout << mesh.numberOfNodes();
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
fout << mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces();
fout << mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
(mesh.degree() > 2) * mesh.numberOfCells();
} else {
throw NotImplementedError("unexpected mesh type");
}
......@@ -521,7 +522,26 @@ VTKWriter::_write(const MeshType& mesh,
item_value_variant);
}
if constexpr (is_polynomial_mesh_v<MeshType>) {
Array<int> node_number(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces());
size_t number_of_inner_nodes = 0;
if (mesh.degree() > 2) {
const auto cell_type = mesh.connectivity().cellType();
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
switch (cell_type[cell_id]) {
case CellType::Triangle: {
number_of_inner_nodes += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
break;
}
case CellType::Quadrangle: {
number_of_inner_nodes += (mesh.degree() - 2) * (mesh.degree() - 2);
break;
}
default: {
throw NormalError("invalid cell type for polynomal mesh of degree > 2");
}
}
}
}
Array<int> node_number(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() + number_of_inner_nodes);
node_number.fill(0);
_write_array(fout, "node_number", node_number, serialize_data_list);
}
......@@ -545,7 +565,8 @@ VTKWriter::_write(const MeshType& mesh,
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
const NodeValue<const Rd>& xr = mesh.xr();
const FaceArray<const Rd>& xl = mesh.xl();
Array<TinyVector<3>> positions(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces());
Array<TinyVector<3>> positions(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
(mesh.degree() > 2) * mesh.numberOfCells());
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
......@@ -570,6 +591,25 @@ VTKWriter::_write(const MeshType& mesh,
}
}
});
if (mesh.degree() > 2) {
const size_t cell_begining = mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i_node = 0; i_node < number_of_node_per_face; ++i_node) {
const size_t index = cell_begining + cell_id;
for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
positions[index][i] = xj[cell_id][i];
}
for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
positions[index][i] = 0;
}
}
});
}
_write_array(fout, "Positions", positions, serialize_data_list);
} else {
throw NotImplementedError("unexpected mesh type");
......@@ -596,25 +636,31 @@ VTKWriter::_write(const MeshType& mesh,
_write_array(fout, "offsets", offsets, serialize_data_list);
}
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
if (mesh.degree() > 2) {
throw NotImplementedError("only parabolic meshes are supported");
}
const size_t degree = mesh.degree();
const size_t number_of_vertices = mesh.numberOfNodes();
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const auto& cell_type = mesh.connectivity().cellType();
Array<unsigned int> offsets(mesh.numberOfCells());
unsigned int offset = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const auto& cell_nodes = cell_to_node_matrix[cell_id];
offset += degree * cell_nodes.size();
if (degree > 2) {
if (cell_type[cell_id] == CellType::Triangle) {
offset += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
} else if (cell_type[cell_id] == CellType::Quadrangle) {
offset += (mesh.degree() - 1) * (mesh.degree() - 1);
}
}
offsets[cell_id] = offset;
}
_write_array(fout, "offsets", offsets, serialize_data_list);
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
const auto& cell_type = mesh.connectivity().cellType();
const auto& cell_face_is_reversed_matrix = mesh.connectivity().cellFaceIsReversed();
Array<NodeId::base_type> nodes{offset};
size_t index = 0;
......@@ -624,14 +670,51 @@ VTKWriter::_write(const MeshType& mesh,
for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
nodes[index++] = cell_nodes[i_node];
}
if (degree == 2) {
if (degree > 1) {
auto cell_face_is_reversed = cell_face_is_reversed_matrix[cell_id];
if (cell_type[cell_id] == CellType::Triangle) {
for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
nodes[index++] = number_of_vertices + cell_faces[(i_face + 2) % 3];
const size_t i_face_local_number = (i_face + 2) % 3;
if (cell_face_is_reversed[i_face_local_number]) {
for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
nodes[index++] =
number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
}
} else {
for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
nodes[index++] =
number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
}
}
}
if (degree > 2) {
for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 2) / 2; ++i) {
nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
}
}
} else {
for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
nodes[index++] = number_of_vertices + cell_faces[i_face];
if ((i_face >= 2) xor cell_face_is_reversed[i_face]) {
for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
}
} else {
for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
}
}
}
if (degree > 2) {
if (cell_type[cell_id] != CellType::Quadrangle) {
std::ostringstream error_msg;
error_msg << "cell type: " << name(cell_type[cell_id])
<< " is not supported for VTK output of polynomial meshes of degree " << mesh.degree()
<< '\n';
throw NormalError(error_msg.str());
}
for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 1); ++i) {
nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
}
}
}
}
......@@ -673,7 +756,7 @@ VTKWriter::_write(const MeshType& mesh,
} else if (element_degree == 2) {
types[j] = 22;
} else {
throw NotImplementedError("unexpected degree for mesh writer");
types[j] = 69;
}
}
break;
......@@ -687,7 +770,7 @@ VTKWriter::_write(const MeshType& mesh,
} else if (element_degree == 2) {
types[j] = 23;
} else {
throw NotImplementedError("unexpected degree for mesh writer");
types[j] = 70;
}
}
break;
......@@ -696,12 +779,12 @@ VTKWriter::_write(const MeshType& mesh,
if constexpr (is_polygonal_mesh_v<MeshType>) {
types[j] = 7;
} else if constexpr (is_polynomial_mesh_v<MeshType>) {
if (element_degree == 2) {
if (element_degree == 1) {
types[j] = 7;
} else if (element_degree == 2) {
types[j] = 36;
} else {
throw NotImplementedError("unexpected degree for mesh writer");
throw NormalError("output for polynomial cells of degree > 2 is not available for polygonal cells");
}
}
break;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment