Intrepid
test_02.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
62int main(int argc, char *argv[]) {
63
64 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
65 // This little trick lets us print to std::cout only if
66 // a (dummy) command-line argument is provided.
67 int iprint = argc - 1;
68 Teuchos::RCP<std::ostream> outStream;
69 Teuchos::oblackholestream bhs; // outputs nothing
70 if (iprint > 0)
71 outStream = Teuchos::rcp(&std::cout, false);
72 else
73 outStream = Teuchos::rcp(&bhs, false);
74
75 // Save the format state of the original std::cout.
76 Teuchos::oblackholestream oldFormatState;
77 oldFormatState.copyfmt(std::cout);
78
79 *outStream \
80 << "===============================================================================\n" \
81 << "| |\n" \
82 << "| Unit Test (FunctionSpaceTools) |\n" \
83 << "| |\n" \
84 << "| 1) basic operator transformations and integration in HCURL |\n" \
85 << "| |\n" \
86 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
87 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
88 << "| |\n" \
89 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
90 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
91 << "| |\n" \
92 << "===============================================================================\n";
93
94
95 int errorFlag = 0;
96
97 typedef FunctionSpaceTools fst;
98
99 *outStream \
100 << "\n"
101 << "===============================================================================\n"\
102 << "| TEST 1: correctness of math operations |\n"\
103 << "===============================================================================\n";
104
105 outStream->precision(20);
106
107 try {
108 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >(); // cell type: hex
109
110 /* Related to cubature. */
111 DefaultCubatureFactory<double> cubFactory; // create cubature factory
112 int cubDegree = 20; // cubature degree
113 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature
114 int spaceDim = myCub->getDimension(); // get spatial dimension
115 int numCubPoints = myCub->getNumPoints(); // get number of cubature points
116 /* Related to basis. */
117 Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexBasis; // create H-curl basis on a hex
118 int numFields = hexBasis.getCardinality(); // get basis cardinality
119
120 /* Cell geometries and orientations. */
121 int numCells = 4;
122 int numNodes = 8;
123 int numCellData = numCells*numNodes*spaceDim;
124 int numSignData = numCells*numFields;
125
126 double hexnodes[] = {
127 // hex 0 -- affine
128 -1.0, -1.0, -1.0,
129 1.0, -1.0, -1.0,
130 1.0, 1.0, -1.0,
131 -1.0, 1.0, -1.0,
132 -1.0, -1.0, 1.0,
133 1.0, -1.0, 1.0,
134 1.0, 1.0, 1.0,
135 -1.0, 1.0, 1.0,
136 // hex 1 -- affine
137 -3.0, -3.0, 1.0,
138 6.0, 3.0, 1.0,
139 7.0, 8.0, 0.0,
140 -2.0, 2.0, 0.0,
141 -3.0, -3.0, 4.0,
142 6.0, 3.0, 4.0,
143 7.0, 8.0, 3.0,
144 -2.0, 2.0, 3.0,
145 // hex 2 -- affine
146 -3.0, -3.0, 0.0,
147 9.0, 3.0, 0.0,
148 15.0, 6.1, 0.0,
149 3.0, 0.1, 0.0,
150 9.0, 3.0, 0.1,
151 21.0, 9.0, 0.1,
152 27.0, 12.1, 0.1,
153 15.0, 6.1, 0.1,
154 // hex 3 -- nonaffine
155 -2.0, -2.0, 0.0,
156 2.0, -1.0, 0.0,
157 1.0, 6.0, 0.0,
158 -1.0, 1.0, 0.0,
159 0.0, 0.0, 1.0,
160 1.0, 0.0, 1.0,
161 1.0, 1.0, 1.0,
162 0.0, 1.0, 1.0
163 };
164
165 double edgesigns[] = {
166 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
167 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1,
168 -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1,
169 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1
170 };
171
172 /* Computational arrays. */
173 FieldContainer<double> cub_points(numCubPoints, spaceDim);
174 FieldContainer<double> cub_weights(numCubPoints);
175 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
176 FieldContainer<double> field_signs(numCells, numFields);
177 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
178 FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim);
179 FieldContainer<double> jacobian_det(numCells, numCubPoints);
180 FieldContainer<double> weighted_measure(numCells, numCubPoints);
181
182 FieldContainer<double> curl_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
183 FieldContainer<double> transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
184 FieldContainer<double> weighted_transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
185 FieldContainer<double> stiffness_matrices(numCells, numFields, numFields);
186
187 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim);
188 FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
189 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim);
190 FieldContainer<double> mass_matrices(numCells, numFields, numFields);
191
192 /******************* START COMPUTATION ***********************/
193
194 // get cubature points and weights
195 myCub->getCubature(cub_points, cub_weights);
196
197 // fill cell vertex array
198 cell_nodes.setValues(hexnodes, numCellData);
199 // set basis function signs, for each cell
200 field_signs.setValues(edgesigns, numSignData);
201
202 // compute geometric cell information
203 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
204 CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
205 CellTools<double>::setJacobianDet(jacobian_det, jacobian);
206 // compute weighted measure
207 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
208
209 // Computing stiffness matrices:
210 // tabulate curls of basis functions at (reference) cubature points
211 hexBasis.getValues(curl_of_basis_at_cub_points, cub_points, OPERATOR_CURL);
212
213 // transform curls of basis functions
214 fst::HCURLtransformCURL<double>(transformed_curl_of_basis_at_cub_points,
215 jacobian,
216 jacobian_det,
217 curl_of_basis_at_cub_points);
218 // multiply with weighted measure
219 fst::multiplyMeasure<double>(weighted_transformed_curl_of_basis_at_cub_points,
220 weighted_measure,
221 transformed_curl_of_basis_at_cub_points);
222
223 // we can apply the field signs to the basis function arrays, or after the fact, see below
224 fst::applyFieldSigns<double>(transformed_curl_of_basis_at_cub_points, field_signs);
225 fst::applyFieldSigns<double>(weighted_transformed_curl_of_basis_at_cub_points, field_signs);
226
227 // compute stiffness matrices
228 fst::integrate<double>(stiffness_matrices,
229 transformed_curl_of_basis_at_cub_points,
230 weighted_transformed_curl_of_basis_at_cub_points,
231 COMP_CPP);
232 // Computing mass matrices:
233 // tabulate values of basis functions at (reference) cubature points
234 hexBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE);
235
236 // transform values of basis functions
237 fst::HCURLtransformVALUE<double>(transformed_value_of_basis_at_cub_points,
238 jacobian_inv,
239 value_of_basis_at_cub_points);
240
241 // multiply with weighted measure
242 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points,
243 weighted_measure,
244 transformed_value_of_basis_at_cub_points);
245
246 // compute mass matrices
247 fst::integrate<double>(mass_matrices,
248 transformed_value_of_basis_at_cub_points,
249 weighted_transformed_value_of_basis_at_cub_points,
250 COMP_CPP);
251
252 // apply field signs (after the fact, as a post-processing step)
253 fst::applyLeftFieldSigns<double>(mass_matrices, field_signs);
254 fst::applyRightFieldSigns<double>(mass_matrices, field_signs);
255 /******************* STOP COMPUTATION ***********************/
256
257
258 /******************* START COMPARISON ***********************/
259 string basedir = "./testdata";
260 for (int cell_id = 0; cell_id < numCells-1; cell_id++) {
261
262 stringstream namestream;
263 string filename;
264 namestream << basedir << "/mass_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
265 namestream >> filename;
266
267 ifstream massfile(&filename[0]);
268 if (massfile.is_open()) {
269 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
270 errorFlag++;
271 massfile.close();
272 }
273 else {
274 errorFlag = -1;
275 std::cout << "End Result: TEST FAILED\n";
276 return errorFlag;
277 }
278
279 namestream.clear();
280 namestream << basedir << "/stiff_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
281 namestream >> filename;
282
283 ifstream stifffile(&filename[0]);
284 if (stifffile.is_open())
285 {
286 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
287 errorFlag++;
288 stifffile.close();
289 }
290 else {
291 errorFlag = -1;
292 std::cout << "End Result: TEST FAILED\n";
293 return errorFlag;
294 }
295
296 }
297
298 for (int cell_id = 3; cell_id < numCells; cell_id++) {
299
300 stringstream namestream;
301 string filename;
302 namestream << basedir << "/mass_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
303 namestream >> filename;
304
305 ifstream massfile(&filename[0]);
306 if (massfile.is_open()) {
307 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
308 errorFlag++;
309 massfile.close();
310 }
311 else {
312 errorFlag = -1;
313 std::cout << "End Result: TEST FAILED\n";
314 return errorFlag;
315 }
316
317 namestream.clear();
318 namestream << basedir << "/stiff_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
319 namestream >> filename;
320
321 ifstream stifffile(&filename[0]);
322 if (stifffile.is_open())
323 {
324 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
325 errorFlag++;
326 stifffile.close();
327 }
328 else {
329 errorFlag = -1;
330 std::cout << "End Result: TEST FAILED\n";
331 return errorFlag;
332 }
333
334 }
335
336 /******************* STOP COMPARISON ***********************/
337
338 *outStream << "\n";
339 }
340 catch (const std::logic_error & err) {
341 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
342 *outStream << err.what() << '\n';
343 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
344 errorFlag = -1000;
345 };
346
347 if (errorFlag != 0)
348 std::cout << "End Result: TEST FAILED\n";
349 else
350 std::cout << "End Result: TEST PASSED\n";
351
352 // reset format state of std::cout
353 std::cout.copyfmt(oldFormatState);
354
355 return errorFlag;
356}
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.
Header file for the Intrepid::HCURL_HEX_I1_FEM 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...