/usr/include/rheolef/field_expr_v2_value_assembly.h is in librheolef-dev 6.7-6.
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 | #ifndef _RHEO_FIELD_EXPR_V2_VALUE_ASSEMBLY_H
#define _RHEO_FIELD_EXPR_V2_VALUE_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:
I = int_omega expr(x) dx
The integrals are evaluated over each element K of the domain omega
by using a quadrature formulae given by qopt
expr(x) is a nonlinear expression
The integrals over K are transformed on the reference element with
the piola transformation:
F : hat_K ---> K
hat_x |--> x = F(hat_x)
int_K expr(x) dx
= int_{hat_K} expr(F(hat_x) det(DF(hat_x)) d hat_x
= sum_q expr(F(hat_xq)) det(DF(hat_xq)) hat_wq
On each K, the computation needs a loop in q.
The expresion is represented by a tree at compile-time.
The 'expr' is represented by field_expr_terminal<f> : it evaluation gives :
value1(q) = expr(F(hat_xq))
This approch generilizes to an expression tree.
*/
namespace rheolef { namespace details {
template <class T, class M, class Expr, class Result>
void
field_expr_v2_value_assembly (
const geo_basic<T,M>& omega,
const Expr& expr,
const quadrature_option_type& qopt,
Result& result)
{
typedef typename geo_basic<T,M>::size_type size_type;
// ------------------------------------
// 0) initialize the loop
// ------------------------------------
quadrature<T> quad;
check_macro (qopt.get_order() != std::numeric_limits<quadrature_option_type::size_type>::max(),
"integrate: the quadrature formulae order may be chosen");
quad.set_order (qopt.get_order());
quad.set_family (qopt.get_family());
expr.initialize (omega, quad);
expr.template valued_check<Result>();
basis_on_pointset<T> piola_on_quad;
piola_on_quad.set (quad, omega.get_piola_basis());
std::vector<size_type> dis_idx;
std::vector<Result> value;
tensor_basic<T> DF;
std::vector<size_type> dis_inod_K;
size_type d = omega.dimension();
size_type map_d = omega.map_dimension();
// ------------------------------------
// 1) loop on elements
// ------------------------------------
result = Result(0);
for (size_type ie = 0, ne = omega.size(map_d); ie < ne; ie++) {
const geo_element& K = omega.get_geo_element (map_d, ie);
// ------------------------------------
// 1.1) locally compute values
// ------------------------------------
// locally evaluate int_K expr dx
// with loop on a quadrature formulae
// and accumulate in result
value.resize (dis_idx.size());
reference_element hat_K = K.variant();
omega.dis_inod (K, dis_inod_K);
expr.evaluate (K, value);
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 (omega, piola_on_quad, K, dis_inod_K, q, DF);
T det_DF = det_jacobian_piola_transformation (DF, d, map_d);
T wq = weight_coordinate_system (omega, piola_on_quad, K, dis_inod_K, q);
wq *= det_DF;
wq *= (*first_quad).w;
result += value[q]*wq;
}
}
#ifdef _RHEOLEF_HAVE_MPI
if (is_distributed<M>::value) {
result = mpi::all_reduce (omega.geo_element_ownership(0).comm(), result, std::plus<Result>());
}
#endif // _RHEOLEF_HAVE_MPI
}
}}// namespace rheolef::details
#endif // _RHEO_FIELD_EXPR_V2_VALUE_ASSEMBLY_H
|