Intrepid
test_01.cpp
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#include "Teuchos_oblackholestream.hpp"
51#include "Teuchos_RCP.hpp"
52#include "Teuchos_GlobalMPISession.hpp"
53
54#include <random>
55
56using namespace std;
57using 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
73int 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
111 Basis_HGRAD_TET_COMP12_FEM<double, FieldContainer<double> > tetBasis;
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
313 FieldContainer<double> vals;
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}
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::HGRAD_TET_COMP12_FEM class.