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