Intrepid
test_01.cpp
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 
51 #include "Intrepid_CellTools.hpp"
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 
59 using namespace std;
60 using 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 
75 int 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
125  FieldContainer<double> a_2(2);
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 }
virtual int getCardinality() const
Returns cardinality of the basis.
Header file for the Intrepid::CellTools class.
Header file for utility class to provide multidimensional containers.
Defines expert-level interfaces for the evaluation of functions and operators in physical space (supp...
void getValues(ArrayScalar &outputValues, const ArrayScalar &inputPoints, const EOperator operatorType) const
FEM basis evaluation on a reference Tetrahedron cell.
Intrepid utilities.
Header file for the abstract base class Intrepid::DefaultCubatureFactory.
Header file for the Intrepid::FunctionSpaceTools class.
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.
A stateless class for operations on cell data. Provides methods for:
Implementation of the default H(grad)-compatible FEM basis of degree 1 on Tetrahedron cell...