From 75e11993c02768c05d1f062fe1c58f4934473d16 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 16:31:58 +0200 Subject: [PATCH 1/8] make naming consistent with dune-vtk and transform examples into tests --- dune/gmsh4/filereader.hh | 1 - dune/gmsh4/{reader.hh => gmsh4reader.hh} | 2 +- .../{reader.impl.hh => gmsh4reader.impl.hh} | 4 +- dune/gmsh4/gridcreators/common.hh | 1 - .../gridcreators/continuousgridcreator.hh | 27 +- dune/gmsh4/gridcreators/derivedgridcreator.hh | 10 +- .../gridcreators/discontinuousgridcreator.hh | 50 +- .../gmsh4/gridcreators/lagrangegridcreator.hh | 29 +- .../gmsh4/gridcreators/parallelgridcreator.hh | 4 +- dune/gmsh4/gridcreators/serialgridcreator.hh | 40 +- dune/gmsh4/types.cc | 1 - dune/gmsh4/utility/CMakeLists.txt | 1 + dune/gmsh4/utility/filesystem.cc | 2 +- dune/gmsh4/utility/lagrangepoints.hh | 716 ++++-------------- dune/gmsh4/utility/lagrangepoints.impl.hh | 507 +++++++++++++ src/ALUGridTests.cc | 32 +- src/CMakeLists.txt | 27 +- src/FoamGridTests.cc | 20 +- src/gmsh4reader.cc | 8 +- 19 files changed, 746 insertions(+), 736 deletions(-) rename dune/gmsh4/{reader.hh => gmsh4reader.hh} (99%) rename dune/gmsh4/{reader.impl.hh => gmsh4reader.impl.hh} (99%) create mode 100644 dune/gmsh4/utility/lagrangepoints.impl.hh diff --git a/dune/gmsh4/filereader.hh b/dune/gmsh4/filereader.hh index 918ffca..cd9daee 100644 --- a/dune/gmsh4/filereader.hh +++ b/dune/gmsh4/filereader.hh @@ -11,7 +11,6 @@ namespace Dune { namespace Gmsh4 { - template class FileReader { diff --git a/dune/gmsh4/reader.hh b/dune/gmsh4/gmsh4reader.hh similarity index 99% rename from dune/gmsh4/reader.hh rename to dune/gmsh4/gmsh4reader.hh index 46e1817..8be484e 100644 --- a/dune/gmsh4/reader.hh +++ b/dune/gmsh4/gmsh4reader.hh @@ -274,4 +274,4 @@ namespace Dune } // end namespace Dune -#include "reader.impl.hh" +#include "gmsh4reader.impl.hh" diff --git a/dune/gmsh4/reader.impl.hh b/dune/gmsh4/gmsh4reader.impl.hh similarity index 99% rename from dune/gmsh4/reader.impl.hh rename to dune/gmsh4/gmsh4reader.impl.hh index 2cb8890..f7e6954 100644 --- a/dune/gmsh4/reader.impl.hh +++ b/dune/gmsh4/gmsh4reader.impl.hh @@ -10,8 +10,8 @@ #include #include -//Helper-function to deal with endianness in binary msh-files. -void swapBytes(char *array, int size) +// Helper-function to deal with endianness in binary msh-files. +inline void swapBytes (char *array, int size) { char *tmp = new char[size]; memcpy(tmp, array, size); diff --git a/dune/gmsh4/gridcreators/common.hh b/dune/gmsh4/gridcreators/common.hh index 5a6e440..5c563fe 100644 --- a/dune/gmsh4/gridcreators/common.hh +++ b/dune/gmsh4/gridcreators/common.hh @@ -6,7 +6,6 @@ namespace Dune { namespace Gmsh4 { - template using HasInsertVertex = decltype( std::declval().insertVertex(std::declval()...) ); diff --git a/dune/gmsh4/gridcreators/continuousgridcreator.hh b/dune/gmsh4/gridcreators/continuousgridcreator.hh index c6c1b7b..6fc95a7 100644 --- a/dune/gmsh4/gridcreators/continuousgridcreator.hh +++ b/dune/gmsh4/gridcreators/continuousgridcreator.hh @@ -34,8 +34,8 @@ namespace Dune template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { vertexMap_.resize(nodeTagRange.second - nodeTagRange.first + 1); vertexShift_ = nodeTagRange.first; @@ -55,9 +55,9 @@ namespace Dune template void insertElementsImpl (std::size_t numElements, - std::pair elementTagRange, - std::vector const& entityBlocks, - BoundaryEntities const& boundaryEntities) + std::pair elementTagRange, + std::vector const& entityBlocks, + BoundaryEntities const& boundaryEntities) { std::vector connectivity; std::size_t cornerIndex = 0; @@ -69,22 +69,7 @@ namespace Dune auto type = Gmsh4::to_geometry(entityBlock.elementType); Gmsh4::CellType cell{type}; - //this segment was probably added due to a misunderstanding about gmsh's boundary-handling - /*if (entityBlock.entityDim == Grid::dimension-1) { //boundary - if (boundaryEntities.count(entityBlock.entityTag)) { - auto refElem = referenceElement(cell.type()); - connectivity.resize(refElem.size(Grid::dimension-1)); - - for (auto const& element : entityBlock.elements) { - assert(element.nodes.size() >= connectivity.size()); - for (std::size_t j = 0; j < connectivity.size(); ++j) - connectivity[cell.permutation(j)] = vertexMap_[element.nodes[j] - vertexShift_]; - - factory().insertBoundarySegment(connectivity); - } - } - } - else */ + if (entityBlock.entityDim == Grid::dimension) { //element auto refElem = referenceElement(cell.type()); connectivity.resize(refElem.size(Grid::dimension)); diff --git a/dune/gmsh4/gridcreators/derivedgridcreator.hh b/dune/gmsh4/gridcreators/derivedgridcreator.hh index 92042ff..573ec26 100644 --- a/dune/gmsh4/gridcreators/derivedgridcreator.hh +++ b/dune/gmsh4/gridcreators/derivedgridcreator.hh @@ -31,17 +31,17 @@ namespace Dune template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { gridCreator_.insertVertices(numNodes, nodeTagRange, entityBlocks); } template void insertElementsImpl (std::size_t numElements, - std::pair elementTagRange, - std::vector const& entityBlocks, - BoundaryEntities const& boundaryEntities) + std::pair elementTagRange, + std::vector const& entityBlocks, + BoundaryEntities const& boundaryEntities) { gridCreator_.insertElements(numElements, elementTagRange, entityBlocks, boundaryEntities); } diff --git a/dune/gmsh4/gridcreators/discontinuousgridcreator.hh b/dune/gmsh4/gridcreators/discontinuousgridcreator.hh index e192548..58ae02f 100644 --- a/dune/gmsh4/gridcreators/discontinuousgridcreator.hh +++ b/dune/gmsh4/gridcreators/discontinuousgridcreator.hh @@ -47,55 +47,19 @@ namespace Dune template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { - std::cout << "WARNING! Dune::Gmsh4::DiscontinuousGridCreator is not implemented yet!" << std::endl; - // points_ = &points; - // uniquePoints_.clear(); - // std::size_t idx = 0; - - // for (auto const& p : points) { - // auto b = uniquePoints_.emplace(std::make_pair(p,idx)); - // if (b.second) { - // factory().insertVertex(p); - // ++idx; - // } - // } + DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::DiscontinuousGridCreator is not implemented yet!"); } template void insertElementsImpl (std::size_t numElements, - std::pair elementTagRange, - std::vector const& entityBlocks, - BoundaryEntities const& boundaryEntities) + std::pair elementTagRange, + std::vector const& entityBlocks, + BoundaryEntities const& boundaryEntities) { - std::cout << "WARNING! Dune::Gmsh4::DiscontinuousGridCreator is not implemented yet!" << std::endl; - // assert(points_ != nullptr); - // std::size_t idx = 0; - // for (std::size_t i = 0; i < types.size(); ++i) { - // auto type = Vtk::to_geometry(types[i]); - // Vtk::CellType cellType{type}; - - // int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]); - // assert(nNodes > 0); - // std::vector vtk_cell; vtk_cell.reserve(nNodes); - // for (int j = 0; j < nNodes; ++j) { - // std::size_t v_j = connectivity[idx++]; - // std::size_t new_idx = uniquePoints_[(*points_)[v_j]]; - // vtk_cell.push_back(new_idx); - // } - - // if (cellType.noPermutation()) { - // factory().insertElement(type,vtk_cell); - // } else { - // // apply index permutation - // std::vector cell(nNodes); - // for (int j = 0; j < nNodes; ++j) - // cell[j] = vtk_cell[cellType.permutation(j)]; - // factory().insertElement(type,cell); - // } - // } + DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::DiscontinuousGridCreator is not implemented yet!"); } private: diff --git a/dune/gmsh4/gridcreators/lagrangegridcreator.hh b/dune/gmsh4/gridcreators/lagrangegridcreator.hh index 57e07aa..fa1b05e 100644 --- a/dune/gmsh4/gridcreators/lagrangegridcreator.hh +++ b/dune/gmsh4/gridcreators/lagrangegridcreator.hh @@ -21,7 +21,6 @@ namespace Dune { namespace Gmsh4 { - // \brief Create a grid from data that represents higher (lagrange) cells. /** * The grid is created from the first nodes of a cell parametrization, representing @@ -68,8 +67,8 @@ namespace Dune /// Implementation of the interface function `insertVertices()` template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { vertexMap_.resize(nodeTagRange.second - nodeTagRange.first + 1); vertexShift_ = nodeTagRange.first; @@ -94,9 +93,9 @@ namespace Dune /// Implementation of the interface function `insertElements()` template void insertElementsImpl (std::size_t numElements, - std::pair elementTagRange, - std::vector const& entityBlocks, - BoundaryEntities const& boundaryEntities) + std::pair elementTagRange, + std::vector const& entityBlocks, + BoundaryEntities const& boundaryEntities) { const int dim = GridType::dimension; std::vector connectivity; @@ -109,22 +108,7 @@ namespace Dune auto type = Gmsh4::to_geometry(entityBlock.elementType); Gmsh4::CellType cell{type}; - //this segment was probably added due to a misunderstanding about gmsh's boundary-handling - /*if (entityBlock.entityDim == dim - 1) { //boundary - if (boundaryEntities.count(entityBlock.entityTag)) { - auto refElem = referenceElement(type); - connectivity.resize(refElem.size(dim - 1)); - - for (auto const& element : entityBlock.elements) { - assert(element.nodes.size() >= connectivity.size()); - for (std::size_t j = 0; j < connectivity.size(); ++j) - connectivity[cell.permutation(j)] = vertexMap_[element.nodes[j] - vertexShift_]; - - factory().insertBoundarySegment(connectivity); - } - } - } - else */ + if (entityBlock.entityDim == dim) { //element auto refElem = referenceElement(type); connectivity.resize(refElem.size(dim)); @@ -409,5 +393,4 @@ namespace Dune }; } // end namespace Gmsh4 - } // end namespace Dune diff --git a/dune/gmsh4/gridcreators/parallelgridcreator.hh b/dune/gmsh4/gridcreators/parallelgridcreator.hh index d8a9ec8..9a3c14a 100644 --- a/dune/gmsh4/gridcreators/parallelgridcreator.hh +++ b/dune/gmsh4/gridcreators/parallelgridcreator.hh @@ -35,8 +35,8 @@ namespace Dune template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { std::cout << "WARNING! Dune::Gmsh4::ParallelGridCreator is unfinished!" << std::endl; GlobalCoordinate p; diff --git a/dune/gmsh4/gridcreators/serialgridcreator.hh b/dune/gmsh4/gridcreators/serialgridcreator.hh index e0456ba..06b7c1a 100644 --- a/dune/gmsh4/gridcreators/serialgridcreator.hh +++ b/dune/gmsh4/gridcreators/serialgridcreator.hh @@ -29,48 +29,24 @@ namespace Dune template void insertVerticesImpl (std::size_t numNodes, - std::pair nodeTagRange, - std::vector const& entityBlocks) + std::pair nodeTagRange, + std::vector const& entityBlocks) { - std::cout << "WARNING! Dune::Gmsh4::SerialGridCreator is unfinished!" << std::endl; - // shift_.push_back(points_.size()); - // points_.reserve(points_.size() + points.size()); - // points_.insert(points_.end(), points.begin(), points.end()); + DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::SerialGridCreator is unfinished!"); } template void insertElementsImpl (std::size_t numElements, - std::pair elementTagRange, - std::vector const& entityBlocks, - BoundaryEntities const& boundaryEntities) + std::pair elementTagRange, + std::vector const& entityBlocks, + BoundaryEntities const& boundaryEntities) { - std::cout << "WARNING! Dune::Gmsh4::SerialGridCreator is unfinished!" << std::endl; - // types_.reserve(types_.size() + types.size()); - // types_.insert(types_.end(), types.begin(), types.end()); - - // offsets_.reserve(offsets_.size() + offsets.size()); - // std::transform(offsets.begin(), offsets.end(), std::back_inserter(offsets_), - // [shift=offsets_.empty() ? 0 : offsets_.back()](std::int64_t o) { return o + shift; }); - - // connectivity_.reserve(connectivity_.size() + connectivity.size()); - // std::transform(connectivity.begin(), connectivity.end(), std::back_inserter(connectivity_), - // [shift=shift_.back()](std::int64_t idx) { return idx + shift; }); + DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::SerialGridCreator is unfinished!"); } void insertPiecesImpl (std::vector const& pieces) { - std::cout << "WARNING! Dune::Gmsh4::SerialGridCreator is unfinished!" << std::endl; - // if (this->comm().rank() == 0) { - // Gmsh4Reader pieceReader(*this); - // for (std::string const& piece : pieces) { - // pieceReader.readFromFile(piece, false); - // pieceReader.createGrid(false); - // } - - // DiscontinuousGridCreator creator(this->factory()); - // creator.insertVertices(points_, {}); - // creator.insertElements(types_, offsets_, connectivity_); - // } + DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::SerialGridCreator is unfinished!"); } private: diff --git a/dune/gmsh4/types.cc b/dune/gmsh4/types.cc index 895766e..4ad8934 100644 --- a/dune/gmsh4/types.cc +++ b/dune/gmsh4/types.cc @@ -6,7 +6,6 @@ namespace Dune { namespace Gmsh4 { - GeometryType to_geometry (int elementType) { switch (elementType) { diff --git a/dune/gmsh4/utility/CMakeLists.txt b/dune/gmsh4/utility/CMakeLists.txt index 0206b72..96bab44 100644 --- a/dune/gmsh4/utility/CMakeLists.txt +++ b/dune/gmsh4/utility/CMakeLists.txt @@ -5,5 +5,6 @@ dune_add_library("filesystem" OBJECT install(FILES filesystem.hh lagrangepoints.hh + lagrangepoints.impl.hh string.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gmsh4/utility) diff --git a/dune/gmsh4/utility/filesystem.cc b/dune/gmsh4/utility/filesystem.cc index 7ba4f01..08e5700 100644 --- a/dune/gmsh4/utility/filesystem.cc +++ b/dune/gmsh4/utility/filesystem.cc @@ -43,7 +43,7 @@ void Path::split(std::string p) bool relative = true; trim(p); - Dune::Vtk::split(p.begin(), p.end(), separators.begin(), separators.end(), + Dune::Gmsh4::split(p.begin(), p.end(), separators.begin(), separators.end(), [this,&relative](auto first, auto end) { auto token = std::string(first, end); diff --git a/dune/gmsh4/utility/lagrangepoints.hh b/dune/gmsh4/utility/lagrangepoints.hh index 04c77b5..7ea8752 100644 --- a/dune/gmsh4/utility/lagrangepoints.hh +++ b/dune/gmsh4/utility/lagrangepoints.hh @@ -9,601 +9,191 @@ namespace Dune { - -namespace Gmsh4 -{ - -namespace Impl { - // forward declaration - template - class LagrangePointSetBuilder; -} - - -/// \brief A set of lagrange points compatible with the numbering of VTK and Gmsh -/** - * \tparam K Field-type for the coordinates - * \tparam dim Dimension of the coordinates - **/ -template -class LagrangePointSet - : public EmptyPointSet -{ - using Super = EmptyPointSet; - -public: - static const unsigned int dimension = dim; - - LagrangePointSet (std::size_t order) - : Super(order) - { - assert(order > 0); - } - - /// Fill the lagrange points for the given geometry type - void build (GeometryType gt) - { - assert(gt.dim() == dimension); - builder_(gt, order(), points_); - } - - /// Fill the lagrange points for the given topology type `Topology` - template - bool build () - { - build(GeometryType(Topology{})); - return true; - } - - /// Returns whether the point set support the given topology type `Topology` and can - /// generate point for the given order. - template - static bool supports (std::size_t order) - { - return true; - } - - using Super::order; - -private: - using Super::points_; - Impl::LagrangePointSetBuilder builder_; -}; - - -namespace Impl { - -// Build for lagrange point sets in different dimensions -// Specialized for dim=1,2,3 -template -class LagrangePointSetBuilder -{ -public: - template - void operator()(GeometryType, unsigned int, Points& points) const - { - DUNE_THROW(Dune::NotImplemented, - "Lagrange points not yet implemented for this GeometryType."); - } -}; - -/** - * The implementation of the point set builder is directly derived from VTK. - * Modification are a change in data-types and naming scheme. Additionally - * a LocalKey is created to reflect the concept of a Dune PointSet. - * - * Included is the license of the BSD 3-clause License included in the VTK - * source code from 2020/04/13 in commit b90dad558ce28f6d321420e4a6b17e23f5227a1c - * of git repository https://gitlab.kitware.com/vtk/vtk. - * - Program: Visualization Toolkit - Module: Copyright.txt - - Copyright (c) 1993-2015 Ken Martin, Will Schroeder, Bill Lorensen - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, - this list of conditions and the following disclaimer. - - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - - * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names - of any contributors may be used to endorse or promote products derived - from this software without specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS'' - AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR - ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR - SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER - CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, - OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE - OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - * - **/ - -// Lagrange points on point geometries -template -class LagrangePointSetBuilder -{ - static constexpr int dim = 0; - using LP = LagrangePoint; - using Vec = typename LP::Vector; - using Key = LocalKey; - -public: - template - void operator()(GeometryType gt, int /*order*/, Points& points) const - { - assert(gt.isVertex()); - points.push_back(LP{Vec{},Key{0,0,0}}); - } -}; - - -// Lagrange points on line geometries -template -class LagrangePointSetBuilder -{ - static constexpr int dim = 1; - using LP = LagrangePoint; - using Vec = typename LP::Vector; - using Key = LocalKey; - -public: - template - void operator()(GeometryType gt, int order, Points& points) const - { - assert(gt.isLine()); - - // Vertex nodes - points.push_back(LP{Vec{0.0}, Key{0,dim,0}}); - points.push_back(LP{Vec{1.0}, Key{1,dim,0}}); - - if (order > 1) { - // Inner nodes - Vec p{0.0}; - for (unsigned int i = 0; i < order-1; i++) - { - p[0] += 1.0 / order; - points.push_back(LP{p,Key{0,dim-1,i}}); - } - } - } -}; - - -// Lagrange points on 2d geometries -template -class LagrangePointSetBuilder -{ - static constexpr int dim = 2; - using LP = LagrangePoint; - using Vec = typename LP::Vector; - using Key = LocalKey; - - friend class LagrangePointSetBuilder; - -public: - template - void operator()(GeometryType gt, int order, Points& points) const - { - std::size_t nPoints = numLagrangePoints(gt.id(), dim, order); - - if (gt.isTriangle()) - buildTriangle(nPoints, order, points); - else if (gt.isQuadrilateral()) - buildQuad(nPoints, order, points); - else { - DUNE_THROW(Dune::NotImplemented, - "Lagrange points not yet implemented for this GeometryType."); - } - - assert(points.size() == nPoints); - } - -private: - // Construct the point set in a triangle element. - // Loop from the outside to the inside - template - void buildTriangle (std::size_t nPoints, int order, Points& points) const + namespace Gmsh4 { - points.reserve(nPoints); - - const int nVertexDOFs = 3; - const int nEdgeDOFs = 3 * (order-1); - - static const unsigned int vertexPerm[3] = {0,1,2}; - static const unsigned int edgePerm[3] = {0,2,1}; - - auto calcKey = [&](int p) -> Key + namespace Impl { - if (p < nVertexDOFs) { - return Key{vertexPerm[p], dim, 0}; - } - else if (p < nVertexDOFs+nEdgeDOFs) { - unsigned int entityIndex = (p - nVertexDOFs) / (order-1); - unsigned int index = (p - nVertexDOFs) % (order-1); - return Key{edgePerm[entityIndex], dim-1, index}; - } - else { - unsigned int index = p - (nVertexDOFs + nEdgeDOFs); - return Key{0, dim-2, index}; - } - }; - - std::array bindex; - - double order_d = double(order); - for (std::size_t p = 0; p < nPoints; ++p) { - barycentricIndex(p, bindex, order); - Vec point{bindex[0] / order_d, bindex[1] / order_d}; - points.push_back(LP{point, calcKey(p)}); + // forward declaration + template + class LagrangePointSetBuilder; } - } - // "Barycentric index" is a triplet of integers, each running from 0 to - // . It is the index of a point on the triangle in barycentric - // coordinates. - static void barycentricIndex (int index, std::array& bindex, int order) - { - int max = order; - int min = 0; - - // scope into the correct triangle - while (index != 0 && index >= 3 * order) - { - index -= 3 * order; - max -= 2; - min++; - order -= 3; - } - // vertex DOFs - if (index < 3) + /// \brief A set of lagrange points compatible with the numbering of VTK and Gmsh + /** + * \tparam K Field-type for the coordinates + * \tparam dim Dimension of the coordinates + **/ + template + class LagrangePointSet + : public EmptyPointSet { - bindex[index] = bindex[(index + 1) % 3] = min; - bindex[(index + 2) % 3] = max; - } - // edge DOFs - else - { - index -= 3; - int dim = index / (order - 1); - int offset = (index - dim * (order - 1)); - bindex[(dim + 1) % 3] = min; - bindex[(dim + 2) % 3] = (max - 1) - offset; - bindex[dim] = (min + 1) + offset; - } - } - + using Super = EmptyPointSet; - // Construct the point set in the quad element - // 1. build equispaced points with index tuple (i,j) - // 2. map index tuple to DOF index and LocalKey - template - void buildQuad(std::size_t nPoints, int order, Points& points) const - { - points.resize(nPoints); - - std::array orders{order, order}; - std::array nodes{Vec{0., 0.}, Vec{1., 0.}, Vec{1., 1.}, Vec{0., 1.}}; - - for (int n = 0; n <= orders[1]; ++n) { - for (int m = 0; m <= orders[0]; ++m) { - // int idx = pointIndexFromIJ(m,n,orders); - - const double r = double(m) / orders[0]; - const double s = double(n) / orders[1]; - Vec p = (1.0 - r) * (nodes[3] * s + nodes[0] * (1.0 - s)) + - r * (nodes[2] * s + nodes[1] * (1.0 - s)); + public: + static const unsigned int dimension = dim; - auto [idx,key] = calcQuadKey(m,n,orders); - points[idx] = LP{p, key}; - // points[idx] = LP{p, calcQuadKey(n,m,orders)}; + LagrangePointSet (std::size_t order) + : Super(order) + { + assert(order > 0); } - } - } - - // Obtain the VTK DOF index of the node (i,j) in the quad element - // and construct a LocalKey - static std::pair calcQuadKey (int i, int j, std::array order) - { - const bool ibdy = (i == 0 || i == order[0]); - const bool jbdy = (j == 0 || j == order[1]); - const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); // How many boundaries do we lie on at once? - int dof = 0; - unsigned int entityIndex = 0; - unsigned int index = 0; - - if (nbdy == 2) // Vertex DOF - { - dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)); - entityIndex = (j ? (i ? 3 : 2) : (i ? 1 : 0)); - return std::make_pair(dof,Key{entityIndex, dim, 0}); - } - - int offset = 4; - if (nbdy == 1) // Edge DOF - { - if (!ibdy) { - dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset; - entityIndex = j ? 3 : 2; - index = i-1; - } - else if (!jbdy) { - dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset; - entityIndex = i ? 1 : 0; - index = j-1; + /// Fill the lagrange points for the given geometry type + void build (GeometryType gt) + { + assert(gt.dim() == dimension); + builder_(gt, order(), points_); } - return std::make_pair(dof, Key{entityIndex, dim-1, index}); - } - offset += 2 * (order[0]-1 + order[1]-1); - - // nbdy == 0: Face DOF - dof = offset + (i - 1) + (order[0]-1) * ((j - 1)); - Key innerKey = LagrangePointSetBuilder::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second; - return std::make_pair(dof, Key{0, dim-2, innerKey.index()}); - } -}; - - -// Lagrange points on 3d geometries -template -class LagrangePointSetBuilder -{ - static constexpr int dim = 3; - using LP = LagrangePoint; - using Vec = typename LP::Vector; - using Key = LocalKey; - -public: - template - void operator() (GeometryType gt, unsigned int order, Points& points) const - { - std::size_t nPoints = numLagrangePoints(gt.id(), dim, order); - - if (gt.isTetrahedron()) - buildTetra(nPoints, order, points); - else if (gt.isHexahedron()) - buildHex(nPoints, order, points); - else { - DUNE_THROW(Dune::NotImplemented, - "Lagrange points not yet implemented for this GeometryType."); - } + /// Fill the lagrange points for the given topology type `Topology` + template + bool build () + { + build(GeometryType(Topology{})); + return true; + } - assert(points.size() == nPoints); - } + /// Returns whether the point set support the given topology type `Topology` and can + /// generate point for the given order. + template + static bool supports (std::size_t order) + { + return true; + } -private: - // Construct the point set in the tetrahedron element - // 1. construct barycentric (index) coordinates - // 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index - template - void buildTetra (std::size_t nPoints, int order, Points& points) const - { - points.reserve(nPoints); + using Super::order; - const int nVertexDOFs = 4; - const int nEdgeDOFs = 6 * (order-1); - const int nFaceDOFs = 4 * (order-1)*(order-2)/2; + private: + using Super::points_; + Impl::LagrangePointSetBuilder builder_; + }; - static const unsigned int vertexPerm[4] = {0,1,2,3}; - static const unsigned int edgePerm[6] = {0,2,1,3,4,5}; - static const unsigned int facePerm[4] = {1,2,0,3}; - auto calcKey = [&](int p) -> Key + namespace Impl { - if (p < nVertexDOFs) { - return Key{vertexPerm[p], dim, 0}; - } - else if (p < nVertexDOFs+nEdgeDOFs) { - unsigned int entityIndex = (p - nVertexDOFs) / (order-1); - unsigned int index = (p - nVertexDOFs) % (order-1); - return Key{edgePerm[entityIndex], dim-1, index}; - } - else if (p < nVertexDOFs+nEdgeDOFs+nFaceDOFs) { - unsigned int index = (p - (nVertexDOFs + nEdgeDOFs)) % ((order-1)*(order-2)/2); - unsigned int entityIndex = (p - (nVertexDOFs + nEdgeDOFs)) / ((order-1)*(order-2)/2); - return Key{facePerm[entityIndex], dim-2, index}; - } - else { - unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs); - return Key{0, dim-3, index}; - } - }; + // Build for lagrange point sets in different dimensions + // Specialized for dim=1,2,3 + template + class LagrangePointSetBuilder + { + public: + template + void operator()(GeometryType, unsigned int, Points& points) const + { + DUNE_THROW(Dune::NotImplemented, + "Lagrange points not yet implemented for this GeometryType."); + } + }; - std::array bindex; - double order_d = double(order); - for (std::size_t p = 0; p < nPoints; ++p) { - barycentricIndex(p, bindex, order); - Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d}; - points.push_back(LP{point, calcKey(p)}); - } - } + // Lagrange points on point geometries + template + class LagrangePointSetBuilder + { + static constexpr int dim = 0; + using LP = LagrangePoint; + using Vec = typename LP::Vector; + using Key = LocalKey; + + public: + template + void operator()(GeometryType gt, int /*order*/, Points& points) const; + }; - // "Barycentric index" is a set of 4 integers, each running from 0 to - // . It is the index of a point in the tetrahedron in barycentric - // coordinates. - static void barycentricIndex (std::size_t p, std::array& bindex, int order) - { - const int nVertexDOFs = 4; - const int nEdgeDOFs = 6 * (order-1); - static const int edgeVertices[6][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}}; - static const int linearVertices[4][4] = {{0,0,0,1}, {1,0,0,0}, {0,1,0,0}, {0,0,1,0}}; - static const int vertexMaxCoords[4] = {3,0,1,2}; - static const int faceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}}; - static const int faceMinCoord[4] = {1,3,0,2}; + // Lagrange points on line geometries + template + class LagrangePointSetBuilder + { + static constexpr int dim = 1; + using LP = LagrangePoint; + using Vec = typename LP::Vector; + using Key = LocalKey; - int max = order; - int min = 0; + public: + template + void operator()(GeometryType gt, int order, Points& points) const; + }; - // scope into the correct tetra - while (p >= 2 * (order * order + 1) && p != 0 && order > 3) - { - p -= 2 * (order * order + 1); - max -= 3; - min++; - order -= 4; - } - // vertex DOFs - if (p < nVertexDOFs) - { - for (int coord = 0; coord < 4; ++coord) - bindex[coord] = (coord == vertexMaxCoords[p] ? max : min); - } - // edge DOFs - else if (p < nVertexDOFs+nEdgeDOFs) - { - int edgeId = (p - nVertexDOFs) / (order-1); - int vertexId = (p - nVertexDOFs) % (order-1); - for (int coord = 0; coord < 4; ++coord) + // Lagrange points on 2d geometries + template + class LagrangePointSetBuilder { - bindex[coord] = min + - (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) + - linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId)); - } - } - // face DOFs - else - { - int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2); - int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2); - - std::array projectedBIndex; - if (order == 3) - projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0; - else - LagrangePointSetBuilder::barycentricIndex(vertexId, projectedBIndex, order-3); + static constexpr int dim = 2; + using LP = LagrangePoint; + using Vec = typename LP::Vector; + using Key = LocalKey; - for (int i = 0; i < 3; i++) - bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]); + friend class LagrangePointSetBuilder; - bindex[faceMinCoord[faceId]] = min; - } - } - -private: - // Construct the point set in the heyhedral element - // 1. build equispaced points with index tuple (i,j,k) - // 2. map index tuple to DOF index and LocalKey - template - void buildHex (std::size_t nPoints, int order, Points& points) const - { - points.resize(nPoints); - - std::array orders{order, order, order}; - std::array nodes{Vec{0., 0., 0.}, Vec{1., 0., 0.}, Vec{1., 1., 0.}, Vec{0., 1., 0.}, - Vec{0., 0., 1.}, Vec{1., 0., 1.}, Vec{1., 1., 1.}, Vec{0., 1., 1.}}; - - for (int p = 0; p <= orders[2]; ++p) { - for (int n = 0; n <= orders[1]; ++n) { - for (int m = 0; m <= orders[0]; ++m) { - const double r = double(m) / orders[0]; - const double s = double(n) / orders[1]; - const double t = double(p) / orders[2]; - Vec point = (1.0-r) * ((nodes[3] * (1.0-t) + nodes[7] * t) * s + (nodes[0] * (1.0-t) + nodes[4] * t) * (1.0-s)) + - r * ((nodes[2] * (1.0-t) + nodes[6] * t) * s + (nodes[1] * (1.0-t) + nodes[5] * t) * (1.0-s)); - - auto [idx,key] = calcHexKey(m,n,p,orders); - points[idx] = LP{point, key}; - } - } - } - } + public: + template + void operator()(GeometryType gt, int order, Points& points) const; - // Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element - static std::pair calcHexKey (int i, int j, int k, std::array order) - { - const bool ibdy = (i == 0 || i == order[0]); - const bool jbdy = (j == 0 || j == order[1]); - const bool kbdy = (k == 0 || k == order[2]); - const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); // How many boundaries do we lie on at once? + private: // implementation details - int dof = 0; - unsigned int entityIndex = 0; - unsigned int index = 0; + // Construct the point set in a triangle element. + // Loop from the outside to the inside + template + void buildTriangle (std::size_t nPoints, int order, Points& points) const; - if (nbdy == 3) // Vertex DOF - { - dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0); - entityIndex = (i ? 1 : 0) + (j ? 2 : 0) + (k ? 4 : 0); - return std::make_pair(dof, Key{entityIndex, dim, 0}); - } + // "Barycentric index" is a triplet of integers, each running from 0 to + // . It is the index of a point on the triangle in barycentric + // coordinates. + static void barycentricIndex (int index, std::array& bindex, int order); - int offset = 8; - if (nbdy == 2) // Edge DOF - { - entityIndex = (k ? 8 : 4); - if (!ibdy) - { // On i axis - dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset; - index = i; - entityIndex += (i ? 1 : 0); - } - else if (!jbdy) - { // On j axis - dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset; - index = j; - entityIndex += (j ? 3 : 2); - } - else - { // !kbdy, On k axis - offset += 4 * (order[0]-1) + 4 * (order[1]-1); - dof = (k - 1) + (order[2]-1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset; - index = k; - entityIndex = (i ? 1 : 0) + (j ? 2 : 0); - } - return std::make_pair(dof, Key{entityIndex, dim-1, index}); - } + // Construct the point set in the quad element + // 1. build equispaced points with index tuple (i,j) + // 2. map index tuple to DOF index and LocalKey + template + void buildQuad(std::size_t nPoints, int order, Points& points) const; - offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1); - if (nbdy == 1) // Face DOF - { - Key faceKey; - if (ibdy) // On i-normal face - { - dof = (j - 1) + ((order[1]-1) * (k - 1)) + (i ? (order[1]-1) * (order[2]-1) : 0) + offset; - entityIndex = (i ? 1 : 0); - faceKey = LagrangePointSetBuilder::calcQuadKey(j-1,k-1,{order[1]-2, order[2]-2}).second; - } - else { - offset += 2 * (order[1] - 1) * (order[2] - 1); - if (jbdy) // On j-normal face - { - dof = (i - 1) + ((order[0]-1) * (k - 1)) + (j ? (order[2]-1) * (order[0]-1) : 0) + offset; - entityIndex = (j ? 3 : 2); - faceKey = LagrangePointSetBuilder::calcQuadKey(i-1,k-1,{order[0]-2, order[2]-2}).second; - } - else - { // kbdy, On k-normal face - offset += 2 * (order[2]-1) * (order[0]-1); - dof = (i - 1) + ((order[0]-1) * (j - 1)) + (k ? (order[0]-1) * (order[1]-1) : 0) + offset; - entityIndex = (k ? 5 : 4); - faceKey = LagrangePointSetBuilder::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second; - } - } - return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()}); - } + // Obtain the VTK DOF index of the node (i,j) in the quad element + // and construct a LocalKey + static std::pair calcQuadKey (int i, int j, std::array order); + }; - // nbdy == 0: Body DOF - offset += 2 * ((order[1]-1) * (order[2]-1) + (order[2]-1) * (order[0]-1) + (order[0]-1) * (order[1]-1)); - dof = offset + (i - 1) + (order[0]-1) * ((j - 1) + (order[1]-1) * ((k - 1))); - Key innerKey = LagrangePointSetBuilder::calcHexKey(i-1,j-1,k-1,{order[0]-2, order[1]-2, order[2]-2}).second; - return std::make_pair(dof, Key{0, dim-3, innerKey.index()}); - } -}; -}}} // end namespace Dune::Gmsh4::Impl + // Lagrange points on 3d geometries + template + class LagrangePointSetBuilder + { + static constexpr int dim = 3; + using LP = LagrangePoint; + using Vec = typename LP::Vector; + using Key = LocalKey; + + public: + template + void operator() (GeometryType gt, unsigned int order, Points& points) const; + + private: // implementation details + + // Construct the point set in the tetrahedron element + // 1. construct barycentric (index) coordinates + // 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index + template + void buildTetra (std::size_t nPoints, int order, Points& points) const; + + // "Barycentric index" is a set of 4 integers, each running from 0 to + // . It is the index of a point in the tetrahedron in barycentric + // coordinates. + static void barycentricIndex (std::size_t p, std::array& bindex, int order); + + // Construct the point set in the heyhedral element + // 1. build equispaced points with index tuple (i,j,k) + // 2. map index tuple to DOF index and LocalKey + template + void buildHex (std::size_t nPoints, int order, Points& points) const; + + // Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element + static std::pair calcHexKey (int i, int j, int k, std::array order); + }; + + } // end namespace Impl + } // end namespace Gmsh4 +} // end namespace Dune + +#include diff --git a/dune/gmsh4/utility/lagrangepoints.impl.hh b/dune/gmsh4/utility/lagrangepoints.impl.hh new file mode 100644 index 0000000..a66e9cf --- /dev/null +++ b/dune/gmsh4/utility/lagrangepoints.impl.hh @@ -0,0 +1,507 @@ +#pragma once + +#include +#include + +#include +#include +#include + +namespace Dune { +namespace Gmsh4 { +namespace Impl { + +/** + * The implementation of the point set builder is directly derived from VTK. + * Modification are a change in data-types and naming scheme. Additionally + * a LocalKey is created to reflect the concept of a Dune PointSet. + * + * Included is the license of the BSD 3-clause License included in the VTK + * source code from 2020/04/13 in commit b90dad558ce28f6d321420e4a6b17e23f5227a1c + * of git repository https://gitlab.kitware.com/vtk/vtk. + * + Program: Visualization Toolkit + Module: Copyright.txt + + Copyright (c) 1993-2015 Ken Martin, Will Schroeder, Bill Lorensen + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names + of any contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS'' + AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + **/ + + +template + template +void LagrangePointSetBuilder::operator() (GeometryType gt, int /*order*/, Points& points) const +{ + assert(gt.isVertex()); + points.push_back(LP{Vec{},Key{0,0,0}}); +} + + +template + template +void LagrangePointSetBuilder::operator() (GeometryType gt, int order, Points& points) const +{ + assert(gt.isLine()); + + // Vertex nodes + points.push_back(LP{Vec{0.0}, Key{0,dim,0}}); + points.push_back(LP{Vec{1.0}, Key{1,dim,0}}); + + if (order > 1) { + // Inner nodes + Vec p{0.0}; + for (unsigned int i = 0; i < order-1; i++) + { + p[0] += 1.0 / order; + points.push_back(LP{p,Key{0,dim-1,i}}); + } + } +} + + +template + template +void LagrangePointSetBuilder::operator() (GeometryType gt, int order, Points& points) const +{ + std::size_t nPoints = numLagrangePoints(gt.id(), dim, order); + + if (gt.isTriangle()) + buildTriangle(nPoints, order, points); + else if (gt.isQuadrilateral()) + buildQuad(nPoints, order, points); + else { + DUNE_THROW(Dune::NotImplemented, + "Lagrange points not yet implemented for this GeometryType."); + } + + assert(points.size() == nPoints); +} + + +// Construct the point set in a triangle element. +// Loop from the outside to the inside +template + template +void LagrangePointSetBuilder::buildTriangle (std::size_t nPoints, int order, Points& points) const +{ + points.reserve(nPoints); + + const int nVertexDOFs = 3; + const int nEdgeDOFs = 3 * (order-1); + + static const unsigned int vertexPerm[3] = {0,1,2}; + static const unsigned int edgePerm[3] = {0,2,1}; + + auto calcKey = [&](int p) -> Key + { + if (p < nVertexDOFs) { + return Key{vertexPerm[p], dim, 0}; + } + else if (p < nVertexDOFs+nEdgeDOFs) { + unsigned int entityIndex = (p - nVertexDOFs) / (order-1); + unsigned int index = (p - nVertexDOFs) % (order-1); + return Key{edgePerm[entityIndex], dim-1, index}; + } + else { + unsigned int index = p - (nVertexDOFs + nEdgeDOFs); + return Key{0, dim-2, index}; + } + }; + + std::array bindex; + + double order_d = double(order); + for (std::size_t p = 0; p < nPoints; ++p) { + barycentricIndex(p, bindex, order); + Vec point{bindex[0] / order_d, bindex[1] / order_d}; + points.push_back(LP{point, calcKey(p)}); + } +} + + +// "Barycentric index" is a triplet of integers, each running from 0 to +// . It is the index of a point on the triangle in barycentric +// coordinates. +template +void LagrangePointSetBuilder::barycentricIndex (int index, std::array& bindex, int order) +{ + int max = order; + int min = 0; + + // scope into the correct triangle + while (index != 0 && index >= 3 * order) + { + index -= 3 * order; + max -= 2; + min++; + order -= 3; + } + + // vertex DOFs + if (index < 3) + { + bindex[index] = bindex[(index + 1) % 3] = min; + bindex[(index + 2) % 3] = max; + } + // edge DOFs + else + { + index -= 3; + int dim = index / (order - 1); + int offset = (index - dim * (order - 1)); + bindex[(dim + 1) % 3] = min; + bindex[(dim + 2) % 3] = (max - 1) - offset; + bindex[dim] = (min + 1) + offset; + } +} + + +// Construct the point set in the quad element +// 1. build equispaced points with index tuple (i,j) +// 2. map index tuple to DOF index and LocalKey +template + template +void LagrangePointSetBuilder::buildQuad(std::size_t nPoints, int order, Points& points) const +{ + points.resize(nPoints); + + std::array orders{order, order}; + std::array nodes{Vec{0., 0.}, Vec{1., 0.}, Vec{1., 1.}, Vec{0., 1.}}; + + for (int n = 0; n <= orders[1]; ++n) { + for (int m = 0; m <= orders[0]; ++m) { + // int idx = pointIndexFromIJ(m,n,orders); + + const double r = double(m) / orders[0]; + const double s = double(n) / orders[1]; + Vec p = (1.0 - r) * (nodes[3] * s + nodes[0] * (1.0 - s)) + + r * (nodes[2] * s + nodes[1] * (1.0 - s)); + + auto [idx,key] = calcQuadKey(m,n,orders); + points[idx] = LP{p, key}; + // points[idx] = LP{p, calcQuadKey(n,m,orders)}; + } + } +} + + +// Obtain the VTK DOF index of the node (i,j) in the quad element +// and construct a LocalKey +template +std::pair::Key> +LagrangePointSetBuilder::calcQuadKey (int i, int j, std::array order) +{ + const bool ibdy = (i == 0 || i == order[0]); + const bool jbdy = (j == 0 || j == order[1]); + const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); // How many boundaries do we lie on at once? + + int dof = 0; + unsigned int entityIndex = 0; + unsigned int index = 0; + + if (nbdy == 2) // Vertex DOF + { + dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)); + entityIndex = (j ? (i ? 3 : 2) : (i ? 1 : 0)); + return std::make_pair(dof,Key{entityIndex, dim, 0}); + } + + int offset = 4; + if (nbdy == 1) // Edge DOF + { + if (!ibdy) { + dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset; + entityIndex = j ? 3 : 2; + index = i-1; + } + else if (!jbdy) { + dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset; + entityIndex = i ? 1 : 0; + index = j-1; + } + return std::make_pair(dof, Key{entityIndex, dim-1, index}); + } + + offset += 2 * (order[0]-1 + order[1]-1); + + // nbdy == 0: Face DOF + dof = offset + (i - 1) + (order[0]-1) * ((j - 1)); + Key innerKey = LagrangePointSetBuilder::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second; + return std::make_pair(dof, Key{0, dim-2, innerKey.index()}); +} + + +// Lagrange points on 3d geometries +template + template +void LagrangePointSetBuilder::operator() (GeometryType gt, unsigned int order, Points& points) const +{ + std::size_t nPoints = numLagrangePoints(gt.id(), dim, order); + + if (gt.isTetrahedron()) + buildTetra(nPoints, order, points); + else if (gt.isHexahedron()) + buildHex(nPoints, order, points); + else { + DUNE_THROW(Dune::NotImplemented, + "Lagrange points not yet implemented for this GeometryType."); + } + + assert(points.size() == nPoints); +} + + +// Construct the point set in the tetrahedron element +// 1. construct barycentric (index) coordinates +// 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index +template + template +void LagrangePointSetBuilder::buildTetra (std::size_t nPoints, int order, Points& points) const +{ + points.reserve(nPoints); + + const int nVertexDOFs = 4; + const int nEdgeDOFs = 6 * (order-1); + const int nFaceDOFs = 4 * (order-1)*(order-2)/2; + + static const unsigned int vertexPerm[4] = {0,1,2,3}; + static const unsigned int edgePerm[6] = {0,2,1,3,4,5}; + static const unsigned int facePerm[4] = {1,2,0,3}; + + auto calcKey = [&](int p) -> Key + { + if (p < nVertexDOFs) { + return Key{vertexPerm[p], dim, 0}; + } + else if (p < nVertexDOFs+nEdgeDOFs) { + unsigned int entityIndex = (p - nVertexDOFs) / (order-1); + unsigned int index = (p - nVertexDOFs) % (order-1); + return Key{edgePerm[entityIndex], dim-1, index}; + } + else if (p < nVertexDOFs+nEdgeDOFs+nFaceDOFs) { + unsigned int index = (p - (nVertexDOFs + nEdgeDOFs)) % ((order-1)*(order-2)/2); + unsigned int entityIndex = (p - (nVertexDOFs + nEdgeDOFs)) / ((order-1)*(order-2)/2); + return Key{facePerm[entityIndex], dim-2, index}; + } + else { + unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs); + return Key{0, dim-3, index}; + } + }; + + std::array bindex; + + double order_d = double(order); + for (std::size_t p = 0; p < nPoints; ++p) { + barycentricIndex(p, bindex, order); + Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d}; + points.push_back(LP{point, calcKey(p)}); + } +} + + +// "Barycentric index" is a set of 4 integers, each running from 0 to +// . It is the index of a point in the tetrahedron in barycentric +// coordinates. +template +void LagrangePointSetBuilder::barycentricIndex (std::size_t p, std::array& bindex, int order) +{ + const int nVertexDOFs = 4; + const int nEdgeDOFs = 6 * (order-1); + + static const int edgeVertices[6][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}}; + static const int linearVertices[4][4] = {{0,0,0,1}, {1,0,0,0}, {0,1,0,0}, {0,0,1,0}}; + static const int vertexMaxCoords[4] = {3,0,1,2}; + static const int faceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}}; + static const int faceMinCoord[4] = {1,3,0,2}; + + int max = order; + int min = 0; + + // scope into the correct tetra + while (p >= 2 * (order * order + 1) && p != 0 && order > 3) + { + p -= 2 * (order * order + 1); + max -= 3; + min++; + order -= 4; + } + + // vertex DOFs + if (p < nVertexDOFs) + { + for (int coord = 0; coord < 4; ++coord) + bindex[coord] = (coord == vertexMaxCoords[p] ? max : min); + } + // edge DOFs + else if (p < nVertexDOFs+nEdgeDOFs) + { + int edgeId = (p - nVertexDOFs) / (order-1); + int vertexId = (p - nVertexDOFs) % (order-1); + for (int coord = 0; coord < 4; ++coord) + { + bindex[coord] = min + + (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) + + linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId)); + } + } + // face DOFs + else + { + int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2); + int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2); + + std::array projectedBIndex; + if (order == 3) + projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0; + else + LagrangePointSetBuilder::barycentricIndex(vertexId, projectedBIndex, order-3); + + for (int i = 0; i < 3; i++) + bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]); + + bindex[faceMinCoord[faceId]] = min; + } +} + + +// Construct the point set in the hexahedral element +// 1. build equispaced points with index tuple (i,j,k) +// 2. map index tuple to DOF index and LocalKey +template + template +void LagrangePointSetBuilder::buildHex (std::size_t nPoints, int order, Points& points) const +{ + points.resize(nPoints); + + std::array orders{order, order, order}; + std::array nodes{Vec{0., 0., 0.}, Vec{1., 0., 0.}, Vec{1., 1., 0.}, Vec{0., 1., 0.}, + Vec{0., 0., 1.}, Vec{1., 0., 1.}, Vec{1., 1., 1.}, Vec{0., 1., 1.}}; + + for (int p = 0; p <= orders[2]; ++p) { + for (int n = 0; n <= orders[1]; ++n) { + for (int m = 0; m <= orders[0]; ++m) { + const double r = double(m) / orders[0]; + const double s = double(n) / orders[1]; + const double t = double(p) / orders[2]; + Vec point = (1.0-r) * ((nodes[3] * (1.0-t) + nodes[7] * t) * s + (nodes[0] * (1.0-t) + nodes[4] * t) * (1.0-s)) + + r * ((nodes[2] * (1.0-t) + nodes[6] * t) * s + (nodes[1] * (1.0-t) + nodes[5] * t) * (1.0-s)); + + auto [idx,key] = calcHexKey(m,n,p,orders); + points[idx] = LP{point, key}; + } + } + } +} + + +// Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element +template +std::pair::Key> +LagrangePointSetBuilder::calcHexKey (int i, int j, int k, std::array order) +{ + const bool ibdy = (i == 0 || i == order[0]); + const bool jbdy = (j == 0 || j == order[1]); + const bool kbdy = (k == 0 || k == order[2]); + const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); // How many boundaries do we lie on at once? + + int dof = 0; + unsigned int entityIndex = 0; + unsigned int index = 0; + + if (nbdy == 3) // Vertex DOF + { + dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0); + entityIndex = (i ? 1 : 0) + (j ? 2 : 0) + (k ? 4 : 0); + return std::make_pair(dof, Key{entityIndex, dim, 0}); + } + + int offset = 8; + if (nbdy == 2) // Edge DOF + { + entityIndex = (k ? 8 : 4); + if (!ibdy) + { // On i axis + dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset; + index = i; + entityIndex += (i ? 1 : 0); + } + else if (!jbdy) + { // On j axis + dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset; + index = j; + entityIndex += (j ? 3 : 2); + } + else + { // !kbdy, On k axis + offset += 4 * (order[0]-1) + 4 * (order[1]-1); + dof = (k - 1) + (order[2]-1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset; + index = k; + entityIndex = (i ? 1 : 0) + (j ? 2 : 0); + } + return std::make_pair(dof, Key{entityIndex, dim-1, index}); + } + + offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1); + if (nbdy == 1) // Face DOF + { + Key faceKey; + if (ibdy) // On i-normal face + { + dof = (j - 1) + ((order[1]-1) * (k - 1)) + (i ? (order[1]-1) * (order[2]-1) : 0) + offset; + entityIndex = (i ? 1 : 0); + faceKey = LagrangePointSetBuilder::calcQuadKey(j-1,k-1,{order[1]-2, order[2]-2}).second; + } + else { + offset += 2 * (order[1] - 1) * (order[2] - 1); + if (jbdy) // On j-normal face + { + dof = (i - 1) + ((order[0]-1) * (k - 1)) + (j ? (order[2]-1) * (order[0]-1) : 0) + offset; + entityIndex = (j ? 3 : 2); + faceKey = LagrangePointSetBuilder::calcQuadKey(i-1,k-1,{order[0]-2, order[2]-2}).second; + } + else + { // kbdy, On k-normal face + offset += 2 * (order[2]-1) * (order[0]-1); + dof = (i - 1) + ((order[0]-1) * (j - 1)) + (k ? (order[0]-1) * (order[1]-1) : 0) + offset; + entityIndex = (k ? 5 : 4); + faceKey = LagrangePointSetBuilder::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second; + } + } + return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()}); + } + + // nbdy == 0: Body DOF + offset += 2 * ((order[1]-1) * (order[2]-1) + (order[2]-1) * (order[0]-1) + (order[0]-1) * (order[1]-1)); + dof = offset + (i - 1) + (order[0]-1) * ((j - 1) + (order[1]-1) * ((k - 1))); + Key innerKey = LagrangePointSetBuilder::calcHexKey(i-1,j-1,k-1,{order[0]-2, order[1]-2, order[2]-2}).second; + return std::make_pair(dof, Key{0, dim-3, innerKey.index()}); +} + +}}} // end namespace Dune::Gmsh4::Impl diff --git a/src/ALUGridTests.cc b/src/ALUGridTests.cc index 1d7e920..27da190 100644 --- a/src/ALUGridTests.cc +++ b/src/ALUGridTests.cc @@ -10,12 +10,12 @@ #include -#include +#include #include -#include +#include #include -#include +#include #include #include @@ -57,7 +57,7 @@ int main(int argc, char** argv) { //test 02: triangle-pair: gmsh4 -> vtk (flat, order 3) std::cout << "\ntriangle-pair: gmsh4 -> vtk (flat, order 3)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/triangles_3d_order3_subset.msh"); + GRID_PATH "/triangles_3d_order3_subset.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 3); @@ -71,7 +71,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/triangles_3d_order3_subset.msh"); + reader.read(GRID_PATH "/triangles_3d_order3_subset.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -88,7 +88,7 @@ int main(int argc, char** argv) { //test 04: sphere: gmsh4 -> vtk (flat, order 4) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, ascii)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4.msh"); + GRID_PATH "/sphere_order4.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -102,7 +102,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4.msh"); + reader.read(GRID_PATH "/sphere_order4.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -119,7 +119,7 @@ int main(int argc, char** argv) { //test 06: sphere: gmsh4 -> vtk (flat, order 1, binary) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 1, binary)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order1_binary.msh"); + GRID_PATH "/sphere_order1_binary.msh"); auto& grid = *gridPtr; VtkUnstructuredGridWriter vtkWriter(grid.leafGridView(), Vtk::ASCII); @@ -129,7 +129,7 @@ int main(int argc, char** argv) { //test 07: sphere: gmsh4 -> vtk (flat, order 4, binary) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, binary)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4_binary.msh"); + GRID_PATH "/sphere_order4_binary.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -143,7 +143,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4_binary.msh"); + reader.read(GRID_PATH "/sphere_order4_binary.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -160,7 +160,7 @@ int main(int argc, char** argv) { //test 09: sphere: gmsh4 -> vtk (flat, order 1, quadrilaterals) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 1, ascii, quadrilaterals)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order1_quad.msh"); + GRID_PATH "/sphere_order1_quad.msh"); auto& grid = *gridPtr; VtkUnstructuredGridWriter vtkWriter(grid.leafGridView(), Vtk::ASCII); @@ -170,7 +170,7 @@ int main(int argc, char** argv) { //test 10: sphere: gmsh4 -> vtk (flat, order 4, quadrilaterals) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, ascii, quadrilaterals)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4_quad.msh"); + GRID_PATH "/sphere_order4_quad.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -184,7 +184,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4_quad.msh"); + reader.read(GRID_PATH "/sphere_order4_quad.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -201,7 +201,7 @@ int main(int argc, char** argv) { //test 12: sphere: gmsh4 -> vtk (flat, order 1, binary, quadrilaterals) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 1, binary, quadrilaterals)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order1_quad_binary.msh"); + GRID_PATH "/sphere_order1_quad_binary.msh"); auto& grid = *gridPtr; VtkUnstructuredGridWriter vtkWriter(grid.leafGridView(), Vtk::ASCII); @@ -211,7 +211,7 @@ int main(int argc, char** argv) { //test 13: sphere: gmsh4 -> vtk (flat, order 4, binary, quadrilaterals) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, binary, quadrilaterals)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4_quad_binary.msh"); + GRID_PATH "/sphere_order4_quad_binary.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -225,7 +225,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4_quad_binary.msh"); + reader.read(GRID_PATH "/sphere_order4_quad_binary.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2370f99..52a0bd3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,17 +1,24 @@ +set(GRID_PATH "${CMAKE_CURRENT_SOURCE_DIR}/meshes/") + if(dune-alugrid_FOUND) add_executable("gmsh4reader" gmsh4reader.cc) target_link_dune_default_libraries("gmsh4reader") target_link_libraries(gmsh4reader dunegmsh4) + target_compile_definitions(gmsh4reader PRIVATE "GRID_PATH=\"${GRID_PATH}\"") endif() -if(dune-alugrid_FOUND AND dune-vtk_FOUND AND dune-curvedsurfacegrid_FOUND) - add_executable("ALUGridTests" ALUGridTests.cc) - target_link_dune_default_libraries("ALUGridTests") - target_link_libraries(ALUGridTests dunegmsh4) -endif() +dune_add_test(SOURCES ALUGridTests.cc + LINK_LIBRARIES dunegmsh4 + COMPILE_DEFINITIONS "GRID_PATH=\"${GRID_PATH}\"" + CMAKE_GUARD + dune-alugrid_FOUND + dune-vtk_FOUND + dune-curvedsurfacegrid_FOUND) -if(dune-foamgrid_FOUND AND dune-vtk_FOUND AND dune-curvedsurfacegrid_FOUND) - add_executable("FoamGridTests" FoamGridTests.cc) - target_link_dune_default_libraries("FoamGridTests") - target_link_libraries(FoamGridTests dunegmsh4) -endif() +dune_add_test(SOURCES FoamGridTests.cc + LINK_LIBRARIES dunegmsh4 + COMPILE_DEFINITIONS "GRID_PATH=\"${GRID_PATH}\"" + CMAKE_GUARD + dune-foamgrid_FOUND + dune-vtk_FOUND + dune-curvedsurfacegrid_FOUND) diff --git a/src/FoamGridTests.cc b/src/FoamGridTests.cc index 1841bc0..88cd669 100644 --- a/src/FoamGridTests.cc +++ b/src/FoamGridTests.cc @@ -10,12 +10,12 @@ #include -#include +#include #include -#include +#include #include -#include +#include #include #include @@ -56,7 +56,7 @@ int main(int argc, char** argv) { //test 02: triangle-pair: gmsh4 -> vtk (flat, order 3) std::cout << "\ntriangle-pair: gmsh4 -> vtk (flat, order 3)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/triangles_3d_order3_subset.msh"); + GRID_PATH "/triangles_3d_order3_subset.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 3); @@ -70,7 +70,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/triangles_3d_order3_subset.msh"); + reader.read(GRID_PATH "/triangles_3d_order3_subset.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -87,7 +87,7 @@ int main(int argc, char** argv) { //test 04: sphere: gmsh4 -> vtk (flat, order 4) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, ascii)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4.msh"); + GRID_PATH "/sphere_order4.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -101,7 +101,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4.msh"); + reader.read(GRID_PATH "/sphere_order4.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; @@ -118,7 +118,7 @@ int main(int argc, char** argv) { //test 06: sphere: gmsh4 -> vtk (flat, order 1, binary) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 1, binary)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order1_binary.msh"); + GRID_PATH "/sphere_order1_binary.msh"); auto& grid = *gridPtr; VtkUnstructuredGridWriter vtkWriter(grid.leafGridView(), Vtk::ASCII); @@ -128,7 +128,7 @@ int main(int argc, char** argv) { //test 07: sphere: gmsh4 -> vtk (flat, order 4, binary) std::cout << "\nsphere: gmsh4 -> vtk (flat, order 4, binary)" << std::endl; std::unique_ptr gridPtr = Gmsh4Reader::createGridFromFile( - "dune-gmsh4/src/meshes/sphere_order4_binary.msh"); + GRID_PATH "/sphere_order4_binary.msh"); auto& grid = *gridPtr; Vtk::LagrangeDataCollector dataCollector(grid.leafGridView(), 4); @@ -142,7 +142,7 @@ int main(int argc, char** argv) Gmsh4::LagrangeGridCreator creator(factory); Gmsh4Reader reader(creator); - reader.read("dune-gmsh4/src/meshes/sphere_order4_binary.msh"); + reader.read(GRID_PATH "/sphere_order4_binary.msh"); std::unique_ptr gridPtr = factory.createGrid(); auto& grid = *gridPtr; diff --git a/src/gmsh4reader.cc b/src/gmsh4reader.cc index 0c6c176..3d97f92 100644 --- a/src/gmsh4reader.cc +++ b/src/gmsh4reader.cc @@ -16,7 +16,7 @@ #include #include -#include +#include #include #include @@ -34,7 +34,7 @@ int main(int argc, char** argv) using GridType = Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>; using GridView = typename GridType::LeafGridView; { - auto gridPtr = Gmsh4Reader::createGridFromFile("src/meshes/square.msh"); + auto gridPtr = Gmsh4Reader::createGridFromFile(GRID_PATH "/square.msh"); auto& grid = *gridPtr; VTKWriter vtkWriter(grid.leafGridView()); @@ -42,7 +42,7 @@ int main(int argc, char** argv) } { - auto gridPtr = Gmsh4Reader::createGridFromFile("src/meshes/square_part1_topology.pro"); + auto gridPtr = Gmsh4Reader::createGridFromFile(GRID_PATH "/square_part1_topology.pro"); auto& grid = *gridPtr; VTKWriter vtkWriter(grid.leafGridView()); @@ -50,7 +50,7 @@ int main(int argc, char** argv) } { - auto gridPtr = Gmsh4Reader::createGridFromFile("src/meshes/square_part2.msh"); + auto gridPtr = Gmsh4Reader::createGridFromFile(GRID_PATH "/square_part2.msh"); auto& grid = *gridPtr; VTKWriter vtkWriter(grid.leafGridView()); -- GitLab From cd6e6cb5064847c540f95e9ffd3d463323cd05d6 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 16:35:48 +0200 Subject: [PATCH 2/8] added gitlab-ci configuration --- .gitlab-ci.yml | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 .gitlab-ci.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000..613ec5d --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,17 @@ +--- + +include: + - "https://gitlab.mn.tu-dresden.de/amdis/ci-config/raw/master/config/common/master.yml" + - "https://gitlab.mn.tu-dresden.de/amdis/ci-config/raw/master/jobs/common/master.yml" + +before_script: + - . /duneci/bin/duneci-init-job + - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git + - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git + - duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git + - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git + - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git + - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git + - duneci-install-module https://gitlab.mn.tu-dresden.de/iwr/dune-curvedsurfacegrid.git -- GitLab From 1ce74eb2c259b0bd5eb558853263168562726489 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 16:49:03 +0200 Subject: [PATCH 3/8] update gitlab-ci configuration --- .gitlab-ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 613ec5d..d60b451 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -8,6 +8,7 @@ before_script: - . /duneci/bin/duneci-init-job - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git - duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git + - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git - duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git -- GitLab From e871279061c4f1da2fb5d731fc07b9847d3fb942 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 17:54:52 +0200 Subject: [PATCH 4/8] correct bug of overloaded by temporary --- dune/gmsh4/gridcreators/lagrangegridcreator.hh | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/dune/gmsh4/gridcreators/lagrangegridcreator.hh b/dune/gmsh4/gridcreators/lagrangegridcreator.hh index fa1b05e..5072c41 100644 --- a/dune/gmsh4/gridcreators/lagrangegridcreator.hh +++ b/dune/gmsh4/gridcreators/lagrangegridcreator.hh @@ -245,7 +245,11 @@ namespace Dune return LocalFunction{gridCreator}; } - friend LocalFunction localFunction (LagrangeGridCreator&& gridCreator) = delete; + friend LocalFunction localFunction (LagrangeGridCreator&& gridCreator) + { + DUNE_THROW(Dune::Exception, "Cannot pass temporary LagrangeGridCreator to localFunction(). Pass an lvalue-reference instead."); + return LocalFunction{gridCreator}; + } struct EntitySet { -- GitLab From b0bbd8e96f1d062f8450a5e42e4cd24c498c6a31 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 18:10:13 +0200 Subject: [PATCH 5/8] change alugrid branch to 2.7 --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index d60b451..e430e02 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -11,7 +11,7 @@ before_script: - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git - - duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git + - duneci-install-module --branch releases/2.7 https://gitlab.dune-project.org/extensions/dune-alugrid.git - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git -- GitLab From 3d1aae1ba904f7e4dba3005914f4c1e7042c8d91 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 18:30:46 +0200 Subject: [PATCH 6/8] fix ALUGrid bug --- src/ALUGridTests.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/ALUGridTests.cc b/src/ALUGridTests.cc index 27da190..d38ad41 100644 --- a/src/ALUGridTests.cc +++ b/src/ALUGridTests.cc @@ -5,6 +5,10 @@ # include "config.h" #endif +// workaround for bug in ALUGrid +#include +namespace Dune { using std::shared_ptr; } + #include // An initializer of MPI #include // We use exceptions -- GitLab From 2bea84402dbd9620cb8dc3f0c648533c7516d09d Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Mon, 17 Aug 2020 18:42:12 +0200 Subject: [PATCH 7/8] deactivate ALUGrid in ci tests --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e430e02..a9f4c3e 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -11,7 +11,7 @@ before_script: - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git - - duneci-install-module --branch releases/2.7 https://gitlab.dune-project.org/extensions/dune-alugrid.git + - # duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git -- GitLab From 26a6e23505c54f92b47bef6b84ebb11bf3c54262 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 6 Sep 2020 14:54:22 +0200 Subject: [PATCH 8/8] remove comment from gitlabci yml --- .gitlab-ci.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a9f4c3e..e0d302a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -11,7 +11,6 @@ before_script: - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git - - # duneci-install-module https://gitlab.dune-project.org/extensions/dune-alugrid.git - duneci-install-module https://gitlab.dune-project.org/extensions/dune-foamgrid.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-vtk.git - duneci-install-module https://gitlab.mn.tu-dresden.de/spraetor/dune-curvedgeometry.git -- GitLab