Panzer  Version of the Day
Panzer_BasisValues2.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #include "PanzerDiscFE_config.hpp"
44 #include "Panzer_Traits.hpp"
45 #include "Panzer_BasisValues2.hpp"
46 
48 
49 #include "Intrepid2_FunctionSpaceTools.hpp"
50 
51 
52 namespace panzer {
53 
54 template <typename Scalar>
56 evaluateValues(const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
57  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
58  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
59  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
60 {
61  PHX::MDField<Scalar,Cell,IP> weighted_measure;
62  PHX::MDField<Scalar,Cell,NODE,Dim> vertex_coordinates;
63  build_weighted = false;
64  evaluateValues(cub_points,jac,jac_det,jac_inv,weighted_measure,vertex_coordinates,false);
65 }
66 
67 template <typename Scalar>
69 evaluateValues(const PHX::MDField<Scalar,IP,Dim,void,void,void,void,void,void> & cub_points,
70  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
71  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
72  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv,
73  const PHX::MDField<Scalar,Cell,IP> & weighted_measure,
74  const PHX::MDField<Scalar,Cell,NODE,Dim> & vertex_coordinates,
75  bool use_vertex_coordinates)
76 {
77  MDFieldArrayFactory af("",ddims_,true);
78 
79  int num_dim = basis_layout->dimension();
80 
81  // currently this just copies on the basis objects are converted
82  // in intrepid
83  evaluateReferenceValues(cub_points,compute_derivatives,use_vertex_coordinates);
84 
85  PureBasis::EElementSpace elmtspace = getElementSpace();
86  if(elmtspace==PureBasis::CONST ||
87  elmtspace==PureBasis::HGRAD) {
88  Intrepid2::FunctionSpaceTools::
89  HGRADtransformVALUE<Scalar>(basis_scalar,
90  basis_ref_scalar);
91 
92  if(build_weighted) {
93  Intrepid2::FunctionSpaceTools::
94  multiplyMeasure<Scalar>(weighted_basis_scalar,
95  weighted_measure,
96  basis_scalar);
97  }
98  }
99  else if(elmtspace==PureBasis::HCURL) {
100  Intrepid2::FunctionSpaceTools::
101  HCURLtransformVALUE<Scalar>(basis_vector,
102  jac_inv,
103  basis_ref_vector);
104 
105  if(build_weighted) {
106  Intrepid2::FunctionSpaceTools::
107  multiplyMeasure<Scalar>(weighted_basis_vector,
108  weighted_measure,
109  basis_vector);
110  }
111  }
112  else if(elmtspace==PureBasis::HDIV)
113  {
114  Intrepid2::FunctionSpaceTools::
115  HDIVtransformVALUE<Scalar>(basis_vector,
116  jac,
117  jac_det,
118  basis_ref_vector);
119 
120  if(build_weighted) {
121  Intrepid2::FunctionSpaceTools::
122  multiplyMeasure<Scalar>(weighted_basis_vector,
123  weighted_measure,
124  basis_vector);
125  }
126  }
127  else { TEUCHOS_ASSERT(false); }
128 
129  if(elmtspace==PureBasis::HGRAD && compute_derivatives) {
130  Intrepid2::FunctionSpaceTools::
131  HGRADtransformGRAD<Scalar>(grad_basis,
132  jac_inv,
133  grad_basis_ref);
134 
135  if(build_weighted) {
136  Intrepid2::FunctionSpaceTools::
137  multiplyMeasure<Scalar>(weighted_grad_basis,
138  weighted_measure,
139  grad_basis);
140  }
141  }
142  else if(elmtspace==PureBasis::HCURL && num_dim==2 && compute_derivatives) {
143  Intrepid2::FunctionSpaceTools::
144  HDIVtransformDIV<Scalar>(curl_basis_scalar,
145  jac_det, // note only volume deformation is needed!
146  // this relates directly to this being in
147  // the divergence space in 2D!
148  curl_basis_ref_scalar);
149 
150  if(build_weighted) {
151  Intrepid2::FunctionSpaceTools::
152  multiplyMeasure<Scalar>(weighted_curl_basis_scalar,
153  weighted_measure,
154  curl_basis_scalar);
155  }
156  }
157  else if(elmtspace==PureBasis::HCURL && num_dim==3 && compute_derivatives) {
158  Intrepid2::FunctionSpaceTools::
159  HCURLtransformCURL<Scalar>(curl_basis_vector,
160  jac,
161  jac_det,
162  curl_basis_ref_vector);
163 
164  if(build_weighted) {
165  Intrepid2::FunctionSpaceTools::
166  multiplyMeasure<Scalar>(weighted_curl_basis_vector,
167  weighted_measure,
168  curl_basis_vector);
169  }
170  }
171  else if(elmtspace==PureBasis::HDIV && compute_derivatives) {
172  Intrepid2::FunctionSpaceTools::
173  HDIVtransformDIV<Scalar>(div_basis,
174  jac_det,
175  div_basis_ref);
176 
177  if(build_weighted) {
178  Intrepid2::FunctionSpaceTools::
179  multiplyMeasure<Scalar>(weighted_div_basis,
180  weighted_measure,
181  div_basis);
182  }
183  }
184 
185  // If basis supports coordinate values at basis points, then
186  // compute these values
187  if(use_vertex_coordinates) {
188  Teuchos::RCP<Intrepid2::DofCoordsInterface<ArrayDynamic> > coords
189  = Teuchos::rcp_dynamic_cast<Intrepid2::DofCoordsInterface<ArrayDynamic> >(intrepid_basis);
190  if (!Teuchos::is_null(coords)) {
191 /*
192  ArrayDynamic dyn_basis_coordinates_ref = af.buildArray<Scalar,BASIS,Dim>("basis_coordinates_ref",basis_coordinates_ref.dimension(0),basis_coordinates_ref.dimension(1));
193  coords->getDofCoords(dyn_basis_coordinates_ref);
194 
195  // fill in basis coordinates
196  for (size_type i = 0; i < basis_coordinates_ref.dimension(0); ++i)
197  for (size_type j = 0; j < basis_coordinates_ref.dimension(1); ++j)
198  basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
199 */
200 
201  Intrepid2::CellTools<Scalar> cell_tools;
202  cell_tools.mapToPhysicalFrame(basis_coordinates,
203  basis_coordinates_ref,
204  vertex_coordinates,
205  intrepid_basis->getBaseCellTopology());
206  }
207  }
208 }
209 
210 template <typename Scalar>
212 evaluateValuesCV(const PHX::MDField<Scalar,Cell,IP,Dim,void,void,void,void,void> & cell_cub_points,
213  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac,
214  const PHX::MDField<Scalar,Cell,IP,void,void,void,void,void,void> & jac_det,
215  const PHX::MDField<Scalar,Cell,IP,Dim,Dim,void,void,void,void> & jac_inv)
216 {
217  MDFieldArrayFactory af("",ddims_,true);
218 
219  int num_ip = basis_layout->numPoints();
220  int num_card = basis_layout->cardinality();
221  int num_dim = basis_layout->dimension();
222 
223  size_type num_cells = jac.dimension(0);
224 
225  PureBasis::EElementSpace elmtspace = getElementSpace();
226  ArrayDynamic dyn_cub_points = af.buildArray<Scalar,IP,Dim>("dyn_cub_points", num_ip, num_dim);
227 
228  // Integration points are located on physical cells rather than reference cells,
229  // so we evaluate the basis in a loop over cells.
230  for (size_type icell = 0; icell < num_cells; ++icell)
231  {
232  for (int ip = 0; ip < num_ip; ++ip)
233  for (int d = 0; d < num_dim; ++d)
234  dyn_cub_points(ip,d) = cell_cub_points(icell,ip,d);
235 
236  if(elmtspace==PureBasis::CONST) {
237  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_ip);
238 
239  intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
240  Intrepid2::OPERATOR_VALUE);
241 
242  // transform values method just transfers values to array with cell index - no need to call
243  for (int b = 0; b < num_card; ++b)
244  for (int ip = 0; ip < num_ip; ++ip)
245  basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
246 
247  }
248  if(elmtspace==PureBasis::HGRAD) {
249  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_ip);
250 
251  intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
252  Intrepid2::OPERATOR_VALUE);
253 
254  // transform values method just transfers values to array with cell index - no need to call
255  for (int b = 0; b < num_card; ++b)
256  for (int ip = 0; ip < num_ip; ++ip)
257  basis_scalar(icell,b,ip) = dyn_basis_ref_scalar(b,ip);
258 
259  if(compute_derivatives) {
260 
261  int one_cell = 1;
262  ArrayDynamic dyn_grad_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_grad_basis_ref",num_card,num_ip,num_dim);
263  ArrayDynamic dyn_grad_basis = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_grad_basis",one_cell,num_card,num_ip,num_dim);
264  ArrayDynamic dyn_jac_inv = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim);
265 
266  intrepid_basis->getValues(dyn_grad_basis_ref, dyn_cub_points,
267  Intrepid2::OPERATOR_GRAD);
268 
269  int cellInd = 0;
270  for (int ip = 0; ip < num_ip; ++ip)
271  for (int d1 = 0; d1 < num_dim; ++d1)
272  for (int d2 = 0; d2 < num_dim; ++d2)
273  dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
274 
275  Intrepid2::FunctionSpaceTools::HGRADtransformGRAD<Scalar>(dyn_grad_basis,
276  dyn_jac_inv,
277  dyn_grad_basis_ref);
278 
279  for (int b = 0; b < num_card; ++b)
280  for (int ip = 0; ip < num_ip; ++ip)
281  for (int d = 0; d < num_dim; ++d)
282  grad_basis(icell,b,ip,d) = dyn_grad_basis(0,b,ip,d);
283 
284  }
285  }
286  else if(elmtspace==PureBasis::HCURL) {
287  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_ip,num_dim);
288 
289  intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
290  Intrepid2::OPERATOR_VALUE);
291 
292  int one_cell = 1;
293  ArrayDynamic dyn_basis_vector = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_basis_vector",one_cell,num_card,num_ip,num_dim);
294  ArrayDynamic dyn_jac_inv = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac_inv",one_cell,num_ip,num_dim,num_dim);
295 
296  int cellInd = 0;
297  for (int ip = 0; ip < num_ip; ++ip)
298  for (int d1 = 0; d1 < num_dim; ++d1)
299  for (int d2 = 0; d2 < num_dim; ++d2)
300  dyn_jac_inv(cellInd,ip,d1,d2) = jac_inv(icell,ip,d1,d2);
301 
302  Intrepid2::FunctionSpaceTools::HCURLtransformVALUE<Scalar>(dyn_basis_vector,
303  dyn_jac_inv,
304  dyn_basis_ref_vector);
305 
306  for (int b = 0; b < num_card; ++b)
307  for (int ip = 0; ip < num_ip; ++ip)
308  for (int d = 0; d < num_dim; ++d)
309  basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
310 
311  if(compute_derivatives && num_dim ==2) {
312 
313  int one_cell = 1;
314  ArrayDynamic dyn_curl_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_ip);
315  ArrayDynamic dyn_curl_basis_scalar = af.buildArray<Scalar,Cell,BASIS,IP>("dyn_curl_basis_scalar",one_cell,num_card,num_ip);
316  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
317 
318  intrepid_basis->getValues(dyn_curl_basis_ref_scalar, dyn_cub_points,
319  Intrepid2::OPERATOR_CURL);
320 
321  int cellInd = 0;
322  for (int ip = 0; ip < num_ip; ++ip)
323  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
324 
325  Intrepid2::FunctionSpaceTools::HDIVtransformDIV<Scalar>(dyn_curl_basis_scalar,
326  dyn_jac_det,
327  dyn_curl_basis_ref_scalar);
328 
329  for (int b = 0; b < num_card; ++b)
330  for (int ip = 0; ip < num_ip; ++ip)
331  curl_basis_scalar(icell,b,ip) = dyn_curl_basis_scalar(0,b,ip);
332 
333  }
334  if(compute_derivatives && num_dim ==3) {
335 
336  int one_cell = 1;
337  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_ip,num_dim);
338  ArrayDynamic dyn_curl_basis = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_curl_basis_vector",one_cell,num_card,num_ip,num_dim);
339  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
340  ArrayDynamic dyn_jac = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac",one_cell,num_ip,num_dim,num_dim);
341 
342  intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
343  Intrepid2::OPERATOR_CURL);
344 
345  int cellInd = 0;
346  for (int ip = 0; ip < num_ip; ++ip)
347  {
348  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
349  for (int d1 = 0; d1 < num_dim; ++d1)
350  for (int d2 = 0; d2 < num_dim; ++d2)
351  dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2);
352  }
353 
354  Intrepid2::FunctionSpaceTools::HCURLtransformCURL<Scalar>(dyn_curl_basis,
355  dyn_jac,
356  dyn_jac_det,
357  dyn_curl_basis_ref);
358 
359  for (int b = 0; b < num_card; ++b)
360  for (int ip = 0; ip < num_ip; ++ip)
361  for (int d = 0; d < num_dim; ++d)
362  curl_basis_vector(icell,b,ip,d) = dyn_curl_basis(0,b,ip,d);
363 
364  }
365 
366  }
367  else if(elmtspace==PureBasis::HDIV) {
368 
369  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_ip,num_dim);
370 
371  intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
372  Intrepid2::OPERATOR_VALUE);
373 
374  int one_cell= 1;
375  ArrayDynamic dyn_basis_vector = af.buildArray<Scalar,Cell,BASIS,IP,Dim>("dyn_basis_vector",one_cell,num_card,num_ip,num_dim);
376  ArrayDynamic dyn_jac = af.buildArray<Scalar,Cell,IP,Dim,Dim>("dyn_jac",one_cell,num_ip,num_dim,num_dim);
377  ArrayDynamic dyn_jac_det = af.buildArray<Scalar,Cell,IP>("dyn_jac_det",one_cell,num_ip);
378 
379  int cellInd = 0;
380  for (int ip = 0; ip < num_ip; ++ip)
381  {
382  dyn_jac_det(cellInd,ip) = jac_det(icell,ip);
383  for (int d1 = 0; d1 < num_dim; ++d1)
384  for (int d2 = 0; d2 < num_dim; ++d2)
385  dyn_jac(cellInd,ip,d1,d2) = jac(icell,ip,d1,d2);
386  }
387 
388  Intrepid2::FunctionSpaceTools::HDIVtransformVALUE<Scalar>(dyn_basis_vector,
389  dyn_jac,dyn_jac_det,
390  dyn_basis_ref_vector);
391 
392  for (int b = 0; b < num_card; ++b)
393  for (int ip = 0; ip < num_ip; ++ip)
394  for (int d = 0; d < num_dim; ++d)
395  basis_vector(icell,b,ip,d) = dyn_basis_vector(0,b,ip,d);
396 
397  if(compute_derivatives) {
398 
399  ArrayDynamic dyn_div_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_ip);
400  ArrayDynamic dyn_div_basis = af.buildArray<Scalar,Cell,BASIS,IP>("dyn_div_basis_scalar",one_cell,num_card,num_ip);
401 
402  intrepid_basis->getValues(dyn_div_basis_ref, dyn_cub_points,
403  Intrepid2::OPERATOR_DIV);
404 
405  Intrepid2::FunctionSpaceTools::HDIVtransformDIV<Scalar>(dyn_div_basis,
406  dyn_jac_det,
407  dyn_div_basis_ref);
408 
409  for (int b = 0; b < num_card; ++b)
410  for (int ip = 0; ip < num_ip; ++ip)
411  div_basis(icell,b,ip) = dyn_div_basis(0,b,ip);
412 
413  }
414 
415  }
416  else { TEUCHOS_ASSERT(false); }
417 
418  } // cell loop
419 
420 }
421 
422 template <typename Scalar>
424 evaluateReferenceValues(const PHX::MDField<Scalar,IP,Dim> & cub_points,bool compute_derivatives,bool use_vertex_coordinates)
425 {
426  MDFieldArrayFactory af("",ddims_,true);
427 
428  int num_quad = basis_layout->numPoints();
429  int num_dim = basis_layout->dimension();
430  int num_card = basis_layout->cardinality();
431 
432  ArrayDynamic dyn_cub_points = af.buildArray<Scalar,IP,Dim>("dyn_cub_points", num_quad,num_dim);
433 
434  for (int ip = 0; ip < num_quad; ++ip)
435  for (int d = 0; d < num_dim; ++d)
436  dyn_cub_points(ip,d) = cub_points(ip,d);
437 
438  PureBasis::EElementSpace elmtspace = getElementSpace();
439  if(elmtspace==PureBasis::HGRAD || elmtspace==PureBasis::CONST) {
440  ArrayDynamic dyn_basis_ref_scalar = af.buildArray<Scalar,BASIS,IP>("dyn_basis_ref_scalar",num_card,num_quad);
441 
442  intrepid_basis->getValues(dyn_basis_ref_scalar, dyn_cub_points,
443  Intrepid2::OPERATOR_VALUE);
444 
445  for (int b = 0; b < num_card; ++b)
446  for (int ip = 0; ip < num_quad; ++ip)
447  basis_ref_scalar(b,ip) = dyn_basis_ref_scalar(b,ip);
448  }
449  else if(elmtspace==PureBasis::HDIV || elmtspace==PureBasis::HCURL) {
450  ArrayDynamic dyn_basis_ref_vector = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
451 
452  intrepid_basis->getValues(dyn_basis_ref_vector, dyn_cub_points,
453  Intrepid2::OPERATOR_VALUE);
454 
455  for (int b = 0; b < num_card; ++b)
456  for (int ip = 0; ip < num_quad; ++ip)
457  for (int d = 0; d < num_dim; ++d)
458  basis_ref_vector(b,ip,d) = dyn_basis_ref_vector(b,ip,d);
459  }
460  else { TEUCHOS_ASSERT(false); }
461 
462  if(elmtspace==PureBasis::HGRAD && compute_derivatives) {
463  ArrayDynamic dyn_grad_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_basis_ref_vector",num_card,num_quad,num_dim);
464 
465  intrepid_basis->getValues(dyn_grad_basis_ref, dyn_cub_points,
466  Intrepid2::OPERATOR_GRAD);
467 
468  for (int b = 0; b < num_card; ++b)
469  for (int ip = 0; ip < num_quad; ++ip)
470  for (int d = 0; d < num_dim; ++d)
471  grad_basis_ref(b,ip,d) = dyn_grad_basis_ref(b,ip,d);
472  }
473  else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==2) {
474  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_curl_basis_ref_scalar",num_card,num_quad);
475 
476  intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
477  Intrepid2::OPERATOR_CURL);
478 
479  for (int b = 0; b < num_card; ++b)
480  for (int ip = 0; ip < num_quad; ++ip)
481  curl_basis_ref_scalar(b,ip) = dyn_curl_basis_ref(b,ip);
482  }
483  else if(elmtspace==PureBasis::HCURL && compute_derivatives && num_dim==3) {
484  ArrayDynamic dyn_curl_basis_ref = af.buildArray<Scalar,BASIS,IP,Dim>("dyn_curl_basis_ref_vector",num_card,num_quad,num_dim);
485 
486  intrepid_basis->getValues(dyn_curl_basis_ref, dyn_cub_points,
487  Intrepid2::OPERATOR_CURL);
488 
489  for (int b = 0; b < num_card; ++b)
490  for (int ip = 0; ip < num_quad; ++ip)
491  for (int d = 0; d < num_dim; ++d)
492  curl_basis_ref_vector(b,ip,d) = dyn_curl_basis_ref(b,ip,d);
493  }
494  else if(elmtspace==PureBasis::HDIV && compute_derivatives) {
495  ArrayDynamic dyn_div_basis_ref = af.buildArray<Scalar,BASIS,IP>("dyn_div_basis_ref_scalar",num_card,num_quad);
496 
497  intrepid_basis->getValues(dyn_div_basis_ref, dyn_cub_points,
498  Intrepid2::OPERATOR_DIV);
499 
500  for (int b = 0; b < num_card; ++b)
501  for (int ip = 0; ip < num_quad; ++ip)
502  div_basis_ref(b,ip) = dyn_div_basis_ref(b,ip);
503  }
504 
505 
506  if(use_vertex_coordinates) {
507  Teuchos::RCP<Intrepid2::DofCoordsInterface<ArrayDynamic> > coords
508  = Teuchos::rcp_dynamic_cast<Intrepid2::DofCoordsInterface<ArrayDynamic> >(intrepid_basis);
509  // This will trip if the basis is not inheirited off DofCoordsInterface and doesnt have the right function,getDofCoords
510  TEUCHOS_ASSERT(elmtspace ==PureBasis::CONST || !Teuchos::is_null(coords));
511  if (!Teuchos::is_null(coords)) {
512  ArrayDynamic dyn_basis_coordinates_ref = af.buildArray<Scalar,BASIS,Dim>("basis_coordinates_ref",basis_coordinates_ref.dimension(0),basis_coordinates_ref.dimension(1));
513  coords->getDofCoords(dyn_basis_coordinates_ref);
514 
515  // fill in basis coordinates
516  for (int i = 0; i < basis_coordinates_ref.extent_int(0); ++i)
517  for (int j = 0; j < basis_coordinates_ref.extent_int(1); ++j)
518  basis_coordinates_ref(i,j) = dyn_basis_coordinates_ref(i,j);
519  }
520  }
521 
522  references_evaluated = true;
523 }
524 
525 // method for applying orientations
526 template <typename Scalar>
528 applyOrientations(const PHX::MDField<Scalar,Cell,BASIS> & orientations)
529 {
530  int num_cell = orientations.dimension(0);
531  int num_basis = orientations.dimension(1);
532  int num_dim = basis_layout->dimension();
533  int num_ip = basis_layout->numPoints();
534  PureBasis::EElementSpace elmtspace = getElementSpace();
535 
536  if(elmtspace==PureBasis::HCURL && num_dim==2) {
537 
538  // setup the orientations for the trial space
539  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
540 
541  for (int c=0; c<num_cell; c++)
542  for (int b=0; b<num_basis; b++)
543  for (int p=0; p<num_ip; p++)
544  for (int d=0; d<num_dim; d++)
545  basis_vector(c, b, p, d) *= orientations(c, b);
546 
547  if(compute_derivatives) {
548  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(curl_basis_scalar,orientations);
549  for (int c=0; c<num_cell; c++)
550  for (int b=0; b<num_basis; b++)
551  for (int p=0; p<num_ip; p++)
552  curl_basis_scalar(c, b, p) *= orientations(c, b);
553  }
554 
555  // setup the orientations for the test space
556  if(build_weighted) {
557  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_basis_vector,orientations);
558 
559  if(compute_derivatives)
560  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_curl_basis_scalar,orientations);
561  }
562  }
563  else if(elmtspace==PureBasis::HCURL && num_dim==3) {
564 
565  // setup the orientations for the trial space
566  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
567 
568  for (int c=0; c<num_cell; c++)
569  for (int b=0; b<num_basis; b++)
570  for (int p=0; p<num_ip; p++)
571  for (int d=0; d<num_dim; d++)
572  basis_vector(c, b, p, d) *= orientations(c, b);
573 
574  if(compute_derivatives) {
575  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(curl_basis_vector,orientations);
576  for (int c=0; c<num_cell; c++)
577  for (int b=0; b<num_basis; b++)
578  for (int p=0; p<num_ip; p++)
579  for (int d=0; d<num_dim; d++)
580  curl_basis_vector(c, b, p,d) *= orientations(c, b);
581  }
582 
583  // setup the orientations for the test space
584  if(build_weighted) {
585  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_basis_vector,orientations);
586 
587  if(compute_derivatives)
588  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_curl_basis_vector,orientations);
589  }
590  }
591  else if(elmtspace==PureBasis::HDIV) {
592  // setup the orientations for the trial space
593  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(basis_vector,orientations);
594 
595  for (int c=0; c<num_cell; c++)
596  for (int b=0; b<num_basis; b++)
597  for (int p=0; p<num_ip; p++)
598  for (int d=0; d<num_dim; d++)
599  basis_vector(c, b, p, d) *= orientations(c, b);
600 
601  if(compute_derivatives) {
602  // Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(div_basis,orientations);
603 
604  for (int c=0; c<num_cell; c++)
605  for (int b=0; b<num_basis; b++)
606  for (int p=0; p<num_ip; p++)
607  div_basis(c, b, p) *= orientations(c, b);
608  }
609 
610  // setup the orientations for the test space
611  if(build_weighted) {
612  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_basis_vector,orientations);
613 
614  if(compute_derivatives)
615  Intrepid2::FunctionSpaceTools::applyFieldSigns<Scalar>(weighted_div_basis,orientations);
616  }
617  }
618 }
619 
620 template <typename Scalar>
622 { return basis_layout->getBasis()->getElementSpace(); }
623 
624 template <typename Scalar>
626 setupArrays(const Teuchos::RCP<const panzer::BasisIRLayout>& layout,
627  bool computeDerivatives)
628 {
629  MDFieldArrayFactory af(prefix,alloc_arrays);
630 
631  compute_derivatives = computeDerivatives;
632  basis_layout = layout;
633  Teuchos::RCP<const panzer::PureBasis> basisDesc = layout->getBasis();
634 
635  // for convience pull out basis and quadrature information
636  int num_quad = layout->numPoints();
637  int dim = basisDesc->dimension();
638  int card = basisDesc->cardinality();
639  int numcells = basisDesc->numCells();
640  panzer::PureBasis::EElementSpace elmtspace = basisDesc->getElementSpace();
641  Teuchos::RCP<const shards::CellTopology> cellTopo = basisDesc->getCellTopology();
642 
643  intrepid_basis = basisDesc->getIntrepid2Basis<Scalar,ArrayDynamic>();
644 
645  // allocate field containers
646  // field sizes defined by http://trilinos.sandia.gov/packages/docs/dev/packages/intrepid/doc/html/basis_page.html#basis_md_array_sec
647 
648  // compute basis fields
649  if(elmtspace==panzer::PureBasis::HGRAD) {
650  // HGRAD is a nodal field
651 
652  // build values
654  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
655  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
656 
657  if(build_weighted)
658  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
659 
660  // build gradients
662 
663  if(compute_derivatives) {
664  grad_basis_ref = af.buildStaticArray<Scalar,BASIS,IP,Dim>("grad_basis_ref",card,num_quad,dim); // F, P, D
665  grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("grad_basis",numcells,card,num_quad,dim);
666 
667  if(build_weighted)
668  weighted_grad_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_grad_basis",numcells,card,num_quad,dim);
669  }
670 
671  // build curl
673 
674  // nothing - HGRAD does not support CURL operation
675  }
676  else if(elmtspace==panzer::PureBasis::HCURL) {
677  // HCURL is a vector field
678 
679  // build values
681 
682  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
683  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
684 
685  if(build_weighted)
686  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
687 
688  // build gradients
690 
691  // nothing - HCURL does not support GRAD operation
692 
693  // build curl
695 
696  if(compute_derivatives) {
697  if(dim==2) {
698  // curl of HCURL basis is not dimension dependent
699  curl_basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("curl_basis_ref",card,num_quad); // F, P
700  curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("curl_basis",numcells,card,num_quad);
701 
702  if(build_weighted)
703  weighted_curl_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_curl_basis",numcells,card,num_quad);
704  }
705  else if(dim==3){
706  curl_basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("curl_basis_ref",card,num_quad,dim); // F, P, D
707  curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("curl_basis",numcells,card,num_quad,dim);
708 
709  if(build_weighted)
710  weighted_curl_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_curl_basis",numcells,card,num_quad,dim);
711  }
712  else { TEUCHOS_ASSERT(false); } // what do I do with 1D?
713  }
714  }
715  else if(elmtspace==panzer::PureBasis::HDIV) {
716  // HDIV is a vector field
717 
718  // build values
720 
721  basis_ref_vector = af.buildStaticArray<Scalar,BASIS,IP,Dim>("basis_ref",card,num_quad,dim); // F, P, D
722  basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("basis",numcells,card,num_quad,dim);
723 
724  if(build_weighted)
725  weighted_basis_vector = af.buildStaticArray<Scalar,Cell,BASIS,IP,Dim>("weighted_basis",numcells,card,num_quad,dim);
726 
727  // build gradients
729 
730  // nothing - HCURL does not support GRAD operation
731 
732  // build curl
734 
735  // nothing - HDIV does not support CURL operation
736 
737  // build div
739 
740  if(compute_derivatives) {
741  div_basis_ref = af.buildStaticArray<Scalar,BASIS,IP>("div_basis_ref",card,num_quad); // F, P
742  div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("div_basis",numcells,card,num_quad);
743 
744  if(build_weighted)
745  weighted_div_basis = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_div_basis",numcells,card,num_quad);
746  }
747  }
748  else if(elmtspace==panzer::PureBasis::CONST) {
749  // CONST is a nodal field
750 
751  // build values
753  basis_ref_scalar = af.buildStaticArray<Scalar,BASIS,IP>("basis_ref",card,num_quad); // F, P
754  basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("basis",numcells,card,num_quad);
755 
756  if(build_weighted)
757  weighted_basis_scalar = af.buildStaticArray<Scalar,Cell,BASIS,IP>("weighted_basis",numcells,card,num_quad);
758 
759  // build gradients
761 
762  // nothing - CONST does not support GRAD operation
763 
764  // build curl
766 
767  // nothing - CONST does not support CURL operation
768 
769  // build div
771 
772  // nothing - CONST does not support DIV operation
773  }
774  else { TEUCHOS_ASSERT(false); }
775 
776  basis_coordinates_ref = af.buildStaticArray<Scalar,BASIS,Dim>("basis_coordinates_ref",card,dim);
777  basis_coordinates = af.buildStaticArray<Scalar,Cell,BASIS,Dim>("basis_coordinates",numcells,card,dim);
778 }
779 
780 // do some explicit instantiation so things compile faster.
781 
782 #define BASIS_VALUES_INSTANTIATION(SCALAR) \
783 template class BasisValues2<SCALAR>;
784 
787 #ifdef Panzer_BUILD_HESSIAN_SUPPORT
788 BASIS_VALUES_INSTANTIATION(panzer::Traits::HessianType)
789 #endif
790 
791 } // namespace panzer
PHX::MDField< Scalar, T0 > buildStaticArray(const std::string &str, int d0) const
void evaluateValues(const PHX::MDField< Scalar, IP, Dim, void, void, void, void, void, void > &cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
void evaluateReferenceValues(const PHX::MDField< Scalar, IP, Dim > &cub_points, bool compute_derivatives, bool use_vertex_coordinates)
ArrayTraits< Scalar, PHX::MDField< Scalar, void, void, void, void, void, void, void, void > >::size_type size_type
PHX::MDField< Scalar > buildArray(const std::string &str, int d0) const
PHX::MDField< double, Cell, BASIS, IP, Dim > weighted_basis_vector
#define BASIS_VALUES_INSTANTIATION(SCALAR)
PANZER_FADTYPE FadType
void setupArrays(const Teuchos::RCP< const panzer::BasisIRLayout > &basis, bool computeDerivatives=true)
Sizes/allocates memory for arrays.
void applyOrientations(const PHX::MDField< Scalar, Cell, BASIS > &orientations)
Method to apply orientations to a basis values container.
void evaluateValuesCV(const PHX::MDField< Scalar, Cell, IP, Dim, void, void, void, void, void > &cell_cub_points, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac, const PHX::MDField< Scalar, Cell, IP, void, void, void, void, void, void > &jac_det, const PHX::MDField< Scalar, Cell, IP, Dim, Dim, void, void, void, void > &jac_inv)
PureBasis::EElementSpace getElementSpace() const
PHX::MDField< Scalar > ArrayDynamic