Intrepid
test_01.cpp
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 #include "Teuchos_oblackholestream.hpp"
51 #include "Teuchos_RCP.hpp"
52 #include "Teuchos_GlobalMPISession.hpp"
53 
54 #include <random>
55 
56 using namespace std;
57 using namespace Intrepid;
58 
59 #define INTREPID_TEST_COMMAND( S , throwCounter, nException ) \
60 { \
61  ++nException; \
62  try { \
63  S ; \
64  } \
65  catch (const std::logic_error & err) { \
66  ++throwCounter; \
67  *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
68  *outStream << err.what() << '\n'; \
69  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
70  }; \
71 }
72 
73 int main(int argc, char *argv[]) {
74 
75  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
76 
77  // This little trick lets us print to std::cout only if
78  // a (dummy) command-line argument is provided.
79  int iprint = argc - 1;
80  Teuchos::RCP<std::ostream> outStream;
81  Teuchos::oblackholestream bhs; // outputs nothing
82  if (iprint > 0)
83  outStream = Teuchos::rcp(&std::cout, false);
84  else
85  outStream = Teuchos::rcp(&bhs, false);
86 
87  // Save the format state of the original std::cout.
88  Teuchos::oblackholestream oldFormatState;
89  oldFormatState.copyfmt(std::cout);
90 
91  *outStream \
92  << "===============================================================================\n" \
93  << "| |\n" \
94  << "| Unit Test (Basis_HGRAD_TET_COMP12_FEM) |\n" \
95  << "| |\n" \
96  << "| 1) Evaluation of Basis Function Values |\n" \
97  << "| |\n" \
98  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \
99  << "| Denis Ridzal (dridzal@sandia.gov), |\n" \
100  << "| Kara Peterson (kjpeter@sandia.gov), |\n" \
101  << "| Jake Ostien (jtostie@sandia.gov). |\n" \
102  << "| |\n" \
103  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
104  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
105  << "| |\n" \
106  << "===============================================================================\n"\
107  << "| TEST 1: Basis creation, values |\n"\
108  << "===============================================================================\n";
109 
110  // Define basis and error flag
112  int errorFlag = 0;
113 
114  // Initialize throw counter for exception testing
115  //int nException = 0;
116  //int throwCounter = 0;
117 
118  // Define array containing the 10 vertices of the reference TET
119  FieldContainer<double> tetNodes(10, 3);
120  tetNodes(0,0) = 0.0; tetNodes(0,1) = 0.0; tetNodes(0,2) = 0.0;
121  tetNodes(1,0) = 1.0; tetNodes(1,1) = 0.0; tetNodes(1,2) = 0.0;
122  tetNodes(2,0) = 0.0; tetNodes(2,1) = 1.0; tetNodes(2,2) = 0.0;
123  tetNodes(3,0) = 0.0; tetNodes(3,1) = 0.0; tetNodes(3,2) = 1.0;
124  tetNodes(4,0) = 0.5; tetNodes(4,1) = 0.0; tetNodes(4,2) = 0.0;
125  tetNodes(5,0) = 0.5; tetNodes(5,1) = 0.5; tetNodes(5,2) = 0.0;
126  tetNodes(6,0) = 0.0; tetNodes(6,1) = 0.5; tetNodes(6,2) = 0.0;
127  tetNodes(7,0) = 0.0; tetNodes(7,1) = 0.0; tetNodes(7,2) = 0.5;
128  tetNodes(8,0) = 0.5; tetNodes(8,1) = 0.0; tetNodes(8,2) = 0.5;
129  tetNodes(9,0) = 0.0; tetNodes(9,1) = 0.5; tetNodes(9,2) = 0.5;
130  // Define array containing 5 integration points
131  FieldContainer<double> tetPoints(9, 3);
132  // from the 5 point integration
133  tetPoints(0,0) = 0.25; tetPoints(0,1) = 0.25; tetPoints(0,2) = 0.25;
134  tetPoints(1,0) = 0.5; tetPoints(1,1) = (1./6.); tetPoints(1,2) = (1./6.);
135  tetPoints(2,0) = (1./6.); tetPoints(2,1) = 0.5; tetPoints(2,2) = (1./6.);
136  tetPoints(3,0) = (1./6.); tetPoints(3,1) = (1./6.); tetPoints(3,2) = 0.5;
137  tetPoints(4,0) = (1./6.); tetPoints(4,1) = (1./6.); tetPoints(4,2) = (1./6.);
138  // from the 4 point integration
139  tetPoints(5,0) = 0.1381966011250105151795413165634361882280;
140  tetPoints(5,1) = 0.1381966011250105151795413165634361882280;
141  tetPoints(5,2) = 0.1381966011250105151795413165634361882280;
142 
143  tetPoints(6,0) = 0.5854101966249684544613760503096914353161;
144  tetPoints(6,1) = 0.1381966011250105151795413165634361882280;
145  tetPoints(6,2) = 0.1381966011250105151795413165634361882280;
146 
147  tetPoints(7,0) = 0.1381966011250105151795413165634361882280;
148  tetPoints(7,1) = 0.5854101966249684544613760503096914353161;
149  tetPoints(7,2) = 0.1381966011250105151795413165634361882280;
150 
151  tetPoints(8,0) = 0.1381966011250105151795413165634361882280;
152  tetPoints(8,1) = 0.1381966011250105151795413165634361882280;
153  tetPoints(8,2) = 0.5854101966249684544613760503096914353161;
154 
155  // output precision
156  outStream -> precision(20);
157 
158  // VALUE: Each row gives the 10 correct basis set values at an evaluation point
159  double nodalBasisValues[] = {
160  // first 4 vertices
161  1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
162  0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
163  0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
164  0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
165  // second 6 vertices
166  0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
167  0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
168  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
169  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
170  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
171  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0
172  };
173  double pointBasisValues[] = {
174  // pt 0 {0.25, 0.25, 0.25}
175  0.0, 0.0, 0.0, 0.0, 1./6., 1./6., 1./6., 1./6., 1./6., 1./6.,
176  // pt 1 {0.5, 1/6, 1/6}
177  0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3., 0.0,
178  // pt 2 {1/6, 0.5, 0.1/6}
179  0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 0.0, 0.0, 1./3.,
180  // pt 3 {1/6, 1/6, 0.5}
181  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1./3., 1./3., 1./3.,
182  // pt 4 {1/6, 1/6, 1/6}
183  0.0, 0.0, 0.0, 0.0, 1./3., 0.0, 1./3., 1./3., 0.0, 0.0,
184  // pt 5
185  0.170820393249936908922752100619382870632, 0.0, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0,
186  // pt 6
187  0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.0,
188  // pt 7
189  0.0, 0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.0, 0.0, 0.276393202250021030359082633126872376456,
190  // pt 8
191  0.0, 0.0, 0.0, 0.170820393249936908922752100619382870632, 0.0, 0.0, 0.0, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456, 0.276393202250021030359082633126872376456,
192  };
193 
194  // GRAD and D1: each row gives the 3x10 correct values of the gradients of the 10 basis functions
195  double pointBasisGrads[] = {
196  // point 0
197  -1./4., -1./4., -1./4., \
198  1./4., 0.0, 0.0, \
199  0.0, 1./4., 0.0, \
200  0.0, 0.0, 1./4., \
201  0.0, -3./4., -3./4., \
202  3./4., 3./4., 0.0, \
203  -3./4., 0.0, -3./4., \
204  -3./4., -3./4., 0.0, \
205  3./4., 0.0, 3./4., \
206  0.0, 3./4., 3./4., \
207 
208  // point 1
209  -1./24., -1./24., -1./24., \
210  7./8., 0.0, 0.0, \
211  0.0, 1./24., 0.0, \
212  0.0, 0.0, 1./24., \
213  -35./36., -19./12., -19./12., \
214  11./18., 19./12., 0.0, \
215  -17./36., 0.0, -1./3., \
216  -17./36., -1./3., 0.0, \
217  11./18., 0.0, 19./12., \
218  -5./36., 1./3., 1./3., \
219 
220  // point 2
221  -1./24., -1./24., -1./24., \
222  1./24., 0.0, 0.0, \
223  0.0, 7./8., 0.0, \
224  0.0, 0.0, 1./24., \
225  0.0, -17./36., -1./3., \
226  19./12., 11./18., 0.0, \
227  -19./12., -35./36., -19./12., \
228  -1./3., -17./36., 0.0, \
229  1./3., -5./36., 1./3., \
230  0.0, 11./18., 19./12., \
231 
232  // point 3
233  -1./24., -1./24., -1./24., \
234  1./24., 0.0, 0.0, \
235  0.0, 1./24., 0.0, \
236  0.0, 0.0, 7./8., \
237  0.0, -1./3., -17./36., \
238  1./3., 1./3., -5./36., \
239  -1./3., 0.0, -17./36., \
240  -19./12., -19./12., -35./36., \
241  19./12., 0.0, 11./18., \
242  0.0, 19./12., 11./18., \
243 
244  // point 4
245  -7./8., -7./8., -7./8., \
246  1./24., 0.0, 0.0, \
247  0.0, 1./24., 0.0, \
248  0.0, 0.0, 1./24., \
249  35./36., -11./18., -11./18., \
250  17./36., 17./36., 5./36., \
251  -11./18., 35./36., -11./18., \
252  -11./18., -11./18., 35./36., \
253  17./36., 5./36., 17./36., \
254  5./36., 17./36., 17./36., \
255 
256  // point 5
257  -1.088525491562421136153440125774228588290, -1.088525491562421136153440125774228588290, -1.088525491562421136153440125774228588290, \
258  -0.029508497187473712051146708591409529430, 0.0, 0.0, \
259  0.0, -0.029508497187473712051146708591409529430, 0.0, \
260  0.0, 0.0, -0.029508497187473712051146708591409529430, \
261  1.30437298687487732290535130675991113734, -0.563661001875017525299235527605726980380, -0.563661001875017525299235527605726980380, \
262  0.377322003750035050598471055211453960760, 0.377322003750035050598471055211453960760, 0.186338998124982474700764472394273019620, \
263  -0.563661001875017525299235527605726980380, 1.30437298687487732290535130675991113734, -0.563661001875017525299235527605726980380, \
264  -0.563661001875017525299235527605726980380, -0.563661001875017525299235527605726980380, 1.30437298687487732290535130675991113734, \
265  0.377322003750035050598471055211453960760, 0.186338998124982474700764472394273019620, 0.377322003750035050598471055211453960760, \
266  0.186338998124982474700764472394273019620, 0.377322003750035050598471055211453960760, 0.377322003750035050598471055211453960760, \
267 
268  // point 6
269  0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \
270  1.088525491562421136153440125774228588290, 0.0, 0.0, \
271  0.0, -0.029508497187473712051146708591409529430, 0.0, \
272  0.0, 0.0, -0.029508497187473712051146708591409529430, \
273  -1.30437298687487732290535130675991113734, -1.868033988749894848204586834365638117720, -1.868033988749894848204586834365638117720, \
274  0.563661001875017525299235527605726980380, 1.868033988749894848204586834365638117720, 0.0, \
275  -0.377322003750035050598471055211453960760, 0.0, -0.190983005625052575897706582817180941140, \
276  -0.377322003750035050598471055211453960760, -0.190983005625052575897706582817180941140, 0.0, \
277  0.563661001875017525299235527605726980380, 0.0, 1.868033988749894848204586834365638117720, \
278  -0.186338998124982474700764472394273019620, 0.190983005625052575897706582817180941140, 0.19098300562505257589770658281718094114, \
279 
280  // point 7
281  0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \
282  -0.029508497187473712051146708591409529430, 0.0, 0.0, \
283  0.0, 1.088525491562421136153440125774228588290, 0.0, \
284  0.0, 0.0, -0.029508497187473712051146708591409529430, \
285  0.0, -0.377322003750035050598471055211453960760, -0.190983005625052575897706582817180941140, \
286  1.868033988749894848204586834365638117720, 0.563661001875017525299235527605726980380, 0.0, \
287  -1.868033988749894848204586834365638117720, -1.30437298687487732290535130675991113734, -1.868033988749894848204586834365638117720, \
288  -0.190983005625052575897706582817180941140, -0.377322003750035050598471055211453960760, 0.0, \
289  0.190983005625052575897706582817180941140, -0.186338998124982474700764472394273019620, 0.190983005625052575897706582817180941140, \
290  0.0, 0.563661001875017525299235527605726980380, 1.868033988749894848204586834365638117720, \
291 
292  // point 8
293  0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, 0.029508497187473712051146708591409529430, \
294  -0.029508497187473712051146708591409529430, 0.0, 0.0, \
295  0.0, -0.029508497187473712051146708591409529430, 0.0, \
296  0.0, 0.0, 1.088525491562421136153440125774228588290, \
297  0.0, -0.190983005625052575897706582817180941140, -0.377322003750035050598471055211453960760, \
298  0.190983005625052575897706582817180941140, 0.190983005625052575897706582817180941140, -0.186338998124982474700764472394273019620, \
299  -0.190983005625052575897706582817180941140, 0.0, -0.377322003750035050598471055211453960760, \
300  -1.868033988749894848204586834365638117720, -1.868033988749894848204586834365638117720, -1.30437298687487732290535130675991113734,
301  1.868033988749894848204586834365638117720, 0.0, 0.563661001875017525299235527605726980380, \
302  0.0, 1.868033988749894848204586834365638117720, 0.563661001875017525299235527605726980380, \
303  };
304 
305  try{
306 
307  // Dimensions for the output arrays:
308  int numFields = tetBasis.getCardinality();
309  int numNodes = tetNodes.dimension(0);
310  int spaceDim = tetBasis.getBaseCellTopology().getDimension();
311 
312  // Generic array for values, grads, curls, etc. that will be properly sized before each call
314 
315  // Check VALUE of basis functions at nodes: resize vals to rank-2 container:\n";
316  *outStream << " check VALUE of basis functions at nodes\n";
317  vals.resize(numFields, numNodes);
318  tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
319 
320  for (int i = 0; i < numFields; i++) {
321  for (int j = 0; j < numNodes; j++) {
322  int l = i + j * numFields;
323  if (std::abs(vals(i,j) - nodalBasisValues[l]) > INTREPID_TOL) {
324  errorFlag++;
325  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
326 
327  // Output the multi-index of the value where the error is:
328  *outStream << " At multi-index { ";
329  *outStream << i << " ";*outStream << j << " ";
330  *outStream << "} computed value: " << vals(i,j)
331  << " but reference value: " << nodalBasisValues[l] << "\n";
332  }
333  }
334  }
335 
336  // Check VALUE of basis functions at points: resize vals to rank-2 container:\n";
337  *outStream << " check VALUE of basis functions at points\n";
338  int numPoints = tetPoints.dimension(0);
339  vals.resize(numFields, numPoints);
340  vals.initialize(0.0);
341  tetBasis.getValues(vals, tetPoints, OPERATOR_VALUE);
342 
343  for (int i = 0; i < numFields; i++) {
344  for (int j = 0; j < numPoints; j++) {
345  int l = i + j * numFields;
346  if (std::abs(vals(i,j) - pointBasisValues[l]) > INTREPID_TOL) {
347  errorFlag++;
348  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
349 
350  // Output the multi-index of the value where the error is:
351  *outStream << " At multi-index { ";
352  *outStream << i << " ";*outStream << j << " ";
353  *outStream << "} computed value: " << vals(i,j)
354  << " but reference value: " << pointBasisValues[l] << "\n";
355  }
356  }
357  }
358 
359  // Check VALUE of basis functions at random points: resize vals to rank-2 container:\n";
360  *outStream << " check VALUE of basis functions at random points\n";
361  int numRandomPoints = 16384;
362  FieldContainer<double> tetRandomPoints(numRandomPoints, 3);
363  vals.resize(numFields, numRandomPoints);
364  vals.initialize(0.0);
365 
366  int point = 0;
367  int count = 0;
368  std::random_device rd;
369  std::mt19937 gen(rd());
370  std::uniform_real_distribution<> dis(0, 1);
371  while (point < numRandomPoints) {
372 
373  count++;
374  double r = dis(gen);
375  double s = dis(gen);
376  double t = dis(gen);
377  if (r + s + t > 1.0) continue;
378 
379  tetRandomPoints(point, 0) = r;
380  tetRandomPoints(point, 1) = s;
381  tetRandomPoints(point, 2) = t;
382 
383  point++;
384  }
385 
386  tetBasis.getValues(vals, tetRandomPoints, OPERATOR_VALUE);
387 
388  for (int j = 0; j < numRandomPoints; j++) {
389  double sum = 0.0;
390  for (int i = 0; i < numFields; i++) {
391  sum += vals(i,j);
392  }
393  if (std::abs(sum - 1.0) > INTREPID_TOL) {
394  errorFlag++;
395  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
396 
397  // Just indicate something bad happened
398  *outStream << " Composite tet basis functions";
399  *outStream << " are not summing to 1.0\n";
400  *outStream << " sum : " << sum << "\n";
401  }
402  }
403 
404  // Check GRAD of basis functions at points: resize vals to rank-3 container:\n";
405  numPoints = tetPoints.dimension(0);
406  vals.resize(numFields, numPoints, spaceDim);
407  vals.initialize(0.0);
408  tetBasis.getValues(vals, tetPoints, OPERATOR_GRAD);
409 
410  for (int i = 0; i < numFields; i++) {
411  for (int j = 0; j < numPoints; j++) {
412  for (int k = 0; k < spaceDim; k++) {
413  int l = k + i * spaceDim + j * spaceDim * numFields;
414  if (std::abs(vals(i,j,k) - pointBasisGrads[l]) > INTREPID_TOL) {
415  errorFlag++;
416  *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
417 
418  // Output the multi-index of the value where the error is:
419  *outStream << " At multi-index { ";
420  *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
421  *outStream << "} computed grad component: " << vals(i,j,k)
422  << " but reference grad component: " << pointBasisGrads[l] << "\n";
423  }
424  }
425  }
426  }
427  }
428  // Catch unexpected errors
429  catch (const std::logic_error & err) {
430  *outStream << err.what() << "\n\n";
431  errorFlag = -1000;
432  };
433 
434 
435  if (errorFlag != 0)
436  std::cout << "End Result: TEST FAILED\n";
437  else
438  std::cout << "End Result: TEST PASSED\n";
439 
440  // reset format state of std::cout
441  std::cout.copyfmt(oldFormatState);
442 
443  return errorFlag;
444 }
virtual int getCardinality() const
Returns cardinality of the basis.
int dimension(const int whichDim) const
Returns the specified dimension.
Header file for utility class to provide multidimensional containers.
void resize(const int dim0)
Resizes FieldContainer to a rank-1 container with the specified dimension, initialized by 0...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
void initialize(const Scalar value=0)
Initializes a field container by assigning value to all its elements.
Header file for the Intrepid::HGRAD_TET_COMP12_FEM class.
virtual const shards::CellTopology getBaseCellTopology() const
Returns the base cell topology for which the basis is defined. See Shards documentation http://trilin...