Intrepid
test_04.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
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) volume integration on tetrahedra, testing dataIntegral |\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
111 typedef FunctionSpaceTools fst;
112
113 *outStream \
114 << "\n"
115 << "===============================================================================\n"\
116 << "| TEST 1: correctness of cell volumes |\n"\
117 << "===============================================================================\n";
118
119 outStream->precision(20);
120
121 try {
122 shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tet
123
124 /* Related to cubature. */
125 DefaultCubatureFactory<double> cubFactory; // create cubature factory
126 int cubDegree = 0; // cubature degree
127 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature
128 int spaceDim = myCub->getDimension(); // get spatial dimension
129 int numCubPoints = myCub->getNumPoints(); // get number of cubature points
130
131 /* Cell geometries. */
132 int numCells = 4;
133 int numNodes = 4;
134 int numCellData = numCells*numNodes*spaceDim;
135 double tetnodes[] = {
136 // tet 0
137 0.0, 0.0, 0.0,
138 1.0, 0.0, 0.0,
139 0.0, 1.0, 0.0,
140 0.0, 0.0, 1.0,
141 // tet 1
142 4.0, 5.0, 1.0,
143 -6.0, 2.0, 0.0,
144 4.0, -3.0, -1.0,
145 0.0, 2.0, 5.0,
146 // tet 2
147 -6.0, -3.0, 1.0,
148 9.0, 2.0, 1.0,
149 8.9, 2.1, 0.9,
150 8.9, 2.1, 1.1,
151 // tet 3
152 -6.0, -3.0, 1.0,
153 12.0, 3.0, 1.0,
154 2.9, 0.1, 0.9,
155 2.9, 0.1, 1.1
156 };
157
158 /* Analytic volumes. */
159 double tetvols[] = {1.0/6.0, 194.0/3.0, 1.0/15.0, 2.0/25.0};
160
161 /* Computational arrays. */
162 FieldContainer<double> cub_points(numCubPoints, spaceDim);
163 FieldContainer<double> cub_weights(numCubPoints);
164 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
165 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
166 FieldContainer<double> jacobian_det(numCells, numCubPoints);
167 FieldContainer<double> weighted_measure(numCells, numCubPoints);
168 FieldContainer<double> data_one(numCells, numCubPoints);
169 FieldContainer<double> volumes(numCells);
170
171 /******************* START COMPUTATION ***********************/
172
173 // get cubature points and weights
174 myCub->getCubature(cub_points, cub_weights);
175
176 // fill cell vertex array
177 cell_nodes.setValues(tetnodes, numCellData);
178
179 // compute geometric cell information
180 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
181 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
182
183 // compute weighted measure
184 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
185
186 // set data to 1.0
187 for (int cell=0; cell<data_one.dimension(0); cell++) {
188 for (int qp=0; qp<data_one.dimension(1); qp++) {
189 data_one(cell,qp) = 1.0;
190 }
191 }
192
193 // compute volumes
194 fst::integrate<double>(volumes, data_one, weighted_measure, COMP_CPP);
195
196 /******************* STOP COMPUTATION ***********************/
197
198
199 /******************* START COMPARISON ***********************/
200 for (int cell_id = 0; cell_id < numCells; cell_id++) {
201 *outStream << "Volume of cell " << cell_id << " = " << volumes(cell_id) << " vs. Analytic value = " << tetvols[cell_id] << "\n";
202 if (std::fabs(volumes(cell_id)-tetvols[cell_id]) > INTREPID_TOL) {
203 errorFlag++;
204 }
205 }
206 /******************* STOP COMPARISON ***********************/
207
208 *outStream << "\n";
209 }
210 catch (const std::logic_error & err) {
211 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
212 *outStream << err.what() << '\n';
213 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
214 errorFlag = -1000;
215 };
216
217
218 if (errorFlag != 0)
219 std::cout << "End Result: TEST FAILED\n";
220 else
221 std::cout << "End Result: TEST PASSED\n";
222
223 // reset format state of std::cout
224 std::cout.copyfmt(oldFormatState);
225
226 return errorFlag;
227}
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.
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...