← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« 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, Pierre Saramito, Sylvestre Ledru
  • Date: 2012-05-14 14:02:09 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20120514140209-dzbdlidkotyflf9e
Tags: 6.1-1
http://bugs.debian.org/671996
[ Pierre Saramito ]
* New upstream release 6.1 (minor changes):
  - support arbitrarily polynomial order Pk
  - source code supports g++-4.7 (closes: #671996)

[ Sylvestre Ledru ]
* update of the watch file

Show diffs side-by-side

added added

removed removed

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