diff --git a/cpppetsc/tools/generate_xdmf.py b/cpppetsc/tools/generate_xdmf.py
index 34de96cc5116b2545498ff6a6d2522d7bc9d0a68..1c5fca94474ac79fe82ed05d31849831c838b05c 100755
--- a/cpppetsc/tools/generate_xdmf.py
+++ b/cpppetsc/tools/generate_xdmf.py
@@ -79,6 +79,8 @@ def number_of_corners_to_type(
'Triangle'
>>> number_of_corners_to_type(4, 2)
'Quadrilateral'
+ >>> number_of_corners_to_type(6, 2)
+ 'Tri_6'
>>> number_of_corners_to_type(2, 2)
Traceback (most recent call last):
generate_xdmf.UnsupportedElementType
@@ -104,6 +106,7 @@ def number_of_corners_to_type(
2: {
3: "Triangle",
4: "Quadrilateral",
+ 6: "Tri_6",
},
3: {
4: "Tetrahedron",
diff --git a/elements/src/CMakeLists.txt b/elements/src/CMakeLists.txt
index 9b93172e53f14671bc5b3db1b7bf69fd6f93f0b4..3c1af024fc49f745bc9d23899f56d383c12196ca 100644
--- a/elements/src/CMakeLists.txt
+++ b/elements/src/CMakeLists.txt
@@ -70,6 +70,7 @@ add_library(${PROJECT_NAME}-elements
shape/Quad4.cc
shape/Seg2.cc
shape/Tri3.cc
+ shape/Tri6.cc
shape/ShapeBase.cc
shape/ValueTrait.cc
shape/automatic_gradients.cc
diff --git a/elements/src/include/ae108/elements/shape/Tri6.h b/elements/src/include/ae108/elements/shape/Tri6.h
new file mode 100644
index 0000000000000000000000000000000000000000..33cb73bef4dd54c61254d8d0ced5526a496c2950
--- /dev/null
+++ b/elements/src/include/ae108/elements/shape/Tri6.h
@@ -0,0 +1,62 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "ae108/elements/shape/GradientTrait.h"
+#include "ae108/elements/shape/PointTrait.h"
+#include "ae108/elements/shape/ShapeBase.h"
+#include "ae108/elements/shape/ValueTrait.h"
+#include <cstddef>
+
+namespace ae108 {
+namespace elements {
+namespace shape {
+
+// refer to Cook et. al (2002), "Concepts and applications of Finite Element
+// Analysis", 4th ed., p.265
+struct Tri6 : ShapeBase<std::size_t, double, 2, 6> {};
+
+AE108_ELEMENTS_SHAPE_DEFINE_VALUE(Tri6, {{
+ (1 - xi[0] - xi[1]) *
+ (1 - 2 * xi[0] - 2 * xi[1]),
+ xi[0] * (2 * xi[0] - 1),
+ xi[1] * (2 * xi[1] - 1),
+ 4 * xi[0] * (1 - xi[0] - xi[1]),
+ 4 * xi[0] * xi[1],
+ 4 * xi[1] * (1 - xi[0] - xi[1]),
+ }});
+
+AE108_ELEMENTS_SHAPE_DEFINE_GRADIENTS(
+ Tri6, {{
+ {{4. * (xi[0] + xi[1]) - 3., 4. * (xi[0] + xi[1]) - 3.}},
+ {{4. * xi[0] - 1., 0.}},
+ {{0., 4 * xi[1] - 1}},
+ {{4. * (1 - 2 * xi[0] - xi[1]), -4. * xi[0]}},
+ {{4. * xi[1], 4. * xi[0]}},
+ {{-4. * xi[1], 4. * (1 - xi[0] - 2. * xi[1])}},
+ }});
+
+AE108_ELEMENTS_SHAPE_DEFINE_POINTS(Tri6, {{
+ {{0., 0.}},
+ {{1., 0.}},
+ {{0., 1.}},
+ {{0.5, 0.}},
+ {{0.5, 0.5}},
+ {{0., 0.5}},
+ }});
+
+} // namespace shape
+} // namespace elements
+} // namespace ae108
\ No newline at end of file
diff --git a/elements/src/shape/Tri6.cc b/elements/src/shape/Tri6.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a7a67797db492553f0c688b0be2f093e4e5db39d
--- /dev/null
+++ b/elements/src/shape/Tri6.cc
@@ -0,0 +1,15 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/elements/shape/Tri6.h"
\ No newline at end of file
diff --git a/elements/test/CMakeLists.txt b/elements/test/CMakeLists.txt
index 11ff253edf6b70b34ccc10e8312fba744b9587bd..92574cd20de2b0f5b5d60aba1c58cbbeba23b839 100644
--- a/elements/test/CMakeLists.txt
+++ b/elements/test/CMakeLists.txt
@@ -30,6 +30,7 @@ add_executable(${PROJECT_NAME}-ElementsTests
shape/Quad4_Test.cc
shape/Seg2_Test.cc
shape/Tri3_Test.cc
+ shape/Tri6_Test.cc
tensor/Tensor_Test.cc
tensor/as_matrix_of_columns_Test.cc
tensor/as_matrix_of_rows_Test.cc
diff --git a/elements/test/shape/Tri6_Test.cc b/elements/test/shape/Tri6_Test.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d4f9a94d68a96e0e638da7f973a54c96a929ddec
--- /dev/null
+++ b/elements/test/shape/Tri6_Test.cc
@@ -0,0 +1,52 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "Shape_Test.h"
+#include "ae108/elements/shape/Tri6.h"
+#include "ae108/elements/shape/get_points.h"
+#include <gmock/gmock.h>
+
+using testing::DoubleEq;
+using testing::ElementsAre;
+using testing::Test;
+using testing::Types;
+
+namespace ae108 {
+namespace elements {
+namespace shape {
+namespace {
+
+using Configurations = Types<Tri6>;
+
+INSTANTIATE_TYPED_TEST_CASE_P(Tri6_Test, Shape_Test, Configurations);
+
+struct Tri6_Test : Test {
+ using Shape = Tri6;
+};
+
+TEST_F(Tri6_Test, points_are_correct) {
+ const auto points = get_points<Shape>();
+
+ EXPECT_THAT(points, ElementsAre(ElementsAre(DoubleEq(0.), DoubleEq(0.)),
+ ElementsAre(DoubleEq(1.), DoubleEq(0.)),
+ ElementsAre(DoubleEq(0.), DoubleEq(1.)),
+ ElementsAre(DoubleEq(0.5), DoubleEq(0.)),
+ ElementsAre(DoubleEq(0.5), DoubleEq(0.5)),
+ ElementsAre(DoubleEq(0.), DoubleEq(0.5))));
+}
+
+} // namespace
+} // namespace shape
+} // namespace elements
+} // namespace ae108
\ No newline at end of file
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index b5b59853892cca1d74ab7c4624eda0e58ea491e5..5536ee6c6d8d0c26450bd131375b269b9d4a2a99 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -48,3 +48,10 @@ add_executable(${PROJECT_NAME}-Input Input.cc)
target_link_libraries(${PROJECT_NAME}-Input
PRIVATE ae108::cpppetsc
)
+
+add_executable(${PROJECT_NAME}-QuadraticTriangles QuadraticTriangles.cc)
+target_link_libraries(${PROJECT_NAME}-QuadraticTriangles
+ PRIVATE ae108::solve
+ PRIVATE ae108::assembly
+ PRIVATE ae108::elements
+)
diff --git a/examples/QuadraticTriangles.cc b/examples/QuadraticTriangles.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4b8c3b03ef032d53af2d11c39085cac589fcf7e1
--- /dev/null
+++ b/examples/QuadraticTriangles.cc
@@ -0,0 +1,199 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/assembly/Assembler.h"
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/Vector.h"
+#include "ae108/cpppetsc/Viewer.h"
+#include "ae108/cpppetsc/createVectorFromSource.h"
+#include "ae108/cpppetsc/setName.h"
+#include "ae108/cpppetsc/writeToViewer.h"
+#include "ae108/elements/CoreElement.h"
+#include "ae108/elements/embedding/IsoparametricEmbedding.h"
+#include "ae108/elements/integrator/IsoparametricIntegrator.h"
+#include "ae108/elements/materialmodels/Hookean.h"
+#include "ae108/elements/mesh/generate_quadratic_triangle_mesh.h"
+#include "ae108/elements/quadrature/Quadrature.h"
+#include "ae108/elements/shape/Tri6.h"
+#include "ae108/solve/NonlinearSolver.h"
+#include <array>
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Mesh = cpppetsc::Mesh<Policy>;
+using Vector = cpppetsc::Vector<Policy>;
+using Viewer = cpppetsc::Viewer<Policy>;
+using BoundaryCondition = cpppetsc::MeshBoundaryCondition<Mesh>;
+
+namespace embedding = ae108::elements::embedding;
+namespace integrator = ae108::elements::integrator;
+namespace materialmodels = ae108::elements::materialmodels;
+namespace quadrature = ae108::elements::quadrature;
+namespace shape = ae108::elements::shape;
+namespace mesh = ae108::elements::mesh;
+
+// In this example we will simulate a mesh of quadratic triangles that
+// that we will generate later.
+
+// Let's specify the parameters of this mesh.
+
+constexpr auto number_of_vertices_per_element = Mesh::size_type{6};
+constexpr auto dimension = Mesh::size_type{2};
+constexpr auto dof_per_vertex = Mesh::size_type{2};
+
+// Let's specify how the elements behave when deformed.
+
+// First, we select a Hookean (linear elastic) material model in 2D.
+
+using MaterialModel = materialmodels::Hookean<dimension>;
+
+// We'll choose the Tri6 shape functions (6 shape functions) that are defined
+// in 2D reference space.
+
+using Shape = shape::Tri6;
+
+// Since also our physical space is 2D we may use the isoparametric embedding.
+
+using Embedding = embedding::IsoparametricEmbedding<Shape>;
+
+// Let's choose a simple quadrature rule in reference space.
+
+using Quadrature = quadrature::Quadrature<quadrature::QuadratureType::Simplex,
+ dimension, 2 /* order */>;
+
+// We use an "Integrator" to integrate in physical space. Of course, it depends
+// on the selected quadrature rule and the embedding. In the case of the
+// isoparametric embedding, this embedding is specified by the "Shape".
+
+using Integrator = integrator::IsoparametricIntegrator<Shape, Quadrature>;
+
+// Now we are ready to select our element type.
+
+using Element = elements::CoreElement<MaterialModel, Integrator>;
+
+// We will assemble e.g. energy using a collection of elements. This is done by
+// the assembler. (The list DefaultFeaturePlugins contain the features (e.g.
+// energy) that most elements support.)
+
+using Assembler =
+ assembly::Assembler<Element, assembly::DefaultFeaturePlugins, Policy>;
+
+// Our goal is to minimize the energy. This is done by the solver.
+
+using Solver = solve::NonlinearSolver<Assembler>;
+
+int main(int argc, char **argv) {
+ // PETSc must be initialized before using it.
+ Policy::handleError(PetscInitialize(&argc, &argv, NULL, NULL));
+
+ // We use a scope around our computation to make sure everything is cleaned up
+ // before we call PetscFinalize.
+ {
+ // Now we generate a geometry of 4 elements connected in x-direction.
+ // *-----*-----*
+ // | \ | \ |
+ // *-----*-----*
+ const auto geometry =
+ mesh::generate_quadratic_triangle_mesh({{2., 1.}}, {{2, 1}});
+
+ const auto mesh =
+ Mesh::fromConnectivity(dimension, geometry.connectivity(),
+ geometry.number_of_positions(), dof_per_vertex);
+ auto assembler = Assembler();
+
+ const auto model = MaterialModel(1.0, 0.);
+
+ // Depending on whether we use MPI, our mesh may be distributed and not all
+ // elements are present on this computational node.
+
+ // Let's add those elements that are "local" to the assembler.
+
+ for (const auto &element : mesh.localElements()) {
+ const auto vertexIndices = element.vertexIndices();
+ assembler.emplaceElement(
+ element, model,
+ Integrator(Embedding(Embedding::Collection<Embedding::PhysicalPoint>{{
+ geometry.position_of_vertex(vertexIndices.at(0)),
+ geometry.position_of_vertex(vertexIndices.at(1)),
+ geometry.position_of_vertex(vertexIndices.at(2)),
+ geometry.position_of_vertex(vertexIndices.at(3)),
+ geometry.position_of_vertex(vertexIndices.at(4)),
+ geometry.position_of_vertex(vertexIndices.at(5)),
+ }})));
+ }
+
+ // We need to create a solver. We do not use the time, so we can set it to
+ // zero.
+
+ const auto solver = Solver(&mesh);
+ const auto time = Element::Time{0.};
+
+ // Before we can produce meaningful results, we need to specify boundary
+ // conditions. Let's fix the nodes at x=0 and pull on the nodes at x=2.
+
+ std::vector<BoundaryCondition> boundary_conditions;
+ for (const auto &vertex : mesh.localVertices()) {
+ const auto position = geometry.position_of_vertex(vertex.index());
+ const auto tolerance = .1;
+ if (std::abs(position[0] - 0.) < tolerance) {
+ // The displacement in x direction is zero.
+ boundary_conditions.push_back({vertex, 0, 0.});
+ // The displacement in y direction is zero.
+ boundary_conditions.push_back({vertex, 1, 0.});
+ } else if (std::abs(position[0] - 2.) < tolerance) {
+ // The displacement in x direction is .5.
+ boundary_conditions.push_back({vertex, 0, .5});
+ // The displacement in y direction is 0.
+ boundary_conditions.push_back({vertex, 1, 0.});
+ }
+ }
+
+ // We are ready to minimize the energy.
+
+ auto result = solver.computeSolution(
+ boundary_conditions, Vector::fromGlobalMesh(mesh), time, &assembler);
+
+ // As a final step, we write the result to a file.
+
+ // We may view this result in e.g. ParaView by generating an XDMF file with
+ // the `generate_xdmf.py` script in the repository, and opening this file
+ // with ParaView ("XDMF Reader").
+
+ // First we collect the coordinates in a vector.
+ using DataSource = std::function<void(Mesh::size_type, Mesh::value_type *)>;
+ auto coordinates = cpppetsc::createVectorFromSource(
+ mesh, dimension,
+ DataSource(
+ [&](const Mesh::size_type index, Mesh::value_type *const out) {
+ const auto &position = geometry.position_of_vertex(index);
+ std::copy(position.begin(), position.end(), out);
+ }));
+ cpppetsc::setName("coordinates", &coordinates);
+
+ // Now we write the mesh to a file.
+ auto viewer = Viewer::fromHdf5FilePath("quadratic_triangles.ae108",
+ Viewer::Mode::write);
+ cpppetsc::writeToViewer(mesh, coordinates, &viewer);
+
+ // Let's write the result to the file.
+ cpppetsc::setName("result", &result);
+ cpppetsc::writeToViewer(result, &viewer);
+
+ fprintf(stderr, "The data has been written to the file "
+ "\"quadratic_triangles.ae108\".\n");
+ }
+ Policy::handleError(PetscFinalize());
+}