refinedsimplexlocalbasis.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_REFINED_SIMPLEX_LOCALBASIS_HH
#define DUNE_REFINED_SIMPLEX_LOCALBASIS_HH
 
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/common/localbasis.hh>
 
namespace Dune
{
  template<class D, int dim>
  class RefinedSimplexLocalBasis
  {
  protected:
    RefinedSimplexLocalBasis()
    {
      DUNE_THROW(Dune::NotImplemented,"RefinedSimplexLocalBasis not implemented for dim > 3.");
    }
  };
 
  template<class D>
  class RefinedSimplexLocalBasis<D,1>
  {
  protected:
 
    RefinedSimplexLocalBasis() {}
 
    static int getSubElement(const FieldVector<D,1>& global)
    {
      if (global[0] <= 0.5)
        return 0;
      else if (global[0] <= 1.0)
        return 1;
 
      DUNE_THROW(InvalidStateException, "no subelement defined");
    }
 
    static void getSubElement(const FieldVector<D,1>& global,
                              int& subElement,
                              FieldVector<D,1>& local)
    {
      if (global[0] <= 0.5) {
        subElement = 0;
        local[0] = 2.0 * global[0];
        return;
      }
 
      subElement = 1;
      local[0] = 2.0 * global[0] - 1.0;
    }
 
  };
 
 
  template<class D>
  class RefinedSimplexLocalBasis<D,2>
  {
  protected:
 
    RefinedSimplexLocalBasis() {}
 
    static int getSubElement(const FieldVector<D,2>& global)
    {
      if (global[0] + global[1] <= 0.5)
        return 0;
      else if (global[0] >= 0.5)
        return 1;
      else if (global[1] >= 0.5)
        return 2;
 
      return 3;
    }
 
    static void getSubElement(const FieldVector<D,2>& global,
                              int& subElement,
                              FieldVector<D,2>& local)
    {
      if (global[0] + global[1] <= 0.5) {
        subElement = 0;
        local[0] = 2*global[0];
        local[1] = 2*global[1];
        return;
      } else if (global[0] >= 0.5) {
        subElement = 1;
        local[0] = 2*global[0]-1;
        local[1] = 2*global[1];
        return;
      } else if (global[1] >= 0.5) {
        subElement = 2;
        local[0] = 2*global[0];
        local[1] = 2*global[1]-1;
        return;
      }
 
      subElement = 3;
      local[0] = -2 * global[0] + 1;
      local[1] = -2 * global[1] + 1;
 
    }
 
 
  };
 
  template<class D>
  class RefinedSimplexLocalBasis<D,3>
  {
  protected:
 
    RefinedSimplexLocalBasis() {}
 
    static int getSubElement(const FieldVector<D,3>& global)
    {
      if (global[0] + global[1] + global[2] <= 0.5)
        return 0;
      else if (global[0] >= 0.5)
        return 1;
      else if (global[1] >= 0.5)
        return 2;
      else if (global[2] >= 0.5)
        return 3;
      else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] <= 0.5))
        return 4;
      else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] <= 0.5))
        return 5;
      else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] >= 0.5))
        return 6;
      else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] >= 0.5))
        return 7;
 
      DUNE_THROW(InvalidStateException, "no subelement defined");
 
    }
    static void getSubElement(const FieldVector<D,3>& global,
                              int& subElement,
                              FieldVector<D,3>& local)
    {
      if (global[0] + global[1] + global[2] <= 0.5) {
        subElement = 0;
        local = global;
        local *= 2.0;
        return;
      } else if (global[0] >= 0.5) {
        subElement = 1;
        local = global;
        local[0] -= 0.5;
        local *= 2.0;
        return;
      } else if (global[1] >= 0.5) {
        subElement = 2;
        local = global;
        local[1] -= 0.5;
        local *= 2.0;
        return;
      } else if (global[2] >= 0.5) {
        subElement = 3;
        local = global;
        local[2] -= 0.5;
        local *= 2.0;
        return;
      } else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] <= 0.5)) {
        subElement = 4;
        local[0] = 2.0 * global[1];
        local[1] = 2.0 * (0.5 - global[0] - global[1]);
        local[2] = 2.0 * (-0.5 + global[0] + global[1] + global[2]);
        //              Dune::FieldMatrix<double,3,3> A(0.0);
        //              A[0][1] =  2.0;
        //              A[1][0] = -2.0;
        //              A[1][1] = -2.0;
        //              A[2][0] =  2.0;
        //              A[2][1] =  2.0;
        //              A[2][2] =  2.0;
        //              A.mv(global,local);
        //              local[1] += 1.0;
        //              local[2] -= 1.0;
        return;
      } else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] <= 0.5)) {
        subElement = 5;
        local[0] = 2.0 * (0.5 - global[0]);
        local[1] = 2.0 * (0.5 - global[1] - global[2]);
        local[2] = 2.0 * global[2];
        //              Dune::FieldMatrix<double,3,3> A(0.0);
        //              A[0][0] = -2.0;
        //              A[1][1] = -2.0;
        //              A[1][2] = -2.0;
        //              A[2][2] =  2.0;
        //              A.mv(global,local);
        //              local[0] += 1.0;
        //              local[1] += 1.0;
        return;
      } else if ((global[0] + global[1] <= 0.5)and (global[1] + global[2] >= 0.5)) {
        subElement = 6;
        local[0] = 2.0 * (0.5 - global[0] - global[1]);
        local[1] = 2.0 * global[0];
        local[2] = 2.0 * (-0.5 + global[1] + global[2]);
        //              Dune::FieldMatrix<double,3,3> A(0.0);
        //              A[0][0] = -2.0;
        //              A[0][1] = -2.0;
        //              A[1][0] =  2.0;
        //              A[2][1] =  2.0;
        //              A[2][2] =  2.0;
        //              A.mv(global,local);
        //              local[0] += 1.0;
        //              local[2] -= 1.0;
        return;
      } else if ((global[0] + global[1] >= 0.5)and (global[1] + global[2] >= 0.5)) {
        subElement = 7;
        local[0] = 2.0 * (-0.5 + global[1] + global[2]);
        local[1] = 2.0 * (0.5 - global[1]);
        local[2] = 2.0 * (-0.5 + global[0] + global[1]);
        //              Dune::FieldMatrix<double,3,3> A(0.0);
        //              A[0][1] =  2.0;
        //              A[0][2] =  2.0;
        //              A[1][1] = -2.0;
        //              A[2][0] =  2.0;
        //              A[2][1] =  2.0;
        //              A.mv(global,local);
        //              local[0] -= 1.0;
        //              local[1] += 1.0;
        //              local[2] -= 1.0;
        return;
      }
 
      DUNE_THROW(InvalidStateException, "no subelement defined");
 
    }
 
  };
 
 
}
 
#endif