Intrepid
test_06.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
53#include "Intrepid_HGRAD_TET_C1_FEM.hpp"
55#include "Intrepid_Utils.hpp"
56#include "Teuchos_oblackholestream.hpp"
57#include "Teuchos_RCP.hpp"
58#include "Teuchos_GlobalMPISession.hpp"
59
60using namespace std;
61using namespace Intrepid;
62
63#define INTREPID_TEST_COMMAND( S ) \
64{ \
65 try { \
66 S ; \
67 } \
68 catch (const std::logic_error & err) { \
69 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
70 *outStream << err.what() << '\n'; \
71 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
72 }; \
73}
74
75
76int main(int argc, char *argv[]) {
77
78 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
79
80 // This little trick lets us print to std::cout only if
81 // a (dummy) command-line argument is provided.
82 int iprint = argc - 1;
83 Teuchos::RCP<std::ostream> outStream;
84 Teuchos::oblackholestream bhs; // outputs nothing
85 if (iprint > 0)
86 outStream = Teuchos::rcp(&std::cout, false);
87 else
88 outStream = Teuchos::rcp(&bhs, false);
89
90 // Save the format state of the original std::cout.
91 Teuchos::oblackholestream oldFormatState;
92 oldFormatState.copyfmt(std::cout);
93
94 *outStream \
95 << "===============================================================================\n" \
96 << "| |\n" \
97 << "| Unit Test (FunctionSpaceTools) |\n" \
98 << "| |\n" \
99 << "| 1) basic operator transformations and integration in HGRAD |\n" \
100 << "| |\n" \
101 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
102 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
103 << "| |\n" \
104 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
105 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
106 << "| |\n" \
107 << "===============================================================================\n";
108
109
110 int errorFlag = 0;
111#ifdef HAVE_INTREPID_DEBUG
112 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
113 int endThrowNumber = beginThrowNumber + 28;
114#endif
115
116 typedef FunctionSpaceTools fst;
117
118 *outStream \
119 << "\n"
120 << "===============================================================================\n"\
121 << "| TEST 1: exceptions |\n"\
122 << "===============================================================================\n";
123
124 try{
125#ifdef HAVE_INTREPID_DEBUG
126 FieldContainer<double, 1> a_2(2);
127 FieldContainer<double, 2> a_2_2(2, 2);
128 FieldContainer<double, 3> a_2_3(2, 3);
129 FieldContainer<double, 4> a_3_2(3, 2);
130 FieldContainer<double, 5> a_2_2_3(2, 2, 3);
131 FieldContainer<double, 6> a_2_2_3_3(2, 2, 3, 3);
132 FieldContainer<double, 7> a_2_2_2(2, 2, 2);
133 FieldContainer<double, 8> a_2_2_2_3_3(2, 2, 2, 3, 3);
134 FieldContainer<double, 9> a_2_2_2_2_2(2, 2, 2, 2, 2);
135 FieldContainer<double, 10> a_2_2_2_2(2, 2, 2, 2);
136 FieldContainer<double, 11> a_3_2_2_2(3, 2, 2, 2);
137 FieldContainer<double, 12> a_2_3_2_2(2, 3, 2, 2);
138 FieldContainer<double, 13> a_2_2_3_2(2, 2, 3, 2);
139 FieldContainer<double, 14> a_2_2_2_3(2, 2, 2, 3);
140
141 *outStream << "\n >>>>> TESTING computeCellMeasure:\n";
142 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2, a_2) );
143 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2_2, a_2) );
144
145 *outStream << "\n >>>>> TESTING computeFaceMeasure:\n";
146 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
147 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2_2_3_3, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) );
148
149 *outStream << "\n >>>>> TESTING computeEdgeMeasure:\n";
150 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
151 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2_2_2_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) );
152
153 *outStream << "\n >>>>> TESTING integrate:\n";
154 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
155 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
156 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
157 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
158 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
159 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
160 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
161 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
162 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
163 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
164
165 *outStream << "\n >>>>> TESTING operatorIntegral:\n";
166 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2, a_2_2_2, COMP_CPP) );
167 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) );
168 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) );
169 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
170
171 *outStream << "\n >>>>> TESTING functionalIntegral:\n";
172 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_2_3_3, a_2_2_2, COMP_CPP) );
173 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) );
174 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) );
175 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) );
176
177 *outStream << "\n >>>>> TESTING dataIntegral:\n";
178 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2, a_2_2_2, COMP_CPP) );
179 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2, a_2_2, COMP_CPP) );
180 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) );
181 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) );
182
183 *outStream << "\n >>>>> TESTING applyLeftFieldSigns:\n";
184 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2, a_2) );
185 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2) );
186 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_3_2) );
187 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_3) );
188 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_2) );
189
190 *outStream << "\n >>>>> TESTING applyRightFieldSigns:\n";
191 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2, a_2) );
192 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2) );
193 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_3_2) );
194 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_3) );
195 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_2) );
196
197 *outStream << "\n >>>>> TESTING applyFieldSigns:\n";
198 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2, a_2) );
199 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2) );
200 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_3_2) );
201 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2_3) );
202 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2_2_3_3, a_2_2) );
203
204 *outStream << "\n >>>>> TESTING evaluate:\n";
205 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2) );
206 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2_2_3_3) );
207 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2_2, a_2_2_2_3_3) );
208 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_3_2, a_2_2_2_3_3) );
209 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_2_3, a_2_2_2_3_3) );
210 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_3_2_2_2, a_2_2, a_2_2_2_2_2) );
211 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_3_2_2, a_2_2, a_2_2_2_2_2) );
212 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_2, a_2_2, a_2_2_2_2_2) );
213 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_3, a_2_2, a_2_2_2_2_2) );
214 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_2, a_2_2, a_2_2_2_2_2) );
215#endif
216 }
217 catch (const std::logic_error & err) {
218 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
219 *outStream << err.what() << '\n';
220 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
221 errorFlag = -1000;
222 };
223
224#ifdef HAVE_INTREPID_DEBUG
225 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
226 errorFlag++;
227#endif
228
229 *outStream \
230 << "\n"
231 << "===============================================================================\n"\
232 << "| TEST 2: correctness of math operations |\n"\
233 << "===============================================================================\n";
234
235 outStream->precision(20);
236
237 try {
238 // cell type: tet
239 shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >();
240
241 /* Related to cubature. */
242 // create cubature factory
243 DefaultCubatureFactory<double, FieldContainer<double, 1>, FieldContainer<double, 2> > cubFactory;
244 // cubature degree
245 int cubDegree = 2;
246 // create default cubature
247 Teuchos::RCP<Cubature<double, FieldContainer<double, 1>, FieldContainer<double, 2> > > myCub = cubFactory.create(cellType, cubDegree);
248 // get spatial dimension
249 int spaceDim = myCub->getDimension();
250 // get number of cubature points
251 int numCubPoints = myCub->getNumPoints();
252
253 /* Related to basis. */
254 // create tet basis
255 Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double, 1> > tetBasis;
256 // get basis cardinality
257 int numFields = tetBasis.getCardinality();
258
259 /* Cell geometries. */
260 int numCells = 4;
261 int numNodes = 4;
262 int numCellData = numCells*numNodes*spaceDim;
263 double tetnodes[] = {
264 // tet 0
265 0.0, 0.0, 0.0,
266 1.0, 0.0, 0.0,
267 0.0, 1.0, 0.0,
268 0.0, 0.0, 1.0,
269 // tet 1
270 4.0, 5.0, 1.0,
271 -6.0, 2.0, 0.0,
272 4.0, -3.0, -1.0,
273 0.0, 2.0, 5.0,
274 // tet 2
275 -6.0, -3.0, 1.0,
276 9.0, 2.0, 1.0,
277 8.9, 2.1, 0.9,
278 8.9, 2.1, 1.1,
279 // tet 3
280 -6.0, -3.0, 1.0,
281 12.0, 3.0, 1.0,
282 2.9, 0.1, 0.9,
283 2.9, 0.1, 1.1
284 };
285
286 /* Computational arrays. */
287 FieldContainer<double, 1> cub_points(numCubPoints, spaceDim);
288 FieldContainer<double, 2> cub_weights(numCubPoints);
289 FieldContainer<double, 3> cell_nodes(numCells, numNodes, spaceDim);
290 FieldContainer<double, 4> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
291 FieldContainer<double, 5> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
292 FieldContainer<double, 6> jacobian_det(numCells, numCubPoints);
293 FieldContainer<double, 7> weighted_measure(numCells, numCubPoints);
294
295 FieldContainer<double, 1> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
296 FieldContainer<double, 9> transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
297 FieldContainer<double, 10> weighted_transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
298 FieldContainer<double, 11> stiffness_matrices(numCells, numFields, numFields);
299
300 FieldContainer<double, 1> value_of_basis_at_cub_points(numFields, numCubPoints);
301 FieldContainer<double, 13> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints);
302 FieldContainer<double, 14> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints);
303 FieldContainer<double, 15> mass_matrices(numCells, numFields, numFields);
304
305 /******************* START COMPUTATION ***********************/
306
307 // get cubature points and weights
308 myCub->getCubature(cub_points, cub_weights);
309
310 // fill cell vertex array
311 cell_nodes.setValues(tetnodes, numCellData);
312
313 // compute geometric cell information
314 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
315 CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
316 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
317
318 // compute weighted measure
319 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
320
321 // Computing stiffness matrices:
322 // tabulate gradients of basis functions at (reference) cubature points
323 tetBasis.getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD);
324
325 // transform gradients of basis functions
326 fst::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points,
327 jacobian_inv,
328 grad_of_basis_at_cub_points);
329
330 // multiply with weighted measure
331 fst::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points,
332 weighted_measure,
333 transformed_grad_of_basis_at_cub_points);
334
335 // compute stiffness matrices
336 fst::integrate<double>(stiffness_matrices,
337 transformed_grad_of_basis_at_cub_points,
338 weighted_transformed_grad_of_basis_at_cub_points,
339 COMP_CPP);
340
341
342 // Computing mass matrices:
343 // tabulate values of basis functions at (reference) cubature points
344 tetBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
345
346 // transform values of basis functions
347 fst::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
348 value_of_basis_at_cub_points);
349
350 // multiply with weighted measure
351 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
352 weighted_measure,
353 transformed_value_of_basis_at_cub_points);
354
355 // compute mass matrices
356 fst::integrate<double>(mass_matrices,
357 transformed_value_of_basis_at_cub_points,
358 weighted_transformed_value_of_basis_at_cub_points,
359 COMP_CPP);
360
361 /******************* STOP COMPUTATION ***********************/
362
363
364 /******************* START COMPARISON ***********************/
365 string basedir = "./testdata";
366 for (int cell_id = 0; cell_id < numCells; cell_id++) {
367
368 stringstream namestream;
369 string filename;
370 namestream << basedir << "/mass_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat";
371 namestream >> filename;
372
373 ifstream massfile(&filename[0]);
374 if (massfile.is_open()) {
375 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, 0) > 0)
376 errorFlag++;
377 massfile.close();
378 }
379 else {
380 errorFlag = -1;
381 std::cout << "End Result: TEST FAILED\n";
382 return errorFlag;
383 }
384
385 namestream.clear();
386 namestream << basedir << "/stiff_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat";
387 namestream >> filename;
388
389 ifstream stifffile(&filename[0]);
390 if (stifffile.is_open())
391 {
392 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, 0) > 0)
393 errorFlag++;
394 stifffile.close();
395 }
396 else {
397 errorFlag = -1;
398 std::cout << "End Result: TEST FAILED\n";
399 return errorFlag;
400 }
401
402 }
403
404 /******************* STOP COMPARISON ***********************/
405
406 *outStream << "\n";
407 }
408 catch (const std::logic_error & err) {
409 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
410 *outStream << err.what() << '\n';
411 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
412 errorFlag = -1000;
413 };
414
415
416 if (errorFlag != 0)
417 std::cout << "End Result: TEST FAILED\n";
418 else
419 std::cout << "End Result: TEST PASSED\n";
420
421 // reset format state of std::cout
422 std::cout.copyfmt(oldFormatState);
423
424 return errorFlag;
425}
Header file for the Intrepid::CellTools class.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for utility class to provide multidimensional containers.
Header file for the Intrepid::FunctionSpaceTools class.
Intrepid utilities.
A stateless class for operations on cell data. Provides methods for:
static void setJacobianDet(ArrayJacDet &jacobianDet, const ArrayJac &jacobian)
Computes the determinant of the Jacobian matrix DF of the reference-to-physical frame map F.
static void setJacobianInv(ArrayJacInv &jacobianInv, const ArrayJac &jacobian)
Computes the inverse of the Jacobian matrix DF of the reference-to-physical frame map F.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...