← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/lib/geo-tetgen.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:
nfem/lib/geo-tetgen.cc
 
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
//
 
22
// geo tetgen inputs / outputs
 
23
//
 
24
// author:
 
25
//      Pierre.Saramito@imag.fr
 
26
//
 
27
// date: 23 nov 2009
 
28
//
 
29
 
 
30
// ============================================================================
 
31
//  includes
 
32
// ============================================================================
 
33
 
 
34
#include "rheolef/georep.h"
 
35
#include "rheolef/iorheo.h"
 
36
#include "rheolef/rheostream.h" // i/o utility
 
37
#include "rheolef/field.h"
 
38
#include "rheolef/form.h"
 
39
#include "rheolef/tiny_element.h"
 
40
using namespace std;
 
41
 
 
42
 
 
43
// ============================================================================
 
44
// tetgen input (3d)
 
45
// use as input the concatenation of 3 files : .node .ele .face
 
46
// example:
 
47
// cat toto.node toto.elt toto.face | sed -e 's/#.*$/' | 
 
48
//              geo -input-tetgen -upgrade -geo - > toto.geo
 
49
// Note: strip comments lines as "#comments"
 
50
// ============================================================================
 
51
 
 
52
static
 
53
void
 
54
get_nodes_tetgen (
 
55
        istream& s, 
 
56
        georep::iterator_node p, 
 
57
        size_t np,
 
58
        size_t n_vert_attr,
 
59
        size_t n_bdr_mark,
 
60
        point& xmin, 
 
61
        point& xmax)
 
62
{
 
63
  typedef georep::size_type size_type;
 
64
  xmin[0] = xmin[1] = xmin[2] =  numeric_limits<Float>::max();
 
65
  xmax[0] = xmax[1] = xmax[2] = -numeric_limits<Float>::max();
 
66
  size_t vertice_number = 0;
 
67
  for (size_t i = 0; i < np; i++) {
 
68
    // <point #> <x> <y> <z> [attributes] [boundary marker]
 
69
    size_t idx;
 
70
    s >> idx; idx--;
 
71
    check_macro (idx < np, "node index "<<idx+1<<" is out of range 1:"<<np);
 
72
    s >> p[idx][0] >> p[idx][1] >> p[idx][2];
 
73
    // skip 0d domains
 
74
    size_t dummy;
 
75
    for (size_t j = 0; j < n_vert_attr + n_bdr_mark; j++) {
 
76
      s >> dummy;
 
77
    }
 
78
  }
 
79
}
 
80
//
 
81
// get boundary triangles
 
82
//
 
83
static
 
84
void
 
85
get_boundary_triangles_tetgen (istream& is,
 
86
    size_t n_bdr_triangles, bool has_face_attr, size_t n_vert,
 
87
    vector<domain>& domains,
 
88
    vector<size_t>& triangle_domain_tetgen_number)
 
89
{
 
90
  map<size_t, vector<tiny_element> > domain_map;
 
91
  tiny_element S;
 
92
  S.set_type (reference_element::t);
 
93
  for (size_t j = 0; is.good() && j < n_bdr_triangles; j++) {
 
94
    size_t S_idx;
 
95
    is >> S_idx; S_idx--;
 
96
    for (size_t i = 0 ; i < 3; i++) {
 
97
      is >> S[i]; S[i]--;
 
98
      check_macro (S[i] < n_vert, "tetgen input: face " << S_idx+1 << ": "
 
99
        << i+1 << "-th vertice index "
 
100
        << S[i]+1 << " out of range 1:" << n_vert);
 
101
    }
 
102
    size_t r = 0;
 
103
    if (has_face_attr) is >> r;
 
104
    domain_map[r].push_back (S);
 
105
  }
 
106
  check_macro (is, "a problem occurs while loading a tetgen mesh");
 
107
 
 
108
  // copy domain_map into domains (table)
 
109
  triangle_domain_tetgen_number.resize (domain_map.size()); 
 
110
  domains.resize (domain_map.size());
 
111
  vector<domain>::iterator iter_domain = domains.begin();
 
112
  size_t i_dom = 0;
 
113
  for (map<size_t, vector<tiny_element> >::const_iterator iter_domain_map = domain_map.begin();
 
114
      iter_domain_map != domain_map.end(); iter_domain_map++, iter_domain++, i_dom++) {
 
115
      triangle_domain_tetgen_number[i_dom] = (*iter_domain_map).first;
 
116
      const vector<tiny_element>& dom = (*iter_domain_map).second;
 
117
      if (dom.size() == 0) continue;
 
118
      (*iter_domain).set(dom.begin(), dom.size(), 2, "unamed"+itos(i_dom));
 
119
  }
 
120
}
 
121
istream&
 
122
georep::get_tetgen (istream& s, bool check_optional_domain_name) 
 
123
{
 
124
  // tetgen has only boundary edges that are numbered
 
125
  // and have domain numbers 
 
126
  // => this is only a version 1 mesh file format type
 
127
  //    not consistent for P2 elements
 
128
  _version = 1;
 
129
 
 
130
  typedef georep::size_type size_type;
 
131
  check_macro (s.good(), "bad input stream for tetgen.");
 
132
  bool verbose = iorheo::getverbose(s);
 
133
  //
 
134
  // get coordinates
 
135
  // header: <# of points> <dimension (must be 3)> <# of attributes> <# of boundary markers (0 or 1)>
 
136
  //
 
137
  size_type n_vert, n_vert_attr, n_bdr_mark;
 
138
  s >> n_vert >>_dim >> n_vert_attr >> n_bdr_mark;
 
139
  check_macro (s, "bad input stream for tetgen.");
 
140
  check_macro (_dim == 3, "bad dimension "<<_dim<<" for tetgen: expect 3");
 
141
  _x.resize(n_vert);
 
142
  get_nodes_tetgen (s, begin_node(), n_vert, n_vert_attr, n_bdr_mark, _xmin, _xmax);
 
143
  check_macro (s, "bad input stream for tetgen.");
 
144
  _count_geo [0] = _count_element [0] = n_vert;
 
145
  //
 
146
  // get elements
 
147
  // header: <# of tetrahedra> <nodes per tetrahedron> <# of attributes> 
 
148
  // lines : <tetrahedron #> <node> <node> <node> <node> ... [attributes]
 
149
  size_type n_elt, n_vert_per_elt, n_tetra_attr;
 
150
  s >> n_elt >> n_vert_per_elt >> n_tetra_attr;
 
151
  check_macro (n_vert_per_elt == 4, "bad #node per tetra "<<n_vert_per_elt<<" for tetgen: expect 4");
 
152
  resize(n_elt);
 
153
  iterator elt = begin();
 
154
  for (size_t j = 0; j < n_elt; j++) {
 
155
    size_t K_idx;
 
156
    s >> K_idx; K_idx--;
 
157
    check_macro (K_idx < n_elt, "element index "<<K_idx+1<<" is out of range 1:"<<n_elt);
 
158
    geo_element& K = elt[K_idx];
 
159
    K.set_type (n_vert_per_elt, _dim);
 
160
    K.set_index(K_idx) ;
 
161
    for (size_type i = 0 ; i < n_vert_per_elt; i++) {
 
162
      s >> K[i]; K[i]--;
 
163
      check_macro (K[i] < n_vert, "tetgen input: element " << K_idx+1 << ": "
 
164
        << i+1 << "-th vertice index "
 
165
        << K[i]+1 << " out of range 1:" << n_vert);
 
166
    }
 
167
    // skip 3d domains
 
168
    size_t dummy;
 
169
    for (size_t j = 0; j < n_tetra_attr; j++) {
 
170
      s >> dummy;
 
171
    }
 
172
  }
 
173
  // read triangular faces:
 
174
  // First line: <# of faces> <boundary marker (0 or 1)>
 
175
  // Remaining lines list of # of faces:
 
176
  // <face #> <node> <node> <node> [boundary marker]
 
177
  size_t n_face, bdry_mark;
 
178
  s >> n_face >> bdry_mark;
 
179
  bool has_face_attr = (bdry_mark == 1);
 
180
  vector<domain> triangle_domain;
 
181
  vector<size_t> triangle_domain_tetgen_number;
 
182
  get_boundary_triangles_tetgen (s, n_face, has_face_attr, n_vert,
 
183
        triangle_domain, triangle_domain_tetgen_number);
 
184
 
 
185
  // ---------------------------------------------------------------
 
186
  // check extension to optional domain names
 
187
  // ---------------------------------------------------------------
 
188
  if (!check_optional_domain_name) return s;
 
189
  char c;
 
190
  s >> ws >> c; // skip white and grab a char
 
191
  // have "FaceDomainNames", "EdgeDomainNames" or "VerticeDomainNames" ?
 
192
  // tetgen mesh may be followed by field data and such, so be carrefull...
 
193
  while (c == 'E' || c == 'F' || c == 'V') {
 
194
      s.unget(); // put char back
 
195
      if (c == 'F') {
 
196
        if (!scatch(s,"FaceDomainNames")) return s;
 
197
        // ---------------------------------------------------------------
 
198
        // get triangle domain names
 
199
        // ---------------------------------------------------------------
 
200
        size_type n_dom_triangle;
 
201
        s >> n_dom_triangle;
 
202
        if (n_dom_triangle != triangle_domain.size()) {
 
203
           string numbers;
 
204
           for (size_t j = 0; j < triangle_domain_tetgen_number.size(); j++) {
 
205
              numbers = numbers + itos(triangle_domain_tetgen_number[j]) + " ";
 
206
           }
 
207
           warning_macro("geo input: "
 
208
            << n_dom_triangle << " domain names founded while "
 
209
            << triangle_domain.size() << " tetgen triangle domains numbers are defined: "
 
210
            << numbers);
 
211
           cerr << endl;
 
212
           error_macro("HINT: check domain name file (.dmn)");
 
213
        }
 
214
        size_t start_triangle_domain = _domlist.size();
 
215
        _domlist.resize(_domlist.size() + n_dom_triangle);
 
216
        for (size_type k = 0; k < n_dom_triangle; k++) {
 
217
          string name;
 
218
          s >> name;
 
219
          triangle_domain[k].set_name(name);
 
220
          _domlist.at(start_triangle_domain+k) = triangle_domain[k];
 
221
        }
 
222
      } else if (c == 'E') {
 
223
        if (!scatch(s,"EdgeDomainNames")) return s;
 
224
        // ---------------------------------------------------------------
 
225
        // get edge domain names
 
226
        // ---------------------------------------------------------------
 
227
        warning_macro ("EdgeDomainNames skiped");
 
228
        size_type n_dom_edge;
 
229
        s >> n_dom_edge;
 
230
        for (size_type k = 0; k < n_dom_edge; k++) {
 
231
          string name;
 
232
          s >> name;
 
233
        }
 
234
      } else if (c == 'V') {
 
235
        if (!scatch(s,"VerticeDomainNames")) return s;
 
236
        // ---------------------------------------------------------------
 
237
        // get vertice domain names
 
238
        // ---------------------------------------------------------------
 
239
        warning_macro ("VerticeDomainNames skiped");
 
240
        size_type n_dom_vertice;
 
241
        s >> n_dom_vertice;
 
242
        for (size_type k = 0; k < n_dom_vertice; k++) {
 
243
          string name;
 
244
          s >> name;
 
245
        }
 
246
      }
 
247
      s >> ws >> c; // skip white and grab a char
 
248
  }
 
249
  return s;
 
250
}
 
251
// =======================================================================
 
252
// tetgen output
 
253
// =======================================================================
 
254
ostream&
 
255
georep::put_tetgen (ostream& os, bool only_boundary) const
 
256
{
 
257
  fatal_macro ("put_tetgen: not yet!");
 
258
  return os;
 
259
}