← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/ptst/geo_element_curved_ball_tst.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/ptst/geo_element_curved_ball_tst.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
#include "rheolef.h"
 
22
#include "rheolef/geo_element_curved_ball.h"
 
23
using namespace rheolef;
 
24
using namespace std;
 
25
void put (const geo_basic<Float,sequential>& omega, size_t order, bool check_tetra_with_curved_edge) {
 
26
  typedef geo_basic<Float,sequential>::size_type size_type;
 
27
  typedef point_basic<size_type>                 ilat;
 
28
  size_type d = omega.dimension();
 
29
  check_macro (d >= 2, "unexpected dimension "<<d);
 
30
  point center(0,0);
 
31
  Float radius = 1;
 
32
  bool have_quad = false;
 
33
  array<point,sequential> node = omega.get_nodes();
 
34
  cout << setprecision(7)
 
35
       << "set xrange ["<<omega.xmin()[0]<<":"<<omega.xmax()[0]<<"]" << endl
 
36
       << "set yrange ["<<omega.xmin()[1]<<":"<<omega.xmax()[1]<<"]" << endl
 
37
       ;
 
38
  if (d == 2) {
 
39
    cout << "set size square" << endl;
 
40
  } else {
 
41
    cout << "set zrange ["<<omega.xmin()[2]<<":"<<omega.xmax()[2]<<"]" << endl
 
42
         << "set xyplane at -0.1" << endl
 
43
         << "set view equal xyz # equal scales" << endl
 
44
         ;
 
45
  }
 
46
  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
 
47
    const geo_element& K = omega.get_geo_element(d, ie);
 
48
    switch (K.variant()) {
 
49
      case reference_element::t: {
 
50
        center = point(-1,-1); radius = sqrt(Float(5));
 
51
        curved_ball_t F (node[K[0]], node[K[1]], node[K[2]], 1, center, radius);
 
52
        for (size_type i = 1; i < order; i++) {
 
53
          for (size_type j = 1; i+j < order; j++) {
 
54
            size_type loc_inod = reference_element_t::ilat2loc_inod (order, ilat(i, j));
 
55
            point hat_x (1.*i/order, 1.*j/order);
 
56
            point x = F (hat_x);
 
57
            cout << "set label \"(" << i << "," << j << ")\" at " << x[0] << "," << x[1] << endl;
 
58
          }
 
59
        }
 
60
        break;
 
61
      }
 
62
      case reference_element::q: {
 
63
        have_quad = true;
 
64
        curved_ball_q F (node[K[1]], node[K[2]], node[K[3]], node[K[0]]);
 
65
        for (size_type i = 1; i < order; i++) {
 
66
          for (size_type j = 1; j < order; j++) {
 
67
            size_type loc_inod = reference_element_q::ilat2loc_inod (order, ilat(i, j));
 
68
            point hat_x (-1+2.*i/order, -1+2.*j/order);
 
69
            point x = F (hat_x);
 
70
            cout << "set label \"(" << i << "," << j << ")\" at " << x[0] << "," << x[1] << endl;
 
71
          }
 
72
        }
 
73
        break;
 
74
      }
 
75
      case reference_element::T: {
 
76
        //center = point(-1,-1,-1); radius = sqrt(Float(6));
 
77
        curved_ball_T F (node[K[0]], node[K[1]], node[K[2]], node[K[3]], center, radius);
 
78
        if (!check_tetra_with_curved_edge) {
 
79
          F.set_boundary_face (3);
 
80
        } else {
 
81
          F.set_boundary_edge (1);
 
82
        }
 
83
        for (size_type i = 1; i < order; i++) {
 
84
          for (size_type j = 1; i+j < order; j++) {
 
85
            for (size_type k = 1; i+j+k < order; k++) {
 
86
              size_type loc_inod = reference_element_T::ilat2loc_inod (order, ilat(i, j, k));
 
87
              point hat_x (1.*i/order, 1.*j/order, 1.*k/order);
 
88
              point x = F (hat_x);
 
89
              cout << "set label \"(" << i << "," << j << "," << k << ")\" at "
 
90
                   << x[0] << "," << x[1] << "," << x[2] << endl;
 
91
            }
 
92
          }
 
93
        }
 
94
        break;
 
95
      }
 
96
      default:
 
97
        error_macro ("unexpected variant "<<K.variant());
 
98
    }
 
99
  }
 
100
  cout << "r  = " << radius << endl
 
101
       << "cx = " << center[0] << endl
 
102
       << "cy = " << center[1] << endl
 
103
      ;
 
104
  if (d == 3) {
 
105
    cout << "cz = " << center[2] << endl;
 
106
  }
 
107
  if (d == 2) {
 
108
    if (have_quad) {
 
109
      cout << "bdry(x) = (x > 1./sqrt(2.)) ? (cy + sqrt(r**2-(x-cx)**2)) : (1 + ((1. - sqrt(2.))*x))" << endl;
 
110
    } else {
 
111
      cout << "bdry(x) = cy + sqrt(r**2-(x-cx)**2)" << endl;
 
112
    }
 
113
    cout << "plot bdry(x)" << endl;
 
114
  } else { // d=3
 
115
    cout << "bdry(x,y) = cz + sqrt(r**2-(x-cx)**2-(y-cy)**2)" << endl;
 
116
    cout << "splot bdry(x,y)" << endl;
 
117
  }
 
118
  cout << "pause -1 \"<retour>\"" << endl;
 
119
}
 
120
int main(int argc, char**argv) {
 
121
  environment rheolef (argc, argv);
 
122
  geo_basic<Float,sequential> omega (argv[1]);
 
123
  size_t order = (argc > 2) ? atoi(argv[2]) : omega.order();
 
124
  bool check_tetra_with_curved_edge = (argc > 3);
 
125
  put (omega, order, check_tetra_with_curved_edge);
 
126
}