Intrepid
Intrepid_TensorProductSpaceToolsDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50namespace Intrepid {
51 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
52 class ArrayTypeBasis>
53 void TensorProductSpaceTools::evaluate( ArrayTypeOut &vals ,
54 const ArrayTypeCoeffs &coeffs ,
55 const Array<RCP<ArrayTypeBasis> > &bases )
56 {
57 const unsigned space_dim = bases.size();
58
59#ifdef HAVE_INTREPID_DEBUG
60 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
61 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
62
63 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
64 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
65
66 for (unsigned d=0;d<space_dim;d++)
67 {
68 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
69 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
70 }
71
72 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
73 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
74
75 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
76 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
77
78 int product_of_basis_dimensions = 1;
79 int product_of_basis_points = 1;
80 for (unsigned d=0;d<space_dim;d++)
81 {
82 product_of_basis_dimensions *= bases[d]->dimension(0);
83 product_of_basis_points *= bases[d]->dimension(1);
84 }
85
86 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
87 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
88
89 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
90 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
91#endif
92
93
94 switch (space_dim)
95 {
96 case 2:
97 evaluate2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
98 break;
99 case 3:
100 evaluate3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
101 break;
102 }
103
104 }
105
106 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
107 class ArrayTypeBasis>
109 const ArrayTypeCoeffs &coeffs ,
110 const Array<RCP<ArrayTypeBasis> > &bases )
111 {
112 const unsigned space_dim = bases.size();
113
114#ifdef HAVE_INTREPID_DEBUG
115 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
116 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): vals must be rank 3 array." );
117
118 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
119 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): coeffs must be rank 3 array." );
120
121 for (unsigned d=0;d<space_dim;d++)
122 {
123 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
124 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): each tabulated basis must be rank 2 array." );
125 }
126
127 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
128 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of cells for vals and coeffs must match." );
129
130 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
131 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of fields for vals and coeffs must match." );
132
133 int product_of_basis_dimensions = 1;
134 int product_of_basis_points = 1;
135 for (unsigned d=0;d<space_dim;d++)
136 {
137 product_of_basis_dimensions *= bases[d]->dimension(0);
138 product_of_basis_points *= bases[d]->dimension(1);
139 }
140
141 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
142 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of points in vals and bases ." );
143
144 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
145 ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of basis functions in coeffs and bases ." );
146#endif
147
148
149 switch (space_dim)
150 {
151 case 2:
152 evaluateCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
153 break;
154 case 3:
155 evaluateCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
156 break;
157 }
158
159 }
160
161// template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
162// class ArrayTypeBasis>
163// void TensorProductSpaceTools::evaluateCollocated( ArrayTypeOut &vals ,
164// const ArrayTypeCoeffs &coeffs ,
165// const Array<Array<RCP<ArrayTypeBasis> > > &bases )
166// {
167// const unsigned space_dim = bases.size();
168
169// #ifdef HAVE_INTREPID_DEBUG
170// TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
171// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): vals must be rank 3 array." );
172
173// TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
174// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): coeffs must be rank 3 array." );
175
176// for (unsigned d0=0;d0<bases.size();d0++)
177// {
178// for (unsigned d1=1;d1<space_dim;d1++)
179// {
180// TEUCHOS_TEST_FOR_EXCEPTION( bases[d0][d1]->rank() != 2 , std::invalid_argument ,
181// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): each tabulated basis must be rank 2 array." );
182// }
183// }
184
185// TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
186// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of cells for vals and coeffs must match." );
187
188// TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
189// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Number of fields for vals and coeffs must match." );
190
191// int product_of_basis_dimensions = 1;
192// int product_of_basis_points = 1;
193// for (unsigned d=0;d<space_dim;d++)
194// {
195// product_of_basis_dimensions *= bases[0][d]->dimension(0);
196// product_of_basis_points *= bases[0][d]->dimension(1);
197// }
198
199// TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
200// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of points in vals and bases ." );
201
202// TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
203// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): Incompatible number of basis functions in coeffs and bases ." );
204// #endif
205
206// TEUCHOS_TEST_FOR_EXCEPTION( vals.rank(3) != bases.size() , std::invalid_argument ,
207// ">>> ERROR (TensorProductSpaceTools::evaluateCollocated): wrong dimensions for vals");
208
209
210// switch (space_dim)
211// {
212// case 2:
213// evaluateCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
214// break;
215// case 3:
216// evaluateCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases );
217// break;
218// }
219
220// }
221
222 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
223 class ArrayTypeBasis>
225 const ArrayTypeCoeffs &coeffs ,
226 const Array<RCP<ArrayTypeBasis> > &bases ,
227 const Array<RCP<ArrayTypeBasis> > &Dbases )
228 {
229 const unsigned space_dim = bases.size();
230
231#ifdef HAVE_INTREPID_DEBUG
232 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
233 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
234
235 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
236 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
237
238 TEUCHOS_TEST_FOR_EXCEPTION( Dbases.size() != (int)space_dim , std::invalid_argument ,
239 ">>> ERROR (TensorProductSpaceTools::evaluate): bases and Dbases must be same size." );
240
241 for (unsigned d=0;d<space_dim;d++)
242 {
243 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
244 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
245
246 TEUCHOS_TEST_FOR_EXCEPTION( Dbases[d]->rank() != 3 , std::invalid_argument ,
247 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 3 array." );
248 }
249
250 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
251 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
252
253 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
254 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
255
256 int product_of_basis_dimensions = 1;
257 int product_of_basis_points = 1;
258 for (unsigned d=0;d<space_dim;d++)
259 {
260 product_of_basis_dimensions *= bases[d]->dimension(0);
261 product_of_basis_points *= bases[d]->dimension(1);
262 }
263
264 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
265 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
266
267 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
268 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
269
270 for (unsigned d=0;d<space_dim;d++)
271 {
272 for (unsigned i=0;i<2;i++)
273 {
274 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->dimension(i) != Dbases[d]->dimension(i) , std::invalid_argument ,
275 ">>>ERROR (TensorProductSpaceTools::evaluate): bases and Dbases have incompatible shape." );
276 }
277 }
278#endif
279
280 switch (space_dim)
281 {
282 case 2:
283 TensorProductSpaceTools::evaluateGradient2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
284 break;
285 case 3:
286 TensorProductSpaceTools::evaluateGradient3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
287 break;
288 }
289
290 }
291
292 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
293 class ArrayTypeBasis>
295 const ArrayTypeCoeffs &coeffs ,
296 const Array<RCP<ArrayTypeBasis> > &bases ,
297 const Array<RCP<ArrayTypeBasis> > &Dbases )
298 {
299 const unsigned space_dim = bases.size();
300
301#ifdef HAVE_INTREPID_DEBUG
302 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 4 , std::invalid_argument ,
303 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
304
305 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.rank() != 3 , std::invalid_argument ,
306 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
307
308 TEUCHOS_TEST_FOR_EXCEPTION( Dbases.size() != (int)space_dim , std::invalid_argument ,
309 ">>> ERROR (TensorProductSpaceTools::evaluate): bases and Dbases must be same size." );
310
311 for (unsigned d=0;d<space_dim;d++)
312 {
313 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->rank() != 2 , std::invalid_argument ,
314 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
315
316 TEUCHOS_TEST_FOR_EXCEPTION( Dbases[d]->rank() != 3 , std::invalid_argument ,
317 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 3 array." );
318 }
319
320 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != coeffs.dimension(0) , std::invalid_argument,
321 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
322
323 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != coeffs.dimension(1) , std::invalid_argument,
324 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
325
326 int product_of_basis_dimensions = 1;
327 int product_of_basis_points = 1;
328 for (unsigned d=0;d<space_dim;d++)
329 {
330 product_of_basis_dimensions *= bases[d]->dimension(0);
331 product_of_basis_points *= bases[d]->dimension(1);
332 }
333
334 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_points , std::invalid_argument,
335 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
336
337 TEUCHOS_TEST_FOR_EXCEPTION( coeffs.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
338 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in coeffs and bases ." );
339
340 for (unsigned d=0;d<space_dim;d++)
341 {
342 for (unsigned i=0;i<2;i++)
343 {
344 TEUCHOS_TEST_FOR_EXCEPTION( bases[d]->dimension(i) != Dbases[d]->dimension(i) , std::invalid_argument ,
345 ">>>ERROR (TensorProductSpaceTools::evaluate): bases and Dbases have incompatible shape." );
346 }
347 }
348#endif
349
350 switch (space_dim)
351 {
352 case 2:
353 evaluateGradientCollocated2D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
354 break;
355 case 3:
356 evaluateGradientCollocated3D<Scalar,ArrayTypeOut,ArrayTypeCoeffs,ArrayTypeBasis>( vals , coeffs , bases , Dbases );
357 break;
358 }
359
360 }
361
362 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
363 class ArrayTypeBasis, class ArrayTypeWeights>
364 void TensorProductSpaceTools::moments( ArrayTypeOut &vals ,
365 const ArrayTypeData &data ,
366 const Array<RCP<ArrayTypeBasis> > &basisVals ,
367 const Array<RCP<ArrayTypeWeights> > &wts )
368 {
369 const unsigned space_dim = basisVals.size();
370
371#ifdef HAVE_INTREPID_DEBUG
372 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
373 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
374
375 TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 3 , std::invalid_argument ,
376 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
377
378 TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
379 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
380
381 TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
382 ">>> ERROR (TensorProductSpaceTools::evaluate): quadrature weights ill-sized." );
383
384 for (unsigned d=0;d<space_dim;d++)
385 {
386 TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
387 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
388
389 TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
390 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
391 }
392
393 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
394 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
395
396 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != data.dimension(1) , std::invalid_argument,
397 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
398
399 int product_of_basis_dimensions = 1;
400 int product_of_basis_points = 1;
401 for (unsigned d=0;d<space_dim;d++)
402 {
403 product_of_basis_dimensions *= basisVals[d]->dimension(0);
404 product_of_basis_points *= basisVals[d]->dimension(1);
405 }
406
407 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
408 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
409
410 TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
411 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
412
413#endif
414
415 switch (space_dim)
416 {
417 case 2:
418 moments2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
419 break;
420 case 3:
421 moments3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
422 break;
423 }
424
425 }
426
427 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
428 class ArrayTypeBasis, class ArrayTypeWeights>
430 const ArrayTypeData &data ,
431 const Array<RCP<ArrayTypeBasis> > &basisVals ,
432 const Array<RCP<ArrayTypeWeights> > &wts )
433 {
434 const unsigned space_dim = basisVals.size();
435
436#ifdef HAVE_INTREPID_DEBUG
437 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
438 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
439
440 TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 3 , std::invalid_argument ,
441 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 3 array." );
442
443 TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
444 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
445
446 TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
447 ">>> ERROR (TensorProductSpaceTools::evaluate): quadrature weights ill-sized." );
448
449 for (unsigned d=0;d<space_dim;d++)
450 {
451 TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
452 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
453
454 TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
455 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
456 }
457
458 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
459 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
460
461 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(1) != data.dimension(1) , std::invalid_argument,
462 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of fields for vals and coeffs must match." );
463
464 int product_of_basis_dimensions = 1;
465 int product_of_basis_points = 1;
466 for (unsigned d=0;d<space_dim;d++)
467 {
468 product_of_basis_dimensions *= basisVals[d]->dimension(0);
469 product_of_basis_points *= basisVals[d]->dimension(1);
470 }
471
472 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
473 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
474
475 TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
476 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
477
478#endif
479
480 switch (space_dim)
481 {
482 case 2:
483 moments2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
484 break;
485 case 3:
486 moments3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , wts );
487 break;
488 }
489
490 }
491
492 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
493 class ArrayTypeBasis, class ArrayTypeWeights>
494 void TensorProductSpaceTools::momentsGrad( ArrayTypeOut &vals ,
495 const ArrayTypeData &data ,
496 const Array<RCP<ArrayTypeBasis> > &basisVals ,
497 const Array<RCP<ArrayTypeBasis> > &basisDVals ,
498 const Array<RCP<ArrayTypeWeights> > &wts )
499 {
500 const unsigned space_dim = basisVals.size();
501
502#ifdef HAVE_INTREPID_DEBUG
503 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
504 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
505
506 TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 4 , std::invalid_argument ,
507 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 4 array." );
508
509 TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
510 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
511
512 TEUCHOS_TEST_FOR_EXCEPTION( basisDVals.size() != (int)space_dim , std::invalid_argument ,
513 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
514
515 TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
516 ">>> ERROR (TensorProductSpaceTools::evaluate): quadrature weights ill-sized." );
517
518 for (unsigned d=0;d<space_dim;d++)
519 {
520 TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
521 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
522
523 TEUCHOS_TEST_FOR_EXCEPTION( basisDVals[d]->rank() != 3 , std::invalid_argument ,
524 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated derivative basis must be rank 3 array." );
525
526 TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
527 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
528 }
529
530 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
531 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
532
533 int product_of_basis_dimensions = 1;
534 int product_of_basis_points = 1;
535 for (unsigned d=0;d<space_dim;d++)
536 {
537 product_of_basis_dimensions *= basisVals[d]->dimension(0);
538 product_of_basis_points *= basisVals[d]->dimension(1);
539 }
540
541 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
542 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
543
544 TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
545 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
546
547#endif
548
549 switch (space_dim)
550 {
551 case 2:
552 momentsGrad2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
553 break;
554 case 3:
555 momentsGrad3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
556 break;
557 }
558
559 }
560
561 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
562 class ArrayTypeBasis, class ArrayTypeWeights>
564 const ArrayTypeData &data ,
565 const Array<RCP<ArrayTypeBasis> > &basisVals ,
566 const Array<RCP<ArrayTypeBasis> > &basisDVals ,
567 const Array<RCP<ArrayTypeWeights> > &wts )
568 {
569 const unsigned space_dim = basisVals.size();
570
571#ifdef HAVE_INTREPID_DEBUG
572 TEUCHOS_TEST_FOR_EXCEPTION( vals.rank() != 3 , std::invalid_argument ,
573 ">>> ERROR (TensorProductSpaceTools::evaluate): vals must be rank 3 array." );
574
575 TEUCHOS_TEST_FOR_EXCEPTION( data.rank() != 4 , std::invalid_argument ,
576 ">>> ERROR (TensorProductSpaceTools::evaluate): coeffs must be rank 4 array." );
577
578 TEUCHOS_TEST_FOR_EXCEPTION( basisVals.size() != (int)space_dim , std::invalid_argument ,
579 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
580
581 TEUCHOS_TEST_FOR_EXCEPTION( basisDVals.size() != (int)space_dim , std::invalid_argument ,
582 ">>> ERROR (TensorProductSpaceTools::evaluate): bases ill-sized." );
583
584 TEUCHOS_TEST_FOR_EXCEPTION( wts.size() != (int)space_dim , std::invalid_argument ,
585 ">>> ERROR (TensorProductSpaceTools::evaluate): quadrature weights ill-sized." );
586
587 for (unsigned d=0;d<space_dim;d++)
588 {
589 TEUCHOS_TEST_FOR_EXCEPTION( basisVals[d]->rank() != 2 , std::invalid_argument ,
590 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated basis must be rank 2 array." );
591
592 TEUCHOS_TEST_FOR_EXCEPTION( basisDVals[d]->rank() != 3 , std::invalid_argument ,
593 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated derivative basis must be rank 3 array." );
594
595 TEUCHOS_TEST_FOR_EXCEPTION( wts[d]->rank() != 1 , std::invalid_argument ,
596 ">>> ERROR (TensorProductSpaceTools::evaluate): each tabulated Dbasis must be rank 2 array." );
597 }
598
599 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(0) != data.dimension(0) , std::invalid_argument,
600 ">>> ERROR (TensorProductSpaceTools::evaluate): Number of cells for vals and coeffs must match." );
601
602 int product_of_basis_dimensions = 1;
603 int product_of_basis_points = 1;
604 for (unsigned d=0;d<space_dim;d++)
605 {
606 product_of_basis_dimensions *= basisVals[d]->dimension(0);
607 product_of_basis_points *= basisVals[d]->dimension(1);
608 }
609
610 TEUCHOS_TEST_FOR_EXCEPTION( vals.dimension(2) != product_of_basis_dimensions , std::invalid_argument,
611 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of points in vals and bases ." );
612
613 TEUCHOS_TEST_FOR_EXCEPTION( data.dimension(2) != product_of_basis_points , std::invalid_argument,
614 ">>> ERROR (TensorProductSpaceTools::evaluate): Incompatible number of basis functions in data and bases ." );
615
616#endif
617
618 switch (space_dim)
619 {
620 case 2:
621 momentsGradCollocated2D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
622 break;
623 case 3:
624 momentsGradCollocated3D<Scalar, ArrayTypeOut, ArrayTypeData, ArrayTypeBasis, ArrayTypeWeights>( vals , data , basisVals , basisDVals , wts );
625 break;
626 }
627
628 }
629
630
631 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
632 class ArrayTypeBasis>
633 void TensorProductSpaceTools::evaluate2D( ArrayTypeOut &vals ,
634 const ArrayTypeCoeffs &coeffs ,
635 const Array<RCP<ArrayTypeBasis> > &bases )
636 {
637 const int numBfx = bases[0]->dimension(0);
638 const int numBfy = bases[1]->dimension(0);
639
640 const int numPtsx = bases[0]->dimension(1);
641 const int numPtsy = bases[1]->dimension(1);
642
643 const int numCells = vals.dimension(0);
644 const int numFields = vals.dimension(1);
645
646 const ArrayTypeBasis &Phix = *bases[0];
647 const ArrayTypeBasis &Phiy = *bases[1];
648
649 FieldContainer<double> Xi(numCells,numBfy,numPtsx);
650
651 // sum factorization step
652
653 for (int f=0;f<numFields;f++)
654 {
655 for (int cell=0;cell<numCells;cell++)
656 {
657 for (int j=0;j<numBfy;j++)
658 {
659 for (int i=0;i<numBfx;i++)
660 {
661 const int I = j * numBfx + i;
662 for (int k=0;k<numPtsx;k++)
663 {
664 Xi(cell,j,k) += coeffs(cell,f,I) * Phix(i,k);
665 }
666 }
667 }
668 }
669
670 for (int cell=0;cell<numCells;cell++)
671 {
672 for (int kpty=0;kpty<numPtsy;kpty++)
673 {
674 for (int kptx=0;kptx<numPtsx;kptx++)
675 {
676 vals(cell,f,kptx+numPtsx*kpty) = 0.0;
677 }
678 }
679
680 // evaluation step
681 for (int kpty=0;kpty<numPtsy;kpty++)
682 {
683 for (int kptx=0;kptx<numPtsx;kptx++)
684 {
685 const int I=kptx+numPtsx*kpty;
686
687 for (int j=0;j<numBfy;j++)
688 {
689 vals(cell,f,I) += Phiy(j,kpty) * Xi(cell,j,kptx);
690 }
691 }
692 }
693 }
694 }
695
696 return;
697 }
698
699 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
700 class ArrayTypeBasis>
701 void TensorProductSpaceTools::evaluateCollocated2D( ArrayTypeOut &vals ,
702 const ArrayTypeCoeffs &coeffs ,
703 const Array<RCP<ArrayTypeBasis> > &bases )
704 {
705 // just copy coeffs to vals!
706
707 const int numBfx = bases[0]->dimension(0);
708 const int numBfy = bases[1]->dimension(0);
709
710 const int numCells = vals.dimension(0);
711 const int numFields = vals.dimension(1);
712
713 for (int cell=0;cell<numCells;cell++)
714 {
715 for (int f=0;f<numFields;f++)
716 {
717 for (int j=0;j<numBfy;j++)
718 {
719 for (int i=0;i<numBfx;i++)
720 {
721 const int I = j * numBfx + i;
722 vals(cell,f,I) = coeffs(cell,f,I);
723 }
724 }
725 }
726 }
727
728 return;
729 }
730
731 // template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
732 // class ArrayTypeBasis>
733 // void TensorProductSpaceTools::evaluateCollocated2D( ArrayTypeOut &vals ,
734 // const ArrayTypeCoeffs &coeffs ,
735 // const Array<Array<RCP<ArrayTypeBasis> > > &bases )
736 // {
737 // // just copy coeffs to vals!
738 // const int numCells = vals.dimension(0);
739 // const int numFields = vals.dimension(1);
740
741 // const int numBfx = bases[comp][0]->dimension(0);
742 // const int numBfy = bases[comp][1]->dimension(0);
743
744
745 // for (int cell=0;cell<numCells;cell++)
746 // {
747 // for (int f=0;f<numFields;f++)
748 // {
749 // for (int j=0;j<numBfy;j++)
750 // {
751 // for (int i=0;i<numBfx;i++)
752 // { const int I = j * numBfx + i;
753 // vals(cell,f,I,comp) = coeffs(cell,f,I);
754 // }
755 // }
756 // }
757 // }
758
759 // return;
760 // }
761
762 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
763 class ArrayTypeBasis>
764 void TensorProductSpaceTools::evaluate3D( ArrayTypeOut &vals ,
765 const ArrayTypeCoeffs &coeffs ,
766 const Array<RCP<ArrayTypeBasis> > &bases )
767
768 {
769 // checks to do:
770 // vals point dimension is product of sizes of point arrays
771 // vals
772
773 const int numBfx = bases[0]->dimension(0);
774 const int numBfy = bases[1]->dimension(0);
775 const int numBfz = bases[2]->dimension(0);
776
777 const int numPtsx = bases[0]->dimension(1);
778 const int numPtsy = bases[1]->dimension(1);
779 const int numPtsz = bases[2]->dimension(1);
780
781 const int numCells = vals.dimension(0);
782 const int numFields = vals.dimension(1);
783
784 const ArrayTypeBasis &Phix = *bases[0];
785 const ArrayTypeBasis &Phiy = *bases[1];
786 const FieldContainer<double> &Phiz = *bases[2];
787
788 FieldContainer<double> Xi(numCells,
789 numBfz, numBfy , numPtsx);
790
791 FieldContainer<double> Theta(numCells,
792 numBfz , numPtsy, numPtsx );
793
794
795 for (int f=0;f<numFields;f++)
796 {
797
798 // Xi section
799 for (int c=0;c<numCells;c++)
800 {
801 for (int k=0;k<numBfz;k++)
802 {
803 for (int j=0;j<numBfy;j++)
804 {
805 for (int l=0;l<numPtsx;l++)
806 {
807 for (int i=0;i<numBfx;i++)
808 {
809 const int coeffIndex = k*numBfy*numBfx + j * numBfx + i;
810 Xi(c,k,j,l) += Phix(i,l) * coeffs(c,f,coeffIndex);
811 }
812 }
813 }
814 }
815 }
816
817 // Theta section
818 for (int c=0;c<numCells;c++)
819 {
820 for (int k=0;k<numBfz;k++)
821 {
822 for (int m=0;m<numPtsy;m++)
823 {
824 for (int l=0;l<numPtsx;l++)
825 {
826 for (int j=0;j<numBfy;j++)
827 {
828 Theta(c,k,m,l) += Phiy(j,m) * Xi(c,k,j,l);
829 }
830 }
831 }
832 }
833 }
834
835 // final section
836 for (int c=0;c<numCells;c++)
837 {
838 for (int n=0;n<numPtsz;n++)
839 {
840 for (int m=0;m<numPtsy;m++)
841 {
842 for (int l=0;l<numPtsx;l++)
843 {
844 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) = 0.0;
845 for (int k=0;k<numBfz;k++)
846 {
847 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l) += Theta(c,k,m,l) * Phiz(k,n);
848 }
849 }
850 }
851 }
852 }
853 }
854
855 return;
856
857 }
858
859 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
860 class ArrayTypeBasis>
861 void TensorProductSpaceTools::evaluateCollocated3D( ArrayTypeOut &vals ,
862 const ArrayTypeCoeffs &coeffs ,
863 const Array<RCP<ArrayTypeBasis> > &bases )
864
865 {
866 // copy coeffs to vals
867
868 const int numBfx = bases[0]->dimension(0);
869 const int numBfy = bases[1]->dimension(0);
870 const int numBfz = bases[2]->dimension(0);
871
872 const int numCells = vals.dimension(0);
873 const int numFields = vals.dimension(1);
874
875 for (int cell=0;cell<numCells;cell++)
876 {
877 for (int field=0;field<numFields;field++)
878 {
879 for (int k=0;k<numBfz;k++)
880 {
881 for (int j=0;j<numBfy;j++)
882 {
883 for (int i=0;i<numBfx;i++)
884 {
885 const int I = k*numBfy*numBfx + j * numBfx + i;
886 vals(cell,field,I) = coeffs(cell,field,I);
887 }
888 }
889 }
890 }
891 }
892
893 return;
894
895 }
896
897
898 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
899 class ArrayTypeBasis>
900 void TensorProductSpaceTools::evaluateGradient2D( ArrayTypeOut &vals ,
901 const ArrayTypeCoeffs &coeffs ,
902 const Array<RCP<ArrayTypeBasis> > &bases ,
903 const Array<RCP<ArrayTypeBasis> > &Dbases )
904 {
905
906 const int numBfx = bases[0]->dimension(0);
907 const int numBfy = bases[1]->dimension(0);
908 const int numPtsx = bases[0]->dimension(1);
909 const int numPtsy = bases[1]->dimension(1);
910 const int numCells = vals.dimension(0);
911 const int numFields = vals.dimension(1);
912 const ArrayTypeBasis &Phix = *bases[0];
913 const ArrayTypeBasis &Phiy = *bases[1];
914 const ArrayTypeBasis &DPhix = *Dbases[0];
915 const ArrayTypeBasis &DPhiy = *Dbases[1];
916
917 FieldContainer<double> Xi(numBfy,numPtsx);
918
919 for (int f=0;f<numFields;f++)
920 {
921
922 for (int cell=0;cell<numCells;cell++)
923 {
924 // x derivative
925
926 // sum factorization step
927 for (int j=0;j<numBfy;j++)
928 {
929 for (int k=0;k<numPtsx;k++)
930 {
931 Xi(j,k) = 0.0;
932 }
933 }
934
935 for (int j=0;j<numBfy;j++)
936 {
937 for (int i=0;i<numBfx;i++)
938 {
939 const int I = j * numBfx + i;
940 for (int k=0;k<numPtsx;k++)
941 {
942 Xi(j,k) += coeffs(cell,f,I) * DPhix(i,k,0);
943 }
944 }
945 }
946
947 for (int kpty=0;kpty<numPtsy;kpty++)
948 {
949 for (int kptx=0;kptx<numPtsx;kptx++)
950 {
951 vals(cell,f,kptx+numPtsx*kpty,0) = 0.0;
952 }
953 }
954
955 // evaluation step
956 for (int kpty=0;kpty<numPtsy;kpty++)
957 {
958 for (int kptx=0;kptx<numPtsx;kptx++)
959 {
960 const int I=kptx+numPtsx*kpty;
961
962 for (int j=0;j<numBfy;j++)
963 {
964 vals(cell,f,I,0) += Phiy(j,kpty) * Xi(j,kptx);
965 }
966 }
967 }
968
969 // y derivative
970
971 // sum factorization step
972 for (int j=0;j<numBfy;j++)
973 {
974 for (int k=0;k<numPtsx;k++)
975 {
976 Xi(j,k) = 0.0;
977 }
978 }
979
980 for (int j=0;j<numBfy;j++)
981 {
982 for (int i=0;i<numBfx;i++)
983 {
984 const int I = j * numBfx + i;
985 for (int k=0;k<numPtsx;k++)
986 {
987 Xi(j,k) += coeffs(cell,f,I) * Phix(i,k);
988 }
989 }
990 }
991
992 for (int kpty=0;kpty<numPtsy;kpty++)
993 {
994 for (int kptx=0;kptx<numPtsx;kptx++)
995 {
996 vals(cell,f,kptx+numPtsx*kpty,1) = 0.0;
997 }
998 }
999
1000 // evaluation step
1001 for (int kpty=0;kpty<numPtsy;kpty++)
1002 {
1003 for (int kptx=0;kptx<numPtsx;kptx++)
1004 {
1005 const int I=kptx+numPtsx*kpty;
1006
1007 for (int j=0;j<numBfy;j++)
1008 {
1009 vals(cell,f,I,1) += DPhiy(j,kpty,0) * Xi(j,kptx);
1010 }
1011 }
1012 }
1013 }
1014 }
1015 return;
1016 }
1017
1018 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
1019 class ArrayTypeBasis>
1020 void TensorProductSpaceTools::evaluateGradientCollocated2D( ArrayTypeOut &vals ,
1021 const ArrayTypeCoeffs &coeffs ,
1022 const Array<RCP<ArrayTypeBasis> > &bases ,
1023 const Array<RCP<ArrayTypeBasis> > &Dbases )
1024 {
1025
1026 const int numBfx = bases[0]->dimension(0);
1027 const int numBfy = bases[1]->dimension(0);
1028 const int numPtsx = bases[0]->dimension(1);
1029 const int numPtsy = bases[1]->dimension(1);
1030 const int numCells = vals.dimension(0);
1031 const int numFields = vals.dimension(1);
1032 const ArrayTypeBasis &DPhix = *Dbases[0];
1033 const ArrayTypeBasis &DPhiy = *Dbases[1];
1034
1035 for (int cell=0;cell<numCells;cell++)
1036 {
1037 for (int field=0;field<numFields;field++)
1038 {
1039 for (int j=0;j<numPtsy;j++)
1040 {
1041 for (int i=0;i<numPtsx;i++)
1042 {
1043 const int I = j * numPtsx + i;
1044 vals(cell,field,I,0) = 0.0;
1045 vals(cell,field,I,1) = 0.0;
1046 }
1047 }
1048 }
1049 }
1050
1051 // x derivative
1052 for (int cell=0;cell<numCells;cell++)
1053 {
1054 for (int field=0;field<numFields;field++)
1055 {
1056 for (int j=0;j<numPtsy;j++)
1057 {
1058 for (int i=0;i<numPtsx;i++)
1059 {
1060 const int I = j * numPtsx + i;
1061 for (int ell=0;ell<numBfx;ell++)
1062 {
1063 const int Itmp = j*numPtsx + ell;
1064 vals(cell,field,I,0) += coeffs(cell,field,Itmp) * DPhix( ell , i );
1065 }
1066 }
1067 }
1068 }
1069 }
1070
1071 // y derivative
1072 for (int cell=0;cell<numCells;cell++)
1073 {
1074 for (int field=0;field<numFields;field++)
1075 {
1076 for (int j=0;j<numPtsy;j++)
1077 {
1078 for (int i=0;i<numPtsx;i++)
1079 {
1080 const int I = j * numPtsx + i;
1081 for (int m=0;m<numBfy;m++)
1082 {
1083 const int Itmp = m*numPtsx + i;
1084 vals(cell,field,I,1) += coeffs(cell,field,Itmp) * DPhiy( m , j );
1085 }
1086 }
1087 }
1088 }
1089 }
1090
1091 }
1092
1093 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
1094 class ArrayTypeBasis>
1095 void TensorProductSpaceTools::evaluateGradient3D( ArrayTypeOut &vals ,
1096 const ArrayTypeCoeffs &coeffs ,
1097 const Array<RCP<ArrayTypeBasis> > &bases ,
1098 const Array<RCP<ArrayTypeBasis> > &Dbases )
1099
1100 {
1101 // checks to do:
1102 // vals point dimension is product of sizes of point arrays
1103 // vals
1104
1105 const int numBfx = bases[0]->dimension(0);
1106 const int numBfy = bases[1]->dimension(0);
1107 const int numBfz = bases[2]->dimension(0);
1108 const int numPtsx = bases[0]->dimension(1);
1109 const int numPtsy = bases[1]->dimension(1);
1110 const int numPtsz = bases[2]->dimension(1);
1111 const int numCells = vals.dimension(0);
1112 const int numFields = vals.dimension(1);
1113 const ArrayTypeBasis &Phix = *bases[0];
1114 const ArrayTypeBasis &Phiy = *bases[1];
1115 const ArrayTypeBasis &Phiz = *bases[2];
1116 const ArrayTypeBasis &DPhix = *Dbases[0];
1117 const ArrayTypeBasis &DPhiy = *Dbases[1];
1118 const ArrayTypeBasis &DPhiz = *Dbases[2];
1119
1120 FieldContainer<double> Xi(numCells,
1121 numBfz, numBfy , numPtsx , 3) ;
1122
1123 FieldContainer<double> Theta(numCells,
1124 numBfz , numPtsy, numPtsx , 3);
1125
1126 for (int f=0;f<numFields;f++)
1127 {
1128
1129 // Xi section
1130 for (int c=0;c<numCells;c++)
1131 {
1132 for (int k=0;k<numBfz;k++)
1133 {
1134 for (int j=0;j<numBfy;j++)
1135 {
1136 for (int l=0;l<numPtsx;l++)
1137 {
1138 for (int i=0;i<numBfx;i++)
1139 {
1140 const int coeffIndex = k*numBfz*numBfx + j * numBfx + i;
1141 Xi(c,k,j,l,0) += DPhix(i,l,0) * coeffs(c,f,coeffIndex);
1142 Xi(c,k,j,l,1) += Phix(i,l) * coeffs(c,f,coeffIndex);
1143 Xi(c,k,j,l,2) += Phix(i,l) * coeffs(c,f,coeffIndex);
1144 }
1145 }
1146 }
1147 }
1148 }
1149
1150 // Theta section
1151 for (int c=0;c<numCells;c++)
1152 {
1153 for (int k=0;k<numBfz;k++)
1154 {
1155 for (int m=0;m<numPtsy;m++)
1156 {
1157 for (int l=0;l<numPtsx;l++)
1158 {
1159 for (int j=0;j<numBfy;j++)
1160 {
1161 Theta(c,k,m,l,0) += Phiy(j,m) * Xi(c,k,j,l,0);
1162 Theta(c,k,m,l,1) += DPhiy(j,m,0) * Xi(c,k,j,l,1);
1163 Theta(c,k,m,l,2) += Phiy(j,m) * Xi(c,k,j,l,2);
1164 }
1165 }
1166 }
1167 }
1168 }
1169
1170 // final section
1171 for (int c=0;c<numCells;c++)
1172 {
1173 for (int n=0;n<numPtsz;n++)
1174 {
1175 for (int m=0;m<numPtsy;m++)
1176 {
1177 for (int l=0;l<numPtsx;l++)
1178 {
1179 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) = 0.0;
1180 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) = 0.0;
1181 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) = 0.0;
1182 for (int k=0;k<numBfz;k++)
1183 {
1184 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,0) += Theta(c,k,m,l,0) * Phiz(k,n);
1185 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,1) += Theta(c,k,m,l,1) * Phiz(k,n);
1186 vals(c,f,n*numPtsx*numPtsy+m*numPtsx+l,2) += Theta(c,k,m,l,2) * DPhiz(k,n,0);
1187
1188 }
1189 }
1190 }
1191 }
1192 }
1193 }
1194 return;
1195 }
1196
1197 template<class Scalar, class ArrayTypeOut, class ArrayTypeCoeffs,
1198 class ArrayTypeBasis>
1199 void TensorProductSpaceTools::evaluateGradientCollocated3D( ArrayTypeOut &vals ,
1200 const ArrayTypeCoeffs &coeffs ,
1201 const Array<RCP<ArrayTypeBasis> > &bases ,
1202 const Array<RCP<ArrayTypeBasis> > &Dbases )
1203
1204 {
1205 const int numBfx = bases[0]->dimension(0);
1206 const int numBfy = bases[1]->dimension(0);
1207 const int numBfz = bases[2]->dimension(0);
1208 const int numPtsx = bases[0]->dimension(1);
1209 const int numPtsy = bases[1]->dimension(1);
1210 const int numPtsz = bases[2]->dimension(1);
1211 const int numCells = vals.dimension(0);
1212 const int numFields = vals.dimension(1);
1213 const ArrayTypeBasis &Phix = *bases[0];
1214 const ArrayTypeBasis &Phiy = *bases[1];
1215 const ArrayTypeBasis &Phiz = *bases[2];
1216 const ArrayTypeBasis &DPhix = *Dbases[0];
1217 const ArrayTypeBasis &DPhiy = *Dbases[1];
1218 const ArrayTypeBasis &DPhiz = *Dbases[2];
1219
1220 for (int cell=0;cell<numCells;cell++)
1221 {
1222 for (int field=0;field<numFields;field++)
1223 {
1224 for (int k=0;k<numPtsz;k++)
1225 {
1226 for (int j=0;j<numPtsy;j++)
1227 {
1228 for (int i=0;i<numPtsx;i++)
1229 {
1230 const int I = k * numPtsx * numPtsy + j * numPtsx + i;
1231 vals(cell,field,I,0) = 0.0;
1232 vals(cell,field,I,1) = 0.0;
1233 vals(cell,field,I,2) = 0.0;
1234 }
1235 }
1236 }
1237 }
1238 }
1239
1240 // x derivative
1241 for (int cell=0;cell<numCells;cell++)
1242 {
1243 for (int field=0;field<numFields;field++)
1244 {
1245 for (int k=0;k<numPtsz;k++)
1246 {
1247 for (int j=0;j<numPtsy;j++)
1248 {
1249 for (int i=0;i<numPtsx;i++)
1250 {
1251 const int I = k * numPtsx * numPtsy + j * numPtsx + i;
1252 for (int ell=0;ell<numBfx;ell++)
1253 {
1254 const int Itmp = k * numPtsx * numPtsy + j*numPtsx + ell;
1255 vals(cell,field,I,0) += coeffs(cell,field,Itmp) * DPhix( ell , i );
1256 }
1257 }
1258 }
1259 }
1260 }
1261 }
1262
1263 // y derivative
1264 for (int cell=0;cell<numCells;cell++)
1265 {
1266 for (int field=0;field<numFields;field++)
1267 {
1268 for (int k=0;k<numPtsz;k++)
1269 {
1270 for (int j=0;j<numPtsy;j++)
1271 {
1272 for (int i=0;i<numPtsx;i++)
1273 {
1274 const int I = k * numPtsx * numPtsy + j * numPtsx + i;
1275 for (int m=0;m<numBfy;m++)
1276 {
1277 const int Itmp = k * numPtsx * numPtsy + m * numPtsx + i;
1278 vals(cell,field,I,1) += coeffs(cell,field,Itmp) * DPhiy( m , j );
1279 }
1280 }
1281 }
1282 }
1283 }
1284 }
1285
1286
1287 // z derivative
1288 for (int cell=0;cell<numCells;cell++)
1289 {
1290 for (int field=0;field<numFields;field++)
1291 {
1292 for (int k=0;k<numPtsz;k++)
1293 {
1294 for (int j=0;j<numPtsy;j++)
1295 {
1296 for (int i=0;i<numPtsx;i++)
1297 {
1298 const int I = k * numPtsx * numPtsy + j * numPtsx + i;
1299 for (int n=0;n<numBfz;n++)
1300 {
1301 const int Itmp = n * numPtsx * numPtsy + j * numPtsx + i;
1302 vals(cell,field,I,2) += coeffs(cell,field,Itmp) * DPhiz( n , k );
1303 }
1304 }
1305 }
1306 }
1307 }
1308 }
1309
1310
1311
1312 return;
1313 }
1314
1315 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1316 class ArrayTypeBasis, class ArrayTypeWeights>
1317 void TensorProductSpaceTools::moments2D( ArrayTypeOut &vals ,
1318 const ArrayTypeData &data ,
1319 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1320 const Array<RCP<ArrayTypeWeights> > &wts )
1321 {
1322 const int numBfx = basisVals[0]->dimension(0);
1323 const int numBfy = basisVals[1]->dimension(0);
1324 const int numPtsx = basisVals[0]->dimension(1);
1325 const int numPtsy = basisVals[1]->dimension(1);
1326 const int numCells = vals.dimension(0);
1327 const int numFields = vals.dimension(1);
1328 const ArrayTypeBasis &Phix = *basisVals[0];
1329 const ArrayTypeBasis &Phiy = *basisVals[1];
1330
1331 FieldContainer<double> Xi(numBfx,numPtsy);
1332
1333 const ArrayTypeWeights &wtsx = *wts[0];
1334 const ArrayTypeWeights &wtsy = *wts[1];
1335
1336 for (int f=0;f<numFields;f++)
1337 {
1338 for (int cell=0;cell<numCells;cell++)
1339 {
1340 // sum factorization step
1341 for (int i=0;i<numBfx;i++)
1342 {
1343 for (int k=0;k<numPtsy;k++)
1344 {
1345 Xi(i,k) = 0.0;
1346 }
1347 }
1348
1349 for (int i=0;i<numBfx;i++)
1350 {
1351 for (int l=0;l<numPtsy;l++)
1352 {
1353 for (int k=0;k<numPtsx;k++)
1354 {
1355 Xi(i,l) += wtsx(k) * data(cell,f,l*numPtsx+k) * Phix(i,k);
1356 }
1357 }
1358 }
1359
1360 for (int j=0;j<numBfy;j++)
1361 {
1362 for (int i=0;i<numBfx;i++)
1363 {
1364 vals(cell,f,j*numBfx+i) = 0.0;
1365 }
1366 }
1367
1368 // evaluate moments with sum factorization
1369 for (int j=0;j<numBfy;j++)
1370 {
1371 for (int i=0;i<numBfx;i++)
1372 {
1373 for (int l=0;l<numPtsy;l++)
1374 {
1375 vals(cell,f,j*numBfx+i) += wtsy(l) * Phiy(j,l) * Xi(i,l);
1376 }
1377 }
1378 }
1379 }
1380 }
1381 return;
1382
1383 }
1384
1385 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1386 class ArrayTypeBasis, class ArrayTypeWeights>
1387 void TensorProductSpaceTools::momentsCollocated2D( ArrayTypeOut &vals ,
1388 const ArrayTypeData &data ,
1389 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1390 const Array<RCP<ArrayTypeWeights> > &wts )
1391 {
1392 const int numBfx = basisVals[0]->dimension(0);
1393 const int numBfy = basisVals[1]->dimension(0);
1394 const int numPtsy = basisVals[1]->dimension(1);
1395 const int numCells = vals.dimension(0);
1396 const int numFields = vals.dimension(1);
1397
1398 FieldContainer<double> Xi(numBfx,numPtsy);
1399
1400 const ArrayTypeWeights &wtsx = *wts[0];
1401 const ArrayTypeWeights &wtsy = *wts[1];
1402
1403 for (int cell=0;cell<numCells;cell++)
1404 {
1405 for (int field=0;field<numFields;field++)
1406 {
1407 for (int i=0;i<numBfx;i++)
1408 {
1409 for (int j=0;j<numBfy;j++)
1410 {
1411 const int I = numBfy * i + j;
1412 vals(cell,field,I) = wtsx(i) * wtsy(j) * data(cell,field,I);
1413 }
1414 }
1415 }
1416 }
1417
1418 return;
1419
1420 }
1421
1422 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1423 class ArrayTypeBasis, class ArrayTypeWeights>
1424 void TensorProductSpaceTools::moments3D( ArrayTypeOut &vals ,
1425 const ArrayTypeData &data ,
1426 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1427 const Array<RCP<ArrayTypeWeights> > &wts )
1428 {
1429 const int numBfx = basisVals[0]->dimension(0);
1430 const int numBfy = basisVals[1]->dimension(0);
1431 const int numBfz = basisVals[2]->dimension(0);
1432
1433 const int numPtsx = basisVals[0]->dimension(1);
1434 const int numPtsy = basisVals[1]->dimension(1);
1435 const int numPtsz = basisVals[2]->dimension(1);
1436
1437 const int numCells = vals.dimension(0);
1438 const int numFields = vals.dimension(1);
1439
1440 const ArrayTypeBasis &Phix = *basisVals[0];
1441 const ArrayTypeBasis &Phiy = *basisVals[1];
1442 const ArrayTypeBasis &Phiz = *basisVals[2];
1443
1444 const ArrayTypeWeights &Wtx = *wts[0];
1445 const ArrayTypeWeights &Wty = *wts[1];
1446 const ArrayTypeWeights &Wtz = *wts[2];
1447
1448 FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy);
1449 FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz);
1450
1451 for (int f=0;f<numFields;f++)
1452 {
1453 // Xi phase
1454 for (int c=0;c<numCells;c++)
1455 {
1456 for (int i=0;i<numBfx;i++)
1457 {
1458 for (int n=0;n<numPtsz;n++)
1459 {
1460 for (int m=0;m<numPtsy;m++)
1461 {
1462 for (int l=0;l<numPtsx;l++)
1463 {
1464 Xi(c,i,n,m) += Wtx(l) * Phix(i,l) * data(c,f,n*numPtsy*numPtsx+m*numPtsx+l);
1465 }
1466 }
1467 }
1468 }
1469 }
1470
1471 // Theta phase
1472 for (int c=0;c<numCells;c++)
1473 {
1474 for (int j=0;j<numBfy;j++)
1475 {
1476 for (int i=0;i<numBfx;i++)
1477 {
1478 for (int n=0;n<numPtsz;n++)
1479 {
1480 for (int m=0;m<numPtsy;m++)
1481 {
1482 Theta(c,j,i,n) += Wty(m) * Phiy(j,m) * Xi(c,i,n,m);
1483 }
1484 }
1485 }
1486 }
1487 }
1488
1489 // Final phase
1490 for (int c=0;c<numCells;c++)
1491 {
1492 for (int k=0;k<numBfz;k++)
1493 {
1494 for (int j=0;j<numBfy;j++)
1495 {
1496 for (int i=0;i<numBfx;i++)
1497 {
1498 const int momIdx = k*numBfx*numBfy+j*numBfx+i;
1499 for (int n=0;n<numPtsz;n++)
1500 {
1501 vals(c,f,momIdx) += Wtz(n) * Phiz(k,n) * Theta(c,j,i,n);
1502 }
1503 }
1504 }
1505 }
1506 }
1507
1508 }
1509 return;
1510
1511 }
1512
1513 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1514 class ArrayTypeBasis, class ArrayTypeWeights>
1515 void TensorProductSpaceTools::momentsCollocated3D( ArrayTypeOut &vals ,
1516 const ArrayTypeData &data ,
1517 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1518 const Array<RCP<ArrayTypeWeights> > &wts )
1519 {
1520 const int numBfx = basisVals[0]->dimension(0);
1521 const int numBfy = basisVals[1]->dimension(0);
1522 const int numBfz = basisVals[2]->dimension(0);
1523
1524 const int numCells = vals.dimension(0);
1525 const int numFields = vals.dimension(1);
1526
1527 const ArrayTypeWeights &Wtx = *wts[0];
1528 const ArrayTypeWeights &Wty = *wts[1];
1529 const ArrayTypeWeights &Wtz = *wts[2];
1530
1531 for (int cell=0;cell<numCells;cell++)
1532 {
1533 for (int field=0;field<numFields;field++)
1534 {
1535 for (int k=0;k<numBfz;k++)
1536 {
1537 for (int j=0;j<numBfy;j++)
1538 {
1539 for (int i=0;i<numBfx;i++)
1540 {
1541 const int I = k * numBfy * numBfx + j * numBfx + i;
1542 vals(cell,field,I) = Wtx(i) * Wty(j) * Wtz(k) * data(cell,field,I);
1543 }
1544 }
1545 }
1546 }
1547 }
1548
1549 return;
1550
1551 }
1552
1553 // data is now (C,P,D)
1554 // want to compute the moments against the gradients of the basis
1555 // functions.
1556
1557 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1558 class ArrayTypeBasis, class ArrayTypeWeights>
1559 void TensorProductSpaceTools::momentsGrad2D( ArrayTypeOut &vals ,
1560 const ArrayTypeData &data ,
1561 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1562 const Array<RCP<ArrayTypeBasis> > &Dbases ,
1563 const Array<RCP<ArrayTypeWeights> > &wts )
1564 {
1565
1566 const int numBfx = basisVals[0]->dimension(0);
1567 const int numBfy = basisVals[1]->dimension(0);
1568 const int numPtsx = basisVals[0]->dimension(1);
1569 const int numPtsy = basisVals[1]->dimension(1);
1570 const int numCells = vals.dimension(0);
1571 const int numFields = vals.dimension(1);
1572 const ArrayTypeBasis &Phix = *basisVals[0];
1573 const ArrayTypeBasis &Phiy = *basisVals[1];
1574 const ArrayTypeBasis &DPhix = *Dbases[0];
1575 const ArrayTypeBasis &DPhiy = *Dbases[1];
1576 const ArrayTypeWeights &wtsx = *wts[0];
1577 const ArrayTypeWeights &wtsy = *wts[1];
1578
1579 FieldContainer<double> Xi(numBfx,numPtsy,2);
1580
1581 for (int f=0;f<numFields;f++)
1582 {
1583 // Xi phase
1584 for (int c=0;c<numCells;c++)
1585 {
1586 for (int i=0;i<numBfx;i++)
1587 {
1588 for (int m=0;m<numPtsy;m++)
1589 {
1590 for (int l=0;l<numPtsx;l++)
1591 {
1592 Xi(i,m,0) += wtsx(l) * data(c,f,m*numPtsy+l,0) * DPhix(i,l);
1593 Xi(i,m,1) += wtsx(l) * data(c,f,m*numPtsy+l,1) * Phix(i,l);
1594 }
1595 }
1596 }
1597 }
1598
1599 // main phase
1600 for (int c=0;c<numCells;c++)
1601 {
1602 for (int j=0;j<numBfy;j++)
1603 {
1604 for (int i=0;i<numBfx;i++)
1605 {
1606 const int bfIdx = j*numBfx+i;
1607 vals(c,f,bfIdx) = 0.0;
1608 for (int m=0;m<numPtsy;m++)
1609 {
1610 vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,0) * Phiy(j,m);
1611 vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,1) * DPhiy(j,m);
1612 }
1613 }
1614 }
1615 }
1616 }
1617 return;
1618
1619 }
1620
1621 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1622 class ArrayTypeBasis, class ArrayTypeWeights>
1623 void TensorProductSpaceTools::momentsGradCollocated2D( ArrayTypeOut &vals ,
1624 const ArrayTypeData &data ,
1625 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1626 const Array<RCP<ArrayTypeBasis> > &Dbases ,
1627 const Array<RCP<ArrayTypeWeights> > &wts )
1628 {
1629
1630 const int numBfx = basisVals[0]->dimension(0);
1631 const int numBfy = basisVals[1]->dimension(0);
1632 const int numPtsx = basisVals[0]->dimension(1);
1633 const int numPtsy = basisVals[1]->dimension(1);
1634 const int numCells = vals.dimension(0);
1635 const int numFields = vals.dimension(1);
1636 const ArrayTypeBasis &Phix = *basisVals[0];
1637 const ArrayTypeBasis &Phiy = *basisVals[1];
1638 const ArrayTypeBasis &DPhix = *Dbases[0];
1639 const ArrayTypeBasis &DPhiy = *Dbases[1];
1640 const ArrayTypeWeights &wtsx = *wts[0];
1641 const ArrayTypeWeights &wtsy = *wts[1];
1642
1643 for (int cell=0;cell<numCells;cell++)
1644 {
1645 for (int field=0;field<numFields;field++)
1646 {
1647 for (int j=0;j<numBfy;j++)
1648 {
1649 for (int i=0;i<numBfx;i++)
1650 {
1651 const int I=j*numBfx+i;
1652 vals(cell,field,I) = 0.0;
1653 }
1654 }
1655 }
1656 }
1657
1658 for (int cell=0;cell<numCells;cell++)
1659 {
1660 for (int field=0;field<numFields;field++)
1661 {
1662 for (int j=0;j<numBfy;j++)
1663 {
1664 for (int i=0;i<numBfx;i++)
1665 {
1666 for (int m=0;m<numPtsx;m++)
1667 {
1668 const int Itmp = j * numBfy + m;
1669 vals(cell,field,Itmp) += wtsx(m) * wtsy(j) * DPhix(i,m);
1670 }
1671 }
1672 }
1673 }
1674 }
1675
1676 for (int cell=0;cell<numCells;cell++)
1677 {
1678 for (int field=0;field<numFields;field++)
1679 {
1680 for (int j=0;j<numBfy;j++)
1681 {
1682 for (int i=0;i<numBfx;i++)
1683 {
1684 for (int n=0;n<numPtsy;n++)
1685 {
1686 const int Itmp = n * numBfy + i;
1687 vals(cell,field,Itmp) += wtsx(i) * wtsy(n) * DPhiy(j,n);
1688 }
1689 }
1690 }
1691 }
1692 }
1693
1694 }
1695
1696 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1697 class ArrayTypeBasis, class ArrayTypeWeights>
1698 void TensorProductSpaceTools::momentsGrad3D( ArrayTypeOut &vals ,
1699 const ArrayTypeData &data ,
1700 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1701 const Array<RCP<ArrayTypeBasis> > &basisDVals ,
1702 const Array<RCP<ArrayTypeWeights> > &wts )
1703 {
1704
1705 const int numBfx = basisVals[0]->dimension(0);
1706 const int numBfy = basisVals[1]->dimension(0);
1707 const int numBfz = basisVals[2]->dimension(0);
1708 const int numPtsx = basisVals[0]->dimension(1);
1709 const int numPtsy = basisVals[1]->dimension(1);
1710 const int numPtsz = basisVals[2]->dimension(1);
1711 const int numCells = vals.dimension(0);
1712 const int numFields = vals.dimension(1);
1713 const ArrayTypeBasis &Phix = *basisVals[0];
1714 const ArrayTypeBasis &Phiy = *basisVals[1];
1715 const ArrayTypeBasis &Phiz = *basisVals[2];
1716 const ArrayTypeBasis &DPhix = *basisDVals[0];
1717 const ArrayTypeBasis &DPhiy = *basisDVals[1];
1718 const ArrayTypeBasis &DPhiz = *basisDVals[2];
1719 const ArrayTypeWeights &wtsx = *wts[0];
1720 const ArrayTypeWeights &wtsy = *wts[1];
1721 const ArrayTypeWeights &wtsz = *wts[2];
1722
1723 FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy,3);
1724 FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz,3);
1725
1726 // Xi phase
1727 for (int f=0;f<numFields;f++)
1728 {
1729 for (int c=0;c<numCells;c++)
1730 {
1731 for (int i=0;i<numBfx;i++)
1732 {
1733 for (int n=0;n<numPtsz;n++)
1734 {
1735 for (int m=0;m<numPtsy;m++)
1736 {
1737 for (int l=0;l<numPtsx;l++)
1738 {
1739 const int dataIdx = n * numPtsy * numPtsx + m * numPtsx + l;
1740 Xi(c,i,n,m,0) += wtsx(l) * DPhix(i,l,0) * data(c,f,dataIdx,0);
1741 Xi(c,i,n,m,1) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,1);
1742 Xi(c,i,n,m,2) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,2);
1743 }
1744 }
1745 }
1746 }
1747 }
1748
1749 // Theta phase
1750 for (int c=0;c<numCells;c++)
1751 {
1752 for (int j=0;j<numBfy;j++)
1753 {
1754 for (int i=0;i<numBfx;i++)
1755 {
1756 for (int n=0;n<numPtsz;n++)
1757 {
1758 for (int m=0;m<numPtsy;m++)
1759 {
1760 Theta(c,j,i,n,0) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,0);
1761 Theta(c,j,i,n,1) += wtsy(j) * DPhiy(j,m,0) * Xi(c,i,n,m,1);
1762 Theta(c,j,i,n,2) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,2);
1763 }
1764 }
1765 }
1766 }
1767 }
1768
1769 // last phase
1770 for (int c=0;c<numCells;c++)
1771 {
1772 for (int k=0;k<numBfz;k++)
1773 {
1774 for (int j=0;j<numBfy;j++)
1775 {
1776 for (int i=0;i<numBfx;i++)
1777 {
1778 const int momIdx = k * numBfx * numBfy + j * numBfx + i;
1779 for (int n=0;n<numPtsz;n++)
1780 {
1781 vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,0) * Phiz(k,n);
1782 vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,1) * Phiz(k,n);
1783 vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,2) * DPhiz(k,n,0);
1784 }
1785 }
1786 }
1787 }
1788 }
1789 }
1790
1791}
1792
1793 template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1794 class ArrayTypeBasis, class ArrayTypeWeights>
1795 void TensorProductSpaceTools::momentsGradCollocated3D( ArrayTypeOut &vals ,
1796 const ArrayTypeData &data ,
1797 const Array<RCP<ArrayTypeBasis> > &basisVals ,
1798 const Array<RCP<ArrayTypeBasis> > &basisDVals ,
1799 const Array<RCP<ArrayTypeWeights> > &wts )
1800 {
1801
1802 const int numBfx = basisVals[0]->dimension(0);
1803 const int numBfy = basisVals[1]->dimension(0);
1804 const int numBfz = basisVals[2]->dimension(0);
1805 const int numPtsx = basisVals[0]->dimension(1);
1806 const int numPtsy = basisVals[1]->dimension(1);
1807 const int numPtsz = basisVals[2]->dimension(1);
1808 const int numCells = vals.dimension(0);
1809 const int numFields = vals.dimension(1);
1810 const ArrayTypeBasis &Phix = *basisVals[0];
1811 const ArrayTypeBasis &Phiy = *basisVals[1];
1812 const ArrayTypeBasis &Phiz = *basisVals[2];
1813 const ArrayTypeBasis &DPhix = *basisDVals[0];
1814 const ArrayTypeBasis &DPhiy = *basisDVals[1];
1815 const ArrayTypeBasis &DPhiz = *basisDVals[2];
1816 const ArrayTypeWeights &wtsx = *wts[0];
1817 const ArrayTypeWeights &wtsy = *wts[1];
1818 const ArrayTypeWeights &wtsz = *wts[2];
1819
1820 for (int cell=0;cell<numCells;cell++)
1821 {
1822 for (int field=0;field<numFields;field++)
1823 {
1824 // x component of data versus x derivative of bases
1825 for (int k=0;k<numBfz;k++)
1826 {
1827 for (int j=0;j<numBfy;j++)
1828 {
1829 for (int i=0;i<numBfx;i++)
1830 {
1831 const int I = numBfy * numBfx * k + numBfy * j + i;
1832 for (int ell=0;ell<numPtsx;ell++)
1833 vals(cell,field,I) += wtsx(ell) * wtsy(j) * wtsz(k) * DPhix( i , ell );
1834 }
1835 }
1836 }
1837 // y component of data versus y derivative of bases
1838 for (int k=0;k<numBfz;k++)
1839 {
1840 for (int j=0;j<numBfy;j++)
1841 {
1842 for (int i=0;i<numBfx;i++)
1843 {
1844 const int I = numBfy * numBfx * k + numBfy * j + i;
1845 for (int m=0;m<numPtsy;m++)
1846 vals(cell,field,I) += wtsx(i) * wtsy(m) * wtsz(k) * DPhiy( j , m );
1847 }
1848 }
1849 }
1850 // z component of data versus z derivative of bases
1851 for (int k=0;k<numBfz;k++)
1852 {
1853 for (int j=0;j<numBfy;j++)
1854 {
1855 for (int i=0;i<numBfx;i++)
1856 {
1857 const int I = numBfy * numBfx * k + numBfy * j + i;
1858 for (int n=0;n<numPtsz;n++)
1859 vals(cell,field,I) += wtsx(i) * wtsy(j) * wtsz(n) * DPhiz( k , n );
1860 }
1861 }
1862 }
1863 }
1864 }
1865
1866}
1867
1868
1869
1870} // end namespace Intrepid
static void evaluateCollocated(ArrayTypeOut &vals, const ArrayTypeCoeffs &coeffs, const Array< RCP< ArrayTypeBasis > > &bases)
Computes point values of a set of polynomials expressed in a tensor product basis at output points....
static void momentsGrad(ArrayTypeOut &vals, const ArrayTypeData &data, const Array< RCP< ArrayTypeBasis > > &basisVals, const Array< RCP< ArrayTypeBasis > > &basisDVals, const Array< RCP< ArrayTypeWeights > > &wts)
Computes the moments of a collection of F1 data integrated against a list of functions tabulated at p...
static void evaluateGradient(ArrayTypeOut &vals, const ArrayTypeCoeffs &coeffs, const Array< RCP< ArrayTypeBasis > > &bases, const Array< RCP< ArrayTypeBasis > > &Dbases)
Given a polynomial expressed in a tensor product basis, evaluates the gradient at a tensor product of...
static void momentsGradCollocated(ArrayTypeOut &vals, const ArrayTypeData &data, const Array< RCP< ArrayTypeBasis > > &basisVals, const Array< RCP< ArrayTypeBasis > > &basisDVals, const Array< RCP< ArrayTypeWeights > > &wts)
Computes the moments of a collection of F1 data integrated against a list of functions tabulated at p...
static void moments(ArrayTypeOut &vals, const ArrayTypeData &data, const Array< RCP< ArrayTypeBasis > > &basisVals, const Array< RCP< ArrayTypeWeights > > &wts)
Computes the moments of a set of data integrated against a basis tabulated at points.
static void momentsCollocated(ArrayTypeOut &vals, const ArrayTypeData &data, const Array< RCP< ArrayTypeBasis > > &basisVals, const Array< RCP< ArrayTypeWeights > > &wts)
Computes the moments of a set of data integrated against a basis tabulated at points,...
static void evaluate(ArrayTypeOut &vals, const ArrayTypeCoeffs &coeffs, const Array< RCP< ArrayTypeBasis > > &bases)
Computes point values of a set of polynomials expressed in a tensor product basis at output points....
static void evaluateGradientCollocated(ArrayTypeOut &vals, const ArrayTypeCoeffs &coeffs, const Array< RCP< ArrayTypeBasis > > &bases, const Array< RCP< ArrayTypeBasis > > &Dbases)
Given a polynomial expressed in a tensor product basis, evaluates the gradient at a tensor product of...