Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for Serendipity (order 2) coordinate elements #3662

Merged
merged 6 commits into from
Mar 15, 2025
+60 −41
Merged

Add support for Serendipity (order 2) coordinate elements #3662

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
105caec
Add 20 node hex and 8 node brick as coordinate element.
jorgensd Mar 12, 2025
9319770
Simplify
garth-wells Mar 14, 2025
08f6226
Simplify
garth-wells Mar 14, 2025
d9e81f0
Revert
garth-wells Mar 15, 2025
73021ea
Update
garth-wells Mar 15, 2025
cb13a2b
Merge branch 'main' into dokken/alternative-coordinate-element
garth-wells Mar 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 21 additions & 19 deletions cpp/dolfinx/io/cells.cpp
View file
Original file line number Diff line number Diff line change
@@ -5,11 +5,14 @@
// SPDX-License-Identifier: LGPL-3.0-or-later

#include "cells.h"
#include <array>
#include <cstdint>
#include <dolfinx/common/log.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/cell_types.h>
#include <numeric>
#include <stdexcept>
#include <vector>

using namespace dolfinx;
namespace
@@ -414,38 +417,40 @@ std::vector<std::uint16_t> vtk_pyramid(int num_nodes)
//-----------------------------------------------------------------------------
std::vector<std::uint16_t> vtk_quadrilateral(int num_nodes)
{
const int n = cell_degree(mesh::CellType::quadrilateral, num_nodes);
std::vector<std::uint16_t> map(num_nodes);
std::vector<std::uint16_t> map;
map.reserve(num_nodes);

// Vertices
map[0] = 0;
map[1] = 1;
map[2] = 3;
map[3] = 2;

int j = 4;
map.insert(map.begin(), {0, 1, 3, 2});

const int edge_nodes = n - 1;
// Special handling for second-order serendipity
constexpr mesh::CellType cell = mesh::CellType::quadrilateral;
int edge_nodes = num_nodes == 8 ? 1 : cell_degree(cell, num_nodes) - 1;

// Edges
for (int k = 0; k < edge_nodes; ++k)
map[j++] = 4 + k;
map.push_back(4 + k);
for (int k = 0; k < edge_nodes; ++k)
map[j++] = 4 + 2 * edge_nodes + k;
map.push_back(4 + 2 * edge_nodes + k);
for (int k = 0; k < edge_nodes; ++k)
map[j++] = 4 + 3 * edge_nodes + k;
map.push_back(4 + 3 * edge_nodes + k);
for (int k = 0; k < edge_nodes; ++k)
map[j++] = 4 + edge_nodes + k;
map.push_back(4 + edge_nodes + k);

// Face
for (int k = 0; k < edge_nodes * edge_nodes; ++k)
map[j++] = 4 + edge_nodes * 4 + k;
map.push_back(4 + edge_nodes * 4 + k);

return map;
}
//-----------------------------------------------------------------------------
std::vector<std::uint16_t> vtk_hexahedron(int num_nodes)
{
std::uint16_t n = cell_degree(mesh::CellType::hexahedron, num_nodes);
// Special handling for second order serendipity
constexpr mesh::CellType cell = mesh::CellType::hexahedron;
int edge_nodes = num_nodes == 20 ? 1 : cell_degree(cell, num_nodes) - 1;
int face_nodes = num_nodes == 20 ? 0 : edge_nodes * edge_nodes;
int volume_nodes = face_nodes * edge_nodes;

std::vector<std::uint16_t> map(num_nodes);

@@ -462,16 +467,14 @@ std::vector<std::uint16_t> vtk_hexahedron(int num_nodes)
// Edges
int j = 8;
int base = 8;
const int edge_nodes = n - 1;
const std::vector<int> edges = {0, 3, 5, 1, 8, 10, 11, 9, 2, 4, 7, 6};
const std::array edges = {0, 3, 5, 1, 8, 10, 11, 9, 2, 4, 7, 6};
for (int e : edges)
{
for (int i = 0; i < edge_nodes; ++i)
map[j++] = base + edge_nodes * e + i;
}
base += 12 * edge_nodes;

const int face_nodes = edge_nodes * edge_nodes;
const std::vector<int> faces = {2, 3, 1, 4, 0, 5};
for (int f : faces)
{
@@ -480,7 +483,6 @@ std::vector<std::uint16_t> vtk_hexahedron(int num_nodes)
}
base += 6 * face_nodes;

const int volume_nodes = face_nodes * edge_nodes;
for (int i = 0; i < volume_nodes; ++i)
map[j++] = base + i;

30 changes: 17 additions & 13 deletions cpp/dolfinx/io/xdmf_utils.cpp
View file
Original file line number Diff line number Diff line change
@@ -38,9 +38,11 @@ xdmf_utils::get_cell_type(const pugi::xml_node& topology_node)
{"tetrahedron", {"tetrahedron", 1}},
{"tetrahedron_10", {"tetrahedron", 2}},
{"quadrilateral", {"quadrilateral", 1}},
{"quadrilateral_8", {"quadrilateral", 2}},
{"quadrilateral_9", {"quadrilateral", 2}},
{"quadrilateral_16", {"quadrilateral", 3}},
{"hexahedron", {"hexahedron", 1}},
{"hexahedron_20", {"hexahedron", 2}},
{"wedge", {"prism", 1}},
{"hexahedron_27", {"hexahedron", 2}}};

@@ -169,19 +171,21 @@ std::int64_t xdmf_utils::get_num_cells(const pugi::xml_node& topology_node)
std::string xdmf_utils::vtk_cell_type_str(mesh::CellType cell_type,
int num_nodes)
{
static const std::map<mesh::CellType, std::map<int, std::string>> vtk_map = {
{mesh::CellType::point, {{1, "PolyVertex"}}},
{mesh::CellType::interval, {{2, "PolyLine"}, {3, "Edge_3"}}},
{mesh::CellType::triangle,
{{3, "Triangle"}, {6, "Triangle_6"}, {10, "Triangle_10"}}},
{mesh::CellType::quadrilateral,
{{4, "Quadrilateral"},
{9, "Quadrilateral_9"},
{16, "Quadrilateral_16"}}},
{mesh::CellType::prism, {{6, "Wedge"}}},
{mesh::CellType::tetrahedron,
{{4, "Tetrahedron"}, {10, "Tetrahedron_10"}, {20, "Tetrahedron_20"}}},
{mesh::CellType::hexahedron, {{8, "Hexahedron"}, {27, "Hexahedron_27"}}}};
static const std::map<mesh::CellType, std::map<int, std::string>> vtk_map
= {{mesh::CellType::point, {{1, "PolyVertex"}}},
{mesh::CellType::interval, {{2, "PolyLine"}, {3, "Edge_3"}}},
{mesh::CellType::triangle,
{{3, "Triangle"}, {6, "Triangle_6"}, {10, "Triangle_10"}}},
{mesh::CellType::quadrilateral,
{{4, "Quadrilateral"},
{8, "Quadrilateral_8"},
{9, "Quadrilateral_9"},
{16, "Quadrilateral_16"}}},
{mesh::CellType::prism, {{6, "Wedge"}}},
{mesh::CellType::tetrahedron,
{{4, "Tetrahedron"}, {10, "Tetrahedron_10"}, {20, "Tetrahedron_20"}}},
{mesh::CellType::hexahedron,
{{8, "Hexahedron"}, {27, "Hexahedron_27"}, {20, "Hexahedron_20"}}}};

// Get cell family
auto cell = vtk_map.find(cell_type);
31 changes: 22 additions & 9 deletions python/dolfinx/io/utils.py
View file
Original file line number Diff line number Diff line change
@@ -21,7 +21,7 @@
from dolfinx.cpp.io import perm_gmsh as cell_perm_gmsh
from dolfinx.cpp.io import perm_vtk as cell_perm_vtk
from dolfinx.fem import Function
from dolfinx.mesh import Geometry, GhostMode, Mesh, MeshTags
from dolfinx.mesh import CellType, Geometry, GhostMode, Mesh, MeshTags

__all__ = ["VTKFile", "XDMFFile", "cell_perm_gmsh", "cell_perm_vtk", "distribute_entity_data"]

@@ -187,21 +187,34 @@ def read_mesh(
cells = super().read_topology_data(name, xpath)
x = super().read_geometry_data(name, xpath)

# Build the mesh
cmap = _cpp.fem.CoordinateElement_float64(cell_shape, cell_degree)
msh = _cpp.mesh.create_mesh(
self.comm, cells, cmap, x, _cpp.mesh.create_cell_partitioner(ghost_mode)
)
msh.name = name
domain = ufl.Mesh(
basix.ufl.element(
# Get coordinate element, special handling for second order serendipity.
num_nodes_per_cell = cells.shape[1]
if (cell_shape == CellType.quadrilateral and num_nodes_per_cell == 8) or (
cell_shape == CellType.hexahedron and num_nodes_per_cell == 20
):
el = basix.ufl.element(
basix.ElementFamily.serendipity,
cell_shape.name,
2,
)
cmap = _cpp.fem.CoordinateElement_float64(el.basix_element._e)
basix_el = basix.ufl.blocked_element(el, shape=(x.shape[1],))
else:
basix_el = basix.ufl.element(
"Lagrange",
cell_shape.name,
cell_degree,
basix.LagrangeVariant.unset,
shape=(x.shape[1],),
)
cmap = _cpp.fem.CoordinateElement_float64(cell_shape, cell_degree)

# Build the mesh
msh = _cpp.mesh.create_mesh(
self.comm, cells, cmap, x, _cpp.mesh.create_cell_partitioner(ghost_mode)
)
msh.name = name
domain = ufl.Mesh(basix_el)
return Mesh(msh, domain)

def read_meshtags(