Intrepid
Intrepid_FunctionSpaceToolsInPlaceDef.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 
52  template<class Scalar, class ArrayType>
54  {
55  return;
56  }
57 
58  template<class Scalar, class ArrayType>
60  {
61  return;
62  }
63 
64  template<class Scalar, class ArrayType, class ArrayTypeJac>
66  const ArrayTypeJac & jacobianInverse,
67  const char transpose)
68  {
69  FieldContainer<Scalar> tmp(inOutVals.dimension(3));
70 
71  // test for transpose direction, one loop nest for each direction
72  if (transpose == 'T')
73  {
74  for (int c=0;c<inOutVals.dimension(0);c++)
75  {
76  for (int f=0;f<inOutVals.dimension(1);f++)
77  {
78  for (int p=0;p<inOutVals.dimension(2);p++)
79  {
80  for (int d1=0;d1<inOutVals.dimension(3);d1++)
81  {
82  tmp(d1) = 0.0;
83  for (int d2=0;d2<inOutVals.dimension(3);d2++)
84  {
85  tmp(d1) += jacobianInverse(c,p,d2,d1) * inOutVals(c,f,p,d2);
86  }
87  }
88  for (int d1=0;d1<inOutVals.dimension(3);d1++)
89  {
90  inOutVals(c,f,p,d1) = tmp(d1);
91  }
92  }
93  }
94  }
95  }
96  else if (transpose == 'N')
97  {
98  for (int c=0;c<inOutVals.dimension(0);c++)
99  {
100  for (int f=0;f<inOutVals.dimension(1);f++)
101  {
102  for (int p=0;p<inOutVals.dimension(2);p++)
103  {
104  for (int d1=0;d1<inOutVals.dimension(3);d1++)
105  {
106  tmp(d1) = 0.0;
107  for (int d2=0;d2<inOutVals.dimension(3);d2++)
108  {
109  tmp(d1) += jacobianInverse(c,p,d1,d2) * inOutVals(c,f,p,d2);
110  }
111  }
112  for (int d1=0;d1<inOutVals.dimension(3);d1++)
113  {
114  inOutVals(c,f,p,d1) = tmp(d1);
115  }
116  }
117  }
118  }
119  }
120  else
121  {
122  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
123  "Intrepid:FunctionSpaceToolsInPlace::HGRADtransformGRAD::Unknown transpose type" );
124  }
125  }
126 
127  template<class Scalar, class ArrayType, class ArrayTypeJac>
129  const ArrayTypeJac & jacobianInverse,
130  const char transpose)
131  {
132  char t_loc;
133  if (transpose == 'T')
134  {
135  t_loc = 'N';
136  }
137  else if (transpose == 'N')
138  {
139  t_loc = 'T';
140  }
141  else
142  {
143  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
144  "Intrepid:FunctionSpaceToolsInPlace::HGRADtransformGRADDual::Unknown transpose type" );
145  }
146  FunctionSpaceToolsInPlace::HGRADtransformGRAD<Scalar,ArrayType,ArrayTypeJac>( inOutVals,
147  jacobianInverse,
148  t_loc);
149  }
150 
151 
152  template<class Scalar, class ArrayType, class ArrayTypeJac>
154  const ArrayTypeJac & jacobianInverse,
155  const char transpose)
156  {
157  FunctionSpaceToolsInPlace::HGRADtransformGRAD<Scalar,ArrayType,ArrayTypeJac>( inOutVals , jacobianInverse, transpose );
158  }
159 
160  template<class Scalar, class ArrayType, class ArrayTypeJac>
162  const ArrayTypeJac & jacobianInverse,
163  const char transpose)
164  {
165  char t_loc;
166  if (transpose == 'T')
167  {
168  t_loc = 'N';
169  }
170  else if (transpose == 'N')
171  {
172  t_loc = 'T';
173  }
174  else
175  {
176  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
177  "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformVALUEDual::Unknown transpose type" );
178  }
179  FunctionSpaceToolsInPlace::HCURLtransformVALUEDual<Scalar,ArrayType,ArrayTypeJac>(inOutVals,
180  jacobianInverse,
181  t_loc);
182 
183  }
184 
185  template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet>
187  const ArrayTypeJac & jacobian,
188  const ArrayTypeDet & jacobianDet,
189  const char transpose)
190  {
191  FieldContainer<Scalar> tmp(inOutVals.dimension(3));
192 
193  // test for transpose direction, one loop nest for each direction
194  if (transpose == 'T')
195  {
196  for (int c=0;c<inOutVals.dimension(0);c++)
197  {
198  for (int f=0;f<inOutVals.dimension(1);f++)
199  {
200  for (int p=0;p<inOutVals.dimension(2);p++)
201  {
202  for (int d1=0;d1<inOutVals.dimension(3);d1++)
203  {
204  tmp(d1) = 0.0;
205  for (int d2=0;d2<inOutVals.dimension(3);d2++)
206  {
207  tmp(d1) += jacobian(c,p,d2,d1) * inOutVals(c,f,p,d2);
208  }
209  }
210  for (int d1=0;d1<inOutVals.dimension(3);d1++)
211  {
212  inOutVals(c,f,p,d1) = tmp(d1) / jacobianDet(c,p);
213  }
214  }
215  }
216  }
217  }
218  else if (transpose == 'N')
219  {
220  for (int c=0;c<inOutVals.dimension(0);c++)
221  {
222  for (int f=0;f<inOutVals.dimension(1);f++)
223  {
224  for (int p=0;p<inOutVals.dimension(2);p++)
225  {
226  for (int d1=0;d1<inOutVals.dimension(3);d1++)
227  {
228  tmp(d1) = 0.0;
229  for (int d2=0;d2<inOutVals.dimension(3);d2++)
230  {
231  tmp(d1) += jacobian(c,p,d1,d2) * inOutVals(c,f,p,d2);
232  }
233  }
234  for (int d1=0;d1<inOutVals.dimension(3);d1++)
235  {
236  inOutVals(c,f,p,d1) = tmp(d1) / jacobianDet(c,p);
237  }
238  }
239  }
240  }
241  }
242  else
243  {
244  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
245  "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformCURL::Unknown transpose type" );
246  }
247 
248  }
249 
250  template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet>
252  const ArrayTypeJac & jacobian,
253  const ArrayTypeDet & jacobianDet,
254  const char transpose)
255  {
256  char t_loc;
257  if (transpose == 'T')
258  {
259  t_loc = 'N';
260  }
261  else if (transpose == 'N')
262  {
263  t_loc = 'T';
264  }
265  else
266  {
267  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
268  "Intrepid:FunctionSpaceToolsInPlace::HCURLtransformCURLDual::Unknown transpose type" );
269  }
270  FunctionSpaceToolsInPlace::HCURLtransformCURLDual<Scalar,ArrayType,ArrayTypeJac>(inOutVals,
271  jacobian,
272  jacobianDet,
273  t_loc);
274  }
275 
276  template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet>
278  const ArrayTypeJac & jacobian,
279  const ArrayTypeDet & jacobianDet,
280  const char transpose)
281  {
282  FunctionSpaceToolsInPlace::HCURLtransformCURL<Scalar,ArrayType,ArrayTypeJac,ArrayTypeDet>( inOutVals,jacobian,jacobianDet,transpose);
283  }
284 
285  template<class Scalar, class ArrayType, class ArrayTypeJac, class ArrayTypeDet>
287  const ArrayTypeJac & jacobian,
288  const ArrayTypeDet & jacobianDet,
289  const char transpose)
290  {
291  char t_loc;
292  if (transpose == 'T')
293  {
294  t_loc = 'N';
295  }
296  else if (transpose == 'N')
297  {
298  t_loc = 'T';
299  }
300  else
301  {
302  TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument ,
303  "FunctionSpaceToolsInPlace::HDIVtransformVALUEDual: invalid transpose character");
304  }
305  FunctionSpaceToolsInPlace::HDIVtransformVALUE<Scalar,ArrayType,ArrayTypeJac,ArrayTypeDet>( inOutVals,
306  jacobian,
307  jacobianDet ,
308  t_loc );
309  }
310 
311  template<class Scalar, class ArrayType, class ArrayTypeDet>
313  const ArrayTypeDet & jacobianDet)
314  {
315  for (int c=0;c<inOutVals.dimension(0);c++)
316  {
317  for (int f=0;f<inOutVals.dimension(1);f++)
318  {
319  for (int p=0;p<inOutVals.dimension(2);p++)
320  {
321  inOutVals(c,f,p) /= jacobianDet(c,p);
322  }
323  }
324  }
325  }
326 
327  template<class Scalar, class ArrayType, class ArrayTypeDet>
329  const ArrayTypeDet & jacobianDet)
330  {
331  FunctionSpaceToolsInPlace::HDIVtransformDIV<Scalar,ArrayType,ArrayTypeDet>( inOutVals ,
332  jacobianDet );
333  }
334 
335  template<class Scalar, class ArrayType, class ArrayTypeDet>
337  const ArrayTypeDet & jacobianDet)
338  {
339  FunctionSpaceToolsInPlace::HDIVtransformDIV<Scalar,ArrayType,ArrayTypeDet>( inOutVals,jacobianDet);
340  }
341 
342  template<class Scalar, class ArrayType, class ArrayTypeDet>
344  const ArrayTypeDet & jacobianDet)
345  {
346  FunctionSpaceToolsInPlace::HVOLtransformVALUEDual( inOutVals ,jacobianDet );
347  }
348 
349 
350 
351  template<class Scalar, class ArrayType, class ArrayTypeMeasure>
352  void FunctionSpaceToolsInPlace::multiplyMeasure(ArrayType & inOutVals,
353  const ArrayTypeMeasure & inMeasure)
354  {
355  if (inOutVals.rank() == 2) // inOutVals is (C,P)
356  {
357  for (int c=0;c<inOutVals.dimension(0);c++)
358  {
359  for (int p=0;p<inOutVals.dimension(0);p++)
360  {
361  inOutVals(c,p) *= inMeasure(c,p);
362  }
363  }
364  }
365  else if (inOutVals.rank() == 3) // inOutVals is (C,F,P)
366  {
367  for (int c=0;c<inOutVals.dimension(0);c++)
368  {
369  for (int f=0;f<inOutVals.dimension(1);f++)
370  {
371  for (int p=0;p<inOutVals.dimension(2);p++)
372  {
373  inOutVals(c,f,p) *= inMeasure(c,p);
374  }
375  }
376  }
377  }
378  else if (inOutVals.rank() == 4) // inOutVals is (C,F,P)
379  {
380  for (int c=0;c<inOutVals.dimension(0);c++)
381  {
382  for (int f=0;f<inOutVals.dimension(1);f++)
383  {
384  for (int p=0;p<inOutVals.dimension(2);p++)
385  {
386  for (int d=0;d<inOutVals.dimension(3);d++)
387  {
388  inOutVals(c,f,p,d) *= inMeasure(c,p);
389  }
390  }
391  }
392  }
393  }
394  }
395 
396 } // end namespace Intrepid
static void HCURLtransformCURLDual(ArrayType &outVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const char transpose= 'N')
Applies the dual of the HCURLtransformCURL transformation.
static void HVOLtransformVALUEDual(ArrayType &inOutVals, const ArrayTypeDet &jacobianDet)
Applies the dual of HVOLtransformVALUE.
static void HCURLtransformVALUE(ArrayType &inOutVals, const ArrayTypeJac &jacobianInverse, const char transpose= 'T')
Transformation of a (vector) value field in the H-curl space, defined at points on a reference cell...
static void HCURLtransformVALUEDual(ArrayType &outVals, const ArrayTypeJac &jacobianInverse, const char transpose= 'T')
Applies the dual of the HCURLtransformVALUE transformation.
static void HDIVtransformVALUEDual(ArrayType &outVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const char transpose= 'N')
Applies the dual of HDIVtransformVALUE.
static void HGRADtransformGRAD(ArrayType &inOutVals, const ArrayTypeJac &jacobianInverse, const char transpose= 'T')
Transformation of a gradient field in the H-grad space, defined at points on a reference cell...
static void HDIVtransformDIV(ArrayType &inOutVals, const ArrayTypeDet &jacobianDet)
Transformation of a divergence field in the H-div space, defined at points on a reference cell...
static void HCURLtransformCURL(ArrayType &inOutVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const char transpose= 'N')
Transformation of a curl field in the H-curl space, defined at points on a reference cell...
static void HGRADtransformVALUE(ArrayType &inOutVals)
Transformation of a (scalar) value field in the H-grad space, defined at points on a reference cell...
static void HVOLtransformVALUE(ArrayType &inOutVals, const ArrayTypeDet &jacobianDet)
Transformation of a (scalar) value field in the H-vol space, defined at points on a reference cell...
static void HDIVtransformDIVDual(ArrayType &inOutVals, const ArrayTypeDet &jacobianDet)
Applies the dual of HDIVtransformDIV, which is the same.
static void HGRADtransformGRADDual(ArrayType &inOutVals, const ArrayTypeJac &jacobianInverse, const char transpose= 'T')
Applies the transpose of the HGRADtransformGRAD to the data.
static void HDIVtransformVALUE(ArrayType &inOutVals, const ArrayTypeJac &jacobian, const ArrayTypeDet &jacobianDet, const char transpose= 'N')
Transformation of a (vector) value field in the H-div space, defined at points on a reference cell...
static void HGRADtransformVALUEDual(ArrayType &inOutVals)
Since there is no matrix involved, this is the same transformation as HGRADtransformVALUE.