Mercurial > fem-fenics-eugenio

comparison src/Mesh.cc @ 268:61830a4f9ab9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Improve formatting
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Thu, 14 Aug 2014 12:26:55 +0200
parents a61fc34334ca
children
comparison
equal deleted inserted replaced
267:53039ac90368 268:61830a4f9ab9
27 std::size_t const); 27 std::size_t const);
28 octave_value 28 octave_value
29 compute_cell_markers (dolfin::Mesh const &, Array <octave_idx_type> const &, 29 compute_cell_markers (dolfin::Mesh const &, Array <octave_idx_type> const &,
30 std::size_t const); 30 std::size_t const);
31 31
32 DEFUN_DLD (Mesh, args, nargout,"-*- texinfo -*-\n\ 32 DEFUN_DLD (Mesh, args, nargout, "-*- texinfo -*-\n\
33 @deftypefn {Function File} {[@var{mesh_out}, \ 33 @deftypefn {Function File} {[@var{mesh_out}, \
34 @var{facet_markers}, @var{cell_markers}]} = \ 34 @var{facet_markers}, @var{cell_markers}]} = \
35 Mesh (@var{mesh_in}) \n\ 35 Mesh (@var{mesh_in}) \n\
36 Construct a mesh from file or from (p, e, t) format.\n\ 36 Construct a mesh from file or from (p, e, t) format.\n\
37 The @var{mesh_in} should be either\n\ 37 The @var{mesh_in} should be either\n\
53 { 53 {
54 int nargin = args.length (); 54 int nargin = args.length ();
55 octave_value_list retval; 55 octave_value_list retval;
56 56
57 if (nargin < 1 || nargin > 1 || nargout > 3) 57 if (nargin < 1 || nargin > 1 || nargout > 3)
58 print_usage (); 58 { print_usage (); }
59 else 59 else
60 { 60 {
61 if (!error_state) 61 if (!error_state)
62 { 62 {
63 63
100 dolfin::Mesh const & msh = msh_ov.get_msh (); 100 dolfin::Mesh const & msh = msh_ov.get_msh ();
101 std::size_t const D = p.rows (); 101 std::size_t const D = p.rows ();
102 102
103 retval(1) = compute_facet_markers (msh, e, D); 103 retval(1) = compute_facet_markers (msh, e, D);
104 if (nargout == 3) 104 if (nargout == 3)
105 retval(2) = compute_cell_markers (msh, t, D); 105 { retval(2) = compute_cell_markers (msh, t, D); }
106 } 106 }
107 } 107 }
108 } 108 }
109 109
110 else 110 else
111 error ("Mesh: the argument you provide is invalid"); 111 { error ("Mesh: the argument you provide is invalid"); }
112 112
113 } 113 }
114 } 114 }
115 return retval; 115 return retval;
116 } 116 }
117 117
118 mesh::mesh (Array<double>& p, 118 mesh::mesh (Array<double> & p,
119 Array<octave_idx_type>& e, 119 Array<octave_idx_type> & e,
120 Array<octave_idx_type>& t) 120 Array<octave_idx_type> & t)
121 { 121 {
122 std::size_t D = p.rows (); 122 std::size_t D = p.rows ();
123 if (D < 2 || D > 3) 123 if (D < 2 || D > 3)
124 error ("Mesh constructor: only 2D or 3D meshes are supported"); 124 { error ("Mesh constructor: only 2D or 3D meshes are supported"); }
125 else 125 else
126 { 126 {
127 dolfin::MeshEditor editor; 127 dolfin::MeshEditor editor;
128 SHARED_PTR <dolfin::Mesh> msh (new dolfin::Mesh ()); 128 SHARED_PTR <dolfin::Mesh> msh (new dolfin::Mesh ());
129 editor.open (*msh, D, D); 129 editor.open (*msh, D, D);
162 162
163 editor.close (); 163 editor.close ();
164 164
165 // store information associated with e 165 // store information associated with e
166 msh->init (D - 1); 166 msh->init (D - 1);
167 dolfin::MeshValueCollection<std::size_t> facet(*msh, D - 1); 167 dolfin::MeshValueCollection<std::size_t> facet (*msh, D - 1);
168 std::size_t num_side_edges = e.cols (); 168 std::size_t num_side_edges = e.cols ();
169 169
170 unsigned const size = 170 unsigned const size =
171 #ifdef LATEST_DOLFIN 171 #ifdef LATEST_DOLFIN
172 dolfin::MPI::size (MPI_COMM_WORLD); 172 dolfin::MPI::size (MPI_COMM_WORLD);
221 && (*f).entities(0)[2] == e.xelem (1, i) - 1 221 && (*f).entities(0)[2] == e.xelem (1, i) - 1
222 || (*f).entities(0)[0] == e.xelem (2, i) - 1 222 || (*f).entities(0)[0] == e.xelem (2, i) - 1
223 && (*f).entities(0)[1] == e.xelem (1, i) - 1 223 && (*f).entities(0)[1] == e.xelem (1, i) - 1
224 && (*f).entities(0)[2] == e.xelem (0, i) - 1) 224 && (*f).entities(0)[2] == e.xelem (0, i) - 1)
225 { 225 {
226 std::pair <std::size_t, std::size_t> 226 std::pair <std::size_t, std::size_t>
227 idxvl ((*f).index (), e.xelem (9, i)); 227 idxvl ((*f).index (), e.xelem (9, i));
228 msh->domains ().set_marker (idxvl, D - 1); 228 msh->domains ().set_marker (idxvl, D - 1);
229 break; 229 break;
230 } 230 }
231 } 231 }
398 octave_value 398 octave_value
399 compute_facet_markers (dolfin::Mesh const & _msh, 399 compute_facet_markers (dolfin::Mesh const & _msh,
400 Array <octave_idx_type> const & e, 400 Array <octave_idx_type> const & e,
401 std::size_t const D) 401 std::size_t const D)
402 { 402 {
403 dolfin::MeshFunction <std::size_t> facet (_msh, D-1, 403 dolfin::MeshFunction <std::size_t> facet (_msh, D - 1,
404 std::numeric_limits <std::size_t> 404 std::numeric_limits <std::size_t>
405 ::max ()); 405 ::max ());
406 std::size_t const num_side_edges = e.cols (); 406 std::size_t const num_side_edges = e.cols ();
407 std::vector <std::size_t> const & global_vertices = 407 std::vector <std::size_t> const & global_vertices =
408 _msh.topology ().global_indices (0); 408 _msh.topology ().global_indices (0);
409 std::vector <std::size_t> const & global_facets = 409 std::vector <std::size_t> const & global_facets =
410 _msh.topology ().global_indices (D-1); 410 _msh.topology ().global_indices (D - 1);
411 411
412 if (D == 2) 412 if (D == 2)
413 { 413 {
414 for (std::size_t i = 0; i < num_side_edges; ++i) 414 for (std::size_t i = 0; i < num_side_edges; ++i)
415 { 415 {
416 std::size_t local_vertex = 0; 416 std::size_t local_vertex = 0;
417 std::size_t const e_index = e.xelem (0, i) - 1; 417 std::size_t const e_index = e.xelem (0, i) - 1;
418 while (local_vertex < global_vertices.size () && 418 while (local_vertex < global_vertices.size () &&
419 e_index != global_vertices[local_vertex]) 419 e_index != global_vertices[local_vertex])
420 ++local_vertex; 420 { ++local_vertex; }
421 421
422 if (local_vertex < global_vertices.size ()) 422 if (local_vertex < global_vertices.size ())
423 { 423 {
424 dolfin::Vertex v (_msh, local_vertex); 424 dolfin::Vertex v (_msh, local_vertex);
425 for (dolfin::FacetIterator f (v); ! f.end (); ++f) 425 for (dolfin::FacetIterator f (v); ! f.end (); ++f)
448 { 448 {
449 std::size_t local_vertex = 0; 449 std::size_t local_vertex = 0;
450 std::size_t const e_index = e.xelem (0, i) - 1; 450 std::size_t const e_index = e.xelem (0, i) - 1;
451 while (local_vertex < global_vertices.size () && 451 while (local_vertex < global_vertices.size () &&
452 e_index != global_vertices[local_vertex]) 452 e_index != global_vertices[local_vertex])
453 ++local_vertex; 453 { ++local_vertex; }
454 454
455 if (local_vertex < global_vertices.size ()) 455 if (local_vertex < global_vertices.size ())
456 { 456 {
457 dolfin::Vertex v (_msh, local_vertex); 457 dolfin::Vertex v (_msh, local_vertex);
458 for (dolfin::FacetIterator f (v); ! f.end (); ++f) 458 for (dolfin::FacetIterator f (v); ! f.end (); ++f)
462 std::size_t const & vertex1 = 462 std::size_t const & vertex1 =
463 global_vertices[f->entities(0)[1]]; 463 global_vertices[f->entities(0)[1]];
464 std::size_t const & vertex2 = 464 std::size_t const & vertex2 =
465 global_vertices[f->entities(0)[2]]; 465 global_vertices[f->entities(0)[2]];
466 466
467 if (vertex0 == e(0, i) - 1 467 if (vertex0 == e.xelem (0, i) - 1
468 && vertex1 == e.xelem (1, i) - 1 468 && vertex1 == e.xelem (1, i) - 1
469 && vertex2 == e.xelem (2, i) - 1 469 && vertex2 == e.xelem (2, i) - 1
470 || vertex0 == e.xelem (0, i) - 1 470 || vertex0 == e.xelem (0, i) - 1
471 && vertex1 == e.xelem (2, i) - 1 471 && vertex1 == e.xelem (2, i) - 1
472 && vertex2 == e.xelem (1, i) - 1 472 && vertex2 == e.xelem (1, i) - 1
515 { 515 {
516 std::size_t local_vertex = 0; 516 std::size_t local_vertex = 0;
517 std::size_t const t_index = t.xelem (0, i) - 1; 517 std::size_t const t_index = t.xelem (0, i) - 1;
518 while (local_vertex < global_vertices.size () && 518 while (local_vertex < global_vertices.size () &&
519 t_index != global_vertices[local_vertex]) 519 t_index != global_vertices[local_vertex])
520 ++local_vertex; 520 { ++local_vertex; }
521 521
522 if (local_vertex < global_vertices.size ()) 522 if (local_vertex < global_vertices.size ())
523 { 523 {
524 dolfin::Vertex v (_msh, local_vertex); 524 dolfin::Vertex v (_msh, local_vertex);
525 for (dolfin::CellIterator f (v); ! f.end (); ++f) 525 for (dolfin::CellIterator f (v); ! f.end (); ++f)
564 { 564 {
565 std::size_t local_vertex = 0; 565 std::size_t local_vertex = 0;
566 std::size_t const t_index = t.xelem (0, i) - 1; 566 std::size_t const t_index = t.xelem (0, i) - 1;
567 while (local_vertex < global_vertices.size () && 567 while (local_vertex < global_vertices.size () &&
568 t_index != global_vertices[local_vertex]) 568 t_index != global_vertices[local_vertex])
569 ++local_vertex; 569 { ++local_vertex; }
570 570
571 if (local_vertex < global_vertices.size ()) 571 if (local_vertex < global_vertices.size ())
572 { 572 {
573 dolfin::Vertex v (_msh, local_vertex); 573 dolfin::Vertex v (_msh, local_vertex);
574 for (dolfin::CellIterator f (v); ! f.end (); ++f) 574 for (dolfin::CellIterator f (v); ! f.end (); ++f)
581 global_vertices[f->entities(0)[2]]; 581 global_vertices[f->entities(0)[2]];
582 std::size_t const & vertex3 = 582 std::size_t const & vertex3 =
583 global_vertices[f->entities(0)[3]]; 583 global_vertices[f->entities(0)[3]];
584 584
585 if (vertex0 == t.xelem (0, i) - 1 585 if (vertex0 == t.xelem (0, i) - 1
586 && vertex1 == t.xelem (1, i) - 1 586 && vertex1 == t.xelem (1, i) - 1
587 && vertex2 == t.xelem (2, i) - 1 587 && vertex2 == t.xelem (2, i) - 1
588 && vertex3 == t.xelem (3, i) - 1 588 && vertex3 == t.xelem (3, i) - 1
589 || vertex0 == t.xelem (0, i) - 1 589 || vertex0 == t.xelem (0, i) - 1
590 && vertex1 == t.xelem (1, i) - 1 590 && vertex1 == t.xelem (1, i) - 1
591 && vertex2 == t.xelem (3, i) - 1 591 && vertex2 == t.xelem (3, i) - 1
592 && vertex3 == t.xelem (2, i) - 1 592 && vertex3 == t.xelem (2, i) - 1
593 || vertex0 == t.xelem (0, i) - 1 593 || vertex0 == t.xelem (0, i) - 1
594 && vertex1 == t.xelem (2, i) - 1 594 && vertex1 == t.xelem (2, i) - 1
595 && vertex2 == t.xelem (1, i) - 1 595 && vertex2 == t.xelem (1, i) - 1
596 && vertex3 == t.xelem (3, i) - 1 596 && vertex3 == t.xelem (3, i) - 1
597 || vertex0 == t.xelem (0, i) - 1 597 || vertex0 == t.xelem (0, i) - 1
598 && vertex1 == t.xelem (2, i) - 1 598 && vertex1 == t.xelem (2, i) - 1
599 && vertex2 == t.xelem (3, i) - 1 599 && vertex2 == t.xelem (3, i) - 1
600 && vertex3 == t.xelem (1, i) - 1 600 && vertex3 == t.xelem (1, i) - 1
601 || vertex0 == t.xelem (0, i) - 1 601 || vertex0 == t.xelem (0, i) - 1
602 && vertex1 == t.xelem (3, i) - 1 602 && vertex1 == t.xelem (3, i) - 1
603 && vertex2 == t.xelem (1, i) - 1 603 && vertex2 == t.xelem (1, i) - 1
604 && vertex3 == t.xelem (2, i) - 1 604 && vertex3 == t.xelem (2, i) - 1
605 || vertex0 == t.xelem (0, i) - 1 605 || vertex0 == t.xelem (0, i) - 1
606 && vertex1 == t.xelem (3, i) - 1 606 && vertex1 == t.xelem (3, i) - 1
607 && vertex2 == t.xelem (2, i) - 1 607 && vertex2 == t.xelem (2, i) - 1
608 && vertex3 == t.xelem (1, i) - 1 608 && vertex3 == t.xelem (1, i) - 1
609 609
610 || vertex0 == t.xelem (1, i) - 1 610 || vertex0 == t.xelem (1, i) - 1
611 && vertex1 == t.xelem (0, i) - 1 611 && vertex1 == t.xelem (0, i) - 1
612 && vertex2 == t.xelem (2, i) - 1 612 && vertex2 == t.xelem (2, i) - 1
613 && vertex3 == t.xelem (3, i) - 1 613 && vertex3 == t.xelem (3, i) - 1
614 || vertex0 == t.xelem (1, i) - 1 614 || vertex0 == t.xelem (1, i) - 1
615 && vertex1 == t.xelem (0, i) - 1 615 && vertex1 == t.xelem (0, i) - 1
616 && vertex2 == t.xelem (3, i) - 1 616 && vertex2 == t.xelem (3, i) - 1
617 && vertex3 == t.xelem (2, i) - 1 617 && vertex3 == t.xelem (2, i) - 1
618 || vertex0 == t.xelem (1, i) - 1 618 || vertex0 == t.xelem (1, i) - 1
619 && vertex1 == t.xelem (2, i) - 1 619 && vertex1 == t.xelem (2, i) - 1
620 && vertex2 == t.xelem (0, i) - 1 620 && vertex2 == t.xelem (0, i) - 1
621 && vertex3 == t.xelem (3, i) - 1 621 && vertex3 == t.xelem (3, i) - 1
622 || vertex0 == t.xelem (1, i) - 1 622 || vertex0 == t.xelem (1, i) - 1
623 && vertex1 == t.xelem (2, i) - 1 623 && vertex1 == t.xelem (2, i) - 1
624 && vertex2 == t.xelem (3, i) - 1 624 && vertex2 == t.xelem (3, i) - 1
625 && vertex3 == t.xelem (0, i) - 1 625 && vertex3 == t.xelem (0, i) - 1
626 || vertex0 == t.xelem (1, i) - 1 626 || vertex0 == t.xelem (1, i) - 1
627 && vertex1 == t.xelem (3, i) - 1 627 && vertex1 == t.xelem (3, i) - 1
628 && vertex2 == t.xelem (0, i) - 1 628 && vertex2 == t.xelem (0, i) - 1
629 && vertex3 == t.xelem (2, i) - 1 629 && vertex3 == t.xelem (2, i) - 1
630 || vertex0 == t.xelem (1, i) - 1 630 || vertex0 == t.xelem (1, i) - 1
631 && vertex1 == t.xelem (3, i) - 1 631 && vertex1 == t.xelem (3, i) - 1
632 && vertex2 == t.xelem (2, i) - 1 632 && vertex2 == t.xelem (2, i) - 1
633 && vertex3 == t.xelem (0, i) - 1 633 && vertex3 == t.xelem (0, i) - 1
634 634
635 || vertex0 == t.xelem (2, i) - 1 635 || vertex0 == t.xelem (2, i) - 1
636 && vertex1 == t.xelem (0, i) - 1 636 && vertex1 == t.xelem (0, i) - 1
637 && vertex2 == t.xelem (1, i) - 1 637 && vertex2 == t.xelem (1, i) - 1
638 && vertex3 == t.xelem (3, i) - 1 638 && vertex3 == t.xelem (3, i) - 1
639 || vertex0 == t.xelem (2, i) - 1 639 || vertex0 == t.xelem (2, i) - 1
640 && vertex1 == t.xelem (0, i) - 1 640 && vertex1 == t.xelem (0, i) - 1
641 && vertex2 == t.xelem (3, i) - 1 641 && vertex2 == t.xelem (3, i) - 1
642 && vertex3 == t.xelem (1, i) - 1 642 && vertex3 == t.xelem (1, i) - 1
643 || vertex0 == t.xelem (2, i) - 1 643 || vertex0 == t.xelem (2, i) - 1
644 && vertex1 == t.xelem (1, i) - 1 644 && vertex1 == t.xelem (1, i) - 1
645 && vertex2 == t.xelem (0, i) - 1 645 && vertex2 == t.xelem (0, i) - 1
646 && vertex3 == t.xelem (3, i) - 1 646 && vertex3 == t.xelem (3, i) - 1
647 || vertex0 == t.xelem (2, i) - 1 647 || vertex0 == t.xelem (2, i) - 1
648 && vertex1 == t.xelem (1, i) - 1 648 && vertex1 == t.xelem (1, i) - 1
649 && vertex2 == t.xelem (3, i) - 1 649 && vertex2 == t.xelem (3, i) - 1
650 && vertex3 == t.xelem (0, i) - 1 650 && vertex3 == t.xelem (0, i) - 1
651 || vertex0 == t.xelem (2, i) - 1 651 || vertex0 == t.xelem (2, i) - 1
652 && vertex1 == t.xelem (3, i) - 1 652 && vertex1 == t.xelem (3, i) - 1
653 && vertex2 == t.xelem (0, i) - 1 653 && vertex2 == t.xelem (0, i) - 1
654 && vertex3 == t.xelem (1, i) - 1 654 && vertex3 == t.xelem (1, i) - 1
655 || vertex0 == t.xelem (2, i) - 1 655 || vertex0 == t.xelem (2, i) - 1
656 && vertex1 == t.xelem (3, i) - 1 656 && vertex1 == t.xelem (3, i) - 1
657 && vertex2 == t.xelem (1, i) - 1 657 && vertex2 == t.xelem (1, i) - 1
658 && vertex3 == t.xelem (0, i) - 1 658 && vertex3 == t.xelem (0, i) - 1
659 659
660 || vertex0 == t.xelem (3, i) - 1 660 || vertex0 == t.xelem (3, i) - 1
661 && vertex1 == t.xelem (0, i) - 1 661 && vertex1 == t.xelem (0, i) - 1
662 && vertex2 == t.xelem (1, i) - 1 662 && vertex2 == t.xelem (1, i) - 1
663 && vertex3 == t.xelem (2, i) - 1 663 && vertex3 == t.xelem (2, i) - 1
664 || vertex0 == t.xelem (3, i) - 1 664 || vertex0 == t.xelem (3, i) - 1
665 && vertex1 == t.xelem (0, i) - 1 665 && vertex1 == t.xelem (0, i) - 1
666 && vertex2 == t.xelem (2, i) - 1 666 && vertex2 == t.xelem (2, i) - 1
667 && vertex3 == t.xelem (1, i) - 1 667 && vertex3 == t.xelem (1, i) - 1
668 || vertex0 == t.xelem (3, i) - 1 668 || vertex0 == t.xelem (3, i) - 1
669 && vertex1 == t.xelem (1, i) - 1 669 && vertex1 == t.xelem (1, i) - 1
670 && vertex2 == t.xelem (0, i) - 1 670 && vertex2 == t.xelem (0, i) - 1
671 && vertex3 == t.xelem (2, i) - 1 671 && vertex3 == t.xelem (2, i) - 1
672 || vertex0 == t.xelem (3, i) - 1 672 || vertex0 == t.xelem (3, i) - 1
673 && vertex1 == t.xelem (1, i) - 1 673 && vertex1 == t.xelem (1, i) - 1
674 && vertex2 == t.xelem (2, i) - 1 674 && vertex2 == t.xelem (2, i) - 1
675 && vertex3 == t.xelem (0, i) - 1 675 && vertex3 == t.xelem (0, i) - 1
676 || vertex0 == t.xelem (3, i) - 1 676 || vertex0 == t.xelem (3, i) - 1
677 && vertex1 == t.xelem (2, i) - 1 677 && vertex1 == t.xelem (2, i) - 1
678 && vertex2 == t.xelem (0, i) - 1 678 && vertex2 == t.xelem (0, i) - 1
679 && vertex3 == t.xelem (1, i) - 1 679 && vertex3 == t.xelem (1, i) - 1
680 || vertex0 == t.xelem (3, i) - 1 680 || vertex0 == t.xelem (3, i) - 1
681 && vertex1 == t.xelem (2, i) - 1 681 && vertex1 == t.xelem (2, i) - 1
682 && vertex2 == t.xelem (1, i) - 1 682 && vertex2 == t.xelem (1, i) - 1
683 && vertex3 == t.xelem (0, i) - 1) 683 && vertex3 == t.xelem (0, i) - 1)
684 { 684 {
685 cell[*f] = t.xelem (4, i); 685 cell[*f] = t.xelem (4, i);
686 break; 686 break;
687 } 687 }
688 } 688 }