dune-dec issueshttps://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues2017-09-26T15:07:46Zhttps://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/11Refactoring of half-edge grid implementation2017-09-26T15:07:46ZGhost UserRefactoring of half-edge grid implementationExtract common functionality and algorithms from half-edge grid implementations to define a more abstract interface and reduce the complexity of `indexSet`.
Add adapter for specific implementation if needed.Extract common functionality and algorithms from half-edge grid implementations to define a more abstract interface and reduce the complexity of `indexSet`.
Add adapter for specific implementation if needed.https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/10Add tests for the correctness of the halfe-edge grid's internal structure.2017-09-26T13:08:50ZGhost UserAdd tests for the correctness of the halfe-edge grid's internal structure.There are some basic relations of half edges all implementations of a half edge grid should comply to.
If not all basic operations are implemented, some test should probably be deactivated for that specific implementation.
### Correctn...There are some basic relations of half edges all implementations of a half edge grid should comply to.
If not all basic operations are implemented, some test should probably be deactivated for that specific implementation.
### Correctness of a single basic operation
- [ ] `next()` should be **cyclic** and **consistent**:
- [ ] For interior half edges the cycle length should be equal to
the number of edges of the corresponding element type.
All elements of such a cycle should be interior half edges belonging to the same face.
- [ ] For half edges on the boundary all elements of the cycle should lie on the boundary, too.
- [ ] `previous()` should be **cyclic** and **consistent**:
- [ ] For interior half edges the cycle length should be equal to the number of edges of the corresponding element type.
All elements of such a cycle should be interior half edges belonging to the same face.
- [ ] For half edges on the boundary all elements of the cycle should lie on the boundary, too.
- [ ] `opposite()` should be an **involution** and **consistent**:
- [ ] `opposite() of opposite() == id`
- [ ] A half edge and its `opposite()` belong to the same edge.
### Correctness of combinations of basic operations
- [ ] `next()` and `previous() are the **inverse** of each other:
- [ ] If both `next()` and `previous()` are defined for a half edge:
`next() of previous() == id`
`previous() of next() == id`
- [ ] An iteration around a vertex should **terminate** and be **consistent**:
- [ ] Whether the vertex belonging to a half edge is its target or source
and independently of the specific combination of basic operations used to get the next half edge during this iteration, the iteration should return to the first half edge.
- [ ] The number of steps for an interior vertex should be at least 3.
The number of steps for an vertex on the boundary should be at least 2.
- [ ] All half edges in such an iteration should have the same vertex.https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/9Bug in init() of default halfedge implementation2017-10-24T06:11:04ZGhost UserBug in init() of default halfedge implementationSomething is going wrong with the `init()` function.
The cause is most likely located in `reorient()`.
## Simple example
### Halfedge grid created from simple grid
```
using WorldVector = Dune::FieldVector<float_type, dow>;
WorldVect...Something is going wrong with the `init()` function.
The cause is most likely located in `reorient()`.
## Simple example
### Halfedge grid created from simple grid
```
using WorldVector = Dune::FieldVector<float_type, dow>;
WorldVector x0 = { 0.0, 0.0, 0.0 };
WorldVector x1 = { 1.0, 0.0, 0.0 };
WorldVector x2 = { 0.0, 1.0, 0.0 };
WorldVector x3 = { 0.0, 0.0, 1.0 };
std::vector<std::vector<index_type>> cells = { {0, 1, 2}
, {0, 3, 1}
, {1, 3, 2}
// , {0, 2, 3} }; // same orientation
, {0, 3, 2} }; // different orientation
std::vector<WorldVector> coords = { x0, x1, x2, x3 };
GridBase base_grid {cells, coords};
Grid grid {base_grid};
```
The grid `grid` has correct incidences of edges to vertices/faces and of faces to vertices/edges.
However, the incidences of vertices to edges and faces are wrong:
### Query incidences of vertices
```
for (auto const& vertex : vertices(grid.gridView())) {
msg_("vertex[", vertex.index(), "]:");
msg_("\nvertices: ");
for (auto const& v : vertices(vertex))
msg_(v.index(), " ");
msg_("\nedges: ");
for (auto const& e : edges(vertex))
msg_(e.index(), " ");
msg_("\nfaces: ");
for (auto const& f : faces(vertex))
msg_(f.index(), " ");
msg();
}
```
### Output
```
vertex[0]:
vertices: 0
edges: 1 3 4 5 1 2 4 0
faces: 3 1 2 3 0 2 1 0
vertex[1]:
vertices: 1
edges: 2 4 0 1 3 4 5 1
faces: 2 1 0 3 1 2 3 0
vertex[2]:
vertices: 2
edges: 5 1 2 4 0 1 3 4
faces: 3 0 2 1 0 3 1 2
vertex[3]:
vertices: 3
edges: 3 5 2 0
faces: 3 2 0 1
```https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/8benchmark and optimize Half-Edge implementation2017-09-26T12:17:04ZPraetorius, Simonbenchmark and optimize Half-Edge implementationThe current implementation of half-edge datastructure is based on index-redirections, i.e.
```c++
struct HalfEdge
{
using IndexType = unsigned long int;
// an index representing unset entities or neighbours
static constexpr...The current implementation of half-edge datastructure is based on index-redirections, i.e.
```c++
struct HalfEdge
{
using IndexType = unsigned long int;
// an index representing unset entities or neighbours
static constexpr IndexType invalid = IndexType(-1);
// half edges
IndexType next_ = invalid;
IndexType previous_ = invalid;
IndexType opposite_ = invalid;
// entities
IndexType vertex_ = invalid;
IndexType edge_ = invalid;
IndexType face_ = invalid;
}
```
Can we reduce some of the fields? Is an implementation based on pointers a better choice?
```c++
struct HalfEdge {
HalfEdge* next;
HalfEdge* previous;
HalfEdge* opposite;
EntitySeed vertex;
EntitySeed edge;
EntitySeed face;
};
```
where `EntitySeed` might be a `Dune::EntitySeed` type to create real entities from.
What about rewriting an array of structs to a struct of arrays?
A slightly different datastructure is explained in Tyler *J. Alumbaugh and Xiangmin Jiao*: "Compact Array-Based Mesh Data Structures"
**TODO:** This needs to be benchmarked. Several implementations should be compared.https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/5Weighted triangulation2017-09-26T12:17:04ZPraetorius, SimonWeighted triangulationRecently we have found a way to allow non-circumcentric well-centered grid. The idea is based on a vertex weight that shift the circumcenter towards the barycenter. There is a prove of concept available in the `examples/` directory.
Wh...Recently we have found a way to allow non-circumcentric well-centered grid. The idea is based on a vertex weight that shift the circumcenter towards the barycenter. There is a prove of concept available in the `examples/` directory.
What needs to be implemented is a way of assigning new coordinates to the cached circumcenters and to update the cached dual volumes. Also, there must be found a way to implement the shift of the circumcenter
```
c* = c^old - 1/2 * grad(w)
```
with `grad(w)` constant on the element, and `grad` the surface gradient in the tangent plane of the entity.
There is an explicit formula for an edge `e_ij = (x_i, x_j)`:
```
c*_ij = x_i + (|e_ij|^2 + w_i - w_j) / (2 |e_ij|^2) * e_ij
```
Additionally we need a test for well-centeredness of a weighted triangulation and a local weight calculation (maybe based on a gauss-seidel relaxation of a laplace-beltrami equation, see "*Weighted Triangulations for Geometry Processing*" (de Geos, F. and Memari, P. and Mullen, P. and Desbrun, M.)https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/4Linear-Algebra backend2017-09-26T12:17:04ZPraetorius, SimonLinear-Algebra backendThe is a backend for Eigen3 implemented and for a minimal own data-structure. This needs to be extended for more libraries, e.g. MTL4 and PETSc. For PETSc already a basic implementation is available (for the old dec-pde library) that pro...The is a backend for Eigen3 implemented and for a minimal own data-structure. This needs to be extended for more libraries, e.g. MTL4 and PETSc. For PETSc already a basic implementation is available (for the old dec-pde library) that probably can easily be adapted.
There are currently just a few requirements for the backend. Basically it must be possible to run the CG or BCGS method with the implemented backend on the iterative solvers in `linear_algebra/krylov`.
The backend interface should be documented somewhere.https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/3Musical isomorphisms2017-09-26T12:17:04ZPraetorius, SimonMusical isomorphismsIn the file `FlatSharp.hpp` and `MusicalIsomorphism.hpp` there are some flat and sharp operations already implemented. This needs to be extended to any source data, i.e. for vector-fields given on vertices, edges or cells.
Additionally...In the file `FlatSharp.hpp` and `MusicalIsomorphism.hpp` there are some flat and sharp operations already implemented. This needs to be extended to any source data, i.e. for vector-fields given on vertices, edges or cells.
Additionally an implementation for analytically given vector fields must be provided.
Also the inverse transformation sharp should be implemented, so that the resulting vector can be defined on any entity.
Add tests for the implementation and a comparison for the various methods.https://gitlab.math.tu-dresden.de/spraetor/dune-dec/-/issues/1Circumcenter calculation2017-09-26T12:17:04ZPraetorius, SimonCircumcenter calculationCurrently the calculation of circumcenters is based on a least-squares problem:
For all edges of the entity, formulate pairwise condition
```
|c - x_i|^2 = |x - x_j|^2, were (x_i, x_j) forms an edge
```
with `c = lambda * x`, i.e. `c...Currently the calculation of circumcenters is based on a least-squares problem:
For all edges of the entity, formulate pairwise condition
```
|c - x_i|^2 = |x - x_j|^2, were (x_i, x_j) forms an edge
```
with `c = lambda * x`, i.e. `c` given in barycentric coordinates. Additionally require `sum(lambda) = 1`.
The overdetermined system is solved with a qr-method with full pivoting, using Eigen, or row-pivoting using an own implementation. Both implementation fail for some cases, i.e. produce `inf` or `NaN`.
This needs to be improved! Find a more stable and maybe more simple way of calculating the circumcenter for an (simplcial) entity of any dimension.