HArD::Core3D
Hybrid Arbitrary Degree::Core 3D - Library to implement 3D schemes with vertex, edge, face and cell polynomials as unknowns
Loading...
Searching...
No Matches
IntegrateTripleProduct.hpp
Go to the documentation of this file.
1#ifndef INTEGRATE_TRIPLE_PRODUCT_HPP
2#define INTEGRATE_TRIPLE_PRODUCT_HPP
3
4#include <GMpoly_cell.hpp>
5
6namespace HArDCore3D {
7
13typedef boost::multi_array<double, 3> Scalar3Tensor;
14typedef Eigen::Map<Eigen::MatrixXd, Eigen::Unaligned, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> MatrixSlice;
15typedef Eigen::Map<Eigen::VectorXd, Eigen::Unaligned, Eigen::InnerStride<Eigen::Dynamic>> VectorSlice;
16
17//----------------------------------------------------------------------
18//----------------------------------------------------------------------
19// FUNCTIONS TO SLICE SCALAR 3_TENSORS
20//----------------------------------------------------------------------
21//----------------------------------------------------------------------
22
23//
25
28MatrixSlice slice(
29 Scalar3Tensor & tensor,
30 size_t fixed_dim,
31 size_t index
32 )
33{
34 size_t size_dim1 = tensor.shape()[0];
35 size_t size_dim2 = tensor.shape()[1];
36 size_t size_dim3 = tensor.shape()[2];
37
38 size_t start;
39 size_t num_rows;
40 size_t num_cols;
41 size_t column_stride; // no. elements between two column entries
42 size_t row_stride; // no. elements between two row entries
43
44 if (fixed_dim == 0 && index < size_dim1){
45 start = index*size_dim2*size_dim3;
46 num_rows = size_dim2;
47 num_cols = size_dim3;
48 column_stride = 1;
49 row_stride = size_dim3;
50 } else if (fixed_dim == 1 && index < size_dim2){
51 start = index*size_dim3;
52 num_rows = size_dim1;
53 num_cols = size_dim3;
54 column_stride = 1;
55 row_stride = size_dim2*size_dim3;
56 } else if (fixed_dim == 2 && index < size_dim3){
57 // Note: this map can be slow in calculations possibly due to the the spacing of the data (non-contiguous in either index), try using a copy if tests are slow.
58 start = index;
59 num_rows = size_dim1;
60 num_cols = size_dim2;
61 column_stride = size_dim3;
62 row_stride = size_dim2*size_dim3;
63 } else{
64 std::cerr << "[IntegrateTripleProduct] ERROR: Multiarray slicing dimension or index out of range" << std::endl;
65 exit(1);
66 }
67return MatrixSlice(tensor.data() + start, num_rows, num_cols, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>(column_stride, row_stride));
68};
69
71
72VectorSlice slice(
73 Scalar3Tensor & tensor,
74 size_t fixed_dim1,
75 size_t index1,
76 size_t fixed_dim2,
77 size_t index2
78 )
79{
80 size_t size_dim1 = tensor.shape()[0];
81 size_t size_dim2 = tensor.shape()[1];
82 size_t size_dim3 = tensor.shape()[2];
83
84 size_t start;
85 size_t num_elements;
86 size_t stride; // no. elements between two entries
87
88 if (fixed_dim1 == 0 && index1 < size_dim1 && fixed_dim2 == 1 && index2 < size_dim2){
89 start = index1*size_dim2*size_dim3 + index2*size_dim3;
90 num_elements = size_dim3;
91 stride = 1;
92 } else if (fixed_dim1 == 1 && index1 < size_dim2 && fixed_dim2 == 0 && index2 < size_dim1){
93 start = index2*size_dim2*size_dim3 + index1*size_dim3;
94 num_elements = size_dim3;
95 stride = 1;
96 } else if (fixed_dim1 == 0 && index1 < size_dim1 && fixed_dim2 == 2 && index2 < size_dim3){
97 start = index1*size_dim2*size_dim3 + index2;
98 num_elements = size_dim2;
99 stride = size_dim3;
100 } else if (fixed_dim1 == 2 && index1 < size_dim3 && fixed_dim2 == 0 && index2 < size_dim1){
101 start = index1 + index2*size_dim2*size_dim3;
102 num_elements = size_dim2;
103 stride = size_dim3;
104 } else if (fixed_dim1 == 1 && index1 < size_dim2 && fixed_dim2 == 2 && index2 < size_dim3){
105 start = index1*size_dim3 + index2;
106 num_elements = size_dim1;
107 stride = size_dim2*size_dim3;
108 } else if (fixed_dim1 == 2 && index1 < size_dim3 && fixed_dim2 == 1 && index2 < size_dim2) {
109 start = index1 + index2*size_dim3;
110 num_elements = size_dim1;
111 stride = size_dim2*size_dim3;
112 }else{
113 std::cerr << "[IntegrateTripleProduct] ERROR: Multiarray slicing dimension or index out of range" << std::endl;
114 exit(1);
115 }
116 return VectorSlice(tensor.data()+start, num_elements, Eigen::InnerStride<>(stride));
117};
118
119//----------------------------------------------------------------------
120//----------------------------------------------------------------------
121// FUNCTIONS TO COMPUTE TRIPLE INTEGRALS
122//----------------------------------------------------------------------
123//----------------------------------------------------------------------
124
125//
126/* BASIC ONES: Monomial, tensorized, and generic template */
127
129
133template<size_t N>
134Scalar3Tensor tripleInt(
135 const Cell& T,
136 const MonomialScalarBasisCell & basis1,
137 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & basis2,
138 MonomialCellIntegralsType mono_int_map = {}
139 )
140 {
141 size_t dim1 = basis1.dimension();
142 size_t dim2 = basis2.dimension();
143 size_t anc_dim2 = basis2.ancestor().dimension();
144
145 size_t totaldegree = basis1.max_degree()+basis2.ancestor().max_degree()+basis2.ancestor().max_degree();
146 Scalar3Tensor scalar_integrals(boost::extents[dim1][anc_dim2][anc_dim2]);
147
148 // Obtain integration data from IntegrateCellMonomials
149 MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
150
151 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
152 std::fill_n(integrals.data(), integrals.num_elements(), 0);
153
154 for (size_t i = 0; i < dim1; i++) {
155 for (size_t j = 0; j < anc_dim2; j++) {
156 for (size_t k = 0; k <= j; k++) {
157 scalar_integrals[i][j][k] = intmap.at(basis1.powers(i) + basis2.ancestor().powers(j) + basis2.ancestor().powers(k));
158 scalar_integrals[i][k][j] = scalar_integrals[i][j][k];
159 }
160 }
161 }
162
163 for (size_t i = 0; i < N; i++){
164 integrals[boost::indices[boost::multi_array_types::index_range(0,dim1)][boost::multi_array_types::index_range(i*anc_dim2,i*anc_dim2+anc_dim2)][boost::multi_array_types::index_range(i*anc_dim2,i*anc_dim2+anc_dim2)]] = scalar_integrals;
165 }
166
167 return integrals;
168 };
169
171
174template<size_t N>
175Scalar3Tensor tripleInt(
176 const Cell& T,
177 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & tens_family1,
178 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & tens_family2,
179 MonomialCellIntegralsType mono_int_map = {}
180 )
181 {
182 size_t dim1 = tens_family1.dimension();
183 size_t anc_dim1 = tens_family1.ancestor().dimension();
184 size_t dim2 = tens_family2.dimension();
185 size_t anc_dim2 = tens_family2.ancestor().dimension();
186 size_t totaldegree = tens_family1.ancestor().max_degree() + tens_family2.ancestor().max_degree() + tens_family2.ancestor().max_degree();
187
188 // Obtain integration data from IntegrateCellMonomials
189 MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
190
191 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
192 std::fill_n(integrals.data(), integrals.num_elements(), 0);
193
194 Eigen::Matrix3i id = Eigen::Matrix3i::Identity();
195
196 // Loop over the basis vectors: e_i . (e_j x e_k)
197 for (size_t i = 0; i < N; i++){
198 for (size_t j = 0; j < N; j++){
199 for (size_t k = 0; k < j; k++){
200 int sign = id.col(i).dot(id.col(j).cross(id.col(k)));
201 if (sign){
202 // Loop over the monomials
203 for (size_t sub_i = 0; sub_i < anc_dim1; sub_i++){
204 for (size_t sub_j = 0; sub_j < anc_dim2; sub_j++){
205 for (size_t sub_k = 0; sub_k < anc_dim2; sub_k++){
206 integrals[i*anc_dim1 + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] =
207 sign * intmap.at(
208 tens_family1.ancestor().powers(sub_i)
209 + tens_family2.ancestor().powers(sub_j)
210 + tens_family2.ancestor().powers(sub_k)
211 );
212 if (sub_j != sub_k){
213 integrals[i*anc_dim1 + sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[i*anc_dim1 + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
214 }
215 integrals[i*anc_dim1 + sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[i*anc_dim1 + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
216 integrals[i*anc_dim1 + sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[i*anc_dim1 + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
217 }// for sub_k
218 }// for sub_j
219 }// for sub_i
220 };
221 }// for k
222 }// for j
223 }// for i
224 return integrals;
225 };
226
228
231template<size_t N>
232Scalar3Tensor tripleInt(
233 const Cell& T,
234 const GradientBasis<MonomialScalarBasisCell> & grad_basis,
235 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & tens_family,
236 MonomialCellIntegralsType mono_int_map = {}
237 )
238 {
239 size_t dim1 = grad_basis.dimension();
240 size_t dim2 = tens_family.dimension();
241 size_t anc_dim2 = tens_family.ancestor().dimension();
242 size_t totaldegree = grad_basis.ancestor().max_degree()-1 + tens_family.ancestor().max_degree() + tens_family.ancestor().max_degree();
243
244 // Obtain integration data from IntegrateCellMonomials
245 MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
246
247 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
248 std::fill_n(integrals.data(), integrals.num_elements(), 0);
249
250 Eigen::Matrix3i id = Eigen::Matrix3i::Identity();
251
252 // Loop over the basis vectors: e_i . (e_j x e_k)
253 for (size_t i = 0; i < N; i++){
254 for (size_t j = 0; j < N; j++){
255 for (size_t k = 0; k < j; k++){
256 int sign = id.col(i).dot(id.col(j).cross(id.col(k)));
257 if (sign){
258 // Loop over the monomials
259 for (size_t sub_i = 0; sub_i < dim1; sub_i++){
260 if (grad_basis.ancestor().powers(sub_i)[i]>0){
261 for (size_t sub_j = 0; sub_j < anc_dim2; sub_j++){
262 for (size_t sub_k = 0; sub_k < anc_dim2; sub_k++){
263 integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] =
264 sign * grad_basis.ancestor().powers(sub_i)[i]
265 * intmap.at(
266 grad_basis.ancestor().powers(sub_i)-id.col(i)
267 + tens_family.ancestor().powers(sub_j)
268 + tens_family.ancestor().powers(sub_k)
269 ) / T.diam();
270 if (sub_j != sub_k){
271 integrals[sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
272 }
273 integrals[sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
274 integrals[sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
275 }// for sub_k
276 }// for sub_j
277 };
278 }// for sub_i
279 };
280 }// for k
281 }// for j
282 }// for i
283 return integrals;
284 };
285
287
290template<size_t N>
291Scalar3Tensor tripleInt(
292 const Cell& T,
293 const GolyComplBasisCell & golycompl_basis,
294 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & tens_family,
295 MonomialCellIntegralsType mono_int_map = {}
296 )
297 {
298 size_t dim1 = golycompl_basis.dimension();
299 size_t dimPkmo = golycompl_basis.dimPkmo();
300 size_t dimPkmo2 = dim1 - 2*dimPkmo;
301 size_t dim2 = tens_family.dimension();
302
303 size_t anc_dim2 = tens_family.ancestor().dimension();
304 size_t totaldegree = golycompl_basis.max_degree()+ tens_family.ancestor().max_degree() + tens_family.ancestor().max_degree();
305
306 // Obtain integration data from IntegrateCellMonomials
307 MonomialCellIntegralsType intmap = CheckIntegralsDegree(T, totaldegree, mono_int_map);
308
309 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
310 std::fill_n(integrals.data(), integrals.num_elements(), 0);
311
312 Eigen::Matrix3i id = Eigen::Matrix3i::Identity();
313
314 // Loop over the basis vectors: (e_j x e_k)
315 for (size_t j = 0; j < N; j++){
316 for (size_t k = 0; k < N; k++){
317 Eigen::Vector3i direction = id.col(j).cross(id.col(k));
318 int sign = direction.sum();
319 // Loop over the basis/monomials
320 for (size_t sub_j = 0; sub_j < anc_dim2; sub_j++){
321 for (size_t sub_k = 0; sub_k <= sub_j; sub_k++){
322 if (direction==id.col(0)){
323 // Second direction of golycompl
324 for (size_t sub_i = 0; sub_i < dimPkmo; sub_i++){
325 integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = -sign * intmap.at(golycompl_basis.powers(dimPkmo + sub_i) + id.col(2) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
326 integrals[dimPkmo + sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
327 integrals[dimPkmo + sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
328 integrals[dimPkmo + sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
329 }
330 // Third direction of golycompl
331 for (size_t sub_i = 0; sub_i < dimPkmo2; sub_i++){
332 integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = sign * intmap.at(golycompl_basis.powers(2*dimPkmo + sub_i) + id.col(1) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
333 integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
334 integrals[2*dimPkmo + sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
335 integrals[2*dimPkmo + sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
336 }
337 }else if (direction==id.col(1)){
338 // First direction of golycompl
339 for (size_t sub_i = 0; sub_i < dimPkmo; sub_i++){
340 integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = sign * intmap.at(golycompl_basis.powers(sub_i)+ id.col(2) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
341 integrals[sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
342 integrals[sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
343 integrals[sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
344 }
345 // Third direction of golycompl
346 for (size_t sub_i = 0; sub_i < dimPkmo2; sub_i++){
347 integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = -sign * intmap.at(golycompl_basis.powers(2*dimPkmo + sub_i) + id.col(0) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
348 integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
349 integrals[2*dimPkmo + sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
350 integrals[2*dimPkmo + sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[2*dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
351 }
352 }else if (direction==id.col(2)){
353 // First direction of golycompl
354 for (size_t sub_i = 0; sub_i < dimPkmo; sub_i++){
355 integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = -sign * intmap.at(golycompl_basis.powers(sub_i)+ id.col(1) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
356 integrals[sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
357 integrals[sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
358 integrals[sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
359 }
360 // Second direction of golycompl
361 for (size_t sub_i = 0; sub_i < dimPkmo; sub_i++){
362 integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k] = sign * intmap.at(golycompl_basis.powers(dimPkmo + sub_i) + id.col(0) + tens_family.ancestor().powers(sub_j) + tens_family.ancestor().powers(sub_k));
363 integrals[dimPkmo + sub_i][j*anc_dim2 + sub_k][k*anc_dim2 + sub_j] = integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
364 integrals[dimPkmo + sub_i][k*anc_dim2 + sub_j][j*anc_dim2 + sub_k] = -integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
365 integrals[dimPkmo + sub_i][k*anc_dim2 + sub_k][j*anc_dim2 + sub_j] = -integrals[dimPkmo + sub_i][j*anc_dim2 + sub_j][k*anc_dim2 + sub_k];
366 }
367 } //else
368 }// for sub_k
369 }// for sub_j
370 }// for k
371 }// for j
372 return integrals;
373 };
374
375//----------------------------------------------------------------------
376// Family, Shift, Tensorized
377
379template<typename ScalarBasisType1, typename ScalarBasisType2, size_t N>
380Scalar3Tensor tripleInt(
381 const Cell& T,
382 const GradientBasis<ShiftedBasis<ScalarBasisType1>> & grad_shift_basis,
383 const TensorizedVectorFamily<Family<ScalarBasisType2>, N> & basis2,
384 MonomialCellIntegralsType mono_int_map = {}
385 )
386 {
387 size_t dim2 = basis2.dimension();
388 size_t unshifted_dim1 = grad_shift_basis.ancestor().ancestor().dimension();
389 size_t shift = grad_shift_basis.ancestor().shift();
390
391 GradientBasis<ScalarBasisType1> grad_basis(grad_shift_basis.ancestor().ancestor());
392
393 Scalar3Tensor anc_integrals = tripleInt(T, grad_basis, basis2, mono_int_map);
394
395 Scalar3Tensor integrals = anc_integrals[boost::indices[boost::multi_array_types::index_range(shift,unshifted_dim1)][boost::multi_array_types::index_range(0,dim2)][boost::multi_array_types::index_range(0,dim2)]];
396
397 return integrals;
398 };
399
400
402template<typename ScalarBasisType, typename BasisType>
403Scalar3Tensor tripleInt(
404 const Cell& T,
405 const GradientBasis<ShiftedBasis<ScalarBasisType>> & grad_shift_basis,
406 const BasisType & basis2,
407 MonomialCellIntegralsType mono_int_map = {}
408 )
409 {
410 size_t dim2 = basis2.dimension();
411 size_t unshifted_dim1 = grad_shift_basis.ancestor().ancestor().dimension();
412 size_t shift = grad_shift_basis.ancestor().shift();
413
414 GradientBasis<ScalarBasisType> grad_basis(grad_shift_basis.ancestor().ancestor());
415
416 Scalar3Tensor anc_integrals = tripleInt(T, grad_basis, basis2, mono_int_map);
417
418 Scalar3Tensor integrals = anc_integrals[boost::indices[boost::multi_array_types::index_range(shift,unshifted_dim1)][boost::multi_array_types::index_range(0,dim2)][boost::multi_array_types::index_range(0,dim2)]];
419
420 return integrals;
421 };
422
423
425template<typename BasisType, typename ScalarBasisType, size_t N>
426Scalar3Tensor tripleInt(
427 const Cell& T,
428 const BasisType & basis1,
429 const TensorizedVectorFamily<Family<ScalarBasisType>, N> & basis2,
430 MonomialCellIntegralsType mono_int_map = {}
431 )
432 {
433 size_t dim1 = basis1.dimension();
434 size_t dim2 = basis2.dimension();
435 size_t anc_dim2 = basis2.ancestor().dimension();
436 size_t anc_anc_dim2 = basis2.ancestor().ancestor().dimension();
437
438 Eigen::MatrixXd M2 = basis2.ancestor().matrix();
439
440 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
441 std::fill_n(integrals.data(), integrals.num_elements(), 0);
442
443 TensorizedVectorFamily<ScalarBasisType, N> red_basis2(basis2.ancestor().ancestor());
444
445 Scalar3Tensor tensor_integrals = tripleInt(T, basis1, red_basis2, mono_int_map);
446
447 for (size_t j = 0; j < N; j++){
448 for (size_t k = 0; k < N; k++){
449 Scalar3Tensor anc_scalar_integrals = tensor_integrals[boost::indices[boost::multi_array_types::index_range(0,dim1)][boost::multi_array_types::index_range(j*anc_anc_dim2,(j+1)*anc_anc_dim2)][boost::multi_array_types::index_range(k*anc_anc_dim2,(k+1)*anc_anc_dim2)]];
450 for (size_t i = 0; i < dim1; i++){
451 MatrixSlice int_jk = slice(anc_scalar_integrals, 0, i);
452 for (size_t sub_j = 0; sub_j < anc_dim2; sub_j++){
453 for (size_t sub_k = 0; sub_k <= sub_j; sub_k++){
454 integrals[i][anc_dim2*j + sub_j][anc_dim2*k + sub_k] = M2.row(sub_j) * int_jk * M2.row(sub_k).transpose();
455 integrals[i][anc_dim2*j + sub_k][anc_dim2*k + sub_j] = integrals[i][anc_dim2*j + sub_j][anc_dim2*k + sub_k];
456 }
457 }
458 }
459 }
460 }
461 return integrals;
462 };
463
464
466template<typename BasisType, typename ScalarBasisType, size_t N>
467Scalar3Tensor tripleInt(
468 const Cell& T,
469 const Family<BasisType> & basis1,
470 const Family<TensorizedVectorFamily<ScalarBasisType, N>> & basis2,
471 MonomialCellIntegralsType mono_int_map = {}
472 )
473 {
474 size_t dim1 = basis1.dimension();
475 size_t dim2 = basis2.dimension();
476 size_t anc_dim1 = basis1.ancestor().dimension();
477 size_t anc_dim2 = basis2.ancestor().dimension();
478
479 Eigen::MatrixXd M1 = basis1.matrix();
480 Eigen::MatrixXd M2 = basis2.matrix();
481
482 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
483 std::fill_n(integrals.data(), integrals.num_elements(), 0);
484
485 Scalar3Tensor anc_integrals = tripleInt(T, basis1.ancestor(), basis2.ancestor(), mono_int_map);
486
487 for (size_t i = 0; i < dim1; i++){
488 Eigen::MatrixXd int_jk = Eigen::MatrixXd::Zero(anc_dim2, anc_dim2);
489 for (size_t l = 0; l < anc_dim1; l++){
490 int_jk += M1(i, l) * slice(anc_integrals, 0, l);
491 }
492
493 for (size_t j = 0; j < dim2; j++){
494 for (size_t k = 0; k < dim2; k++){
495 integrals[i][j][k] = M2.row(j) * int_jk * M2.row(k).transpose();
496 }
497 }
498 }
499 return integrals;
500 };
501
503template<typename BasisType, size_t N>
504Scalar3Tensor tripleInt(
505 const Cell& T,
506 const Family<BasisType> & basis1,
507 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & basis2,
508 MonomialCellIntegralsType mono_int_map = {}
509 )
510 {
511 size_t dim1 = basis1.dimension();
512 size_t dim2 = basis2.dimension();
513 size_t anc_dim1 = basis1.ancestor().dimension();
514
515 Eigen::MatrixXd M1 = basis1.matrix();
516
517 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
518 std::fill_n(integrals.data(), integrals.num_elements(), 0);
519
520 Scalar3Tensor scalar_integrals = tripleInt(T, basis1.ancestor(), basis2, mono_int_map);
521
522 for (size_t i = 0; i < dim1; i++){
523 MatrixSlice int_jk = slice(integrals, 0, i);
524 for (size_t l = 0; l < anc_dim1; l++){
525 int_jk += M1(i, l) * slice(scalar_integrals, 0, l);
526 }
527 }
528 return integrals;
529 };
530
532template<typename BasisType1, typename BasisType2>
533Scalar3Tensor tripleInt(
534 const Cell& T,
535 const BasisType1 & basis1,
536 const Family<BasisType2> & basis2,
537 MonomialCellIntegralsType mono_int_map = {}
538 )
539 {
540 size_t dim1 = basis1.dimension();
541 size_t dim2 = basis2.dimension();
542
543 Eigen::MatrixXd M2 = basis2.matrix();
544
545 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
546 std::fill_n(integrals.data(), integrals.num_elements(), 0);
547
548 Scalar3Tensor scalar_integrals = tripleInt(T, basis1, basis2.ancestor(), mono_int_map);
549
550 for (size_t i = 0; i < dim1; i++){
551 MatrixSlice int_jk = slice(scalar_integrals, 0, i);
552 for (size_t j = 0; j < dim2; j++){
553 for (size_t k = 0; k <= j; k++){
554 integrals[i][j][k] = M2.row(j) * int_jk * M2.row(k).transpose();
555 integrals[i][k][j] = integrals[i][j][k];
556 }
557 }
558 }
559 return integrals;
560 };
561
563template<size_t N, typename ScalarBasisType>
564Scalar3Tensor tripleInt(
565 const Cell& T,
566 const TensorizedVectorFamily<Family<ScalarBasisType>, N> & tens_family1,
567 const TensorizedVectorFamily<MonomialScalarBasisCell, N> & tens_family2,
568 MonomialCellIntegralsType mono_int_map = {}
569 )
570 {
571 size_t dim1 = tens_family1.dimension();
572 size_t anc_dim1 = tens_family1.ancestor().dimension();
573 size_t anc_anc_dim1 = tens_family1.ancestor().ancestor().dimension();
574 size_t dim2 = tens_family2.dimension();
575 size_t anc_dim2 = tens_family2.ancestor().dimension();
576
577 Eigen::MatrixXd M1 = tens_family1.ancestor().matrix();
578
579 Scalar3Tensor integrals(boost::extents[dim1][dim2][dim2]);
580 std::fill_n(integrals.data(), integrals.num_elements(), 0);
581
582 TensorizedVectorFamily<ScalarBasisType, N> red_basis1(tens_family1.ancestor().ancestor());
583
584 Scalar3Tensor tensor_integrals = tripleInt(T, red_basis1, tens_family2, mono_int_map);
585 Eigen::Matrix3i id = Eigen::Matrix3i::Identity();
586 // Loop over the basis vectors: e_i . (e_j x e_k)
587 for (size_t i = 0; i < N; i++){
588 for (size_t j = 0; j < N; j++){
589 for (size_t k = 0; k < j; k++){
590 int sign = id.col(i).dot(id.col(j).cross(id.col(k)));
591 if (sign){
592 Scalar3Tensor anc_scalar_integrals = tensor_integrals[boost::indices[boost::multi_array_types::index_range(i*anc_anc_dim1,(i+1)*anc_anc_dim1)][boost::multi_array_types::index_range(j*anc_dim2,(j+1)*anc_dim2)][boost::multi_array_types::index_range(k*anc_dim2,(k+1)*anc_dim2)]];
593 for (size_t sub_j = 0; sub_j < anc_dim2; sub_j++){
594 Eigen::MatrixXd fam_int_jk = M1 * slice(anc_scalar_integrals, 1, sub_j);
595 for (size_t sub_i = 0; sub_i < anc_dim1; sub_i++){
596 for (size_t sub_k = 0; sub_k < anc_dim2; sub_k++){
597 integrals[i*anc_dim1+sub_i][j*anc_dim2+sub_j][k*anc_dim2+sub_k] = fam_int_jk(sub_i, sub_k);
598 integrals[i*anc_dim1+sub_i][k*anc_dim2+sub_j][j*anc_dim2+sub_k] = -fam_int_jk(sub_i, sub_k);
599 }
600 }
601 }
602 };
603 }// for k
604 }// for j
605 }// for i
606 return integrals;
607 };
608}
609
610#endif
GMpoly_cell.hpp
HArDCore3D::Family
Family of functions expressed as linear combination of the functions of a given basis.
Definition basis.hpp:389
HArDCore3D::GolyComplBasisCell
Basis for the complement G^{c,k}(T) in P^k(T)^3 of the range of grad.
Definition basis.hpp:1531
HArDCore3D::GradientBasis
Basis for the space of gradients of polynomials.
Definition basis.hpp:1256
HArDCore3D::MonomialScalarBasisCell
Scalar monomial basis on a cell.
Definition basis.hpp:155
HArDCore3D::ShiftedBasis
Generate a basis where the function indices are shifted.
Definition basis.hpp:1040
HArDCore3D::TensorizedVectorFamily
Vector family obtained by tensorization of a scalar family.
Definition basis.hpp:610
HArDCore3D::Matrix
@ Matrix
Definition basis.hpp:67
HArDCore3D::Scalar3Tensor
boost::multi_array< double, 3 > Scalar3Tensor
Definition IntegrateTripleProduct.hpp:13
HArDCore3D::tripleInt
Scalar3Tensor tripleInt(const Cell &T, const MonomialScalarBasisCell &basis1, const TensorizedVectorFamily< MonomialScalarBasisCell, N > &basis2, MonomialCellIntegralsType mono_int_map={})
Computes the triple integral product of a scalar times the dot product of two vectors - basis1(basis2...
Definition IntegrateTripleProduct.hpp:134
HArDCore3D::CheckIntegralsDegree
MonomialCellIntegralsType CheckIntegralsDegree(const Cell &T, const size_t degree, const MonomialCellIntegralsType &mono_int_map={})
Checks if the degree of an existing list of monomial integrals is sufficient, other re-compute and re...
Definition GMpoly_cell.cpp:160
HArDCore3D::VectorSlice
Eigen::Map< Eigen::VectorXd, Eigen::Unaligned, Eigen::InnerStride< Eigen::Dynamic > > VectorSlice
Definition IntegrateTripleProduct.hpp:15
HArDCore3D::slice
MatrixSlice slice(Scalar3Tensor &tensor, size_t fixed_dim, size_t index)
Function to slice a 3-tensor with respect to one index (returns a 2-tensor)
Definition IntegrateTripleProduct.hpp:28
HArDCore3D::MatrixSlice
Eigen::Map< Eigen::MatrixXd, Eigen::Unaligned, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > MatrixSlice
Definition IntegrateTripleProduct.hpp:14
HArDCore3D::MonomialCellIntegralsType
std::unordered_map< VectorZd, double, VecHash > MonomialCellIntegralsType
Type for list of integrals of monomials.
Definition GMpoly_cell.hpp:49
HArDCore3D
Definition ddr-magnetostatics.hpp:41
Generated by doxygen 1.9.8