This file is indexed.

/usr/include/rheolef/field_vf_assembly.h is in librheolef-dev 6.5-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#ifndef _RHEO_FIELD_VF_ASSEMBLY_H
#define _RHEO_FIELD_VF_ASSEMBLY_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
#include "rheolef/field.h"
#include "rheolef/test.h"
#include "rheolef/quadrature.h"
/*
   let:
     l(v) = int_domain expr(v) dx
  
   The integrals are evaluated over each element K of the domain
   by using a quadrature formulae given by qopt
  
   expr(v) is a linear expression with respect to the test-function v

   The test-function v is replaced by each of the basis function of 
   the corresponding finite element space Xh: (phi_i), i=0..dim(Xh)-1
 
   The integrals over K are transformed on the reference element with
   the piola transformation:
     F : hat_K ---> K
         hat_x |--> x = F(hat_x)
 
   exemples:
   1) expr(v) = v
    int_K phi_i(x) dx 
     = int_{hat_K} hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
     = sum_q hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq

    The value(q,i) = (hat_phi_i(hat_xq)) are evaluated on time for all over the 
    reference element hat_K and for the given quadrature formulae by:
  	expr.initialize (dom, quad);
    This expression is represented by the 'test' class (see test.h)

   2) expr(v) = f*v
    int_K f(x)*phi_i(x) dx 
     = int_{hat_K} f(F(hat_x)*hat_phi_i(hat_x) det(DF(hat_x)) d hat_x
     = sum_q f(F(hat_xq))*hat_phi_i(hat_xq) det(DF(hat_xq)) hat_wq
 
    On each K, the computation needs two imbricated loops in (q,i).
    The expresion is represented by a tree at compile-time.
    The first subexpr 'f' is represented by field_vf_expr_terminal<f> : it evaluates : 
      value1(q,i) = f(F(hat_xq)) : the values depends only of q and are independent of i.
    They could be evaluated juste before the 'i' loop.
    As the field_vf_expr_terminal<> is general and can handle also more complex expressions
    involving fields with various approximations, all the values on K are evaluated 
    in a specific 'q' loop by
      subexpr1.element_initialize (K);
    The second subexpr 'phi_i' is represdented by the 'test' class (see before).
      value2(q,i) = (hat_phi_i(hat_xq))
    The '*' performs on the fly the product
      value(q,i) = value1(q,i)*value2(q,i)
  
   3) expr(v) = dot(f,grad(v)) dx
    The 'f' function is here vector-valued.
    The expresion is represented by a tree at compile-time :
      dot -> f 
          -> grad -> v
    The 'grad' node returns  
      value(q,i) = trans(inv(DF(hat_wq))*grad_phi_i(hat_xq) that is vector-valued
    The grad_phi values are obtained by a grad_value(q,i) method on the 'test' class.
    The 'dot' performs on the fly the product
      value(q,i) = ot (value1(q,i), value2(q,i))

   This approch generilize for an expression tree.
*/
namespace rheolef {

// ====================================================================
// compute dis_idof on K
// ====================================================================
template <class T, class M>
void
assembly_dis_idof (
  const space_basic<T,M>&              X,
  const geo_basic<T,M>&                dom,
  const geo_element&                   K,
  std::vector<geo_element::size_type>& dis_idx);

// ====================================================================
// common implementation for integration on a band or an usual domain
// ====================================================================
template <class T, class M>
template <class Expr>
void
field_basic<T,M>::assembly_internal (
    const geo_basic<T,M>&         dom,
    const geo_basic<T,M>&         band,
    const band_basic<T,M>&        gh,
    const Expr&                   expr,
    const quadrature_option_type& qopt,
    bool                          is_on_band)
{
  // ------------------------------------
  // 0) initialize the loop
  // ------------------------------------
  const bool ignore_sys_coord = false; // TODO: available as an option via fopt
  const space_basic<T,M>& Xh = expr.get_vf_space();
  const geo_basic<T,M>& omega     = Xh.get_geo();
  const geo_basic<T,M>& bgd_omega = Xh.get_geo().get_background_geo();
  check_macro (band.get_background_geo() == bgd_omega,
	"assembly: incompatible integration domain "<<dom.name() << " and test function based domain "
	<< Xh.get_geo().name());
  resize (Xh, 0);

  quadrature<T> quad;
  if (qopt.get_order() != 0 && qopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max()) {
    quad.set_order (qopt.get_order());
  } else {
    size_type k = Xh.degree();
    if (dom.coordinate_system() != space_constant::cartesian) k++; // multiplies by a 'r' weight
    size_type quad_order = 2*k+1; // see Raviart & Thomas, Masson, 1988, page 130, theorem 5.3.1
    quad.set_order (quad_order);
  }
  quad.set_family (qopt.get_family());
  if (is_on_band) {
    expr.initialize (gh, quad, ignore_sys_coord);
  } else {
    expr.initialize (dom, quad, ignore_sys_coord);
  }
  expr.template valued_check<T>();
  basis_on_pointset<T> piola_on_quad;
  piola_on_quad.set (quad, Xh.get_geo().get_piola_basis());
  vec<T,M>& u = set_u();
  vec<T,M>& b = set_b();
  std::vector<size_type> dis_idx;
  std::vector<T> lk, value;
  tensor_basic<T> DF;
  std::vector<size_type> dis_inod_K;
  size_type d     = dom.dimension();
  size_type map_d = dom.map_dimension();
  bool is_dg = (!Xh.get_numbering().is_continuous() && map_d < omega.map_dimension());
  // ------------------------------------
  // 1) loop on elements
  // ------------------------------------
  for (size_type ie = 0, ne = dom.size(map_d); ie < ne; ie++) {
    const geo_element& K = dom.get_geo_element (map_d, ie);
    // ------------------------------------
    // 1.1) locally compute dofs
    // ------------------------------------
    if (!is_on_band) {
      assembly_dis_idof (Xh, dom, K, dis_idx);
    } else { // on a band
      size_type L_ie = gh.sid_ie2bnd_ie (ie);
      const geo_element& L = band [L_ie];
      Xh.dis_idof (L, dis_idx);
    }
    // ------------------------------------
    // 1.2) locally compute values
    // ------------------------------------
    // locally evaluate int_K expr dx
    // with loop on a quadrature formulae 
    // and accumulate in lk
    value.resize (dis_idx.size());
    lk.resize    (dis_idx.size());
    fill (lk.begin(), lk.end(), T(0));
    expr.element_initialize (K);
    reference_element hat_K = K.variant();
    dom.dis_inod (K, dis_inod_K);
    typename quadrature<T>::const_iterator
      first_quad = quad.begin(hat_K),
      last_quad  = quad.end  (hat_K);
    for (size_type q = 0; first_quad != last_quad; ++first_quad, ++q) {
      jacobian_piola_transformation (dom, piola_on_quad, K, dis_inod_K, q, DF);
      T det_DF = det_jacobian_piola_transformation (DF, d, map_d);
      T wq = 1;
      if (!ignore_sys_coord) {
        wq *= weight_coordinate_system (dom, piola_on_quad, K, dis_inod_K, q);
      }
      wq *= det_DF;
      wq *= (*first_quad).w;
      expr.basis_evaluate (hat_K, q, value);
      for (size_type loc_idof = 0, loc_ndof = value.size(); loc_idof < loc_ndof; ++loc_idof) {
        lk [loc_idof] += value[loc_idof]*wq;
      }
    }
    // -----------------------------------------
    // 1.3) assembly local lk in global field lh
    // -----------------------------------------
    check_macro (dis_idx.size() == lk.size(),
      "incompatible sizes dis_idx.size="<<dis_idx.size()<<", lk.size="<<lk.size());
    for (size_type loc_idof = 0, loc_ndof = lk.size(); loc_idof < loc_ndof; loc_idof++) {
      const T& value = lk [loc_idof];
      size_type dis_idof = dis_idx [loc_idof];
      size_type dis_iub  = Xh.dis_iub(dis_idof);
      if (Xh.dis_is_blocked(dis_idof)) b.dis_entry(dis_iub) += value;
      else                             u.dis_entry(dis_iub) += value;
    }
  }
  // -----------------------------------------
  // 2) finalize the assembly process
  // -----------------------------------------
  u.dis_entry_assembly(set_add_op<T,T>());
  b.dis_entry_assembly(set_add_op<T,T>());
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
    const geo_basic<T,M>&         dom,
    const Expr&                   expr,
    const quadrature_option_type& qopt)
{
  assembly_internal (dom, dom, band_basic<T,M>(), expr, qopt, false);
}
template <class T, class M>
template <class Expr>
inline
void
field_basic<T,M>::assembly (
    const band_basic<T,M>&        gh,
    const Expr&                   expr,
    const quadrature_option_type& qopt)
{
  assembly_internal (gh.level_set(), gh.band(), gh, expr, qopt, true);
}

}// namespace rheolef
#endif // _RHEO_FIELD_VF_ASSEMBLY_H