~ubuntu-branches/ubuntu/trusty/rheolef/trusty-proposed

« back to all changes in this revision

Viewing changes to nfem/basis/P2_numbering.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
ImportĀ upstreamĀ versionĀ 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
#include "P2_numbering.h"
 
22
using namespace std;
 
23
 
 
24
/* For mopre info: Dhatt et Touzot 2eme ed p110 
 
25
 * WARNING:
 
26
 * The node numbering correspond to the basis numbering.
 
27
 * The node numbering convention is different from Dhatt et Touzot:
 
28
 * we use P. L. Georges et H. Borouchaki
 
29
 *        "Triangulation de Delaunay et maillages",
 
30
 *        Hermes, 1997
 
31
 *   - first the vertices
 
32
 *   - then the middle of the edges
 
33
 *
 
34
 * WARNING2: the transformation is still linear
 
35
 */
 
36
 
 
37
/*
 
38
 * global numbering statement: uses symbolic values (do not uses as variable !):
 
39
 *
 
40
 *    K_idx                 : the current element index in the mesh
 
41
 *    K.face(i)             : the i-th element face index in the mesh
 
42
 *    K.edge(i)             : the i-th element edge index in the mesh
 
43
 *    K[i]                  : the i-th element vertice index in the mesh
 
44
 *    mesh_n_geo (dim)      : total number of vertices(dim=0), edges(1), faces(2), volume(3)
 
45
 *                            in the mesh
 
46
 *    mesh_n_element (type) : total number of element of type 'p', 't', 'q', 'T', 'P', 'H'
 
47
 *                            in the mesh
 
48
 */
 
49
std::string
 
50
numbering_P2::name() const
 
51
{
 
52
  return "P2";
 
53
}
 
54
numbering_P2::size_type
 
55
numbering_P2::idof (
 
56
        const size_type*      mesh_n_geo,
 
57
        const size_type*      mesh_n_element,
 
58
        const geo_element&    K,
 
59
        size_type             i_dof_local) const
 
60
{
 
61
        // all vertices-based dof are numbered first
 
62
        // then all edge-based dof
 
63
        // and then all dof based on interior nodes (e.g. quadrangles)
 
64
 
 
65
        // for mixed triangle-quadrangle mesh : assume triangle are numbered first
 
66
        switch (K.type()) {
 
67
          case reference_element::p: 
 
68
            return K.index();
 
69
          case reference_element::e: 
 
70
            if (i_dof_local < 2) {
 
71
              return K [i_dof_local];
 
72
            } else {
 
73
              return mesh_n_geo[0] + K.index();
 
74
            }
 
75
          case reference_element::t: 
 
76
            if (i_dof_local < 3) {
 
77
              return K [i_dof_local];
 
78
            } else {
 
79
              return mesh_n_geo[0] + K.edge(i_dof_local - 3);
 
80
            }
 
81
          case reference_element::q: 
 
82
            if (i_dof_local < 4) {
 
83
              return K [i_dof_local];
 
84
            } else if (i_dof_local < 8) {
 
85
              return mesh_n_geo[0] + K.edge(i_dof_local - 4);
 
86
            } else {
 
87
              return mesh_n_geo[0] + mesh_n_geo[1]
 
88
                  + K.index() - mesh_n_element[reference_element::t];
 
89
            }
 
90
          case reference_element::T: 
 
91
            if (i_dof_local < 4) {
 
92
              return K [i_dof_local];
 
93
            } else {
 
94
              return mesh_n_geo[0] + K.edge(i_dof_local - 4);
 
95
            }
 
96
          case reference_element::P: 
 
97
          case reference_element::H:
 
98
          default:
 
99
            error_macro ("unsupported P2 element on `" << K.name() << "'");
 
100
            return 0;
 
101
        }
 
102
}
 
103
void
 
104
numbering_P2::idof (
 
105
        const size_type*      mesh_n_geo,
 
106
        const size_type*      mesh_n_element,
 
107
        const geo_element&    K, 
 
108
        vector<size_type>&    i_dof) const
 
109
{
 
110
  for (size_type i_dof_local = 0; i_dof_local < K.size(); i_dof_local++)
 
111
    i_dof[i_dof_local] 
 
112
     = numbering_P2::idof (mesh_n_geo, mesh_n_element, K, i_dof_local);
 
113
}
 
114
numbering_P2::size_type
 
115
numbering_P2::ndof (
 
116
              size_type  mesh_map_dimension,
 
117
        const size_type* mesh_n_geo,
 
118
        const size_type* mesh_n_element) const
 
119
{
 
120
    switch (mesh_map_dimension) {
 
121
      case 0:
 
122
        return mesh_n_geo[0];
 
123
      case 1:
 
124
        return mesh_n_geo[0] + mesh_n_geo[1];
 
125
      case 2:
 
126
        // number of vertices + nb of edges in the mesh
 
127
        // + number of quadrangle (extra interior node)
 
128
        return mesh_n_geo[0]
 
129
          + mesh_n_geo[1]
 
130
          + mesh_n_element[reference_element::q];
 
131
      case 3:
 
132
        return mesh_n_geo[0]
 
133
          + mesh_n_geo[1]
 
134
          + mesh_n_element[reference_element::P]
 
135
          + mesh_n_element[reference_element::H];
 
136
      default:
 
137
        fatal_macro ("unsupported P2 approximation in `" << mesh_map_dimension << "'");
 
138
        return 0;
 
139
    }
 
140
}