Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
msh2geo.cc
Go to the documentation of this file.
1
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// the msh2geo unix command
22// author: Pierre.Saramito@imag.fr
23
24namespace rheolef {
122} // namespace rheolef
123
124#include "rheolef/compiler.h" // after index_set to avoid boost
125#include "scatch.icc"
126#include "index_set_header.icc"
127#include "index_set_body.icc"
128#include <array>
129#include <cstring>
130
131namespace rheolef {
132
133template<class T>
134struct point_basic: std::array<T,3> {
135 typedef std::array<T,3> base;
136 point_basic(T x0=T(), T x1=T(), T x2=T())
137 : base() {
138 base::operator[](0) = x0;
139 base::operator[](1) = x1;
140 base::operator[](2) = x2;
141 }
142};
144
145namespace edge {
146#include "edge.icc"
147} // namespace edge
148
149namespace triangle {
150#include "triangle.icc"
151} // namespace triangle
152
153namespace quadrangle {
154#include "quadrangle.icc"
155} // namespace quadrangle
156
157namespace tetrahedron {
158#include "tetrahedron.icc"
159} // namespace tetrahedron
160
161namespace prism {
162#include "prism.icc"
163} // namespace prism
164
165namespace hexahedron {
166#include "hexahedron.icc"
167} // namespace hexahedron
168
169} // namespace rheolef
170
171#include "reference_element_aux.icc"
172
173namespace rheolef {
174
175// TODO: move somewhere in reference_element_xxx.cc
176const char
178 'p',
179 'e',
180 't',
181 'q',
182 'T',
183 'P',
184 'H'
185};
186} // namespace rheolef
187
188namespace rheolef {
189
191typedef int shift_type;
192
193// side (edge & face): element with reduced node list when order > 1, for sides (edges or faces)
194struct geo_element_side: std::array<size_t,4> {
195 void setname (char name) { _variant = reference_element_variant(name); }
196 void setindex (size_t index) { _index = index; }
197 char name() const { return _reference_element_name [_variant]; }
198 size_t index() const { return _index; }
199 size_t variant() const { return _variant; }
200 size_t dimension()const { return reference_element_dimension_by_variant(variant()); }
201 size_t n_vertex() const { return reference_element_n_vertex (variant()); }
202 size_t size() const { return n_vertex(); }
203 geo_element_side()
204 : array<size_t,4>(), _variant(-1), _index(-1)
205 {}
206protected:
207 size_t _variant;
208 size_t _index;
209};
210// side index with orient and shift for geo_element
212 typedef int orientation_type;
213 typedef int shift_type;
214 void setindex (size_t index) { _index = index; }
215 void setorientation (size_t orient) { _orient = orient; }
216 void setshift (size_t shift) { _shift = shift; }
217 size_t index() const { return _index; }
218 int orientation() const { return _orient; }
219 size_t shift() const { return _shift; }
222protected:
223 size_t _index;
226};
227// element with full node list when order > 1
228struct geo_element: vector<size_t> {
229 void setname (char name) { _name = name; }
230 void setorder (size_t order) { _order = order; }
231 void setindex (size_t index) { _index = index; }
233 geo_element_indirect& edge_indirect (size_t i) { return _edge[i]; }
234 geo_element_indirect& face_indirect (size_t i) { return _face[i]; }
235 char name() const { return _name; }
236 size_t order() const { return _order; }
237 size_t index() const { return _index; }
238 size_t gmshtype() const { return _gmshtype; }
239 size_t variant() const { return reference_element_variant (name()); }
240 size_t dimension()const { return reference_element_dimension_by_name(name()); }
241 size_t n_vertex() const { return reference_element_n_vertex (variant()); }
242 const geo_element_indirect& edge_indirect (size_t i) const { return _edge[i]; }
243 const geo_element_indirect& face_indirect (size_t i) const { return _face[i]; }
244 size_t n_edge() const { return _reference_element_n_edge_by_variant [variant()]; }
245 size_t edge (size_t i) const { return _edge[i].index(); }
246 size_t edge_local_vertex(size_t iedg, size_t edg_iv) const;
247 size_t face (size_t i) const { return _face[i].index(); }
248 size_t n_face() const { return _reference_element_n_face_by_variant [variant()]; }
249 size_t face_n_vertex(size_t loc_ifac) const;
250 size_t face_local_vertex(size_t iedg, size_t edg_iv) const;
252 : vector<size_t>(), _name('z'), _order(-1), _index(-1), _gmshtype(-1), _edge(), _face()
253 {}
254protected:
255 char _name; // TODO: store variant instead of name
256 size_t _order;
257 size_t _index;
258 size_t _gmshtype;
259 array<geo_element_indirect,12> _edge;
260 array<geo_element_indirect,6> _face;
261};
262size_t
263geo_element::edge_local_vertex (size_t iedg, size_t edg_iv) const
264{
265 switch (variant()) {
266 case reference_element__e: return vector<size_t>::operator[] (edg_iv);
267 case reference_element__t: return triangle::edge [iedg][edg_iv%2];
268 case reference_element__q: return quadrangle::edge [iedg][edg_iv%2];
269 case reference_element__T: return tetrahedron::edge [iedg][edg_iv%2];
270 case reference_element__P: return prism::edge [iedg][edg_iv%2];
271 case reference_element__H: return hexahedron::edge [iedg][edg_iv%2];
272 default: error_macro ("invalid variant"); return 0;
273 }
274}
275size_t
276geo_element::face_local_vertex (size_t loc_ifac, size_t fac_iv) const
277{
278 switch (variant()) {
279 case reference_element__t:
280 case reference_element__q: return vector<size_t>::operator[] (fac_iv);
281 case reference_element__T: return tetrahedron::face [loc_ifac][fac_iv%3];
282 case reference_element__P: return prism::face [loc_ifac][fac_iv%4];
283 case reference_element__H: return hexahedron::face [loc_ifac][fac_iv%4];
284 default: error_macro ("invalid variant"); return 0;
285 }
286}
287size_t
288geo_element::face_n_vertex(size_t loc_ifac) const
289{
290 switch (variant()) {
291 case reference_element__t: return 3;
292 case reference_element__q: return 4;
293 case reference_element__T: return 3;
294 case reference_element__H: return 4;
295 case reference_element__P: return (loc_ifac < 2) ? 3 : 4;
296 default: error_macro ("invalid variant"); return 0;
297 }
298}
299
300} // namespace rheolef
301
302#include "msh2geo_defs.icc"
303#include "msh2geo_node_renum.icc"
304#include "geo_element_aux.icc"
305
306using namespace std;
307using namespace rheolef;
308// necessite de numeroter les aretes
309// puis de le stoker dans ts les triangles
310// le num arete et son orient
311// => doit etre appeler apres la num des aretes
312void
314 size_t order,
315 const array<size_t,reference_element__max_variant>& size_by_variant,
316 const array<size_t,reference_element__max_variant+1>& first_by_variant,
317 const array<size_t,reference_element__max_variant>& loc_ndof_by_variant,
318 const geo_element& K,
319 vector<size_t>& dis_idof1)
320{
321 dis_idof1.resize (reference_element_n_node (K.variant(), order));
322 std::fill (dis_idof1.begin(), dis_idof1.end(), std::numeric_limits<size_t>::max());
323
324 for (size_t subgeo_variant = 0, variant = K.variant(); subgeo_variant <= variant; subgeo_variant++) {
325 size_t subgeo_dim = reference_element_dimension_by_variant (subgeo_variant);
326 size_t loc_sub_ndof = loc_ndof_by_variant [subgeo_variant];
327 for (size_t first_loc_idof = reference_element_first_inod_by_variant (variant, order, subgeo_variant),
328 last_loc_idof = reference_element_last_inod_by_variant (variant, order, subgeo_variant),
329 loc_idof = first_loc_idof;
330 loc_idof < last_loc_idof; loc_idof++) {
331 // 1) local loc_igev on subgeo
332 size_t loc_igev = (loc_idof - first_loc_idof) / loc_sub_ndof;
333 size_t loc_igev_j = (loc_idof - first_loc_idof) % loc_sub_ndof;
334 // 2) then compute ige; computation depends upon the subgeo dimension:
335 size_t ige;
336 switch (subgeo_dim) {
337 case 0: {
338 // convert node numbering to vertex numbering, for geo order > 1
339 size_t loc_inod = loc_idof; // only one dof per vertex
340 ige = K [loc_inod];
341 break;
342 }
343 case 1: {
344 loc_igev_j = geo_element_fix_edge_indirect (K, loc_igev, loc_igev_j, order);
345 size_t loc_ige = loc_igev;
346 ige = K.edge_indirect(loc_ige).index();
347 break;
348 }
349 case 2: {
350 size_t loc_ige = loc_igev;
351 if (subgeo_variant == reference_element__t) {
352 loc_igev_j = geo_element_fix_triangle_indirect (K, loc_igev, loc_igev_j, order);
353 } else {
354 size_t loc_ntri = (K.variant() == reference_element__P) ? 2 : 0;
355 loc_ige += loc_ntri;
356 loc_igev_j = geo_element_fix_quadrangle_indirect (K, loc_ige, loc_igev_j, order);
357 }
358 ige = K.face(loc_ige);
359 break;
360 }
361 case 3: {
362 ige = K.index();
363 break;
364 }
365 default: {
366 error_macro ("unexpected subgeo_dim="<<subgeo_dim);
367 }
368 }
369 // 3) then compute igev, by variant
370 size_t igev = ige;
371 for (size_t prev_subgeo_variant = reference_element_first_variant_by_dimension(subgeo_dim);
372 prev_subgeo_variant < subgeo_variant;
373 prev_subgeo_variant++) {
374 size_t nprev = size_by_variant [prev_subgeo_variant];
375 assert_macro (ige >= nprev, "invalid index");
376 igev -= nprev;
377 }
378 size_t dis_igev = igev;
379 // 4) finally compute dis_idof
380 size_t dis_idof = 0;
381 for (size_t prev_subgeo_variant = 0;
382 prev_subgeo_variant < subgeo_variant;
383 prev_subgeo_variant++) {
384 dis_idof += size_by_variant [prev_subgeo_variant]
385 *loc_ndof_by_variant [prev_subgeo_variant];
386 }
387 dis_idof += dis_igev
388 *loc_ndof_by_variant [subgeo_variant]
389 + loc_igev_j;
390 assert_macro (loc_idof < dis_idof1.size(), "msh2geo: invalid index");
391 dis_idof1 [loc_idof] = dis_idof;
392 }
393 }
394}
395void
397 ostream& out,
398 size_t dim,
399 size_t dom_dim,
400 const map<size_t,list<size_t> >& domain_map,
401 map<size_t, string>& phys,
402 const vector<geo_element>& element)
403{
404 if (dim == dom_dim && domain_map.size() == 1) return;
405 for (map<size_t,list<size_t> >::const_iterator
406 first = domain_map.begin(),
407 last = domain_map.end(); first != last; first++) {
408 size_t dom_idx = (*first).first;
409 const list<size_t>& dom = (*first).second;
410 string dom_name = phys[dom_idx];
411 if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
412 out << "domain" << endl
413 << dom_name << endl
414 << "1 " << dom_dim << " " << dom.size() << endl;
415 for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
416 variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
417 for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
418 size_t ie = *fe;
419 if (element[ie].variant() != variant) continue;
420 out << element[ie].name() << "\t";
421 if (element[ie].order() > 1) {
422 out << "p" << element[ie].order() << " ";
423 }
424 for (size_t j = 0; j < element[ie].size(); j++) {
425 out << element[ie][j] << " ";
426 }
427 out << endl;
428 }
429 }
430 out << endl;
431 }
432}
433void
435 ostream& out,
436 size_t dim,
437 size_t dom_dim,
438 const map<size_t,list<size_t> >& domain_map,
439 map<size_t, string>& phys,
440 const vector<geo_element>& element,
441 const vector<geo_element_side>& edge,
442 const vector<index_set>& edge_ball,
443 const vector<geo_element_side>& face,
444 const vector<index_set>& face_ball)
445{
446 if (dim == dom_dim && domain_map.size() == 1) return;
447 for (map<size_t,list<size_t> >::const_iterator
448 first = domain_map.begin(),
449 last = domain_map.end(); first != last; first++) {
450 size_t dom_idx = (*first).first;
451 const list<size_t>& dom = (*first).second;
452 string dom_name = phys[dom_idx];
453 if (dom_name == "") dom_name = "unnamed" + std::to_string(dom_idx);
454 out << "domain" << endl
455 << dom_name << endl
456 << "2 " << dom_dim << " " << dom.size() << endl;
457 for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
458 variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
459 for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
460 size_t ie = *fe;
461 if (element[ie].variant() != variant) continue;
462 size_t isid = 0;
463 orientation_type orient = 1;
464 if (element[ie].name() == 'p') {
465 isid = element[ie][0];
466 orient = 1;
467 } else if (element[ie].name() == 'e') {
468 size_t inod0 = element[ie][0];
469 size_t inod1 = element[ie][1];
470 index_set iedge_set = edge_ball[inod0];
471 iedge_set.inplace_intersection (edge_ball[inod1]);
472 check_macro (iedge_set.size() == 1,
473 "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
474 isid = *(iedge_set.begin());
475 orient = (inod0 == edge[isid][0]) ? 1 : -1;
476 } else if (element[ie].dimension() == 2) {
477 size_t inod0 = element[ie][0];
478 index_set iface_set = face_ball[inod0];
479 for (size_t j = 1; j < element[ie].n_vertex(); ++j) {
480 size_t inodj = element[ie][j];
481 iface_set.inplace_intersection (face_ball[inodj]);
482 }
483 check_macro (iface_set.size() == 1,
484 "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
485 isid = *(iface_set.begin());
486 shift_type shift;
487 if (element[ie].name() == 't') {
488 geo_element_get_orientation_and_shift (
489 face[isid],
490 element[ie][0],
491 element[ie][1],
492 element[ie][2],
493 orient, shift);
494 } else {
495 geo_element_get_orientation_and_shift (
496 face[isid],
497 element[ie][0],
498 element[ie][1],
499 element[ie][2],
500 element[ie][3],
501 orient, shift);
502 }
503 } else { // 3d domain
504 isid = element[ie].index();
505 orient = 1;
506 }
507 out << orient*int(isid) << endl;
508 }
509 }
510 out << endl;
511 }
512}
513void msh2geo (istream& in, ostream& out, string sys_coord_name, bool do_upgrade)
514{
515 // ----------------------------------
516 // 1. input gmsh
517 // ----------------------------------
518 // 1.0. preambule
519 check_macro (scatch(in,"$MeshFormat",true),
520 "msh2geo: input stream does not contains a gmsh mesh file ($MeshFormat not found).");
521 double gmsh_fmt_version;
522 size_t file_type, float_data_size;
523 in >> gmsh_fmt_version >> file_type >> float_data_size;
524 check_macro (gmsh_fmt_version <= 2.2,
525 "gmsh format version " << gmsh_fmt_version << " founded ; expect version 2.2 (HINT: use gmsh -format msh2)");
526 check_macro (file_type == 0, "msh2geo: unsupported gmsh non-ascii format");
527 check_macro (scatch(in,"$EndMeshFormat",true), "msh2geo: gmsh input error: $EndMeshFormat not found.");
528 //
529 // 1.1 optional domain names
530 //
531 bool full_match = true;
532 bool partial_match = !full_match;
533 check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
534 size_t nphys;
535 map<size_t, string> phys;
536 string label;
537 in >> label;
538 size_t n_names = 0;
539 if (label == "PhysicalNames") {
540 in >> nphys;
541 for (size_t i = 0; i < nphys; i++) {
542 string name;
543 size_t dom_dim, dom_idx;
544 in >> dom_dim >> dom_idx;
545 // get name:
546 char c;
547 in >> std::ws >> c;
548 if (c != '"') name.push_back(c);
549 do {
550 in.get(c);
551 name.push_back(c);
552 } while (c != '"' && c != '\n');
553 // strip '"' in name
554 size_t start = 0, end = name.length();
555 if (name[start] == '"') start++;
556 if (name[end-1] == '"') end--;
557 name = name.substr(start,end-start);
558 // rename spaces and tabs
559 for (size_t i = 0; i < name.size(); i++) {
560 if (name[i] == ' ' || name[i] == '\t') name[i] = '_';
561 }
562 phys[dom_idx] = name.substr(start,end-start);
563 }
564 check_macro (scatch(in,"$EndPhysicalNames",true), "msh2geo: gmsh input error ($EndPhysicalNames not found).");
565 check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
566 in >> label;
567 }
568 //
569 // 1.2. nodes
570 //
571 if (label != "Nodes") {
572 check_macro (scatch(in,"$Nodes",true), "msh2geo: $Nodes not found in gmsh file");
573 }
574 size_t nnode;
575 in >> nnode;
576 vector<point> node (nnode);
577 double infty = numeric_limits<double>::max();
578 point xmin ( infty, infty, infty);
579 point xmax (-infty, -infty, -infty);
580 for (size_t k = 0; k < node.size(); k++) {
581 size_t dummy;
582 in >> dummy;
583 for (size_t i = 0; i < 3; ++i) {
584 in >> node[k][i];
585 }
586 for (size_t j = 0 ; j < 3; j++) {
587 xmin[j] = min(node[k][j], xmin[j]);
588 xmax[j] = max(node[k][j], xmax[j]);
589 }
590 }
591 //
592 // dimension is deduced from bounding box
593 //
594 size_t dim = 3;
595 if (xmax[2] == xmin[2]) {
596 dim = 2;
597 if (xmax[1] == xmin[1]) {
598 dim = 1;
599 if (xmax[0] == xmin[0]) dim = 0;
600 }
601 }
602 //
603 // 1.3. elements
604 //
605 check_macro (scatch(in,"$Elements",true), "msh2geo: $Elements not found in gmsh file");
606 size_t n_element;
607 in >> n_element;
608 vector<geo_element> element (n_element);
609 vector<size_t> node_subgeo_variant (node.size(), std::numeric_limits<size_t>::max());
610 array<map<size_t, list<size_t> >,4> domain_map; // domain_map[dim][idom][ie]
611 size_t map_dim = 0;
612 size_t order = 0;
613 for (size_t i = 0; i < element.size(); i++) {
614 size_t id, dummy, gmshtype;
615 in >> id >> gmshtype;
616 size_t n_tag_gmsh;
617 in >> n_tag_gmsh;
618 size_t domain_idx = 0;
619 for (size_t j = 0 ; j < n_tag_gmsh; j++) {
620 // the first tag is the physical domain index
621 // the second tag is the object index, defined for all elements
622 // the third is zero (in all examples)
623 size_t tag_dummy;
624 in >> tag_dummy;
625 if (j == 0) {
626 domain_idx = tag_dummy;
627 }
628 }
629 check_macro (gmshtype < gmshtype_max,
630 "msh2geo: element #" << id << ": unexpected gmsh type '" << gmshtype << "'");
631 check_macro (gmsh_table[gmshtype].supported,
632 "msh2geo: element #" << id << ": unsupported gmsh type '" << gmshtype << "'");
633
634 element[i].setname (gmsh_table[gmshtype].name);
635 element[i].setorder(gmsh_table[gmshtype].order);
636 element[i].setgmshtype(gmshtype);
637 size_t nv = gmsh_table[gmshtype].nv;
638 size_t nn = gmsh_table[gmshtype].nn_tot;
639 size_t dim_i = element[i].dimension();
640 size_t variant = element[i].variant();
641 map_dim = max(map_dim,dim_i);
642 if (order == 0) {
643 order = element[i].order();
644 } else {
645 check_macro (order == element[i].order(), "msh2geo: unexpected order="<<element[i].order());
646 }
647 element[i].resize(nn);
648 for (size_t j = 0; j < nn; j++) {
649 in >> element[i][j];
650 element[i][j]--;
651 }
652 for (size_t subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) { // set node dimension
653 for (size_t loc_inod = reference_element_first_inod_by_variant(variant,element[i].order(),subgeo_variant),
654 loc_nnod = reference_element_last_inod_by_variant(variant,element[i].order(),subgeo_variant);
655 loc_inod < loc_nnod; loc_inod++) {
656 node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
657 }
658 }
659 // from gmsh to rheolef/geo local node renumbering: element[i][j] are modified
660 msh2geo_node_renum (element[i], element[i].name(), element[i].order());
661 if (domain_idx != 0) {
662 domain_map[dim_i][domain_idx].push_back(i);
663 }
664 }
665 array<size_t,reference_element__max_variant> loc_ndof_by_variant;
666 reference_element_init_local_nnode_by_variant (order, loc_ndof_by_variant);
667 // set dis_ie: by increasing variant order, for map_dim element only
668 size_t last_index = 0;
669 for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
670 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
671 for (size_t ie = 0; ie < element.size(); ++ie) {
672 if (element[ie].variant() != variant) continue;
673 element[ie].setindex (last_index);
674 switch (element[ie].dimension()) {
675 case 0: element[ie].resize(1); element[ie][0] = last_index; break;
676 case 1: element[ie].edge_indirect(0).setindex(last_index); break;
677 case 2: element[ie].face_indirect(0).setindex(last_index); break;
678 default: break;
679 }
680 last_index++;
681 }
682 }
683 // -------------------------------------------------
684 // 2. node reordering, first pass
685 // -------------------------------------------------
686 // permut: first vertex, then edge, faces and inernal nodes, by increasing subgeo_dim
687 vector<size_t> old2new_inode (node.size(), std::numeric_limits<size_t>::max());
688 if (true || !do_upgrade) {
689 // 2.1. when no upgrade
690 size_t new_inode = 0;
691 for (size_t subgeo_variant = 0; subgeo_variant < reference_element_last_variant_by_dimension(map_dim); subgeo_variant++) {
692 for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
693 if (node_subgeo_variant [old_inode] != subgeo_variant) continue;
694 old2new_inode[old_inode] = new_inode++;
695 }
696 }
697 } else {
698 // 2.2. when no upgrade: will be renumbered after side creation
699 for (size_t inode = 0; inode < node.size(); inode++) {
700 old2new_inode[inode] = inode;
701 }
702 }
703 for (size_t inode = 0; inode < node.size(); inode++) {
704 check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
705 }
706 // inv permut
707 vector<size_t> new2old_inode (node.size(), std::numeric_limits<size_t>::max());
708 for (size_t inode = 0; inode < node.size(); inode++) {
709 new2old_inode[old2new_inode[inode]] = inode;
710 }
711 for (size_t inode = 0; inode < node.size(); inode++) {
712 check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
713 }
714 // for convenience, apply the new numbering to nodes and elements
715 {
716 vector<point> new_node (node.size());
717 for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
718 new_node [old2new_inode[old_inode]] = node [old_inode];
719 }
720 for (size_t inode = 0; inode < node.size(); inode++) {
721 node [inode] = new_node [inode];
722 }
723 }
724 for (size_t ie = 0; ie < element.size(); ++ie) {
725 for (size_t i = 0; i < element[ie].size(); ++i) {
726 element[ie][i] = old2new_inode[element[ie][i]];
727 }
728 }
729 // ------------------------------------------------------------------------
730 // 3.1) compute all edges
731 // ------------------------------------------------------------------------
732 vector<index_set> edge_ball (node.size());
733 vector<geo_element_side> edge;
734 if (do_upgrade && map_dim >= 2) {
735 list <geo_element_side> edge_list;
736 size_t nedg = 0;
737 for (size_t ie = 0; ie < element.size(); ++ie) {
738 if (element[ie].dimension() != map_dim) continue;
739 for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
740 size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
741 size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
742 size_t inod0 = element[ie][iv0];
743 size_t inod1 = element[ie][iv1];
744 index_set iedge_set = edge_ball[inod0];
745 iedge_set.inplace_intersection (edge_ball[inod1]);
746 check_macro (iedge_set.size() <= 1,
747 "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
748 if (iedge_set.size() == 1) continue; // edge already exists
749 edge_ball[inod0].insert (nedg);
750 edge_ball[inod1].insert (nedg);
751 geo_element_side new_edge;
752 new_edge.setname('e');
753 new_edge[0] = inod0;
754 new_edge[1] = inod1;
755 new_edge.setindex(nedg);
756 edge_list.push_back (new_edge);
757 nedg++;
758 }
759 }
760 edge.resize (edge_list.size());
761 std::copy (edge_list.begin(), edge_list.end(), edge.begin());
762 }
763 // propagate edge numbering inside elements
764 if (do_upgrade && map_dim >= 2) {
765 for (size_t ie = 0; ie < element.size(); ++ie) {
766 if (element[ie].dimension() != map_dim) continue;
767 for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
768 size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
769 size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
770 size_t inod0 = element[ie][iv0];
771 size_t inod1 = element[ie][iv1];
772 index_set iedge_set = edge_ball[inod0];
773 iedge_set.inplace_intersection (edge_ball[inod1]);
774 check_macro (iedge_set.size() == 1,
775 "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
776 size_t iedge = *(iedge_set.begin());
777 element[ie].edge_indirect(iedg).setindex (iedge);
778 element[ie].edge_indirect(iedg).setorientation((edge[iedge][0] == inod0) ? 1 : -1);
779 }
780 }
781 }
782 // ------------------------------------------------------------------------
783 // 3.2) compute all faces
784 // ------------------------------------------------------------------------
785 vector<index_set> face_ball (node.size());
786 vector<geo_element_side> face;
787 if (do_upgrade && map_dim >= 3) {
788 list <geo_element_side> face_list;
789 size_t nfac = 0;
790 for (size_t ie = 0; ie < element.size(); ++ie) {
791 if (element[ie].dimension() != map_dim) continue;
792 for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
793 size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
794 size_t inod0 = element[ie][iv0];
795 index_set iface_set = face_ball[inod0];
796 for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
797 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
798 size_t inodj = element[ie][ivj];
799 iface_set.inplace_intersection (face_ball[inodj]);
800 }
801 check_macro (iface_set.size() <= 1,
802 "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
803 if (iface_set.size() == 1) continue; // face already exists
804 geo_element_side new_face;
805 char face_name = (element[ie].face_n_vertex(loc_ifac) == 3) ? 't' : 'q';
806 new_face.setname(face_name);
807 new_face.setindex(nfac);
808 for (size_t j = 0; j < element[ie].face_n_vertex(loc_ifac); ++j) {
809 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
810 size_t inodj = element[ie][ivj];
811 face_ball[inodj].insert (nfac);
812 new_face[j] = inodj;
813 }
814 face_list.push_back (new_face);
815 nfac++;
816 }
817 }
818 face.resize (face_list.size());
819 std::copy (face_list.begin(), face_list.end(), face.begin());
820 }
821 // propagate face numbering inside elements
822 if (do_upgrade && map_dim >= 3) {
823 for (size_t ie = 0; ie < element.size(); ++ie) {
824 if (element[ie].dimension() != map_dim) continue;
825 for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
826 size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
827 size_t inod0 = element[ie][iv0];
828 index_set iface_set = face_ball[inod0];
829 for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
830 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
831 size_t inodj = element[ie][ivj];
832 iface_set.inplace_intersection (face_ball[inodj]);
833 }
834 check_macro (iface_set.size() == 1,
835 "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
836 size_t iface = *(iface_set.begin());
837 orientation_type orient;
838 shift_type shift;
839 if (element[ie].face_n_vertex(loc_ifac) == 3) {
840 geo_element_get_orientation_and_shift (
841 face[iface],
842 element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
843 element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
844 element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
845 orient, shift);
846 } else {
847 geo_element_get_orientation_and_shift (
848 face[iface],
849 element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
850 element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
851 element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
852 element[ie][element[ie].face_local_vertex(loc_ifac, 3)],
853 orient, shift);
854 }
855 element[ie].face_indirect(loc_ifac).setindex (iface);
856 element[ie].face_indirect(loc_ifac).setorientation (orient);
857 element[ie].face_indirect(loc_ifac).setshift (shift);
858 }
859 }
860 }
861 // ------------------------------------------------------------------------
862 // 3.4) count all elements
863 // ------------------------------------------------------------------------
864 array<size_t,reference_element__max_variant> size_by_variant;
865 array<size_t,reference_element__max_variant+1> first_by_variant;
866 size_by_variant.fill (0);
867 first_by_variant.fill (0);
868 if (true) {
869 size_t n_vertex = 0;
870 for (size_t inode = 0; inode < node.size(); inode++) {
871 if (node_subgeo_variant [inode] == reference_element__p) n_vertex++;
872 }
873 size_by_variant [reference_element__p] = n_vertex;
874 if (map_dim >= 2 && edge.size() != 0) {
875 size_by_variant [reference_element__e] = edge.size();
876 }
877 if (map_dim >= 3 && face.size() != 0) {
878 for (size_t loc_ifac = 0; loc_ifac < face.size(); ++loc_ifac) {
879 size_by_variant [face[loc_ifac].variant()]++;
880 }
881 }
882 for (size_t ie = 0; ie < element.size(); ++ie) {
883 if (element[ie].dimension() != map_dim) continue;
884 size_t variant = element[ie].variant();
885 size_by_variant [variant]++;
886 }
887 for (size_t dim = 0; dim <= 3; ++dim) {
888 for (size_t variant = reference_element_first_variant_by_dimension(dim);
889 variant < reference_element_last_variant_by_dimension(dim); variant++) {
890 first_by_variant [variant+1] = size_by_variant [variant];
891 }
892 }
893 }
894 // -------------------------------------------------
895 // 4. node reordering, second pass
896 // -------------------------------------------------
897 // need edges & faces with index & orient
898 if (do_upgrade) {
899 std::fill (old2new_inode.begin(), old2new_inode.end(), std::numeric_limits<size_t>::max());
900 std::fill (new2old_inode.begin(), new2old_inode.end(), std::numeric_limits<size_t>::max());
901 for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
902 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
903 for (size_t ie = 0; ie < element.size(); ++ie) {
904 if (element[ie].variant() != variant) continue;
905 const geo_element& old_K = element[ie];
906 std::vector<size_t> new_K;
908 order, size_by_variant, first_by_variant, loc_ndof_by_variant, old_K, new_K);
909 for (size_t loc_inod = 0; loc_inod < old_K.size(); loc_inod++) {
910 size_t old_inod = old_K [loc_inod];
911 size_t new_inod = new_K [loc_inod];
912 if (old2new_inode [old_inod] != std::numeric_limits<size_t>::max()) continue; // already done
913 old2new_inode [old_inod] = new_inod;
914 }
915 }
916 }
917 } else {
918 // when no upgrade:
919 for (size_t inode = 0; inode < node.size(); inode++) {
920 old2new_inode[inode] = inode;
921 }
922 }
923 for (size_t inode = 0; inode < node.size(); inode++) {
924 check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
925 }
926 // inv permut
927 for (size_t inode = 0; inode < node.size(); inode++) {
928 new2old_inode[old2new_inode[inode]] = inode;
929 }
930 for (size_t inode = 0; inode < node.size(); inode++) {
931 check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
932 }
933#ifdef TODO
934#endif // TODO
935 // ----------------------------------
936 // 5. output geo
937 // ----------------------------------
938 // 5.1. output the mesh
939 size_t version = 4;
940 static const char* geo_variant_name [reference_element__max_variant] = {
941 "points",
942 "edges",
943 "triangles",
944 "quadrangles",
945 "tetrahedra",
946 "prisms",
947 "hexahedra"
948 };
949 out << setprecision(numeric_limits<double>::digits10)
950 << "#!geo" << endl
951 << endl
952 << "mesh" << endl
953 << version << endl
954 << "header" << endl
955 << " dimension\t" << dim << endl;
956 if (sys_coord_name != "cartesian") {
957 out << " coordinate_system " << sys_coord_name << endl;
958 }
959 if (order != 1) {
960 out << " order\t\t" << order << endl;
961 }
962 out << " nodes\t\t" << node.size() << endl;
963
964 if (map_dim > 0) {
965 for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
966 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
967 if (size_by_variant[variant] > 0) {
968 out << " " << geo_variant_name [variant] << "\t" << size_by_variant[variant] << endl;
969 }
970 }
971 }
972 if (map_dim >= 3 && size_by_variant[reference_element__t] != 0) {
973 out << " triangles\t" << size_by_variant[reference_element__t] << endl;
974 }
975 if (map_dim >= 3 && size_by_variant[reference_element__q] != 0) {
976 out << " quadrangles\t" << size_by_variant[reference_element__q] << endl;
977 }
978 if (map_dim >= 2 && size_by_variant[reference_element__e] != 0) {
979 out << " edges\t" << size_by_variant[reference_element__e] << endl;
980 }
981 out << "end header" << endl
982 << endl;
983 // nodes:
984 for (size_t inode = 0; inode < node.size(); inode++) {
985 for (size_t i = 0; i < dim; ++i) {
986 out << node[new2old_inode[inode]][i];
987 if (i+1 != dim) out << " ";
988 }
989 out << endl;
990 }
991 out << endl;
992 // elements:
993 for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
994 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
995 for (size_t ie = 0; ie < element.size(); ie++) {
996 if (element[ie].variant() != variant) continue;
997 if (!do_upgrade) {
998 if (element[ie].name() != 'e' || element[ie].order() > 1) {
999 out << element[ie].name() << "\t";
1000 }
1001 if (element[ie].order() > 1) {
1002 out << "p" << element[ie].order() << " ";
1003 }
1004 for (size_t iloc = 0, nloc = element[ie].size(); iloc < nloc; iloc++) {
1005 out << element[ie][iloc];
1006 if (iloc+1 != nloc) out << " ";
1007 }
1008 } else {
1009 // upgrade: truncate inodes from pk to p1
1010 out << element[ie].name() << "\t";
1011 for (size_t iloc = 0, nloc = element[ie].n_vertex(); iloc < nloc; iloc++) {
1012 out << element[ie][iloc];
1013 if (iloc+1 != nloc) out << " ";
1014 }
1015 }
1016 out << endl;
1017 }
1018 }
1019 out << endl;
1020 // faces: no-empty when upgrade & map_dim >= 3
1021 for (size_t variant = reference_element_first_variant_by_dimension(2);
1022 variant < reference_element_last_variant_by_dimension(2); variant++) {
1023 for (size_t ie = 0; ie < face.size(); ie++) {
1024 if (face[ie].variant() != variant) continue;
1025 out << face[ie].name() << "\t";
1026 for (size_t iloc = 0, nloc = face[ie].n_vertex(); iloc < nloc; iloc++) {
1027 out << face[ie][iloc];
1028 if (iloc+1 != nloc) out << " ";
1029 }
1030 out << endl;
1031 }
1032 }
1033 // edges: no-empty when upgrade & map_dim >= 2
1034 for (size_t ie = 0, ne = edge.size(); ie < ne; ++ie) {
1035 out << "e\t" << old2new_inode[edge[ie][0]]
1036 << " " << old2new_inode[edge[ie][1]] << endl;
1037 }
1038 out << endl;
1039 //
1040 // 5.2. output domains
1041 //
1042 for (size_t d = 0; d <= 3; ++d) {
1043 if (!do_upgrade) {
1044 put_domain_noupgrade (out, map_dim, d, domain_map[d], phys, element);
1045 } else {
1046 put_domain_upgrade (out, map_dim, d, domain_map[d], phys, element,
1047 edge, edge_ball, face, face_ball);
1048 }
1049 }
1050}
1051void usage()
1052{
1053 cerr << "msh2geo: usage:" << endl
1054 << "msh2geo [-[no]upgrade][-rz|-zr] in.msh > out.geo" << endl
1055 << "msh2geo [-[no]upgrade][-rz|-zr] < in.msh > out.geo" << endl;
1056 exit (1);
1057}
1058int main (int argc, char**argv) {
1059 // ----------------------------
1060 // scan the command line
1061 // ----------------------------
1062 string input_filename = "";
1063 std::string sys_coord_name = "cartesian";
1064 bool do_upgrade = true; // upgrade: put sides & truncate pk inodes to p1
1065 for (int i = 1; i < argc; i++) {
1066 if (strcmp (argv[i], "-rz") == 0) sys_coord_name = "rz";
1067 else if (strcmp (argv[i], "-zr") == 0) sys_coord_name = "zr";
1068 else if (strcmp (argv[i], "-cartesian") == 0) sys_coord_name = "cartesian";
1069 else if (strcmp (argv[i], "-upgrade") == 0) do_upgrade = true;
1070 else if (strcmp (argv[i], "-noupgrade") == 0) do_upgrade = false;
1071 else if (argv[i][0] == '-') {
1072 cerr << "field: unknown option `" << argv[i] << "'" << endl;
1073 usage();
1074 } else {
1075 input_filename = argv[i];
1076 }
1077 }
1078 if (input_filename == "") {
1079 msh2geo (cin, cout, sys_coord_name, do_upgrade);
1080 } else {
1081 ifstream fin (input_filename.c_str());
1082 if (!fin.good()) {
1083 cerr << "msh2geo: unable to read file \""<<input_filename<<"\"" << endl; exit (1);
1084 }
1085 msh2geo (fin, cout, sys_coord_name, do_upgrade);
1086 }
1087}
see the edge page for the full documentation
see the hexahedron page for the full documentation
see the point page for the full documentation
see the prism page for the full documentation
see the quadrangle page for the full documentation
void setorientation(size_t orient)
Definition msh2geo.cc:215
void setindex(size_t index)
Definition msh2geo.cc:214
void setshift(size_t shift)
Definition msh2geo.cc:216
see the geo_element page for the full documentation
size_t face(size_t i) const
Definition msh2geo.cc:247
array< geo_element_indirect, 12 > _edge
Definition msh2geo.cc:259
const geo_element_indirect & face_indirect(size_t i) const
Definition msh2geo.cc:243
size_t variant() const
Definition msh2geo.cc:239
void setindex(size_t index)
Definition msh2geo.cc:231
size_t edge(size_t i) const
Definition msh2geo.cc:245
size_type face(size_type i) const
geo_element_indirect & edge_indirect(size_t i)
Definition msh2geo.cc:233
size_type size() const
size_t index() const
Definition msh2geo.cc:237
void setorder(size_t order)
Definition msh2geo.cc:230
size_t dimension() const
Definition msh2geo.cc:240
const geo_element_indirect & edge_indirect(size_t i) const
Definition msh2geo.cc:242
void setname(char name)
Definition msh2geo.cc:229
size_t edge_local_vertex(size_t iedg, size_t edg_iv) const
Definition msh2geo.cc:263
size_t order() const
Definition msh2geo.cc:236
const geo_element_indirect & edge_indirect(size_type i) const
size_t n_face() const
Definition msh2geo.cc:248
void setgmshtype(size_t gmshtype)
Definition msh2geo.cc:232
array< geo_element_indirect, 6 > _face
Definition msh2geo.cc:260
size_t face_local_vertex(size_t iedg, size_t edg_iv) const
Definition msh2geo.cc:276
size_t face_n_vertex(size_t loc_ifac) const
Definition msh2geo.cc:288
variant_type variant() const
geo_element_indirect & face_indirect(size_t i)
Definition msh2geo.cc:234
size_t n_edge() const
Definition msh2geo.cc:244
size_t gmshtype() const
Definition msh2geo.cc:238
size_t n_vertex() const
Definition msh2geo.cc:241
size_type order() const
void inplace_intersection(const index_set &b)
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition msh2geo.cc:136
std::array< T, 3 > base
Definition msh2geo.cc:135
point_basic< Float > point
Definition point.h:163
int shift_type
Definition msh2geo.cc:191
int orientation_type
Definition msh2geo.cc:190
constexpr size_t reference_element__max_variant
const char _reference_element_name[reference_element__max_variant]
see the tetrahedron page for the full documentation
see the triangle page for the full documentation
#define assert_macro(ok_condition, message)
Definition dis_macros.h:113
#define error_macro(message)
Definition dis_macros.h:49
edge - reference element
const size_t dimension
Definition edge.icc:64
const size_t n_vertex
Definition edge.icc:66
int main()
Definition field2bb.cc:58
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
hexahedron - reference element
const size_t face[n_face][4]
void numbering_Pk_dis_idof(size_t order, const array< size_t, reference_element__max_variant > &size_by_variant, const array< size_t, reference_element__max_variant+1 > &first_by_variant, const array< size_t, reference_element__max_variant > &loc_ndof_by_variant, const geo_element &K, vector< size_t > &dis_idof1)
Definition msh2geo.cc:313
void usage()
Definition msh2geo.cc:1051
void put_domain_noupgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element)
Definition msh2geo.cc:396
void msh2geo(istream &in, ostream &out, string sys_coord_name, bool do_upgrade)
Definition msh2geo.cc:513
void put_domain_upgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element, const vector< geo_element_side > &edge, const vector< index_set > &edge_ball, const vector< geo_element_side > &face, const vector< index_set > &face_ball)
Definition msh2geo.cc:434
This file is part of Rheolef.
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition scatch.icc:44
STL namespace.
prism - reference element
quadrangle - reference element
tetrahedron - reference element
triangle - reference element