Intrepid
test_05.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 "Shards_Array.hpp"
58#include "Intrepid_Utils.hpp"
59#include "Teuchos_oblackholestream.hpp"
60#include "Teuchos_RCP.hpp"
61#include "Teuchos_GlobalMPISession.hpp"
62
63using namespace std;
64using namespace Intrepid;
65
66SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell )
67SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell )
68
69SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field )
70SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field )
71
72SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point )
73SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point )
74
75SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim )
76SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim )
77
78SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Node )
79SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Node )
80
81
82typedef FieldContainer<double> FC;
83
84
85int main(int argc, char *argv[]) {
86
87 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
88
89 // This little trick lets us print to std::cout only if
90 // a (dummy) command-line argument is provided.
91 int iprint = argc - 1;
92 Teuchos::RCP<std::ostream> outStream;
93 Teuchos::oblackholestream bhs; // outputs nothing
94 if (iprint > 0)
95 outStream = Teuchos::rcp(&std::cout, false);
96 else
97 outStream = Teuchos::rcp(&bhs, false);
98
99 // Save the format state of the original std::cout.
100 Teuchos::oblackholestream oldFormatState;
101 oldFormatState.copyfmt(std::cout);
102
103 *outStream \
104 << "===============================================================================\n" \
105 << "| |\n" \
106 << "| Unit Test (FunctionSpaceTools) |\n" \
107 << "| |\n" \
108 << "| 1) basic operator transformations and integration in HDIV |\n" \
109 << "| via shards::Array wrappers |\n" \
110 << "| |\n" \
111 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
112 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
113 << "| |\n" \
114 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
115 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
116 << "| |\n" \
117 << "===============================================================================\n";
118
119
120 int errorFlag = 0;
121
122 typedef FunctionSpaceTools fst;
123
124 *outStream \
125 << "\n"
126 << "===============================================================================\n"\
127 << "| TEST 1: correctness of math operations |\n"\
128 << "===============================================================================\n";
129
130 outStream->precision(20);
131
132 try {
133 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >(); // cell type: hex
134
135 /* Related to cubature. */
136 DefaultCubatureFactory<double> cubFactory; // create cubature factory
137 int cubDegree = 20; // cubature degree
138 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature
139 int spaceDim = myCub->getDimension(); // get spatial dimension
140 int numCubPoints = myCub->getNumPoints(); // get number of cubature points
141
142 /* Related to basis. */
143 Basis_HDIV_HEX_I1_FEM<double, FieldContainer<double> > hexBasis; // create H-div basis on a hex
144 int numFields = hexBasis.getCardinality(); // get basis cardinality
145
146 /* Cell geometries and orientations. */
147 int numCells = 4;
148 int numNodes = 8;
149 int numCellData = numCells*numNodes*spaceDim;
150 int numSignData = numCells*numFields;
151
152 double hexnodes[] = {
153 // hex 0 -- affine
154 -1.0, -1.0, -1.0,
155 1.0, -1.0, -1.0,
156 1.0, 1.0, -1.0,
157 -1.0, 1.0, -1.0,
158 -1.0, -1.0, 1.0,
159 1.0, -1.0, 1.0,
160 1.0, 1.0, 1.0,
161 -1.0, 1.0, 1.0,
162 // hex 1 -- affine
163 -3.0, -3.0, 1.0,
164 6.0, 3.0, 1.0,
165 7.0, 8.0, 0.0,
166 -2.0, 2.0, 0.0,
167 -3.0, -3.0, 4.0,
168 6.0, 3.0, 4.0,
169 7.0, 8.0, 3.0,
170 -2.0, 2.0, 3.0,
171 // hex 2 -- affine
172 -3.0, -3.0, 0.0,
173 9.0, 3.0, 0.0,
174 15.0, 6.1, 0.0,
175 3.0, 0.1, 0.0,
176 9.0, 3.0, 0.1,
177 21.0, 9.0, 0.1,
178 27.0, 12.1, 0.1,
179 15.0, 6.1, 0.1,
180 // hex 3 -- nonaffine
181 -2.0, -2.0, 0.0,
182 2.0, -1.0, 0.0,
183 1.0, 6.0, 0.0,
184 -1.0, 1.0, 0.0,
185 0.0, 0.0, 1.0,
186 1.0, 0.0, 1.0,
187 1.0, 1.0, 1.0,
188 0.0, 1.0, 1.0
189 };
190
191 short facesigns[] = {
192 1, 1, 1, 1, 1, 1,
193 1, -1, 1, -1, 1, -1,
194 -1, -1, 1, 1, -1, 1,
195 -1, -1, 1, 1, -1, -1
196 };
197
198 /* Computational arrays. */
199 // First allocate one very large work space.
200 Teuchos::Array<double> work_space(numCubPoints*spaceDim +
201 numCubPoints +
202 numCells*numNodes*spaceDim +
203 numCells*numCubPoints*spaceDim*spaceDim +
204 2*numCells*numCubPoints +
205 numFields*numCubPoints +
206 2*numCells*numFields*numCubPoints +
207 numCells*numFields*numFields +
208 numFields*numCubPoints*spaceDim +
209 2*numCells*numFields*numCubPoints*spaceDim +
210 numCells*numFields*numFields
211 );
212
213 int offset = 0;
214 shards::Array<double,shards::NaturalOrder,Point,Dim> cub_points(&work_space[offset], numCubPoints, spaceDim);
215 FC cub_points_FC(cub_points);
216 offset += numCubPoints*spaceDim;
217 shards::Array<double,shards::NaturalOrder,Point> cub_weights(&work_space[offset], numCubPoints);
218 FC cub_weights_FC(cub_weights);
219 offset += numCubPoints;
220 shards::Array<double,shards::NaturalOrder,Cell,Node,Dim> cell_nodes(&work_space[offset], numCells, numNodes, spaceDim);
221 FC cell_nodes_FC(cell_nodes);
222 offset += numCells*numNodes*spaceDim;
223 FieldContainer<short> field_signs(numCells, numFields);
224 shards::Array<double,shards::NaturalOrder,Cell,Point,Dim,Dim> jacobian(&work_space[offset], numCells, numCubPoints, spaceDim, spaceDim);
225 FC jacobian_FC(jacobian);
226 offset += numCells*numCubPoints*spaceDim*spaceDim;
227 //shards::Array<double,shards::NaturalOrder> jacobian_inv(&work_space[offset]numCells, numCubPoints, spaceDim, spaceDim);
228 shards::Array<double,shards::NaturalOrder,Cell,Point> jacobian_det(&work_space[offset], numCells, numCubPoints);
229 FC jacobian_det_FC(jacobian_det);
230 offset += numCells*numCubPoints;
231 shards::Array<double,shards::NaturalOrder,Cell,Point> weighted_measure(&work_space[offset], numCells, numCubPoints);
232 FC weighted_measure_FC(weighted_measure);
233 offset += numCells*numCubPoints;
234
235 shards::Array<double,shards::NaturalOrder,Field,Point> div_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints);
236 FC div_of_basis_at_cub_points_FC(div_of_basis_at_cub_points);
237 offset += numFields*numCubPoints;
238 shards::Array<double,shards::NaturalOrder,Cell,Field,Point>
239 transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints);
240 FC transformed_div_of_basis_at_cub_points_FC(transformed_div_of_basis_at_cub_points);
241 offset += numCells*numFields*numCubPoints;
242 shards::Array<double,shards::NaturalOrder,Cell,Field,Point>
243 weighted_transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints);
244 FC weighted_transformed_div_of_basis_at_cub_points_FC(weighted_transformed_div_of_basis_at_cub_points);
245 offset += numCells*numFields*numCubPoints;
246 shards::Array<double,shards::NaturalOrder,Cell,Field,Field> stiffness_matrices(&work_space[offset], numCells, numFields, numFields);
247 FC stiffness_matrices_FC(stiffness_matrices);
248 offset += numCells*numFields*numFields;
249
250 shards::Array<double,shards::NaturalOrder,Field,Point,Dim> value_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints, spaceDim);
251 FC value_of_basis_at_cub_points_FC(value_of_basis_at_cub_points);
252 offset += numFields*numCubPoints*spaceDim;
253 shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim>
254 transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim);
255 FC transformed_value_of_basis_at_cub_points_FC(transformed_value_of_basis_at_cub_points);
256 offset += numCells*numFields*numCubPoints*spaceDim;
257 shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim>
258 weighted_transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim);
259 FC weighted_transformed_value_of_basis_at_cub_points_FC(weighted_transformed_value_of_basis_at_cub_points);
260 offset += numCells*numFields*numCubPoints*spaceDim;
261 shards::Array<double,shards::NaturalOrder,Cell,Field,Field> mass_matrices(&work_space[offset], numCells, numFields, numFields);
262 FC mass_matrices_FC(mass_matrices);
263
264
265 /******************* START COMPUTATION ***********************/
266
267 // get cubature points and weights
268 myCub->getCubature(cub_points_FC, cub_weights_FC);
269
270 // fill cell vertex array
271 cell_nodes_FC.setValues(hexnodes, numCellData);
272
273 // set basis function signs, for each cell
274 field_signs.setValues(facesigns, numSignData);
275
276 // compute geometric cell information
277 CellTools<double>::setJacobian(jacobian_FC, cub_points_FC, cell_nodes_FC, cellType);
278 //CellTools<double>::setJacobianInv(jacobian_inv, jacobian);
279 CellTools<double>::setJacobianDet(jacobian_det_FC, jacobian_FC);
280
281 // compute weighted measure
282 fst::computeCellMeasure<double>(weighted_measure_FC, jacobian_det_FC, cub_weights_FC);
283
284
285 // Computing stiffness matrices:
286 // tabulate divergences of basis functions at (reference) cubature points
287 hexBasis.getValues(div_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_DIV);
288
289 // transform divergences of basis functions
290 fst::HDIVtransformDIV<double>(transformed_div_of_basis_at_cub_points_FC,
291 jacobian_det_FC,
292 div_of_basis_at_cub_points_FC);
293
294 // multiply with weighted measure
295 fst::multiplyMeasure<double>(weighted_transformed_div_of_basis_at_cub_points_FC,
296 weighted_measure_FC,
297 transformed_div_of_basis_at_cub_points_FC);
298
299 // compute stiffness matrices
300 fst::integrate<double>(stiffness_matrices_FC,
301 transformed_div_of_basis_at_cub_points_FC,
302 weighted_transformed_div_of_basis_at_cub_points_FC,
303 COMP_CPP);
304
305 // apply field signs
306 fst::applyLeftFieldSigns<double>(stiffness_matrices_FC, field_signs);
307 fst::applyRightFieldSigns<double>(stiffness_matrices_FC, field_signs);
308
309
310 // Computing mass matrices:
311 // tabulate values of basis functions at (reference) cubature points
312 hexBasis.getValues(value_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_VALUE);
313
314 // transform values of basis functions
315 fst::HDIVtransformVALUE<double>(transformed_value_of_basis_at_cub_points_FC,
316 jacobian_FC,
317 jacobian_det_FC,
318 value_of_basis_at_cub_points_FC);
319
320 // multiply with weighted measure
321 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_FC,
322 weighted_measure_FC,
323 transformed_value_of_basis_at_cub_points_FC);
324
325 // compute mass matrices
326 fst::integrate<double>(mass_matrices_FC,
327 transformed_value_of_basis_at_cub_points_FC,
328 weighted_transformed_value_of_basis_at_cub_points_FC,
329 COMP_CPP);
330
331 // apply field signs
332 fst::applyLeftFieldSigns<double>(mass_matrices_FC, field_signs);
333 fst::applyRightFieldSigns<double>(mass_matrices_FC, field_signs);
334
335 /******************* STOP COMPUTATION ***********************/
336
337
338 /******************* START COMPARISON ***********************/
339 string basedir = "./testdata";
340 for (int cell_id = 0; cell_id < numCells-1; cell_id++) {
341
342 stringstream namestream;
343 string filename;
344 namestream << basedir << "/mass_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
345 namestream >> filename;
346
347 ifstream massfile(&filename[0]);
348 if (massfile.is_open()) {
349 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0)
350 errorFlag++;
351 massfile.close();
352 }
353 else {
354 errorFlag = -1;
355 std::cout << "End Result: TEST FAILED\n";
356 return errorFlag;
357 }
358
359 namestream.clear();
360 namestream << basedir << "/stiff_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
361 namestream >> filename;
362
363 ifstream stifffile(&filename[0]);
364 if (stifffile.is_open())
365 {
366 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0)
367 errorFlag++;
368 stifffile.close();
369 }
370 else {
371 errorFlag = -1;
372 std::cout << "End Result: TEST FAILED\n";
373 return errorFlag;
374 }
375
376 }
377 for (int cell_id = 3; cell_id < numCells; cell_id++) {
378
379 stringstream namestream;
380 string filename;
381 namestream << basedir << "/mass_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
382 namestream >> filename;
383
384 ifstream massfile(&filename[0]);
385 if (massfile.is_open()) {
386 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
387 errorFlag++;
388 massfile.close();
389 }
390 else {
391 errorFlag = -1;
392 std::cout << "End Result: TEST FAILED\n";
393 return errorFlag;
394 }
395
396 namestream.clear();
397 namestream << basedir << "/stiff_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat";
398 namestream >> filename;
399
400 ifstream stifffile(&filename[0]);
401 if (stifffile.is_open())
402 {
403 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0)
404 errorFlag++;
405 stifffile.close();
406 }
407 else {
408 errorFlag = -1;
409 std::cout << "End Result: TEST FAILED\n";
410 return errorFlag;
411 }
412
413 }
414 /******************* STOP COMPARISON ***********************/
415
416 *outStream << "\n";
417 }
418 catch (const std::logic_error & err) {
419 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
420 *outStream << err.what() << '\n';
421 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
422 errorFlag = -1000;
423 };
424
425
426 if (errorFlag != 0)
427 std::cout << "End Result: TEST FAILED\n";
428 else
429 std::cout << "End Result: TEST PASSED\n";
430
431 // reset format state of std::cout
432 std::cout.copyfmt(oldFormatState);
433
434 return errorFlag;
435}
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.
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.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...