Skip to content
GitLab
Commit 75e11993 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

make naming consistent with dune-vtk and transform examples into tests

parent eaed36c2
Expand all Hide whitespace changes
Inline Side-by-side
dune/gmsh4/filereader.hh
View file @ 75e11993
......@@ -11,7 +11,6 @@ namespace Dune
{
namespace Gmsh4
{
template <class Grid, class FilerReaderImp>
class FileReader
{
......
dune/gmsh4/reader.hh → dune/gmsh4/gmsh4reader.hh
View file @ 75e11993
......@@ -274,4 +274,4 @@ namespace Dune
} // end namespace Dune
#include "reader.impl.hh"
#include "gmsh4reader.impl.hh"
dune/gmsh4/reader.impl.hh → dune/gmsh4/gmsh4reader.impl.hh
View file @ 75e11993
......@@ -10,8 +10,8 @@
#include <dune/gmsh4/utility/filesystem.hh>
#include <dune/gmsh4/utility/string.hh>
//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);
......
dune/gmsh4/gridcreators/common.hh
View file @ 75e11993
......@@ -6,7 +6,6 @@ namespace Dune
{
namespace Gmsh4
{
template <class Factory, class... Args>
using HasInsertVertex = decltype( std::declval<Factory>().insertVertex(std::declval<Args>()...) );
......
dune/gmsh4/gridcreators/continuousgridcreator.hh
View file @ 75e11993
......@@ -34,8 +34,8 @@ namespace Dune
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
{
vertexMap_.resize(nodeTagRange.second - nodeTagRange.first + 1);
vertexShift_ = nodeTagRange.first;
......@@ -55,9 +55,9 @@ namespace Dune
template <class ElementAttributes, class BoundaryEntities>
void insertElementsImpl (std::size_t numElements,
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
{
std::vector<unsigned int> 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<double,Grid::dimension-1>(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<double,Grid::dimension>(cell.type());
connectivity.resize(refElem.size(Grid::dimension));
......
dune/gmsh4/gridcreators/derivedgridcreator.hh
View file @ 75e11993
......@@ -31,17 +31,17 @@ namespace Dune
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
{
gridCreator_.insertVertices(numNodes, nodeTagRange, entityBlocks);
}
template <class ElementAttributes, class BoundaryEntities>
void insertElementsImpl (std::size_t numElements,
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
{
gridCreator_.insertElements(numElements, elementTagRange, entityBlocks, boundaryEntities);
}
......
dune/gmsh4/gridcreators/discontinuousgridcreator.hh
View file @ 75e11993
......@@ -47,55 +47,19 @@ namespace Dune
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> 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 <class ElementAttributes, class BoundaryEntities>
void insertElementsImpl (std::size_t numElements,
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> 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<unsigned int> 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<unsigned int> 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:
......
dune/gmsh4/gridcreators/lagrangegridcreator.hh
View file @ 75e11993
......@@ -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 <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> 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 <class ElementAttributes, class BoundaryEntities>
void insertElementsImpl (std::size_t numElements,
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
{
const int dim = GridType::dimension;
std::vector<unsigned int> 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<double, dim - 1>(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<double, dim>(type);
connectivity.resize(refElem.size(dim));
......@@ -409,5 +393,4 @@ namespace Dune
};
} // end namespace Gmsh4
} // end namespace Dune
dune/gmsh4/gridcreators/parallelgridcreator.hh
View file @ 75e11993
......@@ -35,8 +35,8 @@ namespace Dune
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
{
std::cout << "WARNING! Dune::Gmsh4::ParallelGridCreator is unfinished!" << std::endl;
GlobalCoordinate p;
......
dune/gmsh4/gridcreators/serialgridcreator.hh
View file @ 75e11993
......@@ -29,48 +29,24 @@ namespace Dune
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> const& entityBlocks)
std::pair<std::size_t,std::size_t> nodeTagRange,
std::vector<NodeAttributes> 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 <class ElementAttributes, class BoundaryEntities>
void insertElementsImpl (std::size_t numElements,
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> const& entityBlocks,
BoundaryEntities const& boundaryEntities)
std::pair<std::size_t,std::size_t> elementTagRange,
std::vector<ElementAttributes> 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<std::string> const& pieces)
{
std::cout << "WARNING! Dune::Gmsh4::SerialGridCreator is unfinished!" << std::endl;
// if (this->comm().rank() == 0) {
// Gmsh4Reader<Grid, Self> pieceReader(*this);
// for (std::string const& piece : pieces) {
// pieceReader.readFromFile(piece, false);
// pieceReader.createGrid(false);
// }
// DiscontinuousGridCreator<Grid> creator(this->factory());
// creator.insertVertices(points_, {});
// creator.insertElements(types_, offsets_, connectivity_);
// }
DUNE_THROW(Dune::NotImplemented, "Dune::Gmsh4::SerialGridCreator is unfinished!");
}
private:
......
dune/gmsh4/types.cc
View file @ 75e11993
......@@ -6,7 +6,6 @@ namespace Dune
{
namespace Gmsh4
{
GeometryType to_geometry (int elementType)
{
switch (elementType) {
......
dune/gmsh4/utility/CMakeLists.txt
View file @ 75e11993
......@@ -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)
dune/gmsh4/utility/filesystem.cc
View file @ 75e11993
......@@ -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);
......
dune/gmsh4/utility/lagrangepoints.hh
View file @ 75e11993
This diff is collapsed.
dune/gmsh4/utility/lagrangepoints.impl.hh 0 → 100644
View file @ 75e11993
#pragma once
#include <cassert>
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
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 <class K>
template <class Points>
void LagrangePointSetBuilder<K,0>::operator() (GeometryType gt, int /*order*/, Points& points) const
{
assert(gt.isVertex());
points.push_back(LP{Vec{},Key{0,0,0}});
}
template <class K>
template <class Points>
void LagrangePointSetBuilder<K,1>::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 <class K>
template <class Points>
void LagrangePointSetBuilder<K,2>::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 <class K>
template <class Points>
void LagrangePointSetBuilder<K,2>::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<int,3> 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
// <Order>. It is the index of a point on the triangle in barycentric
// coordinates.
template <class K>
void LagrangePointSetBuilder<K,2>::barycentricIndex (int index, std::array<int,3>& 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 <class K>
template <class Points>
void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints, int order, Points& points) const
{
points.resize(nPoints);
std::array<int,2> orders{order, order};
std::array<Vec,4> 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 <class K>
std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
LagrangePointSetBuilder<K,2>::calcQuadKey (int i, int j, std::array<int,2> 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<K,2>::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 K>