test06_borderInTriangle.cc

#include "config.h"
#include <iostream>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/alugrid/grid.hh>
#include <dune/alugrid/dgf.hh>
#include <iostream>
#include <array>
#include <typeinfo>

/*
for every triangle checks if corners are above, below or inbetween the edges --> if corners are on different sides border within trinagle

*/

using namespace Dune;
constexpr int dim = 2; // Grid and world dimension


double distance(FieldVector<double, dim> p1, FieldVector<double, dim> p2){
    return std::abs((p1 - p2).two_norm());
}






void horizontalFlow(std::shared_ptr<Dune::ALUGrid< dim, dim, Dune::simplex, Dune::conforming>> grid, double width, FieldVector<double, dim> pos, double minSizeFactor = 0.5 ){

    //TODO check if pos is at edge of map --> special case

    typedef Dune::ALUGrid< dim, dim, Dune::simplex, Dune::conforming > Grid;
    using GridView = Grid::LeafGridView;
    const GridView gridView = grid->leafGridView();

    double epsilon = 0.001;
    double minSize =  minSizeFactor*width;
    std::cout << width << " " << minSizeFactor << std::endl;
    std::cout << minSize << std::endl;  

    double yBorderU = pos[1] + width/2;
    double yBorderD = pos[1] - width/2;

    double xRangeL = pos[0] - 0.5;
    double xRangeR = pos[0] + 0.5;
    
    
    int c=0;
    while( c<100){

        c++;
        
        for (const auto& element : elements(gridView)){

            bool between = false;
            bool upper = false;
            bool under = false;

            bool continueFor = false;
            int xOut = 0;

            for (int i = 0; i < 3; i++){
                if(element.geometry().corner(i)[0] < xRangeL || element.geometry().corner(i)[0] > xRangeR){
                    xOut ++;
                }
                if (element.geometry().corner(i)[1] > yBorderU) upper = true;
                else if (element.geometry().corner(i)[1] < yBorderD) under = true;
                else if (std::abs(element.geometry().corner(i)[1] - yBorderD) < 1e-5 ||std::abs( element.geometry().corner(i)[1] - yBorderU) < 1e-5){ //border on edge of triangle
                    grid->mark(-1, element); 
                    continueFor = true;
                    break;
                }
                else between = true;
            }

            if(xOut == 3 || continueFor || between + upper + under <= 1) continue; // all 3 corners outside x boundry/border onborder not within triangle
    

            
            double dist01 = distance(element.geometry().corner(0), element.geometry().corner(1));
            double dist02 = distance(element.geometry().corner(0), element.geometry().corner(2));
            double dist12 = distance(element.geometry().corner(1), element.geometry().corner(2));
            
            if (std::min({dist01, dist02}) < minSize) continue; //small enough

           
            Dune::FieldVector<double, 2> v;
            double yCenterHeight;

            if (dist01 == dist02){ //dist12
                v = element.geometry().corner(2) - element.geometry().corner(1);
                if (element.geometry().corner(0)[1] == element.geometry().corner(1)[1]){
                    yCenterHeight = element.geometry().corner(0)[1] + (element.geometry().corner(2)[1] - element.geometry().corner(1)[1]) / 2;
                }
                else yCenterHeight = element.geometry().corner(0)[1] + (element.geometry().corner(1)[1] - element.geometry().corner(2)[1]) / 2;
                }
            if ( dist01 > dist02){  //calc hypothenuse v; y value middle of height of triangle
                v = element.geometry().corner(0) - element.geometry().corner(1);
                if(element.geometry().corner(2)[1] == element.geometry().corner(0)[1]) {
                    yCenterHeight = element.geometry().corner(2)[1] + (element.geometry().corner(1)[1] - element.geometry().corner(0)[1]) / 2;}
                else yCenterHeight = element.geometry().corner(2)[1] + (element.geometry().corner(0)[1] - element.geometry().corner(1)[1]) / 2;
                }
            else {
                v = element.geometry().corner(0) - element.geometry().corner(2);
                if (element.geometry().corner(1)[1] == element.geometry().corner(2)[1]){
                    yCenterHeight = element.geometry().corner(1)[1] + (element.geometry().corner(0)[1] - element.geometry().corner(2)[1]) / 2;}
                else yCenterHeight = element.geometry().corner(1)[1] + (element.geometry().corner(2)[1] - element.geometry().corner(0)[1]) / 2;
                }
           

            Dune::FieldVector<double, 2> u1{1.0, 0.0};
            Dune::FieldVector<double, 2> u2{0, 1.0};
            double dotProduct1 = v[0] * u1[0] + v[1] * u1[1];
            double dotProduct2 = v[0] * u2[0] + v[1] * u2[1];
            if (std::abs(dotProduct1) < 1e-8 || std::abs(dotProduct2) < 1e-8){
                grid->mark(1, element);//hypothenuse vertical or horizontal
                continue;
            }


            double epsilon3 = 1e-2;
            if(std::abs(yBorderU - yCenterHeight) < epsilon3 || std::abs(yBorderD - yCenterHeight) < epsilon3) continue; // border close to middle of triangle (in right angle)

            grid->mark(1,element);
        }
       
        grid->preAdapt();
        grid->adapt();
        grid->postAdapt();            
    }  
}



// { grid_setup_begin }
int main(int argc, char *argv[])
{
    // Set up MPI if available
    MPIHelper::instance(argc, argv);

    typedef Dune::ALUGrid< dim, dim, Dune::simplex, Dune::conforming > Grid;
    

    // Start with a structured grid
    const std::array<unsigned, dim> n = {6, 6};
    const FieldVector<double, dim> lower = {0, 0};
    const FieldVector<double, dim> upper = {6, 6};

    std::shared_ptr<Grid> grid = StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, n);

    using GridView = Grid::LeafGridView;
    const GridView gridView = grid->leafGridView();

  

    horizontalFlow(grid, 0.1, FieldVector<double, dim> {{4, 1.5}});
    horizontalFlow(grid, 0.3, FieldVector<double, dim> {{3, 1.5}});
    

    // Write grid to file
    VTKWriter<GridView> vtkWriter(gridView);
    vtkWriter.write("refined_grid_borderInTriangle");
    
}