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 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
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.5