improvegrid.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_IMPROVEGRID_HH
#define DUNE_FUFEM_IMPROVEGRID_HH
 
#include <cmath>
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/istl/bvector.hh>
 
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/prolongboundarypatch.hh>
 
namespace Dune::Fufem::Impl {
 
    template <class T>
    Dune::FieldVector<T, 3> crossProduct(const Dune::FieldVector<T, 3>& a, const Dune::FieldVector<T, 3>& b)
    {
        Dune::FieldVector<T, 3> r;
        r[0] = a[1] * b[2] - a[2] * b[1];
        r[1] = a[2] * b[0] - a[0] * b[2];
        r[2] = a[0] * b[1] - a[1] * b[0];
        return r;
    }
 
} // end namespace Dune::Fufem::Impl
 
template <class EntityType>
double aspectRatio(const EntityType& entity)
{
    if (!entity.type().isTetrahedron())
        DUNE_THROW(Dune::NotImplemented, "aspect ratio only for simplices!");
 
    const int dim = EntityType::dimension;
 
    namespace Impl = Dune::Fufem::Impl;
    typename EntityType::Geometry geometry = entity.geometry();
 
    const Dune::FieldVector<double,dim> a = geometry.corner(0);
    const Dune::FieldVector<double,dim> b = geometry.corner(1);
    const Dune::FieldVector<double,dim> c = geometry.corner(2);
    const Dune::FieldVector<double,dim> d = geometry.corner(3);
 
    Dune::FieldVector<double,dim> ba = b - a;
    Dune::FieldVector<double,dim> ca = c - a;
    Dune::FieldVector<double,dim> da = d - a;
 
    Dune::FieldVector<double,dim> circumRadiusNumerator0 = Impl::crossProduct(ba,ca);
    circumRadiusNumerator0 *= da.two_norm2();
    Dune::FieldVector<double,dim> circumRadiusNumerator1 = Impl::crossProduct(da,ba);
    circumRadiusNumerator1 *= ca.two_norm2();
    Dune::FieldVector<double,dim> circumRadiusNumerator2 = Impl::crossProduct(ca,da);
    circumRadiusNumerator2 *= ba.two_norm2();
 
    Dune::FieldVector<double,dim> circumRadiusNumerator = circumRadiusNumerator0
        + circumRadiusNumerator1 + circumRadiusNumerator2;
 
    double circumRadius = circumRadiusNumerator.two_norm() / (2*(ba*Impl::crossProduct(ca,da)));
 
    double volume = Impl::crossProduct(ba,ca) * da / 6;
 
    double area = Impl::crossProduct(ba,ca).two_norm()/2;
    area += Impl::crossProduct(da,ba).two_norm()/2;
    area += Impl::crossProduct(da,ca).two_norm()/2;
    area += Impl::crossProduct(Dune::FieldVector<double,dim>(c-b),Dune::FieldVector<double,dim>(d-b)).two_norm()/2;
 
    double inRadius = 3*volume/area;
 
    //  May be negative due to inexact arithmetic
    return std::abs(circumRadius / inRadius);
 
}
 
template <class GridType>
int improveGridLevel(GridType& grid, double threshold, int level)
{
    const int dim = GridType::dimension;
    const auto& levelView = grid.levelGridView(level);
    const auto& indexSet = levelView.indexSet();
 
    Dune::BitSetVector<1> isNewOnThisLevel(indexSet.size(dim), false);
    Dune::BitSetVector<1> isTetMidPoint(indexSet.size(dim), false);
    std::vector<Dune::FieldVector<double, dim> > bisectionPosition(indexSet.size(dim));
 
    // ///////////////////////////////////////////////////////////////
    //   Create the whole grid boundary as a boundary patch
    // ///////////////////////////////////////////////////////////////
 
    BoundaryPatch<typename GridType::LevelGridView> boundary(levelView, true);
 
    // /////////////////////////////////////////////////////////////
    //   We need to know which vertices is new on the maxlevel and which
    //   vertices are descendants of vertices on coarser levels.
    // /////////////////////////////////////////////////////////////
 
    double edgeMidPoints[6][3] = {{0.5,0,0},
                                  {0.5,0.5,0},
                                  {0,0.5,0},
                                  {0,0,0.5},
                                  {0.5,0,0.5},
                                  {0,0.5,0.5}};
 
    for (const auto& e : elements(levelView)) {
 
        const auto& father = e.father();
 
        const auto geometryInFather = e.geometryInFather();
        for (size_t i=0; i<e.subEntities(dim); i++) {
 
            const auto& positionInFather = geometryInFather.corner(i);
 
            for (int j=0; j<6; j++) {
                if ( std::abs(positionInFather[0]-edgeMidPoints[j][0]) < 0.1
                     && std::abs(positionInFather[1]-edgeMidPoints[j][1]) < 0.1
                     && std::abs(positionInFather[2]-edgeMidPoints[j][2]) < 0.1) {
 
                    isNewOnThisLevel[indexSet.subIndex(e,i,dim)][0] = true;
                    bisectionPosition[indexSet.subIndex(e,i,dim)] = father.geometry().global(positionInFather);
                    break;
                }
            }
            // tetraeder mid points must be moved, too
            if ( std::abs(positionInFather[0]-0.25) < 0.1
                 && std::abs(positionInFather[1]-0.25) < 0.1
                 && std::abs(positionInFather[2]-0.25) < 0.1) {
 
                isNewOnThisLevel[indexSet.subIndex(e,i,dim)][0] = true;
                isTetMidPoint[indexSet.subIndex(e,i,dim)][0] = true;
                bisectionPosition[indexSet.subIndex(e,i,dim)] = father.geometry().global(positionInFather);
            }
        }
    }
 
    // ///////////////////////////////////////////////////////////////
    //   Mark the vertices of all 'bad' elements
    // ///////////////////////////////////////////////////////////////
    Dune::BitSetVector<1> badVertices(indexSet.size(dim),false);
    for (const auto& e : elements(levelView)){
 
        Dune::FieldVector<double,dim> dummy(0.5);
 
        double indicator = aspectRatio(e);
        if (indicator > threshold || e.geometry().jacobianInverseTransposed(dummy).determinant() < 0) {
 
            // This element is of bad quality.  Mark all its vertices
            for (size_t i=0; i<e.subEntities(dim); i++){
                if (boundary.containsVertex(indexSet.subIndex(e,i,dim)) || isTetMidPoint[indexSet.subIndex(e,i,dim)][0]) {
                    badVertices[indexSet.subIndex(e,i,dim)][0] = true;
                }
            }
        }
    }
 
    // Loop over all new vertices
    for (const auto& v : vertices(levelView)) {
 
        // Get the index
        int vIdx = indexSet.index(v);
 
        if (!isNewOnThisLevel[vIdx][0] || !badVertices[vIdx][0])
            continue;
 
        grid.setPosition(v, bisectionPosition[vIdx]);
    }
 
    return badVertices.count();
}
 
template <class GridType>
void improveGrid(GridType& grid, double threshold, int maxIter=10)
{
  const int dim = GridType::dimension;
  int maxLevel = grid.maxLevel();
 
  if (maxLevel==0)
      return;
 
  // Only works in 3d
  if (dim!=3)
      DUNE_THROW(Dune::GridError, "Grid improvement only for 3d grids!");
 
  // By moving some vertices to improve certain elements other elements can worsen.
  for(int i=0;i<maxIter;++i) {
 
    int nBadVertices = 0;
    // Could be, that new nodes were added on some corser level, too. So we must check all levels.
    for(int level=maxLevel; level>0; --level) {
      nBadVertices += improveGridLevel(grid, threshold, level);
    }
    std::cout << nBadVertices << " bad vertices were detected in grid improvement iteration " << i << "." << std::endl;
    if(nBadVertices==0) break;
  }
}
 
 
template <class GridView>
void printAspectRatios(const GridView& gridView)
{
    double max(0);
    double min(std::numeric_limits<double>::max());
 
    for (const auto& e : elements(gridView)) {
 
        double indicator = aspectRatio(e);
 
        max = std::max(max ,indicator);
        min = std::min(min ,indicator);
    }
    std::cout<<"Aspect ratios range from "<<min<<" -- "<<max<<std::endl;
}
 
#endif