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 
50 namespace 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>
108  void TensorProductSpaceTools::evaluateCollocated( ArrayTypeOut &vals ,
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>
224  void TensorProductSpaceTools::evaluateGradient( ArrayTypeOut &vals ,
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>
294  void TensorProductSpaceTools::evaluateGradientCollocated( ArrayTypeOut &vals ,
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>
429  void TensorProductSpaceTools::momentsCollocated( ArrayTypeOut &vals ,
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>
563  void TensorProductSpaceTools::momentsGradCollocated( ArrayTypeOut &vals ,
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 numPtsx = basisVals[0]->dimension(1);
1395  const int numPtsy = basisVals[1]->dimension(1);
1396  const int numCells = vals.dimension(0);
1397  const int numFields = vals.dimension(1);
1398 
1399  FieldContainer<double> Xi(numBfx,numPtsy);
1400 
1401  const ArrayTypeWeights &wtsx = *wts[0];
1402  const ArrayTypeWeights &wtsy = *wts[1];
1403 
1404  for (int cell=0;cell<numCells;cell++)
1405  {
1406  for (int field=0;field<numFields;field++)
1407  {
1408  for (int i=0;i<numBfx;i++)
1409  {
1410  for (int j=0;j<numBfy;j++)
1411  {
1412  const int I = numBfy * i + j;
1413  vals(cell,field,I) = wtsx(i) * wtsy(j) * data(cell,field,I);
1414  }
1415  }
1416  }
1417  }
1418 
1419  return;
1420 
1421  }
1422 
1423  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1424  class ArrayTypeBasis, class ArrayTypeWeights>
1425  void TensorProductSpaceTools::moments3D( ArrayTypeOut &vals ,
1426  const ArrayTypeData &data ,
1427  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1428  const Array<RCP<ArrayTypeWeights> > &wts )
1429  {
1430  const int numBfx = basisVals[0]->dimension(0);
1431  const int numBfy = basisVals[1]->dimension(0);
1432  const int numBfz = basisVals[2]->dimension(0);
1433 
1434  const int numPtsx = basisVals[0]->dimension(1);
1435  const int numPtsy = basisVals[1]->dimension(1);
1436  const int numPtsz = basisVals[2]->dimension(1);
1437 
1438  const int numCells = vals.dimension(0);
1439  const int numFields = vals.dimension(1);
1440 
1441  const ArrayTypeBasis &Phix = *basisVals[0];
1442  const ArrayTypeBasis &Phiy = *basisVals[1];
1443  const ArrayTypeBasis &Phiz = *basisVals[2];
1444 
1445  const ArrayTypeWeights &Wtx = *wts[0];
1446  const ArrayTypeWeights &Wty = *wts[1];
1447  const ArrayTypeWeights &Wtz = *wts[2];
1448 
1449  FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy);
1450  FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz);
1451 
1452  for (int f=0;f<numFields;f++)
1453  {
1454  // Xi phase
1455  for (int c=0;c<numCells;c++)
1456  {
1457  for (int i=0;i<numBfx;i++)
1458  {
1459  for (int n=0;n<numPtsz;n++)
1460  {
1461  for (int m=0;m<numPtsy;m++)
1462  {
1463  for (int l=0;l<numPtsx;l++)
1464  {
1465  Xi(c,i,n,m) += Wtx(l) * Phix(i,l) * data(c,f,n*numPtsy*numPtsx+m*numPtsx+l);
1466  }
1467  }
1468  }
1469  }
1470  }
1471 
1472  // Theta phase
1473  for (int c=0;c<numCells;c++)
1474  {
1475  for (int j=0;j<numBfy;j++)
1476  {
1477  for (int i=0;i<numBfx;i++)
1478  {
1479  for (int n=0;n<numPtsz;n++)
1480  {
1481  for (int m=0;m<numPtsy;m++)
1482  {
1483  Theta(c,j,i,n) += Wty(m) * Phiy(j,m) * Xi(c,i,n,m);
1484  }
1485  }
1486  }
1487  }
1488  }
1489 
1490  // Final phase
1491  for (int c=0;c<numCells;c++)
1492  {
1493  for (int k=0;k<numBfz;k++)
1494  {
1495  for (int j=0;j<numBfy;j++)
1496  {
1497  for (int i=0;i<numBfx;i++)
1498  {
1499  const int momIdx = k*numBfx*numBfy+j*numBfx+i;
1500  for (int n=0;n<numPtsz;n++)
1501  {
1502  vals(c,f,momIdx) += Wtz(n) * Phiz(k,n) * Theta(c,j,i,n);
1503  }
1504  }
1505  }
1506  }
1507  }
1508 
1509  }
1510  return;
1511 
1512  }
1513 
1514  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1515  class ArrayTypeBasis, class ArrayTypeWeights>
1516  void TensorProductSpaceTools::momentsCollocated3D( ArrayTypeOut &vals ,
1517  const ArrayTypeData &data ,
1518  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1519  const Array<RCP<ArrayTypeWeights> > &wts )
1520  {
1521  const int numBfx = basisVals[0]->dimension(0);
1522  const int numBfy = basisVals[1]->dimension(0);
1523  const int numBfz = basisVals[2]->dimension(0);
1524 
1525  const int numPtsx = basisVals[0]->dimension(1);
1526  const int numPtsy = basisVals[1]->dimension(1);
1527  const int numPtsz = basisVals[2]->dimension(1);
1528 
1529  const int numCells = vals.dimension(0);
1530  const int numFields = vals.dimension(1);
1531 
1532  const ArrayTypeWeights &Wtx = *wts[0];
1533  const ArrayTypeWeights &Wty = *wts[1];
1534  const ArrayTypeWeights &Wtz = *wts[2];
1535 
1536  for (int cell=0;cell<numCells;cell++)
1537  {
1538  for (int field=0;field<numFields;field++)
1539  {
1540  for (int k=0;k<numBfz;k++)
1541  {
1542  for (int j=0;j<numBfy;j++)
1543  {
1544  for (int i=0;i<numBfx;i++)
1545  {
1546  const int I = k * numBfy * numBfx + j * numBfx + i;
1547  vals(cell,field,I) = Wtx(i) * Wty(j) * Wtz(k) * data(cell,field,I);
1548  }
1549  }
1550  }
1551  }
1552  }
1553 
1554  return;
1555 
1556  }
1557 
1558  // data is now (C,P,D)
1559  // want to compute the moments against the gradients of the basis
1560  // functions.
1561 
1562  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1563  class ArrayTypeBasis, class ArrayTypeWeights>
1564  void TensorProductSpaceTools::momentsGrad2D( ArrayTypeOut &vals ,
1565  const ArrayTypeData &data ,
1566  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1567  const Array<RCP<ArrayTypeBasis> > &Dbases ,
1568  const Array<RCP<ArrayTypeWeights> > &wts )
1569  {
1570 
1571  const int numBfx = basisVals[0]->dimension(0);
1572  const int numBfy = basisVals[1]->dimension(0);
1573  const int numPtsx = basisVals[0]->dimension(1);
1574  const int numPtsy = basisVals[1]->dimension(1);
1575  const int numCells = vals.dimension(0);
1576  const int numFields = vals.dimension(1);
1577  const ArrayTypeBasis &Phix = *basisVals[0];
1578  const ArrayTypeBasis &Phiy = *basisVals[1];
1579  const ArrayTypeBasis &DPhix = *Dbases[0];
1580  const ArrayTypeBasis &DPhiy = *Dbases[1];
1581  const ArrayTypeWeights &wtsx = *wts[0];
1582  const ArrayTypeWeights &wtsy = *wts[1];
1583 
1584  FieldContainer<double> Xi(numBfx,numPtsy,2);
1585 
1586  for (int f=0;f<numFields;f++)
1587  {
1588  // Xi phase
1589  for (int c=0;c<numCells;c++)
1590  {
1591  for (int i=0;i<numBfx;i++)
1592  {
1593  for (int m=0;m<numPtsy;m++)
1594  {
1595  for (int l=0;l<numPtsx;l++)
1596  {
1597  Xi(i,m,0) += wtsx(l) * data(c,f,m*numPtsy+l,0) * DPhix(i,l);
1598  Xi(i,m,1) += wtsx(l) * data(c,f,m*numPtsy+l,1) * Phix(i,l);
1599  }
1600  }
1601  }
1602  }
1603 
1604  // main phase
1605  for (int c=0;c<numCells;c++)
1606  {
1607  for (int j=0;j<numBfy;j++)
1608  {
1609  for (int i=0;i<numBfx;i++)
1610  {
1611  const int bfIdx = j*numBfx+i;
1612  vals(c,f,bfIdx) = 0.0;
1613  for (int m=0;m<numPtsy;m++)
1614  {
1615  vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,0) * Phiy(j,m);
1616  vals(c,f,bfIdx) += wtsy(m) * Xi(i,m,1) * DPhiy(j,m);
1617  }
1618  }
1619  }
1620  }
1621  }
1622  return;
1623 
1624  }
1625 
1626  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1627  class ArrayTypeBasis, class ArrayTypeWeights>
1628  void TensorProductSpaceTools::momentsGradCollocated2D( ArrayTypeOut &vals ,
1629  const ArrayTypeData &data ,
1630  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1631  const Array<RCP<ArrayTypeBasis> > &Dbases ,
1632  const Array<RCP<ArrayTypeWeights> > &wts )
1633  {
1634 
1635  const int numBfx = basisVals[0]->dimension(0);
1636  const int numBfy = basisVals[1]->dimension(0);
1637  const int numPtsx = basisVals[0]->dimension(1);
1638  const int numPtsy = basisVals[1]->dimension(1);
1639  const int numCells = vals.dimension(0);
1640  const int numFields = vals.dimension(1);
1641  const ArrayTypeBasis &Phix = *basisVals[0];
1642  const ArrayTypeBasis &Phiy = *basisVals[1];
1643  const ArrayTypeBasis &DPhix = *Dbases[0];
1644  const ArrayTypeBasis &DPhiy = *Dbases[1];
1645  const ArrayTypeWeights &wtsx = *wts[0];
1646  const ArrayTypeWeights &wtsy = *wts[1];
1647 
1648  for (int cell=0;cell<numCells;cell++)
1649  {
1650  for (int field=0;field<numFields;field++)
1651  {
1652  for (int j=0;j<numBfy;j++)
1653  {
1654  for (int i=0;i<numBfx;i++)
1655  {
1656  const int I=j*numBfx+i;
1657  vals(cell,field,I) = 0.0;
1658  }
1659  }
1660  }
1661  }
1662 
1663  for (int cell=0;cell<numCells;cell++)
1664  {
1665  for (int field=0;field<numFields;field++)
1666  {
1667  for (int j=0;j<numBfy;j++)
1668  {
1669  for (int i=0;i<numBfx;i++)
1670  {
1671  const int I=j*numBfx+i;
1672  for (int m=0;m<numPtsx;m++)
1673  {
1674  const int Itmp = j * numBfy + m;
1675  vals(cell,field,Itmp) += wtsx(m) * wtsy(j) * DPhix(i,m);
1676  }
1677  }
1678  }
1679  }
1680  }
1681 
1682  for (int cell=0;cell<numCells;cell++)
1683  {
1684  for (int field=0;field<numFields;field++)
1685  {
1686  for (int j=0;j<numBfy;j++)
1687  {
1688  for (int i=0;i<numBfx;i++)
1689  {
1690  const int I=j*numBfx+i;
1691  for (int n=0;n<numPtsy;n++)
1692  {
1693  const int Itmp = n * numBfy + i;
1694  vals(cell,field,Itmp) += wtsx(i) * wtsy(n) * DPhiy(j,n);
1695  }
1696  }
1697  }
1698  }
1699  }
1700 
1701  }
1702 
1703  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1704  class ArrayTypeBasis, class ArrayTypeWeights>
1705  void TensorProductSpaceTools::momentsGrad3D( ArrayTypeOut &vals ,
1706  const ArrayTypeData &data ,
1707  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1708  const Array<RCP<ArrayTypeBasis> > &basisDVals ,
1709  const Array<RCP<ArrayTypeWeights> > &wts )
1710  {
1711 
1712  const int numBfx = basisVals[0]->dimension(0);
1713  const int numBfy = basisVals[1]->dimension(0);
1714  const int numBfz = basisVals[2]->dimension(0);
1715  const int numPtsx = basisVals[0]->dimension(1);
1716  const int numPtsy = basisVals[1]->dimension(1);
1717  const int numPtsz = basisVals[2]->dimension(1);
1718  const int numCells = vals.dimension(0);
1719  const int numFields = vals.dimension(1);
1720  const ArrayTypeBasis &Phix = *basisVals[0];
1721  const ArrayTypeBasis &Phiy = *basisVals[1];
1722  const ArrayTypeBasis &Phiz = *basisVals[2];
1723  const ArrayTypeBasis &DPhix = *basisDVals[0];
1724  const ArrayTypeBasis &DPhiy = *basisDVals[1];
1725  const ArrayTypeBasis &DPhiz = *basisDVals[2];
1726  const ArrayTypeWeights &wtsx = *wts[0];
1727  const ArrayTypeWeights &wtsy = *wts[1];
1728  const ArrayTypeWeights &wtsz = *wts[2];
1729 
1730  FieldContainer<double> Xi(numCells,numBfx,numPtsz,numPtsy,3);
1731  FieldContainer<double> Theta(numCells,numBfy,numBfx,numPtsz,3);
1732 
1733  // Xi phase
1734  for (int f=0;f<numFields;f++)
1735  {
1736  for (int c=0;c<numCells;c++)
1737  {
1738  for (int i=0;i<numBfx;i++)
1739  {
1740  for (int n=0;n<numPtsz;n++)
1741  {
1742  for (int m=0;m<numPtsy;m++)
1743  {
1744  for (int l=0;l<numPtsx;l++)
1745  {
1746  const int dataIdx = n * numPtsy * numPtsx + m * numPtsx + l;
1747  Xi(c,i,n,m,0) += wtsx(l) * DPhix(i,l,0) * data(c,f,dataIdx,0);
1748  Xi(c,i,n,m,1) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,1);
1749  Xi(c,i,n,m,2) += wtsx(l) * Phix(i,l) * data(c,f,dataIdx,2);
1750  }
1751  }
1752  }
1753  }
1754  }
1755 
1756  // Theta phase
1757  for (int c=0;c<numCells;c++)
1758  {
1759  for (int j=0;j<numBfy;j++)
1760  {
1761  for (int i=0;i<numBfx;i++)
1762  {
1763  for (int n=0;n<numPtsz;n++)
1764  {
1765  for (int m=0;m<numPtsy;m++)
1766  {
1767  Theta(c,j,i,n,0) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,0);
1768  Theta(c,j,i,n,1) += wtsy(j) * DPhiy(j,m,0) * Xi(c,i,n,m,1);
1769  Theta(c,j,i,n,2) += wtsy(j) * Phiy(j,m) * Xi(c,i,n,m,2);
1770  }
1771  }
1772  }
1773  }
1774  }
1775 
1776  // last phase
1777  for (int c=0;c<numCells;c++)
1778  {
1779  for (int k=0;k<numBfz;k++)
1780  {
1781  for (int j=0;j<numBfy;j++)
1782  {
1783  for (int i=0;i<numBfx;i++)
1784  {
1785  const int momIdx = k * numBfx * numBfy + j * numBfx + i;
1786  for (int n=0;n<numPtsz;n++)
1787  {
1788  vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,0) * Phiz(k,n);
1789  vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,1) * Phiz(k,n);
1790  vals(c,f,momIdx) += wtsz(n) * Theta(c,j,i,n,2) * DPhiz(k,n,0);
1791  }
1792  }
1793  }
1794  }
1795  }
1796  }
1797 
1798 }
1799 
1800  template<class Scalar, class ArrayTypeOut, class ArrayTypeData,
1801  class ArrayTypeBasis, class ArrayTypeWeights>
1802  void TensorProductSpaceTools::momentsGradCollocated3D( ArrayTypeOut &vals ,
1803  const ArrayTypeData &data ,
1804  const Array<RCP<ArrayTypeBasis> > &basisVals ,
1805  const Array<RCP<ArrayTypeBasis> > &basisDVals ,
1806  const Array<RCP<ArrayTypeWeights> > &wts )
1807  {
1808 
1809  const int numBfx = basisVals[0]->dimension(0);
1810  const int numBfy = basisVals[1]->dimension(0);
1811  const int numBfz = basisVals[2]->dimension(0);
1812  const int numPtsx = basisVals[0]->dimension(1);
1813  const int numPtsy = basisVals[1]->dimension(1);
1814  const int numPtsz = basisVals[2]->dimension(1);
1815  const int numCells = vals.dimension(0);
1816  const int numFields = vals.dimension(1);
1817  const ArrayTypeBasis &Phix = *basisVals[0];
1818  const ArrayTypeBasis &Phiy = *basisVals[1];
1819  const ArrayTypeBasis &Phiz = *basisVals[2];
1820  const ArrayTypeBasis &DPhix = *basisDVals[0];
1821  const ArrayTypeBasis &DPhiy = *basisDVals[1];
1822  const ArrayTypeBasis &DPhiz = *basisDVals[2];
1823  const ArrayTypeWeights &wtsx = *wts[0];
1824  const ArrayTypeWeights &wtsy = *wts[1];
1825  const ArrayTypeWeights &wtsz = *wts[2];
1826 
1827  for (int cell=0;cell<numCells;cell++)
1828  {
1829  for (int field=0;field<numFields;field++)
1830  {
1831  // x component of data versus x derivative of bases
1832  for (int k=0;k<numBfz;k++)
1833  {
1834  for (int j=0;j<numBfy;j++)
1835  {
1836  for (int i=0;i<numBfx;i++)
1837  {
1838  const int I = numBfy * numBfx * k + numBfy * j + i;
1839  for (int ell=0;ell<numPtsx;ell++)
1840  {
1841  const int Itmp = numBfy * numBfx * k + numBfy * j + ell;
1842  vals(cell,field,I) += wtsx(ell) * wtsy(j) * wtsz(k) * DPhix( i , ell );
1843  }
1844  }
1845  }
1846  }
1847  // y component of data versus y derivative of bases
1848  for (int k=0;k<numBfz;k++)
1849  {
1850  for (int j=0;j<numBfy;j++)
1851  {
1852  for (int i=0;i<numBfx;i++)
1853  {
1854  const int I = numBfy * numBfx * k + numBfy * j + i;
1855  for (int m=0;m<numPtsy;m++)
1856  {
1857  const int Itmp = numBfy * numBfx * k + numBfy * m + i;
1858  vals(cell,field,I) += wtsx(i) * wtsy(m) * wtsz(k) * DPhiy( j , m );
1859  }
1860  }
1861  }
1862  }
1863  // z component of data versus z derivative of bases
1864  for (int k=0;k<numBfz;k++)
1865  {
1866  for (int j=0;j<numBfy;j++)
1867  {
1868  for (int i=0;i<numBfx;i++)
1869  {
1870  const int I = numBfy * numBfx * k + numBfy * j + i;
1871  for (int n=0;n<numPtsz;n++)
1872  {
1873  const int Itmp = numBfy * numBfx * n + numBfy * j + i;
1874  vals(cell,field,I) += wtsx(i) * wtsy(j) * wtsz(n) * DPhiz( k , n );
1875  }
1876  }
1877  }
1878  }
1879  }
1880  }
1881 
1882 }
1883 
1884 
1885 
1886 } // end namespace Intrepid
Intrepid::TensorProductSpaceTools::evaluateGradient
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...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:224
Intrepid::TensorProductSpaceTools::momentsGradCollocated
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...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:563
Intrepid::TensorProductSpaceTools::moments
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.
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:364
Intrepid::TensorProductSpaceTools::momentsGrad
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...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:494
Intrepid::TensorProductSpaceTools::evaluate
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...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:53
Intrepid::TensorProductSpaceTools::momentsCollocated
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, assuming that the basis nodes and integration points coincide.
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:429
Intrepid::TensorProductSpaceTools::evaluateCollocated
static void evaluateCollocated(ArrayTypeOut &vals, const ArrayTypeCoeffs &coeffs, const Array< RCP< ArrayTypeBasis > > &bases)
Computes point values of a set of array-valued polynomials expressed in a tensor product basis at out...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:108
Intrepid::TensorProductSpaceTools::evaluateGradientCollocated
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...
Definition: Intrepid_TensorProductSpaceToolsDef.hpp:294
Generated by   doxygen 1.8.13