Intrepid2
Intrepid2_OrientationToolsDefMatrixData.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 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 Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
43 
48 #ifndef __INTREPID2_ORIENTATIONTOOLS_DEF_MATRIX_DATA_HPP__
49 #define __INTREPID2_ORIENTATIONTOOLS_DEF_MATRIX_DATA_HPP__
50 
51 // disable clang warnings
52 #if defined (__clang__) && !defined (__INTEL_COMPILER)
53 #pragma clang system_header
54 #endif
55 
56 namespace Intrepid2 {
57 
58  template<typename SpT>
59  template<typename BasisPtrType>
60  typename OrientationTools<SpT>::CoeffMatrixDataViewType
61  OrientationTools<SpT>::createCoeffMatrixInternal(BasisPtrType basis) {
62  const ordinal_type order(basis->getDegree());
63  const std::string name(basis->getName());
64  CoeffMatrixDataViewType matData;
65 
66  auto ordinalToTag = basis->getAllDofTags();
67  auto tagToOrdinal = basis->getAllDofOrdinal();
68 
69  //
70  // High order HGRAD Elements
71  //
72 
73  if (name == "Intrepid2_HGRAD_QUAD_Cn_FEM") {
74  if (order >1) {
75  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 4, numOrts = 2;
76  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HGRAD_QUAD_Cn_FEM",
77  numEdges,
78  numOrts,
79  matDim,
80  matDim);
81 
82  init_HGRAD_QUAD_Cn_FEM(matData, order);
83  } else {
84  // add dummy
85  matData = CoeffMatrixDataViewType();
86  }
87  }
88  else if (name == "Intrepid2_HGRAD_HEX_Cn_FEM") {
89  if (order > 1) {
90  const ordinal_type matDim = ordinalToTag(tagToOrdinal(2, 0, 0), 3), numSubcells = 18, numOrts = 8;
91  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HGRAD_HEX_Cn_FEM",
92  numSubcells,
93  numOrts,
94  matDim,
95  matDim);
96 
97  init_HGRAD_HEX_Cn_FEM(matData, order);
98  } else {
99  // add dummy
100  matData = CoeffMatrixDataViewType();
101  }
102  }
103  else if (name == "Intrepid2_HGRAD_TRI_Cn_FEM") {
104  if (order > 1) {
105  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 3, numOrts = 2;
106  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HGRAD_TRI_Cn_FEM",
107  numEdges,
108  numOrts,
109  matDim,
110  matDim);
111 
112  init_HGRAD_TRI_Cn_FEM(matData, order);
113  } else {
114  // add dummy
115  matData = CoeffMatrixDataViewType();
116  }
117  }
118  else if (name == "Intrepid2_HGRAD_TET_Cn_FEM") {
119  const ordinal_type
120  matDim = max(ordinalToTag(tagToOrdinal(1, 0, 0), 3),
121  ordinalToTag(tagToOrdinal(2, 0, 0), 3)),
122  numSubcells = 10, numOrts = 6;
123  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HGRAD_TET_Cn_FEM",
124  numSubcells,
125  numOrts,
126  matDim,
127  matDim);
128 
129  init_HGRAD_TET_Cn_FEM(matData, order);
130  }
131 
132  //
133  // High order HCURL Elements
134  //
135 
136  else if (name == "Intrepid2_HCURL_QUAD_In_FEM") {
137  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 4, numOrts = 2;
138  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_QUAD_In_FEM",
139  numEdges,
140  numOrts,
141  matDim,
142  matDim);
143 
144  init_HCURL_QUAD_In_FEM(matData, order);
145  }
146  else if (name == "Intrepid2_HCURL_HEX_In_FEM") {
147  // if order == 1, there is no face dofs
148  const ordinal_type matDim = ordinalToTag(tagToOrdinal((order < 2 ? 1 : 2), 0, 0), 3), numSubcells = 18, numOrts = 8;
149  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_HEX_In_FEM",
150  numSubcells,
151  numOrts,
152  matDim,
153  matDim);
154 
155  init_HCURL_HEX_In_FEM(matData, order);
156  }
157  else if (name == "Intrepid2_HCURL_TRI_In_FEM") {
158  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 3, numOrts = 2;
159  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_TRI_In_FEM",
160  numEdges,
161  numOrts,
162  matDim,
163  matDim);
164 
165  init_HCURL_TRI_In_FEM(matData, order);
166  }
167  else if (name == "Intrepid2_HCURL_TET_In_FEM") {
168  const ordinal_type matDim = ordinalToTag(tagToOrdinal((order < 2 ? 1 : 2), 0, 0), 3), numSubcells = 10, numOrts = 6;
169  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_TET_In_FEM",
170  numSubcells,
171  numOrts,
172  matDim,
173  matDim);
174 
175  init_HCURL_TET_In_FEM(matData, order);
176  }
177 
178  //
179  // High order HDIV Elements
180  //
181 
182  else if (name == "Intrepid2_HDIV_QUAD_In_FEM") {
183  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 4, numOrts = 2;
184  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_QUAD_In_FEM",
185  numEdges,
186  numOrts,
187  matDim,
188  matDim);
189 
190  init_HDIV_QUAD_In_FEM(matData, order);
191  }
192  else if (name == "Intrepid2_HDIV_HEX_In_FEM") {
193  const ordinal_type matDim = ordinalToTag(tagToOrdinal(2, 0, 0), 3), numSubcells = 6, numOrts = 8;
194  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_HEX_In_FEM",
195  numSubcells,
196  numOrts,
197  matDim,
198  matDim);
199 
200  init_HDIV_HEX_In_FEM(matData, order);
201  }
202  else if (name == "Intrepid2_HDIV_TRI_In_FEM") {
203  const ordinal_type matDim = ordinalToTag(tagToOrdinal(1, 0, 0), 3), numEdges = 3, numOrts = 2;
204  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_TRI_In_FEM",
205  numEdges,
206  numOrts,
207  matDim,
208  matDim);
209 
210  init_HDIV_TRI_In_FEM(matData, order);
211  }
212  else if (name == "Intrepid2_HDIV_TET_In_FEM") {
213  const ordinal_type matDim = ordinalToTag(tagToOrdinal(2, 0, 0), 3), numSubcells = 4, numOrts = 6;
214  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_TET_In_FEM",
215  numSubcells,
216  numOrts,
217  matDim,
218  matDim);
219 
220  init_HDIV_TET_In_FEM(matData, order);
221  }
222 
223  //
224  // 2D H(Curl/Div) I1 Elements
225  //
226 
227  else if (name == "Intrepid2_HCURL_QUAD_I1_FEM" ||
228  name == "Intrepid2_HDIV_QUAD_I1_FEM") {
229  const ordinal_type matDim = 1, numEdges = 4, numOrts = 2;
230  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_QUAD_I1_FEM",
231  numEdges,
232  numOrts,
233  matDim,
234  matDim);
235  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId)
236  init_EDGE_ELEMENT_I1_FEM(matData, edgeId);
237  }
238  else if (name == "Intrepid2_HCURL_TRI_I1_FEM" ||
239  name == "Intrepid2_HDIV_TRI_I1_FEM") {
240  const ordinal_type matDim = 1, numEdges = 3, numOrts = 2;
241  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_TRI_I1_FEM",
242  numEdges,
243  numOrts,
244  matDim,
245  matDim);
246  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId)
247  init_EDGE_ELEMENT_I1_FEM(matData, edgeId);
248  }
249 
250  //
251  // 3D H(Curl) I1 Elements
252  //
253 
254  else if (name == "Intrepid2_HCURL_HEX_I1_FEM") {
255  const ordinal_type matDim = 1, numEdges = 12, numOrts = 2;
256  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_HEX_I1_FEM",
257  numEdges,
258  numOrts,
259  matDim,
260  matDim);
261  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId)
262  init_EDGE_ELEMENT_I1_FEM(matData, edgeId);
263  }
264  else if (name == "Intrepid2_HCURL_TET_I1_FEM") {
265  const ordinal_type matDim = 1, numEdges = 6, numOrts = 2;
266  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_TET_I1_FEM",
267  numEdges,
268  numOrts,
269  matDim,
270  matDim);
271  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId)
272  init_EDGE_ELEMENT_I1_FEM(matData, edgeId);
273  }
274  else if (name == "Intrepid2_HCURL_WEDGE_I1_FEM") {
275  const ordinal_type matDim = 1, numEdges = 9, numOrts = 2;
276  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HCURL_WEDGE_I1_FEM",
277  numEdges,
278  numOrts,
279  matDim,
280  matDim);
281  for (ordinal_type edgeId=0;edgeId<numEdges;++edgeId)
282  init_EDGE_ELEMENT_I1_FEM(matData, edgeId);
283  }
284 
285  //
286  // 3D H(Div) I1 Elements
287  //
288 
289  else if (name == "Intrepid2_HDIV_HEX_I1_FEM") {
290  const ordinal_type matDim = 1, numFaces = 6, numOrts = 8;
291  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_HEX_I1_FEM",
292  numFaces,
293  numOrts,
294  matDim,
295  matDim);
296  for (ordinal_type faceId=0;faceId<numFaces;++faceId)
297  init_QUAD_FACE_ELEMENT_I1_FEM(matData, faceId);
298  }
299  else if (name == "Intrepid2_HDIV_TET_I1_FEM") {
300  const ordinal_type matDim = 1, numFaces = 4, numOrts = 6;
301  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_TET_I1_FEM",
302  numFaces,
303  numOrts,
304  matDim,
305  matDim);
306  for (ordinal_type faceId=0;faceId<numFaces;++faceId)
307  init_TRI_FACE_ELEMENT_I1_FEM(matData, faceId);
308  }
309  else if (name == "Intrepid2_HDIV_WEDGE_I1_FEM") {
310  const ordinal_type matDim = 1, numFaces = 5, numOrts = 8;
311  matData = CoeffMatrixDataViewType("Orientation::CoeffMatrix::Intrepid2_HDIV_WEDGE_I1_FEM",
312  numFaces,
313  numOrts,
314  matDim,
315  matDim);
316  ordinal_type faceId = 0;
317  for ( ;faceId<3;++faceId)
318  init_QUAD_FACE_ELEMENT_I1_FEM(matData, faceId);
319  for ( ;faceId<numFaces;++faceId)
320  init_TRI_FACE_ELEMENT_I1_FEM(matData, faceId);
321  }
322  return matData;
323  }
324 
325  //
326  // Quad elements
327  //
328 
329  template<typename SpT>
330  void
331  OrientationTools<SpT>::
332  init_HGRAD_QUAD_Cn_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
333  const ordinal_type order) {
334  if (order > 1) {
335  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> lineBasis(order);
336  Basis_HGRAD_QUAD_Cn_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
337 
338  const ordinal_type numEdge = 4, numOrt = 2;
339  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
340  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
341  auto mat = Kokkos::subview(matData,
342  edgeId, edgeOrt,
343  Kokkos::ALL(), Kokkos::ALL());
345  lineBasis, cellBasis,
346  edgeId, edgeOrt);
347  }
348  }
349  }
350 
351  template<typename SpT>
352  void
353  OrientationTools<SpT>::
354  init_HCURL_QUAD_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
355  const ordinal_type order) {
356  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1, POINTTYPE_GAUSS);
357  Basis_HCURL_QUAD_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
358 
359  const ordinal_type numEdge = 4, numOrt = 2;
360  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
361  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
362  auto mat = Kokkos::subview(matData,
363  edgeId, edgeOrt,
364  Kokkos::ALL(), Kokkos::ALL());
366  bubbleBasis, cellBasis,
367  edgeId, edgeOrt);
368  }
369  }
370 
371  template<typename SpT>
372  void
373  OrientationTools<SpT>::
374  init_HDIV_QUAD_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
375  const ordinal_type order) {
376  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1, POINTTYPE_GAUSS);
377  Basis_HDIV_QUAD_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
378 
379  const ordinal_type numEdge = 4, numOrt = 2;
380  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
381  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
382  auto mat = Kokkos::subview(matData,
383  edgeId, edgeOrt,
384  Kokkos::ALL(), Kokkos::ALL());
386  bubbleBasis, cellBasis,
387  edgeId, edgeOrt);
388  }
389  }
390 
391  //
392  // Hexahedral elements
393  //
394 
395  template<typename SpT>
396  void
397  OrientationTools<SpT>::
398  init_HGRAD_HEX_Cn_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
399  const ordinal_type order) {
400  if (order > 1) {
401  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> lineBasis(order);
402  Basis_HGRAD_QUAD_Cn_FEM<Kokkos::DefaultHostExecutionSpace> quadBasis(order);
403  Basis_HGRAD_HEX_Cn_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
404 
405  const ordinal_type numEdge = 12, numFace = 6;
406  {
407  const ordinal_type numOrt = 2;
408  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
409  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
410  auto mat = Kokkos::subview(matData,
411  edgeId, edgeOrt,
412  Kokkos::ALL(), Kokkos::ALL());
414  lineBasis, cellBasis,
415  edgeId, edgeOrt);
416  }
417  }
418  {
419  const ordinal_type numOrt = 8;
420  for (ordinal_type faceId=0;faceId<numFace;++faceId)
421  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
422  auto mat = Kokkos::subview(matData,
423  numEdge+faceId, faceOrt,
424  Kokkos::ALL(), Kokkos::ALL());
426  quadBasis, cellBasis,
427  faceId, faceOrt);
428  }
429  }
430  }
431  }
432 
433  template<typename SpT>
434  void
435  OrientationTools<SpT>::
436  init_HCURL_HEX_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
437  const ordinal_type order) {
438  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1, POINTTYPE_GAUSS);
439  Basis_HCURL_QUAD_In_FEM<Kokkos::DefaultHostExecutionSpace> quadBasis(order);
440  Basis_HCURL_HEX_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
441 
442  const ordinal_type numEdge = 12, numFace = 6;
443  {
444  const ordinal_type numOrt = 2;
445  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
446  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
447  auto mat = Kokkos::subview(matData,
448  edgeId, edgeOrt,
449  Kokkos::ALL(), Kokkos::ALL());
451  bubbleBasis, cellBasis,
452  edgeId, edgeOrt);
453  }
454  }
455  if (order > 1) {
456  const ordinal_type numOrt = 8;
457  for (ordinal_type faceId=0;faceId<numFace;++faceId)
458  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
459  auto mat = Kokkos::subview(matData,
460  numEdge+faceId, faceOrt,
461  Kokkos::ALL(), Kokkos::ALL());
463  quadBasis, cellBasis,
464  faceId, faceOrt);
465  }
466  }
467  }
468 
469  template<typename SpT>
470  void
471  OrientationTools<SpT>::
472  init_HDIV_HEX_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
473  const ordinal_type order) {
474  Basis_HGRAD_QUAD_Cn_FEM<Kokkos::DefaultHostExecutionSpace> quadBasis(order-1, POINTTYPE_GAUSS);
475  Basis_HDIV_HEX_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
476 
477  const ordinal_type numFace = 6;
478  {
479  const ordinal_type numOrt = 8;
480  for (ordinal_type faceId=0;faceId<numFace;++faceId)
481  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
482  auto mat = Kokkos::subview(matData,
483  faceId, faceOrt,
484  Kokkos::ALL(), Kokkos::ALL());
486  quadBasis, cellBasis,
487  faceId, faceOrt);
488  }
489  }
490  }
491 
492  //
493  // Triangle elements
494  //
495 
496  template<typename SpT>
497  void
498  OrientationTools<SpT>::
499  init_HGRAD_TRI_Cn_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
500  const ordinal_type order) {
501  if (order > 1) {
502  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> lineBasis(order);
503  Basis_HGRAD_TRI_Cn_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
504 
505  const ordinal_type numEdge = 3, numOrt = 2;
506  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
507  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
508  auto mat = Kokkos::subview(matData,
509  edgeId, edgeOrt,
510  Kokkos::ALL(), Kokkos::ALL());
512  lineBasis, cellBasis,
513  edgeId, edgeOrt);
514  }
515  }
516  }
517 
518 
519  template<typename SpT>
520  void
521  OrientationTools<SpT>::
522  init_HCURL_TRI_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
523  const ordinal_type order) {
524  Basis_HVOL_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1);
525  Basis_HCURL_TRI_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
526 
527  const ordinal_type numEdge = 3, numOrt = 2;
528  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
529  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
530  auto mat = Kokkos::subview(matData,
531  edgeId, edgeOrt,
532  Kokkos::ALL(), Kokkos::ALL());
534  bubbleBasis, cellBasis,
535  edgeId, edgeOrt);
536  }
537  }
538 
539  template<typename SpT>
540  void
541  OrientationTools<SpT>::
542  init_HDIV_TRI_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
543  const ordinal_type order) {
544  Basis_HVOL_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1);
545  Basis_HDIV_TRI_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
546 
547  const ordinal_type numEdge = 3, numOrt = 2;
548  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
549  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
550  auto mat = Kokkos::subview(matData,
551  edgeId, edgeOrt,
552  Kokkos::ALL(), Kokkos::ALL());
554  bubbleBasis, cellBasis,
555  edgeId, edgeOrt);
556  }
557  }
558 
559  //
560  // Tetrahedral elements
561  //
562 
563  template<typename SpT>
564  void
565  OrientationTools<SpT>::
566  init_HGRAD_TET_Cn_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
567  const ordinal_type order) {
568  if (order > 1) {
569  Basis_HGRAD_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> lineBasis(order);
570  Basis_HGRAD_TRI_Cn_FEM<Kokkos::DefaultHostExecutionSpace> triBasis(order);
571  Basis_HGRAD_TET_Cn_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
572 
573  const ordinal_type numEdge = 6, numFace = 4;
574  {
575  const ordinal_type numOrt = 2;
576  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
577  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
578  auto mat = Kokkos::subview(matData,
579  edgeId, edgeOrt,
580  Kokkos::ALL(), Kokkos::ALL());
582  lineBasis, cellBasis,
583  edgeId, edgeOrt);
584  }
585  }
586  if (order > 2) {
587  const ordinal_type numOrt = 6;
588  for (ordinal_type faceId=0;faceId<numFace;++faceId)
589  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
590  auto mat = Kokkos::subview(matData,
591  numEdge+faceId, faceOrt,
592  Kokkos::ALL(), Kokkos::ALL());
594  triBasis, cellBasis,
595  faceId, faceOrt);
596  }
597  }
598  }
599  }
600 
601  template<typename SpT>
602  void
603  OrientationTools<SpT>::
604  init_HCURL_TET_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
605  const ordinal_type order) {
606  Basis_HVOL_LINE_Cn_FEM<Kokkos::DefaultHostExecutionSpace> bubbleBasis(order-1);
607  Basis_HCURL_TRI_In_FEM<Kokkos::DefaultHostExecutionSpace> triBasis(order);
608  Basis_HCURL_TET_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
609 
610  const ordinal_type numEdge = 6, numFace = 4;
611  {
612  const ordinal_type numOrt = 2;
613  for (ordinal_type edgeId=0;edgeId<numEdge;++edgeId)
614  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
615  auto mat = Kokkos::subview(matData,
616  edgeId, edgeOrt,
617  Kokkos::ALL(), Kokkos::ALL());
619  bubbleBasis, cellBasis,
620  edgeId, edgeOrt);
621  }
622  }
623  if (order > 1) {
624  const ordinal_type numOrt = 6;
625  for (ordinal_type faceId=0;faceId<numFace;++faceId)
626  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
627  auto mat = Kokkos::subview(matData,
628  numEdge+faceId, faceOrt,
629  Kokkos::ALL(), Kokkos::ALL());
631  triBasis, cellBasis,
632  faceId, faceOrt);
633  }
634  }
635  }
636 
637  template<typename SpT>
638  void
639  OrientationTools<SpT>::
640  init_HDIV_TET_In_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
641  const ordinal_type order) {
642  Basis_HVOL_TRI_Cn_FEM<Kokkos::DefaultHostExecutionSpace> triBasis(order-1);
643  Basis_HDIV_TET_In_FEM<Kokkos::DefaultHostExecutionSpace> cellBasis(order);
644 
645  const ordinal_type numFace = 4, numOrt = 6;
646  for (ordinal_type faceId=0;faceId<numFace;++faceId)
647  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
648  auto mat = Kokkos::subview(matData,
649  faceId, faceOrt,
650  Kokkos::ALL(), Kokkos::ALL());
652  triBasis, cellBasis,
653  faceId, faceOrt);
654  }
655  }
656 
657  //
658  // Lower order I1 elements
659  //
660 
661  template<typename SpT>
662  void
663  OrientationTools<SpT>::
664  init_EDGE_ELEMENT_I1_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
665  const ordinal_type edgeId) {
666  const ordinal_type numOrt = 2;
667  const double edgeOrtCoeff[2] = { 1.0, -1.0 };
668  for (ordinal_type edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
669  auto mat = Kokkos::subview(matData,
670  edgeId, edgeOrt,
671  Kokkos::ALL(), Kokkos::ALL());
672  mat(0,0) = edgeOrtCoeff[edgeOrt];
673  }
674  }
675 
676  template<typename SpT>
677  void
678  OrientationTools<SpT>::
679  init_TRI_FACE_ELEMENT_I1_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
680  const ordinal_type faceId) {
681  const ordinal_type numOrt = 6;
682  const double faceOrtCoeff[6] = { 1.0, 1.0, 1.0,
683  -1.0, -1.0, -1.0 };
684 
685  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
686  auto mat = Kokkos::subview(matData,
687  faceId, faceOrt,
688  Kokkos::ALL(), Kokkos::ALL());
689  mat(0,0) = faceOrtCoeff[faceOrt];
690  }
691  }
692 
693  template<typename SpT>
694  void
695  OrientationTools<SpT>::
696  init_QUAD_FACE_ELEMENT_I1_FEM(typename OrientationTools<SpT>::CoeffMatrixDataViewType matData,
697  const ordinal_type faceId) {
698  const ordinal_type numOrt = 8;
699  const double faceOrtCoeff[8] = { 1.0, 1.0, 1.0, 1.0,
700  -1.0, -1.0, -1.0, -1.0 };
701 
702  for (ordinal_type faceOrt=0;faceOrt<numOrt;++faceOrt) {
703  auto mat = Kokkos::subview(matData,
704  faceId, faceOrt,
705  Kokkos::ALL(), Kokkos::ALL());
706  mat(0,0) = faceOrtCoeff[faceOrt];
707  }
708  }
709 
710  template<typename SpT>
711  template<typename BasisPtrType>
712  typename OrientationTools<SpT>::CoeffMatrixDataViewType
713  OrientationTools<SpT>::createCoeffMatrix(BasisPtrType basis) {
714 #ifdef HAVE_INTREPID2_DEBUG
715  INTREPID2_TEST_FOR_EXCEPTION( !basis->requireOrientation(), std::invalid_argument,
716  ">>> ERROR (OrientationTools::createCoeffMatrix): basis does not require orientations." );
717 #endif
718  Kokkos::push_finalize_hook( [=] {
719  ortCoeffData=std::map<std::pair<std::string,ordinal_type>, typename OrientationTools<SpT>::CoeffMatrixDataViewType>();
720  });
721 
722  const std::pair<std::string,ordinal_type> key(basis->getName(), basis->getDegree());
723  const auto found = ortCoeffData.find(key);
724 
725  CoeffMatrixDataViewType matData;
726  if (found == ortCoeffData.end()) {
727  matData = createCoeffMatrixInternal(basis);
728  ortCoeffData.insert(std::make_pair(key, matData));
729  } else {
730  matData = found->second;
731  }
732 
733  return matData;
734  }
735 
736  template<typename SpT>
738  ortCoeffData.clear();
739  }
740 }
741 
742 #endif
743 
744 
745 // template<typename SpT>
746 // void
747 // OrientationTools<SpT>::CoeffMatrix::import(const OrientationTools<SpT>::DenseMatrix &b,
748 // const bool transpose) {
749 // #ifdef HAVE_INTREPID2_DEBUG
750 // INTREPID2_TEST_FOR_ABORT( !( NumRows() == b.NumRows() && NumCols() == b.NumCols() ), std::invalid_argument,
751 // ">>> ERROR (Intrepid::Orientation::CoeffMatrix::import): "
752 // "Matrix dimensions are not matched");
753 // #endif
754 // // count size
755 // const SpT eps = 1.0e-8;
756 // const ordinal_type nrows = (transpose ? b.NumCols() : b.NumRows());
757 // const ordinal_type ncols = (transpose ? b.NumRows() : b.NumCols());
758 // size_type nnz = b.countNumNonZeros(eps);
759 // createInternalArrays(nrows, ncols, nnz);
760 
761 // // construct sparse array
762 // nnz = 0;
763 // for (ordinal_type i=0;i<nrows;++i) {
764 // _ap(i) = nnz;
765 // for (ordinal_type j=0;j<ncols;++j) {
766 // const SpT val = (transpose ? b.Value(j,i) : b.Value(i,j));
767 // const SpT val2 = val*val;
768 
769 // // consider it as a nonzero entry
770 // if (val2 > eps) {
771 // _aj(nnz) = j;
772 // _ax(nnz) = val;
773 // ++nnz;
774 // }
775 // }
776 // }
777 // _ap(nrows) = nnz;
778 // }
779 
780 // template<typename SpT>
781 // std::ostream&
782 // OrientationTools<SpT>::CoeffMatrix::showMe(std::ostream &os) const {
783 // std::ofstream prec;
784 // prec.copyfmt(os);
785 
786 // os.precision(3);
787 
788 // os << " -- OrientationTools::CoeffMatrix -- " << std::endl
789 // << " # of Rows = " << _m << std::endl
790 // << " # of Cols = " << _n << std::endl
791 // << std::endl
792 // << " RowPtrArray length = " << _ap.extent(0) << std::endl
793 // << " ColArray length = " << _aj.extent(0) << std::endl
794 // << " ValueArray length = " << _ax.extent(0) << std::endl
795 // << std::endl;
796 
797 // const ordinal_type w = 10;
798 // if (_ap.size() && _aj.size() && _ax.size()) {
799 // os << std::setw(w) << "Row" << " "
800 // << std::setw(w) << "Col" << " "
801 // << std::setw(w) << "Val" << std::endl;
802 // for (ordinal_type i=0;i<_m;++i) {
803 // size_type jbegin = _ap[i], jend = _ap[i+1];
804 // for (ordinal_type j=jbegin;j<jend;++j) {
805 // SpT val = _ax[j];
806 // os << std::setw(w) << i << " "
807 // << std::setw(w) << _aj[j] << " "
808 // << std::setw(w) << val << std::endl;
809 // }
810 // }
811 // }
812 // os.copyfmt(prec);
813 
814 // return os;
815 // }
816 
817 // template<class Scalar>
818 // size_type
819 // OrientationTools<Scalar>::DenseMatrix::countNumNonZeros(const Scalar epsilon) const {
820 // size_type nnz = 0;
821 // for (ordinal_type j=0;j<NumCols();++j) {
822 // for (ordinal_type i=0;i<NumRows();++i) {
823 // const Scalar val = Value(i,j);
824 // nnz += ((val*val) > epsilon);
825 // }
826 // }
827 // return nnz;
828 // }
829 
830 // template<class Scalar>
831 // std::ostream&
832 // OrientationTools<Scalar>::DenseMatrix::showMe(std::ostream &os) const {
833 // std::ofstream prec;
834 // prec.copyfmt(os);
835 
836 // os.precision(3);
837 
838 // os << " -- OrientationTools::DenseMatrix -- " << std::endl
839 // << " # of Rows = " << _m << std::endl
840 // << " # of Cols = " << _n << std::endl
841 // << " Col Stride = " << _cs << std::endl
842 // << " Row Stride = " << _rs << std::endl
843 // << std::endl
844 // << " ValueArray dimensions = " << _a.extent(0) << std::endl
845 // << std::endl;
846 
847 // const ordinal_type w = 10;
848 // if (_a.size()) {
849 // for (ordinal_type i=0;i<_m;++i) {
850 // for (ordinal_type j=0;j<_n;++j) {
851 // const Scalar val = this->Value(i,j);
852 // os << std::setw(w) << val << " ";
853 // }
854 // os << std::endl;
855 // }
856 // }
857 // os.copyfmt(prec);
858 
859 // return os;
860 // }
861 
862 
863 
864 
865 
866 // template<typename SpT>
867 // void
868 // OrientationTools<SpT>::
869 // initHexahedron(Kokkos::View<CoeffMatrixType***,SpT> MatrixData,
870 // const EFunctionSpace space,
871 // const ordinal_type order) {
872 
873 // switch (space) {
874 // case FUNCTION_SPACE_HGRAD: {
875 // Basis_HGRAD_LINE_Cn_FEM<SpT> lineBasis(order);
876 // Basis_HGRAD_QUAD_Cn_FEM<SpT> faceBasis(order);
877 // Basis_HGRAD_HEXA_Cn_FEM<SpT> cellBasis(order);
878 
879 // {
880 // const ordinal_type numEdge = 12, numOrt = 2;
881 // for (auto edgeId=0;edgeId<numEdge;++edgeId)
882 // for (auto edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
883 // const auto C = Impl::OrientationTools::getEdgeCoeffMatrix_HGRAD(lineBasis, cellBasis, edgeId, edgeOrt);
884 // MatrixData(0, edgeId, edgeOrt) = Kokkos::create_mirror_view(typename SpT::memory_space(), C);
885 // Kokkos::deep_copy(MatrixData(0, edgdId, edgeOrt), C);
886 // }
887 // }
888 // {
889 // const ordinal_type numFace = 6, numOrt = 8;
890 // for (auto faceId=0;faceId<numFace;++faceId)
891 // for (auto faceOrt=0;faceOrt<numOrt;++faceOrt) {
892 // const auto C = Impl::OrientationTools::getQuadrilateralCoeffMatrix_HGRAD(faceBasis, cellBasis, faceId, faceOrt);
893 // MatrixData(1, faceId, faceOrt) = Kokkos::create_mirror_view(typename SpT::memory_space(), C);
894 // Kokkos::deep_copy(MatrixData(1, faceId, faceOrt), C);
895 // }
896 // }
897 // break;
898 // }
899 // default: {
900 // INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
901 // ">>> ERROR (Intrepid::OrientationTools::initQuadrilateral): Invalid function space.");
902 // break;
903 // }
904 // }
905 // }
906 
907 // template<typename SpT>
908 // void
909 // OrientationTools<SpT>::
910 // initTetrahedron(Kokkos::View<CoeffMatrixType***,SpT> MatrixData,
911 // const EFunctionSpace space,
912 // const ordinal_type order) {
913 
914 // switch (space) {
915 // case FUNCTION_SPACE_HGRAD: {
916 // Basis_HGRAD_LINE_Cn_FEM<SpT> lineBasis(order);
917 // Basis_HGRAD_TRI_Cn_FEM <SpT> faceBasis(order);
918 // Basis_HGRAD_TET_Cn_FEM <SpT> cellBasis(order);
919 
920 // {
921 // const ordinal_type numEdge = 6, numOrt = 2;
922 // for (auto edgeId=0;edgeId<numEdge;++edgeId)
923 // for (auto edgeOrt=0;edgeOrt<numOrt;++edgeOrt) {
924 // const auto C = Impl::OrientationTools::getEdgeCoeffMatrix_HGRAD(lineBasis, cellBasis, edgeId, edgeOrt);
925 // MatrixData(0, edgeId, edgeOrt) = Kokkos::create_mirror_view(typename SpT::memory_space(), C);
926 // Kokkos::deep_copy(MatrixData(0, edgdId, edgeOrt), C);
927 // }
928 // }
929 // {
930 // const ordinal_type numFace = 4, numOrt = 6;
931 // for (auto faceId=0;faceId<numFace;++faceId)
932 // for (auto faceOrt=0;faceOrt<numOrt;++faceOrt) {
933 // const auto C = Impl::OrientationTools::getTriangleCoeffMatrix_HGRAD(faceBasis, cellBasis, faceId, faceOrt);
934 // MatrixData(1, faceId, faceOrt) = Kokkos::create_mirror_view(typename SpT::memory_space(), C);
935 // Kokkos::deep_copy(MatrixData(1, faceId, faceOrt), C);
936 // }
937 // }
938 // break;
939 // }
940 // default: {
941 // INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
942 // ">>> ERROR (Intrepid::OrientationTools::initQuadrilateral): Invalid function space.");
943 // break;
944 // }
945 // }
946 // }
static void getCoeffMatrix_HGRAD(outputViewType &output, const subcellBasisType subcellBasis, const cellBasisType cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute coefficient matrix for HGRAD by collocating point values.
static void getCoeffMatrix_HCURL(outputViewType &output, const subcellBasisType subcellBasis, const cellBasisType cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute coefficient matrix for HCURL by collocating point values.
static void getCoeffMatrix_HDIV(outputViewType &output, const subcellBasisType subcellBasis, const cellBasisType cellBasis, const ordinal_type subcellId, const ordinal_type subcellOrt)
Compute coefficient matrix for HDIV by collocating point values.
static void clearCoeffMatrix()
Clear coefficient matrix.
static CoeffMatrixDataViewType createCoeffMatrix(BasisPtrType basis)
Create coefficient matrix.