tetrahedroncell.cpp

利用C

// Copyright (C) 2006-2008 Anders Logg.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Johan Hoffman 2006.
// Modified by Garth N. Wells 2006.
// Modified by Kristian Oelgaard 2006.
//
// First added:  2006-06-05
// Last changed: 2008-06-20

#include <algorithm>
#include <dolfin/log/dolfin_log.h>
#include "Cell.h"
#include "MeshEditor.h"
#include "MeshGeometry.h"
#include "Facet.h"
#include "TetrahedronCell.h"
#include "Vertex.h"
#include "GeometricPredicates.h"

using namespace dolfin;

//-----------------------------------------------------------------------------
dolfin::uint TetrahedronCell::dim() const
{
  return 3;
}
//-----------------------------------------------------------------------------
dolfin::uint TetrahedronCell::numEntities(uint dim) const
{
  switch ( dim )
  {
  case 0:
    return 4; // vertices
  case 1:
    return 6; // edges
  case 2:
    return 4; // faces
  case 3:
    return 1; // cells
  default:
    error("Illegal topological dimension %d for tetrahedron.", dim);
  }

  return 0;
}
//-----------------------------------------------------------------------------
dolfin::uint TetrahedronCell::numVertices(uint dim) const
{
  switch ( dim )
  {
  case 0:
    return 1; // vertices
  case 1:
    return 2; // edges
  case 2:
    return 3; // faces
  case 3:
    return 4; // cells
  default:
    error("Illegal topological dimension %d for tetrahedron.", dim);
  }

  return 0;
}
//-----------------------------------------------------------------------------
dolfin::uint TetrahedronCell::orientation(const Cell& cell) const
{
  // This is a trick to be allowed to initialize mesh entities from cell
  Cell& c = const_cast<Cell&>(cell);

  Vertex v0(c.mesh(), c.entities(0)[0]);
  Vertex v1(c.mesh(), c.entities(0)[1]);
  Vertex v2(c.mesh(), c.entities(0)[2]);
  Vertex v3(c.mesh(), c.entities(0)[3]);

  Point p01 = v1.point() - v0.point();
  Point p02 = v2.point() - v0.point();
  Point p03 = v3.point() - v0.point();

  Point n = p01.cross(p02);

  return ( n.dot(p03) < 0.0 ? 1 : 0 );
}
//-----------------------------------------------------------------------------
void TetrahedronCell::createEntities(uint** e, uint dim, const uint* v) const
{
  // We only need to know how to create edges and faces
  switch ( dim )
  {
  case 1:
    // Create the six edges
    e[0][0] = v[2]; e[0][1] = v[3];
    e[1][0] = v[1]; e[1][1] = v[3];
    e[2][0] = v[1]; e[2][1] = v[2];
    e[3][0] = v[0]; e[3][1] = v[3];
    e[4][0] = v[0]; e[4][1] = v[2];
    e[5][0] = v[0]; e[5][1] = v[1];
    break;
  case 2:
    // Create the four faces
    e[0][0] = v[1]; e[0][1] = v[2]; e[0][2] = v[3];
    e[1][0] = v[0]; e[1][1] = v[2]; e[1][2] = v[3];
    e[2][0] = v[0]; e[2][1] = v[1]; e[2][2] = v[3];
    e[3][0] = v[0]; e[3][1] = v[1]; e[3][2] = v[2];
    break;
  default:
    error("Don't know how to create entities of topological dimension %d.", dim);
  }
}
//-----------------------------------------------------------------------------
void TetrahedronCell::orderEntities(Cell& cell) const
{
  // Sort i - j for i > j: 1 - 0, 2 - 0, 2 - 1, 3 - 0, 3 - 1, 3 - 2

  // Get mesh topology
  MeshTopology& topology = cell.mesh().topology();

  // Sort local vertices on edges in ascending order, connectivity 1 - 0
  if ( topology(1, 0).size() > 0 )
  {
    dolfin_assert(topology(3, 1).size() > 0);

    // Get edges
    uint* cell_edges = cell.entities(1);

    // Sort vertices on each edge
    for (uint i = 0; i < 6; i++)
    {
      uint* edge_vertices = topology(1, 0)(cell_edges[i]);
      std::sort(edge_vertices, edge_vertices + 2);
    }
  }

  // Sort local vertices on facets in ascending order, connectivity 2 - 0
  if ( topology(2, 0).size() > 0 )
  {
    dolfin_assert(topology(3, 2).size() > 0);

    // Get facets
    uint* cell_facets = cell.entities(2);

    // Sort vertices on each facet
    for (uint i = 0; i < 4; i++)
    {
      uint* facet_vertices = topology(2, 0)(cell_facets[i]);
      std::sort(facet_vertices, facet_vertices + 3);
    }
  }

  // Sort local edges on local facets after non-incident vertex, connectivity 2 - 1
  if ( topology(2, 1).size() > 0 )
  {
    dolfin_assert(topology(3, 2).size() > 0);
    dolfin_assert(topology(2, 0).size() > 0);
    dolfin_assert(topology(1, 0).size() > 0);

    // Get facet numbers
    uint* cell_facets = cell.entities(2);

    // Loop over facets on cell
    for (uint i = 0; i < 4; i++)
    {
      // For each facet number get the global vertex numbers
      uint* facet_vertices = topology(2, 0)(cell_facets[i]);

      // For each facet number get the global edge number
      uint* cell_edges = topology(2, 1)(cell_facets[i]);

      // Loop over vertices on facet
      uint m = 0;
      for (uint j = 0; j < 3; j++)
      {
        // Loop edges on facet
        for (uint k(m); k < 3; k++)
        {
          // For each edge number get the global vertex numbers
          uint* edge_vertices = topology(1, 0)(cell_edges[k]);

          // Check if the jth vertex of facet i is non-incident on edge k
          if (!std::count(edge_vertices, edge_vertices+2, facet_vertices[j]))
          {
            // Swap facet numbers
            uint tmp = cell_edges[m];
            cell_edges[m] = cell_edges[k];
            cell_edges[k] = tmp;
            m++;
            break;
          }
        }
      }
    }
  }

  // Sort local vertices on cell in ascending order, connectivity 3 - 0
  if ( topology(3, 0).size() > 0 )
  {
    uint* cell_vertices = cell.entities(0);
    std::sort(cell_vertices, cell_vertices+4);
  }

  // Sort local edges on cell after non-incident vertex tuble, connectivity 3-1
  if ( topology(3, 1).size() > 0 )
  {
    dolfin_assert(topology(1, 0).size() > 0);

    // Get cell vertices and edge numbers
    uint* cell_vertices = cell.entities(0);
    uint* cell_edges = cell.entities(1);

    // Loop two vertices on cell as a lexicographical tuple
    // (i, j): (0,1) (0,2) (0,3) (1,2) (1,3) (2,3)
    uint m = 0;
    for (uint i = 0; i < 3; i++)
    {
      for (uint j = i+1; j < 4; j++)
      {
        // Loop edge numbers
        for (uint k = m; k < 6; k++)
        {
          // Get local vertices on edge
          uint* edge_vertices = topology(1, 0)(cell_edges[k]);

          // Check if the ith and jth vertex of the cell are non-incident on edge k
          if (!std::count(edge_vertices, edge_vertices+2, cell_vertices[i]) && \
              !std::count(edge_vertices, edge_vertices+2, cell_vertices[j]) )
          {
            // Swap edge numbers
            uint tmp = cell_edges[m];
            cell_edges[m] = cell_edges[k];
            cell_edges[k] = tmp;
            m++;
            break;
          }
        }
      }
    }
  }

  // Sort local facets on cell after non-incident vertex, connectivity 3 - 2
  if ( topology(3, 2).size() > 0 )
  {
    dolfin_assert(topology(2, 0).size() > 0);

    // Get cell vertices and facet numbers
    uint* cell_vertices = cell.entities(0);
    uint* cell_facets = cell.entities(2);

    // Loop vertices on cell
    for (uint i = 0; i < 4; i++)
    {
      // Loop facets on cell
      for (uint j = i; j < 4; j++)
      {
        uint* facet_vertices = topology(2, 0)(cell_facets[j]);

        // Check if the ith vertex of the cell is non-incident on facet j
        if (!std::count(facet_vertices, facet_vertices+3, cell_vertices[i]))
        {
          // Swap facet numbers
          uint tmp = cell_facets[i];
          cell_facets[i] = cell_facets[j];
          cell_facets[j] = tmp;
          break;
        }
      }
    }
  }
}
//-----------------------------------------------------------------------------
void TetrahedronCell::refineCell(Cell& cell, MeshEditor& editor,
			     uint& current_cell) const
{
  // Get vertices and edges
  const uint* v = cell.entities(0);
  const uint* e = cell.entities(1);
  dolfin_assert(v);
  dolfin_assert(e);

  // Get offset for new vertex indices
  const uint offset = cell.mesh().numVertices();

  // Compute indices for the ten new vertices
  const uint v0 = v[0];
  const uint v1 = v[1];
  const uint v2 = v[2];
  const uint v3 = v[3];
  const uint e0 = offset + e[findEdge(0, cell)];
  const uint e1 = offset + e[findEdge(1, cell)];
  const uint e2 = offset + e[findEdge(2, cell)];
  const uint e3 = offset + e[findEdge(3, cell)];
  const uint e4 = offset + e[findEdge(4, cell)];
  const uint e5 = offset + e[findEdge(5, cell)];

  // Regular refinement: 8 new cells
  editor.addCell(current_cell++, v0, e3, e4, e5);
  editor.addCell(current_cell++, v1, e1, e2, e5);
  editor.addCell(current_cell++, v2, e0, e2, e4);
  editor.addCell(current_cell++, v3, e0, e1, e3);
  editor.addCell(current_cell++, e0, e1, e2, e5);
  editor.addCell(current_cell++, e0, e1, e3, e5);
  editor.addCell(current_cell++, e0, e2, e4, e5);
  editor.addCell(current_cell++, e0, e3, e4, e5);
}
//-----------------------------------------------------------------------------
void TetrahedronCell::refineCellIrregular(Cell& cell, MeshEditor& editor, 
				      uint& current_cell, uint refinement_rule, 
				      uint* marked_edges) const
{
  warning("Not implemented yet.");