← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/tst/form_grad_grad_s_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/tst/form_grad_grad_s_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
//
 
22
// test des formes grad_grad_s et mass_s
 
23
//
 
24
#include "rheolef.h"
 
25
using namespace std;
 
26
using namespace rheolef;
 
27
// --------------------------------------------------------------------------------
 
28
// definition de phi
 
29
// --------------------------------------------------------------------------------
 
30
point origin (0.5,0.5,0);
 
31
point normal (0,1,0);
 
32
 
 
33
Float levelset_function (const point& x) {
 
34
    return dot(normal, x-origin); // definition de la fonction level set//
 
35
}
 
36
// test
 
37
string i_test = "1";
 
38
Float u (const point& x) {
 
39
    if (i_test == "1") return 1;
 
40
    if (i_test == "x") return x[0];
 
41
    if (i_test == "y") return x[1];
 
42
    if (i_test == "z") return x[2];
 
43
    if (i_test == "x2") return sqr(x[0]);
 
44
    if (i_test == "y2") return sqr(x[1]);
 
45
    if (i_test == "z2") return sqr(x[2]);
 
46
    if (i_test == "xy") return x[0]*x[1];
 
47
    if (i_test == "xz") return x[0]*x[2];
 
48
    if (i_test == "yz") return x[1]*x[2];
 
49
    error_macro ("invalid test " << i_test);
 
50
}
 
51
void save_as_matlab (const form& a, string name = "a") {
 
52
    string filename = "matrix-" + name + ".m";
 
53
    ofstream out (filename.c_str());
 
54
    cerr << "! file " << filename << " created." << endl;
 
55
    out << setprecision(15) 
 
56
        << name << " = " << matlab << a.uu << ";" << endl
 
57
        << "vp = eig(" << name <<");" << endl
 
58
        << "printf (\"vp = %26.16g\\n\", vp);" << endl
 
59
        << "d = det(" << name << ")" << endl
 
60
        ;
 
61
    out.close();
 
62
}
 
63
void save_as_matlab (const field& u, string name = "u") {
 
64
    string filename = "vector-" + name + ".m";
 
65
    ofstream out (filename.c_str());
 
66
    cerr << "! file " << filename << " created." << endl;
 
67
    out << setprecision(15) 
 
68
        << name << " = " << matlab << u.u << ";" << endl
 
69
        ;
 
70
    out.close();
 
71
}
 
72
void 
 
73
usage () {
 
74
  cerr << "usage: grad_grad_tst " 
 
75
       << "mesh.geo " 
 
76
       << "-normal x y z" 
 
77
       << "-origin x y z" 
 
78
       << "-a val " 
 
79
       << "-m val " 
 
80
       << "-test poly " 
 
81
       << endl;
 
82
  exit (1);
 
83
}
 
84
int main (int argc, char**argv)
 
85
{
 
86
    Float a_valid = 0;
 
87
    Float m_valid = -1;
 
88
    string filename;
 
89
    for (int i = 1; i < argc; i++) {
 
90
 
 
91
        if (strcmp (argv[i], "-test") == 0)   {
 
92
            if (i+1 == argc) usage();
 
93
            i_test = argv[++i];
 
94
        } else if (strcmp (argv[i], "-a") == 0)   {
 
95
            if (i+1 == argc || !is_float(argv[i+1])) usage();
 
96
            a_valid = to_float (argv[++i]);
 
97
        } else if (strcmp (argv[i], "-m") == 0)   {
 
98
            if (i+1 == argc || !is_float(argv[i+1])) usage();
 
99
            m_valid = to_float (argv[++i]);
 
100
        } else if ((strcmp (argv[i], "-origin") == 0) || (strcmp (argv[i], "-normal") == 0))   {
 
101
            point x;
 
102
            unsigned int io = i;
 
103
            if (i+1 == argc || !is_float(argv[i+1])) {
 
104
                warning_macro ("invalid argument to `" << argv[i] << "'");
 
105
                usage();
 
106
            }
 
107
            x[0] = to_float (argv[++i]);
 
108
            if (i+1 < argc && is_float(argv[i+1])) {
 
109
                x[1] = to_float (argv[++i]);
 
110
                if (i+1 < argc && is_float(argv[i+1])) {
 
111
                    x[2] = to_float (argv[++i]);
 
112
                }       
 
113
            }   
 
114
            if (strcmp (argv[io], "-origin") == 0)   {
 
115
               origin = x;
 
116
            } else {
 
117
               normal = x;
 
118
            }
 
119
        } else {
 
120
            // input on file
 
121
            filename = argv[i];
 
122
        }
 
123
    }
 
124
    if (filename == "") usage();
 
125
    geo mesh(filename);
 
126
    space Vh (mesh, "P1");
 
127
    field phi_h_lambda = interpolate(Vh, levelset_function);
 
128
    geo band = banded_level_set (phi_h_lambda);
 
129
    if (mesh.dimension() == 3) {
 
130
        cout << "$ DATA = CONTCURVE" << endl
 
131
             << "% equalscale" << endl
 
132
             << "%contstyle=2 # fill color painting" << endl
 
133
             << "%meshplot" << endl
 
134
             << "%contlabel=false # write labels when using line style contours" << endl;
 
135
    }
 
136
    phi_h_lambda.save("phi");
 
137
    space Bh (band, "P1");
 
138
    field phi_h = interpolate(Bh, levelset_function);
 
139
    form a = form (Bh, Bh, "grad_grad_s", phi_h);
 
140
    form m = form (Bh, Bh, "mass_s",      phi_h);
 
141
    save_as_matlab (a, "a");
 
142
    save_as_matlab (m, "m");
 
143
    save_as_matlab (phi_h, "phi");
 
144
    field uh = interpolate (Bh, u);
 
145
    Float a_tst = a(uh,uh);
 
146
    Float m_tst = m(uh,uh);
 
147
    cerr << setprecision(15);
 
148
    cerr << "m(uh,uh)=" << m_tst << endl;
 
149
    cerr << "a(uh,uh)=" << a_tst << endl;
 
150
    Float tol = 1e-7;
 
151
    bool status = (fabs(a_tst - a_valid) < tol) && (m_valid < 0 || fabs(m_tst - m_valid) < tol) ;
 
152
    if (! status) {
 
153
      if (                fabs(a_tst - a_valid) >= tol) cerr << "ERROR: a(uh,uh)=" << a_tst << " != " << a_valid << endl;
 
154
      if (m_valid >= 0 && fabs(m_tst - m_valid) >= tol) cerr << "ERROR: m(uh,uh)=" << m_tst << " != " << m_valid << endl;
 
155
    }
 
156
    return (status ? 0 : 1); 
 
157
}