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