← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rheolef/saucy-proposed

« back to all changes in this revision

Viewing changes to nfem/form_element/mass_s.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/form_element/mass_s.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
// mass matrix on a surface, as defined by a level set
 
23
//
 
24
#include "mass_s.h"
 
25
#include "rheolef/ublas_matrix_range.h"
 
26
 
 
27
namespace rheolef {
 
28
extern void compute_matrix_m (const geo_element& K, const ublas::vector<point>& x, const ublas::vector<Float>& f, ublas::matrix<Float>& mk);
 
29
extern void compute_matrix_extension (const geo_element& K, const ublas::vector<point>& x, const ublas::vector<Float>& f, ublas::matrix<Float>& mk);
 
30
extern bool belongs_to_band (const geo_element& K, const ublas::vector<Float>& f);
 
31
} // namespace rheolef
 
32
 
 
33
using namespace rheolef;
 
34
using namespace std;
 
35
using namespace ublas;
 
36
 
 
37
void 
 
38
mass_s::operator() (const geo_element& K, ublas::matrix<Float>& m) const
 
39
{
 
40
    size_type nloc = K.size();
 
41
    ublas::vector<point> x (nloc);
 
42
    ublas::vector<Float> f (nloc);
 
43
    tiny_vector<size_t> bgd_idof (nloc);
 
44
    const field& phi_h = _wh;
 
45
    space Bh;
 
46
    switch (get_XY_type()) {
 
47
      case Bh_Bh: // Bh is used only to get the mesh from, so no need to check whether first space is P0 or P1
 
48
      case Bh_Wh: Bh = get_first_space(); break;
 
49
      case Wh_Bh:
 
50
      default:    Bh = get_second_space(); break;
 
51
    }
 
52
    const geo& beta_h   = Bh.get_geo();
 
53
    const geo& lambda_h = Bh.get_global_geo();
 
54
    phi_h.get_space().set_global_dof (K, bgd_idof);
 
55
    for (size_t iloc = 0; iloc < nloc; iloc++) {
 
56
      x[iloc] = lambda_h.vertex(K[iloc]);
 
57
      f[iloc] = phi_h.at (bgd_idof[iloc]);
 
58
    }
 
59
    if (!belongs_to_band(K, f)) {
 
60
        // empty m matrix => no assembly peformed.
 
61
        m.resize(0,0);
 
62
        return;
 
63
    }
 
64
    if (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() &&
 
65
        get_second_space().dimension() == get_second_space().get_geo().map_dimension()) {
 
66
 
 
67
      // form ("mass_s", Bh, Bh, phi_h):
 
68
      m.resize (nloc, nloc);
 
69
      compute_matrix_m (K, x, f, m);
 
70
 
 
71
      // then aggregate matrix in case it's a P0 matrix which is wanted:
 
72
      if (get_second_space().get_approx()=="P0")
 
73
       {
 
74
         ublas::matrix<Float> ma (1, nloc, 0);
 
75
         for (size_t j=0; j<nloc; j++) 
 
76
          { 
 
77
//          ma(0, j)=0;
 
78
            for (size_t i=0; i<nloc; i++) ma(0, j) += m(i,j);
 
79
          }
 
80
         m.resize(1, nloc);
 
81
         for (size_t j=0; j<nloc; j++) m(0,j)=ma(0,j);
 
82
       }
 
83
 
 
84
      if (get_first_space().get_approx()=="P0")
 
85
       {
 
86
         size_t ni = m.size1();
 
87
         ublas::matrix<Float> ma (ni, 1, 0);
 
88
         for (size_t i=0; i<ni; i++) 
 
89
          { 
 
90
//          ma(0, j)=0;
 
91
            for (size_t j=0; j<nloc; j++) ma(i, 0) += m(i,j);
 
92
          }
 
93
         m.resize(ni, 1);
 
94
         for (size_t i=0; i<ni; i++) m(i,0)=ma(i,0);
 
95
       }
 
96
 
 
97
    } else if (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() +1 &&
 
98
               get_second_space().dimension() == get_second_space().get_geo().map_dimension()) {
 
99
 
 
100
      // form ("mass_s", Wh, Bh, phi_h): extension
 
101
      m.resize (nloc, nloc-1);
 
102
      compute_matrix_extension (K, x, f, m);
 
103
 
 
104
    } else if (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() &&
 
105
               get_second_space().dimension() == get_second_space().get_geo().map_dimension()+1) {
 
106
 
 
107
      // form ("mass_s", Bh, Wh, phi_h): projection
 
108
      ublas::matrix<Float> mt (nloc, nloc-1);
 
109
      compute_matrix_extension (K, x, f, mt);
 
110
      m.resize (nloc-1, nloc);
 
111
      for (size_t i = 0; i < nloc-1; i++) {
 
112
      for (size_t j = 0; j < nloc;   j++) {
 
113
        m(i,j) = mt(j,i);
 
114
      }}
 
115
    } else {
 
116
      error_macro ("unexpected dimensions and map_dimensions: not yet implemented.");
 
117
    }
 
118
}
 
119
void
 
120
mass_s::check_after_initialize () const
 
121
{
 
122
  // suppose also that multi-component spaces are homogeneous,
 
123
  // i.e. that all components have the same approx
 
124
  check_macro (get_first_space().n_component() == 1,
 
125
    "incompatible spaces for the `mass_s' form.");
 
126
  check_macro (get_first_space().n_component() == get_second_space().n_component(),
 
127
    "incompatible spaces for the `mass_s' form.");
 
128
  check_macro (get_first_space().get_approx() == "P1" || get_first_space().get_approx() == "P0", 
 
129
        "`mass_s' form: P0 or P1 approx required");
 
130
  check_macro (get_second_space().get_approx() == "P1" || get_second_space().get_approx() == "P0", 
 
131
        "`mass_s' form: P0 or P1 approx required");
 
132
  check_macro (is_weighted(), "`mass_s' form may specify the level set function");
 
133
 
 
134
  if        (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() +1 &&
 
135
             get_second_space().dimension() == get_second_space().get_geo().map_dimension()) {
 
136
       _XY_switch = Wh_Bh; // extension
 
137
  } else if (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() &&
 
138
             get_second_space().dimension() == get_second_space().get_geo().map_dimension()+1) {
 
139
       _XY_switch = Bh_Wh; // prolongement
 
140
  } else if (get_first_space().dimension()  == get_first_space().get_geo().map_dimension() &&
 
141
             get_second_space().dimension() == get_second_space().get_geo().map_dimension()) {
 
142
       _XY_switch = Bh_Bh; // classic banded level set
 
143
  } else {
 
144
       error_macro ("form on a banded level set: at least one of the two spaces may be a band");
 
145
  }
 
146
}